function hs=fractal(N,axiom,gen,cns,scal)
% Function to generate fractal geometries.
% Works in two ways: 1) Generation of one of the three
% predefined geometries or 2) general use, where axiom
% and generator must be specified.
%
% Pre-defined use of the function:
% -------------------------------------------------------------------------
% N = the level to be generated.
% axiom = geometry type:
% 1: Sierpinski's pyramide
% 2: Sierpinski's triangle
% 3: Menger's sponge
% EXAMPLES:
% fractal(2,1); % show the Sierpinski's triangle with 2 levels.
% fractal(0,3); % show the axiom of the Menger's sponge.
% fractal(1,3); % show the generator of the Menger's sponge.
%
% The remaining arguments are only used for the general use of the
% function.
%
% General use of the function: NOT WELL TESTED !!!
% -------------------------------------------------------------------------
% N = the level to be generated.
% axiom = is the coordinates of the faces (NOT just the corner coordinates).
% generator = is the coordinates of the pattern faces.
% cns = is the number of coordinates in the axiom that belongs to each face.
% (The length of cns must equal the number of faces).
% scal = the scaling factor.
% hs = handles to surfaces.
% -------------------------------------------------------------------------
% Copyright (C) Jakob L. Laugesen, laugesen@fysik.dtu.dk / jakob@sensory.dk,
% Chaos Group and BioSim, Dep. of Physics, DTU, Denmark.
% Created: 16/11/09
% Version 1
% -------------------------------------------------------------------------
if nargin==2
if (axiom==1 || axiom==2) && N>6
disp('At most six levels for Sierpinski''s triangle and pyramide.');
return;
end
if (axiom==3) && N>3
disp('At most three levels for Menger''s sponge');
return;
end
end
axiomc=[0,0,0];hs=[];N=N-1;
if nargin==2
if axiom==1 % pyramide
scal=2;
h=sqrt(3)/2;
axiom=[-0.5 -0.5 -h/2;...
0.5 -0.5 -h/2;...
0.5 0.5 -h/2;...
-0.5 0.5 -h/2;...
-0.5 -0.5 -h/2;...
0.5 -0.5 -h/2;...
0 0 h/2;...
0.5 -0.5 -h/2;...
0.5 0.5 -h/2;...
0 0 h/2;...
0.5 0.5 -h/2;...
-0.5 0.5 -h/2;...
0 0 h/2;...
-0.5 -0.5 -h/2;...
-0.5 0.5 -h/2;...
0 0 h/2];
r=size(axiom,1);
gen=[axiom/2+repmat([-0.25 -0.25 -h/4],r,1); axiom/2+repmat([-0.25 0.25 -h/4],r,1);...
axiom/2+repmat([0.25 -0.25 -h/4],r,1); axiom/2+repmat([0.25 0.25 -h/4],r,1);...
axiom/2+repmat([0 0 h/4],r,1)];
cns=[4 3 3 3 3];
elseif axiom==2 % triangle
scal=2;
h1=1/(2*cos(pi/6));
h2=sqrt(3)/2-h1;
h3=sqrt(1-h1^2);
axiom=[-0.5 -h2 -h3/2;...
0.5 -h2 -h3/2;...
0.5 -h2 -h3/2;...
-0.5 -h2 -h3/2;...
0.5 -h2 -h3/2;...
0.0 0.0 h3/2;...
-0.5 -h2 -h3/2;...
0.0 h1 -h3/2;...
0.0 0.0 h3/2;...
0.5 -h2 -h3/2;...
0.0 h1 -h3/2;...
0.0 0.0 h3/2];
gen=[axiom/2+repmat([-0.25 -h2/2 -h3/4],12,1); axiom/2+repmat([0.25 -h2/2 -h3/4],12,1);...
axiom/2+repmat([0 h1/2 -h3/4],12,1); axiom/2+repmat([0 0 h3/4],12,1)];
cns=[3 3 3 3];
elseif axiom==3 % menger sponge
scal=3;
e=1/scal;
axiom=[ -0.5 -0.5 -0.5;...
0.5 -0.5 -0.5;...
0.5 0.5 -0.5;...
-0.5 0.5 -0.5;...
-0.5 -0.5 0.5;...
0.5 -0.5 0.5;...
0.5 0.5 0.5;...
-0.5 0.5 0.5;...
0.5 -0.5 -0.5;...
0.5 0.5 -0.5;...
0.5 0.5 0.5;...
0.5 -0.5 0.5;...
-0.5 -0.5 -0.5;...
-0.5 0.5 -0.5;...
-0.5 0.5 0.5;...
-0.5 -0.5 0.5;...
0.5 0.5 -0.5;...
0.5 0.5 0.5;...
-0.5 0.5 0.5;...
-0.5 0.5 -0.5;...
0.5 -0.5 -0.5;...
0.5 -0.5 0.5;...
-0.5 -0.5 0.5;...
-0.5 -0.5 -0.5];
D=size(axiom,1);
gen=[axiom/scal+repmat([-e -e -e],D,1); axiom/scal+repmat([0 -e -e],D,1); axiom/scal+repmat([e -e -e],D,1);...
axiom/scal+repmat([-e -e 0],D,1); axiom/scal+repmat([e -e 0],D,1);...
axiom/scal+repmat([-e -e e],D,1); axiom/scal+repmat([0 -e e],D,1); axiom/scal+repmat([e -e e],D,1);...
axiom/scal+repmat([-e 0 -e],D,1); axiom/scal+repmat([e 0 -e],D,1);...
axiom/scal+repmat([-e 0 e],D,1); axiom/scal+repmat([e 0 e],D,1);...
axiom/scal+repmat([-e e -e],D,1); axiom/scal+repmat([0 e -e],D,1); axiom/scal+repmat([e e -e],D,1);...
axiom/scal+repmat([-e e 0],D,1); axiom/scal+repmat([e e 0],D,1);...
axiom/scal+repmat([-e e e],D,1); axiom/scal+repmat([0 e e],D,1); axiom/scal+repmat([e e e],D,1)];
cns=[4 4 4 4 4 4];
end
end
genc=[];
for n=1:size(axiom,1):size(gen,1)
genc=[genc;mean(gen(n:n+size(axiom,1)-1,:))];
end
xyz1=gen;cen1=axiomc;
if N==-1
xyz2=axiom;
elseif N==0
xyz2=gen;
else
for n=1:N
cen2=[];
xyz2=[];
for t=1:size(cen1,1)
centers=repmat(cen1(t,:),size(genc,1),1)+genc/(scal^(n-1));
cen2=[cen2;centers];
xyz2=[];
if n==N
for tt=1:size(cen2,1)
temp=repmat(cen2(tt,:),size(gen,1),1)+gen/(scal^n);
xyz2=[xyz2;temp];
end
end
end
cen1=cen2;
end
end
figure('Renderer','zbuffer');
view(3);
for n=1:size(axiom,1):length(xyz2)
[h]=filltri(xyz2(n:n+size(axiom,1)-1,:),cns);
hs=[hs;h];
try
x=cell2mat(get(hs,'XData')');
y=cell2mat(get(hs,'YData')');
z=cell2mat(get(hs,'ZData')');
ee=unique(([x' y' z']),'rows');
x=ee(:,1:4)';
y=ee(:,5:8)';
z=ee(:,9:12)';
delete(hs(2:end));
set(hs(1),'XData',x);
set(hs(1),'YData',y);
set(hs(1),'ZData',z);
set(hs(1),'LineStyle','none');
hs=hs(1);
end
end
% T=[0.1822 -0.9833 0 0.4005;...
% 0.0343 0.0064 0.9994 -0.5200;...
% 0.9827 0.1821 -0.0349 8.0953;...
% 0 0 0 1.0000];
% T=[ -1.1511 -0.0905 -0.0000 0.6208;...
% 0.0339 -0.4312 0.9272 -0.2649;...
% 0.0840 -1.0673 -0.3746 10.2533;...
% 0 0 0 1.0000];
%
% T=[ 0.3067 -1.0269 0 0.3601;...
% 1.1598 0.2598 0.2079 -0.8138;...
% 0.2465 0.0552 -0.9781 9.8579;...
% 0 0 0 1.0000 ];
% T=[ 0.0416 -0.3407 0 0.1495;...
% 0.0301 0.0023 0.3415 -0.1869;...
% 0.4309 0.0328 -0.0239 3.0308;...
% 0 0 0 0.0010]*1000;
% T=[ -0.3007 0.9537 -0.0000 -0.3265;...
% -0.3261 -0.1028 0.9397 -0.2554;...
% -0.8962 -0.2826 -0.3420 9.4206;...
% 0 0 0 1.0000];
T=[ 1.1273 -0.2500 -0.0000 -0.4387;...
0.0174 0.0786 0.9976 -0.5468;...
0.2494 1.1246 -0.0697 8.9222;...
0 0 0 1.0000]
view(T);
light('Position',[0 0.1 0],'Color',[1 0 0],'Style','infinite');
set(gca,'XColor','white');
set(gca,'YColor','white');
set(gca,'ZColor','white');
set(gca,'XTick',[]);
set(gca,'YTick',[]);
set(gca,'ZTick',[]);
axis equal;
function [hs]=filltri(x,cns)
hold on;
hs=[];
co=1;
for n=1:length(cns)
k=co:co+cns(n)-1;
h=fill3(x(k,1),x(k,2),x(k,3),'cyan');
hs=[hs;h];
co=co+cns(n);
end
