Games103 数学与刚体动力学

GAMES103

Introduction

如果想要在每一个时间点更新物体状态，这就是动画。所以动画的本质就是更新状态，比如内容、位置、姿势、外观等，而$\Delta t$可以和帧率不同步。物理模拟分成四个方向：

• 刚体 Rigidbody
• 衣服和头发 Cloth, Hair
• 软体 Soft Bodies
• 流体 Fluids

前三种我们一般用Mesh来模拟，因为形态比较固定。不过，Rigid Body的Fracture（破裂）可能用Particle，这样可以不需要更新Mesh；也有用Grid来模拟Cloth或Hair，来简化碰撞处理。

Fluids的展示方式都不太一样，所以研究方法也有所不同。Smoke常用Particle和Grid，Drops和Waves由于很多时候是整体，所以可以用Mesh做实时、Grid做离线。Splashes的实时模拟比较难，Particle和Grid比较多。

当然，这些东西也不一定是孤立的。图形学中有些方法是Hybrid的，比如MPM是Particle和Grid的结合，所以适合做雪或者粘制物体。也可能是Coupling，比如场景中又有Rigid Bodies也有Cloth。

Math

Vector

Vector : a geometry entity endowed with magnitude and direction.我们常常把Vector的空间记成$p \in R^3$。在paper中，一般用$\mathbf{p}$表示矢量，$p$表示标量。图形学中，做首席和右手系都比较常见。这主要因为屏幕空间，$z$的前后可以变化。

但其实Vector未必要有一个确切的形态，比如我们可以把
$$
\mathbf{p} = \begin{pmatrix}\mathbf{x_0} \\mathbf{x_1} \ \cdots \ \mathbf{x_n} \end{pmatrix}
$$
定义成一个$\mathbb{R}^{3n}$的向量，其中$\mathbf{x_i}\in \mathbb{R}^3$

矢量可以描述各种各样的东西。比如，
$$
\mathbf{p}(t)=\mathbf{p}+t\mathbf{v}
$$
这就是直线参数方程。也可以是
$$
\mathbf{p}(t)=\mathbf{p}+t(\mathbf{q}-\mathbf{p})
$$
如果我们转化一下，可以写成
$$
\mathbf{p}(t)=(1-t)\mathbf{p}+t\mathbf{q}
$$
t称作interpolant。

矢量的另一个概念是norm，最常见的是Euclidean norm:
$$
||\mathbf{p}|| _2 = (|p_x|^2+|p_y|^2+|p_z|^2)^{1/2}
$$
实际上，我们一般把
$$
||\mathbf{p}||=(|p_x|^p+|p_y|^p+|p_z|^p)^{1/p}
$$
记为p-norm。常用的还有1-norm和infinity norm。

内积是一种常见运算：
$$
\langle \mathbf{p},\mathbf{q} \rangle = \mathbf{p}\cdot \mathbf{q} = \mathbf{p}^T\mathbf{q} = ||\mathbf{p}||||\mathbf{q}||\cos \theta
$$
利用点乘可以定义点到线上的投影：

容易求得
$$
\mathbf{s} = \mathbf{o} + s\overline{\mathbf{v}}=\mathbf{o}+(\mathbf{q}-\mathbf{o})^T\overline{\mathbf{v}}\overline{\mathbf{v}}
$$
还可以定义平面：

那么平面上的点有
$$
s = (\mathbf{p}-\mathbf{c})^T\mathbf{n} = 0
$$
同时，$s>0$表明在正法线方向，$s<0$在负，可以用来做碰撞检测。

另一个例子是点与球的碰撞
$$
||\ve{p}(t)-\ve{c}||^2-r^2
$$
进一步展开，得到
$$
\vd{v}{v}t^2+2\vd{(\ve{p}-\ve{c})}{v}t+(\ve{p}-\ve{c})^2-r^2=0
$$
这个二次方程的解的分布就是碰撞。

叉积不太一样，得到的是一个新的向量：
$$
\ve{r}=\vc{p}{q} = \mat{p_yq_z-p_zq_y \ p_zq_x-p_xq_z\p_xq_y-p_yq_x}
$$
可以拿叉积计算三角形面积和法向

那么其法向是
$$
\ve{n} = \frac{\vc{x_{10}}{x_{20}}}{||\vc{x_{10}}{x_{20}}||} = 2A \nv{n}
$$
另一个问题是，给定一个三角形如何判断在三角形内外？这个也非常经典：如果
$$
(\vc{(x_0-p)}{(x_1-p)})^T\ve{n}>0
$$

$$
(\vc{(x_1-p)}{(x_2-p)})^T\ve{n}>0
$$

$$
(\vc{(x_2-p)}{(x_0-p)})^T\ve{n}>0
$$

那么在形内，否则形外。这个值是什么呢？

实际上，
$$
\frac{1}{2}\vc{(x_0-p)}{(x_1-p)}^T\ve{n} = \begin{cases}|A_2| & \text{inside} \ -|A_2| & \text{outside} \end{cases}
$$
如果把$A_1,A_2,A_3$赋予正负，那么$A_0+A_1+A_2=A$；如果让
$$
b_0 = \frac{A_0}{A},b_1 = \frac{A_1}{A},b_2 = \frac{A_2}{A}
$$
那么$b_0+b_1+b_2=1$，这也构成了一种三角形插值的方法。这个算法称为Barycentric Weights，当然一种加速是用Scanline Conversion.

进一步推广，四面体呢？四面体有四个顶点，它的体积是什么？是
$$
V = \frac{1}{6}\ve{x_{30}}\cdot \ve{x_{10}}\times \ve{x_{20}} = \frac{1}{6}\det{\ve{x_1}&\ve{x_2}&\ve{x_3}&\ve{x_0} \ 1 & 1 & 1 & 1}
$$

体积也可以有正有负，这个结果和$\ve x_1 \ve x_2 \ve x_0$得到的法线方向有关。有了体积的定义，我们就可以定义四个四面体，每个小四面体对应一个面，那么
$$
\cases{V_0 = \operatorname{Vol}(\ve x_3, \ve x_2, \ve x_1, \ve p) \ V_1 = \operatorname{Vol}(\ve x_2, \ve x_3, \ve x_0, \ve p)\ V_2 = \operatorname{Vol}(\ve x_1, \ve x_0, \ve x_3, \ve p) \ V_3 = \operatorname{Vol}(\ve x_0, \ve x_1, \ve x_2, \ve p)}
$$
$\ve p$在体积内，等价于$V_0 ,V_1,V_2,V_3>0$。同样，我们可以定义Barycentric Weight，即
$$
b_i = \frac{V_i}{V}
$$
一个典型的例子是判断点和三角形有没有发生碰撞。

首先，考虑什么时候碰撞：
$$
(\ve p(t)-\ve x_0) \cdot \ve x_{10} \times \ve x_{20} = 0
$$
那么
$$
t = \frac{(\ve p - \ve x_0)\cdot \ve x_{10}\times \ve x_{20}}{\ve v \cdot \ve x_{10} \times \ve x_{20}}
$$
然后我们只需要检查$\ve p(t)$是否在形内。

Matrices

矩阵可以认为是把实数用行列描述。
$$
\ve A = \mat{a_{00} & a_{01} & a_{02} \ & \cdots \ a_{20} & a_{21} & a_{22}} = \mat{\ve a_0 & \ve a_1 & \ve a_2} \in \mathbb{R}^{3\times 3}
$$
Orthogonality（正交）矩阵是一个比较重要的东西：
$$
\ve A = \mat{\ve a_0& \ve a_1& \ve a_2}
$$
其中
$$
\vd{a_i}{a_j} = \begin{cases}1 & i=j \ 0 & i\ne j \end{cases}
$$
对正交矩阵有
$$
\ve A^T = \ve A^{-1}
$$
我们可以用旋转去描述正交矩阵。假设
$$
\ve A \mat{\ve x \ \ve y \ \ve z } = \mat{\ve u \ \ve v \ \ve w}
$$
所以
$$
\ve A =\mat{\ve u \ \ve v \ \ve w}
$$
由于三条向量是局部坐标，所以$\ve A$是一个正交矩阵。

另一个比较重要的概念是singular value decomposition（奇异值分解）。什么是奇异值分解？对于矩阵$\ve A$，那么
$$
\ve A = \ve U \ve D \ve V^T
$$
其中$\ve D$是diagonal（对角）matrix，$\ve U,V$是orthogonal matrix.有一个很有趣的解释：

任何一种变换都可以拆解成三部分：

• 旋转到某个方向
• 沿坐标轴缩放
• 再一次旋转

因此，singular value decomposition的本质是scaling。

Eigenvalue Decomposition（特征值分解）是Singular value decomposition的一个特例。对symmetric matrix $\ve A$：
$$
\ve A = \ve U \ve D \ve U^T = \ve U \ve D \ve U^{-1}
$$
此时，
$$
\ve A \ve u_i = \ve U \ve D \ve U^T \ve u_i = \ve U \ve D \mat{\cdots \ 0 \ 1 \ 0 \ \cdots } = \ve U \mat{\cdots \ 0 \ d_i \ 0 \ \cdots } = d_i \ve u_i
$$
$d_i$记作eigen value，$\ve u_i$记作eigen vector

实际上，asymmetric矩阵也可以分解，只不过结果可能是虚数。

一种重要的矩阵是Symmetric Positive Definiteness(s.p.d. 正定矩阵)：

• $\ve v^T \ve A \ve v > 0$，对$\forall \ve v \ne 0$
• $\ve v^T \ve A \ve v \ge 0$，对$\forall \ve v \ne 0$

这种定义可以这么考虑：由于
$$
\ve v^T \operatorname{diag}(d_0,d_1,\cdots d_n) \ve v > 0
$$
而如果存在正定矩阵$\ve U$，
$$
\ve v^T (\ve U \ve D \ve U^T)\ve v = (\ve U^T \ve v)^T \ve D (\ve U ^T \ve v)
$$
而对于非零向量，$\ve U^T \ve v \ne \ve 0$，所以$\ve {UDU^T}$是正定矩阵。这就得到了上面的定义。

我们一般用一个性质来判断正定性：对角占优矩阵一定正定，即
$$
a_{ii} > \sum_{i \ne j} |a_{ij}|,,, \text{for }\forall i
$$
正定矩阵一定可逆。

【举例】已知$A$为s.p.d，证明
$$
\ve B = \mat{\ve A & -\ve A \ -\ve A& \ve A }
$$
正定。

证明：想要证明正定性，只需要证明对于$\forall \ve x, \ve y$有
$$
\mat{\ve x^T & \ve y^T}\ve B \mat{\ve x^T \ \ve y^T} = \ve x^T\ve A(\ve x - \ve y)-\ve y^T\ve A(\ve x - \ve y) = (\ve x - \ve y)^T\ve A(\ve x - \ve y)
$$
由于$\ve A$为s.p.d，所以上式一定成立。这就证明了$\ve B$的正定性。

我们还需要引入线性系统。一个线性问题可以写成这样的形式：
$$
\ve A \ve x = \ve b
$$
比较简单的想法是计算$\ve A^{-1}$，但由于逆很难计算，并且对于稀疏矩阵来说，其逆矩阵可能不稀疏，这可能也带来庞大的内存开销。这就引入了直接法和迭代法。

直接法的核心原理是LU分解，我们总是能把一个矩阵分解成下三角和上三角矩阵：
$$
\ve A = \ve L \ve U
$$
具体分解过程从略，但在求解之后可以通过计算和回带迅速求解。先求解$\ve L\ve y = \ve b$：
$$
\mat{l_{00} & & \ l_{10} & l_{11} \ \cdots & \cdots & \cdots} \mat{y_0\\cdots\y_n}=\mat{b_0\\cdots\b_n}
$$
求出$\ve y$后继续求解$\ve U \ve x = \ve y$
$$
\mat{\cdots & \cdots & \cdots \ & u_{n-1,n-1} & u_{n-1,n} \ & & u_{n,n} } \mat{x_0 \ \cdots \ x_n} = \mat{y_0 \ \cdots \ y_n}
$$
直接法有几个特征：

• $\ve A$是稀疏矩阵，但$\ve L, \ve U$不一定很稀疏，其稀疏性取决于行列之间的关系。所以我们往往先进行重排列。
• 它包含有分解和求解两部分，对于固定矩阵$A$分解只需要计算一次
• 并行性比较差，有Intel MKL PARDISO库可以进行并行

对于迭代法，迭代求解器有形式
$$
\ve x^{[k+1]} = \ve x^{[k]} + \alpha \ve M^{-1}(\ve b - \ve A \ve x ^{[k]})
$$
其中$\alpha$为松弛系数，$\ve M$为迭代矩阵，$\ve b - \ve A \ve x ^{[k]}$为残差。

如果我们进行变形，可以写成
$$
\ve b-\ve A \ve x^{[k+1]} = \ve b-\ve A \ve x^{[k]}-\alpha\ve A\ve M^{-1}(\ve b-\ve A \ve x^{[k]})
$$
提取公因式，得到
$$
\ve b - \ve A \ve x ^{[k+1]} = (\ve I - \alpha \ve A \ve M^{-1})^{k+1}(\ve b-\ve A \ve x^{[0]})
$$
当且仅当$\rho(\ve I - \alpha \ve A \ve M^{-1})<1$，迭代法收敛。

对于$\ve M$来说，有两种选取方式：

• $\ve M = \operatorname{diag}(\ve A)$，则为Jacobi迭代法
• $\ve M = \operatorname{lower}(\ve A)$，则为G-S迭代法

一般来说，迭代法实现简单，跑起来很快，并且容易并行。但同时，它可能需要收敛条件，并且求精确解比较慢。

Tensor Calculus

首先引入一阶导：
$$
\dd f = \pd{f}{x}\dd x+\pd{f}{y}\dd y + \pd{f}{z}\dd z
$$
我们换一种写法，
$$
\pd{f}{\ve x} = \mat{\pd{f}{x}&\pd{f}{y}&\pd{f}{z}}
$$
也可以写成梯度的形式
$$
\grad f(\ve x) =\mat{\pd{f}{x}&\pd{f}{y}&\pd{f}{z}}^T
$$
接下来，我们进行扩展：如果$\ve f(\ve x) = \mat{f(\ve x)& g(\ve x) & h(\ve x)}^T$，则引出三个比较重要的概念：

（1）Jacobian矩阵
$$
J(\ve x) = \pd{\ve f}{\ve x} = \mat{\pd{f}{x}&\pd{f}{y}&\pd{f}{z} \\pd{g}{x}&\pd{g}{y}&\pd{g}{z} \\pd{h}{x}&\pd{h}{y}&\pd{h}{z} \ }
$$
（2）Divergence
$$
\grad \cdot \ve f = \pd{f}x + \pd{g}y+\pd{h}z
$$
（3）Curl
$$
\grad \times \ve f = \mat{\pd{h}y - \pd{g}z \ \pd{f}z-\pd{h}x\\pd{g}x-\pd{f}y}
$$
进一步的，若$f(\ve x)\in \mathbb{R}$，有

（1）Hessian矩阵
$$
\ve H = \ve J(\grad f(\ve x)) = \mat{
\pd{^2f}{x^2}&\pd{^2f}{x\partial y}&\pd{^2f}{x\partial z}\
\pd{^2f}{x\partial y}&\pd{^2f}{y^2}&\pd{^2f}{y\partial z}\
\pd{^2f}{x\partial z}&\pd{^2f}{y\partial z}&\pd{^2f}{z^2}}
$$
（2）Laplacian算子
$$
\grad \cdot \grad f(\ve x) = \grad^2 f(\ve x) = \Delta f(\ve x) = \pd{^2f}{x^2}+\pd{^2f}{y^2}+\pd{^2f}{z^2}
$$
Hessian矩阵常用于Taylor展开中，对$f(\ve x)$有
$$
f(\ve x) = f(\ve x_0) + \grad f(\ve x_0)\cdot (\ve x - \ve x_0) + \frac{1}{2}(\ve x - \ve x_0)^T \ve H (\ve x - \ve x_0)+\cdots
$$
当然，这里可以看到若$\ve H$是s.p.d，那么$f(\ve x)$的二阶导是正数。

【举例】求解
$$
\pd{||\ve x||}{\ve x}
$$
解：
$$
\pd{||\ve x||}{\ve x} = \pd{(\ve x^T \ve x)^{1/2}}{\ve x} = \frac{1}{2(\ve x^T\ve x)^{1/2}} \pd{(x^2+y^2+z^2)}{\ve x} = \frac{\ve x^T}{||\ve x||}
$$
所以原方向是向量长度变换最短的方向。

考虑有一根弹簧，原长$L$，能量$E(\ve x) = \frac{k}{2}(||\ve x||-L)^2$，我们知道
$$
\ve f(\ve x) = -\grad E(\ve x) = -k(||\ve x|| - L) \left(\pd{||x||}{\ve x}\right)^T
$$
所以
$$
\ve f(\ve x) =-k(||\ve x|| - L) \ve x^0
$$
这就是胡克定律。我们进一步，求Hessian矩阵：
$$
\ve H(\ve x) = -\pd{\ve f(\ve x)}{\ve x} = k\frac{\ve x \ve x ^T}{||\ve x||^2}+k(||\ve x||-L)\frac{1}{||\ve x||} - k(||\ve x||-L)\frac{\ve x}{||\ve x||^2}\frac{\ve x^T}{||\ve x||}
$$
整理一下得到一个很奇怪的玩意
$$
\ve H(\ve x) = k\frac{\ve x \ve x ^T}{||\ve x||^2}+k\left(1-\frac{L}{||\ve x||}\right)\left(\ve I - \frac{\ve x \ve x ^T}{||\ve x||^2}\right)
$$
这就是Hessian矩阵，我们称为Tangent Stiffness。

现在我们再进行一下推导，如果弹簧有$\ve x_0, \ve x_1$两个点固定，$\ve x_{01}=\ve x_1 - \ve x_0$，那么$E(\ve x) = \frac{k}{2}(||\ve x_{01}||-L)^2$。此时，
$$
\ve f(\ve x) = -\grad E(\ve x) = \mat{-\grad_0 E(\ve x) \ -\grad_1 E(\ve x)} = \mat{\ve f_e \ -\ve f_e}
$$

$$
\ve H(\ve x) = \mat{\ve H_e & -\ve H_e \ -\ve H_e & \ve H_e}
$$

其中$\ve f_e, \ve H_e$是之前所求出的结果。

Rigid Body Dynamics

物体的运动

没有形变的物体叫做刚体。几乎所有的引擎都提供了刚体模拟系统，它所做的事情是更新描述物体状态的变量。考虑时刻$t_0t_1\cdots t_n$，对应状态量$\ve s_0 \ve s_1 \cdots \ve s_n$，状态之间的变换是由Update Function所完成。

假设兔子的中心是原点，经过某一个旋转，使用矩阵$\ve R$来描述。那么$\ve R \ve r_i$就是变化后的位置。再平移$\ve x$，兔子就到了$\ve x + \ve R \ve r_i$。

接下来考虑物体的平移。物体从$\ve x^{[0]}$移到$\ve x^{[n]}$，那么有
$$
\begin{cases}
\ve v(t^{[1]}) = \ve v(t^{[0]}) + M^{-1} \int_{t^{[0]}}^{t^{[1]}} \ve f(\ve x(t), \ve v(t), t) \dd t \
\ve x(t^{[1]}) = \ve x(t^{[0]}) + \int_{t^{[0]}}^{t^{[1]}} \ve v(t)\dd t
\end{cases}
$$
如何求解积分？一阶精度的显式积分是
$$
\int_{t_0}^{t_1} \ve v(t) \dd t \approx \Delta t \ve v(t_0)
$$
隐式积分是
$$
\int_{t_0}^{t_1} \ve v(t)\dd t \approx \Delta t \ve v(t_1)
$$
这两种方法都是一阶的。而我们只需要把点换到中间，它就变成了二阶：
$$
\int_{t_0}^{t_1} \ve v(t)\dd t \approx \Delta t \ve v(t_{0.5})
$$
证明可以用Taylor展开，这里从略。对于上面的问题，我们一般进行一下混合：
$$
\begin{cases}
\ve v(t^{[1]}) = \ve v(t^{[0]}) + \Delta tM^{-1} \ve f^{[0]} &(\text{Explicit}) \
\ve x(t^{[1]}) = \ve x(t^{[0]}) + \Delta t \ve v^{[1]} & \text{(Implicit)}
\end{cases}
$$
这个方法也称为semi-implicit或leap-frog。

不妨假设物体受到重力和空气阻力影响，对于重力
$$
\ve f = M \ve g
$$
空气阻力
$$
\ve f = -\sigma \ve v
$$
那么从$\ve v^{[0]},\ve x^{[0]}$更新的过程是
$$
\begin{cases}
\ve f_i^{[0]} = \operatorname{force}(\ve x_i^{[0]}, \ve v_i^{[0]}) \
\ve f^{[0]} = \sum \ve f_i^{[0]} \
\ve v^{[1]} = \ve v^{[0]}+\Delta t M^{-1} \ve f^{[0]} \
\ve x^{[1]} = \ve x^{[0]} + \Delta t \ve v^{[1]}
\end{cases}
$$
再考虑旋转，由于在模拟的时候使用矩阵效率会很低，而且矩阵的9个自由度实际上只需要3个，这就导致了冗余性。同时它不直观，而且对时间微分比较难求，所以矩阵方式不是很好的选择。

欧拉角比较直观，旋转由$z-x-y$的顺序构成，但是做模拟的时候存在gimbal lock，时间导数也不直观。所以我们引入Quaternion，把虚数空间和变换联系起来，定义所有的四则运算。

那么可以定义：对$\ve q = [s \ \ve v]$，有
$$
a\ve q = [as \ \ a\ve v]
$$

$$
\ve q_1 \pm \ve q_2 = [s_1 \pm s_2 \ \ve v_1 \pm \ve v_2]
$$

$$
\ve q_1 \times \ve q_2 = [s_1s_2-\ve v_1 \cdot \ve v_2\ \ s_1\ve v_2 + s_2\ve v_1 + \ve v_1 \times \ve v_2]
$$

$$
||\ve q|| = \sqrt{s^2 + \ve v \cdot \ve v}
$$

四元数可以用来表示旋转。假设旋转轴是$\ve v$，角度是$\theta$，约定$||\ve q||=1$，那么旋转就是
$$
\ve q = [\cos \frac{\theta}{2}\ \ \ve v]
$$
则$||\ve v||^2 = \sin^2 \frac{\theta}{2}$。从四元数转换到矩阵也有一个公式，比较麻烦，不敲了。

这样，我们就可以根据Quaternion来维护旋转了。但我们希望知道它的时间的导数，也就是角速度；如果用3D向量$\ve \omega$表示旋转，方向是旋转轴，模长是速度，这个定义和四元数是比较契合的。

除此之外，我们还需要力矩（Torque）和转动惯量（Inertia）。
$$
\ve{\tau}_i = (\ve R \ve r_i) \times \ve f_i, \tau = \sum \tau_i
$$
力矩是造成旋转趋势的因素，由于我们需要同时考虑旋转轴和力，所以我们用距离和力的叉积来表示旋转的大小，也就是力。$\ve R$是当前的状态参量，表示物体此时本身具有的旋转变换。

对于刚体来说，质量是一个矩阵
$$
\ve {I_{ref}} = \sum m_i(\ve r_i^T \ve r_i \ve E - \ve r_i \ve r_i^T), \ve I = \ve R \ve {I_{ref}} \ve R^T
$$
考虑下面的两只兔子：

给兔子带两只重砝码，他们二者的$\ve\tau$是相同的，但旋转并不相同，所以我们需要考虑到旋转方式的影响。这里要注意，如果兔子本身发生变换$\ve{R}$，则
$$
\ve I = \sum m_i(\ve r_i^T \ve R^T \ve R \ve r_i \ve E - \ve R \ve r_i \ve r_i^T \ve R^T)
$$
考虑到$\ve R^T \ve R = \ve E$，可以把单位矩阵拿到外边，有
$$
\ve I = \sum m_i(\ve R\ve r_i^T \ve r_i \ve E\ve R^T - \ve R \ve r_i \ve r_i^T \ve R^T)
$$
提取公因式，
$$
\ve I = \ve R \ve {I_{ref}} \ve R^T
$$
这样大大简化了计算。因此对于旋转来说
$$
\ve{\omega^{[1]}} = \omega^{[0]} + \Delta t(\ve I^{[0]})^{-1}\tau^{[0]}
$$

$$
\ve q^{[1]} = \ve q^{[0]} +\left[0 \ \ \frac{\Delta t}{2}\omega^{[1]}\right] \times \ve q^{[0]}
$$

我们就定义了全部的状态。

二者之间其实是有相互对应关系的。$\ve f$和$\ve \tau$是相互对应的，$\ve f^{[0]}$和$\boldsymbol\tau$是对应的，$\ve I$和$M$是对应的。

物体的碰撞

Particle Collision Detection and Response

定义signed distance function$\phi(\ve x)$是一个有符号的距离函数。下图表示了它的定义。

因此，一些简单的物体就写出距离函数。

那么如何组合物体呢？我们就需要用到布尔操作。比如下面的物体：

它有两个面和一个圆。当
$$
\phi_0(\ve x) < 0 \and \phi_1(\ve x)<0 \and \phi_2(\ve x) < 0
$$
那么物体在内侧。此时，
$$
\phi(\ve x) = \max (\phi_1(\ve x), \phi_2(\ve x), \phi_3(\ve x))
$$

对于并集，如果
$$
\phi_0(\ve x) < 0 \or \phi_1(\ve x) < 0
$$
那么物体在内测。此时，对于物体内侧，只能认为
$$
\phi(\ve x) \approx \min(\phi_0(\ve x), \phi_1(\ve x))
$$
这个近似对于物体边缘比较接近。对于物体外侧，则
$$
\phi(\ve x) = \min(\phi_0(\ve x), \phi_1(\ve x))
$$
对于碰撞，我们一般有两种模型。一种叫做penalty，也就是在下一帧直接应用力。这个力的方向是$\ve N = \grad \phi(\ve x)$，大小是$-\phi(\ve x)$，即
$$
\ve f = -k \phi(\ve x) \ve N
$$
当然，由于我们检测的时候可能发生碰撞，所以一般设置一个域值$\varepsilon$，即对于$\phi(\ve x) < \varepsilon$，有
$$
\ve f = k(\varepsilon - \phi(\ve x)) \ve N
$$
$k$是决定物体大小的因素，$k$太小物体会一直运动，$k$太大物体会直接飞出去也就是overshooting问题。一种想法是把$k$处理成关于距离的函数，那么
$$
\ve f = \rho \frac{1}{\phi(\ve x)} \ve N
$$
称作log-barrier function。但这个东西仍然避免不了overshooting，而且在物体穿透之后，会越陷越深，反而会往里走。所以这种方式假设$\phi(\ve x)<0$不会发生，因此需要非常小的步长。

log-barrier有一些改进，比如和buffer结合。这个方法的另一个问题是不容易和摩擦结合。

另一种叫做Impulse方法，Penalty的力是滞后一帧的，而Impulse则是在当前帧直接更新速度和位置。当collision发生的时候，我们直接让
$$
\ve {x^{new}} = \ve x + |\phi(\ve x)|\ve N = \ve x - \phi(\ve x) \grad \phi(\ve x)
$$
那么如何考虑速度和摩擦？我们希望考虑速度和力的方向。如果$\vd{\ve v}{\ve N}<0$，那么就需要进行更新。首先进行正交分解，
$$
\ve {v_N} = (\ve v \cdot \ve N) \ve N, \ve{v_T} = \ve v - \ve{v_N}
$$
此时分别进行更新，
$$
\ve{v_N^{new}} = -\mu_N \ve{v_N}, \ve{v_T^{new}} = a\ve{v_T}
$$
然后在加起来。根据Coulomb’s law，
$$
||\ve {v_T^{new}} - \ve {v_T}|| \le \mu_T ||\ve {v_N^{new}} - \ve {v_N}||
$$
进行一下变换
$$
(1-a)||\ve{v_T}|| \le \mu_T(1+\mu_N)||\ve{v_N}||
$$
所以
$$
a = \max(1-\mu_T(1+\mu_N)||\ve {v_N}||/||\ve{v_T}||, 0)
$$
二者分别对应了动摩擦和静摩擦。

Rigidbody Collision Detection and Response

对于一个物体上的点，有
$$
\ve x_i = \ve x + \ve {Rr_i}
$$
然后讨论$\phi(x)$大小即可。对于碰撞响应，我们知道原本
$$
\ve x_i = \ve x + \ve R \ve r_i, \ve v_i = \ve v + \ve{\omega}\times \ve R\ve r_i
$$
但问题在于，我们并不存在具体点的状态参量，而应该修改整个物体的。简单的，我们只去修改物体本身的速度参量。假设有冲量$\ve j$，那么
$$
\ve {v^{new}} = \ve v + \frac{1}{M}\ve j, \ve{\omega^{new}} = \ve \omega + \ve I^{-1}(\ve R \ve r_i \times \ve j)
$$
那么更新到每一个点，
$$
\ve {v_i^{new}} = \ve {v_{new}} + \ve{\omega^{new}\times Rr_i}
$$
由于$\ve v_i = \ve v + \ve \omega \times \ve {Rr_i}$，所以
$$
\ve {v_i^{new}} = \ve {v}_i + \frac{1}{M}\ve j - (\ve R\ve r_i)\times (\ve I^{-1}(\ve R\ve r_i\times \ve j))
$$
为了简化掉叉积，我们引入Cross Matrix：$\ve r \times \ve q$相当于
$$
\ve r \times \ve q =\mat{0 & -r_z & r_y \ r_z & 0 & -r_x \ -r_y & r_x & 0}\mat{q_x \ q_y \ q_z} = \ve r^\ve q
$$
所以
$$
\ve {v_i^{new}} = \ve {v}_i + \frac{1}{M}\ve j - (\ve R\ve r_i)^ (\ve I^{-1}((\ve R\ve r_i)^* \ve j))
$$
进一步简化有
$$
\ve {v_i^{new}} = \ve v_i + \ve K\ve j, \ve K = \frac{1}{M}\ve E - (\ve {Rr}_i)^* \ve I^{-1} (\ve{Rr}_i)^*
$$
因此，只要根据Impulse的方法计算出来$\ve v, \ve {v^{new}}$的结果，由于$\ve K$已知，就可以反求出来$\ve j$。

所以我们构造出如下算法：

// detection
For each vertex x_i = x + Rr_i
if phi(x_i) < 0
	done
v_i = v + omega.cross(R*r_i)
if v_i.dot(N) < 0
	done
// impulse
calculate v_Ni, v_Ti, a
calculate v_Ninew, v_Tinew, v_inew
// calculate impulse j   
K = E_3 / M - M_Rri * inv(I) * M_Rri
j = inv(K) * (v_inew - vi)
// calculate v, omega
v = v + 1/M * j
omega = omega + inv(I)*((R*r_i).cross(j))

这里有一些细节：

• 对于多个点的碰撞，采用位置平均值
• 如果停不下来，可能是因为重力，导致物体表面的抖动（oscillation）。这里可以在物体附近衰减$\mu_N$。

当多点接触，彼此之间有影响，这里不多展开。

Shape Matching

Shape Matching的思路就是，对于一个物体，假设每个点都有自己的速度，对多个质点分别进行模拟。接下来，再让它变回刚体，满足约束。

假设对于一系列的物体，我们希望找到$\ve c, \ve R$。所以我们建立了一个优化问题：寻找质心点的位置，使得收的距离越来越近。即
$$
{\ve c, \ve R } = \operatorname{argmin} \sum_i \frac{1}{2}||\ve c + \ve R\ve r_i -\ve y_i||^2
$$
不妨假设$\ve R$为任意矩阵，那么上述问题转换为
$$
{\ve c, \ve A } = \operatorname{argmin} \sum_i \frac{1}{2}||\ve c + \ve A\ve r_i -\ve y_i||^2
$$
求导，
$$
\pd{f}{\ve c} = \sum_i \ve c + \ve A\ve r_i - \ve y_i=\sum_i \ve c - \ve y_i = 0
$$
这里是基于$\sum \ve A\ve r_i = 0$。因此有
$$
c = \frac{1}{N}\sum_i \ve y_i
$$
再对$\ve A$求导，有
$$
\pd{f}{\ve A} = \sum_i (\ve c + \ve A\ve r_i - \ve y_i)\ve r_i^T = 0
$$
得到
$$
\ve A = \left(\sum_i(\ve y_i - \ve c)\ve r_i^T\right) \left(\sum_i \ve r_i \ve r_i^T\right)^{-1}
$$
由于任意矩阵一定可以分解成旋转和缩放变换，即
$$
\ve A = \ve U \ve D \ve V^T
$$
考虑到我们可以在做最终的旋转$\ve U$之前，先转$\ve V$转回去，再施以旋转$\ve{UV}^T$，即
$$
\ve A = (\ve{UV}^T)(\ve{VDV}^T)
$$

可以看到，前几步实际上只是在本地对某个轴进行了形变，即$\ve S = \ve {VDV}^T$，所以我们知道可以把任意矩阵分解成旋转与缩放矩阵。我们称之为Polar Decomposition，即
$$
\ve A = \ve R \ve S
$$
取$\ve R$为状态参量，则得到的${\ve c, \ve R}$。

因此，Shape Matching的伪代码如下：

// update each vertex
For every vertex
	fi = Force(xi, vi)
	vi = vi + dt / mi * fi
	yi = xi + dt * vi
// estimate state
c = 1 / N * sum(yi)
A = (sum(yi-c)*trans(ri))*inv(sum(ri * trans(ri)))
R = polar(A)
// update each vertex with state
vi = (c + R * ri - xi) / dt
xi = c + R * ri

这个方法没有物理参量，由于基于点，很容易和其它东西结合一起。但是这个问题比较难以保证约束，比如摩擦等，可能需要迭代来反复约束。当碰撞、摩擦接触不频繁或是要求不精确，可以考虑这种方法。