Poincare-disk hyperbolic tiling (Escher / Coxeter)

Math essence

Poincare disk: ds^2 = 4|dz|^2 / (1-|z|^2)^2. The (p, q) tiling is generated by R = rotation by 2*pi/p and S = reflection across the line through the origin at angle theta = pi/p + pi/q - pi/2; S(z) = exp(2i*theta) * conj(z). For each sampled point z, a coset expansion iterates S R^k (up to 2p steps) until the point lands in the fundamental triangle, tracking the tile index. Hyperboloid lift: (x, y, z) = (2x, 2y, 1 + x^2 + y^2) / (1 - x^2 - y^2) embeds H^2 in R^3. Project orthographically with pitch/yaw and splat into a frame buffer.

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/fractals_3d_webgl/HyperbolicTilingPoincare3D/reference.py

from __future__ import annotations

from typing import ClassVar

import numpy as np


def _mobius_add(z: np.ndarray, a: np.ndarray) -> np.ndarray:
    """Mobius (hyperbolic) addition in the open unit disk: z (+) a.

    z and a are complex 1D arrays of the same shape.  The result is
    also a complex array; |z (+)| < 1 whenever |z| < 1 and |a| < 1.
    This is the standard Poincare-disk addition used in Escher's
    Circle Limit etchings and in the M.C. Escher / H.S.M. Coxeter
    correspondence on hyperbolic tessellations.
    """
    za = z * a
    num = z + a
    den = 1.0 + za.conj()
    return num / den


def _tile_reflect(z: np.ndarray, p: int, q: int) -> np.ndarray:
    """Reflect z (a complex array, |z|<1) across one side of the
    (p, q) hyperbolic triangle with one vertex at the origin.

    The triangle has interior angle pi/p at the origin and pi/q at
    each of the two finite vertices on the unit circle.  The
    generators of the (p, q, 2) triangle group are
        R = rotation by 2*pi/p
        S = inversion across the line through 0 and e^{i theta}
    with theta = pi/p + pi/q - pi/2 (so the third angle is pi/2 at
    the finite vertices).  Inversion across a line through the origin
    is just complex conjugation after rotation:
        S(z) = exp(2 i theta) * conj(z)
    A "coset expansion" tile: reflect z against S, then again against
    S R^k until we find a representative inside the fundamental
    triangle; report how many reflections it took.
    """
    theta = np.pi / p + np.pi / q - np.pi / 2.0
    e2i = np.exp(2j * theta)
    n_reflect = np.zeros(z.shape, dtype=np.int32)
    # Naive coset expansion: up to 2*p iterations
    max_iter = 2 * p
    inside = np.ones(z.shape, dtype=bool)
    for k in range(max_iter):
        if not np.any(inside):
            break
        z_work = z[inside]
        rot = np.exp(2j * np.pi * k / p)
        z2 = e2i * np.conj(z_work)  # S
        z3 = z2 * rot  # then rotate
        # Are we inside the wedge [0, 2*pi/p] with theta < arg < pi-theta?
        arg = np.angle(z3)
        ok = (arg >= -theta) & (arg <= 2 * np.pi / p + theta)
        if np.any(ok):
            z[inside] = np.where(ok, z3, z[inside])
            n_reflect[inside] += np.where(ok, 1, 0)
            inside[inside] = ~ok
    return n_reflect