C++ 中的快速行进法

引言

快速行进法 (FMM) 是一种计算方法，它在应用 eikonal 方程时已被证明具有巨大的优势。eikonal 方程用于各种涉及波传播的应用，如 计算机视觉、水力学甚至医学成像。例如，Sethian J.A. 在 1996 年提出的一些新颖方法，包括从某个初始曲线以恒定速度向外传播的旅行时间。更具体地说，FMM 是一种近似技术，它模仿图结构中已知的 Dijkstra 最短路径算法，但将其转移到连续表面上。这种情况在火焰传播建模问题中很常见。

在本文中，我们将介绍快速行进法及其深层机制，解释该方法是什么以及如何在 C++ 中实现它。

问题陈述

设 Ω 表示域，设前沿由边界 Γ⊂Ω 定义，并假设存在以某种速度 F(x) 向外的传播，其中速度函数 F(x) 可以是非均匀的。目标是确定 eikonal 类型方程在区域中的每个点 x∈Ω 处的前沿到达时间 T(x)，该方程如下定义：

其中

• T(x) 是前沿在点 x 的到达时间，
• F(x) 是点 x 处的速度函数，
• ∇T(x) 是到达时间的梯度。

然而，如何找到一个逻辑上、快速且最优的数值方法来处理这个方程，尤其是在广阔的域上，是一个问题。

快速行进法 (FMM)

• FMM 是一种有效地求解 eikonal 方程的数值技术。它特别适用于速度函数 F(x) 非负且平滑变化的场合。
• FMM 指的是一种精确求解 eikonal 方程的技术。它在速度函数 F(x) 始终为正且具有连续导数的情况下特别有用。
• FMM 假设前沿以稳定的方式从边界向外辐射，并能够通过使用堆结构优先附加前沿来扩展前沿，从到达时间最小的点开始。这类似于 Dijkstra 最短路径模型，其中不是先最小化节点之间的距离，而是首先快速达到波前表面。

FMM 的过程

1. 初始化： 到达时间 T(x) 在起始边界处初始化为 0，其他所有位置初始化为无穷大。
2. 标记网格点
   
   • 活动点： 这些是已计算过的点。
   • 试探点： 这些是活动点的直接邻域，预计将在其中计算到达时间。
   • 远点： 这些是尚未到达的点；因此，它们的到达时间仍为无穷大。
3. 更新规则： 在每个时间点，选择一个试探点，该试探点的到达时间最接近已计算邻居的试探点，并根据 eikonal 模型调整其到达时间。
4. 终止： 当确定所有点都已处理完到达时间时，算法结束。

程序 1

下面是快速行进法的一个基本 C++ 实现。为简单起见，我们假设网格是均匀的，并且速度函数 F(x)=1 是恒定的。

输出

 
0 1 2 3 4 
1 2 3 4 5 
2 3 4 5 6 
3 4 5 6 7 
4 5 6 7 8   

说明

• 网格表示： 域被结构化为二维网格，其中每个点由到达时间 T(x,y) 表示。最初，所有点的到达时间都定义为无穷大，源点除外，源点为零。
• 优先队列： 使用最小堆（优先队列）来传播前沿。前沿始终按正确的顺序处理，以便首先考虑具有最小值的点的到达时间。
• Eikonal 更新： 在此过程中，通过处理周围邻近点并根据 eikonal 方程更新它们的到达时间来实现。
• 传播： 该方法中的过程扩展到网格上的所有可用点，执行点处理，并在此时，用邻居的到达时间更新每个点。

复杂度分析

• 快速行进法的时间复杂度可以达到 O (Nlog N)，其中 N 指的是网格点的数量。这是因为每个点都处理一次，并且优先队列操作是 O(log N) 的。
• 快速行进法的空间复杂度可以概括为 O(N)，其中 N 是网格点的总数（即，对于二维网格 n×m，或三维网格 n×m×p）。这种线性空间复杂度是 FMM 的主要特点之一。因此，它非常适合解决大规模问题。

程序 2

输出

 
Arrival times of the front at each grid point:
0 INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF INF 
INF INF INF INF INF INF INF INF INF 0   

说明

1. 速度函数： 速度函数 F(x,y) 是一个变量，可以从正弦和余弦的公式中推导出来。可以使用任何其他函数来代替它，以产生更复杂的传播行为。
2. 模块化设计
   
   • 代码被模块化为 FastMarching 类，其中包含算法的基本操作。
   • speedFunction 是独立的，这使得修改网格上每个点的速度变得简单。
3. 由于对 solveEikonal 函数进行了修改，该函数计算相对于不同速度的到达时间，因此提供了一个求解 eikonal 方程的数值方案。
4. eikonal 方程的二次求解器： 该方法通过一种更通用的方法估计到达时间，该方法使用对角线和正交邻居进行计算。solveEikonal 函数包含一个计算例程，该例程实现二次求解器以获得每次试探点的总时间。
5. 多源： 实现支持多个源，在 initializeSources 方法中对其进行初始化。也就是说，一个波阵面可能从多个前沿完全展开。
6. 优先队列： 这确保了前沿从已知到达时间最小的点扩展，如同 Dijkstra 算法一样，因此波前仍然从已知到达时间最小的点快速传播。

复杂度分析

• 时间复杂度： O(NlogN)，其中 N 是网格点的总数（即，对于 2D 网格 N=n×m）。这是因为使用了优先队列（最小堆），其中每个网格点被处理一次，并且堆操作（插入/提取）需要 O(logN)。
• 空间复杂度： O(N)，其中 N 是网格点的总数。

FMM 的应用

• 计算机视觉： 形状建模、形状分割和对象跟踪。
• 医学成像： 植入物，其中波通过不同密度的组织传播。
• 测地线距离计算： 它涉及在地形中寻找最短路径，其中旅行速度取决于深度或高度，并且还受其他因素的控制。
• 机器人技术： 路径规划，使得机器人的速度在所有地形上不是恒定的。
• 流体动力学： 连续速度变化流体中的前沿传播建模。

结论

总之，快速行进法 作为求解 eikonal 方程的算法，可以通过模拟波前传播过程。该方法高效且准确，因此非常适合需要对广阔域中的到达时间进行快速计算的应用。eikonal 方程有许多应用；C++ 实现中提供的就是其中一个例子。它简单且高效。然而，通过更复杂的速度函数或不同的网格几何形状，也可以解决更高级的问题。