原文:Introduction to 3D Game Programming with DirectX 12 学习笔记之 --- 第二十二章：四元数（QUATERNIONS）
学习目标
- 回顾复数，以及复数相乘如何在平面上表达旋转；
- 理解四元数以及它的运算；
- 理解单位四元数如何表达3D旋转；
- 学习如何转换旋转变量的表达；
- 学习如何对单位四元数线性差值，并且理解它等价于几何上的3D角度差值；
- 熟悉DirectX Math库中的四元数类和操作。
1 回顾复数
四元数可以看做是一个复数，所以我们先要回顾一下复数。本章的主要目标是展示一个复数P乘以一个单位复数，得到的是旋转后的P。
1.1 定义
有很多种方法介绍复数，我们采用下面的方法，将它想象成一个2D点或者向量。
一对排序的实数z = (a, b)是一个复数。第一个部分叫实数部分，第二个部分是虚数部分。并且加减乘除定义如下：
并且很容易证明实数的算术性质对复数也依然有效（交换律，结合律，分配率）（练习1）；
如果一个复数形式是(x, 0)，那么它就是通过实数x定义，然后写成(x, 0)；那么任何实数都可以是一个复数，其虚数部分是0；观察一个实数乘以一个复数x(a, b) = (x, 0)(a, b) = (xa, xb) = (a, b)(x, 0) = (a, b)x，这个形式可以让我们回忆起变量-向量的乘法。
我们定义单位虚数i = (0, 1)，然后使用定义的复数乘法，i^2 = (0, 1)(0, 1) = (−1, 0) = −1，所以i是方程x^2 = −1的解。
复数z = (a, b)的复数共轭表示为z‾\overline{z}z，并且z‾\overline{z}z= (a, −b)；一个简单的记住复数除法公式的方法是：分子和分母乘以分母的共轭，这样分母就成为了一个实数：
下面展示一个复数(a, b)可以写为a + ib，我们有a = (a, 0), b = (b, 0)和i = (0, 1)，所以：
使用a + ib形式，我们可以重做加减乘除：
并且在这种形式下，z = a + ib的共轭复数为a - ib。
1.2 几何解释
我们将复数a + ib = (a, b)理解为几何上的2D向量或者点（在复平面）；复数相加匹配到向量的相加：
复数的绝对值或者长度由向量的长度来定义：
我们将长度为1的复数称之为单位复数：
1.3 极坐标表示和旋转
因为复数可以看做是2D复平面的点或者向量，所以我们也可以将它表示到极坐标：
后面的等式就是复数a + ib的极坐标表示。
令z1 = r1(cosθ1 + isinθ1), z2 = (cosθ2 + isinθ2)，那么：
这里使用了三角定理：
所以几何上，z1z2的乘积表示长度为r1r2的向量旋转了θ1 + θ2角度；如果其中r2为1，那么乘积就表示将z1旋转了θ2角度。如上图，所以复数乘以单位复数，表示把前面的复数旋转。
2 四元数代数
2.1 基本运算的定义
四个有序实数q = (x, y, z, w) = (q1, q2, q3, q4)是一个四元数，通常简写为q = (u, w) = (x, y, z, w)，并且我们称u = (x, y, z)为虚数部分，w为实数部分。那么加减乘除定义如下：
乘法的定义看起来会比较奇怪，但是这些运算是定义，所以我们可以定义为我们想要的形式，并且这些形式很有用。矩阵的乘法定义一开始也看起来很奇怪，但是结果是它们很有用。
令p = (u, p4) = (p1, p2, p3, p4)并且q = (v, q4) = (q1, q2, q3, q4)，那么u × v =(p2q3 – p3q2, p3q1 – p1q3, p1q2 – p2q1)并且u·v = p1q1 + p2q2 + p3q3。那么组件形式下，四元数的乘积r = pq是：
可以写成矩阵乘法形式：
如果你偏爱行向量形式，取其转置矩阵：
2.2 特殊乘积
令i = (1, 0, 0, 0), j = (0, 1, 0, 0), k = (0, 0, 1, 0)为四元数，那么我们有一些特殊乘积，会让我们回忆起叉积：
这些等式是直接从四元数乘积公式得来的，比如：
2.3 特性
四元数的乘法不具备交换律，下图证明ij = −ji。四元数乘积具有结合律；所以四元数可以联想为矩阵的乘积，具有结合律，不具有交换律。四元数e = (0, 0, 0, 1)用以乘法单位：
四元数具有乘法分配律p(q + r) = pq + pr 和 (q + r)p = qp + rp。
2.4 转化
我们将实数、向量与四元数通过下面的方式来关联。令s是实数x，u = (x, y, z)是向量，那么：
可以说任意实数是一个有0向量的四元数，任意向量是一个具有0实数的四元数；另外单位四元数为1 = (0, 0, 0, 1)；一个四元数具有0实数称之为纯四元数（pure quaternion）。
使用四元数乘法的定义，一个实数乘以四元数是标量相乘，具有交换律：
2.5 共轭和标范数
四元数q = (q1, q2, q3, q4) = (u, q4)的共轭由q定义：
也就是直接将虚数部位取反；相比于共轭复数，它具有下面特性：
其中q + q和qq* = q*q等于实数。
四元数的范数（长度），定义为：
范数为1的四元数为单位四元数，范数具有下面的特性：
特性2表示2个单位四元数相乘，依然是单位四元数；如果||p|| = 1，那么||pq|| = ||q||。
共轭和范数的特性可以直接通过定义推导出来，比如：
2.6 逆（Inverses）
因为矩阵和四元数乘法不具有交换律，所以不能直接定义除法运算。但是每个非0四元数具有逆，令q = (q1, q2, q3, q4) = (u, q4)是一个非0的四元数，那么它的逆通过q−1q^{-1}q−1来定义：
很容易证明它是复数的逆：
可以看出如果q是单位四元数，那么∣∣q∣∣2=1||q||^2 = 1∣∣q∣∣2=1，并且q−1=q∗q^{-1} = q^*q−1=q∗。
并且符合下面的特性：
2.7 极坐标表达
如果q = (q1, q2, q3, q4) = (u, q4)是单位四元数，那么：
上图表示，对于θ∈[0, π]，q4 = cosθ，根据三角定义sin2θ + cos2θ = 1：
所以
现在求单位向量：
所以u = sinθn，现在我们可以写出单位四元数q = (u, q4)的极坐标表达，其中n是单位向量：
如果我们将−θ带入等式中的θ，只是取反了向量部分：
下节将会介绍，n表示旋转的轴向。
3 单位四元数和旋转
3.1 旋转运算
令q = (u, w)是一个单位四元数并且令v是一个3D点或者向量，然后我们认为v是一个纯四元数p = (v, 0)。当q是一个单位四元数时，q−1=q∗q^{−1} = q^*q−1=q∗。回顾四元数乘法公式：
现在考虑下面的乘法：
对其长度稍作简化，我们把实数部分和向量部分分开计算。我们做下面的符号替换：
实数部分：
其中u · (v × u) = 0，因为根据叉积的定义，(v × u)是和u正交的。
虚数部分：
其中对u × (u × v)应用了乘法定义：a × (b × c) = (a · c)b − (a · b)c。所以：
计算的结果是一个向量或者点，其实数部分为0。所以随后的等式中，我们放弃实数部分。
因为q是一个单位四元数，所以可以写为：
带入上面公式后：
为了进一步简化，我们带入三角定义：
对比第三章的旋转公式，我们发现它和旋转公式基本一致，它将向量v沿着n轴旋转2θ度。
所以我们定义四元数旋转运算：
所以如果你要沿着n轴旋转θ度，那么你可以构建对于的旋转四元数：
然后应用到旋转公式中Rq(v)R_q(v)Rq(v)。
3.2 四元数旋转运算到矩阵
令q = (u, w) = (q1, q2, q3, q4)是一个单位四元数，根据之前的公式，可以得到：
上面公式的三个部分可以分别写出矩阵形式：
将它们相加：
根据单位四元数的特性（各分组件平方的和为1），做下面的简化：
最后矩阵可以写为：
3.3 矩阵到四元数旋转运算
给出一个旋转矩阵：
我们希望找到四元数q = (q1, q2, q3, q4)，我们的策略是，设置矩阵如下：
然后求解q1, q2, q3, q4；
首先将对角线上的元素相加（最终一个矩阵）：
然后组合对角相反的元素来求解q1, q2, q3：
如果q4 = 0，那上面这些公式就无意义，所以我们要找到R的最大对角元素来除，并且选择矩阵元素的其他组合。加入R11是最大的对角：
如果假设R22 或者 R33为最大对角线，计算模式类似。
3.4 组成
假设p和q是单位四元数，并且对于旋转运算为Rp 和 Rq，令v′=Rp(v)v^{'} = R_p(v)v′=Rp(v)，那么组合：
因为p和q都是单位四元数，pq的乘积也是单位四元数||pq|| = ||p||||q|| = 1；所以pq也表示旋转；也就是说说得到的旋转为：Rq(Rp(v))R_q(R_p(v))Rq(Rp(v))。
4 四元数的插值
因为四元数是由4个实数组成的，所以几何上，可以把它看成是一个4D向量，单位四元数是4D在单位4D球体表面，除了叉积（只定义了3D向量）。特别的，点积也支持四元数，令p = (u, s)并且q = (v, t)，那么：
其中θ是两个四元数之间的夹角，如果p和q是单位四元数，那么p·q = cosθ，所以点积可以帮助我们考虑2个四元数之间的夹角。
出于动画考虑，我们需要在两个方向之间进行插值，为了插值四元数，我们需要在单位球体上进行弧度差值，所以也需要在单位四元数上差值。为何推导出公式，如下图所示：我们需要在a和b中间差值tθ。我们需要找到权重c1和c2支持p = c1a + c2b，其中||p|| = ||a|| = ||b||。我们对两个未知项创建两个等式：
然后可以导出下面的矩阵：
考虑到上面的矩阵等式Ax = b，其中A是可逆的，所以根据克莱姆法则xi=detAidetAx_i = \frac{detA_i}{detA}xi=detAdetAi，其中AiA_iAi是通过交换A中第i列的向量到b，所以：
根据三角毕达哥斯拉定义和加法公式，我们可以得出：
所以：
并且：
所以我们定义出球体差值公式：
如果将单位四元数看成4D向量的话，我们就可以求解四元数之间的夹角：θ = arccos(a · b)。
如果a和b之间的夹角趋近于0，sinθ趋近于0，那么上面公式中的除法就会引发问题，会导致无限大的结果。这种情况下，对两个四元数进行线性差值，并标准化结果，就是对小θ的一个很好的近似：
观察下图，线性差值是通过将四元数差值投影回单位球体，其结果是一个非线性速率的旋转。所以如果你对大角度使用线性差值的话，旋转的速度会时快时慢。
我们现在支持一个四元数有趣的特性，(sq)= sq并且标量-四元数的乘积是具有交换律的，所以我们可以得出：
我们得出q和-q表示的相同的旋转，也可以通过其他方式来证明，如果
Rq表示围绕n旋转θ，R-q表示围绕-n旋转2π − θ。在几何上，一个在4D单位球体上的单位四元数和它的极坐标相反值−q代表的是相同的方向。下图可以看出，这两个旋转到了相同的位置，只是一个旋转了小角度，另一个旋转了大的角度：
所以b和-b表达了相同的方向，我们有2个选择来差值：slerp(a, b, t) 或者 slerp(a, −b, t)。其中一个是从更小的角度直接旋转；另一个是从更大的角度来旋转。如下图所示，选择哪个旋转基于哪个旋转在单位球体上的弧度：选择小弧度代表选择了更直接的路径，选择更长的弧度代表对物体有额外更多的旋转[Eberly01]。
[Watt92]如果要在单位球面上找到四元数最短旋转弧度，我们可以比较||a – b||2和||a – (−b)||2 = ||a + b||2。如果 ||a + b||2 < ||a – b||2我们就选择-b，因为-b更接近a：
// Linear interpolation (for small theta).
public static Quaternion LerpAndNormalize(Quaternion p, Quaternion q, float s)
{
// Normalize to make sure it is a unit quaternion.
return Normalize((1.0f - s)*p + s*q);
}
public static Quaternion Slerp(Quaternion p, Quaternion q, float s)
{
// Recall that q and -q represent the same orientation, but
// interpolating between the two is different: One will take the
// shortest arc and one will take the long arc. To find
// the shortest arc, compare the magnitude of p-q with the
// magnitude p-(-q) = p+q.
if(LengthSq(p-q) > LengthSq(p+q))
q = -q;
float cosPhi = DotP(p, q);
// For very small angles, use linear interpolation.
if(cosPhi > (1.0f - 0.001))
return LerpAndNormalize(p, q, s);
// Find the angle between the two quaternions.
float phi = (float)Math.Acos(cosPhi);
float sinPhi = (float)Math.Sin(phi);
// Interpolate along the arc formed by the intersection of the 4D
// unit sphere and the plane passing through p, q, and the origin of
// the unit sphere.
return ((float)Math.Sin(phi*(1.0- s))/sinPhi)*p + ((float)Math.Sin(phi*s)/sinPhi)*q;
}
5 DIRECTX MATH四元数函数
DirectX数学库支持四元数。因为四元数的数据是4个实数，所以使用XMVECTOR类型类保存四元数。下面是通用的函数：
// Returns the quaternion dot product Q1·Q2.
XMVECTOR XMQuaternionDot(XMVECTOR Q1, XMVECTOR Q2);
// Returns the identity quaternion (0, 0, 0, 1).
XMVECTOR XMQuaternionIdentity();
// Returns the conjugate of the quaternion Q.
XMVECTOR XMQuaternionConjugate(XMVECTOR Q);
// Returns the norm of the quaternion Q.
XMVECTOR XMQuaternionLength(XMVECTOR Q);
// Normalizes a quaternion by treating it as a 4D vector.
XMVECTOR XMQuaternionNormalize(XMVECTOR Q);
// Computes the quaternion product Q1Q2.
XMVECTOR XMQuaternionMultiply(XMVECTOR Q1, XMVECTOR Q2);
// Returns a quaternions from axis-angle rotation representation.
XMVECTOR XMQuaternionRotationAxis(XMVECTOR Axis, FLOAT Angle);
// Returns a quaternions from axis-angle rotation representation, where the axis
// vector is normalized—this is faster than XMQuaternionRotationAxis.
XMVECTOR XMQuaternionRotationNormal(XMVECTOR NormalAxis,FLOAT Angle);
// Returns a quaternion from a rotation matrix.
XMVECTOR XMQuaternionRotationMatrix(XMMATRIX M);
// Returns a rotation matrix from a unit quaternion.
XMMATRIX XMMatrixRotationQuaternion(XMVECTOR Quaternion);
// Extracts the axis and angle rotation representation from the quaternion Q.
VOID XMQuaternionToAxisAngle(XMVECTOR *pAxis, FLOAT *pAngle, XMVECTOR Q);
// Returns slerp(Q1, Q2, t)
XMVECTOR XMQuaternionSlerp(XMVECTOR Q0, XMVECTOR Q1, FLOAT t);
6 旋转Demo
本章中的Demo，我们在简单的场景中运动一个骷髅头。位置、方形和缩放都做动画。我们用四元数来表达骷髅的方向，然后使用球面差值来对方向差值。使用线性差值对位置和缩放差值。它是对下一章中的角色动画做预热。
我们使用关键帧系统对骷髅做动画：
struct Keyframe
{
Keyframe();
˜Keyframe();
float TimePos;
XMFLOAT3 Translation;
XMFLOAT3 Scale;
XMFLOAT4 RotationQuat;
};
动画是一些列通过实践来排序的关键帧:
struct BoneAnimation
{
float GetStartTime()const;
float GetEndTime()const;
void Interpolate(float t, XMFLOAT4X4& M)const;
std::vector<Keyframe> Keyframes;
};
GetStartTime函数用来返回第一个帧的时间；GetEndTime函数返回最后一个关键帧的时间。它对于动画什么时候结束很有用，我们可以停止动画。
现在有了一个关键帧列表，对于每两个帧之间使用插值计算：
void BoneAnimation::Interpolate(float t, XMFLOAT4X4& M)const
{
// t is before the animation started, so just return the first key frame.
if( t <= Keyframes.front().TimePos )
{
XMVECTOR S = XMLoadFloat3(&Keyframes.front().Scale);
XMVECTOR P = XMLoadFloat3(&Keyframes.front().Translation);
XMVECTOR Q = XMLoadFloat4(&Keyframes.front().RotationQuat);
XMVECTOR zero = XMVectorSet(0.0f, 0.0f, 0.0f, 1.0f);
XMStoreFloat4x4(&M, XMMatrixAffineTransformation(S, zero, Q, P));
}
// t is after the animation ended, so just return the last key frame.
else if( t >= Keyframes.back().TimePos )
{
XMVECTOR S = XMLoadFloat3(&Keyframes.back().Scale);
XMVECTOR P = XMLoadFloat3(&Keyframes.back().Translation);
XMVECTOR Q = XMLoadFloat4(&Keyframes.back().RotationQuat);
XMVECTOR zero = XMVectorSet(0.0f, 0.0f, 0.0f, 1.0f);
XMStoreFloat4x4(&M, XMMatrixAffineTransformation(S, zero, Q, P));
}
// t is between two key frames, so interpolate.
else
{
for(UINT i = 0; i < Keyframes.size()-1; ++i)
{
if( t >= Keyframes[i].TimePos && t <= Keyframes[i+1].TimePos )
{
float lerpPercent = (t - Keyframes[i].TimePos) /
(Keyframes[i+1].TimePos - Keyframes[i].TimePos);
XMVECTOR s0 = XMLoadFloat3(&Keyframes[i].Scale);
XMVECTOR s1 = XMLoadFloat3(&Keyframes[i+1].Scale);
XMVECTOR p0 = XMLoadFloat3(&Keyframes[i].Translation);
XMVECTOR p1 = XMLoadFloat3(&Keyframes[i+1].Translation);
XMVECTOR q0 = XMLoadFloat4(&Keyframes[i].RotationQuat);
XMVECTOR q1 = XMLoadFloat4(&Keyframes[i+1].RotationQuat);
XMVECTOR S = XMVectorLerp(s0, s1, lerpPercent);
XMVECTOR P = XMVectorLerp(p0, p1, lerpPercent);
XMVECTOR Q = XMQuaternionSlerp(q0, q1, lerpPercent);
XMVECTOR zero = XMVectorSet(0.0f, 0.0f, 0.0f, 1.0f);
XMStoreFloat4x4(&M, XMMatrixAffineTransformation(S, zero, Q, P));
break;
}
}
}
下图展示了两个关键帧之间插值的结果：
插值过后，我们构造了变换的矩阵，因为在着色器中我们最终使用矩阵做变换。XMMatrixAffineTransformation函数定义如下：
XMMATRIX XMMatrixAffineTransformation(
XMVECTOR Scaling,
XMVECTOR RotationOrigin,
XMVECTOR RotationQuaternion,
XMVECTOR Translation);
现在我们简单的动画系统已经完成，下一步是定义一些关键帧：
// Member data
float mAnimTimePos = 0.0f;
BoneAnimation mSkullAnimation;
//
// In constructor, define the animation keyframes
//
void QuatApp::DefineSkullAnimation()
{
//
// Define the animation keyframes
//
XMVECTOR q0 = XMQuaternionRotationAxis(XMVectorSet(0.0f, 1.0f, 0.0f, 0.0f),
XMConvertToRadians(30.0f));
XMVECTOR q1 = XMQuaternionRotationAxis(XMVectorSet(1.0f, 1.0f, 2.0f, 0.0f),
XMConvertToRadians(45.0f));
XMVECTOR q2 = XMQuaternionRotationAxis(XMVectorSet(0.0f, 1.0f, 0.0f, 0.0f),
XMConvertToRadians(-30.0f));
XMVECTOR q3 = XMQuaternionRotationAxis(XMVectorSet(1.0f, 0.0f, 0.0f, 0.0f),
XMConvertToRadians(70.0f));
mSkullAnimation.Keyframes.resize(5);
mSkullAnimation.Keyframes[0].TimePos = 0.0f;
mSkullAnimation.Keyframes[0].Translation = XMFLOAT3(-0.0f, 0.0f);
mSkullAnimation.Keyframes[0].Scale = XMFLOAT3(0.25f, 0.25f, 0.25f);
XMStoreFloat4(&mSkullAnimation.Keyframes[0].RotationQuat, q0);
mSkullAnimation.Keyframes[1].TimePos = 2.0f;
mSkullAnimation.Keyframes[1].Translation = XMFLOAT3(0.0f, 2.0f, 10.0f);
mSkullAnimation.Keyframes[1].Scale = XMFLOAT3(0.5f, 0.5f, 0.5f);
XMStoreFloat4(&mSkullAnimation.Keyframes[1].RotationQuat, q1);
mSkullAnimation.Keyframes[2].TimePos = 4.0f;
mSkullAnimation.Keyframes[2].Translation = XMFLOAT3(7.0f, 0.0f, 0.0f);
mSkullAnimation.Keyframes[2].Scale = XMFLOAT3(0.25f, 0.25f, 0.25f);
XMStoreFloat4(&mSkullAnimation.Keyframes[2].RotationQuat, q2);
mSkullAnimation.Keyframes[3].TimePos = 6.0f;
mSkullAnimation.Keyframes[3].Translation = XMFLOAT3(0.0f, 1.0f, -10.0f);
mSkullAnimation.Keyframes[3].Scale = XMFLOAT3(0.5f, 0.5f, 0.5f);
XMStoreFloat4(&mSkullAnimation.Keyframes[3].RotationQuat, q3);
mSkullAnimation.Keyframes[4].TimePos = 8.0f;
mSkullAnimation.Keyframes[4].Translation = XMFLOAT3(-0.0f, 0.0f);
mSkullAnimation.Keyframes[4].Scale = XMFLOAT3(0.25f, 0.25f, 0.25f);
XMStoreFloat4(&mSkullAnimation.Keyframes[4].RotationQuat, q0);
}
最后一步使根据时间进行插值操作：
void QuatApp::UpdateScene(float dt)
{
…
// Increase the time position.
mAnimTimePos += dt;
if(mAnimTimePos >= mSkullAnimation.GetEndTime())
{
// Loop animation back to beginning.
mAnimTimePos = 0.0f;
}
// Get the skull’s world matrix at this time instant.
mSkullAnimation.Interpolate(mAnimTimePos, mSkullWorld);
…
}
现在骷髅的世界矩阵每一帧都根据动画来更新。
7 总结
一个有序的4个实时q = (x, y, z, w) = (q1, q2, q3, q4)是一个四元数，一般都简写成q = (u, w) = (x, y, z, w)，并且我们将u = (x, y, z)称为虚向量部分，w为实数部分，进一步它的加减乘除定义为：
四元数乘法不满足交换律，但是满足结合律，四元数e = (0, 0, 0, 1)用以恒等式。四元数支持乘法分配律p(q + r) = pq + pr和(q + r)p = qp + rp；
我们可以将任意实数写成四元数s = (0, 0, 0, s)，也可以将任意向量转换成四元数u = (u, 0)。实数部分为0的四元数为纯四元数。四元数可以和标量相乘：s(p1, p2, p3, p4) = (sp1, sp2, sp3, sp4) = (p1, p2, p3, p4)s，特殊的地方在于标量和四元数的乘法支持交换律；
共轭四元数和四元数范式的定义；
逆四元数的定义和计算；
单位四元数可以写成极向量表达q = (u, q4)，其中n是单位向量；
如果q是一个单位四元数，那么q = (sinθn,cosθ) for ||n|| = 1 and θ ∈ [0, π]，旋转运算为Rq(v)=qvq−1=qvq∗R_q(v) = qvq^{-1} = qvq^*Rq(v)=qvq−1=qvq∗表示将点/向量围绕n旋转2θ。Rq有矩阵表达，任何旋转矩阵都可以转换成四元数用来表达旋转；
我们可以使用球面插值来对两个用单位四元数表示的方向进行插值。