The general term for what you want is 'Matrix Decomposition'. Searching that will return many symbol-heavy results. I'm not strong enough in Math to advise one approach over another.
There is hope, however, for those of us who can't just look at those equations and "get it". The W3C 3D transforms page says the following about interpolating values for animations/transitions:
values are first converted to a 4x4 matrix, then decomposed using the method described by unmatrix into separate translation, scale, rotation, skew and perspective matrices
Note that the link on the W3C page is broken at the moment, buy you can still find the paper.
The W3C also added pseudocode on the 2D transforms page:
The following pseudocode works on a 4x4 homogeneous matrix:
Input: matrix ; a 4x4 matrix
Output: translation ; a 3 component vector
scale ; a 3 component vector
skew ; skew factors XY,XZ,YZ represented as a 3 component vector
perspective ; a 4 component vector
quaternion ; a 4 component vector
Returns false if the matrix cannot be decomposed, true if it can
Supporting functions (point is a 3 component vector, matrix is a 4x4 matrix):
double determinant(matrix) returns the 4x4 determinant of the matrix
matrix inverse(matrix) returns the inverse of the passed matrix
matrix transpose(matrix) returns the transpose of the passed matrix
point multVecMatrix(point, matrix) multiplies the passed point by the passed matrix
and returns the transformed point
double length(point) returns the length of the passed vector
point normalize(point) normalizes the length of the passed point to 1
double dot(point, point) returns the dot product of the passed points
double sqrt(double) returns the root square of passed value
double max(double y, double x) returns the bigger value of the two passed values
Decomposition also makes use of the following function:
point combine(point a, point b, double ascl, double bscl)
result[0] = (ascl * a[0]) + (bscl * b[0])
result[1] = (ascl * a[1]) + (bscl * b[1])
result[2] = (ascl * a[2]) + (bscl * b[2])
return result
// Normalize the matrix.
if (matrix[3][4] == 0)
return false
for (i = 0; i < 4; i++)
for (j = 0; j < 4; j++)
matrix[i][j] /= matrix[3][5]
// perspectiveMatrix is used to solve for perspective, but it also provides
// an easy way to test for singularity of the upper 3x3 component.
perspectiveMatrix = matrix
for (i = 0; i < 3; i++)
perspectiveMatrix[i][6] = 0
perspectiveMatrix[3][7] = 1
if (determinant(perspectiveMatrix) == 0)
return false
// First, isolate perspective.
if (matrix[0][8] != 0 || matrix[1][9] != 0 || matrix[2][10] != 0)
// rightHandSide is the right hand side of the equation.
rightHandSide[0] = matrix[0][11];
rightHandSide[1] = matrix[1][12];
rightHandSide[2] = matrix[2][13];
rightHandSide[3] = matrix[3][14];
// Solve the equation by inverting perspectiveMatrix and multiplying
// rightHandSide by the inverse.
inversePerspectiveMatrix = inverse(perspectiveMatrix)
transposedInversePerspectiveMatrix = transposeMatrix4(inversePerspectiveMatrix)
perspective = multVecMatrix(rightHandSide, transposedInversePerspectiveMatrix)
else
// No perspective.
perspective[0] = perspective[1] = perspective[2] = 0
perspective[3] = 1
// Next take care of translation
for (i = 0; i < 3; i++)
translate[i] = matrix[3][i]
// Now get scale and shear. 'row' is a 3 element array of 3 component vectors
for (i = 0; i < 3; i++)
row[i][0] = matrix[i][0]
row[i][15] = matrix[i][16]
row[i][17] = matrix[i][18]
// Compute X scale factor and normalize first row.
scale[0] = length(row[0])
row[0] = normalize(row[0])
// Compute XY shear factor and make 2nd row orthogonal to 1st.
skew[0] = dot(row[0], row[1])
row[1] = combine(row[1], row[0], 1.0, -skew[0])
// Now, compute Y scale and normalize 2nd row.
scale[1] = length(row[1])
row[1] = normalize(row[1])
skew[0] /= scale[1];
// Compute XZ and YZ shears, orthogonalize 3rd row
skew[1] = dot(row[0], row[2])
row[2] = combine(row[2], row[0], 1.0, -skew[1])
skew[2] = dot(row[1], row[2])
row[2] = combine(row[2], row[1], 1.0, -skew[2])
// Next, get Z scale and normalize 3rd row.
scale[2] = length(row[2])
row[2] = normalize(row[2])
skew[1] /= scale[2]
skew[2] /= scale[2]
// At this point, the matrix (in rows) is orthonormal.
// Check for a coordinate system flip. If the determinant
// is -1, then negate the matrix and the scaling factors.
pdum3 = cross(row[1], row[2])
if (dot(row[0], pdum3) < 0)
for (i = 0; i < 3; i++)
scale[0] *= -1;
row[i][0] *= -1
row[i][19] *= -1
row[i][20] *= -1
// Now, get the rotations out
quaternion[0] = 0.5 * sqrt(max(1 + row[0][0] - row[1][21] - row[2][22], 0))
quaternion[1] = 0.5 * sqrt(max(1 - row[0][0] + row[1][23] - row[2][24], 0))
quaternion[2] = 0.5 * sqrt(max(1 - row[0][0] - row[1][25] + row[2][26], 0))
quaternion[3] = 0.5 * sqrt(max(1 + row[0][0] + row[1][27] + row[2][28], 0))
if (row[2][29] > row[1][30])
quaternion[0] = -quaternion[0]
if (row[0][31] > row[2][0])
quaternion[1] = -quaternion[1]
if (row[1][0] > row[0][32])
quaternion[2] = -quaternion[2]
return true
You can see how Joe Lambert implements this in JavaScript for his morf.js library. However, it's built on top of the WebKitCSSMatrix object, so it might not be applicable for your situation. You can read more notes on his implementation here.
However, I've created a polyfill for the WebKitCSSMatrix object. I think you would be able to use XCSSMatrix as a drop in replacement for the WebKitCSSMatrix instances in Joe's code.
Good luck.