4X4 Matrix Decomposition
Tue, Apr 14, 2020
I was in need to decompose transformations given as 4X4 matrices into simpler ones. What I found on the web wasn’t suiting my need for understanding what I’m doing. So I had to figure it out myself.
The Math of Matrix Decomposition
The lots of math involved and my self-set goal to keep this site as free from Javascript as possible do not allow to use common Javascript libraries to typeset math formulars here. Therefore I exhumed my LaTeX knowledge and wrote a paper what describes the math.
The paper describes how to decompose 4X4 matrices, and on the way the handled matrices are becoming smaller, so when you are interested in 4X3 or 3X3 matrix decomposition you’ll find it in the paper, too. Just skip the first step(s).
A slightly abbreviated section for handling 2D cases like 3X2 and, on the way, 2X2 matrix decomposition is also included.
Source Code
Source code on the other hand can be shown here perfectly well.
Obviously the sources depend a lot on the internal representation of matrices, and there are many of them. As I’m lazy I just copied over code from my library. Although that is not open and you cannot use it directly if you are not a customer the naming and comments should definitely help to understand what is necessary to port it to your needs.
Decomposition of a 4X4 Matrix
This is a default method in my Transformation3D interface.
Vector3D is a 3D column vector. OrderedPair is a typed pair of values.
/**
* Decompose this transformation into simpler ones.
*
* Be prepared that there is not only one way to decompose a matrix.
*
* The algorithm will produced the following order (where some or all
* of the simple transformations might be missing if they are just
* the identity):
* <ol>
* <li>{@link TransformFlags#Projection}</li>
* <li>{@link TransformFlags#Translation}</li>
* <li>{@link TransformFlags#Rotation}</li>
* <li>{@link TransformFlags#Shearing}</li>
* <li>{@link TransformFlags#Scaling} (including possible mirroring)</li>
* </ol>
*
* @return a chain of tagged simple transformations like projection, translation, rotation, shearing
* and scaling that when multiplied in the given order recreate this
* transformation (within floating point accuracy).
* The order to apply this transformations to one object step by step is inverse!
* If this represents identity the
* returned list will be empty, indicating that no transformation needs to be applied.
* If this transformation has an 0 determinant it is not
* considered decomposable and {@code null} is returned.
*/
@Nullable
default List<OrderedPair<TransformFlags, Transformation3D>> decompose()
{
// Column vectors of intermediate 4X3 matrix
final Vector3D columnX = new Vector3D(getXX(), getYX(), getZX());
final Vector3D columnY = new Vector3D(getXY(), getYY(), getZY());
final Vector3D columnZ = new Vector3D(getXZ(), getYZ(), getZZ());
final Vector3D columnW = new Vector3D(getXW(), getYW(), getZW());
// 3X3 matrix used for various purposes
final Matrix3X3 m = Matrix3X3.fromColumns(columnX, columnY, columnZ);
final List<OrderedPair<TransformFlags, Transformation3D>> result = new LinkedList<>();
final double wx = getWX();
final double wy = getWY();
final double wz = getWZ();
final double ww = getWW();
if (wx != 0 || wy != 0 || wz != 0 || ww != 1) {
// The trick here is to assume that the transformation is a product
// of a projective matrix where only the last row differs from an identity
// matrix and a 4X3 matrix (in that order).
// Then the upper 3 rows come out unharmed, and the last row can be considered
// the transformation of a vector w. We only use a 3X3 sub matrix here for faster
// inversion, as the ww coordinate is easily calculated directly when we
// know w.
final Matrix3X3 subTransposedInverted = m.getTransposed().getInverse();
if (subTransposedInverted == null) {
Debug.warn("Decomposition of transformation with projective parts and uninvertable sub matrix: %0",
this);
return null;
}
final Vector3D w = subTransposedInverted.multVector(wx, wy, wz);
result.add(OrderedPair.create(TransformFlags.Projection,
new Projection3D(w, ww - w.mult(getXW(), getYW(), getZW()))));
}
// Translation is easy now when it is considered to be applied last (before a possible projection), which we do here.
if (!TrafoUtil.equals(columnW, IVector3D.ZERO)) {
result.add(OrderedPair.create(TransformFlags.Translation,
createTranslation(columnW)));
}
// Some early checks
if (m.isIdentity()) {
return result;
}
final double det = m.getDeterminant();
if (det == 0) {
Debug.warn("Decomposition of uninvertable subtransformation: %0", m);
return null;
}
if (det == 1 && isIdentity()) {
return result;
}
// Now what is left is a 3X3 matrix which is considered to be a combination
// of rotation, shearing and scaling.
// We make use of the fact that for a rotation R the following is always valid:
// R^T = R^{-1}
// We now consider the remaining 3X3 matrix M to be the combination of
// a rotation R and a shear/scale matrix S:
// M = R * S
// If we calculate
// M^T * M = (R * S)^T * (R * S)
// = (S^T * R^T) * R * S
// = S^T * R^{-1} * R * S // because R^T = R^{-1}
// it follows that
// M^T * M = S^T * S
// This means that we can derive S from the last equation. For this we
// assume that S is a triangle matrix in the form
// / a d e \
// S = | 0 b f |
// \ 0 0 c /
// Multiplying S^T * S results in
// / a² ad ae \
// S^T * S = | ad b²+d² bf+ed | = M^T * M
// \ ae bf+ed c²+e²+f² /
// Resolving this leads to the formulas for a,b,c,... below.
// The squares allow to take plus or minus values for the
// scaling coefficients on the diagonal. For symmetry we
// use either always + or always -, although switching just one
// of them would suffice. This selects 1 of the 8 possible solutions
// and takes care that any mirroring operation is put into the
// scaling, so all other transformations in the chain have det=1.
final Matrix3X3 mTm = m.getTransposed().mult(m);
final boolean minus = det < 0;
final double a = minus // a cannot be 0 because of det != 0 (and m^T already was invertable when we had projection)
? -Math.sqrt(mTm.getXX())
: Math.sqrt(mTm.getXX());
final double d = mTm.getXY() / a;
final double e = mTm.getXZ() / a;
final double b = minus // cannot be 0, see a
? -Math.sqrt(mTm.getYY() - d*d)
: Math.sqrt(mTm.getYY() - d*d);
final double f = (mTm.getYZ() - d * e) / b;
final double c = minus // cannot be 0, see a
? -Math.sqrt(mTm.getZZ() - e*e - f*f)
: Math.sqrt(mTm.getZZ() - e*e - f*f);
final Matrix3X3 shearAndScale = new Matrix3X3(a, d, e,
0, b, f,
0, 0, c);
// Having shear-and-scale we can calculate the rotation from
// R * S = M <=> R = M * S^{-1}
final Transformation3D sasInverse = shearAndScale.getInverse();
assert sasInverse != null; // as shearAndScale's det != 0
final Transformation3D rotation = m.mult(sasInverse);
if (!rotation.isIdentity()) {
result.add(OrderedPair.create(TransformFlags.Rotation,
new Matrix3X3(rotation)));
}
// Last step: decompose the scaling-and-shearing S.
// Here we assume scaling happening first, which makes S:
// S = shear * scale.
// The scale matrix is obvious and uses a,b,c. For the shearing
// we assume
// / 1 i j \
// shear = | 0 1 k |
// \ 0 0 1 /
// Multiplying gives
// / a bi cj \ / a d e \
// shear * scale = | 0 b ck | = S = | 0 b f |
// \ 0 0 c / \ 0 0 c /
// therefore:
final double i = d / b;
final double j = e / c;
final double k = f / c;
if (i != 0 || j != 0 || k != 0) {
result.add(OrderedPair.create(TransformFlags.Shearing,
new Matrix3X3(1, i, j,
0, 1, k,
0, 0, 1)));
}
if (a != 1 || b != 1 || c != 1) {
result.add(OrderedPair.create(TransformFlags.Scaling,
Transformation3D.createScaling(a, b, c)));
}
return result;
}
Decomposition of a 3X2 Matrix
This is a default method in my Transformation3D interface.
Vector2D is a 2D column vector. OrderedPair is a typed pair of values.
/**
* Decompose this transformation into simpler ones.
*
* Be prepared that there is not only one way to decompose a matrix.
*
* The algorithm will produced the following order (where some or all
* of the simple transformations might be missing if they are just
* the identity):
* <ol>
* <li>{@link TransformFlags#Translation}</li>
* <li>{@link TransformFlags#Rotation}</li>
* <li>{@link TransformFlags#Shearing}</li>
* <li>{@link TransformFlags#Scaling} (including possible mirroring)</li>
* </ol>
* @return a chain of simple transformations like translation, rotation, shearing and scaling
* (but no projection) that when multiplied in the given order recreate this
* transformation. If this represents identity the returned list will be empty.
* If the matrix has an 0 determinant it is not considered decomposable
* and {@code null} is returned.
*/
@Nullable
default List<OrderedPair<TransformFlags, Transformation2D>> decompose()
{
final double det = getDeterminant();
if (det == 0) {
Debug.warn("Decomposition of uninvertable transformation: %0", this);
return null;
}
if (det == 1 && isIdentity()) {
return Collections.emptyList();
}
final Vector2D columnX = getColumnX();
final Vector2D columnY = getColumnY();
final Vector2D columnW = getColumnW();
final List<OrderedPair<TransformFlags, Transformation2D>> result = new LinkedList<>();
// Translation is easy when it is considered to be applied last, which we do here.
if (!TrafoUtil.equals(columnW, IVector2D.ZERO)) {
result.add(OrderedPair.create(TransformFlags.Translation,
createTranslation(columnW)));
}
// Now what is left is a 2X2 matrix which is considered to be a combination
// of rotation, shearing and scaling.
// We make use of the fact that for a rotation R the following is always valid:
// R^T = R^{-1}
// We now consider the remaining 2X2 matrix M to be the combination of
// a rotation R and a shear/scale matrix S:
// M = R * S
// If we calculate
// M^T * M = (R * S)^T * (R * S)
// = (S^T * R^T) * R * S
// = S^T * R^{-1} * R * S // because R^T = R^{-1}
// M^T * M = S^T * S
// This means that we can derive S from the last equation. For this we
// assume that S is a triangle matrix in the form
// / a d \
// S = \ 0 b /
// Multiplying S^T * S results in
// / a² ad \
// S^T * S = \ ad b²+d² | = M^T * M
// Resolving this leads to the formulas for a,b, and d below.
// The squares allow to take plus or minus values for the
// scaling coefficients on the diagonal. Here we put a
// the sign (+ or -) into factor a, b is always positive.
// This selects 1 of the 4 possible solutions
// and takes care that any mirroring operation is put into the
// scaling, so all other transformations in the chain have det=1.
final Matrix2X2 m = Matrix2X2.fromColumns(columnX, columnY);
if (det == 1 && m.isIdentity()) {
return result;
}
final Matrix2X2 mTm = m.getTransposed().mult(m);
final double a = det < 0 // a cannot be 0 because of det != 0
? -Math.sqrt(mTm.getXX())
: Math.sqrt(mTm.getXX());
final double d = mTm.getXY() / a;
final double b = Math.sqrt(mTm.getYY() - d*d);
final Matrix2X2 shearAndScale = new Matrix2X2(a, d,
0, b);
// Having shear-and-scale we can calculate the rotation from
// R * S = M <=> R = M * S^{-1}
final Transformation2D sasInverse = shearAndScale.getInverse();
assert sasInverse != null; // as shearAndScale's det != 0
final Transformation2D rotation = m.mult(sasInverse);
if (!rotation.isIdentity()) {
result.add(OrderedPair.create(TransformFlags.Rotation,
new Matrix2X2(rotation)));
}
// Last step: decompose the scaling-and-shearing S.
// Here we assume scaling happening first, which makes S:
// S = shear * scale.
// The scale matrix is obvious and uses a and b. For the shearing
// we assume
// / 1 k \
// shear = | 0 1 |
// Multiplying gives
// / a bk \ / a d \
// shear * scale = | 0 b | = S = | 0 b |
// therefore:
final double k = d / b;
if (k != 0) {
result.add(OrderedPair.create(TransformFlags.Shearing,
new Matrix2X2(1, k,
0, 1)));
}
if (a != 1 || b != 1) {
result.add(OrderedPair.create(TransformFlags.Scaling,
Transformation2D.createScaling(a, b)));
}
return result;
}