mathr / blog / #

The 43MB video linked above is pretty much generated by these 70 lines of GLSL fragment shader:

#extension GL_EXT_gpu_shader4 : enable

uniform vec4 u; // unit
uniform vec4 v; // unit
uniform vec4 w;
uniform float bigness;
uniform vec3 lightPos;

void main() {
  vec2 z = bigness * gl_TexCoord[0].xy;
  vec4 p = u * z.x + v * z.y + w;
  // nearest sphere
  vec4 o1 = round(p);
  vec4 o2 = o1 + vec4(0.5);
  if (p.x < o1.x) o2.x -= 1.0;
  if (p.y < o1.y) o2.y -= 1.0;
  if (p.z < o1.z) o2.z -= 1.0;
  if (p.w < o1.w) o2.w -= 1.0;
  vec4 d1 = o1 - p;
  vec4 d2 = o2 - p;
  float l1 = dot(d1,d1);
  float l2 = dot(d2,d2);
  vec4 o = o1;
  bool odd = false;
  if (l2 < l1) { o = o2; odd = true; }
  // intersection circle
  float ou = dot(u, o - w);
  float ov = dot(v, o - w);
  vec4 c = u * ou + v * ov + w;
  vec4 oc = o - c;
  float d = dot(oc, oc);
  float r2 = 0.25 - d;
  vec2 dz = z - vec2(ou, ov);
  float dr2 = dot(dz, dz);
  vec4 colour = vec4(0.0, 0.0, 0.0, 0.0);
  float alpha = 0.0;
  float brightness = 0.0;
  float specular = 0.0;
  if (dr2 < r2) {
    // inside circle
    float k = o.x + o.y + o.z + o.w;
    bool even = abs(2.0 * floor(k / 2.0) - k) < 0.001;
    if (odd) {
      if (even) {
        colour = vec4(1.0, 0.7, 0.2, 1.0);
      } else {
        colour = vec4(0.7, 1.0, 0.2, 1.0);
      }
    } else {
      if (even) {
        colour = vec4(1.0, 0.2, 0.7, 1.0);
      } else {
        colour = vec4(0.7, 0.2, 1.0, 1.0);
      }
    }
    // sphere shading
    const float shine = 16.0;
    float height = sqrt(r2 - dr2);
    vec3 pos = vec3(z, height);
    vec3 light = normalize(lightPos - pos);
    vec3 normal = normalize(vec3(dz, height));
    float dlm = max(dot(light, normal), 0.0);
    vec3 reflected = 2.0 * dlm * normal - light;
    specular = pow(max(reflected.z, 0.0), shine);
    brightness = 0.0625 + dlm;
  } else {
    discard;
  }
  gl_FragColor = vec4(mix(colour.rgb * brightness, vec3(1.0), specular), 1.0);
}

In 3D Euclidean (flat) space, you could imagine filling it with cubes of side length 1, much like graph paper in 2D. Then you could fill each cube with a sphere of radius 1/2. In 3D the gap between those spheres (at the corners of each cube) has diameter sqrt(3) - 1 (about 0.732). But, in 4D space a similar hypercube lattice each filled with spheres leaves a gap with diameter sqrt(4) - 1 (exactly 1). This means you can fit another radius 1/2 sphere at each corner gap.

The spheres with integer coordinates (the ones at the corners) can be grouped together, and the ones with integer+1/2 coordinates (the ones in the middle of the cubes) can form another group. Moreover, you can use a parity argument to further subgroup each half: determining whether the sum of the coordinates is odd or even gives four groups of spheres in total, arranged in a hard to imagine 4D chessboard-like pattern.

In 3D space you might imagine taking a 2D slice through the grid of spheres. Depending where and at what angle you slice, you might get a regular square grid of circles and gaps or a more unpredictable pattern of differently sized circles and gaps where you chop through different parts of the spheres. You can do the same in the 4D lattice, which can similarly give a regular square grid of circles or irregular-looking patterns.

So, enough imagining: how to actually visualize it. The first step is to define a 2D slice through 4D space:

P(u,v) = u*U + v*V + W

Here lowercase letters are scalar real numbers, and uppercase letters are vectors. For convenience, assume U and V have length 1 and point in different enough directions. W is the origin of the plane, and U and V can be thought of as the local axes within the 2D plane.

The next step given a point (u,v) in this plane (which will eventually be drawn as a pixel on the screen), is to find the center of the nearest sphere in the lattice. The nearest integer sphere can be found by rounding each coordinate of P(u,v), but the nearest half-integer sphere is a bit more tricky - the method in the shader source code above checks in each dimension if the rounded coordinate is bigger or smaller than unrounded coordinate, and picks the neighbouring half-integer coordinate so that the unrounded coordinate is between the two sphere coordinates. Then compare the actual distance between these sphere centers and P(u,v), and pick the closest.

Now we have the closest sphere centered at O, but we want to find the circle that is formed if our plane is slicing through that sphere. I thought the algebra would be hairier than it was, but luckily symmetry comes to the rescue. The center of the circle P(ou, ov) is the nearest point on the plane to the center of the sphere, which can be found by minimizing the distance function (knowing that the differential of a function is 0 at extrema). Leaving out the workings (and assuming U and V are unit length):

ou = dot(U, O - W) ; ov = dot(V, O - W)

But there might be no circle here at all, if (u,v) happens to fall in a gap. Luckily we know the sphere has radius 1/2, so we can check that ||P(ou, ov) - O|| < 1/2. The next step is to find the radius of the circle, which is quite straightforward by Pythagoras' Theorem - the line from the sphere center to the circle center is perpendicular to the plane of the circle, and the hypotenuse has length 1/2, and we already computed the length of the adjacent side.

Now (unless we bailed out with no circle) we have our starting point (u,v), the center of a circle that contains it (ou, ov), and the radius of the circle r. But we want to make it look like a sphere! So we need to find the height of the sphere's surface above the plane. Then for shiny lighting it's easy to calculate the surface normal (as a sphere is uniformly round, the normal points away from its center). With the surface position and normal, and a light moving around in 3D above our plane, standard Phong lighting can make it look like proper disco balls.