绕任意轴旋转矩阵推导

该文是在学习 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)可得下式：
   v c = ( ∥ v ∥ cos ⁡ α ) a ^ = ( ∥ v ∥ ∥ a ^ ∥ cos ⁡ α ) a ^ = ( v ⋅ a ^ ) a ^

   vc=(‖v‖cos⁡α)a^=(‖v‖‖a^‖cos⁡α)a^=(v⋅a^)a^" role="presentation" style="position: relative;">vc=(∥v∥cosα)a^=(∥v∥∥a^∥cosα)a^=(v⋅a^)a^$$\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}$$
   vc​​=(∥v∥cosα)a^=(∥v∥∥a^∥cosα)a^=(v⋅a^)a^​
   故：
   $${\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$
   现在来求解：
   v 1 ′ = p r o j e c t ( v 1 ′ , v 1 ) + p r o j e c t ( v 1 ′ , v 2 ) = ∥ v 1 ′ ∥ cos ⁡ θ v 1 ∥ v 1 ∥ + ∥ v 1 ′ ∥ sin ⁡ θ v 2 ∥ v 2 ∥ = v 1 cos ⁡ θ + v 2 sin ⁡ θ

   v1′=project(v1′,v1)+project(v1′,v2)=‖v1′‖cos⁡θv1‖v1‖+‖v1′‖sin⁡θv2‖v2‖=v1cos⁡θ+v2sin⁡θ" role="presentation" style="position: relative;">v′1=project(v′1,v1)+project(v′1,v2)=∥v′1∥cosθv1∥v1∥+∥v′1∥sinθv2∥v2∥=v1cosθ+v2sinθ$$\begin{array}{rl}{\mathbf{v}}_{\mathbf{1}}^{\prime }& =project({\mathbf{v}}_{\mathbf{1}}^{\prime },{\mathbf{v}}_{\mathbf{1}})+project({\mathbf{v}}_{\mathbf{1}}^{\prime },{\mathbf{v}}_{\mathbf{2}})\\ & =\Vert {\mathbf{v}}_{\mathbf{1}}^{\prime }\Vert \cos \theta \frac{{\mathbf{v}}_{\mathbf{1}}}{\Vert {\mathbf{v}}_{\mathbf{1}}\Vert }+\Vert {\mathbf{v}}_{\mathbf{1}}^{\prime }\Vert \sin \theta \frac{{\mathbf{v}}_{\mathbf{2}}}{\Vert {\mathbf{v}}_{\mathbf{2}}\Vert }\\ & ={\mathbf{v}}_{\mathbf{1}}\cos \theta +{\mathbf{v}}_{\mathbf{2}}\sin \theta \end{array}$$
   v1′​​=project(v1′​,v1​)+project(v1′​,v2​)=∥v1′​∥cosθ∥v1​∥v1​​+∥v1′​∥sinθ∥v2​∥v2​​=v1​cosθ+v2​sinθ​
   上式中相当于把 ${\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 }$ 的式子：
   v ′ = v c + v 1 ′ = v c + v 1 cos ⁡ θ + v 2 sin ⁡ θ

   v′=vc+v1′=vc+v1cos⁡θ+v2sin⁡θ" role="presentation" style="position: relative;">v′=vc+v′1=vc+v1cosθ+v2sinθ$$\begin{array}{rl}{\mathbf{v}}^{\prime }& ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_{\mathbf{1}}^{\prime }\\ & ={\mathbf{v}}_{\mathbf{c}}+{\mathbf{v}}_{\mathbf{1}}\cos \theta +{\mathbf{v}}_{\mathbf{2}}\sin \theta \end{array}$$
   v′​=vc​+v1′​=vc​+v1​cosθ+v2​sinθ​

---

式子推导完成后，我们需要获得旋转矩阵的一般表达。
假设 $\mathbf{v}=(1,0,0)$ ，$\hat{\mathbf{a}}=(a_x,a_y,a_z)$
v c = ( v ⋅ a ^ ) a ^ = ( a x 2 , a x a y , a x a z )

vc​​=(v⋅a^)a^=(ax2​,ax​ay​,ax​az​)​
v 1 = v − ( v ⋅ a ^ ) a ^ = ( 1 − a x 2 , − a x a y , − a x a z )

v1=v−(v⋅a^)a^=(1−ax2,−axay,−axaz)" role="presentation" style="position: relative;">v1=v−(v⋅a^)a^=(1−a2x,−axay,−axaz)$$\begin{array}{rl}{\mathbf{v}}_{\mathbf{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}$$

v1​​=v−(v⋅a^)a^=(1−ax2​,−ax​ay​,−ax​az​)​
v 2 = v 1 × 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 )

v2=v1×a^=(−axayaz+axayaz,(1−ax2)az+ax2az,−(1−ax2)ay+ax2ay)=(0,az,−ay)" role="presentation" style="position: relative;">v2=v1×a^=(−axayaz+axayaz,(1−a2x)az+a2xaz,−(1−a2x)ay+a2xay)=(0,az,−ay)$$\begin{array}{rl}{\mathbf{v}}_{\mathbf{2}}& ={\mathbf{v}}_{\mathbf{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}$$

v2​​=v1​×a^=(−ax​ay​az​+ax​ay​az​,(1−ax2​)az​+ax2​az​,−(1−ax2​)ay​+ax2​ay​)=(0,az​,−ay​)​

v ′ = v c + v 1 ′ = v c + v 1 cos ⁡ θ + v 2 sin ⁡ θ = ( a x 2 , a x a y , a x a z ) + ( 1 − a x 2 , − a x a y , − a x a z ) cos ⁡ θ + ( 0 , a z , − a y ) sin ⁡ θ = ( a x 2 + cos ⁡ θ − a x 2 cos ⁡ θ , a x a y − a x a y cos ⁡ θ + a z sin ⁡ θ , a x a z − a x a z cos ⁡ θ − a y sin ⁡ θ ) = [ a x 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ , a x a y ( 1 − cos ⁡ θ ) + a z sin ⁡ θ , a x a z ( 1 − cos ⁡ θ ) − a y sin ⁡ θ ]

v′​=vc​+v1′​=vc​+v1​cosθ+v2​sinθ=(ax2​,ax​ay​,ax​az​)+(1−ax2​,−ax​ay​,−ax​az​)cosθ+(0,az​,−ay​)sinθ=(ax2​+cosθ−ax2​cosθ,ax​ay​−ax​ay​cosθ+az​sinθ,ax​az​−ax​az​cosθ−ay​sinθ)=[ax2​(1−cosθ)+cosθ,ax​ay​(1−cosθ)+az​sinθ,ax​az​(1−cosθ)−ay​sinθ]​

设置旋转矩阵为 $\mathbf{R}$ ，则 $\mathbf{Rv}={\mathbf{v}}^{\prime }$
即
[ x a y a z a x b y b z b x c y c z c ] [ 1 0 0 ] = [ a x 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ a x a y ( 1 − cos ⁡ θ ) + a z sin ⁡ θ a x a z ( 1 − cos ⁡ θ ) − a y sin ⁡ θ ]

= ⎣ ⎡​xa​xb​xc​​ya​yb​yc​​za​zb​zc​​⎦ ⎤​⎣ ⎡​100​⎦ ⎤​=⎣ ⎡​ax2​(1−cosθ)+cosθax​ay​(1−cosθ)+az​sinθax​az​(1−cosθ)−ay​sinθ​⎦ ⎤​
通过上式可以得到旋转矩阵的第一列向量，即
R = [ x a y a z a x b y b z b x c y c z c ] = [ a x 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ y a z a a x a y ( 1 − cos ⁡ θ ) + a z sin ⁡ θ y b z b a x a z ( 1 − cos ⁡ θ ) − a y sin ⁡ θ y c z c ] \mathbf{R}=

[xayazaxbybzbxcyczc]" role="presentation" style="position: relative;">⎡⎣⎢xaxbxcyaybyczazbzc⎤⎦⎥$$\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]$$

=

[ax2(1−cos⁡θ)+cos⁡θyazaaxay(1−cos⁡θ)+azsin⁡θybzbaxaz(1−cos⁡θ)−aysin⁡θyczc]" role="presentation" style="position: relative;">⎡⎣⎢a2x(1−cosθ)+cosθaxay(1−cosθ)+azsinθaxaz(1−cosθ)−aysinθyaybyczazbzc⎤⎦⎥$$\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]$$

R=⎣ ⎡​xa​xb​xc​​ya​yb​yc​​za​zb​zc​​⎦ ⎤​=⎣ ⎡​ax2​(1−cosθ)+cosθax​ay​(1−cosθ)+az​sinθax​az​(1−cosθ)−ay​sinθ​ya​yb​yc​​za​zb​zc​​⎦ ⎤​

同上，我们可以设 $\mathbf{v}=(0,1,0),(0,0,1)$ 来得到另外两个列向量的表达式。
最终的矩阵形式如下：
R = [ a x 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ a x a y ( 1 − cos ⁡ θ ) − a z sin ⁡ θ a x a z ( 1 − cos ⁡ θ ) + a y sin ⁡ θ a x a y ( 1 − cos ⁡ θ ) + a z sin ⁡ θ a y 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ a y a z ( 1 − cos ⁡ θ ) − a x sin ⁡ θ a x a z ( 1 − cos ⁡ θ ) − a y sin ⁡ θ a y a z ( 1 − cos ⁡ θ ) + a x sin ⁡ θ a z 2 ( 1 − cos ⁡ θ ) + cos ⁡ θ ] \mathbf{R} =

R=⎣ ⎡​ax2​(1−cosθ)+cosθax​ay​(1−cosθ)+az​sinθax​az​(1−cosθ)−ay​sinθ​ax​ay​(1−cosθ)−az​sinθay2​(1−cosθ)+cosθay​az​(1−cosθ)+ax​sinθ​ax​az​(1−cosθ)+ay​sinθay​az​(1−cosθ)−ax​sinθaz2​(1−cosθ)+cosθ​⎦ ⎤​

---

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));
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23