function Res = KochCube(n,KColorMap,Filename_res)
% This function computes the extension of the quadratic type 1 curve of
% Koch snowflake
% It displays the result in MatLab and also generates a 3D obj file.
% To view the obj file result, you can use any 3D file viewer, like MeshLab
%
% Input :
%  - n [Optional] : Number of recusion levels
%  - KColorMap [Optional] : String containing the name of the colormap to
%        be used for representation
%  - Filename_res [Optional] : String containing the name of the obj file
%        to be created, and also the png file.
% Output :
% - Res : Matrix containing the information to represent the Koch cube
%         structure.
%         Each row corresponds to a square to be plotted in 3D.
% The three first columns correspond to the center a square.
% The fourth column is the size of the each square.
% The fifth column to the seventh column give the normal of the square
% The eighth column corresponds to a color level. It is used for
% representation purpose.
%
% Guillaume Jacquenot
% guillaume dot jacquenot at gmail dot com
% 2010 10 03

% Define default values if varaibles do not exist
if ~(exist('n','var')==1)
    n = 3;
end
if ~(exist('KColorMap','var')==1)
    KColorMap = 'gray';
end
if nargin < 3
    % Generate a filename for results
    Filename_res = ['KochCube_' sprintf('%02d',n) '_' KColorMap ];
    if nargin < 1
        help(mfilename);
    end
end
% Initial pattern. The content of each column is the same as the result
% matrix, whose description is given previously.
Ini = [0.5 0.5 0 1 0 0 1 1];
% Apply a multiplying factor on the size and location of the initial
% pattern to work with integer on vertices.
Ini(1:4) = Ini(1:4) * 3^n;

% Define initial pattern "Ini" as "pat"
pat = Ini;
% Loop on the recursion level
for i=1:n
    % Allocate memory for each temporary result / level
    tmp = zeros(13^i,8);
    % Perform one refinement iteration for each face
    for j=1:size(pat,1)
        tmp((13*(j-1)+1):(13*j),:) = KochCube_OneIter(pat(j,:));
    end
    % Set the temporary result "tmp" as "pat"
    pat = tmp;
end
% Set the result of the final level to be "Res"
Res = pat;

if ~iscolormap(KColorMap)
    error([mfilename ':e0'],...
        'KColorMap has to be a valid colormap. See ');
end

% Write the result to an obj file (Use MeshLab, Blender to visualize the
% result
WriteKochCube(Res,n,KColorMap,Filename_reiis);
% Display the result
Res2 = Res;
% Trick to have the same bounding box whatever the recursion level asks.
Res2(:,1:4) = Res2(:,1:4)*3^(5-n);
Faxis = [0 243 0 243 0 125];
DisplayKochCube(Res2,n,KColorMap,Filename_res,Faxis);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Res = KochCube_OneIter(Ini)
% Performs one iteration
if nargin == 0
    Ini = [4.5 4.5 0 9 0 0 1 1];
end
% Define the size of the resulting matrix
Res          = zeros(13,8);
Normale      = Ini(5:7);
% Compute the central square
Res(1,:)     = Ini;
Res(1,1:3)   = Res(1,1:3) + Normale*Ini(4)/3;
Res(1,8)     = Res(1,8) + 2;
% Compute the 8 squares located around the central square
Res(2:9,:)   = repmat(Ini,8,1);
Res(2:9,1:3) = Res(2:9,1:3)+ ...
               Generate_pattern_matrix(Normale)*Ini(4)/3;
% Compute the 4 squares located around the central square, and whose normal
% will change
Res(10:13,4)   = Ini(4);
Res(10:13,5:7) = Generate_pattern_matrix(Normale,1);
Res(10:13,1:3) = repmat(Ini(1:3),4,1) +...
                    Res(10:13,5:7)*Ini(4)/6 + ...
                    repmat(Ini(5:7),4,1)*Ini(4)/6;
Res(10:13,8)   = Ini(8)+1;
Res(:,4)       = Res(:,4)/3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function M = Generate_pattern_matrix(Normale,FirstHalf)
% Normale :
% The second argument is used to create or not the first half of the
% matrix M (ie the 4st rows).
M = [1,0,0;-1,0,0;0, 1,0; 0,-1,0;...
     1,1,0;-1,1,0;-1,-1,0;1,-1,0;];

if abs(Normale(1))==1
    M(:,3) = M(:,1);
    M(:,1) = 0;
elseif abs(Normale(2))==1
    M(:,3) = M(:,2);
    M(:,2) = 0;
end
if nargin==2 && FirstHalf
    M(5:end,:)= [];
elseif nargin==2 && ~FirstHalf
    M(1:4,:)= [];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H = DisplayKochCube(Res,n,KColorMap,Filename_res,Faxis)
% Display the Koch Cube
H = figure;
hold on
axis equal
[Vertices,Faces,ColorFaces] = Generate3DKochCube(Res,n,KColorMap);
if n>0
    patch('Faces',Faces','Vertices',Vertices,...
          'FaceVertexCData',ColorFaces,'FaceColor','flat',...
          'FaceAlpha',0.8);
else
    patch('Faces',Faces','Vertices',Vertices,...
          'FaceVertexCData',ColorFaces);
end
view(3);
axis tight;
if nargin == 5
    axis(Faxis)
end
axis off;
if isempty(strfind(Filename_res,'.png'))
    Filename_res = [Filename_res '.png'];
end
saveas(gcf,Filename_res,'png');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function WriteKochCube(Res,n,KColorMap,Obj_filename)
if nargin==1
    Obj_filename = 'KochCube.obj';
end

Res                          = sortrows(Res,8);
[VertTotal,Faces,ColorFaces] = Generate3DKochCube(Res,n,KColorMap);

if isempty(strfind(Obj_filename,'.obj'))
    Mtl_filename = [Obj_filename '_material.obj'];
    Obj_filename = [Obj_filename '.obj'];
else
    Mtl_filename = [Obj_filename(1:end-4) '_material.obj'];
end
nColor = max(Res(:,8));
IColor1 = zeros(1,nColor);
IColor2 = zeros(1,nColor);
for i=1:nColor
    IColor1(i) =find(Res(:,8)==i,1,'first');
    IColor2(i) =find(Res(:,8)==i,1,'last');
end
fid = fopen(Mtl_filename, 'wt');
fprintf(fid, '# %s\n',Mtl_filename);
fprintf(fid, '#\n');
fprintf(fid, '# File generated with %s\n',mfilename);
fprintf(fid, '# Date : %s\n',datestr(date, 26));
for i=1:nColor
    fprintf(fid, 'newmtl C%03d\n',i);
    fprintf(fid, 'Ka %f %f %f\n',ColorFaces(IColor1(i),:));
    fprintf(fid, 'Kd %f %f %f\n',ColorFaces(IColor1(i),:));
    fprintf(fid, 'd 0.2000\n');
    fprintf(fid, 'illum 1\n');
end
fclose(fid);

fid = fopen(Obj_filename, 'wt');
fprintf(fid, '# %s\n',Obj_filename);
fprintf(fid, '#\n');
fprintf(fid, '# File generated with %s\n',mfilename);
fprintf(fid, '# Date   : %s\n',datestr(date, 26));
fprintf(fid, '# Author : Guillaume Jacquenot\n');
fprintf(fid, '# Koch_cube\n');
fprintf(fid, 'mtllib %s\n',Mtl_filename);
fprintf(fid, 'v %03d %03d %03d\n', VertTotal');
% Normal of each faces has one of the following values
fprintf(fid, 'vn  1.0  0.0  0.0\n');
fprintf(fid, 'vn  0.0  1.0  0.0\n');
fprintf(fid, 'vn  0.0  0.0  1.0\n');
fprintf(fid, 'vn -1.0  0.0  0.0\n');
fprintf(fid, 'vn  0.0 -1.0  0.0\n');
fprintf(fid, 'vn  0.0  0.0 -1.0\n');

% Evaluate the normal of each
NormTotal    = Res(:,5:7);
[X,Y]        = find(NormTotal);
Norm         = Y.*sign(sum(NormTotal,2));
Norm(Norm<0) = 3-Norm(Norm<0);
Norm         = repmat(Norm,1,4)';

% Merge node and vertices
nRes           = size(Res,1);
Total          = zeros(8,nRes);
Total(1:2:7,:) = Faces;
Total(2:2:8,:) = Norm;

% fprintf(fid, 'f  %03d//%1d  %03d//%1d  %03d//%1d\n',Total(1:6,:));
% fprintf(fid, 'f  %03d//%1d  %03d//%1d  %03d//%1d\n',Total([1:2 5:8],:));

for i=1:nColor
    fprintf(fid, 'g G%03d\n',i);
    fprintf(fid, 'usemtl C%03d\n',i);
    fprintf(fid, 'f  %03d//%1d  %03d//%1d  %03d//%1d\n',...
                                         Total(1:6,IColor1(i):IColor2(i)));
    fprintf(fid, 'f  %03d//%1d  %03d//%1d  %03d//%1d\n',...
                                   Total([1:2 5:8],IColor1(i):IColor2(i)));
end
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Vertices,Faces,ColorFaces] = Generate3DKochCube(Res,n,KColorMap)
if nargin==1
    n         = max(Res(:,8));
    KColorMap = 'summer';
end
nRes     = size(Res,1);
Vertices = zeros(4*nRes,3);
for i=1:nRes
    verts = repmat(Res(i,1:3),4,1) + ...
        Res(i,4)/2 * Generate_pattern_matrix(Res(i,5:7),false);
    Vertices(4*(i-1)+1:4*i,:) = verts;
end
Faces = reshape(1:4*nRes,4,nRes);

ColorM = zeros(2*n+1,3);
eval(['ColorM = flipud(' KColorMap '(2*n+1));']);
ColorFaces = ColorM(Res(:,8),:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  res = iscolormap(cmap)
% This function returns true if 'cmap' is a valid colormap
LCmap = {...
    'autumn'
    'bone'
    'colorcube'
    'cool'
    'copper'
    'flag'
    'gray'
    'hot'
    'hsv'
    'jet'
    'lines'
    'pink'
    'prism'
    'spring'
    'summer'
    'white'
    'winter'
};

res = ~isempty(strmatch(cmap,LCmap,'exact'));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%