Position Based Dynamic (一)

本文简要介绍了目前游戏工业界中软体模拟的一种常用方法Position Based Dynamic的原理以及实现。本文分为两个部分：原理以及实现。原理部分主要解读2005年的论文[1]，简要介绍了Position Based Dynamic的思想和原理。在实现部分，本文给出了一个基于Position Based Dynamic框架下的简单三角形模拟程序并详细解读了其代码。

原理

​ 非刚性物体的模拟在当下的游戏产业日渐成为重要的组成部分 (懂得都懂)。相比传统模拟，游戏产业对求解结果的精确性并没有那么高的要求，相反如何快速高效的求出视觉上具有说服力的结果是游戏相关模拟技术关注的重点。因此，Position Based Dynamic(PBD)应运而生。相比于有限元，Material Point Method等更加准确的方法，PBD更关注数值稳定性以及模拟结果的说服力，这点在PBD如何处理碰撞上表现的尤为突出。相比传统方法中通过力的作用修正碰撞结果，PBD在处理碰撞时往往会先估计碰撞后系统的状态，再通过碰撞结果更新当前系统的状态。

PBD的世界

​ 对于PBD方法，世界是由质点以及约束构成的。质点和中学物理中学到的质点性质相同，其包含位置，速度以及质量三个性质，其中位置和速度会受到外力作用也会受到约束限制。一个约束包含受约束的质点集 ,劲度系数以及约束方程以及条件 。约束方程为 的函数，它可以被看成是系统为了维持稳定所定义的势能函数。具体的其可以分为等式约束以及不等式约束两类，其中，等式约束需要系统满足 ,而不等式约束则需要满足 。对于不等式约束，可以分为两个步骤：1.检测当前状态是否满足 2. 对于不满足1.的约束将其转为 并和其它等式约束一同求解。

​ 至此我们以及得到构建PBD世界的全部元素，质点负责记录世界中每个物体的位置以及状态，外力负责改变质点的运动状态，约束负责制约质点运动以保持物体外形。

PBD的基本步骤

​ PBD的基本求解步骤如算法1所示（ 下同）:

1. 对于每个质点
3. 对每个质点
4. 对每个质点
6. 对每个质点

​ 其中，可以看到系统除了需要记录质点的位置，速度，还有一个额外值 负责记录PBD预测的该质点在下一步迭代中的初始位置。最后该位置会在步骤6中被用于逆向更新质点的速度。对于步骤2.中的dampVelocities,其主要步骤如算法2所示:

1. 计算系统的质心坐标以及速度

2. 计算系统的平均角速度

3. 让各个质点角速度于系统相对一致

其中 , 满足 。该衰减项能让系统不减少整体能量的同时消弭系统的形变。

对于第4步中的generateCollisionConstrains以及第5步中的projectConstraints将在1.4以及1.3节中详细介绍。

投影约束

由于系统所有约束都被转为等式约束 ，因此我们在求解时可以只考虑等式情况。在理想情况下，系统的所有约束一起构成了一个 的非线性方程组，直接求解这个线性方程组是十分困难的。因此，[1]提出用类似G-S法的迭代方式逐步求解方程组。在每个迭代步骤里，PBD需要逐个遍历每个等式约束 并求出其近似解。假设求解前的状态为 ,对约束方程泰勒展开有:

$$C(p)=C(p^0)+\Delta p \cdot \bigtriangledown_p C(p^0)=0 (1)$$

该方程是欠约束的，为此，我们不妨假设 的方向和约束函数梯度方向相同 ,带入(1)中得到:

$$C(p^0) + \lambda (\bigtriangledown_p C(p^0)^T \bigtriangledown_p C(p^0))=0$$ $$\lambda = -\frac{C(p^0)}{\bigtriangledown_p C(p^0)^T \bigtriangledown_p C(p^0)} (2)$$

因此，对于每个质点的位置 我们有:

$$\Delta p_i = -\frac{C(p^0)}{\sum_j \bigtriangledown_{p_j} C(p^0)^T \bigtriangledown_{p_j} C(p^0)} \bigtriangledown_{p_i} C(p^0)$$

对于各个质点质量不相同的系统，作者给出的解决方法是给每个距离乘上权重:

$$\Delta p_i = - \frac{ C(p^0)}{ \sum_j w_j \bigtriangledown_{p_j} C(p^0)^T \bigtriangledown_{p_j} C(p^0)} w_i \bigtriangledown_{p_i} C(p^0) (3)$$

在每次迭代中，PBD会遍历每个约束中并给其中的每个质点 添加一个约束投影距离 。由于这个过程必须在每个约束求解完成后才能开始下一个约束的求解，因此这个过程很难被并行化。[2] 提出通过图着色的方法给约束聚类。他们首先构建了以约束为节点，在两个有公共质点的约束之间创建一条边的计算图。接着，通过图着色算法，同一颜色的约束可以放在同一聚类中，之后便可以以聚类为单位平行求解约束方程。

对于非刚性的约束，其劲度系数 ，因此在计算的偏移过程中，我们需要乘上系数 :

$$p_i^{j + 1} = p_i^j + k \Delta p^j_i (4)$$

生成碰撞约束

在所有约束中，有一种特殊约束为碰撞约束。该约束在一个物体于另一个物体碰撞时产生，保证一个物体的质点不会进入另一个物体的内部。假设物体的碰撞点为 ,碰撞表面的法线为则该约束可以写为 。正如1.1中的讨论结果，我们可以把约束转为等式约束

$$(p_i - q_i)\cdot n_i = 0 (5)$$

带入（3）中求解后 为:

$$\Delta p_i = - \frac{(p_i - q_i) \cdot n_i}{ w_i n_i \cdot n_i} w_in_i= -((p_i - q_i)\cdot n_i) n_i (6)$$

可以看到，该结果有十分明确的几何意义:对于进入碰撞点的质点 将其沿着法线方向 移动其到碰撞点距离 在 方向上的投影个距离。

在 步骤中，对每个顶点，PBD会发射一个从 的射线。此时，可以分为三种情况

1. 若 均在物体表面之外，则约束 被满足
2. 若 经过物体表面，令其与物体表面交点为,法线为，则生成一个约束方程为 劲度系数为的 的约束
3. 若 均在物体表面内，则找到 距离物体表面最近的交点 以及其法线,生成则生成一个约束方程为 劲度系数为的 的约束

对于动态物体，则可以通过同时将两物体纳入PBD系统中模拟来实现。例如假设物体2有顶点 ，则可以找到对应碰撞点 以及法线 ，将其带入公式(5)中有:,将该公式带入(3)中便能求解该约束。

实现

一个利用PBD框架实现的刚性三角形模拟代码如下（这个实现起演示作用，并不高效）：

import taichi as ti

ti.init(arch=ti.gpu)  # 使用GPU

## 约束的类型编号
# 等距离约束，数据编码:(i,j,d,_), i,j为两点位置, d为两点距离，只有当 |xi - xj| == d时，约束被满足
c_s_eq_distance = 1

c_c_base = 1024
# 面碰撞约束，数据编码为(i, qx, qy, nx), i为点位置, (qx, qy)为碰撞点qi坐标，nx为碰撞点法线x分量,ny = sqrt(1 - nx ** 2)
# 只有当 <(xi - qi) ,n> > 0  时约束被满足
c_c_ge_surface = c_c_base + 0
## 约束的类型编号

# damp 系数
k_damp = 1e-2
# 场景中三角形数量
n_triangle = 1
# 场景中质点数量
n_vertex = n_triangle * 3
# 模拟步长
dt = 1e-2

# 重力加速度
g = ti.Vector([0, -1])

# 水平面位置
plane = 0.2

# 场景的约束数组容量
n_max_constraint_cnt = n_triangle + n_vertex * 5
# 投影操作的迭代次数
n_max_solve = 5

# 512x512的画布
pixels = ti.Vector.field(3, dtype=ti.f32, shape=(512, 512))

# 质点状态p,x,v
p = ti.Vector.field(2, dtype=ti.f32, shape=(n_vertex))
x = ti.Vector.field(2, dtype=ti.f32, shape=(n_vertex))
v = ti.Vector.field(2, dtype=ti.f32, shape=(n_vertex))
# 质点颜色
c = ti.Vector.field(3, dtype=ti.f32, shape=(n_vertex))

# 约束数组con以及约束的类型数组con_t
con = ti.Vector.field(4, dtype=ti.f32, shape=(n_max_constraint_cnt))
con_t = ti.field(dtype=ti.f32, shape=(n_max_constraint_cnt))

# 系统过去的总时间
acc_t = ti.field(dtype=ti.f32,shape=1)
# [当前所有约束数量, 当前静态约束数量] con数组中从con_cnt[1]开始便是动态生成的碰撞约束
con_cnt = ti.field(dtype=ti.i32, shape=2)

acc_t[0] = 0
con_cnt[0] = 0
con_cnt[1] = 0

@ti.kernel
def updateVelocity():
    for i in v:
        v[i] = v[i] + dt * g

def dampVelocity():
    xc = ti.Vector([0, 0])
    vc = ti.Vector([0, 0])
    for i in range(n_vertex):
        xc += x[i]
        vc += v[i]
    xc /= n_vertex
    vc /= n_vertex

    L = 0
    I = 0
    for i in range(n_vertex):
        ri = x[i] - xc
        vi = v[i] - vc
        L += ri.x * vi.y - ri.y * vi.x
        I += ri.x * ri.x + ri.y * ri.y

    w = L / I
    
    for i in range(n_vertex):
        vi = v[i]
        ri = x[i] - xc
        v[i] = vi + k_damp * (vc + w * ri - vi)


    
@ti.kernel
def updatePosition():
    for i in v:
        p[i] = x[i] + dt * v[i]

@ti.kernel
def generateConstraints():
    for i in v:
        pi = p[i]
        xi = x[i]
        di = pi - xi

        if di[1] < 0 and pi[1] < plane:
            t = (plane - xi[1]) / di[1]
            qi = t * di + xi
            idx = ti.atomic_add(con_cnt[0], 1)

            con[idx] = ti.Vector([i, qi[0], qi[1], 0])
            con_t[idx] = c_c_ge_surface

def projectConstrainsPositions():
    for ci in range(con_cnt[0]):
        con_type = con_t[ci]
        data = con[ci]
        if con_type == c_s_eq_distance:
            i = int(data[0])
            j = int(data[1])
            d = data[2]

            p1 = p[i]
            p2 = p[j]

            norm_p1_p2 = (p1 - p2).norm()
            n = (p1 - p2) / norm_p1_p2
            del_p1 = - (norm_p1_p2 - d) / 2 * n
            del_p2 = (norm_p1_p2 - d) / 2 * n

            p[i] += del_p1
            p[j] += del_p2
            
        elif con_type == c_c_ge_surface:
            # 面碰撞约束，数据编码为(i, qx, qy, nx)
            i = int(data[0])
            qi = ti.Vector([data[1], data[2]])
            nix = data[3]
            ni = ti.Vector([nix, ti.math.sqrt(1 - nix * nix)])

            pi = p[i]
            pq = pi - qi

            p[i] += - (pq[0] * ni[0] + pq[1] * ni[1]) * ni

@ti.kernel
def updateState():
    for i in v:
        v[i] = (p[i] - x[i]) / dt
        x[i] = p[i]
    

def simulation():
    con_cnt[0] = con_cnt[1]

    updateVelocity()
    dampVelocity()
    updatePosition()
    generateConstraints()

    for _ in range(n_max_solve):
        projectConstrainsPositions()

    updateState()
    acc_t[0] += dt


# 计算任意一点P在三角形(p0,p1,p2)内的重心坐标
@ti.func
def barycentric_coords(p0, p1, p2, P):
    # 计算向量
    v0 = p2 - p0
    v1 = p1 - p0
    v2 = P - p0

    # 计算分母
    denom = v0.dot(v0) * v1.dot(v1) - v0.dot(v1) * v1.dot(v0)

    # 计算重心坐标
    u = (v2.dot(v0) * v1.dot(v1) - v2.dot(v1) * v0.dot(v1)) / denom
    v = (v2.dot(v1) * v0.dot(v0) - v2.dot(v0) * v0.dot(v1)) / denom
    
    return u, v

# 将三角形绘制到屏幕上
@ti.kernel
def paint_triangle():
    for i, j in pixels:  # 遍历所有像素
        # 将像素坐标转换为[0,1]范围内的坐标
        pos = ti.Vector([i / 512, j / 512])

        p0 = x[0]
        p1 = x[1]
        p2 = x[2]

        c0 = c[0]
        c1 = c[1]
        c2 = c[2]

        u,v = barycentric_coords(p0, p1, p2, pos)

        color = ti.Vector([0,0,0])
        if u > 0 and v > 0 and u + v < 1:
            color = c0 + (c1 - c0) * u + (c2 - c0) * v

        # 将颜色存储到像素中
        pixels[i, j] = color

gui = ti.GUI("Triangle", res=512)


# 初始化场景，将三角形的刚性约束添加到场景中
@ti.kernel
def init():
    x[0] = ti.Vector([0.37, 0.5])
    x[1] = ti.Vector([0.4, 0.43])
    x[2] = ti.Vector([0.45, 0.51])

    c[0] = ti.Vector([1, 0, 0])
    c[1] = ti.Vector([1, 0, 0])
    c[2] = ti.Vector([1, 0 , 0])

    for i in v:
        # 等距离约束，数据编码:(i,j,d,_), i,j为两点位置, d为两点距离，只有当 |xi - xj| == d时，约束被满足
        j = (i + 1) % 3
        con[i] = ti.Vector([i, j, (x[i] - x[j]).norm(),0])
        con_t[i] = c_s_eq_distance

        con_cnt[0] += 1
        con_cnt[1] += 1

# 初始化系统
init()

while gui.running:
    # 模拟
    simulation()

    # 绘制模拟结果
    paint_triangle()
    # 将画布内容输出到gui上
    gui.set_image(pixels)
    gui.show()

其运行效果如图所示:

场景描述

整个场景包含一个平面以及一个刚体三角形。三角形包含三个质点，每个质点两两之间存在一个等距离约束：

$$C(p_1,p_2)=|p_1 - p_2| - d_{12} = 0$$

其中 为 在初始状态下的距离。将三角形以及其约束添加到场景中的代码为:

@ti.kernel
def init():
    x[0] = ti.Vector([0.37, 0.5])
    x[1] = ti.Vector([0.4, 0.43])
    x[2] = ti.Vector([0.45, 0.51])

    c[0] = ti.Vector([1, 0, 0])
    c[1] = ti.Vector([1, 0, 0])
    c[2] = ti.Vector([1, 0 , 0])

    for i in v:
        # 等距离约束，数据编码:(i,j,d,_), i,j为两点位置, d为两点距离，只有当 |xi - xj| == d时，约束被满足
        j = (i + 1) % 3
        con[i] = ti.Vector([i, j, (x[i] - x[j]).norm(),0])
        con_t[i] = c_s_eq_distance

        con_cnt[0] += 1
        con_cnt[1] += 1

模拟整体框架

模拟的代码按照算法1给出的思路实现

def simulation():
    con_cnt[0] = con_cnt[1]

    updateVelocity()
    dampVelocity()
    updatePosition()
    generateConstraints()

    for _ in range(n_max_solve):
        projectConstrainsPositions()

    updateState()
    acc_t[0] += dt

在开始模拟前需要把上一帧碰撞检测添加的碰撞约束清空con_cnt[0] = con_cnt[1]

updateVelocity() 对应算法1的第一步。由于系统的唯一外力为重力，因此更新时只需加上重力加速度即可

@ti.kernel
def updateVelocity():
    for i in v:
        v[i] = v[i] + dt * g

dampVelocties() 对应算法1的第2步，其实现在2.3节中讨论

updatePosition() 对应算法1的第3步，其就是简单的根据速度更新位置

@ti.kernel
def updatePosition():
    for i in v:
        p[i] = x[i] + dt * v[i]

generateConstraints() 对应算法1 的第4步，其实现在2.4节中讨论

projectConstrainsPositions() 对应算法1 的第5步，其实现在2.5节中讨论

updateState() 对应算法1 的第6步，根据预测的下一帧位置更新质点位置和速度

@ti.kernel
def updateState():
    for i in v:
        v[i] = (p[i] - x[i]) / dt
        x[i] = p[i]

dampVelocities

dampVelocities按照算法2 的思路实现

根据算法2 的第一步，我们需要计算系统的质心坐标以及速度。由于我们假设当前系统每个质点的质量是1，因此:

xc = ti.Vector([0, 0])
vc = ti.Vector([0, 0])
for i in range(n_vertex):
    xc += x[i]
    vc += v[i]
xc /= n_vertex
vc /= n_vertex

根据算法2 的第二步，我们需要计算系统的平均角速度

对于一个二维平面上的系统，其转动惯量为 ,因此我们可以将公式改写为:

实现的代码为:

L = 0
I = 0
for i in range(n_vertex):
    ri = x[i] - xc
    vi = v[i] - vc
    L += ri.x * vi.y - ri.y * vi.x
    I += ri.x * ri.x + ri.y * ri.y

w = L / I

最后根据算法2 的第三步，我们需要对每个质点减去其超出整个系统平均角速度的部分。

for i in range(n_vertex):
    vi = v[i]
    ri = x[i] - xc
    v[i] = vi + k_damp * (vc + w * ri - vi)

generateConstraints

目前我们系统中只有三角形到平面一种情况，因此我们只需遍历所有顶点，判断其是否和三角形有交点即可。

为了简化计算过程，这里我们假设上一帧的位置 不在当前平面中。因此我们只需考虑 是否与平面有交点。

对于水平面 ,质点与其有交点需要满足两个条件:

1. 的y值低于水平线
2. 质点有向下运动的趋势

因此我们有:

if di[1] < 0 and pi[1] < plane:
	...

对于射线 我们可以用以下方程表示:

其与水平面 的交点满足

解得: ,其中

t = (plane - xi[1]) / di[1]
qi = t * di + xi

最后，我们通过原子加给该约束分配一个位置，并把对应信息写到buffer中:

idx = ti.atomic_add(con_cnt[0], 1)

con[idx] = ti.Vector([i, qi[0], qi[1], 0])
con_t[idx] = c_c_ge_surface

projectConstrainsPositions

目前我们系统中需要处理约束有两类: 等值约束以及碰撞约束:

if con_type == c_s_eq_distance:
   ...
elif con_type == c_c_ge_surface:
   ...

对于碰撞约束，我们通过公式(6)的结论直接求解

# 面碰撞约束，数据编码为(i, qx, qy, nx)
i = int(data[0])
qi = ti.Vector([data[1], data[2]])
nix = data[3]
ni = ti.Vector([nix, ti.math.sqrt(1 - nix * nix)])

pi = p[i]
pq = pi - qi

p[i] += - (pq[0] * ni[0] + pq[1] * ni[1]) * ni

对于等值约束，其约束方程为:

$$C(p_1,p_2) = |p_1 - p_2| - d = 0$$

带入公式(3)中解得:

$$\bigtriangledown_{p_1} C(p_1,p_2) = \frac{p_1 - p_2}{|p_1 - p_2|}, \bigtriangledown_{p_2} C(p_1,p_2)=\frac{p_2 - p_1}{|p_1 - p_2|}$$ $$\bigtriangledown_{p_i} C(p_1,p_2) ^T \bigtriangledown_{p_i} C(p_1,p_2) = 1 (i\in \{1,2\})$$ $$\Delta p_1 = - \frac{|p_1 - p_2| - d}{2} \frac{p_1 - p_2}{|p_1 - p_2|}$$ $$\Delta p_2 = - \frac{|p_1 - p_2| - d}{2} \frac{p_2 - p_1}{|p_1 - p_2|}$$

因此，其可以实现为:

i = int(data[0])
j = int(data[1])
d = data[2]

p1 = p[i]
p2 = p[j]

norm_p1_p2 = (p1 - p2).norm()
n = (p1 - p2) / norm_p1_p2
del_p1 = - (norm_p1_p2 - d) / 2 * n
del_p2 = (norm_p1_p2 - d) / 2 * n

p[i] += del_p1
p[j] += del_p2

参考文献

[1] M. Müller, Heidelberger B, Ratcliff M H A J. Position based dynamics[J]. Journal of Visual Communication and Image Representation, 2007(18).

[2] Fratarcangeli M, Pellacini F. Towards a Massively Parallel Solver for Position Based Dynamics[C]. //Proceedings of SIGRAD 2014.