旋转矩阵

Posted 2021-01-25 wangguihandsome

1. 简介

计算机图形学中的应用非常广泛的变换是一种称为仿射变换的特殊变换，在仿射变换中的基本变换包括平移、旋转、缩放、剪切这几种。本文以及接下来的几篇文章重点介绍一下关于旋转的变换，包括二维旋转变换、三维旋转变换以及它的一些表达方式（旋转矩阵、四元数、欧拉角等）。

2. 绕原点二维旋转

首先要明确旋转在二维中是绕着某一个点进行旋转，三维中是绕着某一个轴进行旋转。二维旋转中最简单的场景是绕着坐标原点进行的旋转，如下图所示：

 

 

如图所示点v 绕 原点旋转θ 角，得到点v’，假设 v点的坐标是(x, y) ，那么可以推导得到 v’点的坐标（x’, y’)(设原点到v的距离是r，原点到v点的向量与x轴的夹角是? )
x=rcos?y=rsin?
x′=rcos(θ+?)y′=rsin(θ+?)
通过三角函数展开得到
x′=rcosθcos??rsinθsin?
y′=rsinθcos?+rcosθsin?
带入x和y表达式得到
x′=xcosθ?ysinθ
y′=xsinθ+ycosθ
写成矩阵的形式是：

 

尽管图示中仅仅表示的是旋转一个锐角θ的情形，但是我们推导中使用的是三角函数的基本定义来计算坐标的，因此当旋转的角度是任意角度（例如大于180度，导致v’点进入到第四象限）结论仍然是成立的。

3. 绕任意点的二维旋转

绕原点的旋转是二维旋转最基本的情况，当我们需要进行绕任意点旋转时，我们可以把这种情况转换到绕原点的旋转，思路如下：
1. 首先将旋转点移动到原点处
2. 执行如2所描述的绕原点的旋转
3. 再将旋转点移回到原来的位置

 

也就是说在处理绕任意点旋转的情况下需要执行两次平移的操作。假设平移的矩阵是T(x,y)，也就是说我们需要得到的坐标 v’=T(x,y)*R*T(-x,-y)（我们使用的是列坐标描述点的坐标，因此是左乘，首先执行T(-x,-y)）

在计算机图形学中，为了统一将平移、旋转、缩放等用矩阵表示，需要引入齐次坐标。（假设使用2x2的矩阵，是没有办法描述平移操作的，只有引入3x3矩阵形式，才能统一描述二维中的平移、旋转、缩放操作。同理必须使用4x4的矩阵才能统一描述三维的变换）。

对于二维平移，如下图所示，P点经过x和y方向的平移到P’点，可以得到：

 

 

x′=x+tx

y′=y+ty
由于引入了齐次坐标，在描述二维坐标的时候，使用（x，y，w）的方式（一般w=1），于是可以写成下面矩阵的形式

按矩阵乘法展开，正好得到上面的表达式。也就是说平移矩阵是

 

如果平移值是（-tx，-ty）那么很明显平移矩阵式

 

 

我们可以把2中描述的旋转矩阵也扩展到3x3的方式，变为：

 

从平移和旋转的矩阵可以看出，3x3矩阵的前2x2部分是和旋转相关的，第三列与平移相关。有了上面的基础之后，我们很容易得出二维中绕任意点旋转的旋转矩阵了，只需要把三个矩阵乘起来即可：

 

4. 三维基本旋转

我们可以把一个旋转转换为绕基本坐标轴的旋转，因此有必要讨论一下绕三个坐标值x、y、z的旋转。
本文在讨论过程中使用的是类似于OpenGL中定义的右手坐标系，同时旋转角度的正负也遵循右手坐标系的约定。如下图所示

 

 

4.1 绕X轴的旋转

在三维场景中，当一个点P(x,y,z)绕x轴旋转θ角得到点P’(x’,y’,z’)。由于是绕x轴进行的旋转，因此x坐标保持不变，y和z组成的yoz（o是坐标原点）平面上进行的是一个二维的旋转，可以参考上图（y轴类似于二维旋转中的x轴，z轴类似于二维旋转中的y轴），于是有：
x′=x
y′=ycosθ?zsinθ
z′=ysinθ+zcosθ
写成（4x4）矩阵的形式

 

4.2 绕Y轴旋转

绕Y轴的旋转和绕X轴的旋转类似，Y坐标保持不变，除Y轴之外，ZOX组成的平面进行一次二维的旋转（Z轴类似于二维旋转的X轴，X轴类似于二维旋转中的Y轴，注意这里是ZOX，而不是XOZ，观察上图中右手系的图片可以很容易了解到这一点），同样有：
x′=zsinθ+xcosθ
y′=y
z′=zcosθ?xsinθ
写成（4x4）矩阵的形式

 

 

4.3 绕Z轴旋转

与上面类似，绕Z轴旋转，Z坐标保持不变，xoy组成的平面内正好进行一次二维旋转（和上面讨论二维旋转的情况完全一样）

 

 

4.4 小结

上面描述了三维变换中绕单一轴旋转的矩阵表达形式，绕三个轴旋转的矩阵很类似，其中绕y轴旋转的矩阵与绕x和z轴旋转的矩阵略有点不同（主要是三个轴向顺序和书写矩阵的方式不一致导致的，绕三个不同坐标旋转轴以及其他二个坐标轴组成平面的顺序是： XYZ(绕x轴） YZX（绕y轴） ZXY（绕z轴），其中绕y轴旋转，其他两个轴是ZX，这和我们书写矩阵按

 

 

 

的方式不一致，而导致看起来绕Y轴旋转的矩阵似乎是和其他两个矩阵不一致。如果我们颠倒写法，将公式写成

 

的方式，那么这三个旋转矩阵看起来在形式上就统一了，都是

 

 

这种表现形式了（左上角都是?sinθ）

5. 绕任意轴的三维旋转

绕任意轴的三维旋转可以使用类似于绕任意点的二维旋转一样，将旋转分解为一些列基本的旋转。绕任意轴旋转如下图所示：

 

 

P点绕向量u旋转θ角，得到点Q，已知P点的坐标和向量u，如何求Q点的坐标。
我们可以把向量u进行一些旋转，让它与z轴重合，之后旋转P到Q就作了一次绕Z轴的三维基本旋转，之后我们再执行反向的旋转，将向量u变回到它原来的方向，也就是说需要进行的操作如下：
1. 将旋转轴u绕x轴旋转至xoz平面
2. 将旋转轴u绕y轴旋转至于z轴重合
3. 绕z轴旋转θ角
4. 执行步骤2的逆过程
5. 执行步骤1的逆过程
原始的旋转轴u如下图所示：

 

第1、2、3步骤如下图所示：

 

 

 

步骤1将向量u旋转至xoz平面的操作是一个绕x轴的旋转操作，步骤2将向量u旋转到与z轴重合，第1、2步骤的示意图如下：

 

作点P在yoz平面的投影点q，q的坐标是（0, b, c)，原点o与q点的连线oq和z轴的夹角就是u绕x轴旋转的角度。通过这次旋转使得u向量旋转到xoz平面（图中的or向量）【步骤1】
过r点作z轴的垂线，or与z轴的夹角为β, 这个角度就是绕Y轴旋转的角度，通过这次旋转使得u向量旋转到与z轴重合【步骤2】

步骤1中绕x轴旋转的是一次基本的绕x轴的三维旋转，按照之前的讨论，旋转矩阵是：

 

 

这里的θ就是图中所示的α角 （注意α角度是绕x旋转的正的角度）
从图中我们还可以得到：

于是旋转矩阵(记作 Rx(α))为：

 

在完成步骤1之后，向量u被变换到了r的位置，我们继续步骤2的操作，绕y轴旋转负的β角（注意：这里的β是负的），经过这次变换之后向量u与z轴完全重合，由于这一步也是执行的一次绕Y轴的基本旋转，旋转矩阵（记作 Ry(?β))为：

 

 

使用?β替换表达式中的θ，此外根据图中描述，我们可以计算得到：

 

 

带入上面的表达式，于是旋转矩阵(记作 Ry(?β))为：

 

在完成前面两个步骤之后，u方向和z轴完全重合，因此执行旋转θ角，执行的是一次绕z轴的基本三维旋转（记作 R(θ)，根据之前的讨论，我们可以得到：

最后两步骤是前面1和2的逆操作，也就是绕Y轴旋转β 和绕X轴旋转?α，这两个矩阵分别记作 Ry(β) 和 Rx(?α)，得到它们的方式很简单，只需要将上面步骤1和步骤2中的角度修改成相反数即可，也就是：

 

 

最终得到 绕任意轴u旋转的旋转矩阵是【因为使用的列向量，因此执行的是左乘（从右往左）】：

MR=Rx(?α)Ry(β)Rz(θ)Ry(?β)Rx(α)=

 

 

(注：式中的(u,v,w)对应上文中向量(a,b,c)，公式我自己笔算过，为了减少编辑公式的时间（使用LaTex编辑太繁琐，因此找了一张公式的图片贴在此处)

如果向量是经过单位化的（单位向量），那么有a2+b2+c2=1，可以简化上述的公式，得到：

 

6绕任意轴旋转

绕任意轴旋转的情况比较复杂，主要分为两种情况，一种是平行于坐标轴的，一种是不平行于坐标轴的，对于平行于坐标轴的，我们首先将旋转轴平移至与坐标轴重合，然后进行旋转，最后再平移回去。

• 将旋转轴平移至与坐标轴重合，对应平移操作
• 旋转，对应操作
• 步骤1的逆过程，对应操作

整个过程就是

对于不平行于坐标轴的，可按如下方法处理。（该方法实际上涵盖了上面的情况）

1. 将旋转轴平移至原点
2. 将旋转轴旋转至YOZ平面
3. 将旋转轴旋转至于Z轴重合
4. 绕Z轴旋转θ度
5. 执行步骤3的逆过程
6. 执行步骤2的逆过程
7. 执行步骤1的逆过程

假设用v1(a1, b2, c2)和v2(a2, b2, c2)来表示旋转轴，θ表示旋转角度。为了方便推导，暂时使用右手系并使用列向量，待得出矩阵后转置一下即可，上面步骤对应的流程图如下。

步骤1是一个平移操作，将v1v2平移至原点，对应的矩阵为

步骤2是一个旋转操作，将p(p = v2 -v1)旋转至XOZ平面，步骤3也是一个旋转操作，将p旋转至与Z轴重合，这两个操作对应的图如下。

做点p在平面YOZ上的投影点q。再过q做Z轴垂线，则r是p绕X轴旋转所得，且旋转角度为α，且

,   

于是旋转矩阵为

现在将r绕Y轴旋转至与Z轴重合，旋转的角度为-beta（方向为顺时针），且

,    

于是得到旋转矩阵为

最后是绕Z轴旋转，对应的矩阵如下

如果旋转轴是过原点的，那么第一步和最后一步的平移操作可以省略，也就是把中间五个矩阵连乘起来，再转置一下，得到下面的绕任意轴旋转的矩阵

即

对应的函数代码如下

 1 void RotateArbitraryAxis(D3DXMATRIX* pOut, D3DXVECTOR3* axis, float theta)
 2 {
 3     D3DXVec3Normalize(axis, axis);
 4     float u = axis->x;
 5     float v = axis->y;
 6     float w = axis->z;
 7 
 8     pOut->m[0][0] = cosf(theta) + (u * u) * (1 - cosf(theta));
 9     pOut->m[0][1] = u * v * (1 - cosf(theta)) + w * sinf(theta);
10     pOut->m[0][2] = u * w * (1 - cosf(theta)) - v * sinf(theta);
11     pOut->m[0][3] = 0;
12 
13     pOut->m[1][0] = u * v * (1 - cosf(theta)) - w * sinf(theta);
14     pOut->m[1][1] = cosf(theta) + v * v * (1 - cosf(theta));
15     pOut->m[1][2] = w * v * (1 - cosf(theta)) + u * sinf(theta);
16     pOut->m[1][3] = 0;
17 
18     pOut->m[2][0] = u * w * (1 - cosf(theta)) + v * sinf(theta);
19     pOut->m[2][1] = v * w * (1 - cosf(theta)) - u * sinf(theta);
20     pOut->m[2][2] = cosf(theta) + w * w * (1 - cosf(theta));
21     pOut->m[2][3] = 0;
22 
23     pOut->m[3][0] = 0;
24     pOut->m[3][1] = 0;
25     pOut->m[3][2] = 0;
26     pOut->m[3][3] = 1;

如果旋转轴是不过原点的，那么第一步和最后一步就不能省略，将所有七个矩阵连乘起来，得到如下变换矩阵

对应如下这个超长的矩阵，在这里(u, v, w) = (a2, b2, c2) - (a1, b1, c1)，且是单位向量，a, b, c分别表示(a1, b1, c1)

将上面的过程写成函数，该函数接受四个参数，第一个参数是一个输出参数，用来保存得到的旋转矩阵，第二个和第三个参数是旋转轴的两个端点，最后一个参数是旋转角度θ，注意，在函数中我们已经将上面的矩阵转置了，因为上面是按照列向量计算的。

 1 void RotateArbitraryLine(D3DXMATRIX* pOut, D3DXVECTOR3* v1, D3DXVECTOR3* v2, float theta)
 2 {
 3     float a = v1->x;
 4     float b = v1->y;
 5     float c = v1->z;
 6 
 7     D3DXVECTOR3 p = *v2 - *v1;
 8     D3DXVec3Normalize(&p, &p);
 9     float u = p.x;
10     float v = p.y;
11     float w = p.z;
12 
13     float uu = u * u;
14     float uv = u * v;
15     float uw = u * w;
16     float vv = v * v;
17     float vw = v * w;
18     float ww = w * w;
19     float au = a * u;
20     float av = a * v;
21     float aw = a * w;
22     float bu = b * u;
23     float bv = b * v;
24     float bw = b * w;
25     float cu = c * u;
26     float cv = c * v;
27     float cw = c * w;
28 
29     float costheta = cosf(theta);
30     float sintheta = sinf(theta);
31 
32     pOut->m[0][0] = uu + (vv + ww) * costheta;
33     pOut->m[0][1] = uv * (1 - costheta) + w * sintheta;
34     pOut->m[0][2] = uw * (1 - costheta) - v * sintheta;
35     pOut->m[0][3] = 0;
36 
37     pOut->m[1][0] = uv * (1 - costheta) - w * sintheta;
38     pOut->m[1][1] = vv + (uu + ww) * costheta;
39     pOut->m[1][2] = vw * (1 - costheta) + u * sintheta;
40     pOut->m[1][3] = 0;
41 
42     pOut->m[2][0] = uw * (1 - costheta) + v * sintheta;
43     pOut->m[2][1] = vw * (1 - costheta) - u * sintheta;
44     pOut->m[2][2] = ww + (uu + vv) * costheta;
45     pOut->m[2][3] = 0;
46 
47     pOut->m[3][0] = (a * (vv + ww) - u * (bv + cw)) * (1 - costheta) + (bw - cv) * sintheta;
48     pOut->m[3][1] = (b * (uu + ww) - v * (au + cw)) * (1 - costheta) + (cu - aw) * sintheta;
49     pOut->m[3][2] = (c * (uu + vv) - w * (au + bv)) * (1 - costheta) + (av - bu) * sintheta;
50     pOut->m[3][3] = 1;
51 }

相关学习拓展

https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle

http://blog.miskcoo.com/2018/03/isometry_and_quaternion