图形学中的矩阵变换

摘要

图形学因为要处理三维中的物体，所以经常要用到矩阵变换，包括基础的模型变换（旋转、平移、缩放变换），以及投影变换、视口变换。这其中有很多有意思的数学知识。网络上虽然有很多介绍矩阵变换的博客，但是很多只介绍了基本的矩阵运算，没有深入介绍四元数等高级用法。因此本文将从模型变换入手，介绍图形学中的插值操作，包括四元数等内容。

模型变换

模型变换是将一个点（模型）变换到另一个点（模型）的操作。一般我们会把三维点用齐次坐标（homogeneous notation）表示，所以变换矩阵是$4\times 4$的。另一个需要注意的点是，矩阵乘法运算一般情况下不满足交换律，所以变换的顺序很重要，一个直观的例子是原点处的物体先平移再绕原点旋转的结果与先绕原点旋转再平移的结果是不一样的。这里我们约定，平移变换$T$、缩放变换$S$、旋转变换$R$的顺序依次是先缩放再旋转最后平移，表示为：
$$M=TRS$$

平移变换

平移变换比较自然的一种表示形式是向量——这也是向量的物理意义。假设一个平移操作$\mathbf{t}=(t_x,t_y,t_x)$，写为矩阵是：
$$\mathbf{T}(\mathbf{t})=T(t_x,t_y,t_z)=\left(\begin{array}{cccc}1& 0& 0& t_x\\ 0& 1& 0& t_y\\ 0& 0& 1& t_x\\ 0& 0& 0& 1\end{array}\right)$$
平移矩阵的逆矩阵也很容易求，将$\mathbf{t}$取负即可：${\mathbf{T}}^{-1}(\mathbf{t})=\mathbf{T}(-\mathbf{t})$。逻辑上比较容易理解。

缩放变换

缩放变换用于将$x$、$y$、$z$轴依次按系数$(s_x,s_y,s_y)$进行缩放，缩放矩阵表示为：
$$\mathbf{S}(\mathbf{s})=\left(\begin{array}{cccc}{s}_x& 0& 0& 0\\ 0& s_y& 0& 0\\ 0& 0& s_z& 0\\ 0& 0& 0& 1\end{array}\right)$$
其逆矩阵${\mathbf{S}}^{-1}(\mathbf{s})=\mathbf{S}(1/s_x,1/s_y,1/s_z)$。

旋转变换

旋转变换相较于前两种变换要复杂一些。我们先看简单的形式，假设旋转轴是$x$、$y$、$z$轴，旋转角度为$\varphi$，那么旋转矩阵可以写为：
$$\begin{gathered} {\mathbf{R}}_x(\varphi )=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \varphi & -\sin \varphi & 0\\ 0& \sin \varphi & \cos \varphi & 0\\ 0& 0& 0& 1\end{array}\right) \\ {\mathbf{R}}_y(\varphi )=\left(\begin{array}{cccc}\cos \varphi & 0& \sin \varphi & 0\\ 0& 1& 0& 0\\ -\sin \varphi & 0& \cos \varphi & 0\\ 0& 0& 0& 1\end{array}\right) \\ {\mathbf{R}}_z(\varphi )=\left(\begin{array}{cccc}\cos \varphi & -\sin \varphi & 0& 0\\ \sin \varphi & \cos \varphi & 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right) \end{gathered}$$
其逆矩阵
$${\mathbf{R}}_i^{-1}(\varphi )={\mathbf{R}}_i(-\varphi )={\mathbf{R}}_i^T(\varphi )$$
也即：
$${\mathbf{R}}_i^T(\varphi ){\mathbf{R}}_i(\varphi )=\mathbf{I}$$

其他一些有趣的性质包括：
$$\begin{array}{rl}\text{det}(\mathbf{R}(\varphi ))& =1\\ {\mathbf{R}}_i(0)& =\mathbf{I}\\ {\mathbf{R}}_i({\varphi }_1){\mathbf{R}}_i({\varphi }_2)& ={\mathbf{R}}_i({\varphi }_1+{\varphi }_2)\\ {\mathbf{R}}_i({\varphi }_1){\mathbf{R}}_i({\varphi }_2)& ={\mathbf{R}}_i({\varphi }_2){\mathbf{R}}_i({\varphi }_1)\end{array}$$

绕任意轴的旋转

绕任意轴的旋转看起来会比较复杂，但是仔细分析一下其实并不困难。推导过程主要分两步，

1. 求旋转前后两个向量的关系；
2. 将这个关系转换为矩阵；

假设我们需要将向量$\mathbf{v}$绕着轴$\mathbf{a}$旋转角度$\theta$，如下图所示。

令${\mathbf{v}}_c$表示$\mathbf{v}$在轴$\mathbf{a}$上的投影，即：
$${\mathbf{v}}_c=\mathbf{a}||\mathbf{v}||\cos \alpha =\mathbf{a}(\mathbf{v}\cdot \mathbf{a})$$
并且令：
$$\begin{gathered} {\mathbf{v}}_1=\mathbf{v}-{\mathbf{v}}_c \\ {\mathbf{v}}_2={\mathbf{v}}_1\times \mathbf{a} \end{gathered}$$
则根据${\mathbf{v}}_1$与${\mathbf{v}}_1^{\prime }$、${\mathbf{v}}_2$与${\mathbf{v}}_2^{\prime }$的关系，可以求得：
$${\mathbf{v}}^{\prime }={\mathbf{v}}_c+{\mathbf{v}}_1\cos \theta +{\mathbf{v}}_2\sin \theta$$
这样就得到了${\mathbf{v}}^{\prime }$与$\mathbf{v}$的关系。

接下来，考虑如何将这个向量关系转换为矩阵关系。思路就是，旋转矩阵的每一行相当于对相应坐标轴$(1,0,0)$，$(0,1,0)$，$(0,0,1)$的旋转操作，我们可以分别计算三个坐标轴基底旋转之后的向量，也即旋转矩阵的对应行。
$${\mathbf{R}}_x(\varphi )=\left(\begin{array}{cccc}\cos \theta +(1-\cos \theta )a_x^2& (1-\cos \theta )a_x{a}_y-a_z\sin \theta & (1-\cos \theta )a_x{a}_z+a_y\sin \theta & 0\\ (1-\cos \theta )a_x{a}_y+a_z\sin \theta & \cos \theta +(1-\cos \theta )a_y^2& (1-\cos \theta )a_y{a}_z-a_x\sin \theta & 0\\ (1-\cos \theta )a_x{a}_z-a_y\sin \theta & (1-\cos \theta )a_y{a}_z+a_x\sin \theta & \cos \theta +(1-\cos \theta )a_z^2& 0\\ 0& 0& 0& 1\end{array}\right)$$
具体的代码实现如下。

Matrix4x4 Rotate(float theta, const Vector3f &axis) {
    Vector3f a = Normalize(axis);
    float sinTheta = std::sin(Radians(theta));
    float cosTheta = std::cos(Radians(theta));
    Matrix4x4 m = Matrix4x4::Identity();
    // Compute rotation of first basis vector
    m.m[0][0] = a.x * a.x + (1 - a.x * a.x) * cosTheta;
    m.m[0][1] = a.x * a.y * (1 - cosTheta) - a.z * sinTheta;
    m.m[0][2] = a.x * a.z * (1 - cosTheta) + a.y * sinTheta;
    m.m[0][3] = 0;
    // Compute rotations of second basis vector
    m.m[1][0] = a.x * a.y * (1 - cosTheta) + a.z * sinTheta;
    m.m[1][1] = a.y * a.y + (1 - a.y * a.y) * cosTheta;
    m.m[1][2] = a.y * a.z * (1 - cosTheta) - a.x * sinTheta;
    m.m[1][3] = 0;
    // Compute rotations of third basis vector
    m.m[2][0] = a.x * a.z * (1 - cosTheta) - a.y * sinTheta;
    m.m[2][1] = a.y * a.z * (1 - cosTheta) + a.x * sinTheta;
    m.m[2][2] = a.z * a.z + (1 - a.z * a.z) * cosTheta;
    m.m[2][3] = 0;
    return m;
}

欧拉变换

由于矩阵乘法不满足交换律，因此规定旋转的先后顺序很重要。这里我们介绍比较常用的欧拉变换（Euler Transform）。

如下图所示，假定视角方向是z轴方向，物体的放置为面对z轴负方向，上方向为y轴方向。定义沿x轴正方向的旋转为pitch，沿y轴正方向的旋转为head，z轴正方向的旋转为roll，那么这三个旋转操作可以以此对应到三个旋转矩阵${\mathbf{R}}_x(p)$，${\mathbf{R}}_y(h)$，${\mathbf{R}}_z(r)$。：

欧拉变换定义了这三个旋转的操作顺序为
$$\begin{array}{rl}\mathbf{R}(h,p,r)& ={\mathbf{R}}_z(r){\mathbf{R}}_x(p){\mathbf{R}}_y(h)\\ \mathbf{R}(h,p,r)& =\left(\begin{array}{cccc}\cos r\cos h-\sin r\sin p\sin h& -\sin r\cos p& \cos r\sin h+\sin r\sin p\cos h& 0\\ \sin r\cos h+\cos r\sin p\sin h& \cos r\cos p& \sin r\sin h-\cos r\sin p\cos h& 0\\ -\cos p\sin h& \sin p& \cos p\cos h& 0\\ 0& 0& 0& 1\end{array}\right)\end{array}$$
有了这个约定，就可以顺利的从旋转矩阵提取旋转角度，或者从旋转角度恢复旋转矩阵了。
$$\begin{array}{rl}h& =\arctan (-r_{20},r_{22})\\ p& =\arcsin (r_{21})\\ r& =\arctan (-r_{01},r_{11})\end{array}$$
其中，$\arctan (y,x)$表示$y/x$的反正切函数，其值域是$[0,2\pi )$。注意，这里与传统的值域为$(-\pi /2,\pi /2)$的反正切函数不一样。因为我们定义的旋转角度是$2\pi$的范围，传统的反正切函数是无法满足的。而且，考虑$x$与$y$同号/异号的情况下，$y/x$的符号始终为正/负，是无法区分该角度所在的象限的。因此传入两个参数求反正切函数是有必要的。Real-Time Rendering将这个改进的反正切函数记为$\mathrm{a}\mathrm{t}\mathrm{a}\mathrm{n}2(\mathrm{y},\mathrm{x})$。同时，C++的标准库也内置了这个函数。

代码实现：

Matrix4x4 Rotate(float h, float p, float r){
	Matrix4x4 rotate;
	float sinh = std::sin(h);
	float cosh = std::cos(h);
	float sinp = std::sin(p);
	float cosp = std::cos(p);
	float sinr = std::sin(r);
	float cosr = std::cos(r);
	rotate(0,0) = cosr*cosh-sinr*sinp*sinh;
	rotate(0,1) = -sinr*cosp;
	rotate(0,2) = cosr*sinh+sinr*sinp*cosh;
	rotate(0,3) = 0;
	rotate(1,0) = sinr*cosh+cosr*sinp*sinh;
	rotate(1,1) = cosr*cosp;
	rotate(1,2) = sinr*sinh-cosr*sinp*cosh;
	rotate(1,3) = 0;
	rotate(2,0) = -cosp*sinh;
	rotate(2,1) = sinp;
	rotate(2,2) = cosp*cosh;
	rotate(2,3) = 0;
	rotate(3,0) = 0;
	rotate(3,1) = 0;
	rotate(3,2) = 0;
	rotate(3,3) = 1;
	return rotate;
}
void RotateDecompose(Matrix4x4 rotate, float& h, float& p, float& r){
	h = std::atan2(-rotate(2,0), rotate(2,2));
	p = std::asin(rotate(2,1));
	r = std::atan2(-rotate(0,1), rotate(1,1));
}

矩阵分解

前面介绍的基本上是如何生成模型变换矩阵，有时候我们也需要将给定的模型变换矩阵分解为平移、缩放、旋转矩阵。简便起见，我们假设给定的变换矩阵$M$仅包含平移$T$、缩放$S$和旋转$R$三个操作。

分离平移操作

平移操作是最容易分解出来的，因为矩阵的第四列的前三项代表的就是平移项。提取出来后将对应的位置置零，得到的矩阵${\mathbf{M}}^{\prime }$就是仅包含缩放和旋转的矩阵。
$$T=(M_{03},M_{13},M_{23})$$

分离旋转操作

旋转操作的分解要相对麻烦一些。但是也有很多数值分析的方法，例如QR分解，奇异值分解（SVD），极分解（Polar Decomposition）。对于我们这个应用场景来说，我们希望获得一个正交矩阵R和一个对角线矩阵，极分解是最好的方法。

极分解的问题描述是这样的，给定一个矩阵$\mathbf{M}$，计算正交系数$\mathbf{R}$使得${\mathbf{R}}^T\mathbf{R}=\mathbf{I}$，且$||\mathbf{R}-\mathbf{M}|{|}_F^2={\sum }_{i,j}(R_{ij}-M_{ij}{)}^2$最小。

具体的推导过程在这里就不展开了，有兴趣的读者可以参考[9]的附录部分。这里直接给出近似的数值解：

• 令${\mathbf{R}}_0=\mathbf{M}$
• 计算${\mathbf{R}}_{i+1}=\frac{1}{2}(\mathbf{R}+{\mathbf{R}}^{-T})$
• 迭代，直到${\mathbf{R}}_{i+1}-{\mathbf{R}}_i\approx 0$

得到的$\mathbf{R}$就是我们要求的旋转矩阵。

需要注意的一点是，用极分解的前提是待分解的矩阵仅包含旋转操作和缩放操作。

分离缩放矩阵

有了前面两步，缩放矩阵只需要用${\mathbf{M}}^{\prime }$左乘${\mathbf{R}}^{-1}$即可：
$$\mathbf{S}={\mathbf{R}}^{-1}{\mathbf{M}}^{\prime }$$

代码实现

完整的矩阵分解代码实现（参考PBRT）：

void Decompose(const Matrix4x4 &m, Vector3f *T, Matrix4x4 *Rout, Matrix4x4 *S) {
    // Extract translation _T_ from transformation matrix
    T->x = m.m[0][3];
    T->y = m.m[1][3];
    T->z = m.m[2][3];

    // Compute new transformation matrix _M_ without translation
    Matrix4x4 M = m;
    for (int i = 0; i < 3; ++i) 
        M.m[i][3] = M.m[3][i] = 0.f;
    M.m[3][3] = 1.f;

    // Extract rotation _R_ from transformation matrix
    Float norm;
    int count = 0;
    Matrix4x4 R = M;
    do {
        // Compute next matrix _Rnext_ in series
        Matrix4x4 Rnext;
        Matrix4x4 Rit = Inverse(Transpose(R));
        for (int i = 0; i < 4; ++i)
            for (int j = 0; j < 4; ++j)
                Rnext.m[i][j] = 0.5f * (R.m[i][j] + Rit.m[i][j]);

        // Compute norm of difference between _R_ and _Rnext_
        norm = 0;
        for (int i = 0; i < 3; ++i) {
            Float n = std::abs(R.m[i][0] - Rnext.m[i][0]) +
                      std::abs(R.m[i][1] - Rnext.m[i][1]) +
                      std::abs(R.m[i][2] - Rnext.m[i][2]);
            norm = std::max(norm, n);
        }
        R = Rnext;
    } while (++count < 100 && norm > .0001);
    *Rout = R;

    // Compute scale _S_ using rotation and original matrix
    *S = Matrix4x4::Mul(Inverse(R), M);
}

插值

了解了基本的几种变换，我们再来讨论一下常用的插值操作。

插值（lerp）可以理解为，求两个状态的中间状态。给定状态A和状态B，以及介于0到1的系数$t=a/(a+b)$，求$t$时刻的状态

C = lerp(A, B, t)

比如对位置插值就是求两个点中间某个点。对平移、缩放、旋转做插值就是在平移、缩放、旋转空间求两个状态的中间状态。

线性插值就是用线性刻度（linear scale）度量$x_1$和$x_2$：

$$\begin{array}{rl}\frac{x_2-x}{x-x_1}& =\frac{b}{a}=\frac{1-t}{t}\\ x& =t{x}_2+(1-t)x_1\end{array}$$

对数插值就是用对数刻度（logarithmic scale）度量$x_1$和$x_2$：

$$\begin{array}{rl}\frac{\log x_2-\log x}{\log x-\log x_1}& =\frac{b}{a}=\frac{1-t}{t}\\ x& =x_2^t{x}_1^{1-t}\end{array}$$

这两种插值方式都是常用的插值方式，后面会用到。

对平移和缩放操作做插值比较直接，一般的线性插值就足够，如果觉得线性插值过渡不够平滑，可以采用样条曲线或者贝塞尔曲线的方式插值，这个话题超出了本文的范围，因此此处略过。这里主要讨论一下对旋转操作的插值。理想情况下，对旋转操作做插值应该具有以下三种性质。

1. 力矩最小（torque minimization）；
2. 角速度固定（constant angular velocity）；
3. 满足交换律（commutivity）；

力矩最小指的是，中间状态划过的轨迹在对应的旋转空间是最短的。这个性质可以类比对点的插值，我们知道两点之间线段最短，因此对两个点的插值轨迹就是连接这两个点的线段。如果把旋转空间投影到平面上，那么中间状态的轨迹也应该是一条线段。

角速度固定指的是如果$t$的变化速度是固定的，那么旋转的速度也应该是固定的。

满足交换律指的是，如果对角度进行连续插值，比如从$A$到$B$再到$C$，那么交换插值的顺序（先从$C$到$B$再到$A$）不影响插值结果，前提是插值系数$t$要保持一致。

注：按照上面的描述[2]，感觉说的更像是结合律。不确定为什么这里被称为交换律。个人理解交换律应该是从$A$到$B$的插值结果等于从$B$到$A$的结果。

直观来讲，力矩最小保证了旋转操作走的是最近路线，角速度固定保证了旋转速度不变，交换律能保证连续插值时的稳定性。不幸的是，目前还不存在同时满足这三个性质的对旋转操作的插值方法。

那么可能的对旋转变换的插值方法有哪些呢？

1. 欧拉角线性插值
2. 旋转矩阵线性插值
3. 四元数线性插值
4. 四元数球面插值
5. 四元数对数插值

前两种方法用到了前面介绍的旋转矩阵和欧拉变换，比较直接，下面给出了具体的实现代码。注意，这里输入的矩阵均是旋转矩阵，从完整的模型矩阵分解得到旋转矩阵可以参考Decompose函数。

Matrix4x4 EulerLerp(Matrix4x4 start, Matrix4x4 end, float t){
    float starth, startp, startr, endh, endp, endr;
    RotateDecompose(start, starth, startp, startr);
    RotateDecompose(end, endh, endp, endr);
    float h = starth + (endh - starth) * ratio;
    float p = startp + (endp - startp) * ratio;
    float r = startr + (endr - startr) * ratio;
    return Rotate(h, p, r);
}

Matrix4x4 MatLerp(Matrix4x4 start, Matrix4x4 end, float t){
    return start + (end - start) * t;
}

后三种涉及到四元数，我们先看一下四元数是什么。

四元数法

四元数（quaternions）早在1843年就被提出，最早是作为一种对复数的拓展。之后1985年被引入图形学，用来表示旋转或描述朝向。四元数可以与旋转的矩阵描述、欧拉角描述相互转换，最重要的是，用四元数做旋转插值有着其他的两种描述方式不具有的优良特性，后面会具体讲到。

四元数$\hat{\mathbf{q}}$定义为：
$$\begin{gathered} \hat{\mathbf{q}}=({\mathbf{q}}_v,q_w)=i{q}_x+j{q}_y+k{q}_z+q_w \\ {i}^2=j^2=k^2=-1,jk=-kj=i,ki=-ik=j,ij=-ji=k \end{gathered}$$
其中，$q_w$表示四元数的实部，${\mathbf{q}}_v$表示四元数的虚部，$i$，$j$，$k$均为虚数单位。

四元数的数学性质

对四元数的加减乘操作与向量和虚数的操作法则类似：
$$\begin{array}{rl}\hat{\mathbf{q}}+\hat{\mathbf{r}}& =({\mathbf{q}}_v,q_w)+({\mathbf{r}}_v,r_w)=({\mathbf{q}}_v+{\mathbf{r}}_v,q_w+r_w)\\ s\hat{\mathbf{q}}& =\hat{\mathbf{q}}s\\ & =(0,s)({\mathbf{q}}_v,q_w)=(s{\mathbf{q}}_v,s{q}_w)\\ \hat{\mathbf{q}}\hat{\mathbf{r}}& =(i{q}_x+j{q}_y+k{q}_z+q_w)(i{r}_x+j{r}_y+k{r}_z+r_w)\\ & =i(q_y{r}_z-q_z{r}_y+r_w{q}_x+q_w{r}_x)\\ & +j(q_z{r}_x-q_x{r}_z+r_w{q}_y+q_w{r}_y)\\ & +k(q_x{r}_y-q_y{r}_x+r_w{q}_z+q_w{r}_z)\\ & +q_w{r}_w-q_x{r}_x-q_y{r}_y-q_z{r}_z\\ & =({\mathbf{q}}_v\times {\mathbf{r}}_v+r_w{\mathbf{q}}_v+q_w{\mathbf{r}}_v,q_w{r}_w-{\mathbf{q}}_v\cdot {\mathbf{r}}_v)\end{array}$$
注意，四元数与标量的乘法满足交换律，但是四元数之间的乘法不满足交换律。这主要是因为向量的叉乘不满足交换律。而四元数之间的乘法是满足结合律的。

共轭操作：
$${\hat{\mathbf{q}}}^{\ast }=({\mathbf{q}}_v,q_w{)}^{\ast }=(-{\mathbf{q}}_v,q_w)$$
取模操作：
$$\begin{array}{rl}||\hat{\mathbf{q}}||=n(\hat{\mathbf{q}})& =\sqrt{\hat{\mathbf{q}}{\hat{\mathbf{q}}}^{\ast }}=\sqrt{{\hat{\mathbf{q}}}^{\ast }\hat{\mathbf{q}}}=\sqrt{{\mathbf{q}}_v\cdot {\mathbf{q}}_v+q_w^2}\\ & =\sqrt{q_x^2+q_y^2+q_z^2+q_w^2}\end{array}$$
单位四元数：
$$\hat{\mathbf{i}}=(0,1)$$
逆运算：
$${\hat{\mathbf{q}}}^{-1}=\frac{1}{n(\hat{\mathbf{q}}{)}^2}{\hat{\mathbf{q}}}^{\ast }$$

四元数与旋转操作

单位四元数，即满足$n(\hat{\mathbf{q}})=1$的四元数$\hat{\mathbf{q}}$，可以写为：
$$\hat{\mathbf{q}}=(\sin \varphi {\mathbf{u}}_q,\cos \varphi )=\sin \varphi {\mathbf{u}}_q+\cos \varphi$$
其中${\mathbf{u}}_q$满足$||{\mathbf{u}}_q|{|}^2=1$。

这样的单位四元数的物理意义是，围绕轴${\mathbf{u}}_q$旋转$2\varphi$角度的旋转操作。

单位四元数也可以利用公式$\cos \varphi +i\sin \varphi =e^{i\varphi }$改写为：
$$\hat{\mathbf{q}}=\sin \varphi {\mathbf{u}}_q+\cos \varphi =e^{\varphi {\mathbf{u}}_q}$$
这样我们可以对单位四元数做一些指数、对数操作：
$$\begin{array}{rl}\log \hat{\mathbf{q}}& =\log e^{\varphi {\mathbf{u}}_q}=\varphi {\mathbf{u}}_q\\ {\hat{\mathbf{q}}}^t& =(\sin \varphi {\mathbf{u}}_q+\cos \varphi {)}^t\\ & =e^{\varphi t{\mathbf{u}}_q}\\ & =\sin (\varphi t){\mathbf{u}}_q+\cos \varphi t\end{array}$$
单位四元数的一个性质是，它的逆等于共轭：
$${\hat{\mathbf{q}}}^{-1}={\hat{\mathbf{q}}}^{\ast }$$
这个性质对于计算逆来说方便了很多。

对于一个齐次坐标表示的点$\mathbf{p}=(p_x,p_y,p_z,p_w{)}^T$，给定单位四元数$\hat{\mathbf{q}}=\sin \varphi {\mathbf{u}}_q+\cos \varphi$，那么围绕$\mathbf{u}$轴旋转$2\varphi$角度的旋转操作可以记为：
$$\hat{\mathbf{q}}\hat{\mathbf{p}}{\hat{\mathbf{q}}}^{-1}$$
其中，$\hat{\mathbf{p}}$表示向量$\mathbf{p}$的四元数形式。

从旋转的角度来说，$\hat{\mathbf{q}}$和$-\hat{\mathbf{q}}$的物理意义是一致的，都代表同一个旋转过程。

四元数与矩阵的转换

至此，我们知道了单位四元数所代表的旋转操作的意义，也知道怎么用矩阵表示绕任意轴的旋转操作，因此可以写出单位四元数$\hat{\mathbf{q}}=(q_x,q_y,q_z,q_w)$与旋转矩阵$M^q$的相互转换公式：
$$M^q=\left(\begin{array}{cccc}1-2(q_y^2+q_z^2)& 2(q_x{q}_y-q_w{q}_z)& 2(q_x{q}_z+q_w{q}_y)& 0\\ 2(q_x{q}_y+q_w{q}_z)& 1-2(q_x^2+q_z^2)& 2(q_y{q}_z-q_w{q}_x)& 0\\ 2(q_x{q}_z-q_w{q}_y)& 2(q_y{q}_z+q_w{q}_x)& 1-2(q_x^2+q_y^2)& 0\\ 0& 0& 0& 1\end{array}\right)$$
从旋转矩阵到四元数：
$$\begin{array}{rl}\text{tr}(M^q)& =4-2(q_x^2+q_y^2+q_z^2)=4\left(1-\frac{q_x^2+q_y^2+q_z^2}{q_x^2+q_y^2+q_z^2+q_w^2}\right)\\ & =\frac{4q_w^2}{q_x^2+q_y^2+q_z^2+q_w^2}=4q_w^2\\ q_w& =\frac{1}{2}\sqrt{\text{tr}(M^q)}\\ q_x& =\frac{m_{21}-m_{12}}{4q_w}\\ q_y& =\frac{m_{02}-m_{20}}{4q_w}\\ q_z& =\frac{m_{10}-m_{01}}{4q_w}\end{array}$$
或
$$\begin{array}{rl}4q_w^2& =\text{tr}(M^q)\\ 4q_x^2& =m_{00}-m_{11}-m_{22}+m_{33}\\ 4q_y^2& =-m_{00}+m_{11}-m_{22}+m_{33}\\ 4q_z^2& =-m_{00}-m_{11}+m_{22}+m_{33}\end{array}$$

四元数的插值方法

对四元数进行插值可以有三种方法，线性插值（记为lerp，如果插值后做了归一化操作，则记为nlerp），球面插值（记为slerp），和对数插值（记为llerp）。
$$\begin{array}{rl}\text{lerp}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)& =(1-t)\hat{\mathbf{q}}+t\hat{\mathbf{r}}\\ \text{slerp}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)& =\frac{\sin (\varphi (1-t))}{\sin \varphi }\hat{\mathbf{q}}+\frac{\sin (\varphi t)}{\sin \varphi }\hat{\mathbf{r}}\\ \cos \varphi & =\hat{\mathbf{q}}\cdot \hat{\mathbf{r}}\\ \text{llerp}(\hat{\mathbf{q}},\hat{\mathbf{r}},t)& ={\hat{\mathbf{q}}}^{1-t}{\hat{\mathbf{r}}}^t\\ & =\sin ({\varphi }_q(1-t)){\mathbf{u}}_q+\cos {\varphi }_q(1-t)+\sin ({\varphi }_{r}t){\mathbf{u}}_r+\cos {\varphi }_{r}t\end{array}$$
几个需要注意的点：

1. 单位四元数进行线性插值以后会变成非单位的四元数，如果继续用它来表示旋转操作的话，需要先归一化。否则插值出来的旋转操作会带有额外的缩放操作。
2. 球面插值看起来跟对数插值有点类似，但是并不是一回事。球面插值的物理意义其实是对角度进行插值，系数$t$正比于插值结果和$\hat{\mathbf{q}}$的夹角。插值公式可以从二维的插值推广得来，具体的推导可以参考[2]。
3. 球面插值的物理意义是计算球面上的测地距离，相比于其他两种方式更“平滑”。下图[4]展示了线性插值和球面插值的区别。图a表示插值的两个端点，图b表示线性插值的四等分位置，平分的是两点连线的线段；图c表示球面插值的四等分位置，平分的是两点连成的圆弧。
   
   具体实现：

Quaternion lerp(const Quaternion &q1, const Quaternion &q2, float t) {
	float scale0 = 1.0f - t;
	float scale1 = ( dot(q1, q2) >= 0.0f ) ? t : -t;
	return scale0 * q1 + scale1 * q2;
}

Quaternion nlerp(const Quaternion &q1, const Quaternion &q2, float t) {
	float scale0 = 1.0f - t;
	float scale1 = ( dot(q1, q2) >= 0.0f ) ? t : -t;
	return normalize(scale0 * q1 + scale1 * q2);
}

Quaternion slerp(const Quaternion &q1, const Quaternion &_q2, float t) {
	Quaternion q2(_q2);

	T cosTheta = dot(q1, q2);
	if (cosTheta < 0) {
		/* Take the short way! */
		q2 = -q2;
		cosTheta = -cosTheta;
	}
	if (cosTheta > .9995f) {
		// Revert to plain linear interpolation
		return normalize(q1 * (1.0f - t) +  q2 * t);
	} else {
		float theta = math::safe_acos(math::clamp(cosTheta, (float) -1.0f, (float) 1.0f));
		float thetap = theta * t;
		Quaternion qperp = normalize(q2 - q1 * cosTheta);
		return q1 * std::cos(thetap) + qperp * std::sin(thetap);
	}
}

Quaternion llerp(const Quaternion &q1, const Quaternion &q2, float t) {
	return (q1.log() * (1.0f - t) + q2.log() * t).exp();
}

几种插值方法的比较

了解了四元数插值方法，我们再回头讨论一下为什么用四元数对旋转操作插值会更好.

下图展示了五种插值方式的实际效果。

• 左上：欧拉角线性插值
• 右上：旋转矩阵线性插值
• 左中：lerp，四元数线性插值（未归一化）
• 右中：nlerp，四元数线性插值
• 左下：slerp，四元数球面插值
• 右下：llerp，四元数对数插值
  

第一种欧拉角线性插值方法，很明显角速度不是匀速的，这是因为在这个例子中，三个欧拉角都有变化，如果将旋转矩阵拆分为欧拉角形式，再做插值，虽然在三个轴上角速度都是匀速的，但是组合成旋转矩阵它的速度就不一定是匀速的了。

第二种矩阵线性插值方法，可以看出旋转过程中有微弱的缩放效果（窗户被拉伸了）。如果对插值的中间结果做矩阵分解，可以看出始终为单位矩阵的缩放矩阵被改变了！这个副作用很明显是我们所不希望的。PBRT中举的例子是说用旋转矩阵线性插值会导致旋转的物体“变模糊”，这个模糊效果就是缩放矩阵被改变导致的。如下图所示，左侧是静止状态，中间是四元数球面插值的效果，右侧是旋转矩阵线性插值的效果。注意右侧的旋转球在极点位置是模糊的，而且整个球体的边缘由于模糊而显得更大。

第三种和第四种四元数线性插值，未归一化和归一化。最大的问题是角速度也不是固定的。而且未归一化的情形略有缩放感，窗户的高度先被压缩后被拉伸。

第五种四元数球面插值，可以说是这五种方法中最“平稳”的一种了。但是它也不是那么完美——不满足交换律。因为球面插值与两个四元数的夹角有关系，因此是无法满足交换律的。不过对于绝大部分场景而言，交换律并不是必须的。比如在这个简单的示例中其实看不出来它的缺陷。

第六种四元数对数插值，不清楚实现的对不对，总感觉镜头往反方向转似乎某处的符号写反了。找到的资料非常有限，也不确定这个方法有什么实际用途。

下表展示了三种四元数方法所具有的性质[2]。

$$\begin{array}{cccc}& \text{力矩最小}& \text{角速度固定}& \text{满足交换律}\\ \text{四元数线性插值}& \text{Yes}& \text{No}& \text{Yes}\\ \text{四元数球面插值}& \text{Yes}& \text{Yes}& \text{NO}\\ \text{四元数对数插值}& \text{NO}& \text{Yes}& \text{Yes}\end{array}$$

比较不幸的是，虽然我们前面讨论了四元数的方法要比欧拉角、旋转矩阵的方法更优秀，但它也不是完美的。没有一个四元数方法满足所有的三个特性。因此对不同的场景，要选取不同的插值方法。

对于动画系统（animation systems），比如物体或者相机的旋转操作，球面插值是最优选择。而在实时游戏当中，线性插值可以作为一个次优备选，至少线性插值的计算效率更高，虽然它并不是匀速的，但是很多时候用户并不会注意到这个细节。

未尽事宜

本文主要介绍了模型变换和四元数的基本概念，并讨论了用单位四元数做旋转插值的几种方案。篇幅所限，有一些话题没有涉及到，但是同样很重要，列举在此：

• 缩放部分：沿任意轴的缩放。
• 剪切变换（shear transform）。
• 插值部分：贝塞尔曲线、样条曲线插值方法。线性插值的缺点是连接处不够平滑，即连续但不可导，而贝塞尔曲线、样条曲线能保证连接处n阶可导。
• 对法向量的变换操作，为保证变换后法向仍然垂直于切平面，法向变换的变换矩阵需要做取逆+转置操作，细节可以参考PBRT或RTR。
• 对球面插值的优化。粗暴的球面插值需要做反三角函数运算，计算代价高。具体的优化措施可以参见[6]。

另外还有几个不太明确的点，同样列举出来，以免误导读者：

• 介绍旋转的插值操作时，提到的三种性质中的结合律，不确定是否正确理解了。找到的资料很有些，而且有些矛盾。可以参考[2][3][4]的介绍。
• 关于四元数对数插值，同样不确定理解是否正确。未找到示例代码，所以实现也可能有问题，主要是因为往反方向转太奇怪了。

欢迎指正！

Reference

[1] Linear and Logarithmic Interpolation

[2] Understanding Slerp, Then Not Using It

[3] Practical Parameterization of Rotations Using the Exponential Map

[4] Quaternions, Interpolation and Animation

[5] Quaternion calculus as a basic tool in computer graphics

[6] Slerping Clock Cycles

[7] Animating Rotation with Quaternion Curves

[8] Math Magician – Lerp, Slerp, and Nlerp

[9] Matrix Animation and Polar Decomposition

[10] Real-time Rendering, 4th edition.

[11] Physically Based Rendering, 3rd edition.