Sierpinski's thing
Anolethron
an algorithm that can be generalized to the construction of a Menger Sponge.
e1 = {1, 0}; e2 = {0, 1}; p1 = {0, 0}; p2 = {1, 0}; p3 = {1, 1}; p4 = {0,
1};
Sierpinski[{p1_, p2_, p3_, p4_}] :=
Delete[Flatten[
Table[{p1 + m e1 + n e2, p2 + m e1 + n e2, p3 + n e2 + m e1,
p4 + m e1 + n e2}, {n, 0, 2}, {m, 0, 2}], 1], 5];
Sierpinski1 = Sierpinski[{p1, p2, p3, p4}]
Sierpinski2[ls_] := Flatten[Map[Sierpinski, ls], 1]
S2 = Sierpinski2[Sierpinski1]
Sierpinski3[n_] := Nest[Sierpinski2, {{p1, p2, p3, p4}}, n]
Sierpinski3[3]
Now, I'm not good enough to think of a much more complicated construction
and the problem is that with this algorithm the lengths of the squares I
construct at each step does not scale down with the level of the carpet I'm
constructing: e.g. He builds 9 squares from the big one at the beginning and
deletes the central one, it's ok. But as I Iterate the process at each
smaller square It builds squares of the same size, so what I get is just a
big black figure. It obviously does this way because in the algorithm
there's no instruction to decrease the size of the base vectors (e1,e2).
Thing is I can't think of a way to give mathematica that instruction inside
the Nest or in the definition of the basic "Sierpinski" function. I need
some help. Thanks in advance. This is the expected result:
http://mathworld.wolfram.com/SierpinskiCarpet.html
Szabolcs
I don't understand completely what you were doing above, but here's a
function that constructs a Sierpinski carpet:
In[1]:= pieces = Complement[
Join@@Table[{i, j}, {i, 0, 2}, {j, 0, 2}],
{{1, 1}}]
Out[1]= {{0,0},{0,1},{0,2},{1,0},{1,2},{2,0},{2,1},{2,2}}
In[2]:= sierp[cornerPt_, sideLen_, n_] :=
sierp[cornerPt + #1*(sideLen/3), sideLen/3, n-1] & /@ pieces
In[3]:= sierp[cornerPt_, sideLen_, 0] :=
Rectangle[cornerPt, cornerPt + sideLen*{1, 1}]
In[4]:= Graphics[sierp[{0, 0}, 1, 5], AspectRatio -> Automatic]//Show
Could someone please explain why is this SO MUCH slower in Mathematica 6
than in Mathematica 5.2?
Szabolcs
Jean-Marc Gulliet
Hi,
If you have version 6.0, the following code will draw a nice Sierpinski
Carpet.
rules = {0 -> {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}},
1 -> {{1, 1, 1}, {1, 0, 1}, {1, 1, 1}}};
f[m_: 1] := ArrayFlatten[m /. rules]
drawSerp[n_] := MatrixPlot[Nest[f, 1, n], FrameTicks -> None]
drawSerp[5]
You can see the resulting picture at
http://homepages.nyu.edu/~jmg336/Sierpinski%20Carpet.png
If you do not have access to version 6.0 and are interested by this
approach, let me know so I can tell you how to tweak the code for
version 5.2.
Regards,
Jean-Marc
Jean-Marc Gulliet
For what is worth, I have made some speed tests of your code on version
5.2 and 6.0. The speed difference is striking in this case. The culprit
seems to be the rendering engine in v6. As a summary,
------------------------------------------------------
V5.2 ! V6.0
! Native Engine ! Legacy Engine 5.2*
------------------------------------------------------
Create ! Display ! Create ! Display ! Create ! Display
0.938 ! 0.187 ! 1.703 ! 24.765 ! 0.954 ! 0.234
------------------------------------------------------
Times expressed in second.
Columns "Create" show cpu time.
Columns "Display" show elapsed time.
* See Dimitris's post [1].
V6 appears to be slower than V5.2 even when using the compatibility mode
that emulate (?) the graphic engine of V5.2.
(* Code tested with versions 5.2 and 6.0 *)
In[1]:=
pieces = Complement[Join @@ Table[{i, j}, {i, 0, 2}, {j, 0, 2}], {{1, 1}}];
sierp[cornerPt_, sideLen_, n_] :=
(sierp[cornerPt + #1*(sideLen/3), sideLen/3, n - 1] & ) /@ pieces
sierp[cornerPt_, sideLen_, 0] := Rectangle[cornerPt,
cornerPt + sideLen*{1, 1}]
(* Version 5.2 *)
In[4]:=
Timing[g = Graphics[sierp[{0, 0}, 1, 5], AspectRatio -> Automatic]][[1]]
Out[4]=
0.938 Second
In[5]:=
start = TimeUsed[];
Show[g];
stop = TimeUsed[];
stop - start
[graphic deleted]
Out[8]=
0.187
(* Version 6.0 *)
In[1]:= pieces =
Complement[Join @@ Table[{i, j}, {i, 0, 2}, {j, 0, 2}],
{{1, 1}}];
sierp[cornerPt_, sideLen_, n_] :=
(sierp[cornerPt + #1*(sideLen/3), sideLen/3, n - 1] & ) /@ pieces
sierp[cornerPt_, sideLen_, 0] := Rectangle[cornerPt,
cornerPt + sideLen*{1, 1}]
In[4]:= Timing[
g = Graphics[sierp[{0, 0}, 1, 5], AspectRatio -> Automatic]; ][[
1]]
Out[4]= 1.703
In[5]:= start = TimeUsed[];
g
stop = TimeUsed[];
stop - start
[graphic deleted]
Out[8]= 24.765
(* Switching to legacy graphic engine. See Dimitris's post [1] *)
In[9]:= << Version5`Graphics`
In[10]:= Timing[
g = Graphics[sierp[{0, 0}, 1, 5], AspectRatio -> Automatic]; ][[
1]]
Out[10]= 0.954
In[11]:= start = TimeUsed[];
Show[g];
stop = TimeUsed[];
stop - start
[graphic deleted]
Out[14]= 0.234
In[15]:= << Version6`Graphics`
Regards,
Jean-Marc
[1] "general"
http://groups.google.com/group/comp.soft-sys.math.mathematica/browse_thread/thread/eb2b72709109f82b/?hl=en#
Szabolcs Horvát
> For what is worth, I have made some speed tests of your code on version
> 5.2 and 6.0. The speed difference is striking in this case. The culprit
> seems to be the rendering engine in v6. As a summary,
>
> ------------------------------------------------------
> V5.2 ! V6.0
> ! Native Engine ! Legacy Engine 5.2*
> ------------------------------------------------------
> Create ! Display ! Create ! Display ! Create ! Display
> 0.938 ! 0.187 ! 1.703 ! 24.765 ! 0.954 ! 0.234
> ------------------------------------------------------
It seems that it is not even the rendering itself that takes so long.
With anti-aliasing turned off, the performance does not improve
noticeably. But when the graphics have already been displayed,
resizing is quite fast.
The graphics are probably converted to some other representation
before they are displayed, and the conversion takes a long time.
I experimented a bit more with this:
sierp[ ] is the original function, display time is: -- 21.2 sec
sierp3[ ] is the same as sierp[ ], but produces a flat list. -- 18.8 sec
sierp2[ ] uses Translate[ ] to shift the pieces. -- 75 sec.
sierp4[ ] uses trans[ ] (defined below) to simulate the effect of
Translate[ ]. Creation time is 2.4 sec, so there is no excuse for
sierp2[ ]'s output to be rendered so slowly. It's probably better to
avoid Translate[ ].
I hope that the design of the rendering system allows for some
optimisation here and we shall see a huge improvement in the next
version. This is an unacceptable performance regression from 5.2.
In[1]:= pieces=Complement[Join@@Table[{i,j},{i,0,2},{j,0,2}],{{1,1}}];
sierp[cornerPt_,sideLen_,n_]:=(sierp[cornerPt+#1*(sideLen/3),sideLen/3,n-1]&)/@pieces
sierp[cornerPt_,sideLen_,0]:=Rectangle[cornerPt,cornerPt+sideLen*{1,1}]
In[4]:= gr=Graphics[sierp[{0,0},1,5]];
In[5]:= start=TimeUsed[];
gr
stop=TimeUsed[];
stop-start
Out[8]= 21.25
Here we turn off the anti-aliasing:
In[9]:= start=TimeUsed[];
Style[gr,Antialiasing->False]
stop=TimeUsed[];
stop-start
Out[12]= 21.203
This version producees a flat list of graphics primitives. It seems
that the processing takes a bit less time in this case.
In[13]:= sierp3[cornerPt_,sideLen_,n_]:=Join@@((sierp3[cornerPt+#1*(sideLen/3),sideLen/3,n-1]&)/@pieces)
sierp3[cornerPt_,sideLen_,0]:={Rectangle[cornerPt,cornerPt+sideLen*{1,1}]}
In[15]:= gr3=Graphics[sierp3[{0,0},1,5]];
In[16]:= start=TimeUsed[];
gr3
stop=TimeUsed[];
stop-start
Out[19]= 18.797
With Translate[ ], it is awfully slow ...
In[20]:= sierp2[sideLen_,n_]:=Translate[sierp2[sideLen/3,n-1],sideLen/3#]&/@pieces
sierp2[sideLen_,0]:=Rectangle[{0,0},{sideLen,sideLen}]
In[22]:= gr2=Graphics[sierp2[1,5]];
In[23]:= start=TimeUsed[];
gr2
stop=TimeUsed[];
stop-start
Out[26]= 74.781
Now let's use our own version of Translate[], to see if it is really
necessary for Translate[] to be so slow!
In[27]:= sierp4[sideLen_,n_]:=trans[sierp4[sideLen/3,n-1],sideLen/3#]&/@pieces
sierp4[sideLen_,0]:=Rectangle[{0,0},{sideLen,sideLen}]
In[29]:= trans[elem_,vec_]:=elem/.Rectangle[a_,b_]:>Rectangle[a+vec,b+vec]
In[30]:= Timing[gr4=Graphics[sierp4[1,5]];]
Out[9]= {2.359,Null} (* creation time *)
Szabolcs
Anolethron
Szabolcs
> But how do you generalize it to a menger sponge?
>
>
>
??? That's very starightforward
In[3]:= pieces =
Complement[
Flatten[Table[{i, j, k}, {i, 0, 2}, {j, 0, 2}, {k, 0, 2}],
2], {{1, 1, 1}, {0, 1, 1}, {2, 1, 1}, {1, 0, 1}, {1, 2, 1}, {1, 1,
0}, {1, 1, 2},}]
Out[3]= {{0, 0, 0}, {0, 0, 1}, {0, 0, 2}, {0, 1, 0}, {0, 1, 2}, {0, 2,
0}, {0, 2, 1}, {0, 2, 2}, {1, 0, 0}, {1, 0, 2}, {1, 2, 0}, {1, 2,
2}, {2, 0, 0}, {2, 0, 1}, {2, 0, 2}, {2, 1, 0}, {2, 1, 2}, {2, 2,
0}, {2, 2, 1}, {2, 2, 2}}
In[4]:= menger[cornerPt_, sideLen_, n_] :=
menger[cornerPt + #1*(sideLen/3), sideLen/3, n - 1] & /@ pieces
In[5]:= menger[cornerPt_, sideLen_, 0] :=
Cuboid[cornerPt, cornerPt + sideLen*{1, 1, 1}]
In[9]:= Graphics3D[menger[{0, 0, 0}, 1, 3]]