Poincaré Disk Hyperbolic Tessellation

Math essence

Reflection across a geodesic whose endpoints on the unit circle are (a, b): R(z) = (c*conj(z) + (a-b)*a) / (conj(c)*z + conj((a-b)*a)) where c = (a-b) / (1 - conj(a)*b). Central polygon radius: cosh(rho) = cos(pi/q) / sin(pi/p). Schlafli deficit: (p-2)(q-2) > 4 required for hyperbolic geometry.

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/openprocessing_p5js/HyperbolicTiling/reference.py

from __future__ import annotations

from typing import ClassVar

import numpy as np


def _mobius_reflect_point(z: complex, a: complex, b: complex) -> complex:
    """Reflect complex point ``z`` across the geodesic line through points
    ``a`` and ``b`` (inside the unit disk).

    A Poincaré-disk geodesic is a Euclidean circle orthogonal to the unit
    circle. Given two interior points a, b, the unique such circle passing
    through them has center ``c`` on the perpendicular bisector of ab with
    the orthogonality constraint ``|c|^2 = |c - a|^2 + 1``. Solving:

        Re(conj(a) c) = (|a|^2 + 1) / 2
        Re(conj(b) c) = (|b|^2 + 1) / 2
        Im(conj(a) c - conj(b) c) free

    This is a 2×2 linear system in the real and imaginary parts of c. The
    reflection is then inversion in that circle:

        R(z) = c + |c - a|^2 / conj(z - c)

    Used by Joyce's Hyperbolic Tessellations applet (1995), b5strbal/Escher
    (2020), and Beardon, *The Geometry of Discrete Groups*, §3.
    """
    if abs(a - b) < 1e-12:
        return z
    aa = abs(a) ** 2
    bb = abs(b) ** 2
    rhs_a = (aa + 1.0) * 0.5
    rhs_b = (bb + 1.0) * 0.5
    A = np.array(
        [
            [a.real, -a.imag],
            [b.real, -b.imag],
        ],
        dtype=np.float64,
    )
    rhs = np.array([rhs_a, rhs_b], dtype=np.float64)
    try:
        sol = np.linalg.solve(A, rhs)
    except np.linalg.LinAlgError:
        return z
    c = complex(sol[0], sol[1])
    rho2 = abs(c - a) ** 2
    dz = z - c
    if abs(dz) < 1e-12:
        return z
    return c + rho2 / np.conj(dz)


class HyperbolicTiling:
    """Escher-style regular hyperbolic tessellation in the Poincaré disk.

    Builds the ``{p, q}`` tiling where ``p``-sided regular polygons meet
    ``q`` at every vertex, valid only when the Schlafli angle deficit
    ``(p - 2)(q - 2) > 4``. Implemented as a recursive reflection scheme in