GAMES103课程笔记04-Rigid Body Contacts
这个系列是GAMES103：基于物理的计算机动画入门(GAMES 103: Intro to Physics Based Animation)的同步课程笔记。本课程是计算机动画的入门课程，着重介绍各种基于物理的动画仿真模拟技术。本节课主要介绍刚体的碰撞与接触。
Particle Collision
Signed Distance Function
空间中的曲面可以使用有向距离函数(signed distance function, SDF)来表示。根据距离函数\(\phi(\mathbf{x})\)的符号我们可以表示质点\(\mathbf{x}\)在曲面外侧、曲面内侧或是曲面上：
常见规则曲面的SDF如下图所示：
使用SDF来表示曲面的好处之一在于它可以高效地表示质点与多个曲面的位置关系。假设我们有3个曲面\(\phi_0\)、\(\phi_1\)、\(\phi_2\)，当质点\(\mathbf{x}\)位于3个曲面内部时它到最近曲面的距离为：
\[\phi (\mathbf{x}) = \max \{ \phi_0(\mathbf{x}), \phi_1(\mathbf{x}), \phi_2(\mathbf{x}) \}\] \[\phi_i (\mathbf{x}) \lt 0\]
换句话说我们只需要计算若干个距离函数就可以判断质点是否与曲面发生碰撞，这就为碰撞检测带来了便利。
当然需要说明的是当质点位于曲面并集构成的区域中\((\phi_i (\mathbf{x}) \lt 0)\)时使用上式的计算结果是不精确的：
\[\phi(\mathbf{x}) \approx \max\{\phi_0 (\mathbf{x}), \phi_0 (\mathbf{x}) \}\]
但是如果\(\mathbf{x}\)位于曲面外侧\((\phi_i (\mathbf{x}) \gt 0)\)则可以精确计算质点到曲面的距离：
\[\phi(\mathbf{x}) = \min \{\phi_0 (\mathbf{x}), \phi_0 (\mathbf{x}) \}\]Penalty Method
碰撞仿真的基本方法是利用SDF进行碰撞检测，当质点进入曲面内部时施加一个反方向的力来把它推出去，这样的方法称为penalty method：
施加推力时可以把它想象成弹簧，这样力的大小就和压缩的距离成正比，力的方向为\(\phi(\mathbf{x})\)的梯度方向。显然这样的方法容易造成穿模的问题，因此实际使用时往往还需要设计一个安全距离\(\varepsilon\)，只要质点进入了安全距离内就施加推力：
然而这样方法存在一个问题，当弹簧系数\(k\)比较大时会产生过于巨大的推力。因此实践中也会使用距离的对数来对\(k\)进行近似，这样当距离很小的时候才会产生非常大的推力。但需要说明的是这种方法在质点发生穿模时会产生吸力，在实际应用中必须要调整实践步长\(\Delta t\)来避免这种情况的发生。
Impulse Method
penalty method通过施加推力来模拟质点的碰撞行为，除此之外还可以直接修改质点的运动状态来实现碰撞仿真，这种方法称为impulse method。在impulse mothod中，如果检测到了碰撞首先要更新质点的位置，将它移动到曲面上：
接下来需要修改质点的速度，我们将速度分解为曲面法向和曲面切向两部分\(\mathbf{v}_{\mathbf{N}}\)和\(\mathbf{v}_{\mathbf{T}}\)。然后在两个方向上对速度进行衰减，这里需要注意对法向速度取反。最后再将衰减后的两个速度相加即可得到更新后的质点速度。
Rigid Body Collision
Collision Detection
我们可以把刚体视为质点的集合，对刚体上的每个点进行遍历并分别进行碰撞检测来实现刚体的碰撞检测。尽管这种方法并不完美，当仍然可以得到不错的结果。
Collision Response by Impulse
假设我们检测到了点\(\mathbf{x}_i\)发生碰撞，则该处的位置和速度为：
\[\mathbf{x}_i = \mathbf{x} + \mathbf{R} \mathbf{r}_i\] \[\mathbf{v}_i = \mathbf{v} + \boldsymbol{\omega} \times \mathbf{R} \mathbf{r}_i\]我们可以利用质点的碰撞公式得到\(\mathbf{x}_i\)碰撞后的速度，但需要注意的是\(\mathbf{x}_i\)仍然位于刚体上，整个刚体需要具有相同的速度和角速度。因此我们需要利用质点碰撞后的速度来修正刚体的速度，具体来说我们假定碰撞后刚体受到位于\(\mathbf{x}_i\)处的冲量\(\mathbf{j}\)的作用，根据刚体动力学公式可以得到碰撞后的速度为：
\[\mathbf{v}^\text{new} = \mathbf{v} + \frac{1}{M} \mathbf{j}\] \[\boldsymbol{\omega}^\text{new} = \boldsymbol{\omega} + \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j})\]
此时\(\mathbf{x}_i\)处的速度为：
\[\begin{aligned} \mathbf{v}_i^\text{new} &= \mathbf{v}^\text{new} + \boldsymbol{\omega}^\text{new} \times \mathbf{R} \mathbf{r}_i \\ &= \mathbf{v} + \frac{1}{M} \mathbf{j} + \big(\boldsymbol{\omega} + \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j}) \big) \times \mathbf{R} \mathbf{r}_i \\ &= (\mathbf{v} + \boldsymbol{\omega} \times \mathbf{R} \mathbf{r}_i) + \frac{1}{M} \mathbf{j} + \big( \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j}) \big) \times \mathbf{R} \mathbf{r}_i \\ &= \mathbf{v}_i + \frac{1}{M} \mathbf{j} + \big( \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j}) \big) \times \mathbf{R} \mathbf{r}_i \\ &= \mathbf{v}_i + \frac{1}{M} \mathbf{j} - (\mathbf{R} \mathbf{r}_i) \times \big( \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j}) \big) \end{aligned}\]根据叉乘的性质，向量之间的叉乘可以表示为反对称矩阵和向量的乘法：
\[\mathbf{p} \times \mathbf{q} = \begin{bmatrix} p_y q_z - p_z q_y \\ p_z q_x - p_x q_z \\ p_x q_y - p_y q_x \\ \end{bmatrix} = \begin{bmatrix} 0 & -p_z & p_y \\ p_z & 0 & -p_x \\ -p_y & p_x & 0 \\ \end{bmatrix} \begin{bmatrix} q_x \\ q_y \\ q_z \end{bmatrix} = \mathbf{p}^* \mathbf{q}\] \[\mathbf{p}^* = \begin{bmatrix} 0 & -p_z & p_y \\ p_z & 0 & -p_x \\ -p_y & p_x & 0 \\ \end{bmatrix}\]我们可以把更新后的速度表示为矩阵乘法的形式：
\[\begin{aligned} \mathbf{v}_i^\text{new} &= \mathbf{v}_i + \frac{1}{M} \mathbf{j} - (\mathbf{R} \mathbf{r}_i) \times \big( \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i \times \mathbf{j}) \big) \\ &= \mathbf{v}_i + \frac{1}{M} \mathbf{j} - (\mathbf{R} \mathbf{r}_i)^* \ \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i)^* \ \mathbf{j} \\ &= \mathbf{v}_i + \mathbf{K} \mathbf{j} \end{aligned}\] \[\mathbf{K} = \frac{1}{M} \mathbf{1} - (\mathbf{R} \mathbf{r}_i)^* \ \mathbf{I}^{-1} (\mathbf{R} \mathbf{r}_i)^*\]其中\(\mathbf{1}\)为单位矩阵。上式说明我们可以利用\(\mathbf{x}_i\)点碰撞后的速度反解出作用在刚体上的冲量\(\mathbf{j}\)，然后再利用冲量\(\mathbf{j}\)更新刚体的速度即可实现碰撞的仿真。整个求解过程如下：
对于刚体之间碰撞的情况也可以使用类似的方法进行求解：
Shape Matching
本节课最后介绍了使用shape matching来解决刚体碰撞问题的方法。和前面介绍过的方法相比，shape matching更多地使用了数学优化而非物理定律来处理碰撞问题。shape matching的思路是忽略刚体的约束，将刚体上的每个点视为独立的质点来求解，每个质点的运动状态更新后再施加刚体约束来得到刚体运动后的状态。
具体来说，假设刚体质心的位置\(\mathbf{c}\)和旋转\(\mathbf{R}\)为待求量，则shape matching对应的优化问题可以表示为：
\[\{ \mathbf{c}, \mathbf{R} \} = \arg \min \sum_i \frac{1}{2} \Vert \mathbf{c} + \mathbf{R} \mathbf{r}_i - \mathbf{y}_i \Vert^2\]
我们定义优化的目标函数为\(E = \sum_i \frac{1}{2} \Vert \mathbf{c} + \mathbf{R} \mathbf{r}_i - \mathbf{y}_i \Vert^2\)，首先对质心位置求导得到：
\[\frac{\partial E}{\partial \mathbf{c}} = \sum_i \mathbf{c} + \mathbf{R} \mathbf{r}_i - \mathbf{y}_i = \sum_i \mathbf{c} - \mathbf{y}_i = \mathbf{0}\] \[\mathbf{c} = \frac{1}{N} \sum_i \mathbf{y}_i\]即更新后质心的位置为所有质点的平均值。接下来考虑旋转，我们暂时忽略旋转矩阵的约束只把它当成通常的矩阵\(\mathbf{A}\)进行求导，可以得到：
\[\frac{\partial E}{\partial \mathbf{A}} = \sum_i (\mathbf{c} + \mathbf{A} \mathbf{r}_i - \mathbf{y}_i) \mathbf{r}_i^T = 0\] \[\mathbf{A} = \bigg( \sum_i (\mathbf{y}_i - \mathbf{c}) \mathbf{r}_i^T \bigg) \bigg( \sum_i \mathbf{r}_i \mathbf{r}_i^T \bigg)^{-1}\]然后对\(\mathbf{A}\)进行极分解(polar decomposition)就可以获得所需的旋转矩阵\(\mathbf{R}\)：
\[\mathbf{A} = \mathbf{R} \mathbf{S}\]从SVD的角度来看，极分解的可以理解为在进行第二次旋转前将物体旋转回初始的位置，这样线性变换\(\mathbf{A}\)就可以分解为全局旋转\(\mathbf{R} = \mathbf{U} \mathbf{V}^T\)和本地形变\(\mathbf{S} = \mathbf{V} \mathbf{D} \mathbf{V}^T\)：
总结一下，使用shape matching来求解刚体碰撞问题的流程如下：
shape matching除了可以解决刚体运动外还可以处理布料、弹性体甚至流体的运动仿真。但需要注意的是shape matching很难处理摩擦的问题，因此shape matching比较适合摩擦力可以忽略的场景。