# Generating a Sierpinski carpet

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


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


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


• As a newbie MMA user, I agree this seems most natural. +1 – RunnyKine Mar 25 '13 at 19:12
• The quilt is pretty. – whuber Mar 25 '13 at 19:22
• Very nice! 3D version,carpet[n_] := Nest[ArrayFlatten[{{{#, #, #}, {#, 0, #}, {#, #, #}}, {{#, 0, #}, {0, 0, 0}, {#, 0, #}}, {{#, #, #}, {#, 0, #}, {#, #, #}}}, 3] &, 1, n]; Image3D[carpet[3]] – chyanog Mar 26 '13 at 3:35
• +1 not sure how i missed this one. love the graphics - added to my Pinterest HOF pinterest.com/pin/282389839107715883 – Vitaliy Kaurov Jun 2 '13 at 21:00
• @glS Someone linked to reference.wolfram.com/language/guide/… and checking develop.open.wolframcloud.com/app/view/newNotebook I found that it is live on Open Cloud. – Mr.Wizard Mar 10 '17 at 16:36

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


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


• Shamelessly upvoting. – rcollyer Mar 25 '13 at 16:40
• +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):

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


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


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.

• +1 Once you head down this road to generalization, you might as well implement an IFS. – whuber Mar 26 '13 at 18:42
• @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.

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]


• 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. – chyanog 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]]}]]


• 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
• Not with this procedure, due to the stochastic nature of the method used. – J. M. will be back soon May 2 '15 at 3:20

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


• 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. will be back soon 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. will be back soon Apr 10 '13 at 10:26

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