绕任意轴旋转矩阵推导

该文是在学习 Physically Based Rendering 第2.7.6节绕任意轴旋转时对其公式的推导产生了兴趣。

---

首先，如图所示：

已知条件：
1). $\mathbf{v}$ 是被旋转的向量。
2). $\mathbf{a}$ 是围绕旋转的轴。
3). $\theta$ 是旋转的角度。

解决思路：
通过构建坐标系 $(p,{\mathbf{v}}_1,{\mathbf{v}}_2)$ 获得绕该坐标系旋转的公式并应用到 $(1,0,0),(0,1,0),(0,0,1)$ 的坐标系上即可获得最终应用的旋转矩阵。

假设条件：
1). $\mathbf{a}$ 是单位向量，故之后用 $\hat{\mathbf{a}}$ 表示，即 $\Vert \hat{\mathbf{a}}\Vert =1$。
2). $\alpha$ 为 $\mathbf{v}$ 与 $\hat{\mathbf{a}}$ 的夹角。

求解：

1. 构建向量 ${\mathbf{v}}_1$ ，注意全程需要结合图来理解。
   $${\mathbf{v}}_1=\mathbf{v}-{\mathbf{v}}_{\mathbf{c}}$$
   其中 ${\mathbf{v}}_{\mathbf{c}}=project(\mathbf{v},\hat{\mathbf{a}})$ ，即向量 $\mathbf{v}$ 在向量 $\hat{\mathbf{a}}$ 上的投影，根据假设条件(1)和(2)可得下式：
   $$\begin{array}{rl}{\mathbf{v}}_{\mathbf{c}}& =(\Vert \mathbf{v}\Vert \cos \alpha )\hat{\mathbf{a}}\\ & =(\Vert \mathbf{v}\Vert \Vert \hat{\mathbf{a}}\Vert \cos \alpha )\hat{\mathbf{a}}\\ & =(\mathbf{v}\cdot \hat{\mathbf{a}})\hat{\mathbf{a}}\end{array}$$
   故：
   $${\mathbf{v}}_1=\mathbf{v}-(\mathbf{v}\cdot \hat{\mathbf{a}})\hat{\mathbf{a}}$$

2. 由于 ${\mathbf{v}}_{\mathbf{c}}$ 是投影，故得到的 ${\mathbf{v}}_1\perp \hat{\mathbf{a}}$ ，因此相当于有了两个 basic vector，通过叉乘 我们可以得到第三个基向量 ${\mathbf{v}}_2$ ，注意这里是左手坐标系，因此因该用 ${\mathbf{v}}_1\times \hat{\mathbf{a}}$ 而不是 $\hat{\mathbf{a}}\times {\mathbf{v}}_1$，之前就是犯了这个错误，导致最后推出来的式子不对，区别可以看这篇文章：让人懵圈的左右手坐标系及Unity中的叉积
   $${\mathbf{v}}_2={\mathbf{v}}_1\times \hat{\mathbf{a}}$$
   这里可以通过该式获得一些隐含条件，而原书正是忽略了这些细节导致最后的公式得到的比较突然。
   $$\begin{gathered} \Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\times \hat{\mathbf{a}}\Vert \\ \Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert \Vert \hat{\mathbf{a}}\Vert \sin \beta \end{gathered}$$
   其中如上所述， ${\mathbf{v}}_1\perp \hat{\mathbf{a}}$ ，故 $\beta =\frac{\pi }{2}$ ，那么 $\Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert \Vert \hat{\mathbf{a}}\Vert$ ，又 $\Vert \hat{\mathbf{a}}\Vert =1$ ，则 $\Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert$

3. 现在有了三个基向量，如图，其中我们最终的目标如下：
   $${\mathbf{v}}^{\prime }={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_1^{\prime }$$
   当前需要求解的是 ${\mathbf{v}}_1^{\prime }$
   这里也有一个隐含条件，因为向量 $\mathbf{v}$ 绕 $\hat{\mathbf{a}}$ 轴旋转，且 $\Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert$，因此不管怎么角度怎么旋转，都在一个圆中，即 $\Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert =\Vert {\mathbf{v}}_1^{\prime }\Vert$
   现在来求解：
   $$\begin{array}{rl}{\mathbf{v}}_1^{\prime }& =project({\mathbf{v}}_1^{\prime },{\mathbf{v}}_1)+project({\mathbf{v}}_1^{\prime },{\mathbf{v}}_2)\\ & =\Vert {\mathbf{v}}_1^{\prime }\Vert \cos \theta \frac{{\mathbf{v}}_1}{\Vert {\mathbf{v}}_1\Vert }+\Vert {\mathbf{v}}_1^{\prime }\Vert \sin \theta \frac{{\mathbf{v}}_2}{\Vert {\mathbf{v}}_2\Vert }\\ & ={\mathbf{v}}_1\cos \theta +{\mathbf{v}}_2\sin \theta \end{array}$$
   上式中相当于把 ${\mathbf{v}}_1^{\prime }$ 投影到两个基向量 ${\mathbf{v}}_1,{\mathbf{v}}_2$ 上，他俩相加等于 ${\mathbf{v}}_1^{\prime }$ ，投影的时候由于 ${\mathbf{v}}_1,{\mathbf{v}}_2$ 不是单位向量，故需要归一化，又因为 $\Vert {\mathbf{v}}_2\Vert =\Vert {\mathbf{v}}_1\Vert =\Vert {\mathbf{v}}_1^{\prime }\Vert$ ，故消去 $\Vert {\mathbf{v}}_2\Vert ,\Vert {\mathbf{v}}_1\Vert ,\Vert {\mathbf{v}}_1^{\prime }\Vert$ 得到 ${\mathbf{v}}_1^{\prime }$

4. 有了前面的铺垫，最终得到 ${\mathbf{v}}^{\prime }$ 的式子：
   $$\begin{array}{rl}{\mathbf{v}}^{\prime }& ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_1^{\prime }\\ & ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_1\cos \theta +{\mathbf{v}}_2\sin \theta \end{array}$$

---

式子推导完成后，我们需要获得旋转矩阵的一般表达。
假设 $\mathbf{v}=(1,0,0)$ ，$\hat{\mathbf{a}}=(a_x,a_y,a_z)$
$$\begin{array}{rl}{\mathbf{v}}_{\mathbf{c}}& =(\mathbf{v}\cdot \hat{\mathbf{a}})\hat{\mathbf{a}}\\ & =(a_x^2,a_x{a}_y,a_x{a}_z)\end{array}$$
$$\begin{array}{rl}{\mathbf{v}}_1& =\mathbf{v}-(\mathbf{v}\cdot \hat{\mathbf{a}})\hat{\mathbf{a}}\\ & =(1-a_x^2,-a_x{a}_y,-a_x{a}_z)\end{array}$$
$$\begin{array}{rl}{\mathbf{v}}_2& ={\mathbf{v}}_1\times \hat{\mathbf{a}}\\ & =(-a_x{a}_y{a}_z+a_x{a}_y{a}_z,(1-a_x^2)a_z+a_x^2{a}_z,-(1-a_x^2)a_y+a_x^2{a}_y)\\ & =(0,a_z,-a_y)\end{array}$$

$$\begin{array}{rl}{\mathbf{v}}^{\prime }& ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_1^{\prime }\\ & ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_1\cos \theta +{\mathbf{v}}_2\sin \theta \\ & =(a_x^2,a_x{a}_y,a_x{a}_z)+(1-a_x^2,-a_x{a}_y,-a_x{a}_z)\cos \theta +(0,a_z,-a_y)\sin \theta \\ & =(a_x^2+\cos \theta -a_x^2\cos \theta ,a_x{a}_y-a_x{a}_y\cos \theta +a_z\sin \theta ,a_x{a}_z-a_x{a}_z\cos \theta -a_y\sin \theta )\\ & =[a_x^2(1-\cos \theta )+\cos \theta ,a_x{a}_y(1-\cos \theta )+a_z\sin \theta ,a_x{a}_z(1-\cos \theta )-a_y\sin \theta ]\end{array}$$

设置旋转矩阵为 $\mathbf{R}$ ，则 $\mathbf{Rv}={\mathbf{v}}^{\prime }$
即
$$\left[\begin{array}{ccc}{x}_a& y_a& z_a\\ x_b& y_b& z_b\\ x_c& y_c& z_c\end{array}\right]\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right]=\left[\begin{array}{c}{a}_x^2(1-\cos \theta )+\cos \theta \\ a_x{a}_y(1-\cos \theta )+a_z\sin \theta \\ a_x{a}_z(1-\cos \theta )-a_y\sin \theta \end{array}\right]$$
通过上式可以得到旋转矩阵的第一列向量，即
$$\mathbf{R}=\left[\begin{array}{ccc}{x}_a& y_a& z_a\\ x_b& y_b& z_b\\ x_c& y_c& z_c\end{array}\right]=\left[\begin{array}{ccc}{a}_x^2(1-\cos \theta )+\cos \theta & y_a& z_a\\ a_x{a}_y(1-\cos \theta )+a_z\sin \theta & y_b& z_b\\ a_x{a}_z(1-\cos \theta )-a_y\sin \theta & y_c& z_c\end{array}\right]$$

同上，我们可以设 $\mathbf{v}=(0,1,0),(0,0,1)$ 来得到另外两个列向量的表达式。
最终的矩阵形式如下：
$$\mathbf{R}=\left[\begin{array}{ccc}{a}_x^2(1-\cos \theta )+\cos \theta & a_x{a}_y(1-\cos \theta )-a_z\sin \theta & a_x{a}_z(1-\cos \theta )+a_y\sin \theta \\ a_x{a}_y(1-\cos \theta )+a_z\sin \theta & a_y^2(1-\cos \theta )+\cos \theta & a_y{a}_z(1-\cos \theta )-a_x\sin \theta \\ a_x{a}_z(1-\cos \theta )-a_y\sin \theta & a_y{a}_z(1-\cos \theta )+a_x\sin \theta & a_z^2(1-\cos \theta )+\cos \theta \end{array}\right]$$

---

PBRT源码：

Transform Rotate(Float theta, const Vector3f &axis) {
    
    Vector3f a = Normalize(axis);
    Float sinTheta = std::sin(Radians(theta));
    Float cosTheta = std::cos(Radians(theta));
    Matrix4x4 m;
    // 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 and third basis vectors
    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;

    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 Transform(m, Transpose(m));
}