Shape Analysis课程笔记08-Geodesic Distances Intro & Theory

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

Motivation

在接下来的几节课中我们会介绍如何计算曲面上点之间的距离。这里我们需要区分一下intrinsic和extrinsic的概念，intrinsic是指我们只能在曲面(流形)上进行移动，而extrinsic则允许离开曲面。以下图为例，两个猫爪之间的extrinsic距离非常小，而intrinsic距离则很大。

曲面上的intrinsic距离也称为测地距离(geodesic distance)，它可以理解为限制在曲线上移动时两点之间的最短距离。

计算测地距离并不简单。实际上距离函数可能有大量的局部极小值(local minimum)，这容易导致我们在计算时的效率比较低而且会不那么稳定。

测地距离的另一大特点是它有很多相似的定义。这里需要说明的是尽管这些定义在欧式空间下是等价的，但在更一般的情况下往往是有区别的。

计算测地距离的一种方法是把网格理解为一张图，然后使用图上的最短路径算法来求解。当网格非常稠密时顶点之间的最短路径是测地距离的一个很好的近似，但在大多数情况下测地距离都是严格小于最短路径长度的，而且这种近似方式会产生各种各样的问题。

实际上可以严格证明图上的最短路径不会收敛到测地距离上，而且要计算测地距离还需要一些额外的处理。

Defining Geodesic Distance

False Equivalences

测地距离有很多「等价」的定义，但需要注意的是它们之间是由严格区别的。最常见的定义是全局的最短路径，如下图中棕色路线所示；与之对应的是局部最短路径，由黄色路线表示；而另一个常见的定义则是局部最直的路径，由绿色路线表示。

Derivative Formulas

回忆曲线弧长的定义，我们首先利用全局最短曲线来定义测地线。

由于直接计算弧长比较麻烦，更常用的处理方法是使用曲线的能量$E[\gamma ]=\frac{1}{2}{\int }_a^b\Vert {\gamma }^{\prime }(t){\Vert }_2^{2}dt$来代替弧长作为我们的优化函数。实际上能量$E[\gamma ]$是弧长的一个上界，这里简单证明一下。设有实值函数$f,g:[a,b]\to \mathbb{R}$，我们定义函数的内积为：

$$⟨f,g⟩={\int }_a^{b}f(t)g(t)dt$$

根据Cauchy-Schwarz不等式(Cauchy-Schwarz inequality)有：

$$|⟨f,g⟩|\le \Vert f{\Vert }_2\Vert g{\Vert }_2$$ $$\Vert f{\Vert }_2^2=⟨f,f⟩$$

因此有：

$$\begin{array}{rl}2(b-a)E& =(b-a){\int }_a^b\Vert {\gamma }^{\prime }(t){\Vert }_2^{2}dt\\ & =[{\int }_a^{b}1dt][{\int }_a^b\Vert {\gamma }^{\prime }(t){\Vert }_2^{2}dt]\\ & =⟨1,1⟩\cdot ⟨\Vert {\gamma }^{\prime }{\Vert }_2,\Vert {\gamma }^{\prime }{\Vert }_2⟩\\ & =\Vert 1{\Vert }_2^2\cdot \Vert \Vert {\gamma }^{\prime }{\Vert }_2{\Vert }_2^2\\ & \ge ⟨1,\Vert {\gamma }^{\prime }{\Vert }_2{⟩}^2\\ & =[{\int }_a^b\Vert {\gamma }^{\prime }{\Vert }_{2}dt{]}^2\\ & =L^2[\gamma ]\end{array}$$

当且仅当$\Vert {\gamma }^{\prime }{\Vert }_2$为常数时不等式取等号。

接下来我们引入能量$E[\gamma ]$的微分：

$$\begin{array}{rl}\frac{d}{dt}E[{\gamma }_t]& =\frac{d}{dt}\frac{1}{2}{\int }_a^b\Vert {\gamma }_t^{\prime }(s){\Vert }_2^{2}ds\\ & =\frac{1}{2}{\int }_a^b\frac{d}{dt}\Vert {\gamma }_t^{\prime }(s){\Vert }_2^{2}ds\\ & ={\int }_a^b{\gamma }_t^{\prime }(s)\cdot \frac{d{\gamma }_t^{\prime }}{dt}(s)ds\\ & ={\int }_a^b{\gamma }_t^{\prime }(s)\cdot (\frac{d{\gamma }_t}{dt}{)}^{\prime }(s)ds\end{array}$$

利用分部积分有：

$$\begin{array}{rl}\frac{d}{dt}E[{\gamma }_t]& ={\int }_a^b{\gamma }_t^{\prime }(s)\cdot (\frac{d{\gamma }_t}{dt}{)}^{\prime }(s)ds\\ & =[{\gamma }_t^{\prime }\cdot \frac{d{\gamma }_t}{dt}{]}_a^b-{\int }_a^b{\gamma }_t^{″}(s)\cdot \frac{d{\gamma }_t}{dt}(s)ds\\ & =-{\int }_a^b{\gamma }_t^{″}(s)\cdot \frac{d{\gamma }_t}{dt}(s)ds\end{array}$$

注意到曲线${\gamma }_t$必须依附于曲面(流形)$\mathcal{M}$，因此$\frac{d{\gamma }_t}{dt}(s)$一定位于切平面上，即$\frac{d{\gamma }_t}{dt}(s)\in T_{{\gamma }_t(s)}\mathcal{M}$；同时${\gamma }_t^{″}(s)$则与曲线的法向共线。因此上式等价于先对${\gamma }_t^{″}(s)$利用切平面进行分解，然后把它在切平面上的投影与$\frac{d{\gamma }_t}{dt}(s)$进行内积。

$$\begin{array}{rl}\frac{d}{dt}E[{\gamma }_t]& =-{\int }_a^b{\gamma }_t^{″}(s)\cdot \frac{d{\gamma }_t}{dt}(s)ds\\ & =-{\int }_a^b(\frac{d{\gamma }_t}{dt}(s)\cdot {\text{proj}}_{T_{{\gamma }_t(s)}\mathcal{M}}[{\gamma }_t^{″}(s)])ds\end{array}$$

因此，我们可以把曲线${\gamma }_t$看做能量$E[{\gamma }_t]$的一个local minimizer。显然当${\gamma }_t$为测地线时需要满足$\mathrm{\forall }\mathbf{v}:[a,b]\to {\mathbb{R}}^n$：

$${\int }_a^b(\mathbf{v}\cdot {\text{proj}}_{T_{{\gamma }_t(s)}\mathcal{M}}[{\gamma }_t^{″}(s)])ds=0$$

即测地线$\gamma$需要满足ODE：

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

进一步可以证明此时的测地线$\gamma$的导数为常数：

$$\frac{d}{ds}\Vert {\gamma }^{\prime }(s){\Vert }^2=2{\gamma }^{\prime }(s)\cdot {\gamma }^{″}(s)\equiv 0$$

结合能量$E[\gamma ]$和弧长的关系可以得知$\gamma$具有(局部)最短弧长。

类似地可以证明测地线只有出切平面方向的加速度。

Exponential Map

上面的推导说明我们可以把测地线的计算问题转换为求解ODE。根据初始条件不同会有两种解法：给定端点位置或者给定起点位置以及方向。

利用起点位置和方向我们可以定义曲面(流形)上的指数映射(exponential map)，它表示从起点$\mathbf{p}$出发以$\mathbf{v}$作为方向的测地线。

$${\exp }_{\mathbf{p}}(t\mathbf{v})={\gamma }_{\mathbf{p}}(t)$$

这里我们再重申一下，通过求解ODE得到的测地线是局部的，它们可能不是曲面上两点之间的最短路径。

与之对应的概念是cut locus，它表示从$\mathbf{p}$点出发沿测地线进行移动时不再满足全局最短的点集。

Eikonal Equation

本节课最后介绍了eikonal方程(eikonal equation)，它指出测地距离函数需要满足的条件：

$$\Vert \mathrm{\nabla }u(\mathbf{p}){\Vert }_2=1$$

这里需要注意的是梯度$\mathrm{\nabla }u$是定义在切平面上的，即$\mathrm{\nabla }u(\mathbf{p})\in T_{\mathbf{p}}\mathcal{M}$。

Reference

• Lecture 08: Geodesic distances
• Geodesic Distances: Intro & Theory
• Chapter 6: Geodesic distances