More Circle Magic

Today we revisit the three previous posts about playing with circles with Mathematica Code, as promised.

First, we set up a reference image we want to play with:

width = 3;
height = 2;
spacing = .2;
vectorrectangle = Graphics[{Thickness[.01],
     Line[{{i, 0}, {i, height}}, VertexColors -> {Red, Green}], {i, 0,
       width, spacing}],
     Line[{{0, j}, {width, j}}, VertexColors -> {Blue, Yellow}], {j, 
      0, height, spacing}]


and convert this vector graphics into a bit image:

a = 500 width;
b = 500 height;
bitimage = Image[vectorrectangle, ImageSize -> {a, b}];
rectangleimage = ImageData[bitimage];

Instead of this, you can paste any bitmap into the Image[] function and adjust a and b to the correct dimensions.
ImageData converts the bitmap into an array of RGB values.

Now we do some complex analysis magic, using elliptic functions, in order to map the unit disk to the rectangle.

tau = b/a I;
modul = 1./(ModularLambda[2 tau])^(1/4);
ellipticwidth = (2*EllipticK[1/modul^4])/modul;

DiskToUpperHalfPlane[z_] := (1 - I z)/(-I + z);
UpperHalfPlaneToRectangle[z_] := 
  EllipticF[ArcSin[modul z], 1/modul^4]/modul;
UpperHalfPlaneToNormalizedRectangle[z_] := 
  1/2 + UpperHalfPlaneToRectangle[z]/ellipticwidth;
imre[z_] := {Im[z], Re[z]};


To compute the circular preimage of the rectangle, we create a 750×750 square and map each pixel of the disk within that square to the rectangle, using DiskToNormalizedRectangle, and some rescaling. We record in the square the coordinates of the image pixels, and mark all other pixels in the square that lie not inside the disk as {0,0}.

nn = 750;
scale[{i_, j_}] := 2 (i + I j)/nn - (1 + I);

tab = Transpose[Table[
    If[Abs[scale[{i, j}]] > .99,
     {0, 0},
     imre[a (DiskToNormalizedRectangle[scale[{i, j}]])]], 
	{i, 1, nn}, {j, 1, nn}]];

As the image coordinates are not integers, we interpolate the color of a pixel in the disk as the weighted average of the four adjacent pixels.

WeightedExtract[t_, {x_, y_}] := (1 - FractionalPart[x]) (1 - FractionalPart[y]) t[[Floor[x], Floor[y]]] +
  (1 - FractionalPart[x]) (FractionalPart[y]) t[[Floor[x],  Ceiling[y]]] +
  (FractionalPart[x]) (1 - FractionalPart[y]) t[[Ceiling[x], Floor[y]]] +
  (FractionalPart[x]) (FractionalPart[y]) t[[Ceiling[x], Ceiling[y]]];

colorpick[i_] := 
 If[i == {0, 0}, RGBColor[1, 1, 1], 
  RGBColor @@ WeightedExtract[rectangleimage, i]];

Image[Map[colorpick, tab, {2}]]


To use Black products to create multiple copies of this image within the disk, we first generate a list of the zeroes of the Blasche factor, and then multiply the Black factors. Precomposing the Blaschke function with DiskToNormalizedRectangle makes things pretty.

ww = Table[(1 - 1/n^1.4) E^(.4 I n), {n, 2, 20}];
Blaschke[z_] := Product[
  Abs[ww[[n]]]/ww[[n]] (ww[[n]] - z)/(1 - Conjugate[ww[[n]]] z), {n, 
   1, Length[ww]}];

tab = Transpose[Table[
    If[Abs[scale[{i, j}]] > .99,
     {0, 0},
     imre[a (DiskToNormalizedRectangle[scale[Blaschke[{i, j}]]])]], 
	{i, 1, nn}, {j, 1, nn}]];

Image[Map[colorpick, tab, {2}]]


To get this back into a rectangle, we need the inverse of the DiskToNormalizedRectangle function. This is done with Jacobi’s elliptical functions:

UpperHalfPlaneToDisk[w_] := (I (-I + w))/(I + w);
RectangleToUpperHalfPlane[w_] := JacobiSN[modul w, 1/modul^4]/modul;
NormalizedRectangleToUpperHalfPlane[w_] := -RectangleToUpperHalfPlane[-ellipticwidth (w - 1/2)];

RectangleSelfMap[z_] := UpperHalfPlaneToNormalizedRectangle[

Again, after some bookkeeping

aa = a/2;
bb = b/2;
eps = 0.000001;

scale2[{i_, j_}] := (i + I j)/aa
smash[x_, u_] := If[x  u, u, x]];
smash[{i_, j_}] := {smash[i, b], smash[j, a]};
colorpick2[i_] := colorpick[smash[i]];

tab = Transpose[Table[
    imre[a RectangleSelfMap[(1 - eps) scale2[{i, j}] + eps/2 b/a I ]], {i, 1, aa}, {j, 1, bb}]];
Image[Map[colorpick2, tab, {2}]]

we get our rectangular image back.


Playing with Infinity

After Squaring the Circle and using Blaschke products to squeeze finitely many disks conformally into a single disk, the restless mathematician wants to put even more disks into a disk. Infinitely many, that is. What we need for this is a conformal map that maps the disk infinitely often onto itself. We will cheat a little and employ the exponential function. It maps the left half-plane to the punctured unit disk, and because it is 2π periodic. Each point except the center of the unit disk has infinitely many preimages.

Expo 1

So, my first examples map the disk via a Möbius transformation to the left half plane (that is a bijection), and then via the exponential map back to the (punctured) disk. Missing one point won’t be a big deal – what is a single pixel? Pulling back the default spider web on the left replicates it infinitely often, but in a somewhat disappointing way: The circles of the original web become horocycles that are mapped periodically onto the circles. We just have an infinite hyperbolic wall paper, periodic with respect to a parabolic isometry. We can do slightly better by taking products of several such maps.

Below are products of three such Möbius-exponential maps,

Expo 3

and here the two images that use four factors. I have placed the parabolic fixed points symmetrically but played with the Möbius transformations a little.

Of course one can also use infinite Black products, as long as the zeroes of the factors converge rapidly enough to the boundary circle, and as long one is patient enough to evaluate the infinite products to sufficient accuracy.


In the left example, the zeroes follow a spiral, while in the second example, the zeroes alternate between two spirals that turn the opposite way. The zeroes correspond to the circular holes.

Playing with circular images

After successfully transforming rectangular images into circular ones it is time to do something with them. We have seen already that the one can deform them by shifting one point somewhere else. This is very much like rotating a globe.

But besides these angle preserving symmetries of the disk there are other maps from the disk to itself that also preserve angles but are not anymore 1:1. These are the Blaschke products, written in complex notation as follows.

Product 01

Let’s look at a simple example with just two factors, and choose the a-parameters to be 1/2 and -1/2. Then B(z) maps the double spiderweb on the right to the standard spiderweb on the left:

Blascke2 01

In other words, by taking preimages (or better, by using B(z) to pull back an image…), we can create multiple copies of a circular image within a circular image. The a-parameters designate the locations of the “centers” of the multiple spiderwebs where the strands converge.


For instance, above is a circular image of a Spring wild flower, and to the left its 3-fold mutation. Below are 5-fold mutations with two different choices for the a-parameter.


These images resemble kaleidoscopes, but are improved, because the copies of the original image fit together more evenly (smoothly, and not just by reflection). One can also make the result less symmetrical by choosing the a-parameters less symmetrical. Below the copies of the ferns are places at 120 degree angles but differently far out,


and here we have a large copy of the original budding trillium at the bottom with two smaller copies to the left and right.


Now I need to find somebody who writes an app that implements all this…