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.