Geometría Hiperbólica: Disco de Poincaré

En este post les muestro cómo generar teselaciones hiperbólicas regulares del disco de Poincaré. El código de ejemplo está escrito en C# y utilicé Windows Presentation Foundation (WPF) por su facilidad de generar gráficos vectoriales, hice uso del .Net Framewoek 4 por su soporte para números complejos. Por supuesto, no es difícil transportar todo esto a Java, C++ o cualquier otro lenguaje.

Teselacion {5,4}
Teselación Hiperbólica {5,4} generada en C# y WPF

Modelo de Poincaré para la Geometría Hiperbólica

Se pueden hacer muchos paralelos entre la geometría euclidiana y la hiperbólica. En el modelo de Poincaré, todo el espacio hiperbólico \mathbb{H}^2 está representado dentro de un disco de radio uno.  El borde del disco representa el infinito.  Dentro de este disco se cumplen los postulados de Euclides exceptuando el 5to (el de las paralelas):

  1. Se puede trazar una línea recta que pase por dos puntos.
  2. Se puede prolongar una línea recta indefinidamente a partir de una recta finita.
  3. Se puede trazar una circunferencia con centro y radio dado.
  4. Todos los ángulos rectos son iguales.
  5. Si una línea recta que corta a otras dos rectas forma de un mismo lado con ellas ángulos interiores cuya suma es menor que dos rectos, las dos últimas rectas prolongada indefinidamente se cortan del lado en que la suma de los ángulos es menor que dos rectos.

En \mathbb{H}^2 la suma de los ángulos internos de un triángulo es menor a 180º. Mas sorprendente aún, dos rectas con direcciones distintas pueden ser paralelas. El modelo de Poincaré permite visualizar estos aspectos de la geometría hiperbólica, pero al estar todo el espacio dentro de un disco, las líneas que en realidad son rectas son percibidas como curvas (de ahí que se les llame “Geodésicas”). Y la métrica que nos permite medir distancias dentro del disco de Poincaré no es euclidiana.

Implementación del modelo de Poincaré

Para la implementación de la geometría hiperbólica en el modelo de Poincaré, son útiles los siguientes conceptos:

  • El espacio hiperbólico \mathbb{H}^2 es un disco de radio uno centrado en el origen del plano euclidiano \mathbb{R}^2 llamado disco de Poincaré
  • Los puntos en el espacio hiperbólico \mathbb{H}^2 son puntos en el plano euclidiano que están dentro del disco de Poincaré
  • Las líneas que pasan por dos puntos en \mathbb{H}^2 son círculos euclidianos que pasan por dos puntos en el disco y son ortogonales al disco de Poincaré
  • Las líneas que pasan por el origen (i.e., el centro del disco de Poincaré) son círculos de radio infinito, es decir, son líneas euclidianas.
  • Los ángulos son euclidianos, la medida del ángulo que se forma entre dos geodésicas (líneas hiperbólicas) es el angulo que forman las tangentes de los circulos en el punto en que éstos se intersecan
  • La inversión de un punto en el circulo es una isometría (preserva ángulos y distancias) y es interpretado como la reflexión de un punto en una linea hiperbólica

Construcciones elementales

Punto en \mathbb{H}^2 (también representa puntos en \mathbb{R}^2)

public struct Point
{
     public double X { get; set; }
     public double Y { get; set; }
}

Geodésica (línea) que pasa por los puntos A y B en \mathbb{H}^2

    //Linea hiperbolica
    public class Geodesic
    {
        //Puntos por los que pasa la linea
        public Point A { get; set; }
        public Point B { get; set; }
        //Linea hiperbolica cuando NO pasa por el origen
        public Circle C { get; set; }
        //Linea hiperbolica cuando pasa por el origen
        public Line L { get; set; }
        public bool IsLine { get { return C.HasInfiniteRadius; } }
    }

Inversión del punto p en el circulo C (se puede interpretar como la reflexión del punto p en la geodésica C).

        // Invierte el punto P en el circulo C
        Point PointInversion(Point P, Circle C)
        {
            if (C.HasInfiniteRadius)
                throw new Exception("Inversion en circulos de radio infinito no soportadas");
            double k2 = C.Radius * C.Radius;
            return C.Center + ((P - C.Center) * k2) / (P - C.Center).LengthSquared;
        }
        // Invierte el punto P en la geodesica G
        Point PointInversion(Point P, Geodesic G)
        {
            if (G.IsLine)
                //refleja el punto en la linea euclidiana
                return Geo.L.Reflect(P);
            return PointInversion(P, Geo.C);
        }

Creación de una geodésica que pasa por los puntos P0 y P1.

        //Crea una linea hiperbolica que pasa por puntos p0 y p1
        public Geodesic GetGeodesicLine(Point p0, Point p1)
        {
            //Invierte uno de los puntos con respecto al disco de Poincare
            Point pointInv = PointInversion(p1, PoincareDisk);
            //Calcula el circulo usando tres puntos
            return new Geodesic(A: p0, B: p1, C: new Circle(p0, p1, pointInv));
        }

Estas construcciones son necesarias sólo si nos restringimos a usar el álgebra vectorial. Si en cambio nos decidiéramos a utilizar Álgebra Geométrica Conforme tendríamos todo lo necesario a nuestra disposición, es decir, verdaderas geodésicas (círculos que pueden tener el radio infinito!), inversiones en el circulo, intersecciones entre primitivas, transformaciones conformes (sin la necesidad de usar números complejos ni transformaciones de Möbius, como veremos después), todo integrado en el sistema algebraico!. Además el álgebra geométrica no está ligada limitada al plano 2D (las transformaciones de Möbius si lo están) y se puede extender la geometría a las tres dimensiones sin mucho esfuerzo. En futuros posts estaré mostrando como hacer geometría hiperbólica y euclidiana con Álgebra Geométrica Conforme.

Construcción de Teselaciones Regulares

Una teselación regular es un cubrimiento del espacio con polígonos regulares (polígonos de lados iguales). El plano hiperbólico tiene un número infinito de posibles teselaciones planas, mientras que el plano euclidiano tiene sólo 3 teselaciones posibles (triangulo equilatero, cuadrado y hexágono). El espacio de las posibles teselaciones está dado por dos variables: el número de lados del polígono regular p y el número de polígonos incidentes en un vértice q. El par {p,q} determina la teselación, por ejemplo el siguiente gráfico muestra las teselaciones hiperbólicas {6,4}, {4,6}, {7,3} y {3,7}:

Teselaciones y duales
Arriba: Teselación {6,4} y su dual {4,6} Abajo: Teselación {7,3} y su dual {3,7}

En este sitio se puede ver la tabla de posibles combinaciones y en este otro algunas teselaciones y sus duales. La formula general para saber si una teselación es hiperbólica, esférica o euclidiana es:

(p-2)(q-2) = \left \{ \begin{array}{ll} > 4, & \mbox{Hiperb\'{o}lico} \\ = 4, & \mbox{Euclidiano} \\ < 4, & \mbox{Esf\'{e}rico} \\ \end{array} \right.

La clave para construir estas teselaciones está en crear el polígono central, llamado también región fundamental. Una vez que este polígono está definido, la teselación consiste básicamente en aplicarle transformaciones isométricas, en particular reflejarlo en sus p lados.

Construir el polígono central es bastante simple, en esta página se detallan los pasos:

  1. Calcular el radio del circulo que circunscribe a la región fundamental
  2. Calcular las posiciones de los vértices del polígono central en la circunferencia de la región fundamental
  3. Reflejar el polígono central a través de sus lados, generando N nuevos polígonos.
  4. Para cada lado de un polígono existente que sólo es adyacente a un polígono, reflejar el polígono a través del lado para continuar la teselación
  5. Repetir el paso 4 ad-infititum

Según el paper de Coexeter “The Trigonometry of Hyperbolic Tessellations“, podemos hallar fácilmente un vértice del polígono central usando trigonometría y semejanza de triángulos. Primero tenemos que hallar el centro O y el radio \overline{OC} del circulo que se muestra a continuación:

Donde el ángulo \angle BAC mide \pi/p, el ángulo \angle ABC mide \pi/q y el ángulo \angle ACB mide \pi/2. Aplicando la siguiente formulación (recuerda que en matemáticas el disco de Poincaré tiene radio 1 y está centrado en el origen (0,0)):

\overline{OC} = \left[ \frac{\cos^2 \pi/q} {\sin^2 \pi/p} - 1 \right]^{-1/2} \\ \overline{AO} = \left[ 1 - \frac{\sin^2 \pi/p} {\cos^2 \pi/q} \right]^{-1/2}

Donde el centro del circulo O = (\|\overline{AO}\|, 0) y el radio r = \overline{OC}. Finalmente el vértice B puede ser hallado intersecando el rayo que parte de A con ángulo \pi/p con el círculo. Una vez hallado este vértice, ya podemos hallar el radio del circulo que circunscribe a la región fundamental \|\overline{AB}\| y también las posiciones de los vértices del polígono central en la circunferencia de la región fundamental simplemente añadiendo un vértice en la circunferencia cada 2 \pi/p radianes.

En C# esto se puede codificar sigue:

    public class Edge
    {
        public int v0 { get; set; }
        public int v1 { get; set; }
        public List<int> Neighbors { get; set; }
    }
    public class Poly
    {
        public int level { get; set; }
        public int center { get; set; }
        public List<Edge> edges { get; set; }
    }
        Poly CreateCenterPolygon(int p, int q)
        {
            double s = Math.Sin(Math.PI / p);
            double c = Math.Cos(Math.PI / q);
            // Computo OC y AO segun el articulo de Coexeter
            double OC = 1.0 / Math.Sqrt(((c * c) / (s * s)) - 1.0);
            double AO = 1.0 / Math.Sqrt(1.0 - ((s * s) / (c * c)));

            // Creo el circulo con centro = (AO,0) y radio = OC
            var circle = new Circle(center: new Point(AO, 0), radius: OC);

            // Computo el punto B hallando la interseccion del circulo con
            // la linea que parte desde A en direccion = ( cos(pi/p), sin(pi/p) )
            var lineAB = new Line(o: new Point(0, 0), dir: new Vector(Math.Cos(Math.PI / p), Math.Sin(Math.PI / p)));
            Point B = Geometry.LineCircleIntersection(lineAB, circle);

            // Computo las posiciones de los vertices del polígono central
            // en la circunferencia de la región fundamental añadiendo un
            // vertice en la circunferencia cada 2*pi/p radianes
            double dist = (B - lineAB.O).Length;
            for (int i = 0; i < p; ++i)
            {
                double alpha = (double)i / p * 2.0 * Math.PI;
                Point X = Complex2Vector(dist * new Complex(Math.Cos(alpha), Math.Sin(alpha)));
                Points.Add(X);
            }
            // Añado el punto A a los vertices
            Points.Add(Complex2Vector(new Complex(0, 0))); //centro A

            // Defino el poligono central
            Poly poly = new Poly();
            for (int i = 0; i < p; ++i)
                poly.edges.Add(new Edge(i, (i + 1) % p));
            poly.level = 1;
            poly.center = p;
            // para computar la teselacion dual cada lado guarda sus
            // poligonos adyacentes
            foreach (var e in poly.edges)
                e.Neighbors.Add(poly.center); 

            //retorna el poligono
            return poly;
        }

El metodo Complex2Vector simplemente escala y traslada los vértices desde disco de Poincaré matemático (radio = 1 y centro en (0,0)) hacia el disco de Poincaré en coordenadas de pantalla (representado en la propiedad Disk):

        Point Complex2Vector(Complex c)
        {
            return new Point(c.Real * Disk.Radius + Disk.Center.X, c.Imaginary * Disk.Radius + Disk.Center.Y);
        }

Finalmente queda aplicar reflexiones isométricas al polígono central para poder cubrir todo el plano. Mi estrategia es de fuerza bruta, primero itero por los lados del polígono central y lo reflejo en cada uno de ellos, luego repito lo mismo con los nuevos polígonos reflejados. El algoritmo avanza capa por capa para generar un cubrimiento más uniforme. El código C# es como sigue:

        void ApplyIsometries(Poly centralPoly, int maxLevel)
        {
            //Este metodo es de fuerza bruta.
            HashSet<Edge> edges = new HashSet<Edge>();
            Queue<Poly> Q = new Queue<Poly>();
            Q.Enqueue(centralPoly);
            Polygons.Add(centralPoly);
	        while(!(Q.Count == 0))
	        {
		        Poly current = Q.Dequeue();

		        if(current.level >= maxLevel )
			        continue;

                foreach (var edge in current.edges)
		        {
                    if (!edges.Contains(edge))
			        {
                        // procesa el lado solo una vez
                        edges.Add(edge);
				        int oldPointSize = Points.Count;
                        // Refleja el poligono en el lado
                        Poly poly = ReflectPoly(current, edge);
                        if (oldPointSize != Points.Count || !Polygons.Contains(poly))
				        {
                            // si el poligono no existe lo añado a la lista
                            Polygons.Add(poly);
                            Q.Enqueue(poly);
                        }
                        // para computar la teselacion dual cada lado guarda sus
                        // poligonos adyacentes
                        edge.Neighbors.Add(poly.center);
                    }
		        }
	        }
        }
        private Poly ReflectPoly(Poly current, Edge inversionEdge)
        {
            var geodesic = GetGeodesicLine(inversionEdge);
            List<int> indices = new List<int>();
            foreach(var edge in current.edges)
            {
                if (edge.v0 == inversionEdge.v0 || edge.v0 == inversionEdge.v1)
                    indices.Add(edge.v0);
                else
                    indices.Add(ReflectPoint(geodesic, edge.v0));
            }

            Poly poly = new Poly();
            for (int i = 0; i < indices.Count; ++i)
                poly.edges.Add(new Edge(indices[i], indices[(i + 1) % indices.Count]));
            poly.level = current.level + 1;
            poly.center = ReflectPoint(geodesic, current.center);
            foreach (var edge in poly.edges)
                edge.Neighbors.Add(poly.center);
            return poly;
        }

El método ReflectPoint básicamente realiza una inversión en el circulo para reflejar el vértice, añade el vértice reflejado a la lista Points (si no existía ya) y retorna el indice.

        private int ReflectPoint(Geodesic geodesic, int pointIndex)
        {
            Point newP = PointInversion(Points[pointIndex], geodesic);

            var pt = Points.FindIndex(p => Math.Abs(p.X - newP.X) < 1e-6 && Math.Abs(p.Y - newP.Y) < 1e-6);
            if (pt >= 0)
                return pt;
            var index = Points.Count;
            Points.Add(newP);
            return index;
        }

Con esto se ha terminado la generación de la teselación regular. La clase que encapsula toda la magia se ve algo así:

    public class Tessellation
    {
        public Circle Disk { get; set; }
        public List<Point> Points { get; set; }
        public SortedSet<Poly> Polygons { get; set; }

        public Tessellation(Circle c)
        {
            Disk = c;
            Points = new List<Point>();
            Polygons = new SortedSet<Poly>();
        }

        public void Tessellate(int p, int q, int maxLevel)
        {
            var centralPoly = CreateCenterPolygon(p, q);
            ApplyIsometries(centralPoly, maxLevel);
        }
        ...
    }

Transformaciones de Möbius


Rotación y traslación hiperbólica con transformaciones de Möbius

La forma usual de representar rotaciones, traslaciones, inversiones, y otras transformaciones conformales en el espacio hiperbólico 2D es usando transformaciones de Möbius. Para usar esta técnica tenemos que representar un vector como un número complejo z \in \mathbb{C}. La transformación de Möbius es una función de la forma:

f(z) = \frac{a z + b}{c z + d}

donde los números a, b, c y d son complejos. Podemos expresar esta transformación lineal como una matriz 2×2 utilizando coordenadas homogéneas de la siguiente forma:

Tz = \left[ \begin{array}{cc} a & b \\ c & c \\ \end{array} \right] \left[ \begin{array}{c} z \\ 1 \\ \end{array} \right] = \left[ \begin{array}{c} a z + b \\ c z + d \\ \end{array} \right] = \left[ \begin{array}{c} z' \\ w \\ \end{array} \right]

El resultado final se obtiene al normalizar el resultado tal que w = 1, es decir:

Tz = \frac{1}{w} \left[ \begin{array}{c} z' \\ w \\ \end{array} \right] = \left[ \begin{array}{c} z'/ w \\ 1 \\ \end{array} \right]

De esta forma, toda transformación de Möbius se puede expresar como una matriz, donde la multiplicación de matrices equivale a la composición de transformaciones. De igual forma la inversa de una transformación de Möbius es la inversa de la matriz.

En el espacio hiperbólico de Poincaré, las traslaciones en el origen se pueden expresar de la siguiente forma:

T(z) = \left[ \begin{array}{cc} 1 & b \\ \bar b & 1 \\ \end{array} \right] \left[ \begin{array}{c} z \\ 1 \\ \end{array} \right] = \left[ \begin{array}{c} z + b \\ \bar b z + 1 \\ \end{array} \right]

Donde b es un vector cuya magnitud es menor a 1 y \bar b es el conjugado complejo de b. Las rotaciones en el origen se pueden expresar de la siguiente forma:

R(z) = \left[ \begin{array}{cc} e^{i \theta} & 0 \\ 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} z \\ 1 \\ \end{array} \right] = \left[ \begin{array}{c} e^{i \theta} z \\ 1 \\ \end{array} \right]

Donde e^{i \theta} está definido por la formula de Euler \cos \theta + i \sin \theta. La transformación rígida general dada por T(R(z)) se puede expresar como:

T(R(z)) = \left[ \begin{array}{cc} 1 & b \\ \bar b & 1 \\ \end{array} \right] \left[ \begin{array}{cc} e^{i \theta} & 0 \\ 0 & 1 \\ \end{array} \right] \left[ \begin{array}{c} z \\ 1 \\ \end{array} \right] = \left[ \begin{array}{cc} e^{i \theta} & b \\ \bar b e^{i \theta} & 1 \\ \end{array} \right] \left[ \begin{array}{c} z \\ 1 \\ \end{array} \right]

Implementación de Transformaciones de Möbius

La siguiente implementación está basada en la clase Complex del .Net Framework 4.

    public class Mobius
    {
        ComplexMatrixNxM M;

        public Complex A { get { return M[0, 0]; } set { M[0, 0] = value; } }
        public Complex B { get { return M[0, 1]; } set { M[0, 1] = value; } }
        public Complex C { get { return M[1, 0]; } set { M[1, 0] = value; } }
        public Complex D { get { return M[1, 1]; } set { M[1, 1] = value; } }

        public Mobius(Complex a, Complex b, Complex c, Complex d)
        {
            M = new ComplexMatrixNxM(new Complex[,] { { a, b }, { c, d } });
        }

        public Mobius(ComplexMatrixNxM M)
        {
            Debug.Assert(M != null);
            Debug.Assert(M.cols == 2 && M.rows == 2);
            this.M = M;
        }

        public Complex Apply(Complex z)
        {
            return (A * z + B) / (C * z + D);
        }

        //Rotacion en el origen
        public static Mobius Rotation(double angle)
        {
            return new Mobius(a: Exp(angle), b: Complex.Zero, c: Complex.Zero, d: Complex.One);
        }

        //Traslacion en el origen
        public static Mobius Translation(double x1, double y1)
        {
            return new Mobius(a: Complex.One, b: new Complex(x1, y1), c: Complex.Conjugate(new Complex(x1, y1)), d: Complex.One);
        }

        // Composicion de transformaciones: (F1*F0)(z) = F1(F0(z))
        public static Mobius operator *(Mobius F1, Mobius F0)
        {
            //Multiplicacion de matrices
            return new Mobius(M: F1.M * F0.M);
        }

        // formula de Euler
        public static Complex Exp(double x)
        {
            return new Complex(Math.Cos(x), Math.Sin(x));
        }
    }

Código Fuente

Espero que hayan disfrutado el post. Les dejo el código fuente para que se diviertan.

Advertisements

3 thoughts on “Geometría Hiperbólica: Disco de Poincaré

  1. María says:

    En la entrada “Geometría Hiperbólica: Disco de Poincaré” se ofrece el código en C# para el cálculo de Teselación hiperbólica. Hay un vínculo hacia dropbox, que no redirige a ningún sitio válido. No sé si es posible reestablecer ese vínculo de forma correcta para que vuelva a estar disponible ese código. Gracias!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s