Polychorons

This generates animated 3-dimensional Schlegel diagrams of 4-dimensional objects.

A live WebGL version of the code is at the top of the page.

The object is a tesseract, the 4-dimensional version of a cube. Its corners are pushed out onto a sphere in four dimensions, and then projected down into 3D from one point on that sphere, the same way some maps flatten the globe. Straight edges come out as curves, and everything is drawn with small cubes: an orange cube on each corner and a string of blue ones along each edge.

The tesseract

A polygon is a flat shape bounded by edges, a polyhedron is a solid bounded by flat faces, and a polychoron is the next one up: a 4-dimensional shape bounded by solid cells. The easiest one to think about is the tesseract, or hypercube, which is to a cube what a cube is to a square.

You can build it by counting. A square has its corners at (±1, ±1), four of them. A cube has its corners at (±1, ±1, ±1), eight of them. The tesseract has its corners at (±1, ±1, ±1, ±1), every combination of plus and minus one in four coordinates, which makes 16 vertices.

Two vertices are joined by an edge when they differ in exactly one coordinate. Each vertex has four neighbours, one for each coordinate you can flip, so there are 16 × 4 / 2 = 32 edges. (It also has 24 square faces and 8 cubic cells, but the program only draws the vertices and edges.)

The program makes the vertices by counting from 0 to 15 and using the four bits of the number as the four coordinates, with a 0 bit giving −1 and a 1 bit giving +1:

void defineSchlegelShape()
{
	bool done[256] = {false};

	//  Draw points

	for( int i=0; i<16; i++ )
	{
		int i0 = (i&1) >> 0;
		int i1 = (i&2) >> 1;
		int i2 = (i&4) >> 2;
		int i3 = (i&8) >> 3;

		points[ numPoints++ ] = vec4( 2.0*i0-1.0, 2.0*i1-1.0, 2.0*i2-1.0, 2.0*i3-1.0 );
	}
...

Two numbers that differ in one bit are two vertices that differ in one coordinate, so the test for an edge is that i^j is a single bit: 1, 2, 4 or 8. The program tries every pair of vertices, and a 256-entry done table makes sure each pair is only counted once, which gives the 32 edges:

...
			int k = i^j;

			if( !done[ (j<<4)+i ] )
			{
				if( k == 1 || k == 2 || k == 4 || k == 8 )
				{
					lines[ numLines++ ] = Pair<int>( i, j );
				}

				done[ (j<<4)+i ] = done[ (i<<4)+j ] = true;
			}
...

If you draw the eight vertices with w = −1 as a small cube and the eight with w = +1 as a big cube around it, you get the usual picture of a tesseract: a cube inside a cube, with eight more edges joining each inner corner to the matching outer one. The program's picture comes out in exactly this shape, as we'll see.

Setting up

The window and the OpenGL set-up come from a small Windows OpenGL helper library written by Mark Lord at Gremlin Interactive, in Winmain.cpp and qe_ogl.cpp. WinMain opens a 1600x1200 window called "OpenGL Test Program", asks for a double-buffered RGBA pixel format with 16 bits of colour and 16 bits of depth, and then loops forever, handling any waiting Windows messages and then drawing a frame:

	g_exitflag = FALSE;
	while (! g_exitflag)
	{
		ProcessMessages ();
		GameMain ();
	}

There's no frame limiter, so it draws as fast as the machine allows. Whenever the window changes size, GL_Initialise sets the viewport and a 60 degree perspective with the near plane at 1 and the far plane at 30000, and turns on depth testing. The clear colour is black, and back-face culling is commented out, so both sides of every face get drawn.

	// load a default projection matrix at 60 degree F.O.V.
	glMatrixMode (GL_PROJECTION);
	glLoadIdentity ();
	gluPerspective (60.0, (GLfloat)width/(GLfloat)height, 1.0, 30000.0);
	glMatrixMode (GL_MODELVIEW);
	glLoadIdentity ();
...
	// no blending just yet
	glDisable (GL_BLEND);
  glEnable (GL_DEPTH_TEST);

The far plane is a very long way out, which turns out to be useful: parts of the shape can fly a long way from the centre.

Each frame, GameMain clears the screen and places the camera. It moves everything 15 units away, tilts by r2 and turns by r1 (the two view angles the mouse changes, at 0.2 degrees per unit), and then moves the scene down by 2. So the camera circles around the point (0, 2, 0), which is the height the two shapes are drawn at.

void GameMain()
{
  Initialise();

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	glTranslatef(0.f, 0.f, -15.f);

	glRotatef(-r2*0.2, 1.0, 0.0, 0.0);
	glRotatef(-r1*0.2, 0.0, 1.0, 0.0);

	glTranslatef(0.f, -2.f, 0.f);

	GameDraw();
	SwapBuffers (gl_dc);
}

Initialise only does anything the first time it's called. It builds the tesseract's vertices and edges with defineSchlegelShape and makes the floor texture.

Turning in four dimensions

In 3D we're used to rotating about an axis. Turn a cube about the z axis and x and y change while z stays put. But really the rotation happens in the xy plane: two coordinates mix together and the rest stay still. In 3D there's exactly one axis left over for each plane, so the two ways of describing it are interchangeable.

In 4D they aren't. A turn in the xy plane leaves both z and w alone, so what stays still is a whole plane rather than an axis. And there are more planes to turn in, one for each pair of axes, which in four dimensions makes six: xy, xz, xw, yz, yw and zw.

A turn by an angle r in the plane of two axes is the ordinary 2D rotation, applied to just those two coordinates:

a1 = p1 cos r − p2 sin r
a2 = p1 sin r + p2 cos r

Here's the program's version, vec4::rotate, which takes the two axes as numbers from 0 to 3:

	vec4 rotate( int axis1, int axis2, double r )
	{
		vec4 a = *this;

		if( axis1 >= 0 && axis1 < 4 && axis2 >= 0 && axis2 < 4 )
		{
			a.p[ axis1 ] = p[ axis1 ] * sin(r) + p[ axis2 ] * cos(r);
			a.p[ axis2 ] = p[ axis1 ] * cos(r) - p[ axis2 ] * sin(r);
		}

		return a;
	}

Compare it with the textbook version and you'll see the two lines are the wrong way round: the first coordinate gets what should go into the second, and the second gets what should go into the first. So it's a proper turn by r followed by swapping the two coordinates. At r = 0 it doesn't leave the point alone, it swaps them: (1, 0) becomes (0, 1).

Swapping two coordinates is a mirror, a reflection in the line where they're equal. On its own a mirror would turn the shape inside out, like the difference between a left hand and a right hand. But rotate6 does six of them, one per plane, and the mirrors cancel out in pairs, so the result is still a proper 4D rotation. As r goes up the point still goes smoothly round its plane, just the other way.

	vec4 rotate6( double r[6] )
	{
		return
			rotate( 0, 1, r[0] ).
			rotate( 0, 2, r[1] ).
			rotate( 0, 3, r[2] ).
			rotate( 1, 2, r[3] ).
			rotate( 1, 3, r[4] ).
			rotate( 2, 3, r[5] );
	}

For a tesseract the swaps make no difference to what you see. A swap can be moved past the turns that come after it, which only changes which planes they turn in, so the whole thing works out as a shuffle of the coordinates followed by six ordinary turns. With all six angles at zero you can see the shuffle on its own: (x, y, z, w) becomes (w, z, y, x). The tesseract's vertices are every combination of ±1, so shuffling the coordinates gives the same tesseract back.

The six turns are done one after another, each on the result of the last, so the order matters. Six angles is the right number: 4D has six independent ways to turn, just as 3D has three.

The angles come from the clock. timeGetTime returns milliseconds, and each angle goes once round in its own time:

  double r[6] =
  {
	  2 * time*3.14159 /  7233.0,
	  2 * time*3.14159 /  9215.0,
	  2 * time*3.14159 / 15245.0,
	  2 * time*3.14159 /  6347.0,
	  2 * time*3.14159 / 13746.0,
	  2 * time*3.14159 /  7843.0
  };

The periods don't divide into each other, so the combined motion takes a very long time to repeat.

Three of the planes, xy, xz and yz, don't touch w, and the projection (coming up next) treats x, y and z all alike, so a turn in one of those on its own would just spin the picture like any 3D object. The other three mix w into x, y or z. They're the ones that do the 4D thing: vertices move between the inside of the picture and the outside, so the inner cube swells up through the outer one while the outer one shrinks inside it. With all six chained together and going at once, you see a mixture of both.

Schlegel diagrams

To see a 4D shape we have to get it down to 3D. We could simply drop the w coordinate, but then the inner and outer cubes would land on top of each other. A Schlegel diagram uses a perspective view instead, so that one cell of the shape is seen as the outside and the rest sit inside it. This program gets there in two steps: put every point on a sphere, then project from that sphere into 3D.

Every vertex of the tesseract is √(1 + 1 + 1 + 1) = 2 from the centre, so they all lie on a 4D sphere of radius 2. A sphere in 4D is a 3-dimensional surface, the 3-sphere, in the same way that the surface of a ball is 2-dimensional. vec4::normalise divides a point by its length, which puts it on the 3-sphere of radius 1:

	vec4 normalise()
	{
		double d = sqrt( p[0]*p[0] + p[1]*p[1] + p[2]*p[2] + p[3]*p[3] );

		if( d == 0 )
		{
			return vec4();
		}

		d = 1.0 / d;

		return vec4
		(
			p[0] * d,
			p[1] * d,
			p[2] * d,
			p[3] * d
		);
	}

Then comes the stereographic projection. It's easiest to see in 2D first. Take a circle, pick its top point as the pole, and draw a line from the pole through any other point on the circle. Where the line crosses the horizontal line through the centre is where that point lands.

Points below the centre land inside the circle and points above it land outside. A point on the equator stays where it is and the bottom point lands in the middle. As a point gets close to the pole its line gets flatter and it lands further and further away, and the pole itself would land at infinity.

Working out where it lands is one division. The line from the pole (0, 1) through the point (x, w) goes down 1 − w for every x across, so it reaches w = 0 at x / (1 − w). In 4D it's exactly the same with three coordinates instead of one. The pole is (0, 0, 0, 1) and the point lands at (x, y, z) / (1 − w):

	vec3 stereographicProjection()
	{
		double oo1mz = 1.0 / ( 1.0 - p[3] );

		return vec3
		(
			p[0] * oo1mz,
			p[1] * oo1mz,
			p[2] * oo1mz
		);
	}

(The name oo1mz is "one over one minus z", from the usual 3D version of the formula, where it's z that gets divided out.)

Let's try it on the tesseract. Vertex 15 is (1, 1, 1, 1), which has length 2, so normalising gives (½, ½, ½, ½). Projecting divides (½, ½, ½) by 1 − ½, which gives (1, 1, 1). Vertex 7 is (1, 1, 1, −1), which normalises to (½, ½, ½, −½) and projects to (½, ½, ½) / 1½ = (⅓, ⅓, ⅓).

So the eight vertices with w = +1 make a cube with its corners at ±1, the eight with w = −1 make a cube with its corners at ±⅓, and we get the cube inside a cube, one a third the size of the other. That's the still shape on the left of the demo.

Curved edges

An edge of the tesseract is a straight line in 4D, but the middle of it isn't on the sphere. The middle of the edge from vertex 7 to vertex 15 is (1, 1, 1, 0), which is only √3 from the centre. So drawSchlegalThickLine steps along the straight edge in hundredths, normalises each step onto the sphere, projects it, and draws a small cube there:

void drawSchlegalThickLine( vec4 a, vec4 b )
{
	for( int i=0; i<=100; i++ )
	{
		vec3 p = a.interpolate( b, i*0.01 ).normalise().stereographicProjection();
		double v = getStereographicVolume( p );
		double l = cubeRoot( v ) * 0.03;

		Cube( p, l, Colour( 0.3, 0.6, 0.9 ) );
	}
}

Normalising the points along a straight line between two points on a sphere puts them on an arc of a great circle, the sphere's version of a straight line. Stereographic projection has the handy property that it turns circles into circles, so every edge comes out as an arc of a circle in 3D. The exception is a circle that goes through the pole, which comes out as a straight line.

You can see both in the still shape. The edges of the two cubes bulge outwards: the middle of the edge from vertex 14, (−1, 1, 1, 1), to vertex 15 lands at about (0, 1.37, 1.37) rather than (0, 1, 1). But the eight edges along w, like the one from vertex 7 to vertex 15, lie on great circles through the pole, so they come out as straight spokes from the inner cube's corners to the outer cube's.

Now set it turning. As the shape rotates, its edges sweep all over the sphere, and every so often one passes close to the pole. There 1 − w is tiny, and the division sends the point a very long way out. A point on the sphere at a given w lands √((1 + w) / (1 − w)) from the centre, which is 1 for w = 0, about 14 for w = 0.99 and about 45 for w = 0.999.

With the camera only about 15 units from the centre, that part of the edge swings out past the camera and comes back in from the other side as the edge moves on past the pole. The 101 steps are evenly spaced along the edge in 4D, not in 3D, so a long way out they're far apart and the tube breaks up into a trail of separate cubes.

Sizing the cubes

Stereographic projection magnifies things more the further out they are. Something small at the centre of the picture comes out at half size, on the equator it's full size, and further out it keeps growing. At a point p the length scale is (1 + |p|²) / 2, so the volume scale is that cubed.

The program sizes each cube by how much the projection has magnified the space around it, so the cubes grow with the shape rather than staying a uniform string of beads. getStereographicVolume works out a volume factor and cubeRoot turns it back into a length:

double getStereographicVolume( vec3 a )
{
	double d = 1.0 + a.length2();
	return ( d*d ) / 8.0;
}

double cubeRoot( double a )
{
	return exp( log(a) / 3.0 );
}

The exact volume factor would be (1 + |p|²)³ / 8, whose cube root is exactly the length scale. The code squares instead of cubing, so its cubes grow as (1 + |p|²)^2/3 / 2 rather than (1 + |p|²) / 2. They still get bigger as they go out, just more gently than the space around them.

For the still shape, the inner vertex cubes are at |p|² = ⅓, which gives a half-width of 0.2 × (4/3)^2/3 / 2, about 0.12. The outer ones are at |p|² = 3, which gives 0.2 × 4^2/3 / 2, about 0.25. So the outer corners are about twice the size of the inner ones, where the exact factor would have made them three times.

Near the pole the sizes get big. At 14 units out a vertex cube is about 7 units across, and at 45 units out it's about 32, which with the camera 15 units away fills the screen.

Drawing with cubes

Everything in the picture is a cube. Cube takes two opposite corners and a colour and draws the six faces as triangle fans of four vertices each. There's no lighting. Instead the top and bottom get the full colour, the faces facing along x get 95% of it and the faces facing along z get 90%, which is just enough to show the shape of each cube.

void Cube(GLfloat x1, GLfloat y1, GLfloat z1, GLfloat x2, GLfloat y2, GLfloat z2, GLfloat r, GLfloat g, GLfloat b)
{
  glColor3f(r, g, b);

  glBegin(GL_TRIANGLE_FAN);
    glVertex3f(x1, y1, z1);
    glVertex3f(x1, y1, z2);
    glVertex3f(x2, y1, z2);
    glVertex3f(x2, y1, z1);
  glEnd();
...
  glColor3f(r*0.95, g*0.95, b*0.95);
...

A second version of Cube takes a centre, a half-width and a colour, which is what the drawing functions use. A vertex is an orange cube, Colour( 0.9, 0.6, 0.0 ), with a half-width of 0.2 times the size from the last section, and the edges are made of blue cubes, Colour( 0.3, 0.6, 0.9 ), at 0.03 times:

void drawSchlegalPoint( vec4 a )
{
	vec3 p = a.normalise().stereographicProjection();
	double v = getStereographicVolume( p );
	double l = cubeRoot( v ) * 0.2;

	Cube( p, l, Colour( 0.9, 0.6, 0.0 ) );
}

The cubes always line up with the x, y and z axes, whichever way the edge runs, so a tube is really a stack of overlapping boxes. Near the middle of the picture they're packed tightly enough to look like a solid tube.

That's a lot of cubes. Each shape is 16 + 32 × 101 = 3,248 cubes, so with two shapes a frame is 6,496 cubes, 38,976 faces and 155,904 calls to glVertex3f.

The turning shape is drawn by drawSchlegelShape. It rotates the 16 vertices with rotate6 and draws them, then draws each edge between the rotated vertices. Rotating just the ends is enough, because a rotation keeps straight lines straight: the line between the rotated ends is the rotated edge.

void drawSchlegelShape( double r[6] )
{
	vec4 *tpoints = new vec4[ numPoints ];

	for( int i=0; i<numPoints; i++ )
	{
		tpoints[i] = points[i].rotate6( r );

		drawSchlegalPoint( tpoints[i] );
	}

	for( int i=0; i<numLines; i++ )
	{
		drawSchlegalThickLine( tpoints[ lines[i].a() ], tpoints[ lines[i].b() ] );
	}

	delete[] tpoints;
}

The still shape is drawn by drawSchlegelShape2, which works out the vertices and edges from the bit patterns again every frame and draws them as they are, with no rotation.

The floor

Under the shapes there's a big chequered floor. Its texture is made once in Initialise: 256x256 grey pixels, where each pixel is its x and y XORed together.

  for( j=0; j<256; j++ ) for( i=0; i<256; i++ )
  {
    c = i^j;
    data[ (i+256*j)*3+0 ] = c;
    data[ (i+256*j)*3+1 ] = c;
    data[ (i+256*j)*3+2 ] = c;
  }

XORing the two coordinates makes squares within squares. The top bit of i^j flips every 128 pixels, the next every 64, and so on, so there's a checkerboard at every power of two at once, each one half the brightness of the one above it. It's filtered with GL_LINEAR both ways and has no mipmaps, so in the distance, where lots of texels land on one pixel, the floor shimmers as the view moves.

The floor is a single quad 1000 units across at y = 0, with texture coordinates running up to 200, so the texture repeats 200 times each way, once every 5 units. It's drawn in glColor4f( 0.3, 0.3, 0.3, 0.9 ) with blending on:

  glEnable(GL_BLEND);

	  glEnable(GL_TEXTURE_2D);							// Enable Texture Mapping ( NEW )

    // draw a textured floor
	  glBindTexture(GL_TEXTURE_2D, texture[0]);

    glColor4f( 0.3, 0.3, 0.3, 0.9 );

  //  time = 0;

	double distance_moved = 0;

    glBlendFunc(GL_SRC_ALPHA,GL_ONE_MINUS_SRC_ALPHA);

	  glBegin(GL_QUADS);
		  // Front Face
		  glTexCoord2f(   0.0f + distance_moved*(20.0/100.0),   0.0f ); glVertex3f( -500.0f, 0.0f,  500.0f );
		  glTexCoord2f( 200.0f + distance_moved*(20.0/100.0),   0.0f ); glVertex3f(  500.0f, 0.0f,  500.0f );
		  glTexCoord2f( 200.0f + distance_moved*(20.0/100.0), 200.0f ); glVertex3f(  500.0f, 0.0f, -500.0f );
		  glTexCoord2f(   0.0f + distance_moved*(20.0/100.0), 200.0f ); glVertex3f( -500.0f, 0.0f, -500.0f );
	  glEnd();

	  glDisable(GL_TEXTURE_2D);							// Enable Texture Mapping ( NEW )

  glDisable(GL_BLEND);

OpenGL's default texture environment is GL_MODULATE, which multiplies the texture's colour by the current colour. The texture is RGB, so its alpha counts as 1, and each floor pixel comes out as the texel times 0.3 with an alpha of 0.9. Blending with GL_SRC_ALPHA and GL_ONE_MINUS_SRC_ALPHA then mixes 90% of that with 10% of whatever is already there. (distance_moved would slide the texture along, but it's always 0.)

The comment in GameDraw says the floor is semi-transparent to simulate reflection: a copy of the shape drawn under the floor would show faintly through it. GameDraw has a place for that copy, but the call is commented out, so with only the black background underneath, the floor comes out as the texture at 27% brightness.

Each frame

Putting it together, GameDraw reads the clock, works out the six angles, draws the floor, and then draws the two shapes, each moved up 2 units: the turning one 5 units to the right and the still one 5 units to the left.

  //- Draw the thing above the floor

  glPushMatrix();
  glTranslatef( 5, 2, 0 );
  drawSchlegelShape( r );
  glPopMatrix();



  //- Draw a non animated thing

  glPushMatrix();
  glTranslatef( -5, 2, 0 );
  drawSchlegelShape2();
  glPopMatrix();

The floor is drawn first and writes to the depth buffer, so any part of a shape that swings below y = 0 is hidden by it. Then GameMain calls SwapBuffers to show the frame, and the loop goes round again.

Because the angles come from timeGetTime rather than a frame count, the shape turns at the same speed however fast the machine draws. (GameDraw also works out a moving point mx, my, mz from the clock, and there are functions for working out a camera position, getviewpoint and getscreenpoint, but none of them are used.)

Controls

Holding the left mouse button and dragging turns the view. WM_MOUSEMOVE reads the cursor's screen position, and while the button is down it adds the change since the last move to r1 and r2:

    case WM_MOUSEMOVE:
      POINT pt;
      GetCursorPos(&pt);

      if(lbutton) {
        r1 += (mx-pt.x);
        r2 += (my-pt.y);
      }

      mx = pt.x;
      my = pt.y;

      return 0;

GameMain turns these into degrees at 0.2 per pixel, so a drag of 450 pixels turns the view 90 degrees. There's no limit on the tilt, so you can look up from under the floor or go right over the top.

The WebGL version

The WebGL version at the top of the page is a port of Winmain.cpp, qe_ogl.cpp and test.cpp. The maths is the same, including the swap in rotate and the way the cubes are sized, and it's done in JavaScript numbers, which are doubles just like the original's. The canvas keeps the shape of the original window's drawing area, 1592x1173.

WebGL has no immediate mode, so instead of glVertex3f calls each cube becomes 36 vertices, two triangles for each fan, in a vertex buffer. The turning shape's buffer is rebuilt every frame with the same rotate, normalise and project steps, and the still shape is built once, as it never changes.

The fixed-function pipeline becomes two small shaders, one for the cubes and one for the floor. The cube shader just passes each vertex's colour through. GL_MODULATE is a single multiply in the floor's fragment shader:

    varying vec2 vTexCoord;
    uniform sampler2D uTexture;
    uniform vec4 uColour;
    void main() {
      //  A GL_RGB texture has an alpha of 1, so the alpha comes from glColor4f
      gl_FragColor = vec4(texture2D(uTexture, vTexCoord).rgb, 1.0) * uColour;
    }

A few things are different:

• Escape can't close a browser tab, so the port leaves it out.
• The original asks for 16-bit colour and a 16-bit depth buffer. The browser draws with 8 bits per colour channel, so the colours are smoother than they were on a 16-bit desktop.
• timeGetTime counts milliseconds since Windows started, so the original began wherever the clock happened to be. The port starts its clock at 1,520,000 milliseconds so it always opens the same way, and adding ?t= and another number to the demo's address starts it somewhere else.
• Drags are scaled to the original window's pixels, so a drag across the picture turns it as far as it did on Windows, and touch drags work too.