decompose transformation matrix
参考：https://math.stackexchange.com/questions/237369/given-this-transformation-matrix-how-do-i-decompose-it-into-translation-rotati
写自己的游戏引擎的时候遇到了这个问题，主要有以下几个点没搞清楚：
- decompose transformation matrix的原理是什么，具体应该怎么做
- 什么是shear，什么是trace
- glm::decompose的代码解析
- 正常Transform矩阵的最后一行是[0,0,0,1]，如果不是[0,0,0,1]，那意味着什么？
- 什么是transform matrix里的perspective矩阵
如何分解Transform矩阵
正常的Transform矩阵应该是这么个形式，最后一行的元素为[0,0,0,1]：
如果右下角的值为0，那么该矩阵不是Transform矩阵。如果右下角的值不为1，假设为x，那么可以整个矩阵的元素除以x，把它变成上面这个形式。
提取Translation
这个步骤很简单了，直接取[d, h, l]就是了，没啥好说的
剩下Rotation和Scale部分，看左上角的3X3矩阵即可：
提取Scale
这个步骤也不复杂，需要取前三个列向量的模，即可，如下所示：
提取Rotation
把Scale带来的影响干掉，剩下的矩阵就是旋转矩阵了，如下所示：
准确的说，如果数据正确的话，这个矩阵应该是pure rotation matrix，它有俩特点：
- 矩阵是正交矩阵(正交矩阵的行向量与列向量皆为正交的单位向量，且该矩阵的转置矩阵为其逆矩阵，正交矩阵的行列式为1或-1)
- 矩阵是个特殊的正交矩阵，它的行列式只可能是1
得到一个旋转矩阵还不够，要想在游戏引擎里使用Transform数据，还需要获得对应的欧拉角或者quaternion，所以还需要知道如何把旋转矩阵转换成欧拉角或四元数。
旋转矩阵与欧拉角的互换
参考：https://learnopencv.com/rotation-matrix-to-euler-angles/
参考：http://eecs.qmul.ac.uk/~gslabaugh/publications/euler.pdf
既然是欧拉角，那么先绕哪个轴旋转，这个是自由规定的，一共有A33，六种：X-Y-Z, X-Z-Y, Y-Z-X, Y-X-Z, Z-X-Y, Z-Y-X。
-
绕X轴旋转特定角度的矩阵为：
-
绕Y轴旋转特定角度的矩阵为：
-
绕Z轴旋转特定角度的矩阵为：
按照不同的顺序，三个相乘，得到的结果也不同，比如XYZ顺序的矩阵为：
其他顺序也会有个公式，就不贴出来了，拿这个公式，去跟上面获取的旋转矩阵进行逐个元素的等式运算，就能算出这三个角度了。
旋转矩阵与quaternion的互换
参考：https://automaticaddison.com/how-to-convert-a-quaternion-to-a-rotation-matrix/
很多文章都是直接贴公式了：
不过我还是想研究研究原理，我之前写过一篇文章分析这个问题，可以先看看怎么把quaternion转换成旋转矩阵。
可以思考一下旋转矩阵的组成方式，它本质是从单位矩阵，apply一个旋转得到的，它的3x3的矩阵部分，三条列向量分别为三个正交基。所以从quaternion转换成旋转矩阵的过程就是这个思路，用这个四元数去作用于三条单位向量即可，代码如下：
mat4 quatToMat4(const quat& q)
{
// 这里四元数乘以vec3的运算符被重载了, 具体算法不在本文的讨论范围内
vec3 r = q * vec3(1, 0, 0);
vec3 u = q * vec3(0, 1, 0);
vec3 f = q * vec3(0, 0, 1);
return mat4(r.x, r.y, r.z, 0,
u.x, u.y, u.z, 0,
f.x, f.y, f.z, 0,
0 , 0 , 0 , 1);
}
反过来，要把旋转矩阵转成quaternion，稍微麻烦一些。总之，现在有两组坐标基，一组是正常的xyz对应的坐标基，另外一组就是这三个列向量组成的坐标基，我需要获得一个q，能把第一组转换成第二组。
我可以先算任意坐标轴，旋转到第二组的四元数q1，比如我选择(0,0,1)，那么有：
// newForward为旋转矩阵的第三列
q1 = fromTo((0,0,1), newForward);
此时的q1会作用于原本的三列向量上，但此时只有newForward能够重合，所以还需要保证其他轴在应用旋转后也可以重合，所以有：
vec3 objectUp =q1 * (0,1,0);
q2 = fromTo(objectUp, newUp);//再算一个q2, fromTo会计算最近的RotationPath, 所以newForward不会改变
第三条轴是正交的，只要两轴对好了就可以了，所以最终的四元数就是：
q = q1 * q2;
也有直接套公式的办法，参考http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/，如下图所示：
相关代码示例:
// 最后看看怎么拆成四元数的, 矩阵的trace是右下斜对角线的元素总和
int i, j, k = 0;
T root, trace = Row[0].x + Row[1].y + Row[2].z;
if (trace > static_cast<T>(0))
{
root = sqrt(trace + static_cast<T>(1.0));
Orientation.w = static_cast<T>(0.5) * root;
root = static_cast<T>(0.5) / root;
Orientation.x = root * (Row[1].z - Row[2].y);
Orientation.y = root * (Row[2].x - Row[0].z);
Orientation.z = root * (Row[0].y - Row[1].x);
} // End if > 0
else
{
static int Next[3] = { 1, 2, 0 };
i = 0;
if (Row[1].y > Row[0].x) i = 1;
if (Row[2].z > Row[i][i]) i = 2;
j = Next[i];
k = Next[j];
root = sqrt(Row[i][i] - Row[j][j] - Row[k][k] + static_cast<T>(1.0));
Orientation[i] = static_cast<T>(0.5) * root;
root = static_cast<T>(0.5) / root;
Orientation[j] = root * (Row[i][j] + Row[j][i]);
Orientation[k] = root * (Row[i][k] + Row[k][i]);
Orientation.w = root * (Row[j][k] - Row[k][j]);
} // End if <= 0
shear与trace
参考：https://mathworld.wolfram.com/Shear.html
参考：https://en.wikipedia.org/wiki/Shear_matrix
参考：https://www.youtube.com/watch?v=Z2bR0Mb1Jj0&ab_channel=JoelSperanzaMath
shear翻译为剪切，几何上就是一种切变效果，如下图所示：
把单位矩阵的任意行，或者任意列，与别的行或列相加，得到的矩阵就是shear矩阵，如下图所示就是一个shear矩阵：
矩阵里提取出shear部分
正常代表Transformation的4x4的矩阵是可能含有shear部分的，但是游戏引擎里GameObject的Transform对应的矩阵正常应该是没有这个分量的，由于后面要研究的glm::decompose函数里涉及到了这个，所以提一下。
trace
trace就比较简单了，就是矩阵右下对角线上的元素和，翻译为矩阵的迹，以前线代学过
glm::decompose的代码解析
glm::decompose函数针对的是任意的4X4矩阵，矩阵代表着任意transformation，所以除了基本的RTS，它还可以提取shear和projection的部分，切边部分用vec3表示，projection部分用vec4表示，如下图所示，是一个4X4矩阵各个区间范围对应的分量：
红框对应了切变、Rotation和Scale，绿框对应的Translation，黄框对应的Projection，代码如下，其实很多子模块，前面都分析过了
template<typename T, qualifier Q>
/*GLM_FUNC_QUALIFIER*/ bool decompose(mat<4, 4, T, Q> const& ModelMatrix, vec<3, T, Q>& Scale, qua<T, Q>& Orientation, vec<3, T, Q>& Translation, vec<3, T, Q>& Skew, vec<4, T, Q>& Perspective)
{
// 首先Copy一个输入的Transform矩阵
mat<4, 4, T, Q> LocalMatrix(ModelMatrix);
// Normalize the matrix. 保证矩阵的w为1, 也就是Translation部分的w要为1
if (epsilonEqual(LocalMatrix[3][3], static_cast<T>(0), epsilon<T>()))
return false;// 如果w为0, 那么矩阵有问题
for (length_t i = 0; i < 4; ++i)
for (length_t j = 0; j < 4; ++j)
LocalMatrix[i][j] /= LocalMatrix[3][3];
// 复制Transform矩阵, 叫做PerspectiveMatrix
// perspectiveMatrix is used to solve for perspective, but it also provides
// an easy way to test for singularity of the upper 3x3 component.
mat<4, 4, T, Q> PerspectiveMatrix(LocalMatrix);
// PerspectiveMatrix的最后一行为[0,0,0,1]
for (length_t i = 0; i < 3; i++)
PerspectiveMatrix[i][3] = static_cast<T>(0);
PerspectiveMatrix[3][3] = static_cast<T>(1);
// 正常情况下, PerspectiveMatrix的行列式应该不为0
if (epsilonEqual(determinant(PerspectiveMatrix), static_cast<T>(0), epsilon<T>()))
return false;
// 判断最后一行是否为(0, 0, 0, 1), 如果不是, 则提取出Projection Matrix
// 这里看似是判断矩阵最后一列, 实际上是判断矩阵最后一行, 不为(0, 0, 0, 1)
// glm::mat4在内存里是每列占一行数据的, 所以有点绕, 比如mat3[3][2], 实际上是矩阵第3列第2行的元素
if (
epsilonNotEqual(LocalMatrix[0][3], static_cast<T>(0), epsilon<T>()) ||
epsilonNotEqual(LocalMatrix[1][3], static_cast<T>(0), epsilon<T>()) ||
epsilonNotEqual(LocalMatrix[2][3], static_cast<T>(0), epsilon<T>()))
{
// rightHandSide is the right hand side of the equation.
vec<4, T, Q> RightHandSide;
RightHandSide[0] = LocalMatrix[0][3];
RightHandSide[1] = LocalMatrix[1][3];
RightHandSide[2] = LocalMatrix[2][3];
RightHandSide[3] = LocalMatrix[3][3];
// Solve the equation by inverting PerspectiveMatrix and multiplying
// rightHandSide by the inverse. (This is the easiest way, not
// necessarily the best.)
mat<4, 4, T, Q> InversePerspectiveMatrix = glm::inverse(PerspectiveMatrix);// inverse(PerspectiveMatrix, inversePerspectiveMatrix);
mat<4, 4, T, Q> TransposedInversePerspectiveMatrix = glm::transpose(InversePerspectiveMatrix);// transposeMatrix4(inversePerspectiveMatrix, transposedInversePerspectiveMatrix);
Perspective = TransposedInversePerspectiveMatrix * RightHandSide;
// v4MulPointByMatrix(rightHandSide, transposedInversePerspectiveMatrix, perspectivePoint);
// Clear the perspective partition
LocalMatrix[0][3] = LocalMatrix[1][3] = LocalMatrix[2][3] = static_cast<T>(0);
LocalMatrix[3][3] = static_cast<T>(1);
}
else
{
// 如果矩阵最后一行, 为(0, 0, 0, 1), 正常应该都走这里
// No perspective.
Perspective = vec<4, T, Q>(0, 0, 0, 1);
}
// ================== 1. 平移部分 =====================
// Next take care of translation (easy). 直接取第三列数据的前三个值即可
Translation = vec<3, T, Q>(LocalMatrix[3]);
// 再把矩阵的Translation信息抹掉
LocalMatrix[3] = vec<4, T, Q>(0, 0, 0, LocalMatrix[3].w);
// 创建一个3x3的矩阵, 还是[列号][行号]
vec<3, T, Q> Row[3], Pdum3;
// 2. 创建3x3矩阵
// Now get scale and shear.
for (length_t i = 0; i < 3; ++i)
for (length_t j = 0; j < 3; ++j)
Row[i][j] = LocalMatrix[i][j];
// ================== 2. Scale和切边部分 =====================
// 3. 提取Scale部分, 得到一个3x3的正交矩阵
// Compute X scale factor and normalize first row.
Scale.x = length(Row[0]);// v3Length(Row[0]);
// 矩阵第一列归一化
Row[0] = detail::scale(Row[0], static_cast<T>(1));
// 如果第二列与第一列是正交的, 那么不存在切边
// Compute XY shear factor and make 2nd row orthogonal to 1st.
Skew.z = dot(Row[0], Row[1]);// 正交则Skew.z = 0
// 如果Skew.z=0, 则Row[1]不变, 否则把Row[1]改成与Row[0]正交的向量
Row[1] = detail::combine(Row[1], Row[0], static_cast<T>(1), -Skew.z);
// 去掉切边分量后, 再算scale, 再normalize
Scale.y = length(Row[1]);
Row[1] = detail::scale(Row[1], static_cast<T>(1));
Skew.z /= Scale.y;
// 同上, 如果存在切边, 则Compute XZ and YZ shears, orthogonalize 3rd row.
Skew.y = glm::dot(Row[0], Row[2]);
Row[2] = detail::combine(Row[2], Row[0], static_cast<T>(1), -Skew.y);
Skew.x = glm::dot(Row[1], Row[2]);
Row[2] = detail::combine(Row[2], Row[1], static_cast<T>(1), -Skew.x);
// Next, get Z scale and normalize 3rd row.
Scale.z = length(Row[2]);
Row[2] = detail::scale(Row[2], static_cast<T>(1));
Skew.y /= Scale.z;
Skew.x /= Scale.z;
// 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]); // v3Cross(row[1], row[2], Pdum3);
if (dot(Row[0], Pdum3) < 0)
{
for (length_t i = 0; i < 3; i++)
{
Scale[i] *= static_cast<T>(-1);
Row[i] *= static_cast<T>(-1);
}
}
// ================== 3. Rotation部分 =====================
// Now, get the rotations out, as described in the gem.
// FIXME - Add the ability to return either quaternions (which are
// easier to recompose with) or Euler angles (rx, ry, rz), which
// are easier for authors to deal with. The latter will only be useful
// when we fix https://bugs.webkit.org/show_bug.cgi?id=23799, so I
// will leave the Euler angle code here for now.
// 这里就是按照XYZ的顺序分解得到欧拉角的
// 1. 根据-sin(θ) = mat[0][2](第一列第三个元素), 算出θ值
// ret.rotateY = asin(-Row[0][2]);
// 当cosθ不为0时, 用数组的第二列第三个元素除以数组第三列第三个元素, 直接得到tan(rotateX)
// if (cos(ret.rotateY) != 0)
// {
// ret.rotateX = atan2(Row[1][2], Row[2][2]);
// // 当cosθ不为0时, 根据矩阵第一列算得rotateZ
// ret.rotateZ = atan2(Row[0][1], Row[0][0]);
// }
// else
// {
// ret.rotateX = atan2(-Row[2][0], Row[1][1]);
// ret.rotateZ = 0;
// }
// 最后看看怎么拆成四元数的, 矩阵的trace是右下斜对角线的元素总和
int i, j, k = 0;
T root, trace = Row[0].x + Row[1].y + Row[2].z;
if (trace > static_cast<T>(0))
{
root = sqrt(trace + static_cast<T>(1.0));
Orientation.w = static_cast<T>(0.5) * root;
root = static_cast<T>(0.5) / root;
Orientation.x = root * (Row[1].z - Row[2].y);
Orientation.y = root * (Row[2].x - Row[0].z);
Orientation.z = root * (Row[0].y - Row[1].x);
} // End if > 0
else
{
static int Next[3] = { 1, 2, 0 };
i = 0;
if (Row[1].y > Row[0].x) i = 1;
if (Row[2].z > Row[i][i]) i = 2;
j = Next[i];
k = Next[j];
root = sqrt(Row[i][i] - Row[j][j] - Row[k][k] + static_cast<T>(1.0));
Orientation[i] = static_cast<T>(0.5) * root;
root = static_cast<T>(0.5) / root;
Orientation[j] = root * (Row[i][j] + Row[j][i]);
Orientation[k] = root * (Row[i][k] + Row[k][i]);
Orientation.w = root * (Row[j][k] - Row[k][j]);
} // End if <= 0
return true;
}
