How do I compose a rotation matrix with human readable angles from scratch?
Solution 1:
So the question really is Understanding 4x4 homogenous transform matrices
well without the math behind the only thing that left is geometric representation/meaning which is far better for human abstraction/understanding.
-
So what the 4x4 matrix is?
It is representation of some Cartesian coordinate system and it is composed of:
-
3
basis vectors (one for each axis) red,green,blueSo if the red,green,blue vectors are perpendicular to each other then the coordinate system is orthogonal. If they are also unit vectors then it is orthonormal (like for example unit matrix).
-
origin point gray
-
projection and homogenous side (unmarked bottom rest of the matrix)
This part is there only for enabling rotation and translation at once, therefore point used must be homogenous that means in form
(x,y,z,w=1)
for points and(x,y,z,w=0)
for direction vectors. If it was just(x,y,z)
then the matrix would be3x3
and that is not enough for translation. I will not use any projections they are uneasy to explain geometrically.
This layout is from OpenGL notation there are also transposed representation out there (vectors are rows not columns)
now how to transform any point to/from this coordinate system:
g=M*l; l=Inverse(M)*g;
where:
-
M
is transform matrix -
l
isM
local coordinate system point (LCS) -
g
is global coordinate system point (GCS)
for the transposed version (DirectX) it is:
l=M*g; g=Inverse(M)*l;
That is because transposed orthogonal rotation matrix is also inverse of itself
- for more info see transform matrix anatomy and 3D graphic pipeline
-
-
how to visualize it
Yes you can draw the matrix numbers but they do not make sense at first look especially if the numbers are changing so draw the axises vectors as on image above. Where each axis is a line from
origin
toorigin + line_size*axis_vector
-
how to construct it
Just compute axis vectors and origin and put them inside matrix. To ensure orthogonality exploit cross product (but be careful with order of multiplicants to use the right direction) Here example of getting 3 basis vectors from direction
-
effects
- rotation is done by rotating the axises so you can compute each axis by parametric circle equation ...
- scaling is done by multiplying axises by scale factor
- skewing is just using non perpendicular axises
-
rotation
For most cases the incremental rotation is used. There are two types
-
local rotation
M'=M*rotation_matrix
it rotates around local coordinate axises like you will control plane or car or player ... Most engines/games do not use these and fake it with euler angles instead which is a cheap solution (have many quirks and problems) because most people who using OpenGL do not even know this is possible and rather stack list ofglRotate/glTranslate
calls... -
global rotation
M'=Inverse(Inverse(M)*rotation_matrix)
it rotates around global coordinate system axises.
where
rotation_matrix
is any standard rotation transform matrix.If you have different matrix layout (transposed) then the rotations local and global are computed the other way around ...
You can also compute your
rotation_matrix
from3
angles like:rotation_matrix=rotation_around_x(ax)*rotation_around_y(ay)*rotation_around_z(az);
see Wiki rotation matrices the 3D
Rx,Ry,Rz
fromBasic rotations
are what you need. As you can see they are just unit circle parametric equation really. The order of multiplication change how the angles converge to target position. This is called Euler angles and I do not use it (I integrate step changes instead which has no restrictions if done properly not to mention it is simpler).Anyway if you need you can convert transform matrix into euler angles relatively easily see:
- Is there a way to calculate 3D rotation on X and Y axis from a 4x4 matrix
-
-
glRotate
If you want
glRotate
which is rotation around arbitrary axis not by 3 angles then There is workaround:- create transform matrix
N
for that axis - then transform your matrix
M
to it - rotate
N
by angle - then transform
M
back fromN
to global coordinates
Or you can use Rodrigues_rotation_formula instead
To transform Matrix to/from Matrix in this case just transform axises as points and leave the origin as is but the origin of
N
must be (0,0,0)!!! or the vectors transformed must havew=0
instead. - create transform matrix
-
usage
Transformations are cumulative that means:
-
p'=M1*M2*M3*M4*p;
is the same asM=M1*M2*M3*M4; p'=M*p
So if you have many points to transform then you precompute all transformations to single matrix and use just it. Do not need to multiply points by all subsequent matrices. OK now the concept:
you should have
3
coordinate systems:- camera
C
- world (usually unit matrix)
- object
O
(each object have its own matrix)
so if you have cube with
8
vertexesp0,...,p7
then you have to perform transformation on each point from object local coordinates to camera local coordinates. Some gfx api do some of it so you apply only what you have to so you really need:p(i)'=inverse(C)*unit*M*p(i);
the transforms are cumulative and unit matrix does not change anything so:
Q=inverse(C)*M; p(i)'=Q*p(i);
so before drawing compute
Q
for drawed object then take each pointp(i)
of the object and compute the transformedp(i)'
and draw/use the transformed one ... Thep(i)'
is in local camera coordinate system (x,y of the screen) but there is no perspective there so before drawing you can also add any of the projection matrices and divide byz
cordinate at the end ... The projection is also cumulative so it can be also insideQ
-
[edit1] C++ example
//$$---- Form CPP ----
//---------------------------------------------------------------------------
// apart from math.h include you can ignore this machine generated VCL related code
#include <vcl.h>
#pragma hdrstop
#include "win_main.h"
#include <math.h>
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TMain *Main; // pointer to main window ...
//---------------------------------------------------------------------------
// Here is the important stuff some math first
//---------------------------------------------------------------------------
const double deg=M_PI/180.0;
double divide(double x,double y);
void matrix_mul (double *c,double *a,double *b); // c[16] = a[16] * b[16]
void matrix_mul_vector(double *c,double *a,double *b); // c[ 4] = a[16] * b[ 4]
void matrix_subdet (double *c,double *a); // c[16] = all subdets of a[16]
double matrix_subdet ( double *a,int r,int s);// = subdet(r,s) of a[16]
double matrix_det ( double *a); // = det of a[16]
double matrix_det ( double *a,double *b); // = det of a[16] and subdets b[16]
void matrix_inv (double *c,double *a); // c[16] = a[16] ^ -1
//---------------------------------------------------------------------------
double divide(double x,double y)
{
if (!y) return 0.0;
return x/y;
}
void matrix_mul (double *c,double *a,double *b)
{
double q[16];
q[ 0]=(a[ 0]*b[ 0])+(a[ 1]*b[ 4])+(a[ 2]*b[ 8])+(a[ 3]*b[12]);
q[ 1]=(a[ 0]*b[ 1])+(a[ 1]*b[ 5])+(a[ 2]*b[ 9])+(a[ 3]*b[13]);
q[ 2]=(a[ 0]*b[ 2])+(a[ 1]*b[ 6])+(a[ 2]*b[10])+(a[ 3]*b[14]);
q[ 3]=(a[ 0]*b[ 3])+(a[ 1]*b[ 7])+(a[ 2]*b[11])+(a[ 3]*b[15]);
q[ 4]=(a[ 4]*b[ 0])+(a[ 5]*b[ 4])+(a[ 6]*b[ 8])+(a[ 7]*b[12]);
q[ 5]=(a[ 4]*b[ 1])+(a[ 5]*b[ 5])+(a[ 6]*b[ 9])+(a[ 7]*b[13]);
q[ 6]=(a[ 4]*b[ 2])+(a[ 5]*b[ 6])+(a[ 6]*b[10])+(a[ 7]*b[14]);
q[ 7]=(a[ 4]*b[ 3])+(a[ 5]*b[ 7])+(a[ 6]*b[11])+(a[ 7]*b[15]);
q[ 8]=(a[ 8]*b[ 0])+(a[ 9]*b[ 4])+(a[10]*b[ 8])+(a[11]*b[12]);
q[ 9]=(a[ 8]*b[ 1])+(a[ 9]*b[ 5])+(a[10]*b[ 9])+(a[11]*b[13]);
q[10]=(a[ 8]*b[ 2])+(a[ 9]*b[ 6])+(a[10]*b[10])+(a[11]*b[14]);
q[11]=(a[ 8]*b[ 3])+(a[ 9]*b[ 7])+(a[10]*b[11])+(a[11]*b[15]);
q[12]=(a[12]*b[ 0])+(a[13]*b[ 4])+(a[14]*b[ 8])+(a[15]*b[12]);
q[13]=(a[12]*b[ 1])+(a[13]*b[ 5])+(a[14]*b[ 9])+(a[15]*b[13]);
q[14]=(a[12]*b[ 2])+(a[13]*b[ 6])+(a[14]*b[10])+(a[15]*b[14]);
q[15]=(a[12]*b[ 3])+(a[13]*b[ 7])+(a[14]*b[11])+(a[15]*b[15]);
for(int i=0;i<16;i++) c[i]=q[i];
}
void matrix_mul_vector(double *c,double *a,double *b)
{
double q[3];
q[0]=(a[ 0]*b[0])+(a[ 1]*b[1])+(a[ 2]*b[2])+(a[ 3]);
q[1]=(a[ 4]*b[0])+(a[ 5]*b[1])+(a[ 6]*b[2])+(a[ 7]);
q[2]=(a[ 8]*b[0])+(a[ 9]*b[1])+(a[10]*b[2])+(a[11]);
for(int i=0;i<3;i++) c[i]=q[i];
}
void matrix_subdet (double *c,double *a)
{
double q[16];
int i,j;
for (i=0;i<4;i++)
for (j=0;j<4;j++)
q[j+(i<<2)]=matrix_subdet(a,i,j);
for (i=0;i<16;i++) c[i]=q[i];
}
double matrix_subdet ( double *a,int r,int s)
{
double c,q[9];
int i,j,k;
k=0; // q = sub matrix
for (j=0;j<4;j++)
if (j!=s)
for (i=0;i<4;i++)
if (i!=r)
{
q[k]=a[i+(j<<2)];
k++;
}
c=0;
c+=q[0]*q[4]*q[8];
c+=q[1]*q[5]*q[6];
c+=q[2]*q[3]*q[7];
c-=q[0]*q[5]*q[7];
c-=q[1]*q[3]*q[8];
c-=q[2]*q[4]*q[6];
if (int((r+s)&1)) c=-c; // add signum
return c;
}
double matrix_det ( double *a)
{
double c=0;
c+=a[ 0]*matrix_subdet(a,0,0);
c+=a[ 4]*matrix_subdet(a,0,1);
c+=a[ 8]*matrix_subdet(a,0,2);
c+=a[12]*matrix_subdet(a,0,3);
return c;
}
double matrix_det ( double *a,double *b)
{
double c=0;
c+=a[ 0]*b[ 0];
c+=a[ 4]*b[ 1];
c+=a[ 8]*b[ 2];
c+=a[12]*b[ 3];
return c;
}
void matrix_inv (double *c,double *a)
{
double d[16],D;
matrix_subdet(d,a);
D=matrix_det(a,d);
if (D) D=1.0/D;
for (int i=0;i<16;i++) c[i]=d[i]*D;
}
//---------------------------------------------------------------------------
// now the object representation
//---------------------------------------------------------------------------
const int pnts=8;
double pnt[pnts*3]= // Vertexes for 100x100x100 cube centered at (0,0,0)
{
-100.0,-100.0,-100.0,
-100.0,+100.0,-100.0,
+100.0,+100.0,-100.0,
+100.0,-100.0,-100.0,
-100.0,-100.0,+100.0,
-100.0,+100.0,+100.0,
+100.0,+100.0,+100.0,
+100.0,-100.0,+100.0,
};
const int facs=6;
int fac[facs*4]= // faces (index of point used) no winding rule
{
0,1,2,3,
4,5,6,7,
0,1,5,4,
1,2,6,5,
2,3,7,6,
3,0,4,7,
};
double rep[16]= // 4x4 transform matrix of object (unit from start) at (0,0,+100)
{
1.0,0.0,0.0, 0.0,
0.0,1.0,0.0, 0.0,
0.0,0.0,1.0,100.0,
0.0,0.0,0.0,1.0,
};
double eye[16]= // 4x4 transform matrix of camera at (0,0,-150)
{
1.0,0.0,0.0, 0.0,
0.0,1.0,0.0, 0.0,
0.0,0.0,1.0,-150.0,
0.0,0.0,0.0,1.0,
};
//---------------------------------------------------------------------------
// this is how to draw it
//---------------------------------------------------------------------------
void obj(double *pnt,int pnts,int *fac,int facs,double *rep,double *ieye)
{
// variables for drawing
int i;
double p0[3],p1[3],p2[3],p3[3],m[16],d;
// gfx api variables (change to your stuff) Main is the main form of this application
TCanvas *scr=Main->bmp->Canvas;
double xs2=Main->ClientWidth/2,ys2=Main->ClientHeight/2;
double v=xs2*tan(30.0*deg); // 60 degree viewing angle perspective projection
matrix_mul(m,ieye,rep); // cumulate all needed transforms
for (i=0;i<facs*4;) // go through all faces
{
// convert all points of face
matrix_mul_vector(p0,m,&pnt[fac[i]*3]); i++;
matrix_mul_vector(p1,m,&pnt[fac[i]*3]); i++;
matrix_mul_vector(p2,m,&pnt[fac[i]*3]); i++;
matrix_mul_vector(p3,m,&pnt[fac[i]*3]); i++;
// here goes perspective divide by z coordinate if needed
d=divide(v,p0[2]); p0[0]*=d; p0[1]*=d;
d=divide(v,p1[2]); p1[0]*=d; p1[1]*=d;
d=divide(v,p2[2]); p2[0]*=d; p2[1]*=d;
d=divide(v,p3[2]); p3[0]*=d; p3[1]*=d;
// here is viewport transform (just translate (0,0) to middle of screen in this case
p0[0]+=xs2; p0[1]+=ys2;
p1[0]+=xs2; p1[1]+=ys2;
p2[0]+=xs2; p2[1]+=ys2;
p3[0]+=xs2; p3[1]+=ys2;
// draw quad
// I use VCL GDI TCanvas you use what you have ...
// and wireframe only to keep this simple (no Z buffer,winding culling,...)
scr->Pen->Color=clAqua; // perimeter wireframe
scr->MoveTo(p0[0],p0[1]);
scr->LineTo(p1[0],p1[1]);
scr->LineTo(p2[0],p2[1]);
scr->LineTo(p3[0],p3[1]);
scr->LineTo(p0[0],p0[1]);
// scr->Pen->Color=clBlue; // face cross to visualy check if I correctly generate the fac[]
// scr->MoveTo(p0[0],p0[1]);
// scr->LineTo(p2[0],p2[1]);
// scr->MoveTo(p1[0],p1[1]);
// scr->LineTo(p3[0],p3[1]);
}
}
//---------------------------------------------------------------------------
//---------------------------------------------------------------------------
void TMain::draw()
{
if (!_redraw) return;
bmp->Canvas->Brush->Color=clBlack;
bmp->Canvas->FillRect(TRect(0,0,xs,ys));
// compute inverse of camera need to compute just once for all objects
double ieye[16];
matrix_inv(ieye,eye);
// draw all objects
obj(pnt,pnts,fac,facs,rep,ieye);
Main->Canvas->Draw(0,0,bmp);
_redraw=false;
}
//---------------------------------------------------------------------------
__fastcall TMain::TMain(TComponent* Owner) : TForm(Owner)
{
// window constructor you can ignore this ... (just create a backbuffer bitmap here)
bmp=new Graphics::TBitmap;
bmp->HandleType=bmDIB;
bmp->PixelFormat=pf32bit;
pyx=NULL;
}
//---------------------------------------------------------------------------
void __fastcall TMain::FormDestroy(TObject *Sender)
{
// window destructor release memory ... also ignoe this
if (pyx) delete pyx;
delete bmp;
}
//---------------------------------------------------------------------------
void __fastcall TMain::FormResize(TObject *Sender)
{
// on resize event ... just resize/redraw backbuffer also can ignore this
xs=ClientWidth; xs2=xs>>1;
ys=ClientHeight; ys2=ys>>1;
bmp->Width=xs;
bmp->Height=ys;
if (pyx) delete pyx;
pyx=new int*[ys];
for (int y=0;y<ys;y++) pyx[y]=(int*) bmp->ScanLine[y];
_redraw=true;
}
//---------------------------------------------------------------------------
void __fastcall TMain::FormPaint(TObject *Sender)
{
// repaint event can ignore
_redraw=true;
}
//---------------------------------------------------------------------------
void __fastcall TMain::tim_redrawTimer(TObject *Sender)
{
// timer event to animate the cube ...
_redraw=true;
// rotate the object to see it in motion
double ang,c,s;
ang=5.0*deg; c=cos(ang); s=sin(ang); // rotate baround z by 5 degrees per timer step
double rz[16]= { c, s, 0, 0,
-s, c, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1 };
ang=1.0*deg; c=cos(ang); s=sin(ang); // rotate baround x by 1 degrees per timer step
double rx[16]= { 1, 0, 0, 0,
0, c, s, 0,
0,-s, c, 0,
0, 0, 0, 1 };
matrix_mul(rep,rep,rz);
matrix_mul(rep,rep,rx);
draw();
}
//---------------------------------------------------------------------------
here is how it looks like:
And GIF animation with back face culling:
[notes]
If you have more questions then comment me ...
[Edit2] basic 3D vector operations often needed
If you do not know how to compute vector operations like cross/dot products or absolute value see:
// cross product: W = U x V
W.x=(U.y*V.z)-(U.z*V.y)
W.y=(U.z*V.x)-(U.x*V.z)
W.z=(U.x*V.y)-(U.y*V.x)
// dot product: a = (U.V)
a=U.x*V.x+U.y*V.y+U.z*V.z
// abs of vector a = |U|
a=sqrt((U.x*U.x)+(U.y*U.y)+(U.z*U.z))
here my C++ vector math:
static double vector_tmp[3];
double divide(double x,double y) { if ((y>=-1e-30)&&(y<=+1e-30)) return 0.0; return x/y; }
double* vector_ld(double x,double y,double z) { double *p=vector_tmp; p[0]=x; p[1]=y; p[2]=z; return p;}
double* vector_ld(double *p,double x,double y,double z) { p[0]=x; p[1]=y; p[2]=z; return p;}
void vector_copy(double *c,double *a) { for(int i=0;i<3;i++) c[i]=a[i]; }
void vector_abs(double *c,double *a) { for(int i=0;i<3;i++) c[i]=fabs(a[i]); }
void vector_one(double *c,double *a)
{
double l=divide(1.0,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
c[0]=a[0]*l;
c[1]=a[1]*l;
c[2]=a[2]*l;
}
void vector_len(double *c,double *a,double l)
{
l=divide(l,sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])));
c[0]=a[0]*l;
c[1]=a[1]*l;
c[2]=a[2]*l;
}
void vector_neg(double *c,double *a) { for(int i=0;i<3;i++) c[i]=-a[i]; }
void vector_add(double *c,double *a,double *b) { for(int i=0;i<3;i++) c[i]=a[i]+b[i]; }
void vector_sub(double *c,double *a,double *b) { for(int i=0;i<3;i++) c[i]=a[i]-b[i]; }
void vector_mul(double *c,double *a,double *b) // cross
{
double q[3];
q[0]=(a[1]*b[2])-(a[2]*b[1]);
q[1]=(a[2]*b[0])-(a[0]*b[2]);
q[2]=(a[0]*b[1])-(a[1]*b[0]);
for(int i=0;i<3;i++) c[i]=q[i];
}
void vector_mul(double *c,double *a,double b) { for(int i=0;i<3;i++) c[i]=a[i]*b; }
void vector_mul(double *c,double a,double *b) { for(int i=0;i<3;i++) c[i]=a*b[i]; }
double vector_mul( double *a,double *b) { double c=0; for(int i=0;i<3;i++) c+=a[i]*b[i]; return c; } // dot
double vector_len(double *a) { return sqrt((a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2])); }
double vector_len2(double *a) { return (a[0]*a[0])+(a[1]*a[1])+(a[2]*a[2]); }
[Edit3] local rotations for camera and object control via keyboard
As this has been asked a lot lately here some example answers of mine with demos:
- stationary camera view control (partial pseudo inverse matrix)
- camera and player control (inverse matrix)
- How to preserve accuracy with cumulative transforms over time (full pseudo inverse matrix)
- rotundus style simple OpenGL/C++/VCL player control example