Flag Engine（动画系统）学习笔记（四）——“骨骼姿势”前篇之四元数

2021SC@SDUSC

在上一节中学习了关于骨骼的数据结构，接下来我们会继续学习在动画系统中，骨骼是如何相对某参考系进行收收放，平移和旋转的。而这也被称为骨骼（关节）的姿势。

关节的姿势是通过SQT数据结构，一个4*3或4*4的矩阵表示（缩放，scale；四元数旋转，quaternion；矢量平移（方向，距离（坐标）），translation）。其中旋转是最复杂的一部分。

矩阵旋转
优点：
旋转轴可以是任意向量；
缺点：
旋转其实只需要知道一个向量+一个角度，一共4个值的信息，但矩阵法却使用了16个元素；
而且在做乘法操作时也会增加计算量，造成了空间和时间上的一些浪费；

欧拉旋转
优点：
很容易理解，形象直观；
表示更方便，只需要3个值（分别对应x、y、z轴的旋转角度）；但按我的理解，它还是转换到了3个3*3的矩阵做变换，效率不如四元数；
缺点：
之前提到过这种方法是要按照一个固定的坐标轴的顺序旋转的，因此不同的顺序会造成不同的结果；
会造成万向节锁（Gimbal Lock）的现象。这种现象的发生就是由于上述固定坐标轴旋转顺序造成的。理论上，欧拉旋转可以靠这种顺序让一个物体指到任何一个想要的方向，但如果在旋转中不幸让某些坐标轴重合了就会发生万向节锁，这时就会丢失一个方向上的旋转能力，也就是说在这种状态下我们无论怎么旋转（当然还是要原先的顺序）都不可能得到某些想要的旋转效果，除非我们打破原先的旋转顺序或者同时旋转3个坐标轴。这里有个视频可以直观的理解下；
由于万向节锁的存在，欧拉旋转无法实现球面平滑插值；

四元数旋转
优点：
可以避免万向节锁现象；
只需要一个4维的四元数就可以执行绕任意过原点的向量的旋转，方便快捷，在某些实现下比旋转矩阵效率更高；
可以提供平滑插值；
缺点：
比欧拉旋转稍微复杂了一点点，因为多了一个维度；
理解更困难，不直观；

而在学习SQT数据结构之前，先学习一下四元数旋转。

四元数是对复数的扩充，它使用三个虚部i,j,k它们的关系如下:

i² = j² = k² = –1

四元数形式:

[w  v]或[ w   (x    y     z) ].

四元数定义了一个复数

w + xi + yj + zk

       四元数只能用来代替矩阵保存旋转信息，平移无法代替，旋转我们有 矩阵，欧拉角，四元数，

四元数是里面最复杂的一个但它可以在保证效率的同时减小矩阵1/4的内存占有量，同时又能避免欧拉角的万向锁问题，

无法用语言表达清楚万向锁的现象只有碰到过才能理解。。。

四元数和轴——角对:

设O为旋转角度，N为旋转轴(绕任意轴n旋转O度):

Q = [cos(O/2)  sin(O/2)N] = [ cos(O/2)     ( sin(O/2)Nx     sin(O/2)Ny      sin(O/2)Nz ) ]

负四元数:

-Q = –[w    (x y z)]  = [-w   (-x  -y  -z)] = –[w  v] = [-w -v]

    Q和-Q代表的是相同的角位移，如果我们将O加上360°的倍数，不会改变Q代表的角位移，但它使Q的四个分量都变负了。

因此，3D中的任意角位移都有两种不同的四元数表示方法，它们互为负数。

单位四元数:

几何上存在两个“单位”四元数，它们没有角位移， [1, 0]和[-1,0] 0代表0向量，让我们看看原因：

当O为360°的偶数倍时，有第一种形式，cos(O/2) = 1;

当O为360°的奇数倍时，有第二种形式，cos(0/2) =-1;

在这两种形势下都有sin(O/2) = 0,所以和旋转轴N是无关的,他的意义是，当旋转角O是360°的整数倍时，方位没有改变，并且旋转轴也是关紧要的.

注意:

       数学上只有一个单位四元数[1, 0]，用四元数q乘以单位四元数[1,0]，结果认识q。任意四元数q乘以另一个“几何单位”四元数[-1,0]时得到-q。几何上，因为q和-q代表的角位移相同，可

认为结果是相同。但在数学上，q和-q不相等，所以[-1,0]斌不是“真正”的单位四元数。

四元数的模:

公式:

    || q || = || [w (x y z)] || = √(w² + x² + y² + z²)

    所以

   || q || = √( cos²(O/2) + ( sin(O/2) ||n||)² )

   当n为单位向量时:

             = √( cos²(O/2) + sin²(O/2) * 1 )

    应用三角形公式 sin²(x) + cos²(x) = 1

             = √( cos²(O/2) + sin²(O/2) )

             = √1 = 1

如果为了用四元数来表示方位，我们仅仅使用符合这个规则的单位四元数。非单位四元数要查阅其他资料

四元数共轨和逆:

共轨:

四元数共轨记做 q*可以通过让四元数向量部分变负来获得:

q*  = [w v]* = [w -v]

    = [w (-x –y -z)]

逆:

用q^-1 表示逆

q^-1 = q* / || q ||

四元数乘以自己的逆等于“单位四元数”[1 0]

如果我们使用单位四元数是单位四元数那么就可以得到 q^-1 = q*. 所以单位四元数的共轨等于他的逆

四元数叉乘:

四元数乘法的标准定义:

[W1 V1] *[W2 V2]叉乘后还是四元数

(W1+X1i + Y1j + Z1k)(W2 + X2i + y2j + Z2K)

= [W1W2 - V1·V2  W1V2+W2V1+V2*V1]

叉乘满足结合律，但不满足交换律.

(a*b)*c = a*(b*c)

ab ≠ ba

叉乘的模:

||Q1 * Q2|| = ||Q1|| * ||Q2||

这个结论保证了两个单位四元数相乘还是单位四元数.

四元数逆的乘积:

(a*b)^1 = b^-1  * a^-1

四元数和点：

把一个标准3D点(x,y,z)扩展到四元数空间:P = [0, (x,y,z)]即可，一般情况下它不会是单位四元数

设Q为旋转四元数Q=[cos(O/2),nsin(O/2)],n为单位向量，O为旋转角度,那么P绕n旋转O度就是：

P` = Q * P * Q^-1

四元数的乘法和3D向量旋转的对应关系更多是理论上的，实际上，它几乎和把四元数转换成矩阵形式然后在用矩阵乘以向量时间是一样的.

多次旋转:

P` = B * (A * P * A^-1)B^-1

    = (B * A) * P * (A * B)^-1

这个形式是标准定义但我们不经常使用因为不好看,下面是经常使用的这个影响到了叉乘本身:

红色部分是改变后所影响的

[W1 V1] *[W2 V2]

(W1+X1i + Y1j + Z1k)(W2 + X2i + y2j + Z2K)

= [W1W2 - V1·V2 W1V2+W2V1+V2*V1]

相应的 四元数与3D向量的对应关系也随着改变:

P` = Q^-1 * P * Q

p` = B^-1 * (A^-1 * P * A) * B

     = (A * B)^-1 * P * (A * B)

四元数的差:

我们看：

A * D = B

两边都乘以 A^-1

A^-1 * A = A^-1 * B

因为 A^-1 * A = 1

所以 D = A^-1 * B

这就所谓的 “差”,被定义为一个方位到另一个方位的角位移,顺序是不能错的,从左向右.

四元数的点乘:

点乘那个很简单: Q1 · Q2 = [w1 v1] ·[w2 v2] = w1 · w2 + v1 · v2

点乘结是一个标量这个和向量点乘一样,对于单位四元数 a 和 b 有-1≤  a · b ≤1。

因为 a·b = –( a · –b)，多以b和-b代表的角位移是一样的。

几何解释，四元数点乘a · b的绝对值越大，a 和 b代表的角位移越相似，这一点和 向量点乘类似.

四元数的对数，指数，标量乘运算:

对数：

指数:

标量相乘:

Kq = k[w v] = [kw kv]

从矩阵到欧拉角:

将角位移从矩阵形式转换到欧拉角需要考虑以下几点：

（1）必须清楚矩阵代表什么旋转：物体 -- 惯性还是 惯性 -- 物体，这里讨论使用惯性 -- 物体矩阵的技术，物体 -- 惯性矩阵转换成欧拉角的过程与之类似。

（2）对任意给定角位移，存在无穷多个欧拉角可以用来表示它。因为"别名"问题，这里讨论的技术总是返回"限制欧拉角"，heading和bank的范围±180°，pitch的范围为±90°。

（3）矩阵可能是病态的，我们必须忍受浮点数精度的误差。有些矩阵还包括除旋转外的变换，如缩放、镜像等。这里只讨论工作在旋转矩阵上的变换。   

  //Listing 10.3: Extracting Euler angles from an inertial-to-object rotation matrix// Assume the matrix is stored in these variables:float m11,m12,m13;float m21,m22,m23;float m31,m32,m33;// We will compute the Euler angle values in radians and store them here:float h,p,b;// Extract pitch from m23, being careful for domain errors with asin(). We could have// values slightly out of range due to floating point arithmetic.float sp = –m23;if (sp <= –1.0f) {p = –1.570796f; // –pi/2} else if (sp >= 1.0) {p = 1.570796; // pi/2} else {p = asin(sp);}// Check for the Gimbal lock case, giving a slight tolerance// for numerical imprecisionif (sp > 0.9999f) {// We are looking straight up or down.// Slam bank to zero and just set headingb = 0.0f;h = atan2(–m31, m11);} else {// Compute heading from m13 and m33h = atan2(m13, m33);// Compute bank from m21 and m22b = atan2(m21, m22);}

从四元数到矩阵:

为了将角位移从四元数转换到矩阵形式，可以利用旋转矩阵，它能计算绕任意轴的旋转：

这个矩阵是用n和θ表示的，但四元数的分量是：

w = cos(θ/2)

x = nx sin(θ/2)

y = ny sin(θ/2)

z = nz sin(θ/2)

让我们看看能否将矩阵变形以代入w、x、y、z，矩阵的9个元素都必须这样做。幸运的是，这个矩阵的结构非常好，一旦解出对角线上的一个元素，其他元素就能用同样的方法求出。同样，非对角线元素之间也是彼此类似的。

考虑矩阵对角线上的元素，我们将完整地解出m11，m22和m33解法与之类似：

m11 = nx2(1 - cosθ) + cosθ

我们将从上式的变形开始，变形方法看起来像是在绕圈子，但你马上就能理解这样做的目的：

现在需要消去cosθ项，而代之以包含cosθ/2或sinθ/2的项，因为四元数的元素都是用它们表示的，像以前那样，设α=θ/2，先用α写出cos的倍角公式，再代入θ

// Listing 10.4: Converting a rotation matrix to a quaternion// Input matrix:float m11,m12,m13;float m21,m22,m23;float m31,m32,m33;// Output quaternionfloat w,x,y,z;// Determine which of w, x, y, or z has the largest absolute valuefloat fourWSquaredMinus1 = m11 + m22 + m33;float fourXSquaredMinus1 = m11 – m22 – m33;float fourYSquaredMinus1 = m22 – m11 – m33;float fourZSquaredMinus1 = m33 – m11 – m22;int biggestIndex = 0;float fourBiggestSquaredMinus1 = fourWSquaredMinus1;if (fourXSquaredMinus1 > fourBiggestSquaredMinus1) {fourBiggestSquaredMinus1 = fourXSquaredMinus1;biggestIndex = 1;}if (fourYSquaredMinus1 > fourBiggestSquaredMinus1) {fourBiggestSquaredMinus1 = fourYSquaredMinus1;biggestIndex = 2;}if (fourZSquaredMinus1 > fourBiggestSquaredMinus1) {fourBiggestSquaredMinus1 = fourZSquaredMinus1;biggestIndex = 3;}// Perform square root and divisionfloat biggestVal = sqrt(fourBiggestSquaredMinus1 + 1.0f) * 0.5f;float mult = 0.25f / biggestVal;// Apply table to compute quaternion valuesswitch (biggestIndex) {case 0:w = biggestVal;x = (m23 – m32) * mult;y = (m31 – m13) * mult;z = (m12 – m21) * mult;break;case 1:x = biggestVal;w = (m23 – m32) * mult;y = (m12 + m21) * mult;z = (m31 + m13) * mult;break;case 2:y = biggestVal;w = (m31 – m13) * mult;x = (m12 + m21) * mult;z = (m23 + m32) * mult;break;case 3:z = biggestVal;w = (m12 – m21) * mult;x = (m31 + m13) * mult;y = (m23 + m32) * mult;break;}  

从欧拉角到四元数:

为了将角位移从欧拉角转换到四元数，可以使用从欧拉角构造矩阵类似的方法。先将这三个旋转分别转换为四元数，这是一个简单的运算。再将这三个四元数连接成一个四元数。和矩阵一样，有两种情况需要考虑，第一种是惯性 -- 物体四元数，第二种是物体-- 惯性四元数。因为它们互为共轭关系，所以我们只推导惯性--物体四元数。

设欧拉角为变量h、p、b，设h、p、b分别绕轴y、x、z旋转的四元数。记住，使用负旋转量，因为它们指定坐标系中的旋转角度。

用正确的顺序连接它们得到公式10.24：

（记住，四元数乘法定义是按旋转的顺序从左向右乘。）

物体--惯性四元数是惯性--物体四元数的共轭，见公式10.25：

从四元数到欧拉角:

根据前面的公式发现：

  // Use global variables for input and outputfloat w,x,y,z;float h,p,b;// Extract sin(pitch)float sp = –2.0f * (y*z + w*x);// Check for Gimbal lock, giving slight tolerance for numerical imprecisionif (fabs(sp) > 0.9999f) {// Looking straight up or downp = 1.570796f * sp; // pi/2// Compute heading, slam bank to zeroh = atan2(–x*z – w*y, 0.5f – y*y – z*z);b = 0.0f;} else {// Compute anglesp = asin(sp);h = atan2(x*z – w*y, 0.5f – x*x – y*y);b = atan2(x*y – w*z, 0.5f – x*x – z*z);}// --惯性四元数转换到欧拉角，所用的代码和上面非常类似。只是将x、y、z值变负，因为物体--惯性四元数是惯性--物体四元数的共轭。//Listing 10.6: Converting an object-to-inertial quaternion to Euler angles// Extract sin(pitch)float sp = –2.0f * (y*z – w*x);// Check for Gimbal lock, giving slight tolerance for numerical imprecisionif (fabs(sp) > 0.9999f) {// Looking straight up or downp = 1.570796f * sp; // pi/2// Compute heading, slam bank to zeroh = atan2(–x*z + w*y, 0.5f – y*y – z*z);b = 0.0f;} else {// Compute anglesp = asin(sp);h = atan2(x*z + w*y, 0.5f – x*x – y*y);b = atan2(x*y + w*z, 0.5f – x*x – z*z);}