Generating spherical and hyperbolic tilings in GLSL

Posted on Sat 06 January 2018 in misc

In my previous post, I explored generating hyperbolic, euclidean and spherical tilings using circle inversions. While it produced interesting pictures and animations, I wasn't able to derive a formula to generate particular tilings. There's another way to generate tilings that admits a beautifully simple unification of the three geometries. We will need a bit more machinery, but the result is interesting and approaches the problem in a different way. This post is based nearly entirely on Knighty's shader with some help from notes the author left on fractalforums.

The code below is in the form of fragment shaders. For a shader primer, check out the Book of Shaders. The important thing is that to write a fragment shader, you build a function that takes a location (in our case, of a pixel) and returns a color. All the code here is "live", though changes you make in each box will not cascade down to the subsequent ones.


Triangles in different geometries

To begin, we need to define the triangle that we are going to tile. A triangle can be defined by the angles made at each vertex. In hyperbolic and spherical geometry, three angles are sufficient to define the triangle uniquely (you can't scale triangles without the internal angles also changing). A triangle can tile its geometry if (and only if) all three angles divide \(2\pi\). In other words, around each vertex, you can fit a whole number of triangle copies. The usual way to define a triangle is therefore to pick \(p,q,r\) and set the internal angles to be \(\frac{\pi}{p},\frac{\pi}{q}\) and \(\frac{\pi}{r}\). To tell if \((p,q,r)\) describes a spherical, planar, or hyperbolic triangle, we compare the sum of the interior angles to \(\pi\): if it's less, it's hyperbolic; if greater, it's spherical. Since they all divide \(\pi\), you can simply compare \(\frac{1}{p} + \frac{1}{q} + \frac{1}{r}\) to \(1\).

The next thing we need to define is our model of each geometry. Euclidean and spherical geometry both have natural models: the plane and the sphere. In the previous post, we mentioned the hyperboloid model on the way to working in the Poincaré model, but here we are going to use the hyperboloid model directly. The sphere and plane are both familiar, but we have to inspect the hyperboloid model more closely to understand what's actually going on.


Hyperboloid

The hyperboloid model lives on an infinitely large three-dimensional bowl shape. What you can't see in the diagram is that it doesn't live inside ordinary space at all. This space has a different metric, which means lengths and angles aren't measured quite the same. In ordinary 3D space, the dot product of two vectors \(\vec{a} \cdot \vec{b}\), \(\vec{a}=\left\langle a_0,a_1,a_2\right\rangle,\vec{b}=\left\langle b_0,b_1,b_2\right \rangle\) is \(a_0b_0+a_1b_1+a_2b_2\).

We can use the dot product to define the squared-magnitude (aka length squared) of a vector \(\vec{a}\) :

$$\vec{a} \cdot \vec{a} = a_0^2+a_1^2+a_1^2 = |\vec{a}|^2$$

The hyperboloid model lives in a space where a different dot product applies:

$$\vec{a} \cdot_{\tiny{H}} \vec{b} = a_0b_0+a_1b_1-a_2b_2$$

If we use the dot product to define the magnitude-squared of \(\vec{a}\), we get a different magnitude than in ordinary space. In other words, distances behave differently here.

$$|\vec{a}|^2 = a_0^2+a_1^2-a_1^2$$

Hyperbola

The equation for this hyperboloid is \(x^2 + y^2 - z^2 = 1\). Alternatively, you can define it as vectors \(\vec{a}\) such that \(\vec{a} \cdot_{\tiny{H}} \vec{a} = 1\). In the ordinary 3-D Euclidean metric, \(\vec{a} \cdot \vec{a} = 1\) is a unit sphere. So, this hyperboloid is a sort of "hyperbolic sphere"- it's this space's equivalent of a sphere, the set of all points a unit distance away from the origin. We can draw a 2D version of this easily. Notice that with space = 0, it "degenerates" into two flat planes.

It's important that we are actually flipping the surface so that the two halves of the hyperboloid open vertically up and down, and the Euclidean plane(s) are also in line with the \(xy\) plane. It's points on these surfaces that behave like we want.

In the hyperboloid model of hyperbolic geometry, lines are the intersection of the hyperboloid with planes through the origin ("cutting planes"), which take the form of hyperbolas. These are akin to the great circles, which are lines in spherical geometry, and are also intersections of the sphere with cutting planes. We can also identify lines in Euclidean geometry as intersections of cutting planes with a plane at \(z=1\).

It's obvious that if you take a plane passing through the origin, it cuts the sphere in two identical pieces, and reflecting the sphere across the plane is an isometry: it doesn't change the sphere. This is also true for cutting planes and the hyperboloid... if you reflect it using the hyperbolic dot product, rather than the Euclidean.

First, we identify the plane with the normal vector \(\vec{p}\), which has a unit norm \(\vec{p} \cdot \vec{p} = 1\). To reflect a vector \(\vec{a}\) across a plane with unit normal \(\vec{p}\), we can perform the following operation:

$$\vec{a}' = \vec{a} - 2(\vec{a} \cdot \vec{p})\vec{p}$$

It's easiest to see why this works in the Euclidean case. \(\vec{a} \cdot \vec{p}\) is the distance from \(\vec{a}\) (interpreted as a point) to the plane. \(\vec{a} - (\vec{a} \cdot \vec{p})\vec{p}\) lies on the plane; subtracting again, and you're the same distance away on the other side of the plane. Importantly, this works exactly the same way with points on the vertical hyperboloid and the hyperbolic dot product. Reflecting the hyperboloid across a plane like this is an isometry, just like reflecting a sphere is.

However, to make this work for euclidean tilings we are going to rewrite it in terms of the dot product only. To start, let's try it with the hyperbolic case. If we negate the \(z\) coordinate of \(\vec{p}\), we can build a transformation that reflects properly (negating \(z\) is "safe" because between the two arms of the hyperboloid, the problem is symmetrical).

$$\begin{align*} \vec{p} &= \langle p_0,p_1,p_2 \rangle \\ \vec{q} &= \langle p_0,p_1,-p_2 \rangle \\ \vec{a}' &= \vec{a} - 2(\vec{a} \cdot_{\tiny{H}} \vec{p}) \vec{p} \\ \vec{a}' &= \vec{a} - 2( a_0p_0+a_1p_1-a_2p_2 ) \vec{p} \\ \vec{a}' &= \vec{a} - 2( \vec{a} \cdot \vec{q} ) \vec{p} \\ \vec{a}' &= \vec{a} - 2( \vec{a} \cdot \vec{q} ) \vec{q} \langle1,1,-1\rangle \end{align*}$$

So for \(k \in {1,-1}\) it's clear that \(\vec{a} - 2( \vec{a} \cdot \vec{p} ) \vec{p} \langle1,1,k\rangle\) is the right transformation. But now if we set \(k=0\), it produces a reflection across the intersection of the cutting plane with the plane at \(z=1\) (it definitely maps the plane to itself, because \(z\) is unchanged). Strictly speaking, in the plane it is not quite an isometry: distances end up being scaled as well.

Here it is in action in two dimensions. Switch to the ordinary space by setting space = 1.. See what it does to the red region in each space as you adjust p. Notice that the shape of the surface (circle, hyperboloid or plane) is preserved by the transformation.

So. Now we know how to reflect in our geometries. Now, we need to work out what planes to reflect across.

If two cutting planes intersect with angle \(\theta\), the angle the associated great circles make on the sphere is clearly also \(\theta\), because the planes are always perpendicular to the sphere. This is also true for the angle between two planes and the angle between their associated hyperbolas on the hyperboloid: the cutting planes are actually perpendicular to the hyperboloid, if you measure using the hyperbolic dot product. You can define angles between two cutting planes by the dot product of their respective normal vectors:

$$\cos{\theta} = |\vec{a}||\vec{b}| \vec{a} \cdot \vec{b}$$

This also applies with the hyperbolic dot product, but angles (like distances) stop being exactly what they look like.

So, the problem of finding triangles with a particular \((p,q,r)\) is the same as finding three vectors \(\vec{a}, \vec{b}, \vec{c}\) that define cutting planes that make the correct angles with each other, measured in the correct metric. So:

$$\begin{align*} \vec{a} \cdot \vec{b} & = -\cos{\frac{\pi}{p}}\\ \vec{a} \cdot \vec{c} & = -\cos{\frac{\pi}{q}}\\ \vec{c} \cdot \vec{b} & = -\cos{\frac{\pi}{r}} \end{align*}$$

Choose \(\vec{a}\) and \(\vec{b}\) to lie in the \(xy\) plane making angle \(\frac{\pi}{p}\):

$$\begin{align*} \vec{a} &= \left\langle1,0,0\right\rangle\\ \vec{b} &= \left\langle -\cos{\frac{\pi}{p}},\sin{\frac{\pi}{p}},0\right\rangle \end{align*}$$

Notice that \(a_2=0\) so, regardless of metric,

$$\begin{align*} \vec{a} \cdot \vec{c} &= a_0 c_0 +a_1c_1\\ &= c_0 \end{align*}$$

Next we solve for \(c_1\). \(b_2=0\) so \(\vec{b} \cdot \vec{c}\) is simply (again, regardless of metric):

$$\begin{align*} \vec{b} \cdot \vec{c} &= b_0 c_0 +b_1c_1\\ c_1&=\frac{\vec{b}\cdot \vec{c} -b_0c_0}{b_1}\\ \end{align*}$$

\(\vec{c}\) must be a unit vector, so \(\vec{c} \cdot \vec{c} = 1\). But, we don't know what metric we're in, so we parametrize it with \(k\in\{-1,0,1\}\).

$$\begin{align*} 1&=c_0^2+c_1^2+kc_2^2\\ k{c_2}^2 &= 1-c_0^2-c_1^2\\ \left|k{c_2}^2\right| &= \left|1-c_0^2-c_1^2\right|\\ c_2 &= \sqrt{\left|1-c_0^2-c_1^2\right|} \end{align*}$$

We picked the positive square root for \(c_2\), which is an arbitrary choice in the Euclidean case but in the hyperbolic case, there are actually two arms of the hyperbola (one pointing up, the other down). We want to build the tiling in the upper one and this choice does that for us.

The parameter \(k\) has fallen away, so we actually have a formula that works for either hyperbolic or spherical triangles. For Euclidean ones, \(k=0\) so \(c_2\) can be any constant, and it just acts like a scaling factor.

We can now substitute back in the dot products, to get the final values for the components of \(\vec{c}\)

$$\begin{align*} c_0 &= -\cos{\frac{\pi}{q}}\\ c_1 &= \frac{ -\cos{\frac{\pi}{q}} -\cos{\frac{\pi}{p}}\cos{\frac{\pi}{r}}}{\sin{\frac{\pi}{p}}}\\ c_2 &= \sqrt{\left|\sin^2{\frac{\pi}{r}}-c_1^2\right|}\\ \end{align*}$$

The last thing we need is an equation to project points in the plane to points on the hyperbola or sphere. There are a few choices, but we're going to use a stereographic projection. The idea is to draw a line between \((0,0,-1)\) and a point on the surface, and map that point to the intersection of line with the plane at \(z=0\).

We are going to work in two dimensions to make the algebra much easier. The surfaces have rotational symmetry around the \(z\) axis, so we can work in \((r,z)\) coordinates. \(r\) is the distance away from the \(z\) axis and corresponds to the \(x\) dimension.

$$r,z \mapsto \left( R,Z \right) = \left( \frac{r}{1+z},0\right) $$

We need to invert the transformation to take the plane to the surfaces. We can define the surfaces like this:

$$z^2 + kr^2 = 1, k \in \left(-1,0,1\right)$$

So, we start with \(\frac{r}{1+z}\), square it, substitute in and solve for \(z\):

$$\begin{align*} R^2 &= \frac{r^2}{(1+z)^2}\\ &= \frac{1}{k} \frac{1-z^2}{(1+z)^2} = k \frac{(1+z)(1-z)}{(1+z)^2}\\ &= k \frac{1-z}{1+z} \\ R^2 (1+z) &= k (1-z) \\ R^2 + R^2 z &= k - kz \\ R^2 z + kz &= k - R^2 \\ z &= \frac{k - R^2}{k + R^2} = \frac{k - R^2}{k + R^2} \frac{k}{k}\\ \therefore z &= \frac{1 - kR^2}{1 + kR^2} \end{align*}$$

And we can now solve for \(r\):

$$\begin{align*} R &= \frac{r}{1+z} \\ r &= \left(1+z\right) R \\ r &= \left(1+\frac{1-kR^2}{1+kR^2}\right) R \\ r &= \left(\frac{1+kR^2}{1+kR^2}+\frac{1-kR^2}{1+kR^2}\right) R \\ \therefore r &= \frac{2R}{1+kR^2} \end{align*}$$

So the transformation is

$$ R \mapsto \left(\frac{2R}{1+kR^2}, \frac{1 - kR^2}{1 + kR^2},\right)$$

We can pull this into two dimensions easily. \(R^2 = X^2+Y^2\), and doubling \(R\) doubles both \(X\) and \(Y\), so the final function looks like:

$$ \left(X,Y\right) \mapsto \left(\frac{2X}{1+k\left(X^2+Y^2\right)}, \frac{2Y}{1+k\left(X^2+Y^2\right)}, \frac{1-k\left(X^2+Y^2\right)}{1+k\left(X^2+Y^2\right)}\right)$$

Pulling this all together, and out pops the tilings we want. Adjust p,q,r to see how they behave.

In the hyperbolic case, we can investigate other projections, such as just pulling the point in the plane up to the hyperboloid:

We can also apply a hyperbolic rotation to the hyperbolic tiling, which acts like a translation along the hyperboloid.

That's all for now. This is a lot more laborious than just inverting across circles, but being able to mostly avoid trigonometry has an appeal all its own, I think!