第四章
Transform 变换
变换(transform)是指以点、向量、颜色等实体作为输入,并以某种方式对其进行转换的一种操作。对于计算机图形学从业者而言,熟练掌握变换相关的知识是非常重要的。通过各种变换操作可以对物体、光源和相机进行移动、变形以及设定动画;还可以确保所有的计算都在同一个坐标系下进行,以及使用不同的方式来将物体投影到一个平面上。这里只列举了变换所能完成的部分操作,但是足以证明变换在实时图形学中的重要性,或者可以说,在任何领域图形学中的重要性。
线性变换(linear transform)是指一种仅保留向量加法和标量乘法的变换,具体来说就是:
f
(
x
)
+
f
(
y
)
=
f
(
x
+
y
)
(4.1)
\mathbf{f(x)} + \mathbf{f(y)} = \mathbf{f(x + y)} \tag{4.1}
f(x)+f(y)=f(x+y)(4.1)
k f ( x ) = f ( k x ) (4.2) k \mathbf{f(x)} = \mathbf{f(} k \mathbf{x)} \tag{4.2} kf(x)=f(kx)(4.2)
例如:现在有一个变换 f ( x ) = 5 x \mathbf{f(x)} = 5\mathbf{x} f(x)=5x,它代表将输入向量 x \mathbf{x} x 的每个分量都乘以 5。为了证明这个变换是线性的,它必须满足上述两个条件(即方程4.1和方程4.2):第一个条件是成立的,因为任意两个向量先乘以5再相加,等同于两个向量相加再乘以5;第二个标量乘法的条件也是显然成立的(方程4.2)。这个函数叫做一个缩放变换(scaling transform),因为它改变了一个物体的大小(尺寸)。旋转变换(rotation transform)是另一种线性变换,它将一个向量以原点为中心进行旋转。缩放变换以及旋转变换,以及事实上所有应用于三维向量的线性变换,都可以使用一个 3 × 3 3 \times 3 3×3 矩阵来进行表示;但是, 3 × 3 3 \times 3 3×3 的矩阵尺寸通常是不够的。
函数 f ( x ) = x + ( 7 , 3 , 2 ) \mathbf{f(x)} = \mathbf{x} + (7,3,2) f(x)=x+(7,3,2) 是一个非线性变换,其中的向量 x \mathbf{x} x 是一个三维向量。对两个不同的向量分别执行这个函数,其结果为每个向量各自加上 ( 7 , 3 , 2 ) (7,3,2) (7,3,2)。让一个向量加上另一个固定的向量,意味着完成了一次平移变换(translation),即它以相同的程度对所有输入向量进行了移动,这是一种非常有用的变换类型。还可以将各种各样的变换组合在一起,例如:可以先将一个物体缩放为原来的一半,然后再将其移动到一个不同的位置上。到目前为止,使用这些变换的方式还很简单,而且很难以这种方式( 3 × 3 3 \times 3 3×3矩阵)来将各种变换组合在一起。
可以使用仿射变换(affine transform)的形式来将线性变换和平移变换组合在一起,仿射变换通常存储在一个 4 × 4 4 \times 4 4×4的矩阵中。仿射变换是指先进行一次线性变换,然后再进行一次平移操作的变换。使用齐次符号(homogeneous notation)来表示这样的四维向量,点和方向也会使用同样的方式来进行表示(小写的粗体字母),二者之间的区别是:代表方向的向量,其 w w w分量为0;代表点的向量,其 w w w分量为1。例如:方向向量可以表示为 v = ( v x v y v z 0 ) T \mathbf{v} = \begin{array}{} (v_x & v_y & v_z & 0)^T \end{array} v=(vxvyvz0)T,点可以表示为 v = ( v x v y v z 1 ) T \mathbf{v} = \begin{array}{} (v_x & v_y & v_z & 1)^T \end{array} v=(vxvyvz1)T。
实时渲染中所使用的平移、旋转、缩放、对称和剪切矩阵都是仿射类型的。仿射变换最主要的特征就是它保证了直线的平行性(即两个平行的直线在变换之后仍然是平行的),但是其长度和角度可能会发生一些变化。一个仿射变换也可以表示为一系列独立仿射变换的组合。
本章将从最基本的仿射变换开始说起,章节4.1可以看成是这些简单变换的“参考手册”,在接下来的章节中描述更加复杂的矩阵;然后会对四元数(quaternion)进行讨论和介绍,四元数是一个非常强大的变换工具;再接下来,会讨论有关顶点混合和变形的内容,这是两种用于描述网格动画的方式,它们简单且有效;最后会介绍投影变换矩阵。表4.1总结了本章节中所涉及大部分变换的符号、函数和属性;这里还需要说明一点的是,正交矩阵的逆矩阵(inverse)就是它的转置矩阵(tanspose)。
表4.1 对本章中所讨论的大部分变换的总结。
符号 | 名称 | 特征描述 |
---|---|---|
T ( t ) \mathbf{T(t)} T(t) | 平移矩阵 | 移动一个点;仿射变换 |
R x ( ρ ) \mathbf{R}_x(\rho) Rx(ρ) | 旋转矩阵 | 绕 x x x轴旋转 ρ \rho ρ弧度;绕$y,z $轴旋转的标记也是类似的;正交矩阵 & 仿射变换 |
R \mathbf{R} R | 旋转矩阵 | 任意的旋转矩阵;正交矩阵 & 仿射变换 |
S ( s ) \mathbf{S(s)} S(s) | 缩放矩阵 | 根据传入的向量 s \mathbf{s} s,沿 x , y , z x,y,z x,y,z方向进行缩放;仿射变换 |
H i j ( s ) \mathbf{H}_{ij}(s) Hij(s) | 剪切变换 | 将分量 i i i相对于分量 j j j进行剪切,其中 i , j ∈ { x , y , z } i,j \in \{ x,y,z \} i,j∈{x,y,z};仿射变换 |
E ( h , p , r ) \mathbf{E}(h,p,r) E(h,p,r) | 欧拉变换 | 根据给定的欧拉角(yaw,pitch,roll)调整矩阵的朝向;正交矩阵 & 仿射变换 |
P o ( s ) \mathbf{P}_o(s) Po(s) | 正交投影矩阵 | 平行投影到一个平面或者是体积上;仿射变换 |
P p ( s ) \mathbf{P}_p(s) Pp(s) | 透视投影矩阵 | 透视投影到一个平面或者体积上 |
s l e r p ( p ^ , r ^ , t ) \mathsf{slerp}(\mathbf{\hat{p}, \hat{r}}, t) slerp(p^,r^,t) | 球面插值 | 根据参数 t t t,在四元数 p ^ , r ^ \mathbf{\hat{p}, \hat{r}} p^,r^中创建一个插值四元数 |
这些变换是操控几何物体的基本工具,大部分图形应用程序的接口都允许开发者对这些矩阵进行自定义。有时候还有一些库会对本章节中所讨论的矩阵操作进行封装,以便开发者进行调用;但是,深入理解这些函数调用背后的矩阵操作和变换细节是非常有价值的。了解矩阵在各种函数调用之后做了些什么仅仅是一个开始,深入理解这些矩阵本身的属性,可以帮助在图形学上走得更远。例如:在处理一个正交矩阵的时候就会知道,正交矩阵的逆矩阵实际上就是它的转置矩阵,可以通过获取其转置矩阵的方式,来快速求解正交矩阵的逆矩阵。类似这样的知识可以写出速度更快的代码。
基本变换
平移
可以使用一个平移(translation)矩阵 T \mathbf{T} T来描述从一个位置到另一个位置的变化,这个矩阵需要输入一个向量 t = ( t x , t y , t z ) \mathbf{t} = (t_x, t_y, t_z) t=(tx,ty,tz)来作为参数,从而对物体进行变换。矩阵 T \mathbf{T} T的形式如下:
T ( t ) = T ( t x , t y , t z ) = ( 1 0 0 t x 0 1 0 t y 0 0 1 t z 0 0 0 1 ) (4.3) \mathbf{T(t)} = \mathbf{T}(t_x, t_y, t_z)= \left ( \begin{array}{} 1 & 0 & 0 & t_x \\ 0 & 1 & 0 & t_y \\ 0 & 0 & 1 & t_z \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.3} T(t)=T(tx,ty,tz)= 100001000010txtytz1 (4.3)
图4.1展示了平移变换的一个具体效果。可以很容易的看出,点 p = ( p x , p y , p z , 1 ) \mathbf{p} = (p_x, p_y, p_z, 1) p=(px,py,pz,1)在经过 T ( t ) \mathbf{T(t)} T(t)变换之后,生成了一个新的顶点 p ′ = ( p x + t x , p y + t y , p z + t z , 1 ) \mathbf{p^\prime} = (p_x + t_x, p_y + t_y, p_z + t_z, 1) p′=(px+tx,py+ty,pz+tz,1),很显然这是一个平移变换。请注意,一个向量 v = ( v x , v y , v z , 0 ) \mathbf{v} = (v_x, v_y, v_z, 0) v=(vx,vy,vz,0)和矩阵 T \mathbf{T} T相乘之后,是不会受到影响的,因为一个方向向量是无法被平移的;相比之下,点和向量都会受到其他仿射变换的影响。平移矩阵的逆矩阵 T − 1 ( t ) = T ( − t ) \mathbf{T^{-1}(t)} = \mathbf{T(-t)} T−1(t)=T(−t),即向量 t \mathbf{t} t取反即可;通俗的来理解就是,平移矩阵的逆矩阵就是反向平移相同距离的矩阵。
这里需要提及一下有关矩阵表示的问题。有时候会在计算机图形学中看到另一种有效的矩阵表示方法,这种方法会将平移向量放在变换矩阵的最下面一行(即第四行),DirectX中就使用了这种表示方式。在这种表示方式中,矩阵元素的顺序会被颠倒,即应用程序会以从左往右的方式来读取平移变量;这种向量和矩阵的表示方法被称作行优先(row-major)表示法(也叫做行主序)。而本书中将会采用列优先(column-major)表示法(也叫做列主序),使用哪一种表示方法仅仅是符号上的区别。当采用行主序表示法,那么变换矩阵在内存中进行存储时,代表位移的四个分量会位于所有16个分量的最后四位,其中前三位代表了具体的位移分量,最后一位是1。
旋转
旋转(rotation)变换可以将一个向量(位置或者方向),绕着一个过原点的旋转轴旋转一定的角度。与平移变换一样,旋转变换也是一个刚体变换(rigid-body transform),也就是说,被变换点之间的距离并不会发生改变,并且保持了手性(handedness),即它不会导致物体左右两边相互交换;这两个特性对于计算机图形学中调整物体位置和朝向而言非常有用。方向矩阵(orientation matrix)是一个与相机视角和物体朝向相关的旋转矩阵,它定义了物体在空间中的朝向,即物体向上和向前的方向。
可以很容易地推导出二维空间中的旋转矩阵。假设现在有一个向量 v = ( v x , v y ) \mathbf{v} = (v_x, v_y) v=(vx,vy),可以将其参数化表示为 v = ( v x , v y ) = ( r cos θ , r sin θ ) \mathbf{v} = (v_x, v_y) = (r \cos \theta, r \sin \theta) v=(vx,vy)=(rcosθ,rsinθ)。如果将向量 v \mathbf{v} v顺时针旋转 ϕ \phi ϕ度,那么可以获得一个新向量 v = ( r cos ( θ + ϕ ) , r sin ( θ + ϕ ) ) \mathbf{v} = (r \cos (\theta + \phi), r \sin (\theta + \phi)) v=(rcos(θ+ϕ),rsin(θ+ϕ)),这个变换过程可以写成如下形式:
u = ( r cos ( θ + ϕ ) r sin ( θ + ϕ ) ) = ( r ( cos θ cos ϕ − sin θ sin ϕ ) r ( sin θ cos ϕ + cos θ sin ϕ ) ) = ( cos ϕ − sin ϕ sin ϕ cos ϕ ) ⏟ R ( ϕ ) ( r cos θ r sin θ ) ⏟ v = R ( ϕ ) v , (4.4) \begin{aligned} \mathbf{u} & =\left(\begin{array}{c}r \cos (\theta+\phi) \\ r \sin (\theta+\phi)\end{array}\right)=\left(\begin{array}{l}r(\cos \theta \cos \phi-\sin \theta \sin \phi) \\ r(\sin \theta \cos \phi+\cos \theta \sin \phi)\end{array}\right) \\ & =\underbrace{\left(\begin{array}{rr}\cos \phi & -\sin \phi \\ \sin \phi & \cos \phi\end{array}\right)}_{\mathbf{R}(\phi)} \underbrace{\left(\begin{array}{l}r \cos \theta \\ r \sin \theta\end{array}\right)}_{\mathbf{v}}=\mathbf{R}(\phi) \mathbf{v},\end{aligned}\tag{4.4} u=(rcos(θ+ϕ)rsin(θ+ϕ))=(r(cosθcosϕ−sinθsinϕ)r(sinθcosϕ+cosθsinϕ))=R(ϕ) (cosϕsinϕ−sinϕcosϕ)v (rcosθrsinθ)=R(ϕ)v,(4.4)
在方程4.4的推导过程中,使用了三角函数的二角和差公式来将 cos ( θ + ϕ ) , sin ( θ + ϕ ) \cos (\theta + \phi) ,\sin (\theta + \phi) cos(θ+ϕ),sin(θ+ϕ)进行展开。在三维空间中,经常使用的旋转矩阵是 R x ( ϕ ) \mathbf{R}_x(\phi) Rx(ϕ), R y ( ϕ ) \mathbf{R}_y(\phi) Ry(ϕ), R z ( ϕ ) \mathbf{R}_z(\phi) Rz(ϕ),它们分别代表了绕 x , y , z x,y,z x,y,z轴旋转 ϕ \phi ϕ度,它们的具体形式如下:
R x ( ϕ ) = ( 1 0 0 0 0 cos ϕ − sin ϕ 0 0 sin ϕ cos ϕ 0 0 0 0 1 ) (4.5) \mathbf{R}_x(\phi) = \left ( \begin{array}{} 1 & 0 & 0 & 0 \\ 0 & \cos \phi & -\sin \phi & 0 \\ 0 & \sin \phi & \cos \phi & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.5} Rx(ϕ)= 10000cosϕsinϕ00−sinϕcosϕ00001 (4.5)
R y ( ϕ ) = ( cos ϕ 0 sin ϕ 0 0 1 0 − sin ϕ 0 cos ϕ 0 0 0 0 1 ) (4.6) \mathbf{R}_y(\phi) = \left ( \begin{array}{} \cos \phi & 0 & \sin \phi & 0 \\ 0 & 1 & & 0 \\ -\sin \phi & 0 & \cos \phi & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.6} Ry(ϕ)= cosϕ0−sinϕ00100sinϕcosϕ00001 (4.6)
R z ( ϕ ) = ( cos ϕ − sin ϕ 0 0 sin ϕ cos ϕ 0 0 0 0 1 0 0 0 0 1 ) (4.7) \mathbf{R}_z(\phi) = \left ( \begin{array}{} \cos \phi & -\sin \phi & 0 & 0 \\ \sin \phi & \cos \phi & 0 & 0 \\ 0 & 0& 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.7} Rz(ϕ)= cosϕsinϕ00−sinϕcosϕ0000100001 (4.7)
如果将这个 4 × 4 4 \times 4 4×4矩阵的第四行和第四列删除,那么就可以获得一个 3 × 3 3 \times 3 3×3矩阵。对于一个绕坐标轴旋转 ϕ \phi ϕ度的 3 × 3 3 \times 3 3×3旋转矩阵而言,矩阵的迹(trace:矩阵主对角线上的元素之和,等于矩阵的特征值之和)是一个与坐标轴无关的常量,其计算公式为:
t r ( R ) = 1 + 2 cos ϕ (4.8) \mathsf{tr}\mathbf{(R)} = 1 + 2 \cos \phi \tag{4.8} tr(R)=1+2cosϕ(4.8)
图4.4展示了旋转矩阵的具体作用效果。除了使得物体绕 i i i 轴旋转 ϕ \phi ϕ度之外,旋转矩阵 R i ( ϕ ) \mathbf{R}_i(\phi) Ri(ϕ)的另外一个特征是,它会使得所有位于 i i i轴上的点保持不变。请注意,矩阵 R \mathbf{R} R也可以用来表示绕任意轴进行旋转的旋转矩阵,可以通过对上面提到的三个绕坐标轴旋转的旋转矩阵 R x ( ϕ ) \mathbf{R}_x(\phi) Rx(ϕ), R y ( ϕ ) \mathbf{R}_y(\phi) Ry(ϕ), R z ( ϕ ) \mathbf{R}_z(\phi) Rz(ϕ)进行组合,来获得一个绕任意轴旋转的旋转矩阵。
所有旋转矩阵的行列式值都为1,并且它们都是正交矩阵。这个特征对于任意数量旋转矩阵相乘的结果也同样成立。还可以通过 R i − 1 ( ϕ ) = R i ( − ϕ ) \mathbf{R}_i^{-1}(\phi) = \mathbf{R}_i(-\phi) Ri−1(ϕ)=Ri(−ϕ)的方式,来获得旋转矩阵的逆矩阵,也就是说旋转矩阵的逆矩阵,相当于以相反方向旋转同样的角度。
示例:绕某个点旋转
假设想让一个物体以点 p \mathbf{p} p为中心,绕 z z z轴旋转 ϕ \phi ϕ度,那么其变换过程是什么样的呢?图4.2展示了这个旋转过程。绕某个点进行旋转的意思就是,这个点本身并不会被旋转影响,可以先将物体进行平移,使得点 p \mathbf{p} p和原点重合,这个过程可以通过平移矩阵 T ( − p ) \mathbf{T(-p)} T(−p)来完成;然后可以对物体进行想要的旋转操作,即 R z ( ϕ ) \mathbf{R}_z(\phi) Rz(ϕ);最后再使用反向平移矩阵 T ( p ) \mathbf{T(p)} T(p)来将这个物体平移回原来的位置。将这个三个变换组合在一起,即可获得总的变换矩阵 X \mathbf{X} X:
X = T ( p ) R z ( ϕ ) T ( − p ) (4.9) \mathbf{X} = \mathbf{T(p)} \mathbf{R}_z(\phi) \mathbf{T(-p)} \tag{4.9} X=T(p)Rz(ϕ)T(−p)(4.9)
这里请注意这三个变换的顺序,最先应用的变换矩阵位于方程的最右边,然后接下来的若干个变换依次左乘。
缩放
缩放(scaling)矩阵 S ( s ) = S ( s x , s y , s z ) \mathbf{S(s)} = \mathbf{S} (s_x, s_y, s_z) S(s)=S(sx,sy,sz)可以分别在 x , y , z x,y,z x,y,z方向上,使用缩放因子 s x , s y , s z s_x, s_y, s_z sx,sy,sz来对物体进行缩放,也就是说缩放矩阵可以放大或者缩小物体,缩放因子 s i , i ∈ { x , y , z } s_i, i \in\{ x, y, z \} si,i∈{x,y,z}的值越大,物体在该方向上的尺寸也就会变得越大;如果将缩放矩阵的其中一个分量设置为1,那么物体在该方向上的尺寸将会保持不变。方程4.10展示了缩放矩阵 S \mathbf{S} S的具体形式:
S ( s ) = ( s x 0 0 0 0 s y 0 0 0 0 s z 0 0 0 0 1 ) (4.10) \mathbf{S(s)} = \left ( \begin{array}{} s_x & 0 & 0 & 0 \\ 0 & s_y & 0 & 0 \\ 0 & 0 & s_z & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.10} S(s)= sx0000sy0000sz00001 (4.10)
图4.4展示了缩放矩阵的具体效果。当 s x = s y = s z s_x= s_y= s_z sx=sy=sz的时候,这个缩放操作被称作均匀缩放(uniform),否则会被称作非均匀缩放(nonuniform);有时候也会使用术语各向同性(isotropic)和各向异性(anisotropic)来代指均匀和不均匀。缩放矩阵的逆矩阵可以表示为 S − 1 ( s ) = S ( 1 / s x , 1 / s y , 1 / s z ) \mathbf{S^{-1}(s)} = \mathbf{S} (1/s_x, 1/s_y, 1/s_z) S−1(s)=S(1/sx,1/sy,1/sz),即按照缩放比例进行反向缩放。
S = ( 5 0 0 0 0 5 0 0 0 0 5 0 0 0 0 1 ) , S ′ = ( 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 / 5 ) . (4.11) \mathbf{S} = \left ( \begin{array}{} 5 & 0 & 0 & 0 \\ 0 & 5 & 0 & 0 \\ 0 & 0 & 5 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) , \mathbf{S^{\prime}} = \left ( \begin{array}{} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1/5 \\ \end{array} \right) . \tag{4.11} S= 5000050000500001 ,S′= 1000010000100001/5 .(4.11)
相对于使用矩阵 S \mathbf{S} S来进行均匀缩放,在使用 S ′ \mathbf{S^{\prime}} S′的时候必须遵守齐次坐标的规则。这个缩放矩阵的效率并不高,因为它在齐次坐标过程中引入了除法;即如果矩阵最右下角(位置 ( 3 , 3 ) (3,3) (3,3))上的元素值为1,则不需要进行除法操作。当然了,如果所编写的系统根本没有针对1进行测试和优化(即当点坐标的 w w w分量为1,则跳过齐次坐标除法),那么其实也就没有额外的开销。
如果向量 s \mathbf{s} s中包含1个或者3个为负的分量,那么就获得了一个反射矩阵(reflection matrix),或者叫做镜像矩阵(mirror matrix);如果有两个为负的分量,那么这个矩阵会将物体旋转180度(中心对称)。这里需要注意是,如果一个旋转矩阵和一个反射矩阵相乘,那么生成的结果仍将是一个反射矩阵。例如下列两个矩阵相乘的结果仍然是一个反射矩阵:
( cos ( π / 2 ) sin ( π / 2 ) − sin ( π / 2 ) cos ( π / 2 ) ) ⏟ r o t a t i o n ( 1 0 0 − 1 ) ⏟ r e f l e c t i o n = ( 0 − 1 − 1 0 ) (4.12) \underbrace{ \left( \begin{array}{} \cos (\pi/2) & \sin (\pi/2) \\ -\sin (\pi/2) & \cos (\pi/2) \end{array} \right) }_{\mathsf{rotation}} \underbrace{ \left( \begin{array}{} 1 & 0 \\ 0 & -1 \\ \end{array} \right) }_{\mathsf{reflection}} = \left( \begin{array}{} 0 & -1 \\ -1 & 0 \\ \end{array} \right) \tag{4.12} rotation (cos(π/2)−sin(π/2)sin(π/2)cos(π/2))reflection (100−1)=(0−1−10)(4.12)
通常在检测到一个镜像变换的时候,都会进行一些特殊处理。例如:一个顶点为逆时针顺序定义的三角形,在经过反射矩阵变换之后,其顶点顺序将会变成顺时针;顶点顺序的改变会导致错误的光照效果和背面剔除。可以通过计算左上角 3 × 3 3 \times 3 3×3行列式的值,来判断一个给定的缩放矩阵是否为一个反射矩阵。如果缩放矩阵的行列式为负数,则说明该矩阵是一个反射矩阵。例如:方程4.12中的行列式值为 0 ⋅ 0 − ( − 1 ) ⋅ ( − 1 ) = − 1 0 \cdot 0 - (-1) \cdot (-1) = -1 0⋅0−(−1)⋅(−1)=−1。
示例:按任意方向进行缩放
缩放矩阵 S \mathbf{S} S只能沿坐标轴方向进行缩放,需要使用复合变换来实现在其他方向上的缩放操作。假设现在有一个给定的方向 f \mathbf{f} f,需要首先将这个方向分解为三个标准正交向量 f x , f y , f z \mathbf{f^x, f^y, f^z} fx,fy,fz。然后构建旋转变换矩阵 F \mathbf{F} F,这个矩阵可以用于修改坐标系的基底:
F = ( f x f y f z 0 0 0 0 1 ) (4.13) \mathbf{F} = \left ( \begin{array}{} \mathbf{f^x} & \mathbf{f^y} & \mathbf{f^z} & \mathbf{0} \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.13} F=(fx0fy0fz001)(4.13)
这个变换的核心思想就是,将由 f x , f y , f z \mathbf{f^x, f^y, f^z} fx,fy,fz给出的坐标系与标准坐标轴重合,然后再使用标准的缩放矩阵对其进行缩放,最后再将坐标系反向变换回去即可。这个过程的第一步是乘以矩阵 F \mathbf{F} F的转置矩阵(也就是它的逆矩阵),然后再乘以缩放矩阵(左乘),最后再乘以矩阵 F \mathbf{F} F来将坐标系变换回去,这个过程的总变换如下所示:
X = F S ( s ) F T (4.14) \mathbf{X} = \mathbf{F} \mathbf{S(s)} \mathbf{F^T} \tag{4.14} X=FS(s)FT(4.14)
剪切
剪切(shear)变换是另一种基本变换,它可以用来对整个场景进行扭曲,从而创造出一种迷幻的效果,或者是对单个模型的外观进行扭曲。剪切变换一共包含6个基本矩阵,它们分别是 H x y ( s ) \mathbf{H}_{xy}(s) Hxy(s), H x z ( s ) \mathbf{H}_{xz}(s) Hxz(s), H y x ( s ) \mathbf{H}_{yx}(s) Hyx(s), H y z ( s ) \mathbf{H}_{yz}(s) Hyz(s), H z x ( s ) \mathbf{H}_{zx}(s) Hzx(s), H z y ( s ) \mathbf{H}_{zy}(s) Hzy(s)。其中第一个下标用于表示哪个坐标会被剪切矩阵改变,第二个下标表示将会使用哪个坐标来进行剪切。方程4.15展示了剪切矩阵的一个例子 H x z ( s ) \mathbf{H}_{xz}(s) Hxz(s),从中可以发现,剪切矩阵的两个下标可以用来找到参数 s s s在矩阵中的位置:下标中的 x x x(其索引为0)代表了第0行,下标中的 z z z(其索引为2)代表了第2列,因此可以知道参数 s s s位于剪切矩阵的第0行,第2列,即位置 ( 0 , 2 ) (0,2) (0,2)上:
H x z ( s ) = ( 1 0 s 0 0 1 0 0 0 0 1 0 0 0 0 1 ) (4.15) \mathbf{H_{xz}} (s)= \left ( \begin{array}{} 1 & 0 & s & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.15} Hxz(s)= 10000100s0100001 (4.15)
将这个矩阵和一个顶点相乘,会生成一个新的顶点: ( p x + s p z p y p z ) T \begin{array}{} (p_x + sp_z & p_y & p_z)^T \end{array} (px+spzpypz)T,图4.3生动的展示了一个单位正方形被剪切的过程。通过向相反的方向进行剪切,可以获得 H i j ( s ) \mathbf{H}_{ij}(s) Hij(s)的逆矩阵,即 H i j − 1 ( s ) = H i j ( − s ) \mathbf{H}_{ij}^{-1}(s) = \mathbf{H}_{ij}(-s) Hij−1(s)=Hij(−s)。
还可以使用一种稍微不同的剪切矩阵:
H x y ′ ( s , t ) = ( 1 0 s 0 0 1 t 0 0 0 1 0 0 0 0 1 ) (4.16) \mathbf{H^{\prime}_{xy}} (s,t)= \left ( \begin{array}{} 1 & 0 & s & 0 \\ 0 & 1 & t & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.16} Hxy′(s,t)= 10000100st100001 (4.16)
这里的剪切变换包含两个输入参数,它们代表了这两个坐标( x , y x,y x,y)都会被第三个坐标( z z z)剪切。这个特殊的剪切矩阵可以通过上面一般的剪切矩阵组合而成,即 H x y ′ ( s , t ) = H i k ( s ) H j k ( t ) \mathbf{H^{\prime}_{xy}} (s,t) = \mathbf{H}_{ik}(s) \mathbf{H}_{jk}(t) Hxy′(s,t)=Hik(s)Hjk(t),其中的 k k k代表了 第三个坐标分量的索引。具体使用哪一种形式的剪切矩阵仅仅是一个偏好上的问题。最后需要注意的是,任何剪切矩阵的行列式值都为 ∣ H ∣ = 1 \vert \mathbf{H} \vert =1 ∣H∣=1,这意味着剪切变换是一种体积保持(volume-preserving)的变换,即变换前后,物体的体积并不会发生改变(在三维空间中是体积,在二维空间中则是面积),图4.3中也展示了这个特性。
变换的连接
由于矩阵之间的乘法不具备交换律(noncommutativity),因此矩阵在乘法式子中的顺序十分重要。变换的连接是与顺序相关的。
举一个矩阵顺序相关的例子:假如现在有两个矩阵 S , R \mathbf{S,R} S,R,其中矩阵 S ( 2 , 0.5 , 1 ) \mathbf{S}(2,0.5,1) S(2,0.5,1)是一个缩放变换,它将坐标的 x x x分量缩放为原来的2倍,将 y y y分量缩放为原来的0.5倍。 R z ( π / 6 ) \mathbf{R}_z(\pi / 6) Rz(π/6)是一个旋转矩阵,它绕 z z z轴(在右手坐标系中,这里的 z z z轴指向书页的外面)顺时针旋转了 π / 6 \pi/6 π/6的角度。这两个矩阵可以用两种不同的方式相乘,它们变换的结果是完全不同的,如图4.4所示:
图4.4:上图展示了矩阵相乘的顺序依赖性。在第一行中,先进行了旋转变换 R z ( π / 6 ) \mathbf{R}_z(\pi / 6) Rz(π/6),然后再进行了缩放变换 S ( s ) \mathbf{S(s)} S(s),其中 s = ( 2 , 0.5 , 1 ) \mathbf{s} = (2,0.5,1) s=(2,0.5,1),总变换为 S ( s ) R z ( π / 6 ) \mathbf{S(s)} \mathbf{R}_z(\pi / 6) S(s)Rz(π/6);在第二行中,两个向量交换了相乘的顺序,总变换为 R z ( π / 6 ) S ( s ) \mathbf{R}_z(\pi / 6) \mathbf{S(s)} Rz(π/6)S(s),这两个组合变换的结果完全不同。对于任意的矩阵 M , N \mathbf{M,N} M,N而言, M N ≠ N M \mathbf{MN} \ne \mathbf{NM} MN=NM,即矩阵的乘法不具备交换律。
将一系列矩阵组合成一个独立矩阵有一个好处,那就是可以获得更高的执行效率。例如:想象现在有一个包含几百万顶点的游戏场景,场景中的所有物体都需要进行缩放、旋转和平移变换。这里并不会将所有的顶点都与这三个变换矩阵挨个相乘,因为这样做的效率实在是太低了;会将这三个矩阵连接成一个独立的矩阵,然后对所有顶点都应用这个相同的组合变换矩阵。这个组合矩阵可以写作 C = T R S \mathbf{C = TRS} C=TRS,请注意这里的矩阵顺序,缩放矩阵 S \mathbf{S} S应当首先作用于顶点,因此它出现在组合矩阵的最右侧。这个组合矩阵的顺序意味着 T R S p = ( T ( R ( S p ) ) ) \mathbf{TRSp} = \mathbf{(T(R(Sp)))} TRSp=(T(R(Sp))),其中 p \mathbf{p} p是待变换的顶点。顺便说一句, T R S \mathbf{TRS} TRS顺序是场景图系统中,最为常用的变换组合顺序。
值得注意的是,虽然矩阵连接的结果是与顺序相关的,但是矩阵和矩阵之间可以进行分组计算,也就是说,矩阵乘法是具有结合律的。例如:假设 T R S p \mathbf{TRSp} TRSp的顺序进行计算,并且在计算的过程中,还想顺便计算一下刚体变换 T R \mathbf{TR} TR的值;那么可以将这两个矩阵组合在一起进行计算,即 ( T R ) ( S p ) \mathbf{(TR)(Sp)} (TR)(Sp),这样就可以先获得 T R \mathbf{TR} TR的值,然后再利用这个值去计算 T R S p \mathbf{TRSp} TRSp。
刚体变换
现在想象这样的一个过程,当一个人抓起一个固体物体(比如说从桌子上拿起一支笔),然后将其移动到另一个位置上(可能是衬衫口袋里),在这个过程中,仅仅是物体的位置和朝向发生了变化,其形状和大小并没有受到任何影响。将这样只包含平移和旋转的变换叫做刚体变换(rigid-body transform),它具有保持长度、角度和手性的特点。
任何刚体变换矩阵 X \mathbf{X} X,都可以写成一个平移矩阵 T ( t ) \mathbf{T(t)} T(t)和一个旋转矩阵 R \mathbf{R} R的连接,方程4.17展示了刚体变换矩阵 X \mathbf{X} X的一般形式:
X = T ( t ) R = ( r 00 r 01 r 02 t x r 10 r 11 r 12 t y r 20 r 21 r 22 t z 0 0 0 1 ) (4.17) \mathbf{X} = \mathbf{T(t)R}= \left ( \begin{array}{} r_{00} & r_{01} & r_{02} & t_x \\ r_{10} & r_{11} & r_{12} & t_y \\ r_{20} & r_{21} & r_{22} & t_z \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \tag{4.17} X=T(t)R= r00r10r200r01r11r210r02r12r220txtytz1 (4.17)
刚体变换矩阵 X \mathbf{X} X的逆矩阵可以这样来计算:
X − 1 = ( T ( t ) R ) − 1 = R − 1 T ( t ) − 1 = R T T ( − t ) \mathbf{X^{-1} = (T(t)R)^{-1} = R^{-1}T(t)^{-1} = R^TT(-t)} X−1=(T(t)R)−1=R−1T(t)−1=RTT(−t)
为了计算这个逆矩阵,旋转矩阵 R \mathbf{R} R左上角 3 × 3 3 \times 3 3×3的子矩阵会被转置(旋转矩阵的逆矩阵),平移矩阵 T \mathbf{T} T的参数符号会被取反;然后再将这两个新矩阵调换顺序相乘,便可以获得所需的逆矩阵。计算刚体变换矩阵 X \mathbf{X} X的逆矩阵,另一种方法是将矩阵 R , T \mathbf{R,T} R,T按下列方式进行考虑(符号表示详见方程1.2):
R ˉ = ( r , 0 r , 1 r , 2 ) = ( r 0 , T r 1 , T r 2 , T ) , X = ( R ˉ t 0 T 1 ) , (4.18) \begin{array}{} \mathbf{\bar{R}} = \left ( \begin{array}{} \mathbf{r}_{,0} & \mathbf{r}_{,1} & \mathbf{r}_{,2}\end{array} \right) = \left ( \begin{array}{} \mathbf{r}_{0,} ^T \\[1mm] \mathbf{r}_{1,} ^T \\[1mm] \mathbf{r}_{2,} ^T \end{array} \right) , \\\mathbf{X} = \left ( \begin{array}{} \mathbf{\bar{R}} & \mathbf{t} \\ \mathbf{0} ^T & 1 \end{array} \right) , \end{array} \tag{4.18} Rˉ=(r,0r,1r,2)= r0,Tr1,Tr2,T ,X=(Rˉ0Tt1),(4.18)
其中 r , 0 \mathbf{r}_{,0} r,0代表了旋转矩阵中的第一列(即第一个逗号取值可以是0-2,而第二个下标的值始终为0), r 0 , T \mathbf{r}_{0,} ^T r0,T代表了旋转矩阵中的第一行。还需要注意的是,方程中的 0 \mathbf{0} 0代表一个 3 × 1 3 \times 1 3×1的零向量。在这种表示方式下,方程4.19给出了矩阵 X \mathbf{X} X的逆矩阵的计算过程:
X − 1 = ( r 0 , r 1 , r 2 , − R ˉ T t 0 0 0 1 ) (4.19) \mathbf{X^{-1}} = \left ( \begin{array}{} \mathbf{r}_{0,} & \mathbf{r}_{1,} & \mathbf{r}_{2,} & \mathbf{-\bar{R}^T{}t}\\ 0&0&0&1 \end{array} \right) \tag{4.19} X−1=(r0,0r1,0r2,0−RˉTt1)(4.19)
示例:调整相机的朝向
在图形程序中有一个十分常见的操作,即调整相机的朝向,使其观察某一个指定的点。有一个叫做 g l u L o o k A t ( ) \mathsf{gluLookAt}() gluLookAt()的函数可以完成这个操作(来自OpenGL工具库,简称为GLU)。这里将会具体展示其执行的流程,虽然这个函数现在已经很少用到了,但是这个函数所对应的操作仍然十分常见。假设现在有一个位于点 c \mathbf{c} c的相机,想让这个相机看向位于点 l \mathbf{l} l的物体,并且此时相机指向上方的方向为 u ′ \mathbf{u^{\prime}} u′,如图4.5所示。
图4.5 计算一个变换矩阵,使得位于点 c \mathbf{c} c的相机看向点 1 \mathbf{1} 1,并且向上的向量 u ′ \mathbf{u^{\prime}} u′,需要计算 r , u , v \mathbf{r},\mathbf{u},\mathbf{v} r,u,v
这里需要计算一个包含三个向量 { r , u , v } \{ \mathbf{r,u,v} \} {r,u,v}的基底。首先计算观察向量 v = ( l − c ) / ∥ l − c ∥ \mathbf{v = (l-c) / \Vert l-c \Vert} v=(l−c)/∥l−c∥,即从相机位置指向目标位置的单位向量;然后计算指向右方的向量 r = ( v × u ′ ) / ∥ v × u ′ ∥ \mathbf{r = (v \times u^{\prime}) / \Vert v \times u^{\prime} \Vert} r=(v×u′)/∥v×u′∥;在相机调整朝向之后,向量 u ′ \mathbf{u^{\prime}} u′通常并不会精确地指向相机向上的方向,因此最终指向相机上方的向量,是之前两个向量叉乘的结果,即 u = r × v \mathbf{u = r \times v} u=r×v,这个向量同时垂直于向量 v , r \mathbf{v,r} v,r,并且由于这两个向量都是单位向量,因此计算出的向量 u \mathbf{u} u也是一个单位向量。接下来将构建相机的变换矩阵 M \mathbf{M} M,这里的核心思想是,先将相机平移到坐标原点 ( 0 , 0 , 0 ) (0,0,0) (0,0,0)的位置上;然后再对基底进行变换,使得向量 r \mathbf{r} r指向 ( 1 , 0 , 0 ) (1,0,0) (1,0,0),向量 u \mathbf{u} u指向 ( 0 , 1 , 0 ) (0,1,0) (0,1,0),向量 v \mathbf{v} v指向 ( 0 , 0 , − 1 ) (0,0,-1) (0,0,−1),这个过程可以通过如下变换完成:
M = ( r x r y r z 0 u x u y u z 0 − v x − v y − v z 0 0 0 0 1 ) ⏟ c h a n g e o f b a s i s ( 1 0 0 − t x 0 1 0 − t y 0 0 1 − t z 0 0 0 1 ) ⏟ t r a n s l a t i o n = ( r x r y r z − t ⋅ r u x u y u z − t ⋅ u − v x − v y − v z − t ⋅ v 0 0 0 1 ) (4.20) \begin{array}{} \mathbf{M} &= \underbrace{ \left ( \begin{array}{} r_x & r_y & r_z & 0 \\ u_x & u_y & u_z & 0 \\ -v_x & -v_y & -v_z & 0 \\ 0 & 0 & 0 & 1 \\ \end{array} \right) }_{\mathsf{change\ of \ basis}} \underbrace{ \left ( \begin{array}{} 1 & 0 & 0 & -t_x \\ 0 & 1 & 0 & -t_y \\ 0 & 0 & 1 & -t_z \\ 0 & 0 & 0 & 1 \\ \end{array} \right) }_{\mathsf{translation}} \\[4mm] &=\left ( \begin{array}{} r_x & r_y & r_z & -\mathbf{t \cdot r} \\ u_x & u_y & u_z & -\mathbf{t \cdot u} \\ -v_x & -v_y & -v_z & -\mathbf{t \cdot v} \\ 0 & 0 & 0 & 1 \\ \end{array} \right) \end{array} \tag{4.20} M=change of basis rxux−vx0ryuy−vy0rzuz−vz00001 translation 100001000010−tx−ty−tz1 = rxux−vx0ryuy−vy0rzuz−vz0−t⋅r−t⋅u−t⋅v1 (4.20)
请注意,这里将平移矩阵和基底变换矩阵组合在了一起,第一步进行的应当是平移变换,因此平移矩阵 − t \mathbf{-t} −t会位于总变换矩阵的右侧。有一种方法可以用来帮助记忆向量 r , u , v \mathbf{r,u,v} r,u,v在矩阵中的具体位置:想要让向量 r \mathbf{r} r和向量 ( 1 , 0 , 0 ) (1,0,0) (1,0,0)重合,所以当 ( 1 , 0 , 0 ) (1,0,0) (1,0,0)与基底变换矩阵相乘的时候,其结果矩阵中的第一行,一定是向量 r \mathbf{r} r中的元素;此外,矩阵中第二行和第三行的向量一定会垂直于 r \mathbf{r} r,即 r ⋅ x = 0 \mathbf{r \cdot x} = 0 r⋅x=0。同理,当使用同样的想法来思考向量 u , v \mathbf{u,v} u,v的时候,便可以获得上述的基底变换矩阵。
法线变换
矩阵可以用于对点、线、三角形和其他几何物体进行变换,这些矩阵同样也可以对这些线或者三角形表面的切向量(tangent vector)进行变换,然而有一个重要的几何属性并不能总是使用这些矩阵直接进行变换,即表面法线(以及顶点的光照法线)。图4.6展示了如果使用同样的矩阵同时对几何物体及其法线进行变换的结果。
图4.6 左侧是原始的几何体,它从侧面展示了一个三角形及其表面法线。在中间的图中,使用了一个缩放矩阵来将几何体在 x x x轴方向上缩放为原来的0.5倍,并使用这个缩放矩阵对法线进行了变换,变换后的法线并不是这个三角形的表面法线。右图中则展示了法线正确变换后的结果。
对法线正确的变换方法是:使用原始变换矩阵的伴随矩阵(adjoint)的转置矩阵来对其进行变换,而不是使用原始变换矩阵本身。这里需要知道的是:矩阵的伴随矩阵是始终存在的。法线在经过变换之后,其长度可能会发生变化,因此在变换后通常还需要对法线进行归一化处理。
法线变换的传统方法是,计算原始变换矩阵的逆矩阵的转置,即 ( M − 1 ) T \mathbf{(M^{-1})^T} (M−1)T,这个方法现在也是可以使用的。但是并不需要完整求出这个逆矩阵,这样会做很浪费性能;并且有时候这个逆矩阵可能并不存在,因为逆矩阵是通过伴随矩阵除以矩阵的行列式获得的,但是矩阵的行列式可能会为0。行列式为0的矩阵被称为奇异矩阵(singular matrix),其逆矩阵是不存在的。
进一步说,即使只计算一个 4 × 4 4 \times 4 4×4矩阵的伴随矩阵,也是十分耗时的,而且通常也不需要这样做。由于法线仅仅只是一个向量,因此平移变换对其是没有影响的;此外,大多数模型变换都是仿射变换,这些变换并不会修改齐次坐标的 w w w分量,即它们没有执行投影操作。在这些条件下,只需要计算左上角 3 × 3 3 \times 3 3×3子矩阵的伴随矩阵,即可完成对法线的变换。
实际上,甚至都不需要计算这个伴随矩阵。假设已经知道了这个变换矩阵是完全由平移、旋转和均匀缩放(没有被拉伸或者压缩)操作组合而成的,那么首先平移变换并不会对法线产生影响,其次均匀缩放仅仅会对法线的长度产生影响,剩下的就是一系列旋转变换了,这些旋转变换最终会组合在一起,形成一个总的旋转变换矩阵。上文中说到,可以使用原始变换矩阵的逆转置矩阵,来对法线进行变换;而旋转矩阵本身是一个正交矩阵,其逆矩阵和转置矩阵是一样的,也就是说,对一个旋转矩阵进行两次转置(或者两次求逆)操作之后,得到的矩阵就是最初的旋转矩阵本身。总而言之,在这样的条件下(进行了平移、旋转和均匀缩放操作),可以直接使用模型的变换矩阵来对法线进行变换;但是如果模型变换涉及了非均匀缩放或者投影操作,那么就无法使用这种方法对法线进行变换了。
最后需要考虑是否要对法线进行归一化处理,如果模型变换只涉及平移或者旋转的话,那么法线的长度是不会发生改变的,因此也就不需要进行归一化处理。但是如果变换中涉及了非均匀缩放的话,那么这个均匀缩放系数(如果是已知的,或者是已经被提取出来的话)也可以直接用来对法线进行归一化处理(直接按照均匀缩放系数进行反向缩放即可)。例如:如果知道了一系列的均匀变换最终会导致物体比之前放大5.2倍,并且法线应用了这个变换矩阵的话,可以直接将法线的长度除以5.2即可(即缩小为原来的5.2倍)。除此之外,还可以将这一步放在对法线变换之前,让左上角 3 × 3 3 \times 3 3×3子矩阵除以这个缩放系数,从而构建一个会生成归一化结果的法线变换矩阵。
还有一点需要注意一下,如果表面法线是从变换之后的三角形中计算出来的话(例如使用三角形的边向量进行叉乘,从而获得垂直于三角形表面的法线),那么法线变换的问题就不需要进行考虑了。切向量的本质和法线并不相同,它可以直接使用原始变换矩阵进行变换。
计算逆矩阵
有很多计算都需要使用逆矩阵,例如在不同的坐标系间来回切换的时候。根据一个变换的可用信息不同,可以使用以下三种方法来计算一个逆矩阵:
- 如果某个矩阵是一次很简单的变换,或者是一系列带参简单变换的组合的话,那么可以通过反转“参数”和矩阵次序的方式,来获得这个矩阵的逆矩阵,即进行一系列的反向变换。例如:现在有一个变换矩阵 M = T ( t ) R ( ϕ ) \mathbf{M}=\mathbf{T(t)}\mathbf{R}(\phi) M=T(t)R(ϕ),则其逆矩阵为 M − 1 = R ( − ϕ ) T ( − t ) \mathbf{M^{-1}} = \mathbf{R}(-\phi)\mathbf{T(-t)} M−1=R(−ϕ)T(−t)。这种求取逆矩阵的方式十分简单且准确,这对于渲染一个大世界而言是十分重要的。
- 如果一个矩阵是正交矩阵的话,那么其逆矩阵和转置矩阵是相等的,即 M − 1 = M T \mathbf{M^{-1}} = \mathbf{M^{T}} M−1=MT。旋转矩阵是一个正交矩阵,并且任意数量的旋转矩阵组合在一起仍然是一个旋转矩阵,因此其结果也是正交的。
- 如果这些信息都不知道的话,那么还可以使用伴随矩阵法、Cramer法则、LU分解法、高斯消元法等方法来计算一个矩阵的逆矩阵。伴随矩阵法和Cramer法则通常要更好一些,因为它们所涉及的分支操作较少;在现代的GPU架构中,最好还是避免使用分支结构。
在进行性能优化的时候,还可以对逆矩阵的计算目的进行考虑。例如:如果这个逆矩阵仅仅是用来对向量进行变换的话,那么通常只需要获得左上角 3 × 3 3 \times 3 3×3子矩阵的逆矩阵即可。
特殊的矩阵变换和操作
欧拉变换
欧拉变换可以构建一个旋转矩阵,将自身(相机)或者其他物体指向一个特定的方向,这是一种十分直观的方式,其名字来源于伟大的瑞士数学家 Leonhard Euler (1707–1783)。
图4.7 欧拉变换和改变头部角度(head)、俯仰角(pitch)以及滚转角(roll)之间的关系。相机的默认位置是看向 z z z轴负半轴,同时up方向与 y y y轴正半轴重合。
首先,需要有一个默认的观察方向,通常来说都会让这个方向指向 z z z轴负半轴,并且头部指向 y y y轴正半轴,如图4.7所示(其实就是在应用观察变换矩阵之后,相机的位置和朝向)。欧拉变换是三个旋转矩阵相乘的结果,如图4.7中所示的旋转。欧拉变换通常会使用 E \mathbf{E} E来进行表示,其具体形式如下:
E ( h , p , r ) = R z ( r ) R x ( p ) R y ( h ) (4.21) \mathbf{E}(h,p,r) = \mathbf{R}_z(r) \mathbf{R}_x(p)\mathbf{R}_y(h) \tag{4.21} E(h,p,r)=Rz(r)Rx(p)Ry(h)(4.21)
在欧拉变换中,矩阵的组合方式一共有24种(6个两轴旋转,6个三轴旋转;以及内旋外旋两种方式。 12 × 2 = 24 12\times 2=24 12×2=24),方程4.21所展示的组合方式,是图形学中最为常用的一种组合方式。由于矩阵 E \mathbf{E} E是由一系列旋转矩阵连接组成的,那么矩阵 E \mathbf{E} E本身自然也是一个正交矩阵,因此该矩阵的逆矩阵可以表示为 E − 1 = E T \mathbf{E^{-1}} = \mathbf{E^{T}} E−1=ET = ( R z R x R y ) T = R y T R x T R z T = (\mathbf{R}_z \mathbf{R}_x \mathbf{R}_y)^T = \mathbf{R}_y^T \mathbf{R}_x^T\mathbf{R}_z^T =(RzRxRy)T=RyTRxTRzT,当然,一般直接使用 E T \mathbf{E}^T ET会更加方便。
其中的欧拉角参数 h , p , r h,p,r h,p,r代表了每个方向(head头部,pitch俯仰角,roll滚转角)上绕轴旋转的角度。有时候会将这三个旋转角度都叫做roll,例如头部角度叫做y-roll,俯仰角叫做x-roll等。这几个角度的名称,在不同学科领域中的称呼不太一样,例如在飞行模拟中,头部角度head通常被叫做偏航角(yaw)。
这种变换的方式是非常直观的,因此不使用专业术语也可以很形象对其进行描述。例如:改变head角度会让观察者摇头;改变俯仰角会让观察者点头;改变滚转角会让观察者倾斜头部。这里使用了head,pitch,roll来描述旋转的方向,而不是使用绕 x , y , z x,y,z x,y,z轴。这里需要注意的是,欧拉变换不仅仅可以用来调整相机的方向,还可以用来调整任意物体的朝向;同时,欧拉变换不仅可以用于世界空间中,同样也可以用于局部参考坐标系中。
值得注意的是,在一些欧拉角的表示中,会让 z z z轴指向上方;这虽然会让人感到一些困惑,但是它确实只是一种符号表示上的差异,本质上都是等价的。计算机图形学中,在如何看待和表示世界这个问题上,确实存在着一些分歧,即: y y y轴向上(y-up)还是 z z z轴向上(z-up)。在包括3D打印在内的大部分工业制造领域,都将 z z z轴正方向作为世界空间中向上的方向;而在航空和航海领域中,则将 z z z轴负方向作为世界空间中向上的方向。建筑和GIS领域通常都会使用z-up,因为建筑设计和地图一般都是在二维上进行的( x , y x,y x,y);与多媒体相关的建模系统一般会使用y-up,这与在计算机图形学中的描述是相符的,即相机屏幕也是 y y y轴正方向指向上方。这两种表示方法之间的差异仅仅是旋转 9 0 ∘ 90^{\circ} 90∘而已(或者是一个对称),但是如果不清楚此时使用的是哪一种表示方法的话,可能会导致出现一些问题。在这里,如果没有额外说明的话们将会默认使用y-up。
这里还想指出一点,相机的up方向和世界空间中up方向没有什么特殊关系,可以左右倾斜自己的头部,这时我们眼睛的视野也会相应的倾斜,此时头部的up方向和世界空间中的up方向并不是重合的。再举一个例子:如果现在世界空间使用了y-up,那么相机将会以一个鸟瞰视角俯瞰整个世界;这是因为相机视角向前旋转了 9 0 ∘ 90^{\circ} 90∘,其在世界空间中的up方向变成了 ( 0 , 0 , − 1 ) (0,0,-1) (0,0,−1)。在这种观察角度下,相机的up方向将不再指向世界空间中的正 y y y轴,而是变成了世界空间中的负 z z z轴;但是对于相机而言,自身的y-up仍然是成立的。
欧拉角在小角度变换和调整观察者方向方面十分有用,但是它也有一些严重的限制,即很难将两组欧拉角组合在一起。例如:在两组欧拉角之间进行插值,并不是简单地对每个分量分别进行插值就可以完成的。事实上,两组表示形式完全不同的欧拉角,可能会给出完全相同的方向,因此对这两组欧拉角进行插值的话,中间生成的任何一组欧拉角,在理想情况下都不应当导致物体发生旋转。这也是使用其他方向表示方法(例如四元数)的原因之一,这是十分有价值的。使用欧拉角也会导致一个叫做万向节死锁(gimbal lock)的问题。
从欧拉变换中提取参数
在某些情况下,需要从一个代表欧拉变换的矩阵中,提取出各个方向上所改变的角度,即欧拉变换的参数 h , p , r h,p,r h,p,r。这个过程如下所示:
E ( h , p , r ) = ( e 00 e 01 e 02 e 10 e 11 e 12 e 20 e 21 e 22 ) = R z ( r ) R x ( p ) R y ( h ) (4.22) \mathbf{E}(h,p,r) = \left ( \begin{array}{} e_{00} & e_{01} & e_{02} \\ e_{10} & e_{11} & e_{12} \\ e_{20} & e_{21} & e_{22} \\ \end{array} \right) = \mathbf{R}_z(r) \mathbf{R}_x(p)\mathbf{R}_y(h) \tag{4.22} E(h,p,r)= e00e10e20e01e11e21e02e12e22 =Rz(r)Rx(p)Ry(h)(4.22)
齐次变换矩阵是 4 × 4 4 \times 4 4×4的,这里只使用了左上角 3 × 3 3 \times 3 3×3子矩阵,因为这已经能够提供旋转矩阵的所有必要信息了;也就是说,在完整的 4 × 4 4 \times 4 4×4欧拉变换矩阵中,除了最右下角的元素是1之外,其他剩余的元素均为0。
将三个旋转矩阵相乘,可以获得以下结果:
E = ( cos r cos h − sin r sin p sin h − sin r cos p cos r sin h + sin r sin p cos h sin r cos h + cos r sin p sin h cos r cos p sin r sin h − cos r sin p cos h − cos p sin h sin p cos p cos h ) (4.23) \mathbf{E}= \left ( \begin{array}{} \cos r \cos h - \sin r \sin p \sin h & - \sin r \cos p & \cos r \sin h + \sin r \sin p \cos h \\ \sin r \cos h + \cos r \sin p \sin h & \cos r \cos p & \sin r \sin h - \cos r \sin p \cos h \\ -\cos p \sin h & \sin p & \cos p \cos h \\ \end{array} \right) \tag{4.23} E= cosrcosh−sinrsinpsinhsinrcosh+cosrsinpsinh−cospsinh−sinrcospcosrcospsinpcosrsinh+sinrsinpcoshsinrsinh−cosrsinpcoshcospcosh (4.23)
在方程4.23中,可以很明显的看出 sin p = e 21 \sin p = e_{21} sinp=e21;此外,可以令 e 01 e_{01} e01除以 e 11 e_{11} e11来计算 r r r,令 e 20 e_{20} e20除以 e 22 e_{22} e22来计算 h h h。具体的参数提取方程如下:
e 01 e 11 = − sin r cos r = − tan r , e 20 e 22 = − sin h cos h = − tan h (4.24) \frac{e_{01}}{e_{11}} = \frac{-\sin r}{\cos r} = -\tan r ,\qquad \frac{e_{20}}{e_{22}} = \frac{-\sin h}{\cos h} = -\tan h \tag{4.24} e11e01=cosr−sinr=−tanr,e22e20=cosh−sinh=−tanh(4.24)
也就是说,可以使用 a t a n 2 ( y , x ) \mathsf{atan2(y,x)} atan2(y,x)(包含两个参数的反正切函数)来从矩阵 E \mathbf{E} E中提取欧拉角的参数 h h h(head), p p p(pitch), r r r(roll),即:
h = a t a n 2 ( − e 20 , e 22 ) p = arcsin ( e 21 ) r = a t a n 2 ( − e 01 , e 11 ) (4.25) \begin{align}{} h &= \mathsf{atan2}(-e_{20},e_{22}) \\ p &= \arcsin (e_{21}) \\ r &= \mathsf{atan2}(-e_{01},e_{11}) \\ \end{align} \tag{4.25} hpr=atan2(−e20,e22)=arcsin(e21)=atan2(−e01,e11)(4.25)
但是,这里还需要处理一种特殊情况,即当 cos p = 0 \cos p = 0 cosp=0的时候,会遇到被称为万向节死锁的问题,此时旋转角度 r , h r,h r,h将会围绕着同一个旋转轴进行旋转(尽管它俩的旋转方向可能不同,这取决于旋转角度 p p p是 − π / 2 -\pi/2 −π/2还是 π / 2 \pi/2 π/2),在这种情况下,只需要计算其中任意一个角度即可。如果我们假设 h = 0 h=0 h=0[1769],那么此时矩阵 E \mathbf{E} E为:
E 1 = ( cos r cos h − sin r sin h 0 cos r sin h + sin r cos h sin r cos h + cos r sin h 0 sin r sin h − cos r cos h 0 1 0 ) = ( cos ( r + h ) 0 sin ( r + h ) sin ( r + h ) 0 cos ( r + h ) 0 1 0 ) , p = π / 2 E 2 = ( cos r cos h + sin r sin h 0 cos r sin h − sin r cos h sin r cos h − cos r sin h 0 sin r sin h + cos r cos h 0 − 1 0 ) = ( cos ( r − h ) 0 sin ( r − h ) sin ( r − h ) 0 cos ( r − h ) 0 − 1 0 ) , p = − π / 2 (4.25.1) \begin{align}{} \mathbf{E_1}& = \left ( \begin{array}{} \cos r \cos h - \sin r \sin h & 0 & \cos r \sin h + \sin r \cos h \\ \sin r \cos h + \cos r \sin h & 0 & \sin r \sin h - \cos r \cos h \\ 0 & 1 & 0 \\ \end{array} \right) \\[5mm]&= \left ( \begin{array}{} \cos (r+h) & 0 & \sin (r+h) \\ \sin (r+h) & 0 & \cos (r+h) \\ 0 & 1 & 0 \\ \end{array} \right) ,p = \pi/2 \\[10mm] \mathbf{E_2}&= \left ( \begin{array}{} \cos r \cos h + \sin r \sin h & 0 & \cos r \sin h - \sin r \cos h \\ \sin r \cos h - \cos r \sin h & 0 & \sin r \sin h + \cos r \cos h \\ 0 & -1 & 0 \\ \end{array} \right) \\[5mm]&= \left ( \begin{array}{} \cos (r-h) & 0 & \sin (r-h) \\ \sin (r-h) & 0 & \cos (r-h) \\ 0 & -1 & 0 \\ \end{array} \right) ,p = -\pi/2 \end{align} \tag{4.25.1} E1E2= cosrcosh−sinrsinhsinrcosh+cosrsinh0001cosrsinh+sinrcoshsinrsinh−cosrcosh0 = cos(r+h)sin(r+h)0001sin(r+h)cos(r+h)0 ,p=π/2= cosrcosh+sinrsinhsinrcosh−cosrsinh000−1cosrsinh−sinrcoshsinrsinh+cosrcosh0 = cos(r−h)sin(r−h)000−1sin(r−h)cos(r−h)0 ,p=−π/2(4.25.1)
当 cos p = 0 \cos p = 0 cosp=0的时候,由于方程4.23中的 e 01 , e 11 , e 20 , e 21 e_{01},e_{11},e_{20},e_{21} e01,e11,e20,e21都包含 cos p \cos p cosp项,因此都为0(如方程4.25.1所示),所以就无法使用方程4.25中来解析出 h , r h,r h,r 的值。而且可以在方程 4.25.1 4.25.1 4.25.1中看到,此时方程仅和 r + h r+h r+h或者 r − h r-h r−h有关;这里只需要假设其中任意一个旋转角度为0,然后求出另一个参数即可;例如 h = 0 h = 0 h=0,即 cos h = 1 , sin h = 0 \cos h = 1,\sin h =0 cosh=1,sinh=0,那么此时的欧拉变换矩阵可以写成:
E = ( cos r − sin r cos p sin r sin p sin r cos r cos p − cos r sin p 0 sin p cos p ) (4.26) \mathbf{E}= \left ( \begin{array}{} \cos r & -\sin r \cos p & \sin r \sin p \\ \sin r & \cos r \cos p & - \cos r \sin p \\ 0 & \sin p & \cos p \\ \end{array} \right) \tag{4.26} E= cosrsinr0−sinrcospcosrcospsinpsinrsinp−cosrsinpcosp (4.26)
可以看出,由于旋转角度 p p p并不会对矩阵第一列的值产生任何影响,因此当 cos p = 0 \cos p = 0 cosp=0的时候,可以使用 sin r / cos r = tan r = e 10 / e 00 \sin r / \cos r = \tan r = e_{10} / e_{00} sinr/cosr=tanr=e10/e00,即 r = a t a n 2 ( − e 10 , e 00 ) r = \mathsf{atan2}(-e_{10},e_{00}) r=atan2(−e10,e00)来计算出 r r r。
另外请注意,函数arcsin的定义域是 − π / 2 ≤ p ≤ π / 2 -\pi/2 \le p \le \pi/2 −π/2≤p≤π/2,这意味着如果使用了一个在这个范围之外的 p p p值,来创建一个欧拉变换矩阵 E \mathbf{E} E的话,那么将无法获取最初的 p p p值。参数 h , p , r h,p,r h,p,r的组合并不是唯一的,也就是说有很多组欧拉角的组合,都可以生成效果一样的变换。上述所提到的方法只是一个很简单的版本,它可能会导致数值不稳定(numerical instability)的问题,可以牺牲一些性能来避免这个问题。
当使用欧拉变换的时候,有时会发生一种叫做万向节死锁(gimbal lock)的现象,即在旋转的过程中失去了一个自由度。例如:假设现在按照 x , y , z x,y,z x,y,z的顺序进行变换,然后绕 y y y轴旋转了 π / 2 \pi/2 π/2的角度,即执行了第二个旋转;这个旋转操作会使得局部坐标系中的 z z z轴与原始的 x x x轴重合,最终导致绕 z z z轴旋转的操作是多余的。
在数学上,已经在方程4.25.1和方程4.26中看到了万向节死锁的现象,在方程中我们假设了 cos p = 0 \cos p = 0 cosp=0,即 p = ± π / 2 + 2 π k p = \pm \pi/2 + 2\pi k p=±π/2+2πk,其中 k k k是一个正整数。此时整个矩阵只跟一个角度有关系,如方程4.25.1所示,根据 p p p值的不同,这个角度可能是 r + h r+h r+h或者 r − h r-h r−h。也就是说,原来的参数 r , h r,h r,h代表了两个不同的自由度,而现在只有一个自由度了。
在建模系统中,通常会使用 x , y , z x,y,z x,y,z顺序的欧拉角,每个旋转都对应着一个局部坐标轴。在不同的系统中,通常都会使用不同的欧拉角顺序,例如:在动画中使用 z , x , y z,x,y z,x,y顺序,在动画和物理中也会使用 z , x , z z,x,z z,x,z顺序。所有这些分别指定三个旋转轴的方法都是有效的,刚才提到的最后一个有效顺序 z , x , z z,x,z z,x,z,在某些应用程序中会更加优越,因为只有当绕 x x x轴旋转 π \pi π(旋转半周, 18 0 ∘ 180^{\circ} 180∘)的时候才会出现万向节死锁。没有任何一个序列能够完美的避免万向节死锁问题,但是欧拉角仍然是最为常用的角度表示方法和旋转表示方法,因为动画师们在进行动画制作的时候,更加喜欢使用曲线编辑器来指定某个角度随着时间的变化,这种方法十分直观和便于理解。
示例:约束变换
想象现在拿着一个虚拟的扳手,这个扳手正抓着一个螺栓;为了让这个螺栓安装到位,需要绕着 x x x轴来旋转扳手。现在假设通过操作输入设备(例如鼠标、VR手套、太空球等)给出了一个旋转矩阵,这个旋转矩阵用于控制扳手的运动。但是如果直接将这个旋转矩阵应用到扳手上的话可能会导致一定的问题,因为扳手只应当绕着 x x x轴进行旋转。将这个输入的变换矩阵称为 P \mathbf{P} P,为了限制这个变换矩阵只围绕着 x x x轴进行旋转,需要使用本小节中描述的方法,将变换矩阵中的欧拉角度提取出来,然后创建一个新的变换矩阵 R x ( p ) \mathbf{R}_x(p) Rx(p)即可。这个新的变换矩阵 R x ( p ) \mathbf{R}_x(p) Rx(p),就是希望应用于扳手的变换矩阵,它将会使得扳手只绕着 x x x轴进行旋转(如果此时输入的 P \mathbf{P} P包含绕 x x x轴旋转的变换的话)。
矩阵分解
到目前为止,所做的工作都建立在这样的一个假设上,即已经知道了原始变换矩阵及其变换过程;但是更多的时候,实际上并不清楚这些信息。例如:只知道有一个变换矩阵和被变换的物体,现在的任务是从在这个总变换矩阵中,分解出各种各样的子变换矩阵,这个过程叫做矩阵分解(matrix decomposition)。
矩阵分解有很多用途,例如:
- 提取物体的缩放因子。
- 找到一个指定坐标系中所需要的变换(例如:某些系统和变换不允许使用任意的 4 × 4 4 \times 4 4×4矩阵)。
- 确定一个物体是否只经历了刚体变换。
- 在只有物体的变换矩阵可用的情况下,在动画的关键帧之间进行插值。
- 移除一个旋转变换矩阵中的剪切变换。
在前文中其实已经展示了两个矩阵分解的例子,例如从一个刚体变换中提取出平移矩阵和旋转矩阵;从一个正交矩阵中提取出欧拉角。
在上文中的两个例子中可以看到,从一个变换矩阵中提取出平移矩阵是很简单的,只需要找到 4 × 4 4 \times 4 4×4矩阵中的最后一列元素即可。还可以对变换矩阵的行列式进行检查,如果行列式的值是一个负数,那么就说明这个矩阵包含一个反射变换。而想要分离出旋转、缩放和剪切变换则需要更多的的努力。
绕任意轴旋转
有时候,能够让一个物体绕着任意轴旋转是一件很方便的事情,假设现在有一个归一化的旋转轴 r \mathbf{r} r,希望创建一个矩阵,能够让物体绕着这个旋转轴 r \mathbf{r} r旋转 α \alpha α度。
为了完成这个绕任意轴旋转的操作,首先需要进行一次空间变换,将旋转轴 r \mathbf{r} r与 x x x轴重合。我们可以通过一个旋转矩阵 M \mathbf{M} M来完成,然后当真正需要的旋转变换完成之后,再使用旋转矩阵 M − 1 \mathbf{M}^{-1} M−1将其旋转回最初的位置,整个变换过程如图4.8所示:
图4.8 上图展示了绕任意轴 r \mathbf{r} r进行旋转的变换流程。首先需要找到一组标准正交基 r , s , t \mathbf{r,s,t} r,s,t。然将这个基底与标准基底(坐标轴)重合,即使得旋转轴 r \mathbf{r} r与 x x x轴重合;然后进行绕 x \mathbf{x} x轴旋转 α \alpha α度的操作;最后再反向旋转回来。
为了计算这个旋转矩阵 M \mathbf{M} M,需要找到两个与旋转轴 r \mathbf{r} r正交,且彼此相互正交的轴,从而构建一组标准正交基。其实只要能找到第二个轴 s \mathbf{s} s即可,因为第三个轴 t \mathbf{t} t可以通过前两个轴叉乘的结果获得,即 t = r × s \mathbf{t = r \times s} t=r×s。有一种数值稳定的方法是:找到原始旋转轴 r \mathbf{r} r中的最小分量,然后将其设置为0,再交换剩下两个分量的值,最后再将刚才两个不为0的分量中的任意一个取反即可(这里先设置为0,交换再取反的操作,实际上就是找到与旋转轴 r \mathbf{r} r同一平面的垂直向量)。其数学表达式为:
s ‾ = { ( 0 , − r z , r y ) , if ∣ r x ∣ ≤ ∣ r y ∣ and ∣ r x ∣ ≤ ∣ r z ∣ , ( − r z , 0 , r x ) , if ∣ r y ∣ ≤ ∣ r x ∣ and ∣ r y ∣ ≤ ∣ r z ∣ , ( − r y , r x , 0 ) , if ∣ r z ∣ ≤ ∣ r x ∣ and ∣ r z ∣ ≤ ∣ r y ∣ , s = s ‾ / ∣ ∣ s ‾ ∣ ∣ , t = r × s . (4.27) \begin{array}{l}\overline{\mathbf{s}}=\left\{\begin{array}{ll}\left(0,-r_{z}, r_{y}\right), & \text { if }\left|r_{x}\right| \leq\left|r_{y}\right| \text { and }\left|r_{x}\right| \leq\left|r_{z}\right|, \\ \left(-r_{z}, 0, r_{x}\right), & \text { if }\left|r_{y}\right| \leq\left|r_{x}\right| \text { and }\left|r_{y}\right| \leq\left|r_{z}\right|, \\ \left(-r_{y}, r_{x}, 0\right), & \text { if }\left|r_{z}\right| \leq\left|r_{x}\right| \text { and }\left|r_{z}\right| \leq\left|r_{y}\right|,\end{array}\right. \\[6mm] \mathbf{s}=\overline{\mathbf{s}} /|| \overline{\mathbf{s}}||, \\[1mm] \mathbf{t}=\mathbf{r} \times \mathbf{s} .\end{array} \tag{4.27} s=⎩ ⎨ ⎧(0,−rz,ry),(−rz,0,rx),(−ry,rx,0), if ∣rx∣≤∣ry∣ and ∣rx∣≤∣rz∣, if ∣ry∣≤∣rx∣ and ∣ry∣≤∣rz∣, if ∣rz∣≤∣rx∣ and ∣rz∣≤∣ry∣,s=s/∣∣s∣∣,t=r×s.(4.27)
上述方程保证了 s ˉ \mathbf{\bar{s}} sˉ和 r \mathbf{r} r是正交的(垂直),并且 t \mathbf{t} t也与 r , s \mathbf{r,s} r,s正交,即 ( r , s , t ) \mathbf{(r,s,t)} (r,s,t)是一组标准正交基。Frisvad 提供了一种没有任何分支操作的代码实现,它的速度很快,但是相对精度较低。Max 和Duff等人改进了Frisvad方法的精度。无论使用了哪种方法进行求解,最终都会获得一组标准正交基,从而构建旋转矩阵 M \mathbf{M} M:
M = ( r T s T t T ) (4.28) \mathbf{M}= \left ( \begin{array}{} \mathbf{r}^T \\ \mathbf{s}^T \\ \mathbf{t}^T \\ \end{array} \right) \tag{4.28} M= rTsTtT (4.28)
在经过旋转矩阵 M \mathbf{M} M的变换之后,向量 r \mathbf{r} r与 x x x轴重合,向量 s \mathbf{s} s与 y y y轴重合,向量 t \mathbf{t} t与 z z z轴重合。最终获得了绕归一化旋转轴 r \mathbf{r} r旋转 α \alpha α度的变换矩阵,如下所示:
X = M T R x ( α ) M (4.29) \mathbf{X} = \mathbf{M}^T \mathbf{R}_x(\alpha) \mathbf{M} \tag{4.29} X=MTRx(α)M(4.29)
换句话说,通过变换,先将旋转轴 r \mathbf{r} r与 x x x轴重合(使用旋转矩阵 M \mathbf{M} M),然后绕 x x x轴旋转 α \alpha α度,最后再使用旋转矩阵 M \mathbf{M} M的逆矩阵进行反向变换即可。其中,由于旋转矩阵 M \mathbf{M} M是一个正交矩阵,因此其逆矩阵就是其转置矩阵 M T \mathbf{M}^{T} MT。
Goldman 还提出了另外一种绕任意归一化轴 r \mathbf{r} r旋转 ϕ \phi ϕ度的方法,这里我们直接将最终的变换矩阵展示出来:
R = ( cos ϕ + ( 1 − cos ϕ ) r x 2 ( 1 − cos ϕ ) r x r y − r z sin ϕ ( 1 − cos ϕ ) r x r z + r y sin ϕ ( 1 − cos ϕ ) r x r y + r z sin ϕ cos ϕ + ( 1 − cos ϕ ) r y 2 ( 1 − cos ϕ ) r y r z − r x sin ϕ ( 1 − cos ϕ ) r x r z − r y sin ϕ ( 1 − cos ϕ ) r y r z + r x sin ϕ cos ϕ + ( 1 − cos ϕ ) r z 2 ) (4.30) \mathbf{R} = \left(\begin{array}{ccc} \cos \phi+(1-\cos \phi) r_{x}^{2} & (1-\cos \phi) r_{x} r_{y}-r_{z} \sin \phi & (1-\cos \phi) r_{x} r_{z}+r_{y} \sin \phi \\ (1-\cos \phi) r_{x} r_{y}+r_{z} \sin \phi & \cos \phi+(1-\cos \phi) r_{y}^{2} & (1-\cos \phi) r_{y} r_{z}-r_{x} \sin \phi \\ (1-\cos \phi) r_{x} r_{z}-r_{y} \sin \phi & (1-\cos \phi) r_{y} r_{z}+r_{x} \sin \phi & \cos \phi+(1-\cos \phi) r_{z}^{2} \end{array}\right) \tag{4.30} R= cosϕ+(1−cosϕ)rx2(1−cosϕ)rxry+rzsinϕ(1−cosϕ)rxrz−rysinϕ(1−cosϕ)rxry−rzsinϕcosϕ+(1−cosϕ)ry2(1−cosϕ)ryrz+rxsinϕ(1−cosϕ)rxrz+rysinϕ(1−cosϕ)ryrz−rxsinϕcosϕ+(1−cosϕ)rz2 (4.30)