Is there a way to calculate 3D rotation on X and Y axis from a 4x4 matrix

Question

First of all, I am not a mathematical expert at all. Please be tolerant to my mathematical mistakes and correct me where necessary, I'd love to learn.

I have a cube which is rotating using css animations with transform: matrix3d(4x4). I can also manually rotate the cube, converting user actions to the same matrix3d transformations.

What I want is a rotating cube with css when the user stops interacting that starts from where the user left it. This is something I am successfully doing by getting the cube's transform matrix3d value and using multiplication to set the css's keyframes dynamically.

However when the user starts interacting with the cube, the cube jumps to its last known manual rotation point and continues from there since I can't figure out how to get the rotation on X and Y axis from the 4x4 matrix.

I am currently using the following library, Rematrix, which helps me in the part to go from manual rotation to css rotation, as described above.

I've been looking into articles about Euler, and how to go from Euler to matrixes and visa versa but like I mentioned before, this is where my lack of mathematical knowledge is holding me back I think. I can't seem to figure it out.

As a reference, here are some of the articles I've read to try and solve my problem.

• https://medium.com/@behreajj/3d-rotations-in-processing-vectors-matrices-quaternions-10e2fed5f0a3
• http://www.gregslabaugh.net/publications/euler.pdf
• https://www.learnopencv.com/rotation-matrix-to-euler-angles/
• https://css-tricks.com/get-value-of-css-rotation-through-javascript/

The last source makes the most sense to me but, if I'm correct, is not useful in this case since it is about 2D transformations, and not 3D.

I get the current matrix3d the following way:

const style = getComputedStyle(this.element).transform
const matrix = Rematrix.parse(style)

For manual rotation I use matrix multiplication based on the user's mouse positions (positionY, positionX).

const r1 = Rematrix.rotateX(this.positionY)
const r2 = Rematrix.rotateY(this.positionX)

const transform = [r1, r2].reduce(Rematrix.multiply)

this.element.style[userPrefix.js + 'Transform'] = Rematrix.toString(transform)

Going from manual to css rotation I use the following function:

const setCssAnimationKeyframes = (lastTransform, animationData) => {
  const rotationIncrement = 90

  let matrixes = []

  for (let i = 0; i < 5; i++) {
    const rX = Rematrix.rotateX(rotationIncrement * i)
    const rY = Rematrix.rotateY(rotationIncrement * i)

    const matrix = [lastTransform, rX, rY].reduce(Rematrix.multiply);

    matrixes.push(matrix)
  }

  animationData.innerHTML = `
    @keyframes rotateCube {
      0% {
        transform: ${Rematrix.toString(matrixes[0])};
      }
      25% {
        transform: ${Rematrix.toString(matrixes[1])};
      }
      50% {
        transform: ${Rematrix.toString(matrixes[2])};
      }
      75% {
        transform: ${Rematrix.toString(matrixes[3])}};
      }
      100% {
        transform: ${Rematrix.toString(matrixes[4])};
      }
    }
  `;
}

Please provide answers or comment with any useful information. Although it would most welcome, I do not expect you to provide a fully working code example. Any useful information, in any form, is much appreciated.

Answer 1

First read:

• Understanding 4x4 homogenous transform matrices

as I use terminology from there.

Well I was too lazy to equate the whole stuff for my environment but based on this:

• Computing Euler angles from a rotation matrix

The resulting 3D rotation sub matrix of m for any rotation order will always have these therms:

(+/-)sin(a)
(+/-)sin(b)cos(a)
(+/-)cos(b)cos(a)
(+/-)sin(c)cos(a)
(+/-)cos(c)cos(a)

Only their sign and location will change with transform order and conventions. So to identify them do this:

1. let set some non trivial euler angles first

   their |sin| , |cos| values must be different so none of the 6 values will be the same otherwise this will not work !!!

   I chose these:

   ex = 10 [deg]
   ey = 20 [deg]
   ez = 30 [deg]
   

2. compute rotation matrix m

   so apply 3 euler rotations on unit matrix in their order. In mine setup the resulting matrix looks like this:

   double m[16] = 
    { 
     0.813797652721405, 0.543838143348694,-0.204874128103256, 0, // Xx,Xy,Xz,0.0
    -0.469846308231354, 0.823172926902771, 0.318795770406723, 0, // Yx,Yy,Yz,0.0
     0.342020153999329,-0.163175910711288, 0.925416529178619, 0, // Zx,Zy,Zz,0.0
     0                , 0                , 0                , 1  // Ox,Oy,Oz,1.0
    };
   

   note that I am using OpenGL conventions the basis vectors X,Y,Z and origin O are represented by the lines of matrix and the matrix is direct.

3. identify (+/-)sin(a) therm

   the a can be any of the euler angles so print sin of them all:

   sin(ex) = 0.17364817766693034885171662676931
   sin(ey) = 0.34202014332566873304409961468226
   sin(ez) = 0.5
   

   now see m[8] = sin(ey) so we found our therm... Now we know:

   ey = a = asin(m[8]);
   

4. identify (+/-)???(?)*cos(a) therms

   simply print cos(?)*cos(ey) for the unused angles yet. so if ey is the 20 deg I print 10 and 30 deg ...

   sin(10 deg)*cos(20 deg) = 0.16317591116653482557414168661534
   cos(10 deg)*cos(20 deg) = 0.92541657839832335306523309767123
   sin(30 deg)*cos(20 deg) = 0.46984631039295419202705463866237
   cos(30 deg)*cos(20 deg) = 0.81379768134937369284469321724839
   

   when we look at the m again we can cross match:

   sin(ex)*cos(ey) = 0.16317591116653482557414168661534 = -m[9]
   cos(ex)*cos(ey) = 0.92541657839832335306523309767123 = +m[10]
   sin(ez)*cos(ey) = 0.46984631039295419202705463866237 = -m[4]
   cos(ez)*cos(ey) = 0.81379768134937369284469321724839 = +m[0]
   

   from that we can compute the angles ...

   sin(ex)*cos(ey) = -m[ 9]
   cos(ex)*cos(ey) = +m[10]
   sin(ez)*cos(ey) = -m[ 4]
   cos(ez)*cos(ey) = +m[ 0]
   ------------------------
   sin(ex) = -m[ 9]/cos(ey)
   cos(ex) = +m[10]/cos(ey)
   sin(ez) = -m[ 4]/cos(ey)
   cos(ez) = +m[ 0]/cos(ey)
   

   so finally:

   ---------------------------------------------
   ey = asin(m[8]);
   ex = atan2( -m[ 9]/cos(ey) , +m[10]/cos(ey) )
   ez = atan2( -m[ 4]/cos(ey) , +m[ 0]/cos(ey) )
   ---------------------------------------------
   

And that is it. If you got different layout/conventions/transform order this approach still should work... Only the indexes and signs change. Here small C++/VCL OpenGL example I test this on (X,Y,Z order):

//---------------------------------------------------------------------------
#include <vcl.h>
#include <math.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
bool _redraw=true;                  // need repaint?

//---------------------------------------------------------------------------
double m[16]=               // uniform 4x4 matrix
    {
    1.0,0.0,0.0,0.0,        // Xx,Xy,Xz,0.0
    0.0,1.0,0.0,0.0,        // Yx,Yy,Yz,0.0
    0.0,0.0,1.0,0.0,        // Zx,Zy,Zz,0.0
    0.0,0.0,0.0,1.0         // Ox,Oy,Oz,1.0
    };
double e[3]={0.0,0.0,0.0};  // euler angles x,y,z order
//---------------------------------------------------------------------------
const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
void matrix2euler(double *e,double *m)
    {
    double c;
    e[1]=asin(+m[ 8]);
    c=cos(e[1]); if (fabs(c>1e-20)) c=1.0/c; else c=0.0;
    e[0]=atan2(-m[ 9]*c,m[10]*c);
    e[2]=atan2(-m[ 4]*c,m[ 0]*c);
    }
//---------------------------------------------------------------------------
void gl_draw()
    {
    _redraw=false;
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

    glMatrixMode(GL_PROJECTION);
//  glLoadIdentity();
    glMatrixMode(GL_TEXTURE);
    glLoadIdentity();
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslated(0.0,0.0,-10.0);    // some distance from camera ...

    glDisable(GL_DEPTH_TEST);
    glDisable(GL_TEXTURE_2D);

    int i;
    // draw source matrix:
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glTranslated(-1.0,0.0,0.0); // source matrix on the left
    glMultMatrixd(m);
    glBegin(GL_LINES);
    glColor3f(1.0,0.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3d(1.0,0.0,0.0);
    glColor3f(0.0,1.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3d(0.0,1.0,0.0);
    glColor3f(0.0,0.0,1.0); glVertex3d(0.0,0.0,0.0); glVertex3d(0.0,0.0,1.0);
    glEnd();
    glPopMatrix();

    // draw source matrix:
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glTranslated(m[12],m[13],m[14]);    // source matrix in the middle
    glBegin(GL_LINES);
    glColor3f(1.0,0.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3dv(m+0);
    glColor3f(0.0,1.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3dv(m+4);
    glColor3f(0.0,0.0,1.0); glVertex3d(0.0,0.0,0.0); glVertex3dv(m+8);
    glEnd();
    glPopMatrix();

    // draw euler angles
    matrix2euler(e,m);
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
    glTranslated(+1.0,0.0,0.0); // euler angles on the right
    glRotated(e[0]*rad,1.0,0.0,0.0);
    glRotated(e[1]*rad,0.0,1.0,0.0);
    glRotated(e[2]*rad,0.0,0.0,1.0);
    glBegin(GL_LINES);
    glColor3f(1.0,0.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3d(1.0,0.0,0.0);
    glColor3f(0.0,1.0,0.0); glVertex3d(0.0,0.0,0.0); glVertex3d(0.0,1.0,0.0);
    glColor3f(0.0,0.0,1.0); glVertex3d(0.0,0.0,0.0); glVertex3d(0.0,0.0,1.0);
    glEnd();
    glPopMatrix();

//  glFlush();
    glFinish();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    gl_init(Handle);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glRotated(10.0,1.0,0.0,0.0);
    glRotated(20.0,0.0,1.0,0.0);
    glRotated(30.0,0.0,0.0,1.0);
    glGetDoublev(GL_MODELVIEW_MATRIX,m);
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    gl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
    {
    if (_redraw) gl_draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    gl_resize(ClientWidth,ClientHeight);
    _redraw=true;
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormKeyDown(TObject *Sender, WORD &Key, TShiftState Shift)
    {
//  Caption=Key;
    const double da=5.0;
    if (Key==37){ _redraw=true; glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadMatrixd(m); glRotated(+da,0.0,1.0,0.0); glGetDoublev(GL_MODELVIEW_MATRIX,m); glPopMatrix(); }
    if (Key==39){ _redraw=true; glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadMatrixd(m); glRotated(-da,0.0,1.0,0.0); glGetDoublev(GL_MODELVIEW_MATRIX,m); glPopMatrix(); }
    if (Key==38){ _redraw=true; glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadMatrixd(m); glRotated(+da,1.0,0.0,0.0); glGetDoublev(GL_MODELVIEW_MATRIX,m); glPopMatrix(); }
    if (Key==40){ _redraw=true; glMatrixMode(GL_MODELVIEW); glPushMatrix(); glLoadMatrixd(m); glRotated(-da,1.0,0.0,0.0); glGetDoublev(GL_MODELVIEW_MATRIX,m); glPopMatrix(); }
    }
//---------------------------------------------------------------------------

The only important stuff from it is matrix2euler function converting matrix m to euler angles in x,y,z order. It renders 3 coordinate systems axises. On the left is m used as modelview matrix, in the middle are basis vectors of m using identity modelview and on the right is modelview constructed by the euler angles computed ...

All 3 should match. If the left and middle are not a match then you got different convention of matrix or layout.

Here preview for the (10,20,30) [deg] test case:

It matches even after many rotations (arrow keys)...

The gl_simple.h can be found here:

• complete GL+GLSL+VAO/VBO C++ example

PS. Depending on the platform/environment the computation might need some edge case handling like rounded magnitude for asin bigger than 1, division by zero etc. Also atan2 has its quirks ...

[Edit1] Here the ultimate C++ example which does all this automatically:

//---------------------------------------------------------------------------
enum _euler_cfg_enum
    {
    _euler_cfg_a=0,
    _euler_cfg_b,
    _euler_cfg_c,
    _euler_cfg__sina,
    _euler_cfg_ssina,
    _euler_cfg__sinb_cosa,
    _euler_cfg_ssinb_cosa,
    _euler_cfg__cosb_cosa,
    _euler_cfg_scosb_cosa,
    _euler_cfg__sinc_cosa,
    _euler_cfg_ssinc_cosa,
    _euler_cfg__cosc_cosa,
    _euler_cfg_scosc_cosa,
    _euler_cfgs
    };
//---------------------------------------------------------------------------
void matrix2euler_init(double *e,double *m,int *cfg)    // cross match euler angles e[3] and resulting m[16] transform matrix into cfg[_euler_cfgs]
    {
    int i,j;
    double a,tab[4];
    const double _zero=1e-6;
    for (i=0;i<_euler_cfgs;i++) cfg[i]=-1;      // clear cfg
    // find (+/-)sin(a)
    for (i=0;i<3;i++)                           // test all angles in e[]
        {
        a=sin(e[i]);
        for (j=0;j<16;j++)                      // test all elements in m[]
         if (fabs(fabs(a)-fabs(m[j]))<=_zero)   // find match in |m[j]| = |sin(e[i])|
            {                                   // store configuration
            cfg[_euler_cfg_a]=i;            
            cfg[_euler_cfg__sina]=j;
            cfg[_euler_cfg_ssina]=(a*m[j]<0.0);
            j=-1; break;
            }
        if (j<0){ i=-1; break; }                // stop on match found
        }
    if (i>=0){ cfg[0]=-1; return;   }           // no match !!!
    // find (+/-)???(?)*cos(a)
    a=cos(e[cfg[_euler_cfg_a]]);
    i=0; if (i==cfg[_euler_cfg_a]) i++; tab[0]=sin(e[i])*a; tab[1]=cos(e[i])*a; cfg[_euler_cfg_b]=i;
    i++; if (i==cfg[_euler_cfg_a]) i++; tab[2]=sin(e[i])*a; tab[3]=cos(e[i])*a; cfg[_euler_cfg_c]=i;

    for (i=0;i<4;i++)
        {
        a=tab[i];
        for (j=0;j<16;j++)                      // test all elements in m[]
         if (fabs(fabs(a)-fabs(m[j]))<=_zero)   // find match in |m[j]| = |tab[i]|
            {                                   // store configuration
            cfg[_euler_cfg__sinb_cosa+i+i]=j;
            cfg[_euler_cfg_ssinb_cosa+i+i]=(a*m[j]<0.0);
            j=-1; break;
            }
        if (j>=0){ cfg[0]=-1; return;   }       // no match !!!
        }
    }
//---------------------------------------------------------------------------
void matrix2euler(double *e,double *m,int *cfg) // compute euler angles e[3] from transform matrix m[16] using confing cfg[_euler_cfgs]
    {
    double c;
    //-----angle------         --------------sign--------------     ----------index----------
    e[cfg[_euler_cfg_a]]=asin ((cfg[_euler_cfg_ssina]?-1.0:+1.0) *m[cfg[_euler_cfg__sina     ]]);
    c=cos(e[cfg[_euler_cfg_a]]); if (fabs(c>1e-20)) c=1.0/c; else c=0.0;
    e[cfg[_euler_cfg_b]]=atan2((cfg[_euler_cfg_ssinb_cosa]?-c:+c)*m[cfg[_euler_cfg__sinb_cosa]],
                               (cfg[_euler_cfg_scosb_cosa]?-c:+c)*m[cfg[_euler_cfg__cosb_cosa]]);
    e[cfg[_euler_cfg_c]]=atan2((cfg[_euler_cfg_ssinc_cosa]?-c:+c)*m[cfg[_euler_cfg__sinc_cosa]],
                               (cfg[_euler_cfg_scosc_cosa]?-c:+c)*m[cfg[_euler_cfg__cosc_cosa]]);
    }
//---------------------------------------------------------------------------

Usage:

const double deg=M_PI/180.0;
const double rad=180.0/M_PI;
// variables
double e[3],m[16];
int euler_cfg[_euler_cfgs];
// init angles
e[0]=10.0*deg;
e[1]=20.0*deg;
e[2]=30.0*deg;
// compute coresponding rotation matrix
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
glRotated(e[0]*rad,1.0,0.0,0.0);
glRotated(e[1]*rad,0.0,1.0,0.0);
glRotated(e[2]*rad,0.0,0.0,1.0);
glGetDoublev(GL_MODELVIEW_MATRIX,m);
// cross match e,m -> euler_cfg
matrix2euler_init(e,m,euler_cfg);

// now we can use
matrix2euler(e,m,euler_cfg);

This works for any order of transformation and or convention/layout. the init is called just once and then you can use the conversion for any transform matrix... You can also write your own optimized version based on the euler_cfg results for your environment.