【图形学】刚体的旋转

参考：

• 基于物理的建模与动画 [美] Donald H. House（唐纳德·豪斯），John C. Keyser（约翰·凯泽） 著，叶劲峰 等 译
• 王华民 GAME103 课程
• Baraff, D. (2001). Physically based modeling: Rigid body simulation. SIGGRAPH Course Notes, ACM SIGGRAPH, 2(1), 2-1.

刚体只有两个运动：平移和旋转。
平移非常简单，就不用说了。关键是旋转。

如何描述旋转的大小？用四元数/旋转矩阵/欧拉角
三者都是描述的同一件事情，但是表示方式不同而已。
一般在图形学里就用四元数，因为它是对计算机来说计算最简单的。
（四元数的详细内容可以查看我之前讲的博客：四元数
某个点旋转绕着某个轴旋转，用四元数表示为
$$q=(cos\theta /2,sin\theta /2\hat{u})$$
$\hat{u}=（a,b,c）$确定的是转轴，是个单位向量。$\theta$角度就是旋转的角度。

如何描述旋转的速度？用角速度

$$\omega$$

这个角速度是一个三维的向量。它的长度表示旋转的角速度速率，它的方向表示旋转的转轴。

四元数转换点坐标

用四元数描述旋转，那么刚体上的每个点如何计算其位置呢？
假设任意点
$$\mathbf{r}=(x,y,z)$$
那么经过四元数旋转后就是
$$(0,\mathbf{r})=(0，\mathbf{r})q$$

更新旋转

知道了旋转的大小和速度，那么类比平移，就可以更新下一时刻的旋转

我们知道旋转用四元数q来表示，于是问题就是如何利用角速度更新q。
更新公式如下
$$q=q+\frac{\Delta t}{2}(0,\omega )q$$

注意这里有个1/2是为什么呢？
因为实际上更新q用的是q对时间的导数。而q对时间的导数有如下公式
$$\dot{q}=\frac{1}{2}(0,\omega )q$$
这点和平移有区别。

计算角速度

目前还有个问题：我如何计算角速度呢？

类比速度是由力所造成的加速度驱动的，那么角速度就是由力矩所造成的角加速度驱动的。

我们先要计算力矩，然后计算角加速度，最后更新就得到了角速度。

$$\omega =\omega +\Delta t{\mathbf{a}}_{\mathbf{w}}$$
我们这里角加速度记作${\mathbf{a}}_{\mathbf{w}}$

计算角加速度

那么如何计算角加速度呢？

类比牛顿第二定律，加速度为力
$$\mathbf{a}=\mathbf{F}/m$$

那么
$${\mathbf{a}}_{\mathbf{w}}={\mathbf{I}}^{-1}\tau$$

分子就是力矩，分母是惯性张量。

计算力矩

力矩显然也和力一样，符合叠加原理
$$\tau =\sum {\tau }_{\mathbf{i}}$$

也就是找到物体上每个质点的力矩，叠加后就是整个物体的力矩。为此我们必须先将整个物体离散成许多质点。离散得越多，精度显然就越高。

对于每个质点怎么计算力矩呢？
$${\tau }_{\mathbf{i}}={\mathbf{r}}_{\mathbf{i}}\times {\mathbf{f}}_{\mathbf{i}}$$

这里r_i是矢径，即质点到之心的距离。后面的f_i显然是每个质点的受力。

那么如何计算矢径呢？显然我们必须得先知道质心在哪。

计算质心和矢径

质心就是矢径对质量的平均
$$r_{center}=\frac{\sum r_{0i}{m}_i}{\sum m_i}$$

$r_{0i}$就是初始时刻的时候各个质点的坐标位置。

于是矢径为
$${\mathbf{r}}_{\mathbf{i}}={\mathbf{r}}_{0i}-{\mathbf{r}}_{\mathbf{center}}$$

注意质心和矢径只需要最开始计算一次就够了，因为我们的矢径是相对于质心不变的。这是因为刚体内部各个位置不发生相对位移（也就是不会变形）

计算惯性张量I

惯性张量详细的内容可以查看我之前讲过的一个博客：惯性张量
惯性张量是个二阶张量，可以以矩阵形式表示：

很遗憾的是，惯性张量和质量不同：旋转之后，它是变化的!

（更严谨的说法是在世界坐标下，它是变化的）

这就造成了很多麻烦。如果每次旋转之后都要遍历一遍所有质点计算惯性张量，那么可想而知其麻烦程度和计算量有多大！

幸运的是：我们可以借助旋转矩阵来更新惯性张量。这样，我们只要最开始计算一遍惯性张量，后面就再也不必计算了！每次只需要用旋转矩阵更新就好了。更新的公式如下

$$\mathbf{I}=\mathbf{R}{\mathbf{I}}_0{\mathbf{R}}^{\mathbf{T}}$$

这样只需要做两次矩阵乘法就够了。但是我们在此之前必须先计算旋转矩阵。

这也是为什么我们明明用四元数描述旋转，但是却还是得用旋转矩阵。

实际上，我们用到的其实是惯性张量的逆(上面那个类比牛顿第二定律的公式)。那既然这样，不如我一开始计算的就是逆矩阵，后面更新的也是逆矩阵。

$${\mathbf{I}}^{-1}=\mathbf{R}{\mathbf{I}}_0^{-1}{\mathbf{R}}^{\mathbf{T}}$$

四元数转换旋转矩阵

这里推导有些复杂，我们直接抄Wiki上的公式：

别看这个公式复杂，但是它的每一项只需要计算一次。相较于直接计算惯性张量（每一项都有对所有质点的求和），大大节约了计算量。

它的计算次数为9，另外有两次矩阵相乘（我们姑且认为每个矩阵乘法计算9次），所以总计算量为9+2*9=27次

而直接用惯性张量则每次要计算m*9次（m是质点数，看你离散的程度，越多力矩计算越精确）。假如用100个质点去离散，那么就用900次。

总结

1. 旋转用四元数表示，四元数乘以原来的3维点坐标就是旋转后的点坐标
2. 四元数的更新方式是用角速度时间积分(注意有个1/2，这是四元数对时间求导得到的)
3. 更新角速度的方式是用角加速度。类比牛顿第二定律，惯性张量的逆乘以力矩就是角加速度。为此需要即使算力矩和惯性张量。
4. 计算力矩的方式是矢径乘以各个质点的力，为此要先求出质心才能求出矢径
5. 计算惯性张量较为麻烦，要用一个矩阵遍历所有质点。而且随着旋转，惯性张量总是变的。但是幸运的是可以用旋转矩阵更新惯性张量。只要计算一次，之后不断更新就好了。
6. 为此需要将四元数转换到旋转矩阵。这也是为什么明明采用了四元数，却还是要计算旋转矩阵。这仅仅是为了节约计算量。

（王华民 GAME103 )