新建网页 1
从四元数转换到矩阵
为了将角位移从四元数转换到矩阵形式,可以利用旋转矩阵,它能计算绕任意轴的旋转:
这个矩阵是用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的倍角公式,再代入θ:
上式是正确的,但它和其他的标准形式不同,即:
m11 = 1 - 2y2 - 2z2
实际上,还有其他的形式存在。最著名的一个形式是m11 = w2 + x2 - y2- z2,因为w2 + x2 + y2 + z2 = 1,所以这三种形式是等价的。现在回过头来看看能不能直接导出其他标准形式,第一步,n是单位向量,nx2+ny2 + nz2 = 1,则1 - nx2 = ny2 + nz2。
m11 = 1 - (1 - nx2)(2sin2(θ/2))
= 1 - (ny2 +nz2)(2sin2(θ/2))
= 1 - 2ny2sin2(θ/2) - 2nz2sin2(θ/2)
= 1 - 2y2 - 2z2
元素m22和m33可以用同样的方法求得。
让我们来看看非对角线元素,它们比对角线元素简单一些,以m12为例子:
m12 = nxny(1 - cosθ) + nzsinθ
其他非对角线元素可用同样的方法导出。
最后,给出从四元素构造的完整旋转矩阵,如公式10.23所示:
从矩阵转换到四元数
为了从旋转矩阵中抽出相应的四元数,可以直接利用公式 10.23,检查对角线元素的和(也称作矩阵的轨迹)得到:
通过使轨迹中三个元素中的两个为负,可以用类似的方法求得其他三个元素:
不幸的是,这种方法并不总是能正确工作,因为平方根的结果总是正值。(更加准确地说,没有选择正根还是负根的依据。)但是,q和-q代表相同的方位,我们能任意选择用非负根作为4个分量中的一个并仍能得到正确的四元数,只是不能对四元数的所有4个数都用这种方法。
另一个技巧是检查相对于对角线的对称位置上元素的和与差:
那么应选用四种方法中的哪个呢?似乎最简单的策略是总是先计算同一个分量,如w,然后再计算x、y、z。这伴随着问题,如果w=0,除法就没有定义;如果w非常小,将会出现数值不稳定。Shoemake建议先判断w、x、y、z中哪个最大(不用做平方根运算),根据上面的表,用矩阵对角线计算该元素,再用它计算其他三个。
下面的代码用一种非常直接的方式实现了这个方法。
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 quaternion
float w,x,y,z;
// Determine which of w, x, y, or z has the largest absolute value
float 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 division
float biggestVal = sqrt(fourBiggestSquaredMinus1 + 1.0f) * 0.5f;
float mult = 0.25f / biggestVal;
// Apply table to compute quaternion values
switch (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;
}