• 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);
    }
    
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11

    反过来,要把旋转矩阵转成quaternion,稍微麻烦一些。总之,现在有两组坐标基,一组是正常的xyz对应的坐标基,另外一组就是这三个列向量组成的坐标基,我需要获得一个q,能把第一组转换成第二组。

    我可以先算任意坐标轴,旋转到第二组的四元数q1,比如我选择(0,0,1),那么有:

    // newForward为旋转矩阵的第三列
    q1 = fromTo((0,0,1), newForward);
    
    • 1
    • 2

    此时的q1会作用于原本的三列向量上,但此时只有newForward能够重合,所以还需要保证其他轴在应用旋转后也可以重合,所以有:

    vec3 objectUp =q1 * (0,1,0);
    q2 = fromTo(objectUp, newUp);//再算一个q2, fromTo会计算最近的RotationPath, 所以newForward不会改变
    
    • 1
    • 2

    第三条轴是正交的,只要两轴对好了就可以了,所以最终的四元数就是:

    q = q1 * q2;
    
    • 1

    也有直接套公式的办法,参考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
    
    
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30

    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;
    }
    
    • 1
    • 2
    • 3
    • 4
    • 5
    • 6
    • 7
    • 8
    • 9
    • 10
    • 11
    • 12
    • 13
    • 14
    • 15
    • 16
    • 17
    • 18
    • 19
    • 20
    • 21
    • 22
    • 23
    • 24
    • 25
    • 26
    • 27
    • 28
    • 29
    • 30
    • 31
    • 32
    • 33
    • 34
    • 35
    • 36
    • 37
    • 38
    • 39
    • 40
    • 41
    • 42
    • 43
    • 44
    • 45
    • 46
    • 47
    • 48
    • 49
    • 50
    • 51
    • 52
    • 53
    • 54
    • 55
    • 56
    • 57
    • 58
    • 59
    • 60
    • 61
    • 62
    • 63
    • 64
    • 65
    • 66
    • 67
    • 68
    • 69
    • 70
    • 71
    • 72
    • 73
    • 74
    • 75
    • 76
    • 77
    • 78
    • 79
    • 80
    • 81
    • 82
    • 83
    • 84
    • 85
    • 86
    • 87
    • 88
    • 89
    • 90
    • 91
    • 92
    • 93
    • 94
    • 95
    • 96
    • 97
    • 98
    • 99
    • 100
    • 101
    • 102
    • 103
    • 104
    • 105
    • 106
    • 107
    • 108
    • 109
    • 110
    • 111
    • 112
    • 113
    • 114
    • 115
    • 116
    • 117
    • 118
    • 119
    • 120
    • 121
    • 122
    • 123
    • 124
    • 125
    • 126
    • 127
    • 128
    • 129
    • 130
    • 131
    • 132
    • 133
    • 134
    • 135
    • 136
    • 137
    • 138
    • 139
    • 140
    • 141
    • 142
    • 143
    • 144
    • 145
    • 146
    • 147
    • 148
    • 149
    • 150
    • 151
    • 152
    • 153
    • 154
    • 155
    • 156
    • 157
    • 158
    • 159
    • 160
    • 161
    • 162
    • 163
    • 164
    • 165
    • 166
    • 167
    • 168
    • 169
    • 170
    • 171
    • 172
    • 173
    • 174
    • 175
    • 176
    • 177
    • 178
    • 179
    • 180
    • 181
    • 182
    • 183

  • 相关阅读:
    架构每日一学 1:成为一名架构师,你必须具有“战略意图”
    lua学习笔记
    ElasticSearch三种分页对比
    12、Mybatis框架-1
    MySQL 主从复制与读写分离
    Mybatis一级二级缓存
    机器视觉自动数据标注方法
    计算机网络——物理层(数字传输系统)
    go语言 sync.Map原理
    学习JAVA的二十二天(基础)
  • 原文地址:https://blog.csdn.net/alexhu2010q/article/details/126408366