最新文章

MIT 6.837:Particle System 粒子系统

MIT 6.837:Particle System 粒子系统

在接下来的若干篇博客中,我们将介绍用于实现计算机动画的若干技术。粒子系统是一种基于有限元思想的技术,将复杂的系统转化为有限个粒子模型进行模拟,再利用模拟结果重建出近似的系统模拟结果。

计算机动画的实现方法

关键帧动画

31 - MIT 6.837:Particle System 粒子系统

动画师通过指定动画中的部分帧中物体的状态,再通过插值方法生成状态间的过渡帧。

程序生成动画

将物体的运动描述为算法或函数,例如圆周运动用圆的方程表示、球的弹跳用带阻尼衰减的三角函数表示等。

基于物理的动画

记录物体的物理属性,利用物理规律模拟物体接下来的运动状态。这种方法与现实的物理规律相符,但难以控制其动画效果,模拟过程计算量开销也很大。

动作捕捉

32 - MIT 6.837:Particle System 粒子系统

在真实的物体上粘贴光学标记,再使用高速摄像机拍摄真实世界中的动画,利用这些数据在计算机中重建运动过程。常用于制作需要自然感或细致的动画,但需要在现实中利用设备捕捉动作,成本和难度较高。

粒子系统

模拟物理的过程中,常用的物理元有质点、刚体和变形物体。粒子系统就将要模拟的系统看做一系列质点组成,并模拟这些质点之间的物理规律,从而得到系统结果的近似解。

粒子系统通常生成许多粒子,用空间中的力场描述外界对粒子施加的力,并应用物理学中的规律求解例子状态。这些状态可以包括空间位置、速度、加速度、质量、生命周期等。

粒子的运动

根据经典的牛顿定律,运动涉及的参数包括质量 $m$、位置 $\boldsymbol{x}$、速度 $\boldsymbol{v}$ 和受力 $\boldsymbol{F}$。通过牛顿第二定律,它们建立起了这样的一组关系

$$ \begin{aligned}
& \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{x}(t) = \boldsymbol{v}(t) \\
& \frac{\mathrm{d}}{\mathrm{d}t} \boldsymbol{v}(t) = \frac{1}{m} \boldsymbol{F}(\boldsymbol{x}, \boldsymbol{v}, t)
\end{aligned} $$

这是一组常微分方程(Ordinary Differential Equation,ODE),记 $\boldsymbol{X}(t) = \begin{pmatrix} \boldsymbol{x} \\ \boldsymbol{v} \end{pmatrix}, \boldsymbol{f}(\boldsymbol{X}, t) = \begin{pmatrix} \boldsymbol{v} \\ \dfrac{1}{m} \boldsymbol{F} \end{pmatrix}$,则可以表示为规范化的写法

$$ \frac{\mathrm{d}\boldsymbol{X}}{\mathrm{d}t} = \boldsymbol{f}(\boldsymbol{X}, t) $$

ODE 的解是一个困难的问题,有些 ODE 可以利用通解公式简单求值,但也有许多 ODE 没有解析解。而对于动画中的粒子运动,我们其实只是想知道某些 $t$ 值上的解,因此不需要关注它是否存在解析解,而只需要一个数值解,甚至一个近似解,就可以实现动画效果了。

微分方程的数值解

对于计算机动画来说,我们只需要知道在一些离散时刻状态的值,再将其渲染为帧即可。因此,ODE 的解析解不是必要的,我们只需要有办法得到这些离散时刻的数值解即可。而由于帧间时间间隔很小,利用导数估计下一帧的状态是一个可能的思路。

Euler 方法

欧拉方法直接将时间间隔与导数相乘,即

$$ \begin{aligned}
& t_1 = t_0 + h \\
& \boldsymbol{X}_1 = \boldsymbol{X}_0 + h\boldsymbol{f}(\boldsymbol{X}_0, t_0)
\end{aligned} $$

此处时间间隔 $h$ 取得越小,估计的错误则越小。由于在整个迭代过程中错误会积累,如果 $h$ 的值取的不合适,得到的结果可能与实际结果相去甚远。另外,Euler 方法存在一些无论 $h$ 取得多小都会出错的情况,例如圆周运动。

34 - MIT 6.837:Particle System 粒子系统

Euler 方法的收敛性也无法保证,例如对于下面这个例子:

35 - MIT 6.837:Particle System 粒子系统

这些错误是由于 Euler 方法是一阶估计方法,其错误为 $O(h^2)$,且在整个迭代过程中错误会积累起来。

二阶方法

一阶方法的准确度不佳,因此我们可以使用二阶估计方法,降低 $h$ 取值对估计准确度的影响,使误差缩小到 $O(h^3)$。

36 - MIT 6.837:Particle System 粒子系统

改进 Euler 方法:利用一阶估计终点处导数的值补偿估计误差,令 $\boldsymbol{f}_0 = \boldsymbol{f}(\boldsymbol{X}_0, t_0), \boldsymbol{f}_1 = \boldsymbol{f}(\boldsymbol{X}_0+h\boldsymbol{f}_0, t_0+h)$,则 $\boldsymbol{X}(t_0+h) = \boldsymbol{X}_0 + \dfrac{h}{2}(\boldsymbol{f}_0+\boldsymbol{f}_1)$。

37 - MIT 6.837:Particle System 粒子系统

中点方法:利用步长中点处导数的值补偿误差,令 $\boldsymbol{f}_0 = \boldsymbol{f}(\boldsymbol{X}_0, t_0), \boldsymbol{f}_m = \boldsymbol{f} \left( \boldsymbol{X}_0+\dfrac{h}{2}\boldsymbol{f}_0, t_0+\dfrac{h}{2} \right)$,则 $\boldsymbol{X}(t_0+h) = \boldsymbol{X}_0 + h\boldsymbol{f}_m$。

改进 Euler 方法使用两个全步长估计的均值,而中点方法则是先取半步长导数再估计全步长,两者的结果不一定相同,但具有同阶数的精确度。

粒子动画的流程

伪代码如下:

AnimateParticles(n, y0, t0, tf) {
    y = y0
    t = t0
    DrawParticles(n, y)
    while(t != tf) {
        f = ComputeForces(y, t)
        dydt = AssembleDerivative(y, f)
            //there could be multiple force fields
        {y, t} = ODESolverStep(6n, y, dy/dt)
        DrawParticles(n, y)
    }
}

粒子的外部受力

重力:$\boldsymbol{F} = \begin{pmatrix} 0 \\ 0 \\ -mg \end{pmatrix}$,很轻的粒子(烟尘等)和火焰可以忽略重力作用

阻力:$\boldsymbol{F} = -d \boldsymbol{v}$,可以利用阻力让数值解稳定下来,也可以模拟在粘稠液体中的运动效果

空域力场:力与空间坐标相关,例如风推力/阻力、电磁力、漩涡等,有时需要改变物体的能量则速度也作为参数

相互作用力:与两粒子的空间坐标相关,例如用于模拟液体内部的作用力 $\boldsymbol{F} = \dfrac{k_1}{\left| \boldsymbol{x}^{(i)} – \boldsymbol{x}^{(j)} \right|^m} – \dfrac{k_2}{\left| \boldsymbol{x}^{(i)} – \boldsymbol{x}^{(j)} \right|^n}$(Lennard-Jones 力)

粒子的碰撞

38 - MIT 6.837:Particle System 粒子系统

粒子间的碰撞通常忽略不做处理,一般只做与环境的碰撞检测。碰撞处理大致可以分为两种,一种为模拟真实的反弹,一种只是将其移出物体表面。与平面的求交可以类似之前射线的方法,计算平面方程的值即可。

39 - MIT 6.837:Particle System 粒子系统

在碰撞后,其速度变化规律有多种处理方式。理想反弹时,其切向速度不变,而法向速度反向,也即 $\boldsymbol{v}’ = \boldsymbol{v}_t – \boldsymbol{v}_n$。也可以引入系数 $\varepsilon$ 反应非理想反弹时法向速度的损失,即 $\boldsymbol{v}’ = \boldsymbol{v}_t – \varepsilon \boldsymbol{v}_n$。

粒子的生命周期

利用生命周期,我们可以控制粒子实例的数量,移除超出屏幕或世界范围的粒子等。此外,还可以用生命值来实现颜色、透明度等的改变,实现变色、淡出等效果。

粒子的应用

除粒子模拟火焰、粉尘等效果之外,粒子还可用于实现其他一些效果。

粒子建模

40 - MIT 6.837:Particle System 粒子系统

模拟鸟群、人群

40 1 - MIT 6.837:Particle System 粒子系统
40 - MIT 6.837:Particle System 粒子系统

游戏中的系统和效果

天气、爆炸、草丛等。

MIT 6.837:Sampling & Aliasing 采样与混叠

MIT 6.837:Sampling & Aliasing 采样与混叠

数字图片的采样和量化

真实世界中的场景是连续的,而计算机一般处理离散的数据。

将连续函数变为离散函数的过程称为采样,而将连续的数值变为离散的数值的过程为量化。为了得到适于计算机处理的数据形式,我们必须对原始场景的图片进行采样和量化。

对连续的图片进行采样的过程可以用采样矩阵表示。采样矩阵的基础是二维冲激函数 $\delta(x,y) = \begin{cases} 1, & (x,y)=(0,0) \\ 0, & \text{otherwise} \end{cases}$,可以用冲激函数表示采样矩阵为

$$ \sum_{j=0}^{h-1} \sum_{i=0}^{w-1} \delta(u-i, v-j) $$

将这个矩阵与表示连续图像的连续二维函数相乘,得到二维函数在整点上的颜色值,如下图所示。

.jpg - MIT 6.837:Sampling & Aliasing 采样与混叠

采样与重建

通过傅里叶分析,我们可以知道采样是频域上的周期搬移,而重建是频域上的低通滤波。但由于信号中包含高频,搬移时会发生相邻周期的频谱混叠,而这种混叠将导致重建时无法还原出原始的一个周期,最终体现为图片中的错误与失真。

2 1024x816 - MIT 6.837:Sampling & Aliasing 采样与混叠

为了保证重建中不发生问题,我们需要对信号进行预处理,例如移除其中的高频部分,或增加采样密度。根据采样定理,采样频率应不小于最高频率的两倍,只要预处理能够使信号满足此条件,则不会发生混叠失真。

滤波

在计算机图形学中,滤波器被用于两种情况:重建信号、去除高频部分得到带限信号。这两者都是低通滤波器,通常可以使用同一种滤波器完成。

二维的滤波可以通过矩阵的卷积来计算,因此滤波器可以用卷积核矩阵来表示。

预处理

预处理滤波器用于滤去高频信号,故为低通滤波器,可用覆盖像素面积的比例来实现。例如,光栅化直线的时候,应用预处理滤波器平滑直线边缘的锯齿。

.png - MIT 6.837:Sampling & Aliasing 采样与混叠
1 - MIT 6.837:Sampling & Aliasing 采样与混叠

后处理

后处理坐标系用于对样本滤波,可用对附近像素求平均值实现,也可以实现与预处理滤波器类似的效果。

2 - MIT 6.837:Sampling & Aliasing 采样与混叠

滤波器的种类

理想滤波器:频域低频增益为 1,高频衰减至 0。时域是无穷信号,难以实现。

实际滤波器:存在低频衰减和高频泄露问题,依然可能存在混叠失真。

3 - MIT 6.837:Sampling & Aliasing 采样与混叠

高斯滤波器:CRT 显示器在人眼中成像与之类似。

方框滤波器:时域上是一个门函数,上采样时造成锯齿。对应最近邻插值。

三角滤波器:时域上是一个三角形,易于实现,效果比方框略好。对应双线性插值。

双三次插值:时域上更接近理想滤波器的曲线,效果更好。

纹理贴图反走样

由于纹理贴图是分辨率有限的图片,在一些情况下容易产生走样。利用滤波等原理,我们可以改善这些走样问题。

纹理被放大时,如果使用方框滤波器会产生锯齿,因此可以利用三角滤波器(双线性插值)代替,从而得到平滑的结果,但会带来模糊的效果。

纹理被缩小时,其中的高频分量会导致奇怪的效果(如摩尔纹等),因此可以利用低通滤波器去除其高频分量。可以预处理出多种尺寸的缩小滤波结果,以减少渲染时的计算开销。

MIP Map 是一种保存缩小纹理的方式,它通过预处理出 1/2、1/4、……尺寸的缩小结果,再将其拼在同一张图中。由等比数列求和,新的纹理图片仅为原图片的 2 倍大小。

4 - MIT 6.837:Sampling & Aliasing 采样与混叠

由于每次映射贴图时不一定是等比例缩放,需要沿各边缩放的结果,可以利用各向异性 MIP Map 来处理这种情况。在这种方法中,保存的缩小贴图各边缩小为 1/2、1/4、……等尺寸,图片大小为原图的 4 倍。

4 - MIT 6.837:Sampling & Aliasing 采样与混叠
MIT 6.837:Texture Mapping 纹理映射

MIT 6.837:Texture Mapping 纹理映射

当遇到表面细节复杂的几何体时,利用几何(三角面)来表示这些细节需要巨大的计算开销,从而导致性能瓶颈。一种折中的办法是在平面的纹理上表示这些表面细节,而将纹理图片映射到几何体表面上,大大简化了几何的复杂程度。

二维纹理映射

透视投影后的插值

应用纹理到三角面上通常在投影完成后(fragment shader)进行,此时三角面不在具有空间中的坐标,而只保存了屏幕空间坐标与 Z 缓冲区(或 1/z 值)。因此,纹理的应用必须依赖屏幕空间坐标进行。如果直接对屏幕空间中的坐标进行线性插值,会得到怪异的结果,如下图所示。

50 1024x411 - MIT 6.837:Texture Mapping 纹理映射

在这个例子中,我们把一个方形的纹理映射到两个三角面上,这两个三角面拼成了一个斜的方形面。这两个三角面上的纹理的扭曲不相同。

一种解决此问题的方法是将大三角面细分为小三角面,以让三角形邻接处的扭曲程度接近,减少这种怪异的明显程度。但此方法的效果并不好,我们有其他办法从根本上解决此问题。

这种问题产生的原因是,在透视投影中,屏幕空间上的线性插值在空间中由于 z 坐标原因并不是线性分布的。这个问题在可见性中的 Z 缓冲区话题第一次被提及,而此次我们希望将插值运用在纹理映射中。

51 - MIT 6.837:Texture Mapping 纹理映射

我们可以通过比较在屏幕空间中的线性插值与世界空间中的线性插值,以得到两个插值之间的关系。假设只考虑 xOz 屏幕,则空间中的一条直线被投影到 z=1 上为 x 坐标的一段区间,屏幕空间中 x 坐标即可表达一个点,而世界空间中需要用一个有序对 (x, z) 来表达一个点。

$$
\begin{align}
& 屏幕空间: p(t)=p_1+t(p_2-p_1)=\frac{x_1}{z_1} + t \left( \frac{x_2}{z_2}-\frac{x_1}{z_1} \right) \\
& 世界空间: \begin{bmatrix} x\\z \end{bmatrix}=\begin{bmatrix} x_1\\z_1 \end{bmatrix} + s \left( \begin{bmatrix} x_2\\z_2 \end{bmatrix}-\begin{bmatrix} x_1\\z_1 \end{bmatrix} \right) \\
& 世界空间中点在屏幕空间的投影: p\left( \begin{bmatrix} x\\z \end{bmatrix} \right) = \frac{x_1+s(x_2-x_1)}{z_1+s(z_2-z_1)}
\end{align}
$$

联立屏幕空间与世界空间投影的插值表达式,可得 s 与 t 的关系如下

$$ s=\frac{tz_1}{z_2+t(z_1-z_2)} $$

如果 Z 缓冲区保存的是 1/z 或 w 值,则上式还可以写成

$$ s=\frac{tw_2}{w_1+t(w_2-w_1)} $$

因此,在插值参数为 t 的像素时,我们可以利用上式转换为空间中的插值参数 s,通过 s 的值来获取纹理中的颜色值,这样就可以得到正确的投影效果。

纹理的获取与使用

纹理的获取渠道很多,可以利用现实物体的照片生成纹理,也可自行绘制。如果使用照片生成纹理,直接将物体的三角面与照片中的一部分映射是有问题的,因照片本身也是一种投影结果,投影中会发生拉伸。因此,可以将照片再次投影到同一角度的几何体上,从而得到正常的纹理效果。

52 - MIT 6.837:Texture Mapping 纹理映射

如果使用自己绘制的纹理,则需要指定每一个三角面的纹理,如面数较多,纹理图片的数量也较多,不易管理。此时,可以将一些纹理拼成一张图片,并在模型中指定顶点与纹理图片坐标的对应关系。

52 1 - MIT 6.837:Texture Mapping 纹理映射
52 2 1024x714 - MIT 6.837:Texture Mapping 纹理映射

一种应用纹理颜色的方法是,将 Phong 模型中的漫反射系数替换为纹理图片中的采样值,从而将纹理与光照效果结合起来。

52 - MIT 6.837:Texture Mapping 纹理映射

特殊的纹理

程序化生成纹理

如果能够得到一个函数 $f(x,y,z)$,其输入为要渲染位置的坐标,而输出为颜色,这样的函数也可以看做是纹理。上面基于图片的纹理可以看做是函数值离散化的函数,而如果规定了连续的函数,则纹理可以做到无限大分辨率。

这种纹理通常与噪声算法结合,用于生成一些具有一定随机性同时有一定规律的纹理,如树干、瓷器、动物毛皮花纹等。

53 - MIT 6.837:Texture Mapping 纹理映射

法线贴图

实际的渲染中,我们经常遇到需要在较简单的几何体上实现表面细节的需求。除了常规的材质颜色外,还有光照阴影等效果,这可以通过修改表面的法线的来完成。

一种方法是 Gouraud 渲染,在计算顶点值时,它先计算顶点周围面的平均法向量,再计算颜色,最后在三角面上对顶点颜色进行插值,得到一种较平滑的颜色效果。

另一种方法是 Phong 法线插值,它的插值对象是法线本身,对于每一个渲染的像素插值其法线,再计算颜色值。

54 - MIT 6.837:Texture Mapping 纹理映射

还可以直接在纹理中指定法线方向,利用图片的 RGB 三个通道来表示世界空间中的三维方向向量坐标,直接用纹理中的向量来代替原本的法向量,这样的方法叫做法线贴图。

另一种与法线贴图类似的方法叫凹凸贴图,它指定顶点需要看起来凹凸的程度,通过在此贴图上计算梯度值来改变法向量,形成看起来凹凸的情况。但这种方法无法正确处理自遮挡产生的阴影,因此凹凸程度不宜设置过大。

55 - MIT 6.837:Texture Mapping 纹理映射

置换贴图

置换贴图是凹凸贴图的改进。为了解决凹凸贴图无法正确处理自遮挡的问题,置换贴图应用时将先移动几何体中点的位置,再将新的顶点坐标作为之后阶段渲染的输入。由于置换贴图需要在世界坐标系修改顶点位置,其与上面的方法的应用阶段不同,一种方法是在几何的曲面细分阶段中完成。

56 - MIT 6.837:Texture Mapping 纹理映射

环境贴图

在渲染包含反射项的物体时,反射光线的起点可看成一个新的相机。因此,我们可以预先根据反射光线得到周围环境的图像,再将其直接应用在物体的反射项中,物体较小时可以假设反射线从同一个点发出。

56 - MIT 6.837:Texture Mapping 纹理映射

光照贴图

全局光照的计算开销很大,而场景中的很多物体是静止的。光照贴图通过预先计算那些静止物体的光照强弱,暂时缓存这一部分的光照结果,只需要处理动态物体的光照并合并结果即可。

纹理的混叠

混叠(Aliasing)指的是采样导致频域上不同周期之间的信号叠在一起,导致采样结果发生畸变。

纹理图片的分辨率是有限的,每次获取颜色时要对纹理图片这一离散信号进行采样,如采样率太低则容易发生混叠现象。

57 1024x487 - MIT 6.837:Texture Mapping 纹理映射

对于纹理来说,一种减少混叠的方法是使用 mipmap。