This is what I come up with in C++:
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
#include "list.h"
//---------------------------------------------------------------------------
//--- ND Camera -------------------------------------------------------------
//---------------------------------------------------------------------------
const int ND = 4; // dimensions
double focal_length=2.5; // camera focal length (camera ais axis aligned)
double focus[ND]; // camera position (viewing in -Z direction)
double* perspective3D(double *_p) // camera perspective projection of ND point p to 3D <-1,+1> position on xy plane (Z=0)
{
int i,j;
double p[ND],w,dw;
static double q[3]={0.0,0.0,0.0};
// relative to camera
for (i=0;i<ND;i++) p[i]=_p[i]-focus[i];
// perspective projection
for (w=1.0,i=ND-1;i>=2;i--)
{
if (p[i]<=1e-10)
{ for (i=0;i<3;i++) q[i]=0; return q; } // behind camera
w*=focal_length/p[i];
}
for (j=0;j<2;j++) p[j]*=w;
// conversion to double[3]=(x,y,0.0)
for (i=0;(i<2)&&(i<ND);i++) q[i]=p[i]; q[2]=0.0;
return q;
}
//---------------------------------------------------------------------------
//--- ND Quad mesh ----------------------------------------------------------
//---------------------------------------------------------------------------
struct _quad // Quad face
{
double p[4][ND]; // points
_quad(){}
_quad(_quad& a) { *this=a; }
~_quad(){}
_quad* operator = (const _quad *a) { *this=*a; return this; }
//_quad* operator = (const _quad &a) { ...copy... return this; }
};
//---------------------------------------------------------------------------
void quad_hypercube(List<_quad> &quad,double a) // quad[] = hypercube (2a)^ND
{
if (ND<2) return;
int i0,i1,j,k;
double p[ND];
_quad q;
// set camera position and FOV
for (j=0;j<ND;j++) focus[j]=0.0;
for (j=2;j<ND;j++) focus[j]=-10.0*a/focal_length; // perspective projected axises should have offset
// clear faces
quad.num=0;
// perspective debug
double z,w;
#define add_quad(z,w) { q.p[0][0]=-a; q.p[0][1]=-a; q.p[0][2]=z; q.p[0][3]=w; q.p[1][0]=+a; q.p[1][1]=-a; q.p[1][2]=z; q.p[1][3]=w; q.p[2][0]=+a; q.p[2][1]=+a; q.p[2][2]=z; q.p[2][3]=w; q.p[3][0]=-a; q.p[3][1]=+a; q.p[3][2]=z; q.p[3][3]=w; quad.add(q); }
// iterate through all axis aligned i[0]i1 planes combinations
for (i0=0 ;i0<ND;i0++)
for (i1=i0+1;i1<ND;i1++)
{
// start offset
for (j=0;j<ND;j++) p[j]=-a; p[i0]=0; p[i1]=0;
// iterate all offset combinations
for (;;)
{
// add face
for (j=0;j<ND;j++) q.p[0][j]=p[j]; q.p[0][i0]-=a; q.p[0][i1]-=a;
for (j=0;j<ND;j++) q.p[1][j]=p[j]; q.p[1][i0]+=a; q.p[1][i1]-=a;
for (j=0;j<ND;j++) q.p[2][j]=p[j]; q.p[2][i0]+=a; q.p[2][i1]+=a;
for (j=0;j<ND;j++) q.p[3][j]=p[j]; q.p[3][i0]-=a; q.p[3][i1]+=a;
quad.add(q);
if (ND<=2) break;
// increment offset
for (k=0;;)
{
// first unused axis
while((k==i0)||(k==i1)) k++;
if (k>=ND) break;
p[k]=-p[k];
if (p[k]<0.0) k++;
else break;
}
if (k>=ND) break;
}
}
}
//---------------------------------------------------------------------------
void quad_draw(List<_quad> &quad)
{
for (int i=0;i<quad.num;i++)
{
glBegin(GL_QUADS);
for (int j=0;j<4;j++)
glVertex3dv(perspective3D(quad.dat[i].p[j]));
glEnd();
}
}
//---------------------------------------------------------------------------
to make this simple I used mine dynamic list template so:
List<double> xxx; is the same as double xxx[];
xxx.add(5); adds 5 to end of the list
xxx[7] access array element (safe)
xxx.dat[7] access array element (unsafe but fast direct access)
xxx.num is the actual used size of the array
xxx.reset() clears the array and set xxx.num=0
xxx.allocate(100) preallocate space for 100 items
usage:
// global storage
List<_quad> quad;
// mesh init
quad_hypercube(quad,0.5);
// render
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
// directional (normal shading)
float lightAmbient [4]={0.50,0.50,0.50,1.00}; // rozptylene nesmerove
float lightDiffuse [4]={0.50,0.50,0.50,1.00}; // smerove
float lightDirection[4]={0.00,0.00,+1.0,0.00}; // global smer svetla w=0
glLightfv(GL_LIGHT0,GL_AMBIENT ,lightAmbient );
glLightfv(GL_LIGHT0,GL_DIFFUSE ,lightDiffuse );
glLightfv(GL_LIGHT0,GL_POSITION,lightDirection);
glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_FALSE);
glEnable(GL_COLOR_MATERIAL);
glEnable(GL_LIGHT0);
glDisable(GL_LIGHTING);
glDepthFunc(GL_LEQUAL);
glDisable(GL_CULL_FACE);
glDisable(GL_TEXTURE_2D);
glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
glLineWidth(2.0);
glColor3f(0.5,0.5,0.5);
quad_draw(quad);
glLineWidth(1.0);
glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
glFlush();
SwapBuffers(hdc);
Here results:
As you can see for 4D+ there is some inconsistency but the generated coordinates looks OK so it is most likely some problem with my projection from ND to 2D somewhere (will investigate latter on) as the internal cubes (due to possibly wrong perspective) merges visually with the back face of the other cube. Most likely I would need to change the FOV of camera on per axis basis.
How it works
simply each face of the hypercube is axis aligned QUAD parallel to base plane. So the i0,i1 for loops generate all possible base planes in which I encode the 2D QUAD coordinates.
Then I generate all the combinations of +/-a for the other axises and add QUAD face for each combination. That is all.
To be more clear let assume 3D cube and plane face base plane xy. So i0=0; i1=1; and each face in cube parallel to xy will have the encoded 4 points x,y coordinates encoded by your gray code:
p0(-a,-a,?)
p1(-a,+a,?)
p2(+a,+a,?)
p3(+a,-a,?)
Now in 3D there is just one axis left (z) so for its each combination of -a and +a generate new face so:
p0(-a,-a,-a) // xy face 1
p1(-a,+a,-a)
p2(+a,+a,-a)
p3(+a,-a,-a)
p0(-a,-a,+a) // xy face 2
p1(-a,+a,+a)
p2(+a,+a,+a)
p3(+a,-a,+a)
When done move to next plane i0,i1 combination which is xz
p0(-a,?,-a)
p1(-a,?,+a)
p2(+a,?,+a)
p3(+a,?,-a)
and then generate the combinations again ...
p0(-a,-a,-a) // xz face 1
p1(-a,-a,+a)
p2(+a,-a,+a)
p3(+a,-a,-a)
p0(-a,+a,-a) // xz face 2
p1(-a,+a,+a)
p2(+a,+a,+a)
p3(+a,+a,-a)
and continue until no base plane combinations are left...
In 3D there are 3 base planes (xy,xz,yz) and 1 remainding axsis so 2^1 = 2 parallel faces per plane so 3*2 = 6 faces together.
In 4D you got 6 base planes (xy,xz,xw,yz,yw,zw) and 2 remainding axises give 4 combinations of +/-a meaning 2^2 = 4 parallel faces per base plane leading to 6*4 = 24 faces.
