Shape Analysis课程笔记09-Geodesic Distances Algorithms

这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍计算离散曲面上测地距离的相关算法。

在上一节课中我们介绍了测地距离的概念，简单来说它是指曲面上任意两点之间最短的距离。

Initial Value Problem on Meshes

我们首先来考虑如何在离散曲面上计算局部最短距离的问题。对于测地线它需要满足切平面上加速度为0的条件，即：

\[\text{proj}_{T_{\gamma(s)}\mathcal{M}} [\gamma''(s)] \equiv 0\]

因此我们希望把这个性质直接推广到离散曲面上。具体来说有三种可能情况：

起点和终点在同一三角形上，此时的测地线即为两点的连线；
起点和终点位于具有公共边的两个相邻三角形上，此时可以把两个三角形展开到同一平面上然后相连得到测地线；
起点和终点位于具有公共边的两个相邻三角形上且展开后它们的连线恰好经过三角形的顶点，此时必须要考虑经过该顶点的其它三角形。

对于第三种情况我们必须要考虑顶点的Gauss曲率。

对于\(K < 0\)的情况可能有多条最短路径，而对于\(K < 0\)的情况展开后平面上的连线却不是最短路径。

总而言之，我们很难在离散曲面上保证连续曲面上测地线的所有性质，即使是很多局部上等价的性质在离散的情况下仍然是有区别的。

实际上计算离散曲面测地线时我们很少会使用ODE的角度来进行处理，很多流行的算法都是使用其他的角度来进行设计的。

Fast Marching

另一种计算测地线的方法是使用全局最短路径的定义。我们对于图上的最短路径有很深入的研究，因此希望可以把这些技术推广到网格上来计算测地线。

Dijkstra’s Algorithm

Dijkstra算法(Dijkstra’s algorithm)是计算图上节点之间最短距离的经典算法，它通过BFS来遍历图上的所有节点并且更新起点到当前节点的最短距离。可以证明Dijkstra算法的复杂度为\(O(\vert E \vert + \vert V \vert \log{\vert V \vert})\)。

Fast Marching

当然Dijkstra算法并不能直接用来计算测地线，我们可以对它进行一定的修改来近似计算测地距离。

Dijkstra算法的问题在于相邻三角形两个顶点之间的最短距离不是通过公共顶点的折线，而是穿过公共边的折线。

人们还发现当三角形距离起点比较远时，测地距离的等值线在会在每个三角面片上趋于直线。

因此我们可以利用这样的性质在每个三角面片上使用直线来近似测地线的等值线。

假设已知三角形上两个顶点\(\mathbf{x}_1\)和\(\mathbf{x}_2\)到起点的测地距离，我们的目标是计算第三个顶点\(\mathbf{x}_3\)的测地距离。记三角形上两个已知顶点的测地距离方程为：

\[d_1 = \mathbf{n}^T \mathbf{x}_1 + c\] \[d_2 = \mathbf{n}^T \mathbf{x}_2 + c\]

其中\(\mathbf{x}_i \in \mathbb{R}^2\)，且\(\mathbf{x}_3 = \mathbf{0}\)。因此可以构造出线性方程组：

\[\mathbf{d} = X^T \mathbf{n} + c \cdot \mathbf{1}\] \[\mathbf{d} = \begin{bmatrix} d_1 \\ d_2 \end{bmatrix}, \ X^T = \begin{bmatrix} \mathbf{x}_1^T \\ \mathbf{x}_2^T \end{bmatrix}\]

通过求解这个线性系统可以得到向量\(\mathbf{n}\)：

\[\mathbf{n} = X^{-T} (\mathbf{d} - c \cdot \mathbf{1})\]

结合向量\(\mathbf{n}\)的约束条件可以得到：

\[\begin{aligned} 1 &= \mathbf{n}^T \mathbf{n} \\ &= (\mathbf{d} - c \cdot \mathbf{1})^T X^{-1} X^{-T} (\mathbf{d} - c \cdot \mathbf{1}) \end{aligned}\]

此时上式可以转换为一个关于\(c\)的二次方程：

\[\big[ \mathbf{1}^T (X^T X)^{-1} \mathbf{1} \big] c^2 + \big[ -2 \mathbf{1}^T (X^T X)^{-1} \mathbf{d} \big] c + \big[ \mathbf{d}^T (X^T X)^{-1} \mathbf{d} \big]= 1\]

这样可以得到\(c\)的两个解：

\[c = \frac{-\bar{b} \pm \sqrt{\bar{b}^2 - 4 \bar{a} (\bar{c}-1)}}{2 \bar{a}}\]

其中，\(\bar{a} = \mathbf{1}^T (X^T X)^{-1} \mathbf{1}\)，\(\bar{b} = -2 \mathbf{1}^T (X^T X)^{-1} \mathbf{d}\)，\(\bar{c} = \mathbf{d}^T (X^T X)^{-1} \mathbf{d}\)。

实际上\(c\)的两个解是具有几何意义的，根据\(\mathbf{n}\)的朝向它表示从起点出发测地距离的等值线经过顶点\(\mathbf{x}_3\)和另外两个顶点的相对顺序。如果我们要求测地距离的等值线先经过\(\mathbf{x}_1\)和\(\mathbf{x}_2\)，则需要选择较大的那个值。