How to create a Poincaré disk type kaleidoscope in Mathematica?

Question

Creating a kaleidoscope in Mathematica is not a new topic at all, examples can be found from the links like Wolfram reference and Wolfram demonstrations.

My question is how to create Poincaré disk type of kaleidoscopes in Mathematica like the guys do in Java. Any idea and suggestion would be gratefully appreciated.

Answer 1

Hyperbolic lines in the Poincare disk are arcs of circles which, if extended, would meet the unit circle boundary orthogonally. Tiling the disk requires reflection of a tile's vertices and circular edges in those same circular edges. This corresponds to inversion in a circle.

Begin with a central polygon having p edges, with q such polygons meeting at each vertex. Then reflect this polygon in each of its circular edges to form p new polygons. Nest this iteration to the depth required.

Here is some code (perhaps too much). With these basics you can add variations like colouring hyperbolic polygons (not straight-edged polygons), or marking tiles with lines at the midpoints (forming the dual tiling), or extending to 3D, or truncating the tile corners to form additional shapes.

The function HyperbolicLine draws hyperbolic lines between points P and Q in the Poincare disk.

HyperbolicLine[{{Px_, Py_}, {Qx_, Qy_}}] :=
   If[N[Chop[Px Qy - Py Qx]] =!= 0.,
      Circle[OrthoCentre[{{Px,Py},{Qx,Qy}}],OrthoRadius[{{Px, Py},{Qx, Qy}}],
             OrthoAngles[{{Px, Py}, {Qx, Qy}}]],
      Line[{{Px, Py}, {Qx, Qy}}]]

The three functions within are as follows.

OrthoCentre[{{Px_, Py_}, {Qx_, Qy_}}] := 
   With[{d = 2 Px Qy - 2 Py Qx, p = 1 + Px^2, q = 1 + Qx^2 + Qy^2}, 
      If[N[d] =!= 0.,
         {p Qy + Py^2 Qy - Py q, -p Qx - Py^2 Qx + Px q}/d, 
         ComplexInfinity]]

OrthoRadius[{{Px_, Py_}, {Qx_, Qy_}}] := 
   If[N[Chop[Px Qy - Py Qx]] =!= 0., 
      Sqrt[Total[OrthoCentre[{{Px, Py}, {Qx, Qy}}]^2] - 1],
      Infinity]

OrthoAngles[{{Px_, Py_}, {Qx_, Qy_}}] := 
   Block[{a, b, c = OrthoCentre[{{Px, Py}, {Qx, Qy}}]},
      If[(a = N[Apply[ArcTan, {Px, Py} - c]]) < 0., a = a + 2 π]; 
      If[(b = N[Apply[ArcTan, {Qx, Qy} - c]]) < 0., 
      b = b + 2 π]; {a, b} = Sort[{a, b}]; 
      If[b - a > π, {b, a + 2 π}, {a, b}]]

The Inversion function inverts points, lines, and polygons in a circle.

Inversion[Circle[{Cx_, Cy_}, r_], {Px_, Py_}] := {Cx, Cy} + 
   r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)
Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], {Px_, Py_}] := {Cx, Cy} + 
   r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Line] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Polygon] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], {Ux_, Uy_}] := 
   With[{u = Px - Qx, v = Qy - Py},
      {-Ux (v^2 - u^2) - 2 u v Uy, Uy (v^2 - u^2) - 2 u v Ux}/(u^2 + v^2)]
Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], p_Polygon] :=
   Map[Inversion[Line[{{Px, Py}, {Qx, Qy}}], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_], c_List] := 
   Map[Inversion[Circle[{Cx, Cy}, r], #] &, c]

Invert points of a polygon or list of polygons in a hyperbolic line (which is a line or a circle).

PolygonInvert[p_Polygon] := 
   Map[Inversion[HyperbolicLine[#], p] &, 
       Partition[Join[p[[1]], {p[[1, 1]]}], 2, 1]]
PolygonInvert[p_List] := Flatten[Map[PolygonInvert[#] &, p]]

Join the points of a hyperbolic polygon with straight or hyperbolic lines.

LineRule = Polygon[x_] :> Line[Join[x, {x[[1]]}]];
HyperbolicLineRule = 
   Polygon[x_] :> 
      Map[HyperbolicLine, Partition[Join[x, {x[[1]]}], 2, 1]];

The function CentralPolygon finds the initial polygon to begin the iterations. The first point of the polygon (with p sides) is in the first quadrant at angle π/p with respect to the horizontal. Each vertex is surrounded by q such polygons. The polygon may be rotated counter-clockwise through ϕ radians by using the optional third argument.

CentralPolygon[p_Integer, q_Integer, ϕ_: 0] := 
   With[{r = (Cot[π/p] Cot[π/q] - 1)/Sqrt[Cot[π/p]^2 Cot[π/q]^2 - 1],
        θ = π Range[1, 2 p - 1, 2]/p}, 
        r Map[{{Cos[ϕ], -Sin[ϕ]}, {Sin[ϕ], Cos[ϕ]}}.# &,
             Transpose[{Cos[θ], Sin[θ]}]]]

Delete duplicate polygons, allowing a tolerance.

PolygonUnion[p_Polygon, tol_: 10.^-10] := p
PolygonUnion[p_List, tol_: 10.^-10] := 
   With[{q = p /. Polygon[x_] :> N[Polygon[Round[x, 10.^-10]]]}, 
        DeleteDuplicates[q]]

HyperbolicTessellation tiles polygons with p sides, q of them meeting at each vertex, nesting down to level k. The entire tessellation is rotated counter-clockwise by ϕ radians.

HyperbolicTessellation[p_Integer, q_Integer, ϕ_, k_Integer, t_: 10.^-10] := 
   Map[PolygonUnion[#, t] &, 
       NestList[PolygonInvert, Polygon[CentralPolygon[p,q,ϕ]],k][[{-2, -1}]]] /; k > 0

HyperbolicTessellation[p_Integer, q_Integer, ϕ_, k_Integer, t_: 10.^-10] := 
   Polygon[CentralPolygon[p, q, ϕ]] /; k == 0

And finally,

HyperbolicTessellationGraphics[p_Integer, q_Integer, ϕ_, k_Integer,
      rule_RuleDelayed, opts___] := 
   Graphics[{Circle[{0, 0}, 1], 
      HyperbolicTessellation[p, q, ϕ, k, 10.^-10] /. rule}, opts]

For example,

HyperbolicTessellationGraphics[7, 3, 0., 3, HyperbolicLineRule, PlotLabel -> "{7,3}"]

Graphics[{EdgeForm[{Thin, Black}], 
   Flatten[HyperbolicTessellation[5,4,0.,4]] /. Polygon[x_] :>
      {ColorData["DarkRainbow", Norm[Mean[x]]], Polygon[x]}}]

Block[{poly = Flatten[HyperbolicTessellation[7, 3, 0., 4]]}, 
   Graphics[{EdgeForm[{Thin, Black}], 
      poly /. Polygon[x_] :> Polygon[x, VertexColors -> 
       Map[ColorData["SouthwestColors", Norm[#]^1.2 + RandomReal[{0.,0.3}]]&,x]]}]]

And an example of truncation of a heptagon tiling, left as an exercise...