GAMES103 (二) rigid body

---

游戏引擎中的刚体被认为是不会发生形变的物体，在物理引擎中，我们对刚体的模拟是通过状态函数来实现的，通过确定不同时刻刚体的速度和位置，从而可以得到刚体实时的状态

在物理模拟中，刚体模拟可以拆分为两个部分：平移和旋转

一、平移运动

依托最原始的速度公式
$$v_{i+1}=v_i+a\ast t$$
$$a=\frac{f}{m}$$
其中 $f$ 可以表示为一个关于速度、位移和时间的函数，设置 $t=1$

由于我们目前所只是已知 $t^{[0]}$ 的状态，所以对 $t^{[0]}$ 进行的积分为显式积分，对 $t^{[1]}$ 进行积分为隐式积分，使用泰勒展开得到的误差都在二阶范围内,也就是一阶准确。

将两种方法叠加，可以得到 Mid-Point 方法，三阶误差所以更精确。

- 除了将区域定义为一个矩形，还有方法将区域定义为其它形状，例如梯形。

因此我们最开始关于 $x$、$t$ 的方程可以化作为一个显式一个隐式的方式，将其中 $v^{[1]}$ 写成 $v^{[0.5]}$ ，虽然 $x^{[1]}$的求法没变，但可以看做是对其使用了 Mid-Point 的方法

对于重力一般就是公式：

对于空气阻力如果要求精确，可以使用
$$f_{drag}^{[0]}=-\sigma v^{[0]}$$
不过现在流行的都是直接对速度进行一个衰减，α为衰减系数（老师较为推荐）
$$v^{[1]}=\alpha v^{[0]}$$

二、旋转运动

旋转矩阵的好处与坏处：

– 矩阵中有九个元素，旋转却只有三个自由度（前后、左右、上下）

欧拉角的好处与坏处

为了解决可以表示三维空间中的点，并且可以像复数那样表示加减乘除，因此发明了四元数（三维向量只有加、减、乘没有除）

通过对四元数的使用，可以得到代表绕着 $v$ 轴旋转 $\theta$ 度的四元数 $q$，并可以将四元数转换成矩阵形式：

至于其中的 $q$ 为什么定义为 $cos\frac{\theta }{2}$，个人朴素的理解是绕着 $\theta$ 转动一半的过程中水平方向的作用是不一致的，所以将一个 $\theta$ 的旋转拆成两次。

note: $q*$ 是指 $q$ 的共轭四元数

关于四元数到旋转矩阵的推理（使用 $a^2+b^2+c^2+d^2=1$ 合并）

在旋转中将力称为力矩$(torque)$ $\vec{\tau }$，将质量称为惯性$(inertia)$ $\vec{I}$。

$\vec{\tau }$ 被认为是关于力产生的对旋转的趋势，它与力的方向和旋转半径的方向都是垂直的，由于 $\vec{\tau }$ 的大小与旋转半径和力之间的夹角有关，所以关于 $\vec{\tau }$ 的公式用叉乘表示：

力矩对应于力的旋转，正如力引起线性动量的变化一样。我们知道力$(\vec{F})$可以表示为动量$(\vec{P})$对时间的导数，同样，力矩$(\vec{\tau })$ 可以表示为角动量$(\vec{\omega })$的导数

对于惯性张量矩阵，直接给出定义（$1$ 是单位向量）：

需要注意的是惯性张量描述的是对抗力矩所产生的旋转趋势的力大小，由于力矩具有方向，所以不同方向的惯性张量是不一样的，这与其质量分布有关：

在此科普下这个公式的推导，由质点的角动量公式可知：
$${\vec{L}}_i={\vec{r}}_i\times {\vec{p}}_i$$
依据质点动量和线速度的定义：
$${\vec{p}}_i=m_i\times {\vec{v}}_i$$
$${\vec{v}}_i=\vec{\omega }\times {\vec{r}}_i$$
可以得到最终的 ${\vec{L}}_i$ 表达式：
$${\vec{L}}_i=m_i{\vec{r}}_i\times (\vec{\omega }\times {\vec{r}}_i)$$
使用公式：
$$A\times (B\times C)=B(A\cdot C)-C(A\cdot B)$$
鉴于 ${\vec{r}}_i$ 是一个列向量：
$${\vec{r}}_i=\left[\begin{array}{c}{x}_i\\ y_i\\ z_i\end{array}\right]$$
可以得到：
$${\vec{L}}_i=m_i(\vec{\omega }{\vec{r}}_i^T{\vec{r}}_i-{\vec{r}}_i{\vec{r}}_i^T\vec{\omega })$$
将 $\vec{\omega }$ 提出去，可以得到
$${\vec{L}}_i=[m_i({\vec{r}}_i^T{\vec{r}}_i-{\vec{r}}_i{\vec{r}}_i^T)]\vec{\omega }$$
$${\vec{L}}_i=I_i\vec{\omega }$$
将中括号中的部分提取出来即可得到等式，除此之外，也可推出角动量对时间的导数即为 $\tau$，所以力矩可以类比平移运动中的力
$$\frac{d\vec{L}}{dt}=\tau$$

由于惯性张量矩阵 $I$ 是一个与时间 $t$ 无关的函数，类比于冲量除以质量等于增加的速度：
$$\Delta v=\frac{F\Delta t}{m}$$
因此也可以得到角速度的增量公式为：
$$\Delta \omega =\frac{\tau \Delta t}{I}$$

参考自：
动态分析中的惯性矩阵（张量）基础概念温习
Lecture L26 - 3D Rigid Body Dynamics: The Inertia Tensor

由于$r_i=R{r}_{ref}$，所以可得到对参考坐标系的惯性张量

以下是详细证明：

参考自：Physically Based Modeling Rigid Body Simulation

对于旋转部分，四元数的实数部分代表的是旋转角度，虚数部分是旋转轴方向。由于四元数的导数为：
$$\frac{d\vec{q}}{dt}=\frac{1}{2}\vec{\omega }\vec{q}$$
其中的 $\omega$ 必须转化为 $[0,\omega ]$ 的形式，由于我们四元数累加的意义就是代表积分，所以可以先相乘得到一个对当前旋转的变化表示，再将该变化添加到原来的变化（用四元数 $q^{[i]}$ 表示）。
$${\vec{q}}^{[i+1]}={\vec{q}}^{[i]}+[0,\frac{\vec{\omega }t}{2}]\times {\vec{q}}^{[i]}$$

最后是写代码时候的一些小 $trick$

三、质点碰撞检测

可以用一个距离向量函数 $SDF(SignedDistanceFunction)$ 来表示一个平面，通过符号来表示某个点与该平面的关系。

当负号为负数时，要注意函数 $\varphi (x)$ 越小的表示的距离越大，当然如果表示都在外面，由于未发生碰撞所以就可以不用去关心数值

当检测交集内的点时，$\varphi (x)$ 可能失灵，如下图点与平面的距离既不是 ${\varphi }_0(x)$ 也不是 ${\varphi }_1(x)$，而是两个圆圈的交点。

四、质点碰撞响应

此部分完全照抄了博主写的部分：GAMES103笔记 Lecture4 刚体碰撞（Rigid Body Contacts），如有冒犯，还请联系告知删除

1. Penalty Method (惩罚模型)

检测到质点进入物体，即 $\varphi (x)<0$ 时，给物体施加向外的力 $f=-k\varphi N$ 。

我们希望质点以最快的方向离开物体，所以 $f$ 的方向与 $\nabla \varphi (x)$ 相同。对于平面来说就是法向量的方向。

力的大小和基于胡克定律的弹簧模型是一样的，相当于在物体越过边界后拥有了与距离的绝对值的二次方成正比的弹性势能，所以称该方法为Quadratic Penalty方法。

Quadratic Penalty方法显然存在一个问题，物体进去之后才施加向外的力被弹出来，容易出现穿模现象。

直接有效的解决办法是加个缓冲，只要靠近物体就开始施加向外的力。

加缓冲后穿模的问题得到缓解，但在实践中仍然存在有问题。 $k$ 太小会导致物体速度改变太慢，还是会出现进去以后才出来的情况，而 $k$ 太大的话将导致物体出来速度太大（overshooting）。

可以想到，$k$ 的大小不好调节的本质原因就因为是这是一个类似弹簧的模型，力的大小始终和距离的绝对值成正比，而现实世界显然不是这样的。

Log-Barrier Penalty 方法更接近现实世界的物理规律。

将 SDF 放在分母上，再用一个参数 $\rho$ 来控制力的大小。在距离很小的情况下，受到的力将会增长到非常大的数值。

如果将点在无穷远处的能量看做 $0$ ，将力函数对位移求积分算能量 $E$ 会发现与$log[(\varphi (x)]$ 成正比所以称此方法为 Log-Barrier Penalty 方法。

$$E=-\int fd[\varphi (x)]=-\int \frac{\rho \vec{N}}{\varphi (x)}d[\varphi (x)]=-\rho \vec{N}log[(\varphi (x)]$$

Log-Barrier Penalty 方法假设了 $\varphi (x)$ 永远都不会小于0，虽然很符合物理规律，然而在计算机中我们只能一个时间步一个时间步地计算，如果在一个时间步内物体就穿过了表面，力的方向就反过来了。此外，Log-Barrier Penalty 方法仍然存在overshooting的问题。所以 Log-Barrier Penalty 方法要求时间步步长必须足够小。

Penalty方法的优势是容易实现，碰撞后用力改变物体运动的方式符合真实的物理规律，缺点是需要调节到合适的时间步步长来避免overshooting，而且难以模拟斜着撞过来的物体的水平速度因摩擦力衰减。

2. Impulse Method（冲量模型）

除了从底层模拟物理规律，还可以从表象上来模拟物体的运动。

Penalty 方法施加的力需要经过一小段时间才能起作用，Impulse（冲量）方法在检测到碰撞的瞬间就改变物体的位置和动量（速度）。

检测到进入物体，即 $\varphi (x)<0$ 。首先立即将物体沿着最近的f方向移动到物体表面并更新 $x^{new}$

然后继续检查速度是否指向物体内部，如果 $v\cdot N\ge 0$ ，说明在之前的时间步中已经修改过质点的速度，它已经要出去了，不用做任何处理。反之将速度正交分解再更新 $v_N、v_T$。

其中 $a,{\mu }_N\subset [0,1]$ ，因为我们希望质点离开碰撞的物体，竖直方向上的速度反向要反过来，但因为能量守恒定律，出去的速度的绝对值不能超过进入时的速度。

由于摩擦力的存在，水平方向上的速度需要衰减，依据摩擦力定律（Amontons-Coulomb Friction Laws）：
$$F_f\le {\mu }_T{F}_N$$

可以算出 $a$ 的最大值，但是这个值不能为负的，否则就会让摩擦力方向转向了

Impulse方法的优点是能做到对摩擦力和反弹的精确控制，但对于非质点物体需要做一些额外的处理。

刚体用Impulse方法更多，衣服和弹性体用Penalty方法更多，因为他们对摩擦力和反弹的表示不需要很精确

五、刚体碰撞检测和响应

1、碰撞检测

可以暴力地解决碰撞检测的问题。遍历刚体表面上的每一个点看 SDF 是否小于0，即可判断是否发生碰撞。作业中是用点的平均值来进行
是否可以考虑使用凸包问题来得到最接近撞击面的点

2、 碰撞响应

2.1、Impulse method

首先是可以通过刚体的角速度和线速度来得到碰撞点（相对于刚体）的位置 $x_i$ 和线速度 $v_i$（这里刚体的速度属性其实是刚体质心点的速度属性）

其中碰撞点线速度的公式为什么使用叉乘，可以通过叉乘的方向来理解，碰撞点的线速度方向 $v_i$ 、旋转半径 $r_i$ 以及角速度 $\omega$ 互相垂直。

有了碰撞点的线速度 $v_i$ 就可以得到碰撞点碰撞后的新速度 $v_i^{new}$，它是由刚体的线速度和碰撞点相对于刚体的线速度合成得到的。由我们之前推得的角速度、线速度增量公式：
$$\Delta v=\frac{F\Delta t}{m}$$
$$\Delta \omega =\frac{\tau \Delta t}{I}$$

此时假设有一个冲量 $j$ 作用在我们的碰撞点上，我们假设作用在碰撞点和作用在质心上对整个刚体速度改变的效果一样，就可以得到碰撞点的新速度，也即刚体整体的新线速度 $v^{new}$ 和角速度 ${\omega }^{new}$，从而得到碰撞点最终的速度 $v_i^{new}$，将等式中的 $v$ 和 $\omega \times R{r}_i$ 提取出来得到 $v_i$

将叉乘写成矩阵形式，最终碰撞点在受到冲量 $j$ 撞击前后的速度差可以用冲量与一个矩阵 $K$ 相乘表示出来。

从而可以得到反算出冲量 $j$:
$$j=K^{-1}(v_i^{new}-v_i)$$

从而可以更新刚体整体的速度：

总结一下：$Impulse$ 进行处理碰撞的方法是先判断碰撞点的位置是否进入了 $\varphi (x)$，进入的话就把碰撞点的位置移到 $\varphi (x)$ 表面，再看碰撞点的速度方向是否还在冲着深入 $\varphi (x)$ 的方向前进，如果是的话，就把他的速度进行分解，平行于 $\varphi (x)$的面受到摩擦力减小，垂直的速度进行一个衰减后直接反向。

由于我们之前计算碰撞检测的时候算的都是具体的点，但是碰撞响应是要整个刚体进行响应，所以我们要想出办法，把点的速度改变 $\Delta v_i$ 表示为整个刚体的改变 $\Delta v$。我们可以先通过刚体的角速度 $\omega$、线速度 $v$ 来得到碰撞点碰撞前的速度 $v_i$，然后通过 $v_i$ 算出受到冲量 $j$ 后的新速度 $v_i^{new}$ ，虽然我们不知道这个 $j$ 具体是多少，但 $v_i^{new}$ 的计算不需要用到 $j$，最后用算出来的 $v_i、v_i^{new}$ 来算出冲量 $j$，将 $j$ 直接用于我们刚体的状态改变上（针对质心进行计算）

最后还有个几个小问题。

如果有多个点同时碰撞怎么办，为了减小计算量，只需要计算这些点的位置的平均值，把位置的平均值当做碰撞点做碰撞处理。

物体掉在地面上时在重力和弹力作用下会一直抖动（oscillation），解决办法是检测到物体即将进入静止状态时将反弹系数 $\mu$ 变小。

为什么 $Impulse$ 方法只更新了速度没有更新位置？这是一个非线性问题，之后讲约束的时候回再回来讨论。

2.2、Shape matching

参考论文：Muller et al. 2005. Meshless Deformations Based on Shape Matching. TOG (SIGGRAPH).

Shape Matching 方法先把刚体的每个顶点都当做可以自由移动的顶点，对每个点做单独的碰撞处理。然后通过最小化 MSE（均方误差）的方法求出把它们聚拢回刚体的形状，其中的关键是第二步。

$y_i$ 为聚拢前顶点位置，$c$ 为聚拢后的质心坐标，$R{r}_i$ 为该点以质心做参考点进行约束后相对质心的位置，类似于之前的公式：
$$x_i=x+R{r}_i$$

这里放宽要求让R可以为任意矩阵，记为A（推测这么做是为了方便让 MSE 取到极值）。这里 MSE 的表达式比较简单，直接求导数等于 0 的点就是极值点。其中的 ${\sum }_{i}A{r}_i=A{\sum }_i{r}_i=0$ ，因为质点系中各点相对于质心的向量之和为0。

证明：
$$

$$

$c$ 的表达式比较简单，直觉上也能理解，聚拢前的中心和聚拢后的中心点相同才可以让误差最小。

再对 $A$ 求导，这里求出来的 $A$ 是任意矩形，但我们希望最后刚体能维持原来的形状，还需要做一个 Polar Decomposition ，为何要做这步呢？

奇异值分解可以理解为先把图形旋转为我想要的角度 $(V^T)$，再进行缩放 $D$，最后再转回去 $(U)$，前两个旋转矩阵决定了缩放的角度和大小，与最后一次旋转比较独立。

那么我们也可以直接在本地进行一次缩放，然后再全局旋转为世界坐标系下对应的位置。事实上，由于这个是刚体，我们只想要他的旋转，并不希望他发生形变，所以把它在本地的形变部分刨除，直接用它相对于世界坐标系下的旋转，所以要进行Polar Decomposition。保留旋转矩阵代入，就能得到聚拢后的顶点位置。整体流程如下

具体的算法老师已经在unity中写好了，想了解更多可见 Franca, Leopoldo P. "An algorithm to compute the square root of a 3× 3 positive definite matrix." Computers & Mathematics with Applications 18.5 (1989): 459-466

优点：

容易实现
容易和点的模型结合在一起，包括衣服、软体、基于粒子的流体

缺点：

存在其他约束的情况下，很难保证所有的约束都满足
Shape Matching 适合摩擦力约束不是很重要的情况，比如衣服和扣子的模拟。