Generating Doyle spiral painting

TL;DR Yes, this can be done!

Mathematica graphics Mathematica graphics


If you read the article "Hexagonal circle packings and Doyle spirals" by Leys, you will see that for a choice p and q, we need to find the complex values A, B and r. For that purpose, we can steal this part from the demonstration you linked:

doyle[pi_, qi_] := Module[{p = pi, q = qi, s, t, r},
  r[s_, t_, p_, 
    q_] := (s^2 + s^(2 p/q) - 
      2 s^((p + q)/q) Cos[(2 \[Pi] + p t - q t)/q])/(s + s^(p/q))^2;
  {s, t} = {s, t} /. FindRoot[
     {r[s, t, 0, 1] - r[s, t, p, q] == 0,
      r[s, t, 0, 1] - r[s^(p/q), (p*t + 2. Pi)/q, 0, 1] == 0},
     {{s, 1.5}, {t, 0}}];
  {s*Exp[I*t], s^(p/q)*Exp[I (p*t + 2*Pi)/q], Sqrt[r[s, t, 0, 1]]}
  ]

{a, b, r} = doyle[2, 12]

Now we have the centers for both additional complex circles. Knowing a bit of complex analysis, one understands that for creating all packed circles, we only need to iteratively multiply by e.g. a. So we could write down a function that does this multiplication. I'm keeping complex values until the final visualization, where we use the real part for x and the imaginary part for y:

iterate[a_, b_, j_, n_] := Module[{start = b^j},
  Table[a^i*start, {i, Range[-n, n]}]
]
Graphics[Circle[ReIm[#], r*Abs[#]] & /@ iterate[a, b, 0, 3]]

Mathematica graphics

This shows the $0th$ spiral of packed circles, $3$ circles inward to our base circle and $3$ circles outward. To create the complete plane, we have to create $12$ columns, since q was $12$:

toCirle[z_, r_] := Disk[{Re[z], Im[z]}, Abs[z]*r];
pack = Table[iterate[a, b, j, 5], {j, 12}];
gr = Graphics[{EdgeForm[Black], 
  Map[{RandomColor[], toCirle[#, r] & /@ #} &, pack]}, 
 PlotRange -> {{-10, 10}, {-10, 10}}]

Mathematica graphics

Unfortunately, that is not enough, because the artist chose to use the logarithmic spirals through the circle centers as guide to divide each circle into different parts. In order to do this, we need to go further.

Let us make a cut here and divide the following into small sections where we look at the details. These details will be important for the overall approach

Connection between a, b, r and the circles

As pointed out in the article, the complex numbers a and b are the generators of the circles. This means, all circle centers can be obtained by repeated multiplication. The base circle with the center {1,0} is given by $a^0\cdot b^0$ which is 1 (meaning re=1 and im=0).

Now each multiplication by a or by b gives the center of the circle that is next to the base circle. So 1*a*a=a^ gives the second circle in the direction of a. a^2*b shifts this last circle in the direct of b.

Note, that even a^(-3) is perfectly OK and gives the 3rd circle in the oposite direction. These are the small circles the fill the center. OK, one Manipulate says more than a thousand words. Let us create a dynamic table of all circles in a range for a^i*b^j. Note that, as pointed out in the article, the correct radius for each circle is Abs[a^i*b^j]*r where r is the radius we got from the solution of doyle.

Manipulate[
 Graphics[{EdgeForm[Black], FaceForm[Gray], Table[
    toCirle[a^i*b^j, r],
    {i, i1, i2}, {j, j1, j2}]}, PlotRange -> {{-10, 10}, {-10, 10}}],
 {{i1, 0}, -5, 3, 1},
 {{i2, 0}, i1, 5, 1},
 {{j1, 0}, -5, 3, 1},
 {{j2, 0}, j1, 5, 1}
 ]

Mathematica graphics

Circles and spirals

We have seen that we can go from one circle to any neighbour by increasing (or decreasing) either i or j by 1. But what if we don't make jump from say a^3 to a^4 but a smooth transition? Well, the function for such a thing is easy because a^0=1 and a^1=a, so we can make a function a^3*a^t and let t run from 0 to 1.

Show[
 Graphics[{
   FaceForm[Gray],
   toCirle[#, r] & /@ {a^3, a^4}}
  ],
 ParametricPlot[ReIm[a^3*a^t], {t, 0, 1}, PlotStyle -> White]
 ]

Mathematica graphics

This looks very much like the spirals that were used to divide the circles in the original art. So it seems if we pick out the center any circle next to our base circle, we can create spiral functions that go through the circles. Note that the approach of shifting a spiral to its neighbouring spiral is similar to shifting circles. Here is an example:

Show[
 gr,
 ParametricPlot[Table[ReIm[b^i a^t], {i, 12}], {t, -10, 10}, 
  PlotRange -> {{-10, 10}, {-10, 10}}, PlotStyle -> White],
 ParametricPlot[Table[ReIm[a^j b^t], {j, 5}], {t, -10, 10}, 
  PlotRange -> {{-10, 10}, {-10, 10}}, PlotStyle -> White]
 ]

Mathematica graphics

Spiral functions inside circles

For our later approach, I want to be able to draw the spiral only inside a circle. As we have seen, going from t=0 to t=1 will connect the centers of the circles. This is not what we want. We want values for t that start and end with the circle. Let's make the plot we did earlier again, but use values for t between -1/3 and 1/3

Mathematica graphics

OK, that looks promising. Remember, we know the center of this circle with a^3 and we know its radius with Abs[a^3]*r. What are the bounds where our spiral is exactly on the radius? Let us ask FindRoot:

tb = t /. FindRoot[Abs[1 - a^t] - r, {t, #}] & /@ {-1/3, 1/3}
(* {-0.565183, 0.433533} *)

But wait! I haven't used a^3 at all! Correct. The good thing is that the bounds for the circles apply to each circle of the same spiral. Therefore I'm using the next neighbour of the base circle which is a for FindRoot. Look here:

Show[
 Graphics[{
   FaceForm[Gray],
   toCirle[#, r] & /@ {a^3, a^4}}
  ],
 ParametricPlot[ReIm[a^3*a^t], {t, tb[[1]], tb[[2]]}, 
  PlotStyle -> White]
 ]

Mathematica graphics

What spirals did the artist use?

As it turns out she used the spirals of the following direct neighboring circles of the base circle:

spoints = {a*b^-1, a, b}
(* {1.46301 - 0.54185 I, 1.67036 + 0.343254 I, 0.927594 + 0.578172 I} *)

Let's make a small function that calculates their bounds and returns them with a spiral function. The spiral function will directly incorporate the i and j so that we can easily draw it on every circle we like

spiral[pt_] := Module[{t1, t2},
  {t1, t2} = 
   Block[{t}, 
    t /. FindRoot[Abs[1 - pt^(t)] - r, {t, #}] & /@ {-1/3, 1/3}];
  {t1, t2, Function[{i, j, t}, a^i*b^j*pt^t]}
  ]

Now let's plot these 3 spirals inside our base circle {1,0}

Show[{
  Graphics[Circle[{1, 0}, r]],
  ParametricPlot[ReIm@#3[0, 0, t], {t, #1, #2}] & @@@ spiral /@ spoints
  }]

Mathematica graphics

Now, we can calculate the points of the spirals inside each circle, we have the radius of each circle and through the spiral's start and end points, we have 6 points on each circle.

Creating polygons points for the parts of a circle

For each cake-part of a circle, we can now create a polygon by

  • starting in the center
  • creating points along a spiral outwards to the circle boundary
  • going counterclockwise along the circle to the endpoint of the next spiral
  • create points along this next spiral from the outer point to the center

However, one tiny point is missing. How do we create points along the circle from when we go from one spiral point to the next. That is not as hard as it sounds.

Assume you have two (complex) points that lie on a circle around a center. Then you can subdivide them and create arbitrarily many points between them that all lie on the circle.

circle[z1_, z2_, cent_] := 
 Module[{zz1 = z1 - cent, zz2 = z2 - cent, r},
  r = Abs[Mean[{zz1, zz2}]];
  # + cent & /@ 
   Nest[Riffle[#, 
      Function[zz, With[{m = Mean[zz]}, m/Abs[m]*Abs[zz1]]] /@ 
       Partition[#, 2, 1]] &, {zz1, zz2}, 5]
  ]

Having this, we can create the points for all cake-parts of circle i, j defined by the provided spirals that divide the circle:

createCircleParts[spirals_, i_, j_] := 
 Module[{center, outward, inward},
  outward = Table[#3[i, j, t], {t, 0, #2, #2/10.}] & @@@ spirals;
  inward = Table[#3[i, j, t], {t, 0, #1, #1/10.}] & @@@ spirals;
  center = a^i*b^j;
  {i, j,
   Join[#1[[;; -2]], circle[#1[[-1]], #2[[-1]], center], 
      Reverse[#2[[;; -2]]]] & @@@ 
    Partition[Join[outward, inward, {First[outward]}], 2, 1]}
  ]

The function returns {i,j, {part1, part2, ...}} and we will use i and j later for the coloring as it gives us information about the position of the circle.

To test this function, let us see what happens with the circle i=1 and j=2:

Graphics[{RandomColor[], Polygon[ReIm[#]]} & /@
  Last@createCircleParts[spiral /@ spoints, 1, 2]
 ]

Mathematica graphics

Coloring of circles

For one circle we have the information i, j which encodes the global position and of course we have n cake-parts. An easy way would be to provide a list of colors and select a color depending on the information we have.

I could not really find a pattern in the coloring of the artists image, so lets keep it simple but let us use equivalent colors:

cols = {Black, RGBColor[0.078, 0.71, 0.964], Orange, Red, Darker@Green, Purple};
colorCircleParts[{i_, j_, parts_}, col_List] := 
 Table[{col[[Mod[i + j + n, Length[col]] + 1]], 
   Polygon[ReIm@parts[[n]]]}, {n, Length[parts]}]

Putting everything together

The last thing we need to do is to create a table containing the circles and their parts for a range of i and j values. Then we color the circle parts and display them:

all = Table[colorCircleParts[createCircleParts[spiral /@ spoints, i, j], 
    cols], {i, -5, 6}, {j, 0, 12}];
Graphics[all, PlotRange -> {{-20, 20}, {-20, 20}}]

Mathematica graphics

Aftermath: Getting something close to the artist's work

The webpage of the artist suggests that

The central part of the picture shows the Doyle spiral circle packing P = 2 Q = 12.

That is not true. The values of P and Q define how many circles you need to close one loop. Additionally, the rotation of the circles in the artist's work is clockwise while in mathematics, we usually prefer to go counter-clockwise.

Lucky for us, this is no big deal because to go clockwise we just need to conjugate our complex values a and b. After printing the painting and counting the circles (and paying absolutely no attention to Wjx's comment who already found out that the values are off), I discovered that the painting uses P=3 and Q=8.

Let me show you what that means:

pqPlot[p_, q_] := Module[{a, b, r, c1, c2},
  {a, b, r} = doyle[p, q];
  {a, b} = Conjugate /@ {a, b};
  c1 = toCirle[#, r] & /@ NestList[a*# &, 1, p - 1];
  c2 = toCirle[#, r] & /@ NestList[b*# &, 1, q - 1];
  Graphics[{EdgeForm[Black], FaceForm[LightYellow], c2, 
    FaceForm[LightBlue], c1, FaceForm[LightGreen], EdgeForm[Thick], 
    toCirle[1, r], toCirle[a^p, r]}]
  ]

pqPlot[3, 8]

Mathematica graphics

If you include the inner base circle in your counting, you have 3 circles in the first and 8 circles in the other direction until you reach the outer end circle.

Taking this into account and including some of the colors in the original painting, we can come up with a very close optical copy of what the artist did. I played around with the plot-range to make it fit.

{a, b, r} = doyle[3, 8];
{a, b} = Conjugate /@ {a, b};
spoints = {a*b^-1, a, b};
cols = {GrayLevel[0.1], RGBColor[0.078, 0.71, 0.964], 
   RGBColor[0.95, 0.36, 0.09], RGBColor[0.77, 0.17, 0.12], 
   RGBColor[0.07, 0.6, 0.25], RGBColor[.32, .24, .55]};
range = 5.585;

Graphics[
 Table[colorCircleParts[createCircleParts[spiral /@ spoints, i, j], 
   cols], {i, -5, 2}, {j, 0, 7}],
 PlotRange -> {{-range, range}, {-range, range}}
 ]

Mathematica graphics


(Not an answer but an extended comment.)

Simple mush-ups

In view of question's main part:

Furthermore, can the code be written to generate a parametrized version of the painting so that one could produce and enjoy an unlimited number of pieces of art?

using halirutan's great answer we can do some mush-ups.

Here is an example.

AbsoluteTiming[
 doyleSpiralImages = Flatten@
    Table[(lcols = RandomSample[cols];
      all = 
       Table[colorCircleParts[
         createCircleParts[spiral /@ spoints, i, j], lcols], {i, -5, 
         6}, {j, 0, 12}]; 
      Image[Graphics[all, 
        PlotRange -> {{-20, 20}, {-20, 20}}]]), {20}];
 ]
doyleSpiralImagesBW = 
  ColorConvert[#, "Grayscale"] & /@ doyleSpiralImages;
doyleSpiralImagesBWBin = 
  Flatten@Table[
    Binarize[#, b] & /@ doyleSpiralImagesBW, {b, {0.15, 0.45}}];

(* {31.5697, Null} *)

AbsoluteTiming[
 nc = 16;
 directBlendingImages = 
  Table[ImageAdjust[
    Blend[Colorize[#, 
        ColorFunction -> 
         RandomChoice[{"BrightBands", "FruitPunchColors", 
           "AvocadoColors"}]] & /@ 
      RandomChoice[doyleSpiralImagesBWBin, nc], 
     RandomReal[1, nc]]], {25}];
 ]

(* {42.7025, Null} *)

Multicolumn[ColorNegate /@ directBlendingImages, 6]

enter image description here

Mush-ups with larger collections of images

Ideally, we can use procedures in the spirit of the one above over a large collection of Doyle spiral images. Such images can be made from scratch or obtained from the Web.

For example, if such an image collection is fed to the neural net functions in Mathematica / WL, in principle we will be able to obtain new spiral images by examining the layers.