how should i handle (morphing) 4D objects in opengl?
For starters I would use 4D -> 3D projection instead of cut by hyperplane. The result is not the same but will get you closer to your goal (so you can latter upgrade this to cut). So similarly like in 3D -> 2D conversions used in graphics you got 2 choices one is using perspective projection and second is just ignoring the 4th dimension coordinate while rendering. I will use the latter as it is simpler.
structures
To make this as simple as I can I will use wire-frame instead of BR rendering. So you need to handle 4D mesh (wire-frame). I would use 2 tables:
double pnt[]; // 4D point list (x,y,z,u)
int lin[]; // lines point indexes (i0,i1)first one stores all the vertexes of your mesh and the second hold index pairs of points connected by lines in wire-frame representation.
transforms
If I would only ignore the 4th coordinate then we would not get the desired functionality. So to make the 4th dimension work we need to add 4D transform to orient our mesh in 4D before rendering. So use homogenous transform matrix and lets call ir
rep. In 4D it should be5x5orthonormal matrix with4x4rotation partrot.To make this even easier avoid smooth rotations for now (as in 4D that is not as easy) and compute random rotation
4x4matrix instead. So just set all the cells randomly<-1,+1>. Handle each row as basis vector. To make them orthonormal just make them unit and exploit cross product. For more info see:- How to generate a randomly oriented, high dimension circle in python?
render
simply convert point table by your transform matrix
(x',y',z',u',W) = rep * (x,y,z,u,1)then take the (x
,y,z`) and render ...
Here simple OpenGL/C++ example of 4D hyper cube:
//---------------------------------------------------------------------------
//--- Mesh 4D: ver 0.000 ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _mesh4D_h
#define _mesh4D_h
//---------------------------------------------------------------------------
#include <math.h>
#include "nd_math.h"
#include "list.h"
//---------------------------------------------------------------------------
const double pi = M_PI;
const double pi2 =2.0*M_PI;
const double pipol=0.5*M_PI;
const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
//---------------------------------------------------------------------------
class mesh4D
{
public:
matrix<5> rep; // 4D uniform 5x5 transform matrix
List<double> pnt; // 4D point list (x,y,z,u)
List<int> lin; // lines point indexes (i0,i1)
mesh4D() {}
mesh4D(mesh4D& a) { *this=a; }
~mesh4D() {}
mesh4D* operator = (const mesh4D *a) { *this=*a; return this; }
//mesh4D* operator = (const mesh4D &a) { ...copy... return this; }
void set_randomrep(); // random oriented uniform 4D transform matrix with origin (0,0,0,0)
void set_hypercube(double a);
void draw();
};
//---------------------------------------------------------------------------
void mesh4D::set_randomrep()
{
int i,j;
matrix<4> rot;
rep.unit();
rot.rnd();
rot.orthonormal();
for (i=0;i<4;i++)
for (j=0;j<4;j++)
rep[i][j]=rot[i][j];
}
void mesh4D::set_hypercube(double a)
{
rep.unit(); // reset orientation
pnt.num=0; // clear point list
lin.num=0; // clear line list
pnt.add(-a); pnt.add(-a); pnt.add(-a); pnt.add(-a);
pnt.add(+a); pnt.add(-a); pnt.add(-a); pnt.add(-a);
pnt.add(-a); pnt.add(+a); pnt.add(-a); pnt.add(-a);
pnt.add(+a); pnt.add(+a); pnt.add(-a); pnt.add(-a);
pnt.add(-a); pnt.add(-a); pnt.add(+a); pnt.add(-a);
pnt.add(+a); pnt.add(-a); pnt.add(+a); pnt.add(-a);
pnt.add(-a); pnt.add(+a); pnt.add(+a); pnt.add(-a);
pnt.add(+a); pnt.add(+a); pnt.add(+a); pnt.add(-a);
pnt.add(-a); pnt.add(-a); pnt.add(-a); pnt.add(+a);
pnt.add(+a); pnt.add(-a); pnt.add(-a); pnt.add(+a);
pnt.add(-a); pnt.add(+a); pnt.add(-a); pnt.add(+a);
pnt.add(+a); pnt.add(+a); pnt.add(-a); pnt.add(+a);
pnt.add(-a); pnt.add(-a); pnt.add(+a); pnt.add(+a);
pnt.add(+a); pnt.add(-a); pnt.add(+a); pnt.add(+a);
pnt.add(-a); pnt.add(+a); pnt.add(+a); pnt.add(+a);
pnt.add(+a); pnt.add(+a); pnt.add(+a); pnt.add(+a);
// A0
lin.add( 0+0); lin.add( 0+1);
lin.add( 0+1); lin.add( 0+3);
lin.add( 0+3); lin.add( 0+2);
lin.add( 0+2); lin.add( 0+0);
// A1
lin.add( 4+0); lin.add( 4+1);
lin.add( 4+1); lin.add( 4+3);
lin.add( 4+3); lin.add( 4+2);
lin.add( 4+2); lin.add( 4+0);
// A=A0+A1
lin.add( 0+0); lin.add( 4+0);
lin.add( 0+1); lin.add( 4+1);
lin.add( 0+2); lin.add( 4+2);
lin.add( 0+3); lin.add( 4+3);
// B0
lin.add( 8+0); lin.add( 8+1);
lin.add( 8+1); lin.add( 8+3);
lin.add( 8+3); lin.add( 8+2);
lin.add( 8+2); lin.add( 8+0);
// B1
lin.add(12+0); lin.add(12+1);
lin.add(12+1); lin.add(12+3);
lin.add(12+3); lin.add(12+2);
lin.add(12+2); lin.add(12+0);
// B=B0+B1
lin.add( 8+0); lin.add(12+0);
lin.add( 8+1); lin.add(12+1);
lin.add( 8+2); lin.add(12+2);
lin.add( 8+3); lin.add(12+3);
// hyper cube = A+B
lin.add( 0+0); lin.add( 8+0);
lin.add( 0+1); lin.add( 8+1);
lin.add( 0+2); lin.add( 8+2);
lin.add( 0+3); lin.add( 8+3);
lin.add( 0+4); lin.add( 8+4);
lin.add( 0+5); lin.add( 8+5);
lin.add( 0+6); lin.add( 8+6);
lin.add( 0+7); lin.add( 8+7);
}
//---------------------------------------------------------------------------
void mesh4D::draw()
{
int i,j;
double _zero=1e-3;
vector<5> a,b;
glBegin(GL_LINES);
for (i=0;i<lin.num;)
{
// extrac first point
j=lin[i]*4; i++;
a.a[0]=pnt[j]; j++;
a.a[1]=pnt[j]; j++;
a.a[2]=pnt[j]; j++;
a.a[3]=pnt[j]; j++;
a.a[4]=1.0; // W=1
// extrac second point
j=lin[i]*4; i++;
b.a[0]=pnt[j]; j++;
b.a[1]=pnt[j]; j++;
b.a[2]=pnt[j]; j++;
b.a[3]=pnt[j]; j++;
b.a[4]=1.0; // W=1
// transform
a=rep*a;
b=rep*b;
// render
glVertex3dv(a.a); // use just x,y,z
glVertex3dv(b.a); // use just x,y,z
}
glEnd();
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
I used mine dynamic list.h template so:
List<double> xxx; is the same as double xxx[];xxx.add(5); adds 5 to end of the listxxx[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 arrayxxx.reset() clears the array and set xxx.num=0xxx.allocate(100) preallocate space for 100 items
the nd_math.h is mine lib for N-Dimensional computations. What you need is just 4D,5D vector and 4x4, 5x5 matrix math from linear algebra.
Both libs are a bit big in size and also legal issues prevent me to share their code here.
The usage is simple:
// globals and init
mesh4D mesh
double animx=-50.0,danimx=0.0;
double animy= 0.0,danimy=2.0;
mesh.set_hypercube(0.5);
// render
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluOrtho2D( -2.0, 2.0, -2.0, 2.0 );
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glRotated(animx,1.0,0.0,0.0);
glRotated(animy,0.0,1.0,0.0);
mesh.draw();
glFlush();
SwapBuffers(hdc);
// some timer
animx+=danimx; if (animx>=360.0) animx-=360.0;
animy+=danimy; if (animy>=360.0) animy-=360.0;
call_render_here();
// on key press or mouse wheel or what ever
mesh.set_randomrep();
And here preview for some rep rotations ...
so this way you can render any wire-frame mesh (even BR rendering should works this way).
If you want to upgrade to the cut then you should take each Wire-frame line and compute its intersection with cutting hyperplane. If we chose hyperplane that goes through point
O(0,0,0,u_cut)
and has normal
N(0,0,0,1)
Then the task will simplify a lot. there are 3 options. Lets consider edge line with endpoints A,B:
no intersection
((A.u > u_cut)&&(B.u > u_cut)) || ((A.u < u_cut)&&(B.u < u_cut))just ignore such edge
1 intersection
((A.u >= u_cut)&&(B.u <= u_cut)) || ((A.u <= u_cut)&&(B.u >= u_cut))so compute the intersection via linear interpolation
x = A.x + (B.x-A.x)*(u_cut-A.u)/(B.u-A.u)
y = A.y + (B.y-A.y)*(u_cut-A.u)/(B.u-A.u)
z = A.z + (B.z-A.z)*(u_cut-A.u)/(B.u-A.u)and remember such point and also edge to which it belongs to.
fully inside
(A.u == u_cut)&&(B.u == u_cut)just remember both endpoints and also render this edge.
After all edges are processed this way then you need to analyze the remembered intersection points and create new edges from them based on connectivity info between edges. I did not do this yet so I can not help with this. I would try to connect remembered points sharing the same neighbor but not sure if that is enough in 4D.
For more info take a look at related QAs I found or answer:
- Visualising 4D objects in OpenGL
- 4D to 3D perspective projection
- 4D rotations
- ND vector and square matrix math template
- reper4D
- configurable ND rendering engine + hyperspiral
[Edit1] code with perspective
//---------------------------------------------------------------------------
//--- Mesh 4D: ver 0.001 ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _mesh4D_h
#define _mesh4D_h
//---------------------------------------------------------------------------
#include <math.h>
#include "nd_math.h"
#include "list.h"
//---------------------------------------------------------------------------
const double pi = M_PI;
const double pi2 =2.0*M_PI;
const double pipol=0.5*M_PI;
const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
//---------------------------------------------------------------------------
class mesh4D
{
public:
matrix<5> rep; // 4D uniform 5x5 transform matrix
List<double> pnt; // 4D point list (x,y,z,u)
List<int> lin; // lines point indexes (i0,i1)
mesh4D() {}
mesh4D(mesh4D& a) { *this=a; }
~mesh4D() {}
mesh4D* operator = (const mesh4D *a) { *this=*a; return this; }
//mesh4D* operator = (const mesh4D &a) { ...copy... return this; }
void set_randomrep(); // random oriented uniform 4D transform matrix with origin (0,0,0,0)
void set_hypercube(double a);
void draw();
};
//---------------------------------------------------------------------------
void mesh4D::set_randomrep()
{
int i,j;
matrix<4> rot;
rot.rnd();
rot.orthonormal();
for (i=0;i<4;i++)
for (j=0;j<4;j++)
rep[i][j]=rot[i][j];
}
//---------------------------------------------------------------------------
void mesh4D::set_hypercube(double a)
{
rep.unit(); // reset orientation
rep[0][4]=0.0; // set position
rep[1][4]=0.0;
rep[2][4]=0.0;
rep[3][4]=3.0*a;
pnt.num=0; // clear point list
lin.num=0; // clear line list
pnt.add(-a); pnt.add(-a); pnt.add(-a); pnt.add(-a);
pnt.add(+a); pnt.add(-a); pnt.add(-a); pnt.add(-a);
pnt.add(-a); pnt.add(+a); pnt.add(-a); pnt.add(-a);
pnt.add(+a); pnt.add(+a); pnt.add(-a); pnt.add(-a);
pnt.add(-a); pnt.add(-a); pnt.add(+a); pnt.add(-a);
pnt.add(+a); pnt.add(-a); pnt.add(+a); pnt.add(-a);
pnt.add(-a); pnt.add(+a); pnt.add(+a); pnt.add(-a);
pnt.add(+a); pnt.add(+a); pnt.add(+a); pnt.add(-a);
pnt.add(-a); pnt.add(-a); pnt.add(-a); pnt.add(+a);
pnt.add(+a); pnt.add(-a); pnt.add(-a); pnt.add(+a);
pnt.add(-a); pnt.add(+a); pnt.add(-a); pnt.add(+a);
pnt.add(+a); pnt.add(+a); pnt.add(-a); pnt.add(+a);
pnt.add(-a); pnt.add(-a); pnt.add(+a); pnt.add(+a);
pnt.add(+a); pnt.add(-a); pnt.add(+a); pnt.add(+a);
pnt.add(-a); pnt.add(+a); pnt.add(+a); pnt.add(+a);
pnt.add(+a); pnt.add(+a); pnt.add(+a); pnt.add(+a);
// A0
lin.add( 0+0); lin.add( 0+1);
lin.add( 0+1); lin.add( 0+3);
lin.add( 0+3); lin.add( 0+2);
lin.add( 0+2); lin.add( 0+0);
// A1
lin.add( 4+0); lin.add( 4+1);
lin.add( 4+1); lin.add( 4+3);
lin.add( 4+3); lin.add( 4+2);
lin.add( 4+2); lin.add( 4+0);
// A=A0+A1
lin.add( 0+0); lin.add( 4+0);
lin.add( 0+1); lin.add( 4+1);
lin.add( 0+2); lin.add( 4+2);
lin.add( 0+3); lin.add( 4+3);
// B0
lin.add( 8+0); lin.add( 8+1);
lin.add( 8+1); lin.add( 8+3);
lin.add( 8+3); lin.add( 8+2);
lin.add( 8+2); lin.add( 8+0);
// B1
lin.add(12+0); lin.add(12+1);
lin.add(12+1); lin.add(12+3);
lin.add(12+3); lin.add(12+2);
lin.add(12+2); lin.add(12+0);
// B=B0+B1
lin.add( 8+0); lin.add(12+0);
lin.add( 8+1); lin.add(12+1);
lin.add( 8+2); lin.add(12+2);
lin.add( 8+3); lin.add(12+3);
// hyper cube = A+B
lin.add( 0+0); lin.add( 8+0);
lin.add( 0+1); lin.add( 8+1);
lin.add( 0+2); lin.add( 8+2);
lin.add( 0+3); lin.add( 8+3);
lin.add( 0+4); lin.add( 8+4);
lin.add( 0+5); lin.add( 8+5);
lin.add( 0+6); lin.add( 8+6);
lin.add( 0+7); lin.add( 8+7);
}
//---------------------------------------------------------------------------
void mesh4D::draw()
{
int i,j;
const double _zero=1e-3;
double focal_length=1.0;
vector<5> a,b;
glBegin(GL_LINES);
for (i=0;i<lin.num;)
{
// extrac first point
j=lin[i]*4; i++;
a.a[0]=pnt[j]; j++;
a.a[1]=pnt[j]; j++;
a.a[2]=pnt[j]; j++;
a.a[3]=pnt[j]; j++;
a.a[4]=1.0; // W=1
// extrac second point
j=lin[i]*4; i++;
b.a[0]=pnt[j]; j++;
b.a[1]=pnt[j]; j++;
b.a[2]=pnt[j]; j++;
b.a[3]=pnt[j]; j++;
b.a[4]=1.0; // W=1
// transform
a=rep*a;
b=rep*b;
// perspective: camera projection plane u=0, focus at (0,0,0,-focal_length)
if (a[3]>=0.0) a*=divide(focal_length,a[3]+focal_length); else a.zero();
if (b[3]>=0.0) b*=divide(focal_length,b[3]+focal_length); else b.zero();
// render
glVertex3dv(a.a); // use just x,y,z
glVertex3dv(b.a); // use just x,y,z
}
glEnd();
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
And preview:
[Edit2] solid mesh and cross-section
so I changed the architecture quite a bit. I moved the 4D 5x5 homogenuous transform matrix (reper4D) to separate file and added colors and mesh definition by 4D simplexes (4 point 4 side tetrahedrons). The cut is simply computing the intersection (as described above) of simplex and cutting hyperplane resulting in either 3 points (triangle) , 4 points (tetrahedron) or 0 points. Which can be rendered easily (no need to analyze the connections between edges). For more info see this:
- Calculate 3D cross section of 4D shape (tesseract)
Btw. I think this is how Miegakure works. Here updated code:
//---------------------------------------------------------------------------
//--- Mesh 4D: ver 1.000 ----------------------------------------------------
//---------------------------------------------------------------------------
#ifndef _mesh4D_h
#define _mesh4D_h
//---------------------------------------------------------------------------
#include "list.h"
#include "reper4D.h"
//---------------------------------------------------------------------------
class mesh4D
{
public:
reper4D rep; // 4D uniform 5x5 transform matrix
List<double> pnt; // 4D point list (x,y,z,w)
List<int> lin; // 4D wireframe (i0,i1)
List<int> fac; // 4D simplexes (i0,i1,i2,i3)
List<DWORD> col; // simplex colors (RGB)
mesh4D() {}
mesh4D(mesh4D& a) { *this=a; }
~mesh4D() {}
mesh4D* operator = (const mesh4D *a) { *this=*a; return this; }
//mesh4D* operator = (const mesh4D &a) { ...copy... return this; }
void set_hypercube(double a);
void draw_cut(double w_cut); // render cross section by w=w_cut hyperplane
void draw (double focal_length=-1.0,double w_near=-1.0); // render mesh (focal_length<0) -> no perspective, else perspective view in W+ direction
void draw_wireframe(double focal_length=-1.0,double w_near=-1.0); // render wireframe (focal_length<0) -> no perspective, else perspective view in W+ direction
};
//---------------------------------------------------------------------------
void mesh4D::set_hypercube(double a)
{
const double tab_pnt[]=
{
-a, -a, -a, -a,
+a, -a, -a, -a,
-a, +a, -a, -a,
+a, +a, -a, -a,
-a, -a, +a, -a,
+a, -a, +a, -a,
-a, +a, +a, -a,
+a, +a, +a, -a,
-a, -a, -a, +a,
+a, -a, -a, +a,
-a, +a, -a, +a,
+a, +a, -a, +a,
-a, -a, +a, +a,
+a, -a, +a, +a,
-a, +a, +a, +a,
+a, +a, +a, +a,
};
const int tab_lin[]=
{
// A0
0+0, 0+1,
0+1, 0+3,
0+3, 0+2,
0+2, 0+0,
// A1
4+0, 4+1,
4+1, 4+3,
4+3, 4+2,
4+2, 4+0,
// A=A0+A1
0+0, 4+0,
0+1, 4+1,
0+2, 4+2,
0+3, 4+3,
// B0
8+0, 8+1,
8+1, 8+3,
8+3, 8+2,
8+2, 8+0,
// B1
12+0, 12+1,
12+1, 12+3,
12+3, 12+2,
12+2, 12+0,
// B=B0+B1
8+0, 12+0,
8+1, 12+1,
8+2, 12+2,
8+3, 12+3,
// hyper cube = A+B
0+0, 8+0,
0+1, 8+1,
0+2, 8+2,
0+3, 8+3,
0+4, 8+4,
0+5, 8+5,
0+6, 8+6,
0+7, 8+7,
};
// 5x simplex per cube
#define _cube(a0,a1,a2,a3,a4,a5,a6,a7) a1,a2,a4,a7, a0,a1,a2,a4, a2,a4,a6,a7, a1,a2,a3,a7, a1,a4,a5,a7
// 4D hypercube = 8 cubes
const int tab_fac[]=
{
_cube( 0, 1, 2, 3, 4, 5, 6, 7),
_cube( 0, 1, 2, 3, 8, 9,10,11),
_cube( 4, 5, 6, 7,12,13,14,15),
_cube( 8, 9,10,11,12,13,14,15),
_cube( 0, 1, 4, 5, 8, 9,12,13),
_cube( 0, 2, 4, 6, 8,10,12,14),
_cube( 1, 3, 5, 7, 9,11,13,15),
_cube( 2, 3, 6, 7,10,11,14,15),
};
#undef _cube
const DWORD tab_col[]=
{
// BBGGRR, BBGGRR, BBGGRR, BBGGRR, BBGGRR,
0x00FF0000,0x00FF0000,0x00FF0000,0x00FF0000,0x00FF0000,
0x0000FF00,0x0000FF00,0x0000FF00,0x0000FF00,0x0000FF00,
0x000000FF,0x000000FF,0x000000FF,0x000000FF,0x000000FF,
0x0000FFFF,0x0000FFFF,0x0000FFFF,0x0000FFFF,0x0000FFFF,
0x00FF00FF,0x00FF00FF,0x00FF00FF,0x00FF00FF,0x00FF00FF,
0x00FFFF00,0x00FFFF00,0x00FFFF00,0x00FFFF00,0x00FFFF00,
0x00FFFFFF,0x00FFFFFF,0x00FFFFFF,0x00FFFFFF,0x00FFFFFF,
0x004080FF,0x004080FF,0x004080FF,0x004080FF,0x004080FF,
};
int i,n;
vector<4> p;
rep.reset();
pnt.num=0; for (i=0,n=sizeof(tab_pnt)/sizeof(tab_pnt[0]);i<n;i++) pnt.add(tab_pnt[i]);
lin.num=0; for (i=0,n=sizeof(tab_lin)/sizeof(tab_lin[0]);i<n;i++) lin.add(tab_lin[i]);
fac.num=0; for (i=0,n=sizeof(tab_fac)/sizeof(tab_fac[0]);i<n;i++) fac.add(tab_fac[i]);
col.num=0; for (i=0,n=sizeof(tab_col)/sizeof(tab_col[0]);i<n;i++) col.add(tab_col[i]);
}
//---------------------------------------------------------------------------
void mesh4D::draw_cut(double w_cut)
{
const double _zero=1e-6;
const int edge2[]={0,1,0,2,0,3,1,2,2,3,3,1,-1}; // simplex wireframe i0,i1
const int edge3[]={0,1,2,3,0,1,3,1,2,3,2,0,-1}; // simplex triangles i0,i1,i2
int e,i,j,k,k0,k1,k2,inside[4];
DWORD rgb;
vector<4> p[4],q[4];
vector<3> xyz[4],nor,a,b;
for (i=0;i<fac.num;)
{
rgb=col[i>>2];
// extrac points (x,y,z,w)
for (k=0;k<4;k++)
{
j=fac[i]*4; i++;
p[k].a[0]=pnt[j]; j++;
p[k].a[1]=pnt[j]; j++;
p[k].a[2]=pnt[j]; j++;
p[k].a[3]=pnt[j]; j++;
// transform
rep.l2g(p[k],p[k]);
inside[k]=1;
}
// process edge2 and compute cross section cut intersection points
for (e=0,k=0;edge2[e]>=0;)
{
k0=edge2[e]; e++;
k1=edge2[e]; e++;
// fully inside
if (fabs(p[k0][3]-w_cut)+fabs(p[k1][3]-w_cut)<=_zero)
{
if ((k<4)&&(inside[k0])){ q[k]=p[k0]; k++; inside[k0]=0; }
if ((k<4)&&(inside[k1])){ q[k]=p[k1]; k++; inside[k1]=0; }
continue;
}
// no intersection
if (((p[k0][3]> w_cut)&&(p[k1][3]> w_cut))||((p[k0][3]< w_cut)&&(p[k1][3]< w_cut))) continue;
// 1 intersection
if (k<4)
{
q[k]=p[k1]-p[k0];
q[k]*=divide(w_cut-p[k0][3],p[k1][3]-p[k0][3]);
q[k]+=p[k0];
q[k][3]=w_cut;
k++;
continue;
}
}
// 4D -> 3D vector
for (k0=0;k0<k;k0++) for (k1=0;k1<3;k1++) xyz[k0][k1]=q[k0][k1];
// render triangle
if (k==3)
{
// normal
a=xyz[1]-xyz[0];
b=xyz[2]-xyz[1];
nor.cross(a,b);
nor.unit();
// render
glBegin(GL_TRIANGLES);
glNormal3dv(nor.a);
glColor4ubv((BYTE*)(&rgb));
glVertex3dv(xyz[0].a);
glVertex3dv(xyz[1].a);
glVertex3dv(xyz[2].a);
glEnd();
}
// render simplex
if (k==4)
for (e=0;edge3[e]>=0;)
{
k0=edge3[e]; e++;
k1=edge3[e]; e++;
k2=edge3[e]; e++;
// normal
a=xyz[k1]-xyz[k0];
b=xyz[k2]-xyz[k1];
nor.cross(a,b);
nor.unit();
// render
glBegin(GL_TRIANGLES);
glNormal3dv(nor.a);
glColor4ubv((BYTE*)(&rgb));
glVertex3dv(xyz[k0].a);
glVertex3dv(xyz[k1].a);
glVertex3dv(xyz[k2].a);
glEnd();
}
}
}
//---------------------------------------------------------------------------
void mesh4D::draw(double focal_length,double w_near)
{
const int edge3[]={0,1,2,3,0,1,3,1,2,3,2,0,-1}; // simplex triangles i0,i1,i2
int i,j,k,k0,k1,k2;
DWORD rgb;
vector<4> p;
vector<3> xyz[4],nor,a,b;
// 4D simplexes
glColor3f(0.3,0.3,0.3);
for (i=0;i<fac.num;)
{
rgb=col[i>>2];
// extrac points (x,y,z,w)
for (k=0;k<4;k++)
{
j=fac[i]*4; i++;
p[0]=pnt[j]; j++;
p[1]=pnt[j]; j++;
p[2]=pnt[j]; j++;
p[3]=pnt[j]; j++;
// transform
rep.l2g(p,p);
// perspective projection
if (focal_length>0.0)
{
p[3]-=w_near;
if (p[3]>=0.0) p*=divide(focal_length,p[3]+focal_length); else p.zero();
}
// 4D -> 3D vector
xyz[k].ld(p[0],p[1],p[2]);
}
// render simplex
for (k=0;edge3[k]>=0;)
{
k0=edge3[k]; k++;
k1=edge3[k]; k++;
k2=edge3[k]; k++;
// normal
a=xyz[k1]-xyz[k0];
b=xyz[k2]-xyz[k1];
nor.cross(a,b);
nor.unit();
// render
// glBegin(GL_LINE_LOOP);
glBegin(GL_TRIANGLES);
glNormal3dv(nor.a);
glColor4ubv((BYTE*)(&rgb));
glVertex3dv(xyz[k0].a);
glVertex3dv(xyz[k1].a);
glVertex3dv(xyz[k2].a);
glEnd();
}
}
}
//---------------------------------------------------------------------------
void mesh4D::draw_wireframe(double focal_length,double w_near)
{
int i,j,k;
vector<4> p[4];
// 4D wireframe
glColor3f(1.0,1.0,1.0);
glBegin(GL_LINES);
for (i=0;i<lin.num;)
{
// extrac points (x,y,z,w)
for (k=0;k<2;k++)
{
j=lin[i]*4; i++;
p[k].a[0]=pnt[j]; j++;
p[k].a[1]=pnt[j]; j++;
p[k].a[2]=pnt[j]; j++;
p[k].a[3]=pnt[j]; j++;
// transform
rep.l2g(p[k],p[k]);
// perspective projection
if (focal_length>0.0)
{
p[k][3]-=w_near;
if (p[k][3]>=0.0) p[k]*=divide(focal_length,p[k][3]+focal_length); else p[k].zero();
}
// render
glVertex3dv(p[k].a); // use just x,y,z
}
}
glEnd();
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------
And the preview for cross section render:
The worst part was to define the hypercube as set of simplexes ...
Using OpenGL shaders to plot 4D points in 2D in C++
Since shaders can work with 4 dimensional vectors and 4x4 matrix, we can be smart and use only one vertex shader to do the trick. This shader will take three inputs:
- the point data (a 4 floating vector),
- a selection 4x4 matrix,
- a projection 4x4 matrix.
All the magic is done by the selection matrix, which will map you 4 coordinates vector to planar coordinates. Let's call the point data v, the selection matrix S and P the projection matrix. The output of the vertex shader will be:
P * S * v
By setting correctly the coefficients of S, you will be able to select which component in v should appear in the result. For example, if you want to display YxB, then S will be:
So we have:
P is a standard projection matrix, which will help you to place your graph correctly (offsets and scale).
Following is an example of vertex shader implementation (not tested, may contain mistakes) for OpenGl3:
#version 130
in vec4 vin;
uniform mat4 s;
uniform mat4 p;
out vec4 vout;
void main()
{
vec4 tmp = s * vin;
// Depending on your projection matrix,
// you may force the 4th component to 1.
tmp[3] = 1;
vout = p * tmp;
}
Notes:glUniformMatrix4fv is the function to use to define S and P.
Navigating though the surface of a hypersphere in OpenGL
I did ask the same question on math.stackexchange.com, because the question is too math related and I could possibly get a better answer there. Look here if you are interested in the original answer I've actually used.
Being in 4D, there is not a single vector perpendicular to any rotation, but actually a whole plane of perpendicular vectors, and that was confusing me.
Also, in OpenGL, I soon figured out that I could simplify the computation of the given solution, because I could always assume I am at the 4D projection origin (0,0,0,-1), always facing (0,0,-1,0), then translate a given angle along z and w, and multiply it to the already accumulated GL_MODELVIEW_MATRIX from the previous frame. So I did MODELVIEW <= M x MODELVIEW instead of MODELVIEW <= MODELVIEW x M I would get if I simply called glMultMatrix(M).
How to use 4d rotors
I was able to use Rotors after learning more on the subject thanks to this youtube series about Geometric Algebra
It's really well explained and I recommend it to whoever wants use geometric algebra.
Related Topics
C++ Remove_If on a Vector of Objects
Visual Studio Project Out of Date
Qthread Emits Finished() Signal But Isrunning() Returns True and Isfinished() Returns False
Best Cross-Platform Method to Get Aligned Memory
Is Rebasing Dlls (Or Providing an Appropriate Default Load Address) Worth the Trouble
Multiplying a String by an Int in C++
Using Opengl Glutdisplayfunc Within Class
Window Border Width and Height in Win32 - How to Get It
Linking Boost Library with Boost_Use_Static_Lib Off on Windows
Using Pair as Key in a Map (C++/Stl)
What Is the Use of Volatile Keyword
Copy Constructor in C++ Is Called When Object Is Returned from a Function
Is There a Null Std::Ostream Implementation in C++ or Libraries
Understand Op Registration and Kernel Linking in Tensorflow
How to Install Llvm Toolchain for Eclipse Cdt