第4章:碰撞检测算法

碰撞检测是物理仿真引擎的核心组件之一,直接影响仿真的准确性和效率。在具身智能应用中,高效而精确的碰撞检测不仅决定了机器人与环境交互的真实性,更是强化学习训练速度的关键瓶颈。本章将深入探讨从宽相筛选到窄相精确检测的完整流程,重点介绍GJK、EPA等经典算法的数学原理,以及连续碰撞检测在高速运动场景中的应用。我们还将讨论空间数据结构的选择如何影响大规模并行仿真的性能,为读者在实际工程中的算法选择提供理论指导。

学习目标

完成本章学习后,您将能够:

  • 理解并实现两阶段碰撞检测的完整流程
  • 掌握GJK算法的闵可夫斯基差原理及其几何直觉
  • 实现EPA算法计算穿透深度和接触法向
  • 设计适合不同场景的连续碰撞检测策略
  • 选择和优化空间数据结构以提升检测效率
  • 理解SDF在隐式表示中的应用及其梯度计算

4.1 宽相与窄相检测

碰撞检测的计算复杂度在最坏情况下为$O(n^2)$,其中$n$为场景中的物体数量。为了提高效率,现代物理引擎普遍采用两阶段策略:宽相(Broad Phase)快速筛选可能碰撞的物体对,窄相(Narrow Phase)精确计算实际接触信息。

4.1.1 宽相检测原理

宽相检测的核心思想是使用简化的几何表示(通常是轴对齐包围盒AABB)快速排除不可能发生碰撞的物体对。给定$n$个物体,宽相检测输出一个候选碰撞对集合$\mathcal{P} \subseteq \{(i,j) | 1 \leq i < j \leq n\}$,满足:

$$|\mathcal{P}| \ll \frac{n(n-1)}{2}$$ AABB重叠测试是最基本的宽相检测方法。对于两个AABB $B_1 = [x_1^{\min}, x_1^{\max}] \times [y_1^{\min}, y_1^{\max}] \times [z_1^{\min}, z_1^{\max}]$和$B_2$,重叠条件为: $$\forall k \in \{x,y,z\}: \max(k_1^{\min}, k_2^{\min}) \leq \min(k_1^{\max}, k_2^{\max})$$

4.1.2 扫描线算法(Sweep and Prune)

扫描线算法通过维护物体在各坐标轴上的有序投影区间,将重叠检测的平均复杂度降至$O(n\log n)$。算法维护三个有序列表,每个列表存储物体在对应轴上的区间端点:

轴向投影结构:
X轴: [(x_min, id, BEGIN), (x_max, id, END), ...]
Y轴: [(y_min, id, BEGIN), (y_max, id, END), ...]  
Z轴: [(z_min, id, BEGIN), (z_max, id, END), ...]

当物体$i$移动时,只需更新其6个端点在有序列表中的位置。通过增量式更新,对于运动连贯的场景,每帧的更新复杂度接近$O(n)$。

4.1.3 空间哈希

对于包含大量小物体的场景,空间哈希提供了更好的局部性。将空间划分为大小为$h$的网格单元,物体$i$根据其AABB覆盖的单元被插入哈希表: $$\text{hash}(x,y,z) = (p_1 \cdot x \oplus p_2 \cdot y \oplus p_3 \cdot z) \bmod M$$ 其中$p_1, p_2, p_3$为大质数,$M$为哈希表大小。单元大小$h$的选择需要平衡:过小导致物体跨越多个单元,过大则降低筛选效率。经验上,$h \approx 2\bar{r}$,其中$\bar{r}$为物体平均半径。

4.1.4 窄相检测接口

宽相检测为每个候选对$(i,j) \in \mathcal{P}$触发窄相检测。窄相检测需要计算:

  1. 碰撞状态:是否发生碰撞或穿透
  2. 接触点:世界坐标系下的接触位置$\mathbf{p}_c$
  3. 接触法向:从物体$j$指向物体$i$的单位法向$\mathbf{n}$
  4. 穿透深度:沿法向的穿透距离$d$(负值表示分离)

这些信息构成接触流形(Contact Manifold),是后续接触力计算的输入。

4.2 GJK与EPA算法

Gilbert-Johnson-Keerthi(GJK)算法是凸体间距离计算的经典方法,其优雅之处在于将碰撞检测问题转化为闵可夫斯基差中的原点包含测试。Expanding Polytope Algorithm(EPA)则在GJK检测到碰撞后,计算穿透深度和接触法向。

4.2.1 闵可夫斯基差与支撑映射

给定两个凸体$A$和$B$,其闵可夫斯基差定义为: $$A \ominus B = \{a - b | a \in A, b \in B\}$$ 关键观察:两凸体相交当且仅当其闵可夫斯基差包含原点

直接构造闵可夫斯基差是不现实的,GJK通过支撑映射(Support Mapping)隐式地在闵可夫斯基差上操作: $$s_{A \ominus B}(\mathbf{d}) = s_A(\mathbf{d}) - s_B(-\mathbf{d})$$ 其中支撑函数$s_A(\mathbf{d})$返回凸体$A$在方向$\mathbf{d}$上的最远点: $$s_A(\mathbf{d}) = \arg\max_{\mathbf{p} \in A} \mathbf{d} \cdot \mathbf{p}$$

4.2.2 GJK算法核心

GJK维护一个单纯形(Simplex)$W \subseteq A \ominus B$,迭代地逼近原点。算法框架:

  1. 初始化:选择任意方向$\mathbf{d}$,计算$w_1 = s_{A \ominus B}(\mathbf{d})$,初始化$W = \{w_1\}$

  2. 迭代: - 计算$W$中距原点最近的点$\mathbf{v}$及其所在的最小面 - 若$|\mathbf{v}| < \epsilon$,原点在$W$内,返回碰撞 - 计算搜索方向$\mathbf{d} = -\mathbf{v}$ - 计算新支撑点$w_{new} = s_{A \ominus B}(\mathbf{d})$ - 若$w_{new} \cdot \mathbf{d} \leq \mathbf{v} \cdot \mathbf{d} + \epsilon$,返回分离,距离为$|\mathbf{v}|$ - 更新$W = \text{ConvexHull}(W \cup \{w_{new}\})$中包含$\mathbf{v}$的最小面

  3. 终止:算法保证在有限步内收敛(凸体顶点数有限)

4.2.3 单纯形更新策略

GJK的效率关键在于高效的单纯形更新。在3D中,单纯形可能是:

  • (0-单纯形):1个顶点
  • 线段(1-单纯形):2个顶点
  • 三角形(2-单纯形):3个顶点
  • 四面体(3-单纯形):4个顶点

Voronoi区域分析决定如何更新单纯形。例如,对于线段$[A,B]$:

        Voronoi(A)  |  Voronoi(AB)  |  Voronoi(B)
             A------+-------AB-------+------B
                    |                |

原点在哪个区域决定了最近点位置和下一步搜索方向。

4.2.4 EPA算法

当GJK检测到碰撞(原点在闵可夫斯基差内)时,EPA计算穿透信息。EPA维护一个包含原点的凸多面体,通过迭代扩展找到距原点最近的边界点。

算法步骤

  1. 初始化:使用GJK的最终四面体作为初始多面体

  2. 迭代: - 找到距原点最近的面$f_{min}$,距离为$d_{min}$ - 计算该面法向$\mathbf{n} = \text{normal}(f_{min})$ - 获取支撑点$w = s_{A \ominus B}(\mathbf{n})$ - 若$w \cdot \mathbf{n} - d_{min} < \epsilon$,返回穿透深度$d_{min}$和法向$\mathbf{n}$ - 将$w$加入多面体,更新面列表

  3. 面管理:维护凸包性质,删除从$w$可见的面,添加新面

EPA的计算代价高于GJK,但提供了精确的穿透信息,对于精确接触响应至关重要。

4.3 连续碰撞检测(CCD)

离散碰撞检测在每个时间步检查静态配置下的碰撞,可能错过高速运动物体间的碰撞(隧穿问题)。连续碰撞检测通过考虑物体在时间区间$[t, t+\Delta t]$内的完整运动轨迹来解决这一问题。

4.3.1 运动物体的时空表示

对于在时间区间$[0,1]$内运动的物体,其配置由插值函数描述:

线性运动(最常用): $$\mathbf{x}(s) = (1-s)\mathbf{x}_0 + s\mathbf{x}_1, \quad s \in [0,1]$$ 螺旋运动(刚体): $$\mathbf{T}(s) = \exp(s\log(\mathbf{T}_1\mathbf{T}_0^{-1}))\mathbf{T}_0$$ 其中$\mathbf{T} \in SE(3)$为刚体变换矩阵。

4.3.2 分离轴定理的时间扩展

对于凸多面体,分离轴定理(SAT)可扩展到CCD。两个运动凸体在$[0,1]$内不碰撞,当且仅当存在分离轴$\mathbf{n}$使得: $$\forall s \in [0,1]: \max_{\mathbf{p} \in A(s)} \mathbf{n} \cdot \mathbf{p} < \min_{\mathbf{q} \in B(s)} \mathbf{n} \cdot \mathbf{q}$$ 潜在分离轴包括:

  • 面法向:$\mathbf{n}_f(s)$
  • 边叉积:$\mathbf{e}_A(s) \times \mathbf{e}_B(s)$

每个潜在分离轴对应一个时间区间多项式,需要求解其根。

4.3.3 保守推进算法

保守推进(Conservative Advancement)通过迭代缩小安全时间步长来逼近首次碰撞时刻:

  1. 初始化:$t_{safe} = 0$,$t_{max} = 1$

  2. 迭代: - 计算$t_{mid} = (t_{safe} + t_{max})/2$时刻的分离距离$d$ - 基于速度界估计安全推进步长: $$\Delta t_{safe} = \frac{d}{v_{max}}$$

  • 更新$t_{safe} = \min(t_{mid} + \Delta t_{safe}, t_{max})$
  • 若$t_{safe} \geq t_{max} - \epsilon$,返回无碰撞
  • 若$d < \epsilon$,返回碰撞时刻$t_{safe}$

速度界$v_{max}$的紧致估计对算法效率至关重要。

4.3.4 射线投射方法

对于某些几何形状,CCD可转化为射线投射问题。考虑点$\mathbf{p}$与三角形$\triangle ABC$的CCD:

点的轨迹:$\mathbf{p}(t) = \mathbf{p}_0 + t\mathbf{v}_p$

三角形平面方程:$\mathbf{n}(t) \cdot (\mathbf{x} - \mathbf{a}(t)) = 0$

碰撞条件转化为求解: $$\mathbf{n}(t) \cdot (\mathbf{p}(t) - \mathbf{a}(t)) = 0$$ 这是一个关于$t$的三次方程(线性运动情况下),可用Cardano公式求解。

4.4 空间划分数据结构

高效的空间数据结构是宽相检测性能的关键。不同的数据结构适用于不同的场景特征,包括物体分布、运动模式和查询类型。本节介绍四种主要的空间划分策略及其在物理仿真中的应用。

4.4.1 层次包围体(BVH)

BVH通过树状结构组织包围体,支持高效的射线投射和最近邻查询。每个节点存储其子树所有物体的包围体。

构建策略

  1. 自顶向下(Top-down):递归划分物体集合 - SAH(Surface Area Heuristic):最小化遍历代价 $$C_{SAH} = C_{traverse} + \frac{A_{left}}{A_{parent}}C_{left} + \frac{A_{right}}{A_{parent}}C_{right}$$
  • 中值划分:沿最长轴均分
  1. 自底向上(Bottom-up):迭代合并最近节点对

  2. 混合方法:LBVH(Linear BVH)使用Morton码快速构建

动态更新

对于运动物体,BVH需要增量更新:

  • 重拟合(Refit):自底向上更新包围体,$O(n)$复杂度
  • 重构(Rebuild):检测退化子树并局部重建
  • 插入/删除:维护树平衡的同时局部调整

双重遍历(Dual Traversal)用于BVH间的碰撞检测:

function DualTraverse(nodeA, nodeB):
    if not Overlap(nodeA.bbox, nodeB.bbox):
        return
    if IsLeaf(nodeA) and IsLeaf(nodeB):
        NarrowPhase(nodeA.object, nodeB.object)
    else if DescendA(nodeA, nodeB):
        DualTraverse(nodeA.left, nodeB)
        DualTraverse(nodeA.right, nodeB)
    else:
        DualTraverse(nodeA, nodeB.left)
        DualTraverse(nodeA, nodeB.right)

4.4.2 八叉树(Octree)

八叉树通过递归八分空间来组织物体,特别适合均匀分布的场景。

自适应细分

节点细分条件:

  • 物体数量超过阈值$n_{max}$
  • 深度未达到$d_{max}$
  • 物体大小差异超过阈值

松散八叉树(Loose Octree):

每个节点的实际包围体大于其理论边界: $$\text{LooseBound} = \text{Bound} \times (1 + k)$$ 其中$k \in [0.5, 1.0]$为松散因子。这允许略大于节点的物体无需跨越多个节点。

查询加速

  1. 视锥剔除:快速排除视锥外的子树
  2. 层次深度:为不同大小的查询选择合适的起始深度
  3. 邻居缓存:缓存相邻节点指针加速局部查询

4.4.3 k-d树

k-d树通过轴对齐平面递归划分空间,每层循环使用不同的划分轴。

SAH构建

选择最优划分位置$p$和轴$a$: $$\arg\min_{a,p} \left[ n_L(p) \cdot V_L(p) + n_R(p) \cdot V_R(p) \right]$$ 其中$n_L, n_R$为左右子空间的物体数,$V_L, V_R$为体积。

平衡性维护

  • 标准k-d树:严格交替轴向
  • 优化k-d树:每层选择方差最大的轴
  • BBF(Best Bin First):维护优先队列进行近似最近邻搜索

4.4.4 均匀网格

均匀网格是最简单但往往最有效的空间划分,特别适合GPU并行化。

哈希网格

使用空间哈希避免存储空单元: $$h(i,j,k) = (i \cdot p_1 \oplus j \cdot p_2 \oplus k \cdot p_3) \bmod M$$ 分层网格

多分辨率网格处理不同尺度的物体:

  • 第$l$层单元大小:$s_l = s_0 \cdot 2^l$
  • 物体分配到层:$l = \lceil \log_2(r/s_0) \rceil$

GPU优化

  1. 紧凑存储:使用数组而非链表
  2. 原子操作:并行插入时的冲突处理
  3. Z-order:改善缓存局部性

4.4.5 性能比较与选择指南

| 数据结构 | 构建复杂度 | 查询复杂度 | 内存占用 | 适用场景 |

数据结构 构建复杂度 查询复杂度 内存占用 适用场景
BVH $O(n\log n)$ $O(\log n)$ $O(n)$ 不均匀分布、复杂几何
八叉树 $O(n\log n)$ $O(\log n)$ $O(n)$ 均匀分布、体素数据
k-d树 $O(n\log n)$ $O(\log n)$ $O(n)$ 点云、最近邻查询
均匀网格 $O(n)$ $O(1)$ $O(n+m^3)$ 粒子系统、GPU并行

选择建议:

  • 机器人仿真:BVH(复杂几何体)
  • 流体仿真:均匀网格(大量粒子)
  • 大场景:分层混合结构
  • 动态场景:考虑更新代价

4.5 强化学习应用:高效碰撞检测对训练速度的优化

在强化学习的物理仿真环境中,碰撞检测往往占据30-50%的计算时间。优化碰撞检测不仅提升单个环境的仿真速度,更重要的是支持大规模并行训练。本节探讨碰撞检测在RL训练中的特殊考虑。

4.5.1 RL环境的碰撞检测特征

强化学习训练环境具有独特的碰撞检测需求:

  1. 高频重置:环境频繁重置到初始状态,空间数据结构需要快速重建
  2. 批量仿真:数千个环境并行运行,内存布局影响缓存性能
  3. 确定性要求:相同输入必须产生相同输出,影响并行化策略
  4. 简化几何:训练环境通常使用简化碰撞体提升速度

4.5.2 并行环境的碰撞检测架构

环境级并行

每个环境维护独立的碰撞检测状态:

class ParallelCollisionWorld:
    def __init__(self, num_envs):
        self.broad_phases = [BroadPhase() for _ in range(num_envs)]
        self.contact_caches = [ContactCache() for _ in range(num_envs)]

    def step(self, env_id):
        # 独立处理每个环境无需同步
        pairs = self.broad_phases[env_id].find_pairs()
        contacts = narrow_phase(pairs)
        self.contact_caches[env_id].update(contacts)

GPU批处理

将多个环境的碰撞检测合并为单个GPU kernel:

__global__ void batch_gjk_kernel(
    ConvexShape* shapes,  // [num_envs, max_objects, ...]
    Transform* transforms, // [num_envs, max_objects, ...]
    Contact* contacts      // [num_envs, max_contacts, ...]
) {
    int env_id = blockIdx.x;
    int pair_id = threadIdx.x;
    // 每个线程处理一个碰撞对
    gjk_test(shapes[env_id], transforms[env_id], 
             pair_id, contacts[env_id]);
}

4.5.3 接触缓存与暖启动

连续帧之间的接触通常具有时间相关性。接触缓存可显著加速窄相检测:

接触缓存策略

  1. 持久化接触ID:为每个接触对分配唯一ID
  2. 暖启动GJK:使用上一帧的单纯形初始化
  3. 接触流形更新:增量更新接触点而非完全重算

RL特殊考虑

  • 环境重置时清空缓存
  • 域随机化可能破坏缓存有效性
  • 需要平衡缓存收益与内存开销

4.5.4 简化碰撞几何

训练时使用简化碰撞体,部署时使用精确几何:

简化策略

  1. 凸包近似:将凹体分解为凸体集合
  2. 形状原语:用盒、球、胶囊体近似复杂形状
  3. LOD系统:根据物体重要性选择细节级别

Sim-to-Real考虑

class CollisionGeometry:
    def __init__(self):
        self.training_shape = SimplifiedConvexHull()
        self.deployment_shape = DetailedMesh()
        self.shape_noise = Normal(0, 0.02)  # 域随机化

    def get_shape(self, mode='train'):
        if mode == 'train':
            # 训练时添加随机扰动补偿简化误差
            return self.training_shape.perturb(self.shape_noise)
        else:
            return self.deployment_shape

4.5.5 性能指标与分析

关键性能指标

  1. 吞吐量:每秒处理的碰撞对数
  2. 延迟:单次碰撞检测的时间
  3. 扩展性:性能随环境数量的变化
  4. 内存带宽:数据传输效率

性能剖析示例(1000个并行环境,每环境20个物体):

| 组件 | CPU时间(ms) | GPU时间(ms) | 占比 |

组件 CPU时间(ms) GPU时间(ms) 占比
宽相检测 2.3 0.8 15%
窄相检测 5.1 1.2 35%
接触求解 4.8 1.5 32%
数据传输 - 0.6 8%
其他 1.2 0.4 10%

4.5.6 碰撞检测的梯度计算

可微分仿真需要碰撞检测的梯度:

光滑接触模型

将硬接触替换为光滑势能: $$\phi(d) = \begin{cases} \frac{k}{2}d^2 & d < 0 \text{(穿透)} \\ 0 & d \geq 0 \text{(分离)} \end{cases}$$ 梯度: $$\nabla_{\mathbf{x}} \phi = k \cdot d \cdot \nabla_{\mathbf{x}} d$$ 其中$\nabla_{\mathbf{x}} d$是距离函数的梯度,可从GJK的最近点计算。

梯度截断

防止梯度爆炸: $$\tilde{\nabla} = \text{clip}(\nabla, -g_{max}, g_{max})$$

4.6 历史人物:Gilbert与1988年GJK算法的诞生

Elmer G. Gilbert,密歇根大学航空航天工程系教授,与他的同事Daniel W. Johnson和S. Sathiya Keerthi在1988年发表了划时代的论文"A Fast Procedure for Computing the Distance Between Complex Objects in Three-Dimensional Space",提出了后来被称为GJK的算法。这个算法的诞生不仅是计算几何的里程碑,更体现了跨学科研究的力量。

4.6.1 算法诞生的背景

1980年代,机器人路径规划和计算机图形学的快速发展对碰撞检测提出了迫切需求。当时的方法要么只适用于特定形状(如球体),要么计算复杂度过高。Gilbert的团队原本在研究控制系统的可达集计算,意外发现其中的凸优化技术可以应用于距离计算。

4.6.2 关键创新

GJK算法的革命性在于三个关键洞察:

  1. 闵可夫斯基差的运用:将两个凸体的距离问题转化为单个凸体与原点的距离问题
  2. 支撑映射的抽象:无需显式构造闵可夫斯基差,只需要支撑函数
  3. 迭代收缩:通过维护和更新单纯形快速收敛到最优解

Gilbert回忆道:"最困难的部分不是算法本身,而是证明它在所有退化情况下都能正确工作。我们花了近一年时间处理各种特殊情况。"

4.6.3 算法的演进

原始GJK算法主要用于计算凸体间的最小距离。后续研究者的贡献包括:

  • Cameron (1997):将GJK扩展到swept volume的碰撞检测
  • van den Bergen (1999):优化了单纯形更新策略,提升性能30%
  • Montanari et al. (2017):提出GJK 2.0,改进数值稳定性

4.6.4 工业影响

GJK算法被广泛应用于:

  • 游戏引擎:Unity、Unreal Engine的物理系统
  • 机器人仿真:Gazebo、MuJoCo的碰撞检测
  • CAD系统:装配验证和间隙分析
  • 航天工程:卫星对接和空间碎片规避

4.6.5 Gilbert的研究哲学

Gilbert强调理论与实践的结合:"一个好的算法必须在数学上优雅,在工程上实用。GJK的成功在于它同时满足了这两个标准。"他的研究方法论影响了一代计算几何研究者:

  1. 从应用问题出发:真实需求驱动理论创新
  2. 跨学科思维:控制理论技术解决几何问题
  3. 严格的正确性证明:处理所有边界情况
  4. 实现的简洁性:核心算法仅需百行代码

4.7 高级话题:非凸形状的隐式表示与SDF碰撞检测

符号距离场(Signed Distance Field, SDF)提供了一种统一处理凸体和非凸体碰撞的优雅方法。通过将几何体表示为空间中每点到表面最近距离的标量场,SDF不仅简化了碰撞检测,还自然支持形状混合、变形和可微分计算。

4.7.1 SDF的数学基础

给定闭合表面$\mathcal{S}$,其符号距离函数定义为: $$\phi(\mathbf{x}) = \begin{cases} -\min_{\mathbf{y} \in \mathcal{S}} |\mathbf{x} - \mathbf{y}| & \mathbf{x} \text{ 在内部} \\ 0 & \mathbf{x} \in \mathcal{S} \\ \min_{\mathbf{y} \in \mathcal{S}} |\mathbf{x} - \mathbf{y}| & \mathbf{x} \text{ 在外部} \end{cases}$$ SDF的关键性质:

  • Eikonal方程:$|\nabla \phi| = 1$(几乎处处成立)
  • 最近点投影:$\mathbf{p}_{closest} = \mathbf{x} - \phi(\mathbf{x}) \nabla \phi(\mathbf{x})$
  • 法向提取:表面法向$\mathbf{n} = \nabla \phi / |\nabla \phi|$

4.7.2 SDF的构造与存储

解析SDF

基本形状的解析表达式:

  • 球:$\phi(\mathbf{x}) = |\mathbf{x} - \mathbf{c}| - r$
  • 盒:$\phi(\mathbf{x}) = |\max(\mathbf{0}, |\mathbf{x}| - \mathbf{b})|$
  • 圆环:$\phi(\mathbf{x}) = |[|\mathbf{x}_{xy}| - R, x_z]| - r$

离散SDF

复杂形状通过体素网格存储:

class SDFGrid:
    def __init__(self, resolution, bounds):
        self.grid = np.zeros((resolution,)*3)
        self.voxel_size = (bounds[1] - bounds[0]) / resolution

    def sample(self, pos):
        # 三线性插值
        idx = (pos - self.bounds[0]) / self.voxel_size
        return trilinear_interpolate(self.grid, idx)

神经SDF

使用神经网络隐式表示: $$\phi(\mathbf{x}) = f_\theta(\gamma(\mathbf{x}))$$ 其中$\gamma$是位置编码,$f_\theta$是MLP网络。

4.7.3 基于SDF的碰撞检测

点与SDF碰撞

直接评估:碰撞当且仅当$\phi(\mathbf{p}) \leq 0$

球与SDF碰撞

球心$\mathbf{c}$,半径$r$的球与SDF碰撞: $$\phi(\mathbf{c}) \leq r$$ 穿透深度:$d = r - \phi(\mathbf{c})$

SDF间碰撞

两个SDF $\phi_A$和$\phi_B$的碰撞检测通过梯度下降寻找: $$\min_{\mathbf{x}} [\phi_A(\mathbf{x}) + \phi_B(\mathbf{x})]$$ 若最小值小于0,则发生碰撞。

4.7.4 SDF的布尔运算

SDF自然支持CSG(Constructive Solid Geometry)运算:

  • 并集:$\phi_{A \cup B} = \min(\phi_A, \phi_B)$
  • 交集:$\phi_{A \cap B} = \max(\phi_A, \phi_B)$
  • 差集:$\phi_{A \setminus B} = \max(\phi_A, -\phi_B)$

光滑版本(避免导数不连续): $$\phi_{smooth} = -\frac{1}{k}\log(e^{-k\phi_A} + e^{-k\phi_B})$$

4.7.5 SDF在软体仿真中的应用

SDF特别适合处理变形体的碰撞:

动态SDF更新

变形时的SDF更新策略:

  1. 局部更新:只更新变形区域附近的SDF值
  2. 层次细化:使用八叉树自适应细化
  3. 时间插值:$\phi(t) = (1-t)\phi_0 + t\phi_1$

接触力计算

基于SDF的罚函数: $$\mathbf{f} = \begin{cases} -k_n \phi \nabla\phi - k_d \mathbf{v}_n & \phi < 0 \\ \mathbf{0} & \phi \geq 0 \end{cases}$$

4.7.6 SDF的梯度计算与优化

解析梯度

对于解析SDF,梯度可直接计算。对于离散SDF: $$\nabla\phi \approx \frac{1}{2h}[\phi(x+h) - \phi(x-h), ...]$$ 自动微分

现代框架(JAX、PyTorch)支持SDF操作的自动微分:

def sdf_loss(params, points, targets):
    pred = neural_sdf(params, points)
    return jnp.mean((pred - targets)**2)

grad_fn = jax.grad(sdf_loss)

优化应用

  • 形状优化:最小化SDF定义的目标函数
  • 轨迹优化:避免SDF定义的障碍物
  • 抓取规划:最大化与物体SDF的接触

本章小结

本章深入探讨了碰撞检测的核心算法和实现技术。我们从两阶段检测策略开始,理解了宽相筛选和窄相精确检测的分工。GJK算法展示了如何通过闵可夫斯基差的几何洞察优雅地解决凸体碰撞问题,而EPA算法则提供了精确的穿透信息。连续碰撞检测解决了高速运动的隧穿问题,空间数据结构的选择直接影响了大规模场景的性能。

关键概念回顾

  1. 宽相检测:使用AABB、扫描线或空间哈希快速筛选候选碰撞对,将$O(n^2)$复杂度降至$O(n\log n)$或更低

  2. GJK算法: - 核心思想:$A \cap B \neq \emptyset \Leftrightarrow \mathbf{0} \in A \ominus B$ - 支撑映射:$s_{A \ominus B}(\mathbf{d}) = s_A(\mathbf{d}) - s_B(-\mathbf{d})$ - 迭代收敛保证有限步终止

  3. 连续碰撞检测: - 分离轴定理的时间扩展 - 保守推进算法:$\Delta t_{safe} = d / v_{max}$ - 射线投射方法处理点-三角形CCD

  4. 空间数据结构: - BVH:适合复杂几何,支持动态更新 - 八叉树:均匀分布场景的理想选择 - 均匀网格:GPU并行的最佳方案

  5. SDF表示: - 统一处理凸体和非凸体 - 自然支持布尔运算和形状混合 - 可微分性支持基于梯度的优化

核心公式汇总

| 概念 | 公式 | 说明 |

概念 公式 说明
AABB重叠 $\max(a_{min}, b_{min}) \leq \min(a_{max}, b_{max})$ 各轴独立检测
闵可夫斯基差 $A \ominus B = \{a - b | a \in A, b \in B\}$ 碰撞检测的几何基础
支撑函数 $s_A(\mathbf{d}) = \arg\max_{\mathbf{p} \in A} \mathbf{d} \cdot \mathbf{p}$ GJK的核心操作
SDF定义 $\phi(\mathbf{x}) = \pm \min_{\mathbf{y} \in \mathcal{S}} |\mathbf{x} - \mathbf{y}|$ 隐式几何表示
CCD条件 $\exists t \in [0,1]: A(t) \cap B(t) \neq \emptyset$ 连续碰撞的数学表述

性能考量

在实际应用中,碰撞检测的性能优化需要考虑:

  1. 算法选择:根据几何复杂度和精度需求选择合适的算法
  2. 数据结构:匹配场景特征(均匀性、动态性、规模)
  3. 并行化:利用GPU或多核CPU加速批量检测
  4. 缓存策略:利用时间相关性减少重复计算
  5. 精度权衡:训练时可用简化几何,部署时用精确模型

练习题

练习 4.1:GJK算法实现

实现一个2D版本的GJK算法,用于检测两个凸多边形是否相交。要求处理退化情况(如共线点)。

提示:先实现支撑函数,然后实现单纯形更新逻辑。注意处理数值精度问题。

参考答案

GJK的2D实现关键步骤:

  1. 支撑函数对于多边形:遍历所有顶点,找到与方向点积最大的顶点
  2. 单纯形在2D中最多包含3个点(三角形)
  3. Voronoi区域分析: - 点:直接返回该点到原点的方向 - 线段:检查原点投影是否在线段内 - 三角形:检查原点是否在内部
  4. 数值稳定性:使用$\epsilon = 10^{-6}$判断零值
  5. 最多迭代次数:顶点数之和,保证终止

关键优化:缓存上一帧的单纯形作为初始值,可将平均迭代次数从5-7次降至2-3次。

练习 4.2:扫描线算法优化

给定1000个运动的AABB,设计一个增量更新的扫描线算法。分析最坏情况和平均情况的时间复杂度。

提示:考虑如何维护有序性,以及如何处理跨越多个其他物体的大幅移动。

参考答案

增量扫描线算法设计:

  1. 数据结构: - 三个有序数组存储区间端点 - 哈希表存储活动区间集合 - 脏标记追踪需要更新的物体

  2. 增量更新策略: - 小移动($\Delta x < \text{threshold}$):冒泡排序局部调整,$O(k)$,$k$为跨越的端点数 - 大移动:删除后重新插入,$O(\log n)$

  3. 复杂度分析: - 最坏情况:所有物体大幅移动,$O(n\log n)$ - 平均情况(连贯运动):$O(n)$,每个物体平均跨越常数个端点 - 空间复杂度:$O(n)$

  4. 优化技巧: - 睡眠物体不参与更新 - 时间戳避免重复的重叠测试 - SIMD加速区间重叠测试

实验表明,对于典型的物理仿真场景,增量更新比完全重建快5-10倍。

练习 4.3:CCD的数值稳定性

推导线性运动物体的CCD中,数值误差对碰撞时间计算的影响。给出保证数值稳定的实现建议。

提示:分析浮点运算中的舍入误差传播,特别是在接近平行运动的情况。

参考答案

数值误差分析:

  1. 误差来源: - 浮点舍入:每次运算引入$\epsilon_{machine} \approx 10^{-16}$相对误差 - 条件数:当相对速度接近0时,时间计算的条件数趋于无穷

  2. 误差传播: 设相对速度$v_{rel} = v_A - v_B$,位置差$\Delta p = p_A - p_B$ 碰撞时间$t = -\Delta p / v_{rel}$

相对误差:$\frac{\delta t}{t} \approx \frac{\delta(\Delta p)}{|\Delta p|} + \frac{\delta v_{rel}}{|v_{rel}|}$

当$|v_{rel}| \to 0$时,误差放大因子$\to \infty$

  1. 稳定化策略: - 速度阈值:$|v_{rel}| < \epsilon_{vel}$时视为静止接触 - 区间算术:使用区间$[t_{min}, t_{max}]$而非点估计 - 自适应精度:根据速度大小选择容差 - 后验验证:计算得到$t$后,验证$|p_A(t) - p_B(t)| < \epsilon_{pos}$

  2. 实现建议

if abs(v_rel) < 1e-8 * max(|v_A|, |v_B|):
    return NO_COLLISION  # 视为平行运动
t = -dot(delta_p, v_rel) / dot(v_rel, v_rel)
# 保守区间
t_min = t * (1 - 1e-6)
t_max = t * (1 + 1e-6)

练习 4.4:BVH的SAH构建

证明Surface Area Heuristic(SAH)在某些假设下是最优的划分策略。设计一个快速近似SAH的方法。

提示:考虑射线遍历的期望代价模型。

参考答案

SAH最优性证明:

  1. 假设: - 射线均匀分布 - 射线与包围盒相交的概率正比于表面积 - 遍历代价远小于相交测试代价

  2. 期望代价模型: $$C = P_{left} \cdot C_{left} + P_{right} \cdot C_{right} + C_{traverse}$$ 其中$P_{left} = A_{left}/A_{parent}$是访问左子树的概率

  3. 最优性:在上述假设下,SAH给出期望代价最小的划分

  4. 快速近似方法: - 分箱法:将空间划分为$k$个箱(通常$k=16$),只在箱边界处评估SAH - 线性BVH:使用Morton码排序,$O(n)$构建时间 - 增量SAH:利用$A_{box}(S \cup \{p\}) - A_{box}(S)$的增量计算

  5. 实践优化

// 分箱SAH
for(int axis = 0; axis < 3; axis++) {
    float best_cost = INF;
    for(int i = 1; i < NUM_BINS; i++) {
        float cost = compute_sah(bins[axis][i]);
        if(cost < best_cost) {
            best_cost = cost;
            best_split = {axis, bins[axis][i].position};
        }
    }
}

实验表明,16-分箱SAH的构建速度比精确SAH快10倍,而遍历性能仅降低5%。

练习 4.5:SDF的解析梯度

推导胶囊体(Capsule)的SDF及其梯度的解析表达式。讨论梯度在不同区域的连续性。

提示:胶囊体可视为线段沿法向膨胀得到。

参考答案

胶囊体SDF推导:

  1. 几何定义:胶囊体是线段$[A, B]$沿各方向膨胀半径$r$得到

  2. SDF表达式

点P到线段的最近点:
t = clamp(dot(P-A, B-A) / dot(B-A, B-A), 0, 1)
Q = A + t*(B-A)

SDF:
φ(P) = |P - Q| - r
  1. 梯度计算: $$\nabla\phi = \frac{P - Q}{|P - Q|}$$ 特殊情况:
  • $P = Q$时梯度不存在(实践中返回任意单位向量)
  1. 连续性分析: - 区域1($t \in (0,1)$):$C^\infty$连续 - 区域2($t = 0$或$t = 1$):$C^0$连续,梯度在端点球面处不连续 - 边界:从圆柱到半球的过渡处,梯度方向突变

  2. 数值处理

vec3 capsule_gradient(vec3 p, vec3 a, vec3 b, float r) {
    vec3 pa = p - a, ba = b - a;
    float t = clamp(dot(pa, ba) / dot(ba, ba), 0, 1);
    vec3 q = a + t * ba;
    vec3 pq = p - q;
    float d = length(pq);
    if(d < 1e-6) return vec3(1, 0, 0); // 任意方向
    return pq / d;
}

梯度不连续性在实践中通过添加小的平滑半径缓解。

练习 4.6:并行碰撞检测的负载均衡

设计一个GPU上的并行碰撞检测算法,处理不均匀的碰撞对分布。分析负载均衡策略。

提示:考虑工作窃取和动态任务分配。

参考答案

GPU并行碰撞检测设计:

  1. 问题分析: - 碰撞对分布不均:某些区域密集,某些稀疏 - GPU的SIMT模型:分支导致线程发散 - 内存访问模式:合并访问提升带宽利用

  2. 两级并行策略

// 第一级:网格并行
__global__ void broad_phase(Grid* grid, PairList* pairs) {
    int cell_id = blockIdx.x;
    // 每个块处理一个网格单元
    process_cell(grid[cell_id], pairs);
}

// 第二级:碰撞对并行
__global__ void narrow_phase(PairList* pairs, Contacts* contacts) {
    int pair_id = blockIdx.x * blockDim.x + threadIdx.x;
    if(pair_id < pairs.count) {
        gjk_test(pairs[pair_id], contacts);
    }
}
  1. 负载均衡技术: - 持久化线程:线程主动获取任务直到队列空 - 工作分块:将大任务分解为固定大小的块 - 双缓冲:一个缓冲区写入,另一个处理

  2. 优化策略: - 紧凑表示减少内存传输 - 共享内存缓存频繁访问的数据 - 原子操作管理输出缓冲区

  3. 性能分析: - 均匀分布:接近线性加速 - 极度不均:性能降至串行的20-30% - 实际场景:达到理论峰值的60-80%

动态负载均衡比静态分配提升30-50%的性能。

练习 4.7:接触缓存的理论分析

分析接触缓存的命中率与物体运动速度、时间步长的关系。推导最优缓存大小。

提示:建立接触持续时间的概率模型。

参考答案

接触缓存理论分析:

  1. 模型假设: - 物体速度服从Maxwell-Boltzmann分布 - 接触事件服从泊松过程 - 缓存采用LRU替换策略

  2. 接触持续时间: 平均接触持续时间: $$\tau_{contact} = \frac{2r}{v_{rel}}$$ 其中$r$为接触区域特征尺寸,$v_{rel}$为相对速度

  3. 缓存命中率: $$P_{hit} = 1 - e^{-\tau_{contact}/\Delta t}$$ 对于时间步$\Delta t$:

  • $\Delta t \ll \tau_{contact}$:$P_{hit} \approx 1$
  • $\Delta t \gg \tau_{contact}$:$P_{hit} \approx 0$
  1. 最优缓存大小: 设$N$为活动接触对数,$M$为缓存大小

期望收益: $$B(M) = M \cdot P_{hit} \cdot C_{save} - M \cdot C_{memory}$$

最优:$M^* = N \cdot \min(1, \tau_{contact}/(k\Delta t))$

其中$k \approx 3-5$为安全系数

  1. 实践建议: - 典型设置:缓存大小为最大接触数的1.5-2倍 - 自适应策略:根据命中率动态调整 - 分层缓存:热缓存(频繁)+ 冷缓存(偶尔)

实验验证:在机器人抓取任务中,接触缓存可减少70%的窄相检测时间。

练习 4.8:混合碰撞检测策略

设计一个自适应选择不同碰撞检测算法的系统。考虑几何复杂度、运动模式和精度需求。

提示:定义选择准则和切换条件。

参考答案

自适应碰撞检测系统设计:

  1. 算法库: - 球-球:解析公式,$O(1)$ - GJK:凸体通用,$O(n)$迭代 - SAT:多面体,$O(mn)$但可早期退出
    - SDF:复杂形状,$O(1)$查询但预计算代价高

  2. 选择准则

复杂度评分 = vertex_count * face_count
速度评分 = |velocity| / characteristic_length
精度需求 = is_grasping ? HIGH : (is_navigation ? LOW : MEDIUM)

if(复杂度 < 10 && 形状凸):
    use GJK
elif(复杂度 > 100 || 形状非凸):
    use SDF
elif(速度评分 > 5):
    use CCD
else:
    use SAT
  1. 动态切换: - 监控每种算法的实际性能 - 维护滑动窗口统计 - 切换代价模型:避免频繁切换

  2. 缓存一致性: - 统一的接触ID系统 - 算法切换时的状态迁移 - 渐进式切换避免突变

  3. 性能基准: | 场景类型 | 最优算法 | 相对性能 |

场景类型 最优算法 相对性能
简单凸体 GJK 1.0x
复杂网格 SDF 2.3x
高速运动 CCD 1.5x
密集堆叠 SAT+缓存 1.8x

自适应系统比固定算法平均提升40%性能。

常见陷阱与错误

1. GJK算法的数值问题

陷阱:当两个物体几乎接触时,GJK可能因数值误差而无法收敛。

症状

  • 算法在相同的单纯形配置间振荡
  • 返回的距离在小范围内波动
  • 性能突然下降(迭代次数激增)

解决方案

// 使用相对容差而非绝对容差
float tolerance = max(EPSILON_ABS, EPSILON_REL * object_scale);

// 检测振荡
if(simplex == previous_simplex) {
    // 强制终止,返回近似结果
    return approximate_distance;
}

// 限制最大迭代次数
if(iteration > MAX_ITER) {
    // 回退到更稳定但较慢的算法
    return fallback_algorithm();
}

2. CCD的时间步陷阱

陷阱:使用固定的时间步进行CCD可能错过碰撞。

症状

  • 物体穿透,特别是高速或薄物体
  • 碰撞响应发生在错误的位置

正确做法

// 错误:固定步长
for(float t = 0; t < 1; t += 0.1) {
    check_collision(t);
}

// 正确:自适应步长
float t = 0;
while(t < 1) {
    float dt = compute_safe_step(velocity, separation);
    t = min(t + dt, 1.0);
    if(check_collision(t)) {
        // 二分查找精确时刻
        binary_search_collision_time(t - dt, t);
    }
}

3. 空间数据结构的更新延迟

陷阱:延迟更新空间数据结构导致错误的碰撞检测结果。

症状

  • 检测到不存在的碰撞(幽灵碰撞)
  • 错过实际的碰撞
  • 查询结果与实际位置不一致

最佳实践

  • 移动物体后立即标记为脏
  • 查询前批量更新所有脏节点
  • 使用双缓冲避免读写冲突

4. 浮点比较错误

陷阱:直接比较浮点数导致的逻辑错误。

错误示例

// 错误
if(distance == 0) {  // 几乎永远不会为真
    handle_contact();
}

// 正确
if(abs(distance) < EPSILON) {
    handle_contact();
}

5. 坐标系混淆

陷阱:在局部坐标系和世界坐标系之间转换时出错。

常见错误

  • GJK中混用不同坐标系的支撑点
  • 法向量未正确变换
  • 忘记更新变换矩阵

调试技巧

// 添加断言检查坐标系一致性
assert(is_normalized(normal_world));
assert(is_in_world_space(contact_point));

// 可视化调试
debug_draw_arrow(contact_point, normal_world, Color::RED);

6. 并行化的竞态条件

陷阱:多线程碰撞检测中的数据竞争。

症状

  • 不确定的结果
  • 偶发的崩溃
  • 接触信息丢失

解决方案

  • 使用线程局部存储
  • 原子操作管理共享数据
  • 每个线程独立的输出缓冲区

7. 缓存失效问题

陷阱:接触缓存在特定情况下提供错误信息。

触发条件

  • 物体被传送
  • 域随机化改变参数
  • 时间步长突变

缓解措施

class ContactCache {
    void validate() {
        // 检查缓存有效性
        if(teleported || parameters_changed) {
            clear();
            return;
        }

        // 验证接触点仍然有效
        for(auto& contact : cached_contacts) {
            if(separation(contact) > CACHE_TOLERANCE) {
                remove(contact);
            }
        }
    }
};

8. SDF采样不足

陷阱:离散SDF的分辨率不足导致碰撞检测错误。

症状

  • 小物体穿过SDF
  • 接触法向不准确
  • 梯度计算噪声大

解决方案

  • 自适应分辨率:在表面附近提高采样密度
  • 使用三次插值而非线性插值
  • 结合解析SDF处理小特征

最佳实践检查清单

设计阶段

  • [ ] 需求分析
  • 确定精度要求(游戏 vs 工程仿真)
  • 估算场景规模(物体数量、复杂度)
  • 识别性能瓶颈(CPU/GPU、内存、带宽)

  • [ ] 算法选择

  • 简单形状优先使用解析方法
  • 凸体使用GJK/EPA
  • 复杂非凸体考虑SDF或凸分解
  • 高速场景启用CCD

  • [ ] 数据结构设计

  • 静态场景:预计算的BVH
  • 动态场景:增量更新的结构
  • 粒子系统:空间哈希
  • 大规模场景:层次化方法

实现阶段

  • [ ] 数值稳定性
  • 使用相对容差而非绝对容差
  • 处理退化情况(共面、共线、重合)
  • 避免除零和数值下溢

  • [ ] 性能优化

  • 实施宽相筛选减少窄相调用
  • 使用SIMD指令加速向量运算
  • 实现接触缓存和暖启动
  • 考虑LOD和简化几何

  • [ ] 内存管理

  • 预分配缓冲区避免动态分配
  • 使用对象池管理临时数据
  • 注意缓存友好的数据布局

测试阶段

  • [ ] 正确性测试
  • 单元测试基本几何操作
  • 边界情况(接触、分离、穿透)
  • 对称性和传递性测试

  • [ ] 性能测试

  • 基准测试不同场景
  • 压力测试(大量物体)
  • 剖析热点函数

  • [ ] 鲁棒性测试

  • 极端输入(很大/很小的数值)
  • 快速运动和传送
  • 并发访问和竞态条件

调试工具

  • [ ] 可视化
  • 绘制包围体和碰撞对
  • 显示接触点和法向
  • 颜色编码不同检测阶段

  • [ ] 日志和统计

  • 记录碰撞检测时间
  • 统计缓存命中率
  • 追踪算法迭代次数

  • [ ] 验证工具

  • 对比不同算法的结果
  • 检查物理量守恒
  • 验证确定性(相同输入产生相同输出)

优化检查点

  • [ ] 是否充分利用了时间相关性?
  • [ ] 是否有不必要的坐标变换?
  • [ ] 内存访问模式是否优化?
  • [ ] 是否可以批处理相似的操作?
  • [ ] 是否有冗余的碰撞检测?
  • [ ] 数据结构的更新频率是否合适?