Speeding up RegionPlot for iterated function system fractals

Question

This is my first post here, I have the following basic code for drawing an iterated function system fractal. It works to the third step, but then for the fourth iteration it freezes my old MacBook. I have not included the fourth iteration but it should be clear how to proceed. I reckon it has something to do with the way RegionPlot works. I am looking for an alternative or a way to make it quicker. I have also tried to use DiscretizeRegion to no effect.

ClearAll["Global`*"]
θ = ArcTan[1/2];

sqCreate[{x0_, y0_}, h_] := Rectangle[{x0, y0}, {x0 + h, y0 + h}];
simil[i_, sq_] := Which[
   i == 1, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {0, 1/√5}}]],
   i == 2, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {1/5 + 2/5, 1/√5 - 1/5 + 2/5}}]],
   i == 3, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5 + 2/5, 1/√5 - 1/5 - 1/5}}]],
   i == 4, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5 - 1/5, 1/√5 - 2/5 - 1/5}}]],
   i == 5, 
   TransformedRegion[sq, AffineTransform[{1/√5 RotationMatrix[-θ],
                                         {2/5, 1/√5 - 1/5}}]]
   ];

fullSimil[sq_] := 
  RegionUnion[simil[1, sq], simil[2, sq], simil[3, sq], simil[4, sq], 
   simil[5, sq]];

set0 = sqCreate[{0, 0}, 1];
paint0 = Graphics[{Gray, set0}]

set1 = fullSimil[set0]
paint1 = RegionPlot[set1, PlotRange -> All]

set2 = fullSimil[set1];
paint2 = RegionPlot[set2, PlotRange -> {{0, 1.5}, {-0.5, 1}}, 
  AspectRatio -> Automatic]

set3 = fullSimil[set2];
paint3 = RegionPlot[set3, PlotRange -> {{0, 1.5}, {-1, 1}}, 
  AspectRatio -> Automatic]

Welcome to Mathematica.SE, plebmatician! I suggest the following:

As you receive help, try to give it too, by answering questions in your area of expertise.

Take the tour and check the faqs!

When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge.

Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! — Chris K, Jan 04 '19 at 18:52

Henrik Schumacher · Accepted Answer · 2019-01-04T18:20:29.700

6

This should work a bit faster. The trick is to use BoundaryDiscretizeRegion in each step to convert the internal representation to a single polygonal line. The second ingedient is to use machine precision numbers; these are much faster for numerical computations.

θ = ArcTan[0.5];
A = 1/Sqrt[5] RotationMatrix[-θ];
vecs = Developer`ToPackedArray@N@{{0, 1/Sqrt[5]}, {1/5 + 2/5, 1/Sqrt[5] - 1/5 + 2/5}, {2/5 + 2/5, 1/Sqrt[5] - 1/5 - 1/5}, {2/5 - 1/5, 1/Sqrt[5] - 2/5 - 1/5}, {2/5, 1/Sqrt[5] - 1/5}};

step[sq_] := BoundaryDiscretizeRegion[
  RegionUnion@Table[TransformedRegion[sq, AffineTransform[{A, v}]], {v, vecs}]
  ]

MList = NestList[
  step,
  MeshRegion[N@{{0, 0}, {1, 0}, {1, 1}, {0, 1}}, Polygon[{{1, 2, 3, 4}}]],
  8
  ];

GraphicsGrid[
 Partition[MList, 3],
 ImageSize -> Full
 ]

edited Jan 04 '19 at 18:20

answered Jan 04 '19 at 17:59

Henrik Schumacher

106,770
7
179
309

Thank you very much, amazing how much faster this is. But I don't understand the command "Developer`ToPackedArray", never seen syntax like it before. And what is its purpose, wouldn't the code work the same without it? – plebmatician Jan 04 '19 at 19:14
2

Developer`ToPackedArray just converts the array into a so-called packed array. Packed arrays are usually helpful for numerical performance. But here Developer`ToPackedArray makes probably no difference. – Henrik Schumacher Jan 04 '19 at 19:29

score 2 · Answer 2 · answered Jan 05 '19 at 15:26

With a little refactoring of your code, and with the use of some undocumented functionality:

θ = ArcTan[1/2];
sqCreate[{x0_, y0_}, h_] := Polygon[{{x0, y0}, {x0 + h, y0}, {x0 + h, y0 + h}, {x0, y0 + h}}]

fullSimil[sq_Polygon] := Map[First, Graphics`PolygonUtils`PolygonCombine[
                      Table[Polygon[AffineTransform[{RotationMatrix[-θ]/Sqrt[5], pos}] @@ sq],
                            {pos, {{0, 1/Sqrt[5]}, {3/5, 1/Sqrt[5] + 1/5}, {4/5, 1/Sqrt[5] - 2/5},
                                   {1/5, 1/Sqrt[5] - 3/5}, {2/5, 1/Sqrt[5] - 1/5}}}]]]

Partition[BoundaryDiscretizeGraphics /@ NestList[fullSimil, sqCreate[{0, 0}, 1], 8],
          3] // GraphicsGrid

$fractal meshes$

Speeding up RegionPlot for iterated function system fractals

2 Answers2