Reflex reboot
{3,3,3} 5-cell (simplex)
{4,3,3} 8-cell (hypercube)
{3,3,4} 16-cell (cross polytope)
{3,4,3} 24-cell
{5,3,3} 120-cell
{3,3,5} 600-cell
Explanation
Back in 2009 I posted a short video Reflex preview of truncated 4D polytopes. Today I revisited that old code, reimplementing the core algorithms in OpenGL Shader Language (GLSL) using the FragM environment.
3D polytopes are also known as polyhedra. They can be defined by their Schlaefli symbol, which looks like {4,3} for a cube. This can be understood as having faces made up of squares ({4}), arranged in triangles around each vertex ({3}). This notation extends to 4D, where a hypercube (aka tesseract) has symbol {4,3,3}, made up of cubes {4,3} with {3} around each vertex.
These symbols are all well and good, but to render pictures you need concrete measurements: angles, lengths, etc. In H.S.M. Coxeter's book Regular Polytopes, there is the solution in the form of some recursive equations. Luckily the recursion depth is bounded by the number of dimensions, which is small (either 3 or 4 in my case, though in principle you can go higher). This is important for implementation in GLSL, which disallows recursion. In any case, GLSL doesn't have dynamic length lists, so different functions are needed for different length symbols.
I don't really understand the maths behind it (I don't think I did even when I first wrote the code in Haskell in 2009), but I can describe the implementation. The goal is to find 4 vectors, which are normal to the fundamental region of the tiling of the (hyper)sphere. Given a point anywhere on the sphere, repeated reflections through the planes of the fundamental region can eventually (if you pick the right ones) get you into the fundamental region. Then you can do some signed distance field things to put shapes there, which will be tiled around the whole space when the algorithm is completed.
Starting from the Schlaefli symbol {p,q,r} (doing 4D because it has some tricksy bits, 3D is largely similar), the main task is to find the radii of the sub polytopes ({p,q}, {p}, {q,r}, etc). This is because these radii can be used to calculate the angles of the planes of the fundamental region, using trigonometry. The recursive formula starts here:
float radius(int j, ivec3 p)
{
return sqrt(radius2_j(j, p));
}
Here j will range from 0 to 3 inclusive, and the vector is {p,q,r}. Then radius2_j() evaluates using the radii squared of the relevant subpolytopes, according to j. I think this is an application of Pythagoras' theorem.
float radius2_j(int j, ivec3 p)
{
switch (j)
{
case 0: return radius2_0(p);
case 1: return radius2_0(p) - radius2_0();
case 2: return radius2_0(p) - radius2_0(p.x);
case 3: return radius2_0(p) - radius2_0(p.xy);
}
return 0.0;
}
The function radius2_0() is overloaded for different length inputs, from 0 to 3:
float radius2_0() { return 1.0; }
float radius2_0(int p) { return delta() / delta(p); }
float radius2_0(ivec2 p) { return delta(p.y) / delta(p); }
float radius2_0(ivec3 p) { return delta(p.yz) / delta(p); }
Here it starts to get mysterious, the delta funtion uses trigonometry to find the lengths. I don't know how/why this works, I copy/pasted from the book. Note that it looks recursive at first glance, but in fact each delta calls a different delta(s) with strictly shorter input vectors, so it's just a non-recursive chain of function calls.
float delta()
{
return 1.0;
}
float delta(int p)
{
float s = sin(pi / float(p));
return s * s;
}
float delta(ivec2 p)
{
float c = cos(pi / float(p.x));
return delta(p.y) - delta() * c * c;
}
float delta(ivec3 p)
{
float c = cos(pi / float(p.x));
return delta(p.yz) - delta(p.z) * c * c;
}
Now comes the core function to find the fundamental region: the cosines of the angles are found by ratios of successive radii, the sines are found by Pythagoras' theorem, 3 rotation matrices are constructed, then one axis vector is transformed. Finally these vectors are combined using the 4D cross product (which has 3 inputs), giving the final fundamental region (I'm not sure why the cross products are necessary, but I do know the main property of cross product is that the output is perpendicular to all of the inputs.). Note that some signs need wibbling, either that or permute the order of the inputs.
mat4 fr4(ivec3 pqr)
{
float r0 = radius(0, pqr);
float r1 = radius(1, pqr);
float r2 = radius(2, pqr);
float r3 = radius(3, pqr);
float c1 = r1 / r0;
float c2 = r2 / r1;
float c3 = r3 / r2;
float s1 = sqrt(1.0 - c1 * c1);
float s2 = sqrt(1.0 - c2 * c2);
float s3 = sqrt(1.0 - c3 * c3);
mat4 m1 = mat4(c1, s1, 0, 0, -s1, c1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1);
mat4 m2 = mat4(c2, 0, s2, 0, 0, 1, 0, 0, -s2, 0, c2, 0, 0, 0, 0, 1);
mat4 m3 = mat4(c3, 0, 0, s3, 0, 1, 0, 0, 0, 0, 1, 0, -s3, 0, 0, c3);
vec4 v0 = vec4(1, 0, 0, 0);
vec4 v1 = m1 * v0;
vec4 v2 = m1 * m2 * v0;
vec4 v3 = m1 * m2 * m3 * v0;
return (mat4
( normalize(cross4(v1, v2, v3))
, -normalize(cross4(v0, v2, v3))
, normalize(cross4(v0, v1, v3))
, -normalize(cross4(v0, v1, v2))
));
}
4D cross product can be implemented in terms of 3D determinants:
vec4 cross4(vec4 u, vec4 v, vec4 w)
{
mat3 m0 = mat3(u[1], u[2], u[3], v[1], v[2], v[3], w[1], w[2], w[3]);
mat3 m1 = mat3(u[0], u[2], u[3], v[0], v[2], v[3], w[0], w[2], w[3]);
mat3 m2 = mat3(u[0], u[1], u[3], v[0], v[1], v[3], w[0], w[1], w[3]);
mat3 m3 = mat3(u[0], u[1], u[2], v[0], v[1], v[2], w[0], w[1], w[2]);
return vec4(determinant(m0), -determinant(m1), determinant(m2), -determinant(m3));
}
The projection from 4D to 3D is stereographic, because signed distance fields are implicit we need both directions (one to go from input point in 3D to 4D, then after transformation / tesselation we need to go back to 3D to calculate distances):
vec4 unstereo(vec3 pos)
{
float r = length(pos);
return vec4(2.0 * pos, 1.0 - r * r) / (1.0 + r * r);
}
vec3 stereo(vec4 pos)
{
pos = normalize(pos);
return pos.xyz / (1 - pos.w);
}
float sdistance(vec4 a, vec4 b)
{
return distance(stereo(a), stereo(b));
}
The tiling is done iteratively (the limit of 600 is there because unbounded loops on a GPU are not really advisable, if this limit is set too low, e.g. 100, then visible artifacts occur - it's probably possible to prove an optimal bound somehow):
float poly4(vec4 r)
{
for (int i = 0, j = 0; i < 4 && j < 600; ++j)
{
if (dot(r, FR4[i]) < 0)
{
r = reflect(r, FR4[i]);
i = 0;
}
else
{
i++;
}
}
float de = 1.0 / 0.0;
// signed distance field stuff goes here
return de;
}
The user interface part of the code has some variables for selecting dimension and symmetry group (Schlaefli symbol). It also has 4 sliders for selecting the truncation amount in barycentric coordinates (which makes settings transfer in a meaningful way between polytopes), and 6 checkboxes for enabling different planes (there are 4 ways to choose the first axis and 3 left to choose from for the second axis, divided by 2 because the order doesn't matter).
vec4 s = inverse(transpose(FR4)) * vec4(BX, BY, BZ, BW);
The signed distance field stuff is quite straightforward in the end, though it took a lot of trial and error to get there. The first way I tried was to render tubes for line segments, by projecting the input point 'r' onto the planes of the fundamental region and doing an SDF circle:
float d = sdistance(s, r - FR4[0] * dot(r, FR4[0])) - thickness;
The planes are drawn by projecting 's' onto the cross product of the plane with 'r'. I don't know why this works:
vec4 p = normalize(cross4(FR4[0], FR4[1], r));
float d = sdistance(s, s - p * dot(s, p)) - 0.01;
Finally the DE() function for plugging into the DE-Raytracer.frag that comes with FragM has some animated rotation based on time, and the baseColor() function textures the solid with light and dark circles (actually slices of hyperspheres).
Full code download: reflex.frag.