GAMES103-基于物理的计算机动画入门（1~4，lab1）

前言

本人的GAMES103笔记，以下分节根据知识点来分。同时包含了个人的一点补充。

1. Introduction

课程主页：http://games-cn.org/games103/
B站链接：https://www.bilibili.com/video/BV12Q4y1S73g?from=search&seid=2084315763141666934&spm_id_from=333.337.0.0
课程作业：http://games103.games-cn.org/

课程安排：

理想情况下管线：几何Geometry =》动画Animation =》渲染Rendering

Geometry 的三种表达方式：

• Mesh：分为structured mesh（有一定规律的，因此有些几何算法就可以利用这些规律做文章）和unstructured mesh（没有什么规律的，通用性就更强）
• Point Cloud（点云）：直接给点云是没法用硬件做渲染的，因此有问题比如如何去构造mesh，如何重新采样，如何找近邻等等。
• Volumetric Grid：医学、科学实验很多原始数据都是，比如CT、大气的流体模拟。

Rendering：

• Non-Photorealistic 非真实感渲染
• Photorealistic 真实感渲染
  两种类型：基于光线追踪的和基于传统渲染管线的。

基于物理的计算机动画（physics-based animation）：
本质上是更新每一帧的状态。

两个时刻之间的东西叫做步长。步长可以和帧率一致也可以不一致，完全可以低于帧率。比如渲染一帧用更多的时间步长去模拟，取决于引擎的设计。

四大类：Rigid Bodies、Cloth and Hair、Soft Bodies、Fluids

同样由于 Geometry 有三种表达方式，所以这里的模拟也有三种方式：
mesh在处理形态比较固定的东西上是有很大的好处的，所以很多都是mesh做的。比如刚体、衣服、弹性体都是以mesh为主的。当然也有人研究用别的方式去模拟。


以上仍然有一些未涉及的，比如说：

• Coupling结合、连接，一个场景既有流体又有刚体又有弹性体等，怎样去模拟它们之间的交互。比如水落在衣服上之类的。

• hybrid混合，比如雪、各种粘滞的物体，最大的特点就是既有particle的属性也有grid的属性，所以说用hybrid的方法是最常见最通用的。

这节课涉及到的内容：

任何有形变的都可以视为一个软体。无形变的视为刚体。真正的形变有三个阶段：弹性形变、塑性形变、破碎的形变。

老师的经验面板：

2. Math Background: Vector, Matrix and Tensor Calculus

该节主要为复习和物理模拟相关的数学知识。我之前的笔记：https://zhuanlan.zhihu.com/p/413224014

大多数都非常基础，这里记录几个稍微特别的：

叉乘应用：

计算四面体体积，最底下是 4 * 4 行列式：

四面体重心坐标（通过刚刚的公式得到四个子四面体的体积）：

同样还能用于找交点，当然只是一种解释方式罢了（这里解释为交点和三角形组成的四面体的体积应该为0，然后用刚刚说的四面体体积公式），其实很显然：

例子6就是用叉乘判断点在不在三角形里头，常规操作：

矩阵：

首先是奇异值分解（SVD）：

我们熟知旋转矩阵是正交阵，转置即为逆，这里老师给出了一个很有意思的见解：任何一个线性变换，都可以分解为 旋转=》缩放=》再旋转，之所以搞了两个旋转是因为，缩放只能沿着坐标轴缩放，所以要先旋转来调整位姿，缩放完之后再旋转，其核心是缩放。

图中为二维的例子，但实际上三维也是一样的道理。

为什么说这个解读很有趣呢，因为我从没见过，老师解释的非常好。

正统SVD的描述是：

假设M是一个m×n阶矩阵，其中的元素全部属于域 K，也就是实数域或复数域。如此则存在一个分解使得

其中U是m×m阶酉矩阵；Σ是半正定m×n阶对角矩阵；而V*，即V的共轭转置，是n×n阶酉矩阵。这样的分解就称作M的奇异值分解。Σ对角线上的元素Σi，其中Σi即为M的奇异值。

实际上M未必是方阵，而U和V为酉矩阵（自然便是方阵），并且还可以是复数域（所以扩充出了酉矩阵这样的概念），只不过搬到图形学我们都注重三维（或二维）的实数域了。

特征值分解：

我接着进行补充，这里老师只考虑对称矩阵，实际上上面分解是因为：任何一个实对称矩阵都相似于一个对角阵。因此可以有上面的分解。至于相似矩阵就自行看百度百科了，我放个链接：https://baike.baidu.com/item/%E7%9B%B8%E4%BC%BC%E7%9F%A9%E9%98%B5/10369874?fr=aladdin

对称正定/半正定：

百度百科：https://baike.baidu.com/item/%E6%AD%A3%E5%AE%9A%E7%9F%A9%E9%98%B5


上图这里给了一个结论：如果一个矩阵是对角占优的，那么它就是正定的。

并且一个spd矩阵（对称正定矩阵）一定是可逆的。

接着补充：

一般我们学数学的时候，判断一个矩阵是不是正定的其实常通过看顺序主子式。一个n阶实对称矩阵是正定的当且仅当它的一切顺序主子式均为正。有了这一条，那么自然它就是可逆的（本身的行列式就是最大的主子式，为正，那么自然可逆）。

求解线性问题：

两种方法：直接法和迭代法。

• 直接法：
  主要是通过LU分解。这样我们就可以分两步去求解。
  
  特点：
  如果A是稀疏的，L和U也未必稀疏。我们希望L和U稀疏来减少存储空间，于是可以通过调整顺序（permutation）来做到。
  计算实际上是由两部分构成：LU分解 和 之后的求解。目前而言两部分可以分开来计算，好处是如果矩阵A是固定的，那么可以让分解只做一次，而让求解做多次。

• 迭代法：
  
  
  迭代法实现起来较方便，且不要求精确解的情况下比较有效。并且相对容易并行。

缺点是有一些条件，比如上上张图的最底下，要求spectral radius谱半径要小于1才能收敛，同时要求精确解的话会很慢。

Tensor Calculus：

梯度的本意是一个向量（矢量），表示某一函数在该点处的方向导数沿着该方向取得最大值，即函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大（为该梯度的模）。

散度（divergence）：https://baike.baidu.com/item/%E6%95%A3%E5%BA%A6/8281793?fr=aladdin
旋度（Curl）：https://baike.baidu.com/item/%E6%97%8B%E5%BA%A6/8106439?fr=aladdin

Hessian对称阵，Laplacian就是其对角。

泰勒展开：

粗体的x表示向量。下面的公式即对向量形式的函数做泰勒展开。

可以看到，如果这个H是正定的，由正定的定义，我们知道二阶项就是大于0的，这会带来一些小特性可以优化，不展开了。

一个小测试（把x的转置乘以x看作一个整体，然后求导运用链式法则）：

通过这个小结论，我们进一步可以算出一个弹簧的公式：

如果是两个点的弹簧，则如下：

3. Rigid Body Dynamics

刚体模拟包括其他模拟，本质上就是去更新描述物体状态的一个变量。

因此引出的问题就是如何定义这个状态。

对于一个刚体，由于其无形变，因此状态包含两部分：平移和旋转。

位置很简单，就是一个三维向量。而知道位置还不够，因为实际情况我们还需要知道位置和时间的关系。因此我们还需要速度（位置对时间的导数）。

因为一个力（f函数）本身可能很复杂，可能是时间、速度、位置等的函数，同时由物理学我们知道，物体所受合外力的冲量就是该物体的动量变化量，冲量即力对时间做积分，动量即质量乘以速度，因此可得如下公式：

这里的M就是质量了，就是刚刚说的冲量等于动量变化量的等式，把动量中的质量除过去了而已。然后同样有位置，速度对时间的积分，都是基操了。

随后老师给出了一些积分手段：

分为了三种：显式、隐式和Mid-point，于是进一步公式积分可以写为：

这种结合的方法也叫 semi-implicit 半隐式，也有文献叫做 leapfrog method（像两只青蛙相互跳跃）。

我的补充：
这些都是数值分析的一些常规，进一步还可以了解：牛顿-柯斯特公式、辛普森公式、柯斯特公式、欧拉法（欧拉公式）、龙格-库塔法（R-K法）

力有重力、阻力，如下：

简单来做直接把速度搞成一个衰减的形式来模拟阻力，比较推荐。但不排除有paper用更精确的计算。

有了上面的基础就可以进行状态模拟了，如下：

物理学中质量的英文为mass，我们简单来做就直接把质量和时间步长自己定义（作业里就手动改了尝试一下）。但是可能游戏里 Δt 要和帧率相关。

接下来再谈旋转，我们熟知可以用矩阵来表示，但他也有缺陷：有9个变量却只有3个自由度（有很大的冗余性，即不是所有的矩阵都是旋转矩阵）、很不直观（不能一眼看出如何旋转）、时间微分不容易求（即不容易求旋转速度角速度）。

而像unity的rotation中的x、y、z就是用欧拉角去做的，unity顺序为z-x-y。不过拿欧拉角来模拟也有问题，其一是万向节死锁（gimbal lock），其二是也不好描述时间导数（毕竟有三个旋转，还彼此相关）。推荐我的笔记：https://zhuanlan.zhihu.com/p/413224014

最后便引出了四元数：



具体到unity中，面板上是欧拉角，而内置有四元数方法，两者可以通过赋值之类的操作简单转换：

于是乎，我们就可以用四元数来模拟旋转了，但是还没完，还有时间的导数，即角速度，我们用 w 来表示，其方向就是速度旋转轴：

有了旋转的表示，我们还差力和质量，在旋转中叫做力矩（torque）和转动惯量（inertia）：

关于角速度的更新：

为了避免浮点精度溢出，最好对求得的 q 再normalize一下：

于是这四个变量（速度、位置、角速度、四元数）分别描述了刚体的状态，而质量和inertia都是事先提供的，力和力矩是需要中间计算的。

于是我们的模拟引擎就应该是这样：

模拟中我们只需要每次算出来新的状态：s={v,x,ω,q}，之后再覆盖即可。对上图左边，首先算出每个力，然后求和算出总的力，然后根据牛顿定律得到对速度的更新，然后速度乘以时间步长得到对位置的更新；对于上图右边，就是先有旋转矩阵，然后计算出力矩，接着需要根据reference state的inertia以及当前的旋转矩阵去算出current state的inertia，随后就可以更新角速度和四元数了。

对于扭矩（其实就是一种特殊的力矩），老师也解释了一下：

ri是固定的，而每次的旋转矩阵 R 则是变化的（可以用四元数去描述当前倾斜的状态）。可以从图中直观的看到：扭矩就是那俩向量的叉乘。

对于转动惯量，老师也给了一张很直观的图：

如上图，显然右侧旋转应当更快，因为右侧的质量都集中在轴上，所以产生的旋转更强。即说明了同一个物体产生的旋转不仅和 torque 有关，还和物体本身的状态有关，因此需要用一个矩阵去计算转动惯量：

4. Rigid Body Contacts（Lab 1 ）

这节课讲涉及如下内容：

Particle Collision Detection and Response

首先介绍了带符号距离场sdf，可以参考我的笔记：https://zhuanlan.zhihu.com/p/398656596，此处老师定义的在外侧outside距离函数的值是大于0的。

这里提一下sdf的组合：

那么要是点x在内侧，也就是需要在三条边的距离函数都要在inside侧，即都小于0才行。点x的sdf值就是三者的最大值（因为都是负数，点x的sdf值其实就是到最近的边的距离的负数，图中为Phi2(x)）

对于合并Union的情况，如果两个球合并起来，想要做相交检测，怎么判断这个点在内侧还是外侧呢？显然对于两个球分别的sdf函数，只需要一个sdf值小于0就在内侧了；但是要注意它的值有可能和两个sdf值都毫无关系，比如如图两个sdf值是上图两个红色箭头的模值，而x离边界最近的点却有可能是上图中的红点或是其他的点。因此只能说接近外表面的时候有这样一个近似。当然要是是outside就是两个sdf值的最小值了。

其实也可以认为这就是两个物体，然后分别做碰撞检测，是一样的。

对于一个简单的penalty的模型，检测是否碰撞一样就看sdf值是否小于0（点是否进入了物体内侧），如果发生碰撞（点在物体内侧）一个简单的处理就是把点推出来：

而如上图，定义力一个很简单的方法就是像弹簧一样，进去越深认为力f值越大。如上，方向就是这个sdf函数的gradience，大小就是这个sdf函数值的相反数（因为在内侧sdf值为负），k为这个 penalty 强度。

那么这种方法的弊端就是，要穿透进物体了才施加力，这就势必导致明显的artifact，一个改善措施如下：

即弄一个epsilon，提早进行检测并处理。即只要靠近物体表面就加力把物体推出来，搞了一个小的缓冲空间。

然而这样做还是有点问题，即这个处理和 k 有关，如果 k 太小，那么要是物体穿透的比较剧烈就会有问题。而要是 k 太大，那么模拟的时候则可能导致产生的力过大，简单碰撞都直接被撞飞了，即 overshooting 现象。

避免的思路就是不要用常数 k ，把 k 弄成和距离有关的函数，得到如下的 Log-Barrier 的处理：

所谓 Log-Barrier 其实是对 log 求导就变成了倒数的方式。一般这个强度 barrier strength 是固定的。

但是这样仍然有蛮多问题，一是仍然避免不了 overshooting，二是如果穿透了但是这个力的方向却是相反的，会越陷越深。因此这个方法实际上有一个假设即穿透永远不会发生。因此要避免穿透就只能采用小步长的方法：每次往前挪一点，以确保每次都不发生穿透。因此会带来很多额外的成本。

Penalty Methods 总结：

好处是易实现：加了一个力，然后就可以按正常的模拟方式去计算。同时这里我们可以看出，施加一个力它并不是马上就起作用的，因为力的模拟将在下一个时刻才进行。是有一个滞后性的。

接着介绍 Impulse 的方法：

Impulse，本意就是 脉冲、冲击，物理中意为冲量。即我们不想像 Penalty 一样有滞后性，想要马上起作用，即碰撞检测到之后马上更新速度位置。

想法很简单，把 x 沿着那个 Normal 的方向移动那么多距离就好了：

那么对于速度呢？能否考虑摩擦量呢？

首先要判断速度方向大小，先做个简单的分解（切方向和法方向），然后去进行更新：

这个a和摩擦力有关，我们永远希望越小越好，把切方向的速度衰减掉，但同时a须满足一个库仑定律，并且避免反转速度需要和 0 取最小值。因此计算如上图右下角所示。如果取 0 为静摩擦（说明强到让其停下来，为静摩擦），取左边的值则为动摩擦，实际上反应出摩擦力的两种形式。

Rigid Body Collision Detection and Response

回忆，local space 和 world space（reference 和 current）：

刚刚讲了对每个点的碰撞检测以及响应，这里开始刚体的碰撞检测及响应：

如上，兔子表面每个点可以用质心+一个向量去表示，得到每个点再去遍历做检测即可。

碰撞点的位置速度如下公式：

但是问题是，我们之前讲过，我们的状态仅保存质心的四个变量（速度、位置、角速度、四元数），而对于表面上每一个点我们其实是不保存他们的状态的。因此利用前面impulse的方式去修改他们的值是没有意义的。

解决方法其实很常规，就是把冲量带来的影响直接叠加在质心的速度、角速度上：

然后就可以接着去计算每个点的速度了。

接下来利用叉乘的这个公式：

我们就可以进一步转化当前表达式：

简单意思就是，冲量 -> 造成点的速度的改变 ， 利用上图所示的这些关系把 j 算出来。

理一下整体流程：

首先去对每个顶点计算出每个顶点的位置，然后去检测一下是否有碰撞（根据sdf值），若是有碰撞则进一步去把碰撞点的速度算出来，然后判断顶点是否继续往下掉（是否继续穿透）；如果还要继续穿透就要继续处理，那么这时就要去算 Vi new ，计算方式：先分解出法向和切向，然后去计算关于摩擦的系数a，接着得到质心的新的normal方向和tangent方向的速度，于是就可以接着计算出 Vi new；这里计算 Vi new 那一块是灰的，原因是 Vi new 只是一个中间变量，实际上应该更新线速度和角速度；更新方法是用冲量的方式去更新，先计算出 K ，然后得到冲量。得到冲量后，回到 v 和 ω 去更新。

实现细节：

如果有很多点就用他们位置的平均值；由于重力的影响，会有抖动 oscillation 现象（在表面上，由于重力和所施加的反弹的力，会有反复抖动），从美观的角度可以去做如PPT所示的一个小hack，去衰减这样一个抖动。

事实上，处理这么多点不见得很慢，因为所有点只是都去检测了一遍是否碰撞，真正要处理的只有那些要碰撞的点。

Shape Matching

先假设这些点都不是刚体，简单的说就是先直接对每个点进行模拟，然后再把这些点恢复成一个刚体。因此关键就是恢复成一个刚体的过程：

因此在Shape Matching中，质心位置c和旋转矩阵R我们是不知道的，而模拟后每个点的坐标我们是知道的（记为 yi ，如下图）：

我们要恢复成刚体，自然希望这个范数表达式最小，进一步处理（对c求偏导）整理就得到了上图最下侧表达式。这个表达式很直观，其实也不难想象，质心就直接是这些点坐标的平均值。再对矩阵A求导得到：

得到A之后再分解成旋转矩阵和一个对称矩阵（形变的部分）。这个旋转矩阵就是我们想要的部分，即我们要得的旋转矩阵。至于为什么我们回到SVD中：

如上图下面做一个小的变形，前面三步叠加实际上就是S矩阵，效果是在本地做了一个形变（如图展示的效果就是沿着黑色的轴做出了形变）。

因此总结一下流程：

先对每个点独立地去进行更新模拟，得到每个点的位置y后，再去根据之前所述的公式把质心c和旋转矩阵R搞出来，然后更新速度和位置。

缺点在于很难严格地保证约束，且很难保证所有的约束都满足（可能需要迭代的方式反复约束）。一般当模拟精度不太高的时候（比如模拟衣服和扣子，刚体接触不太频繁）就可以考虑用 Shape Matching 。

Lab1

放一点我个人心得：

由于Unity数学库有点不完善，有一点小trick。

回到第四节课讲的实现细节中：

第一点，要注意用平均值去算的时候，顶点速度v_i 以及 current space 下的R_r_i 检测完都要取平均，带入之后计算。

第二点，我是简单写为restitution *= 0.5f;就可以了。要注意 U_T 和 U_N 要区分开。

由于Unity矩阵库无减法，我是这样做的：

Matrix4x4 K = new Matrix4x4(
	K_minuend.GetColumn(0) - K_subtrahend.GetColumn(0),
	K_minuend.GetColumn(1) - K_subtrahend.GetColumn(1),
	K_minuend.GetColumn(2) - K_subtrahend.GetColumn(2),
	K_minuend.GetColumn(3) - K_subtrahend.GetColumn(3)
);

minuend被减数，这里是被减矩阵；subtrahend减数。即我构造了俩矩阵，之后再直接用new和GetColumn方法构造。

同时矩阵库还没有常数乘以矩阵，我于是这样构造的K_minuend：

Matrix4x4 K_minuend = Matrix4x4.Scale(new Vector3(1.0f / mass, 1.0f / mass, 1.0f / mass));
K_minuend.m33 = 1.0f / mass;

注意m33一定要填上值，我之前就是这里忘记了，导致debug了很久。若不然对之后的Matrix4x4.Inverse会有影响。

还有四元数貌似没有叉乘，我是使用了四元数叉乘公式直接写的，在第三节课引入四元数有讲。

主要参考的PPT：

上图需要在Collision_Impulse函数中一条条对着写，注意要用到平均值去做计算，就需要把 Vi_new 和 后面的 Rr_i 都做替换。