DDG课程笔记5-测地线

这个系列是对CMU离散微分几何课程(CMU CS 15-458/858: Discrete Differential Geometry)的总结和回顾。本节我们介绍曲面上测地线与测地距离的相关内容。

Overview

测地线(geodesics)可以看做是直线在弯曲空间中的推广，基于欧式空间中直线的性质我们可以得到测地线的两条基本性质：

straightest: (局部)没有弯曲；
shortest: (局部)距离最短；

同时，测地线还是一个等距变换下的不变量(isometry invariant)。比如说人身体上任意两点之间的测地线不会因为人体的动作而发生改变。

日常生活中最常见的测地线应该是飞机的航线。对于地球上的任意两点之间，它们之间的最短路线是这两点与球心构成的圆弧。如果我们把航线映射到地图上则可以发现此时的测地线并不是两点之间的直线，这充分体现了欧式空间与曲面的差异。

从测地线的两条基本性质出发我们可以对离散曲面上的测地线做出不同的定义。稍后能看到尽管在光滑曲面上这些定义是相互等价的，但在离散后它们之间可能会有一定的差异。

Shortest Path

Dirichlet Energy on Curves

我们首先把测地线定义为两点之间的最短路径。对于平面上的给定两点存在通过它们的参数化曲线\(\gamma(t): [0, 1] \rightarrow \mathbb{R}^2\)，则它对应的Dirichlet能量为：

\[E_D(\gamma) = \int_0^1 \vert \gamma'(t) \vert^2 dt\]

对\(\gamma(t)\)使用弧长参数进行重参数化，重新参数化后的曲线\(\hat{\gamma}\big( c(t) \big)\)关于弧长参数\(c\)的导数(速度)始终为1，因此：

\[E_D(\gamma) = \int_0^1 \big( c'(t) \big)^2 dt\]

重新参数化后最小化Dirichlet能量等价于最小化\(\big( c'(t) \big)^2\)的积分：

\[\min_\gamma E_D(\gamma) = \min_{\hat{\gamma}} \Bigg( \min_c \int_0^1 \big( c'(t) \big)^2 dt \Bigg)\]

此时最小化曲线\(\gamma(t)\)的Dirichlet能量转换为最小化曲线\(c(t)\)的Dirichlet能量。由于参数\(c\)满足条件\(c(0)=0\)，\(c(1)=L\)，因此Dirichlet能量在\(c(t)\)为直线时取最小值，即

\[\min_\gamma E_D(\gamma) = \min_{\hat{\gamma}} \Bigg( \min_c \int_0^1 \big( c'(t) \big)^2 dt \Bigg) = \min_{\hat{\gamma}} L^2\]

上式说明对于平面曲线，最小化Dirichlet能量等价于最小化曲线长度。显然此时最小化Dirichlet能量的曲线为连接端点的直线。

利用Dirichlet能量的性质还可以将Dirichlet能量写成二阶导数内积的形式：

\[\min_\gamma \int_0^1 \vert \gamma'(t) \vert^2 dt = \min_\gamma -\int_0^1 \langle \gamma(t), \gamma''(t) \rangle dt\]

优化目标取极值的条件为二阶导数恒为0，得到一维Laplace方程：

\[\begin{aligned} \frac{\partial^2}{\partial t^2} \gamma(t) &= 0 \\ \gamma(0) &= p \\ \gamma(1) &= q \end{aligned}\]

我们称满足Laplace方程的函数为调和函数(harmonic functions)，同时上面的推导说明测地线为满足边界条件的调和函数。对于二维平面的情况可以验证直线方程为调和函数，因此测地线为连接端点的直线。

将平面曲线推广到曲面上的曲线，可以得到曲面上的Dirichlet能量为：

\[E_D(\gamma) = \int_0^1 \vert \gamma'(t) \vert^2 dt = \int_0^1 g\big( \gamma'(t), \gamma'(t) \big) dt\]

其中\(g\)为曲面的黎曼度量。显然我们可以通过在曲面上求解Laplace方程来获得测地线，但需要注意的是求解Laplace方程得到的测地线仅为Dirichlet能量的极小值点。实际上满足Laplace方程的曲线可能大多数都是局部最小值点或是鞍点，有时甚至可能会出现没有全局最小的情况。

Discrete Shortest Paths

下面考虑离散曲面上的最短路径。一个最直观的想法是假设最短路径在网格的边上，因此可以使用Dijkstra算法来求解网格上的最短路径。

但遗憾的是这个假设是错误的，网格上的最短路径几乎不会完全位于网格的边上，绝大多数情况下最短路径都会位于网格的三角形面上。

而对于顶点\(i\)附近的最短路径，它们的行为则与\(i\)的angle defect(高斯曲率)有关：

当\(\Omega = 0\)时最短路径会穿过顶点\(i\)；

当\(\Omega > 0\)时最短路径在\(i\)附近的边上且不会穿过顶点\(i\)；

当\(\Omega < 0\)时存在多条穿过顶点\(i\)的最短路径；

因此要计算网格上的最短路径需要将Dijkstra算法进行拓展从而考虑最短路径出现在三角面上的情况。目前来看求解网格上的最短路径所需的计算代价仍是非常大的 — \(O(n^2 \log n)\)。

将光滑曲面与离散曲面上的最短路径进行对比还可以发现对于光滑曲面从同一源点出发的任意两条测地线仅在它们存在包含关系时才会相交，而在离散曲面上的最短路径则不满足这个性质。换句话说，使用最短距离定义的测地线在光滑曲面和离散曲面上并不是完全等价的。

Straightest Path

Geometric Perspective

下面我们从另一个角度来认识测地线。我们可以把测地线定义为没有弯曲的的路径，在平面上看这意味着测地线曲率为0，同时参数化后的曲线加速度为0(曲线运动的速度大小和方向保持不变)。接着把这个定义推广到曲面上可以得到曲面上测地线的性质：测地曲率(geodesic curvature)为0，同时协变导数(covariant derivative)为0。

在曲面与曲率一节中我们介绍过曲线的法曲率(\(\kappa_n\))和测地曲率(\(\kappa_g\))的概念，其中法曲率为曲面自身弯曲带来的曲率，而测地曲率则是曲线在切平面上的局部弯曲。一个典型的例子是球面上的测地线：此时测地线为球面上两点通过球心构成的圆弧，在空间中看测地线是显然存在弯曲(法曲率)的，但从局部来看圆弧上任意位置都是平坦的(测地曲率为0)。

Discrete Straightest Geodesics

将平坦的概念推广到离散曲面上，此时测地线可以定义为曲面上转角\(\kappa_i\)为0的曲线。

对于三角形面上的转角，其定义与平面离散曲线的转角相同；对于穿过三角网边的曲线，定义其转角为将两个三角型面展平后曲线的转角；而对于穿过顶点的曲线则需要注意此时没有办法给出转角的定义。

对于穿过顶点的曲线我们可以通过曲线两侧的角度来度量转角的变化，显然测地线在两侧的角度相等(\(\theta_l = \theta_r\))。

Exponential / Logarithmic Map

考虑在曲面上沿测地线进行移动，定义曲面上\(p\)点的指数映射(exponential map) \(\exp_p: T_p M \rightarrow M\)为在\(p\)点沿\(X\)方向测地线移动\(\vert X \vert\)距离移动到\(q\)点：

\[q = \exp_p (X)\]

类似地，我们定义指数映射的逆变换为对数映射(logarithmic map) \(\log_p: M \rightarrow T_p M\)，它表示从\(q\)点出发沿测地线移动到\(p\)点的最短路径为\(X\)：

\[X = \log_p (q)\]

指数映射和对数映射给出了在曲面上进行移动的法则，同时允许我们通过向量\(X\)来描述曲面上点和点之间的关系。利用指数映射和对数映射我们可以定义曲面上一系列给定点的中心为Karcher均值(Karcher mean)，它表示曲面上到给定点集距离平方和最小的点：

我们可以使用迭代的方法来计算Karcher均值：

在曲面上选择一个随机初始化点\(x\)；
利用对数映射计算\(x\)到每个给定点\(y_i\)的距离向量\(v_i = \log_x(y_i)\)；
计算\(v_i\)的平均值\(v = \frac{1}{n} \sum_i \log_x(y_i)\)；
利用指数映射从\(x\)出发沿\(v\)进行移动得到新的位置\(x = \exp_x(v)\)；
回到第2步直到算法收敛。

在大多数情况下只需少次的迭代即可保证收敛到Karcher均值，但需要主要的是在曲面上给定点集的Karcher均值可能不止一个。

在离散曲面上计算指数映射是比较容易的，利用测地线没有弯曲的定义我们只需要从\(p\)点出发沿给定方向\(u\)移动所需的距离即可。但需要注意的是在离散情况下当顶点的angle defect小于0时可能会出现指数映射不存在的问题，此时从\(p\)点出发沿任意方向移动都无法到达\(q\)点。

同样在离散曲面上计算对数映射也是相对容易的，稍后在Heat Method中会进行更加详细的介绍。

Dynamic Perspective

从动力学的角度上看，测地线可以理解为在曲面上没有发生改变(切向加速度为0)的曲线。但需要注意的是在曲面上的不同位置定义了不同的切平面，因此不能直接对不同位置的切向量进行比较，我们需要想办法去度量不同位置上的切向量。

类似于方向导数，我们定义曲面上的沿自身方向的变化率为协变导数(covariant derivative) \(\nabla_{\dot{\gamma}} \dot{\gamma}\)。根据测地线的定义可以得到测地线上沿自身方向的协变导数为0，\(\nabla_{\dot{\gamma}} \dot{\gamma} = 0\)。

协变导数实际上度量了曲面上向量场的变化率：假设曲面上存在两个向量场\(X\)和\(Y\)，在\(p\)点我们构造一条曲线\(\gamma(t)\)使得\(p\)点的切向量为\(X(p)\)，这样我们可以将曲线上来自\(Y\)的向量表示为\(Y'(t) = Y(\gamma(t))\)。最后我们对\(Y'(t)\)求导并减去曲面的法向量就得到了\(Y\)相对于\(X\)的协变导数\(\nabla_X Y\)。

此外我们也可以只利用曲面自身的局部性质来定义协变导数。对任意函数\(\phi\)，向量场\(X\),\(Y\),\(Z\)，协变导数可由以下性质唯一确定：

\[\nabla_Z (X+Y) = \nabla_Z X + \nabla_Z Y\]
\[\nabla_{X+Y} Z = \nabla_X Z + \nabla_Y Z\]
\[\nabla_{\phi X} Y = \phi \nabla_X Y\]
\[\nabla_X(\phi Y) = (D_X \phi) Y + \phi \nabla_X Y\]
\[D_Z g(X, Y) = g(\nabla_Z X, Y) + g(X, \nabla_Z Y)\]
\[\nabla_X Y - \nabla_Y X = [X, Y]\]

其中\([ \ ]\)为李括号(Lie bracket)，\([X, Y] = \sum_{i, j} (X^j \frac{\partial}{\partial u_j} Y^i -Y^j \frac{\partial}{\partial u_j} X^i) \frac{\partial}{\partial u^i}\)。同时可以证明协变导数由曲面的黎曼度量\(g\)所唯一确定，利用Christoffel记号\(\Gamma\)可以得到：

\[\Gamma_{ij}^p = \frac{1}{2} g^{pj} (g_{ij, k} + g_{jk, i} - g_{ki, j})\]

最后利用测地线协变导数为0的性质得到测地线方程(geodesic equation)：

\[\ddot{\gamma}^k + \Gamma_{ij}^k \dot{\gamma}^i \dot{\gamma}^j = 0\]

Heat Method

综合前面的章节不难发现想要精确地计算曲面上的测地线是比较困难，不过在曲面上计算测地距离实际上是相对容易的，这里介绍一种基于热传导方程的测地距离计算方法。

热传导问题与测地距离间可以通过Varadhan公式(Varadhan’s formula)联系起来：

\[\phi(x, y) = \lim_{t \to 0} \sqrt{-4t \log k_t(x, y)}\]

其中\(x\)表示热源位置，\(y\)为曲面上任意点，\(\phi(x, y)\)为曲面上\(x\)和\(y\)之间的测地距离，\(t\)为传热时间，\(k_t(x, y)\)为热核函数。Varadhan公式说明当传热时间趋于0时\(x\)和\(y\)之间的测地距离由热核函数\(k_t(x, y)\)控制，其物理意义是在极短时间内分子沿测地线进行扩散。

有了Varadhan公式后一个直接的想法是将热核函数带入公式中并进行积分从而求解出测地距离，不过遗憾的是直接使用Varadhan公式在大部分情况下是行不通的。当热核函数存在误差时，直接使用Varadhan公式计算出的测地距离会有非常显著的误差：

尽管Varadhan公式计算出的测地距离存在比较大的误差，我们可以认为测地距离变化的仍是比较准确的。实际上热传导方法计算测地距离的核心就是通过Varadhan公式来估计测地距离的梯度，然后带入Eikonal方程(eikonal equation)来求解出实际的距离场：

\[\vert \nabla \phi \vert = 1\]

Eikonal方程的几何意义为(测地)距离场在曲面上任意位置扩散的速度都是1。综合上面的内容可以得到计算测地距离的算法框架如下：

在曲面上给定点\(x\)上施加一个热源并进行扩散得到温度场；
计算温度场的梯度并进行归一化得到距离场扩散的方向\(X\)；
求解Eikonal方程得到测地距离场\(\phi\)。

Time Discretization

下面介绍如何在离散曲面上使用热传导方法。我们可以在\(x\)上设置温度为1，然后通过隐式迭代求解热传导方程：

\[(I - t \Delta) u^t = u^0\]

对于时间步长\(t\)可以取为网格边长度平均值的平方\(t = h^2\)。

Spatial Discretization

接下来计算温度场\(u^t\)的梯度。在每个三角形上温度场的梯度为：

\[\nabla u = \frac{1}{A_f} \sum_i u_i (N \times e_i)\]

对梯度进行归一化后得到向量场\(X\)，其散度为：

\[\nabla \cdot X = \frac{1}{2} \sum_j \cot \theta_1 (e_1 \cdot X_j) + \cot \theta_2 (e_2 \cdot X_j)\]

最后求解泊松方程即可得到所需距离场：

\[\Delta \phi = \nabla \cdot X\]