Generating a Sierpinski carpet

Question

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.

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["`*"]

Answer 1

Version 11.1 introduces MengerMesh:

MengerMesh[3]

---

---

This seems the most natural to me:

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

ArrayPlot[carpet @ 5, PixelConstrained -> 1]

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"]}
]

---

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] ]

---

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

Graphics3D[Cuboid /@ coords @ carpet3D @ 3]

Answer 2

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}]}]

Here is the 3D version:

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

Answer 3

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]

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

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

Answer 4

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):

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]

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"]

Answer 5

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]

Answer 6

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.

Answer 7

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

This seems also natural.

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]

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]

Answer 8

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]]}]]

Answer 9

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 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]]

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"]]

Answer 10

Wolfram finally (V10.2) has made a function SubstitutionSystem that does the job. These are examples from documentation (see Applications part).

ArrayPlot[
 SubstitutionSystem[{1 -> {{1, 1, 1}, {1, 0, 1}, {1, 1, 1}}, 
 0 -> ConstantArray[0, {3, 3}]}, {{1}}, 5][[-1]]]

And using @Mr.Wizard notation with CrossMatrix from the answer above:

Image3D[
   SubstitutionSystem[{1 -> 1 - CrossMatrix[{1, 1, 1}], 
   0 -> ConstantArray[0, {3, 3, 3}]}, {{{1}}}, {3}][[1]]]