Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to draw a Sierpinski_carpet. I have code that works, but I think there is a more elegant way to do than my way. Maybe I couls use Tuples or Permutations or some similar function to simplify my code.

enter image description here

f[{{x1_, y1_}, {x2_, y2_}}] := Map[Mean, {
    {{{x1, x1, x1}, {y1, y1, y1}}, {{x1, x1, x2}, {y1, y1, y2}}},
    {{{x1, x1, x1}, {y1, y1, y2}}, {{x1, x1, x2}, {y1, y2, y2}}},
    {{{x1, x1, x1}, {y1, y2, y2}}, {{x1, x1, x2}, {y2, y2, y2}}},
    {{{x1, x1, x2}, {y1, y1, y1}}, {{x1, x2, x2}, {y1, y1, y2}}},
    {{{x1, x1, x2}, {y1, y2, y2}}, {{x1, x2, x2}, {y2, y2, y2}}},
    {{{x1, x2, x2}, {y1, y1, y1}}, {{x2, x2, x2}, {y1, y1, y2}}},
    {{{x1, x2, x2}, {y1, y1, y2}}, {{x2, x2, x2}, {y1, y2, y2}}},
    {{{x1, x2, x2}, {y1, y2, y2}}, {{x2, x2, x2}, {y2, y2, y2}}}
    }, {3}];
d = Nest[Join @@ f /@ # &, {{{0., 0.}, {1, 1}}}, 3];
Graphics[Rectangle @@@ d]
Clear["`*"]
share|improve this question

9 Answers 9

up vote 52 down vote accepted

This seems the most natural to me:

Mathematica graphics

carpet[n_] := Nest[ArrayFlatten[{{#, #, #}, {#, 0, #}, {#, #, #}}] &, 1, n]

ArrayPlot[carpet @ 5, PixelConstrained -> 1]

Mathematica graphics

Shorter (in InputForm), but perhaps harder to read and slightly slower, though speed hardly matters given the geometric memory usage:

carpet[n_] := Nest[ArrayFlatten @ ArrayPad[{{0}}, 1, {{#}}] &, 1, n]

Style by level

With a minor change we can increment the values with each fractal level allowing identification such as styling, or other processing.

Wild colors are but a few commands away:

carpet2[n_] := Nest[ArrayFlatten[{{#, #, #}, {#, 0, #}, {#, #, #}}] &[1 + #] &, 1, n]

Table[
  ArrayPlot[carpet2 @ 4, PixelConstrained -> 1, ColorFunction -> color],
  {color, ColorData["Gradients"]}
]

Mathematica graphics


Extension to three dimensions

A Menger sponge courtesy of chyanog, with refinements:

carpet3D[n_] :=
   With[{m = # (1 - CrossMatrix[{1,1,1}])}, Nest[ArrayFlatten[m, 3] &, 1, n]]

Image3D[carpet3D[4], ColorFunction -> GrayLevel]

Mathematica graphics


Element coordinates

If you wish to get coordinates for display with graphics primitives or analysis this can be done efficiently using SparseArray Properties:

coords = SparseArray[#]["NonzeroPositions"] &;

Example usages:

Graphics @ Point @ coords @ carpet @ 4

Mathematica graphics

Graphics3D[Cuboid /@ coords @ carpet3D @ 3]

Mathematica graphics

share|improve this answer
    
As a newbie MMA user, I agree this seems most natural. +1 –  RunnyKine Mar 25 '13 at 19:12
    
@RunnyKine Thanks! :-) –  Mr.Wizard Mar 25 '13 at 19:12
2  
The quilt is pretty. –  whuber Mar 25 '13 at 19:22
    
+1 for having taught me PixelConstrained -> 1 and for the beautiful images. –  Federico Mar 25 '13 at 22:35
3  
Very nice! 3D version,carpet[n_] := Nest[ArrayFlatten[{{{#, #, #}, {#, 0, #}, {#, #, #}}, {{#, 0, #}, {0, 0, 0}, {#, 0, #}}, {{#, #, #}, {#, 0, #}, {#, #, #}}}, 3] &, 1, n]; Image3D[carpet[3]] –  chyaong Mar 26 '13 at 3:35

Assuming you want a vector-based image, it's more efficient to cut holes:

translations = {#, #} & /@ Complement[Tuples[{-1, 0, 1}, 2], {{0, 0}}];
shrink[{{x0_, y0_}, {x1_, y1_}}] := {{2 x0 + x1, 2 y0 + y1}, {x0 + 2 x1, y0 + 2 y1}}/3
children[sq : {{x0_, y0_}, {x1_, y1_}}] := 
  With[{side = x1 - x0, newsq = shrink[sq]},
    (newsq + #) & /@ (side translations)]
gen = NestList[Join @@ children /@ # &, {{{1, 1}, {2, 2}}/3}, 3];
Graphics[Rectangle @@@ Join @@ gen]

For the same generation you need to draw $1/8$ of the rectangles.

ngen = 4;
gen = NestList[Join @@ children /@ # &, {{{1, 1}, {2, 2}}/3}, ngen];
colors = Table[Blend[{RGBColor[0.5, 0.5, 1], White}, i/ngen], {i, 0, ngen}];
Graphics[{Black, Rectangle[{0, 0}, {1, 1}], White, Rectangle @@@ Join @@ gen}]
Graphics[{Red, Rectangle[{0, 0}, {1, 1}], MapThread[Prepend, {Apply[Rectangle, gen, {2}], colors}]}]

Sierpinski 2D

Here is the 3D version:

rule = 0 -> CrossMatrix[{1, 1, 1}];
Graphics3D[Cuboid /@ Position[Nest[ArrayFlatten[# /. rule, 3] &, 0, 3], 0]]

Sierpinski 3D

share|improve this answer

Here are two methods using rules, shamelessly modified from a MathGroup posting (http://forums.wolfram.com/mathgroup/archive/2007/May/msg01356.html).

rules = # -> ArrayPad[{{0}}, 1, #] & /@ {0, 1}
f1[m_] := ArrayFlatten[m /. rules]
drawSerp[n_] := ArrayPlot[Nest[f1, 1, n], Frame -> False]
drawSerp[3]

Mathematica graphics

An alternative cute ASCII implementation, borrowed from http://rosettacode.org/wiki/Sierpinski_carpet#Mathematica

n = 3;
Grid[Nest[ArrayFlatten[# /. rules] &, {{1}}, 
   n] //. {0 -> " ", 1 -> "#"}]

Mathematica graphics

share|improve this answer
    
Shamelessly upvoting. –  rcollyer Mar 25 '13 at 16:40
1  
+1 Using ArrayPlot instead of ColorNegate@Image in the first solution gives a crisper rendering (and happens to be about 50% faster, too). –  whuber Mar 25 '13 at 17:40

After an initial attempt with a Graphics-based solution, it became apparent that raster-based solutions would be far more efficient. Methods based on ArrayPlot work nicely, but I wondered whether image-based manipulations might be the most efficient possible way, given that they would be optimized for precisely the kinds of operations being performed here.

Indeed, the following is an order of magnitude faster than anything I have timed yet, while sharing the expressive clarity of several other answers that have already appeared. Another advantage is that it scales the output resolution to match the depth of the approximation to the carpet.

carpet[n_, white_: 1, black_: 0] := 
 Nest[{ImageAssemble[{{#1, #1, #1}, {#1, #2, #1}, {#1, #1, #1}}], 
       ImageResize[#2, 3 First[ImageDimensions[#2]]]} & @@ # &, 
   Image /@ {{{black}}, {{white}}}, n] // First

It literally pieces all the pixels together, starting with a black pixel (Image[{{black}}]) and a white pixel (Image[{{white}}])--whose colors you may optionally specify as arguments--reassembling them at each stage in the familiar three by three pattern (ImageAssemble) and, preparatory to the next stage, rescaling the central white pixel to match the size (ImageResize). (At the end it throws away the upscaled white image.) Here is carpet[7], a $2187$ by $2187$ image ($0.05$ seconds):

Carpet

It is an easy exercise to modify this to start with any central image (the "focus") instead of just a white pixel. Under prompting by Mr.Wizard (see comments), I offer the sharpest possible solution. To create it, you need to begin with an image whose dimensions are a power of three and downsize it all the way to one pixel, creating a list of images that will serve as the foci of interest:

i = ExampleData[{"TestImage", "Lena"}];             (* Original image *)
n = 3^(k = Floor[Log[3, Min[ImageDimensions[i]]]]); (* Nearest lower power of 3 *)
focus = ImageResize[ImageCrop[i, {n, n}] // ImageAdjust, n/3^#] & /@ Range[k, 0, -1]

An attractive Sierpinski carpet is now particularly simple to make. Here is a general implementation with a default black background

carpet[focus_List, background_: Image[{{0}}]] := 
  Fold[ImageAssemble[{{#1, #1, #1}, {#1, #2, #1}, {#1, #1, #1}}] &, background, focus];

and here is an application to the test image:

carpet[focus]

Lena's carpet

Because this process is so fast, in less than one second we can make carpets of all the example data (assuming they have already been downloaded):

carpet[ExampleData[#], Image[{{1}}]] & /@ ExampleData["TestImage"]

Quilt of carpets

share|improve this answer
    
What's the code to make either of the images in your output? (I think I could figure it out but I'm suggesting that you include it.) –  Mr.Wizard Mar 25 '13 at 19:14
    
@Mr.Wizard I gave the code for the first one :-). –  whuber Mar 25 '13 at 19:17
    
Thanks. v7 users FYI: you'll need black = ColorConvert[black, "RGB"] –  Mr.Wizard Mar 25 '13 at 19:26
    
I observe that this isn't producing pixel-accurate images; perhaps you could calculate the image size necessary and use a different resize algorithm to avoid the blurring. –  Mr.Wizard Mar 25 '13 at 19:28
1  
@Mr.Wizard To achieve a really sharp solution, we need to crop the central image to dimensions which are powers of three. I have modified my answer to show this technique. Thank you for your suggestions! –  whuber Mar 26 '13 at 14:17

One lame method would be to create a replacement rule which replaces one Rectangle's of an graphics which the appropriate 8 others.

f[p_, {min_, max_}] := p/3 max + (1 - p/3) min;
rule = Rectangle[{xmin_, ymin_}, {xmax_, ymax_}] -> 
   With[{expr = Table[If[i =!= 1 || j =!= 1, 
     Rectangle[{f[i, xmin, xmax], f[j, ymin, ymax]}, 
               {f[i + 1, xmin, xmax], f[j + 1, ymin, ymax]}], {}], 
   {j, 0, 2}, {i, 0, 2}]}, Flatten[(expr &)[min, max]]];

Then a simple Nest with one initial rectangle does the job.

Graphics[Nest[# /. rule &, Rectangle[{-1., -1.}, {1., 1.}], 4]]

I spare another black and white image here, because a more beautiful thing can be created when you use the iterative style function which can be found on the Wikipedia page, compile it and save the number of iterations

fillCarpet = Compile[{{pixel, _Integer, 1}},
  Module[{x = pixel[[1]], y = pixel[[2]], result = 1, iter = 0},
   While[x > 0 || y > 0,
    If[Mod[x, 3] === 1 && Mod[y, 3] === 1, result = 0; Break[]];
    x = Quotient[x, 3];
    y = Quotient[y, 3];
    ++iter;
    ];
   {result, iter}
   ], CompilationTarget -> "C", Parallelization -> True, 
  RuntimeAttributes -> {Listable}]

Image@fillCarpet[Table[{i, j}, {j, 0, 3^6-1}, {i, 0, 3^6-1}]]

Here the self-similarity is highlighted through the iteration count. Now we can create various images which have sizes of 3^n-1 for different n

Show[Image@
    fillCarpet[Table[{i, j}, {j, 0, 3^# - 1}, {i, 0, 3^# - 1}]], 
   ImageSize -> 256] & /@ Range[3, 5]

Mathematica graphics Mathematica graphics Mathematica graphics

share|improve this answer

Many of the approaches here use image processing functions and they are blazing fast and very cool. However, there are advantages of a primitives based approach. When studying fractals, sometimes you need vertex information of the approximations, for example. Also, I don't think these image based techniques extend easily to self-similar sets that are not based on a rectangular decomposition. So, here's a reasonably fast approach using graphics primitives.

Clear[step];
step = With[{shifts = N@{
       {0, 0}, {1/3, 0}, {2/3, 0},
       {0, 1/3},                {2/3, 1/3},
       {0, 2/3}, {1/3, 2/3}, {2/3, 2/3}
       }},
   Compile[{{vertices, _Real, 2}},
    Module[{vv},
     vv = vertices/3;
     Table[pt + # & /@ vv, {pt, shifts}]
     ], CompilationTarget -> "C", RuntimeAttributes -> {Listable},
    Parallelization -> True
    ]
   ];

depth = 7;
init = N[{
    {{0, 0}, {1, 0}, {1, 1}, {0, 1}}
    }];
t = AbsoluteTime[];
polygons = Flatten[Nest[step, init, depth], depth]; // AbsoluteTiming
Graphics[Polygon[polygons]]

(* Out: {0.746998, Null} *)

(* Repetitive Sierpinski picture omitted.

(* Put in a new cell and execute all at once for reliable timing. *)
AbsoluteTime[] - t

(* Out: 4.559120 *)

About 10 times slower than whuber's carpet command. Of course, this can never be as fast as the image based stuff, which is manipulated at a lower level.

share|improve this answer
    
+1 Once you head down this road to generalization, you might as well implement an IFS. –  whuber Mar 26 '13 at 18:42
1  
@whuber Almost - have you seen the packages on my Mathematica page: facstaff.unca.edu/mcmcclur/Mathematica.html These include not just image generation code but code to compute Laplacians on self-similar sets and code to compute the dimensions of self-similar sets of finite type, where the open set condition is not satisfied. –  Mark McClure Mar 26 '13 at 18:52
    
@whuber or an L-system. I considered adapting my bf parser to do just that. –  rcollyer Mar 27 '13 at 12:17

Thanks for @Mark McClure. Inspired by him, my original code is simplified.

This seems also natural.

enter image description here

f[v_] := Table[i + j, {i, Drop[Tuples[{0, 1, 2}, 2], {5}]}, {j, v}]/3.;
d = Nest[Join @@ f /@ # &, N@{{{0, 0}, {1, 1}}}, 3];
Graphics[Rectangle @@@ d]

enter image description here

It's also easily generalized to 3D:

f[v_] := Table[ i + j, {i, Select[Tuples[{0, 1, 2}, 3], Count[#, 1] < 2 &]},
  {j, v}]/3.;
d = Nest[Join @@ f /@ # &, N@{{{0, 0, 0}, {1, 1, 1}}}, 3];
Graphics3D[{EdgeForm[], Cuboid /@ d}, Boxed -> False]

enter image description here

share|improve this answer
    
Please see the most recent section of my answer. I believe if you run Timings you will find it is superior. –  Mr.Wizard Mar 27 '13 at 18:32
    
@Mr.Wizard I'm impressed. –  chyaong Mar 28 '13 at 5:03

I might as well... as a variation, here's a chaos game method for generating the carpet:

With[{verts = DeleteCases[Tuples[{-1, 0, 1}, {2}], {0, 0}], n = 1*^6},
   Graphics[{AbsolutePointSize[1/2],
             Point[NestList[(2 RandomChoice[verts] + #)/3 &, RandomReal[{-1, 1}, 2], n]]}]]

chaos game carpet

share|improve this answer
    
Now, the carpet is drawn...and it would be also good for us (mathematicians, or physicists) to obtain the corresponding adj.matrix. However, this can be possible if one makes the program think that the drawn object is a graph. Is there any intelligent way to do this? Sorry, if I am posting too stupid comment...Thank you very much in advance! –  Lady InRed Aug 19 '13 at 9:56

This is more or less equivalent to halirutan's approach, but slightly compacted (adapted from code I wrote ~ 10 years ago, before Tuples[] came along):

Block[{n = 5, pos = Select[Tuples[2 {-1, 0, 1}/3, {2}], (Count[#, 0] < 2) &]},
      Graphics[Nest[Function[g, (g /. v_ /; VectorQ[v, NumericQ] :> v/3 - #) & /@ pos],
                    Rectangle[{-1, -1}, {1, 1}], n]]]

the carpet

The sponge version:

Block[{n = 4, pos = Select[Tuples[2 {-1, 0, 1}/3, {3}], (Count[#, 0] < 2) &]}, 
    Graphics3D[{Directive[EdgeForm[], GrayLevel[1/5], Glow[Gray], Specularity[1/2, 100]], 
                Nest[Function[g3d, (g3d /. v_ /; VectorQ[v, NumericQ] :> v/3 - #) & /@ pos],
                     Cuboid[{-1, -1, -1}, {1, 1, 1}], n]}, Boxed -> False]]

the sponge

A sponge with "less faded" coloring:

Block[{n = 4, pos = Select[Tuples[2 {-1, 0, 1}/3, {3}], (Count[#, 0] < 2) &]}, 
   Graphics3D[{Directive[EdgeForm[], ColorData["Legacy", "DodgerBlue"], 
                         Specularity[3/4, 20]], 
               Nest[Function[g3d, (g3d /. v_ /; VectorQ[v, NumericQ] :> v/3 - #) & /@ pos], 
                    Cuboid[{-1, -1, -1}, {1, 1, 1}], n]},
              Boxed -> False, Lighting -> "Neutral"]]

blue sponge

share|improve this answer
    
In your picture the sponge looks faded, and on my system it renders as a solid gray cube. –  Mr.Wizard Apr 10 '13 at 8:55
    
The faded-ness was more or less deliberate; the coloring directives can be changed if needed. I'll render something differently colored if needed. The "solid gray" was unintentional; it seems I copied a wrong cell from my old notebook. See it now. –  J. M. Apr 10 '13 at 9:50
    
Okay. Both functions (yours and mine) take a long time render at n = 4 but this 3D code takes 1.4 seconds to merely evaluate while mine takes 0.05 second. Likewise your 2D code takes 0.2558 second whereas mine takes 0.0005 second. In what situation will these methods be preferred? –  Mr.Wizard Apr 10 '13 at 10:07
    
@Mr. Wizard, well, it's old code. :) Except for the Tuples[], and Glow[]/Specularity[], everything there was more or less how it was when I wrote it in version 4 (e.g. the original code used SurfaceColor[], and I had rolled my own version of Tuples[] before). I also figured that seeing another replacement-rule based implementation might prove instructive. –  J. M. Apr 10 '13 at 10:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.