50
$\begingroup$

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

10 Answers 10

82
$\begingroup$

Version 11.1 introduces MengerMesh:

MengerMesh[3]

enter image description here

enter image description here



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

enter image description here


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

$\endgroup$
13
  • $\begingroup$ As a newbie MMA user, I agree this seems most natural. +1 $\endgroup$
    – RunnyKine
    Mar 25, 2013 at 19:12
  • 3
    $\begingroup$ The quilt is pretty. $\endgroup$
    – whuber
    Mar 25, 2013 at 19:22
  • 3
    $\begingroup$ Very nice! 3D version,carpet[n_] := Nest[ArrayFlatten[{{{#, #, #}, {#, 0, #}, {#, #, #}}, {{#, 0, #}, {0, 0, 0}, {#, 0, #}}, {{#, #, #}, {#, 0, #}, {#, #, #}}}, 3] &, 1, n]; Image3D[carpet[3]] $\endgroup$
    – chyanog
    Mar 26, 2013 at 3:35
  • 1
    $\begingroup$ +1 not sure how i missed this one. love the graphics - added to my Pinterest HOF pinterest.com/pin/282389839107715883 $\endgroup$ Jun 2, 2013 at 21:00
  • 1
    $\begingroup$ @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. $\endgroup$
    – Mr.Wizard
    Mar 10, 2017 at 16:36
24
$\begingroup$

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

$\endgroup$
18
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ Shamelessly upvoting. $\endgroup$
    – rcollyer
    Mar 25, 2013 at 16:40
  • 1
    $\begingroup$ +1 Using ArrayPlot instead of ColorNegate@Image in the first solution gives a crisper rendering (and happens to be about 50% faster, too). $\endgroup$
    – whuber
    Mar 25, 2013 at 17:40
18
$\begingroup$

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

$\endgroup$
8
  • $\begingroup$ 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.) $\endgroup$
    – Mr.Wizard
    Mar 25, 2013 at 19:14
  • $\begingroup$ @Mr.Wizard I gave the code for the first one :-). $\endgroup$
    – whuber
    Mar 25, 2013 at 19:17
  • $\begingroup$ Thanks. v7 users FYI: you'll need black = ColorConvert[black, "RGB"] $\endgroup$
    – Mr.Wizard
    Mar 25, 2013 at 19:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – Mr.Wizard
    Mar 25, 2013 at 19:28
  • 2
    $\begingroup$ @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! $\endgroup$
    – whuber
    Mar 26, 2013 at 14:17
11
$\begingroup$

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

$\endgroup$
8
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ +1 Once you head down this road to generalization, you might as well implement an IFS. $\endgroup$
    – whuber
    Mar 26, 2013 at 18:42
  • 2
    $\begingroup$ @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. $\endgroup$ Mar 26, 2013 at 18:52
  • $\begingroup$ @whuber or an L-system. I considered adapting my bf parser to do just that. $\endgroup$
    – rcollyer
    Mar 27, 2013 at 12:17
4
$\begingroup$

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

$\endgroup$
3
  • 1
    $\begingroup$ Please see the most recent section of my answer. I believe if you run Timings you will find it is superior. $\endgroup$
    – Mr.Wizard
    Mar 27, 2013 at 18:32
  • $\begingroup$ @Mr.Wizard I'm impressed. $\endgroup$
    – chyanog
    Mar 28, 2013 at 5:03
  • $\begingroup$ @chyanog There seems to be an error in the second piece of code (version 12.1). $\endgroup$ Nov 18, 2020 at 7:12
3
$\begingroup$

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

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$
    – Lady InRed
    Aug 19, 2013 at 9:56
  • $\begingroup$ Not with this procedure, due to the stochastic nature of the method used. $\endgroup$ May 2, 2015 at 3:20
2
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ In your picture the sponge looks faded, and on my system it renders as a solid gray cube. $\endgroup$
    – Mr.Wizard
    Apr 10, 2013 at 8:55
  • $\begingroup$ 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. $\endgroup$ Apr 10, 2013 at 9:50
  • $\begingroup$ 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? $\endgroup$
    – Mr.Wizard
    Apr 10, 2013 at 10:07
  • $\begingroup$ @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. $\endgroup$ Apr 10, 2013 at 10:26
2
$\begingroup$

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

enter image description here

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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