这个系列是MIT 6.838: Shape Analysis的同步课程笔记。本课程会介绍几何方法在图形学、机器学习、计算机视觉、医疗图像以及建筑设计等相关领域的原理和应用。本节主要介绍离散曲面上的Laplace算子。
Discretizing the Laplacian
在上一节课中我们介绍了Laplace算子在不同作用域下的定义和性质。而对于几何处理问题，我们更关心如何在离散曲面上来定义Laplace算子。我们希望离散情况下的Laplace算子能够近似它在连续情况下的各种性质。
回忆\(\mathbb{R}^n\)上Laplace算子可以定义为函数对于坐标轴二阶导数的和，对于曲面的情况可以使用类似的方法进行定义。
实际上数学家已经推导过曲面上的Laplace算子。不过它的表达式非常复杂，在绝大多数情况下我们都会尽量避免直接使用它。
在离散曲面上定义Laplace算子的难点在于它是一个微分算子，对于很多离散几何对象有时我们往往很难去定义各种微分运算。
本节课中我们会利用有限单元法(finite element method, FEM)来推导Laplace算子的近似形式。
Laplacian from FEM
Integration by Parts
在推导离散Laplace算子之前首先简要介绍一下FEM的基本思想。回忆到目前为止我们大量使用了分部积分的技术来处理各种积分问题。
对于Laplace算子而言，通过分部积分可以把二阶微分算子转换为一阶微分(梯度)。在大多数情况下梯度的计算要远比二阶微分要容易得多，有时我们甚至可以利用函数的各种性质来避免一些计算。
Poisson Equation
FEM的本质是把微分方程转换为代数方程进行求解。以泊松方程(Poisson equation)为例，一维泊松方程可以表示为：
\[\Delta f = g\]显然在无穷的函数空间中求解泊松方程是非常困难的。这里我们引入一系列基函数：
\[\phi_1(x), \phi_2(x), \dots \phi_k(x) : \mathcal{M} \rightarrow \mathbb{R}\]利用基函数可以对\(f\)和\(g\)进行分解：
\[f(x) = \sum_k v^k \phi_k(x)\] \[g(x) = \sum_k w^k \phi_k(x)\]这样我们只需要求解系数\(v^k\)和\(w^k\)即可：
\[\sum_k v^k \Delta \phi_k(x) = \sum_k w^k \phi_k(x)\]直接求解上式仍然是非常困难的。这里我们引入测试函数(test function)的概念，放松条件仅要求等式两端与测试函数的内积必须相等。测试函数的一种常见选取方法是选择一个基函数\(\phi_l\)，此时有：
\[\sum_k v^k \int \phi_l(x) \ \Delta \phi_k(x) \ dx = \sum_k w^k \int \phi_l(x) \phi_k(x) \ dx\]通过轮换测试函数\(\phi_l\)可以得到关于系数\(v^k\)和\(w^k\)的线性方程组，其矩阵形式为：
\[L \mathbf{v} = M \mathbf{w}\] \[L_{kl} = \int \phi_l(x) \ \Delta \phi_k(x) \ dx\] \[M_{kl} = \int \phi_l(x) \phi_k(x) \ dx\]这样原始的微分方程就转换为代数方程，我们也就可以使用求解线性系统的相关技术来进行处理。
Galerkin FEM Approach
实际上上面将微分方程转换为代数方程的方法称为Galerkin FEM approach，它的核心步骤包括两步：使用基函数对待求函数进行分解，然后通过函数内积来构造线性系统。
除此之外Galerkin方法也可以从variational calculus的角度来进行理解，泊松方程等价于求解如下优化问题：
\[\min_{f: \mathcal{M} \rightarrow \mathbb{R}} \int \bigg[ \frac{1}{2} \Vert \nabla f \Vert_2^2 - f \cdot g \bigg] \ d x\]通过基函数分解，优化问题可以整理为：
\[\begin{aligned} &\ \min_{v^1, \dots, v^k} \ \int \bigg[ \frac{1}{2} \big\Vert \sum_k v^k \nabla \phi_k \big\Vert_2^2 - \sum_{k, l} v^k w^l \phi_k \phi_l \bigg] \ d x \\ =&\ \min_{v^1, \dots, v^k} \ \int \frac{1}{2} \mathbf{v}^T L \mathbf{v} - \mathbf{v}^T M \mathbf{w} \ d x \\ =&\ \min_{\mathbf{v} \in \mathbb{R}^k} \ \frac{1}{2} \mathbf{v}^T L \mathbf{v} - \mathbf{v}^T M \mathbf{w} \end{aligned}\]而它的解恰为线性系统的解：
\[L \mathbf{v} = M \mathbf{w}\]
需要说明的是Galerkin方法并不是推导FEM唯一的方法，在不同的工程问题中还有其它将微分方程转换为代数方程的方法。
\(L^2\) Dual of a Function
从函数的角度来进行考虑，Galerkin方法实际上是把函数\(f\)转换为它的对偶形式(算子)\(\mathcal{L}_f\)，然后利用测试函数\(g\)来进行求解原始的微分方程。
通过选择合适的测试函数\(g\)我们还可以从对偶形式来恢复函数\(f\)。
当我们使用Laplace算子来构造对偶形式时甚至可以避免直接计算函数的二阶微分。
Laplacian on Triangle Meshes
Hat Function
总结一下，使用Galerkin方法构造和求解FEM问题的核心在于选择合适的函数空间以及测试函数，而且在大多数情况下这两个函数空间是相同的。
在三角网格上最常用的基函数是在给定顶点\(v\)上取1其它顶点取0，然后在相邻三角形上进行线性插值的函数。这样的基函数也称为hat function。
这样三角网格上的任意函数就可以表示为顶点函数值的线性组合：
\[f(\mathbf{x}) = \sum_{v} f^v h_v(\mathbf{x})\]
Piecewise-Linear Elements
接下来考虑hat function的微分。具体来说，对于Laplace算子我们需要计算hat function的梯度以及内积，然后在每个三角形上再进行积分。
Gradient of a Hat Function
我们首先来推导hat function的梯度。假设\(f\)在三角形顶点\(v_1\)上取1并且在另外两个顶点\(v_2\)和\(v_3\)上取0，同时函数值在三角形内部采用线性插值。此时梯度\(\nabla f\)在三角形内部是一个常数，而且函数可以表示为梯度与坐标的线性组合：
\[f(\mathbf{x}) = (\nabla f)^T \mathbf{x} + c\]带入三个顶点可以得到：
\[\nabla f \cdot (\mathbf{v}_1 - \mathbf{v}_2) = 1\] \[\nabla f \cdot (\mathbf{v}_1 - \mathbf{v}_3) = 1\] \[\nabla f \cdot (\mathbf{v}_2 - \mathbf{v}_3) = 0\]这表示梯度\(\nabla f\)的方向垂直于底边\(e_{23}\)。记边\(e_{12}\)和\(e_{13}\)的边长分别为\(l_2\)和\(l_3\)，整理后得到：
\[\Vert \nabla f \Vert_2 \cdot l_3 \cdot \cos (\frac{\pi}{2} - \theta_3) = \Vert \nabla f \Vert_2 \cdot l_3 \cdot \sin (\theta_3) = 1\] \[\Vert \nabla f \Vert_2 = \frac{1}{l_3 \cdot \sin (\theta_3)} = \frac{1}{h}\]即梯度\(\nabla f\)的大小为底边高\(h\)的倒数。记底边向量为\(\mathbf{e}_{23} = \mathbf{v}_3 - \mathbf{v}_2\)，结合\(\nabla f\)的方向可以得到梯度的表达式：
\[\nabla f = \frac{\mathbf{e}_{23}^{\perp}}{\Vert \mathbf{e}_{23}^{\perp} \Vert_2} \cdot \frac{1}{h} = \frac{\mathbf{e}_{23}^{\perp}}{2 A}\]其中\(\mathbf{e}_{23}^{\perp}\)表示对底边向量逆时针转90度，而\(A\)则是三角形的面积。上式说明hat function的梯度可以使用顶点的坐标来进行表示。
Inner Product of Gradients
对于三角形上的内积则需要考虑同一个顶点上的hat function以及不同顶点上的hat function两种情况。对于同一顶点上两个hat function梯度的内积有：
\[\int \nabla f \cdot \nabla f \ d A = A \cdot \Vert \nabla f \Vert_2^2 = \frac{A}{h^2} = \frac{1}{2} (\cot \alpha + \cot \beta)\]
而对于不同顶点上的hat function的情况则有：
\[\begin{aligned} \int \nabla f_\alpha \cdot \nabla f_\beta \ d A &= A \cdot \langle \nabla f_\alpha , \nabla f_\beta \rangle \\ &= A \cdot \langle -\frac{\mathbf{e}_1^{\perp}}{2A} , \frac{\mathbf{e}_2^{\perp}}{2A} \rangle \\ &= -\frac{1}{4A} \langle \mathbf{e}_1 , \mathbf{e}_2 \rangle \\ &= -\frac{1}{4A} l_1 l_2 \cos \theta \\ &= -\frac{1}{2h} l_1 \cos \theta \\ &= -\frac{1}{2} \frac{\cos \theta}{\sin \theta} \\ &= -\frac{1}{2} \cot \theta \end{aligned}\]
把这两种情况整理一下不难发现三角形上梯度的内积只与三角形的内角有关。
对相邻的三角形进行遍历求和可以得到通过相邻顶点来定义的形式：
\[\langle \nabla h_{\mathbf{p}} , \nabla h_{\mathbf{q}} \rangle = -\frac{1}{2} (\cot \theta_1 + \cot \theta_2)\]
同时它和mean curvature normal有相同的形式：
Cotangent Laplacian
这样我们就推导出了三角网格上的Laplace算子。由于它仅与三角形内角的正切有关，因此也称为cotangent Laplacian。
利用cotangent Laplacian就可以在三角网格上求解泊松方程。
需要说明的是这种基于Galerkin方法求解出的结果是原始问题的一个弱解(weak solution)。
从形式来看，cotangent Laplacian实际上就是一个矩阵。
还需要说明的是在求解泊松方程时也要对等式右边进行一些处理。
最常用的方法是引入质量矩阵(mass matrix)。
构造质量矩阵时可以利用顶点的附近三角形面积的\(1/3\)来获得一个对角矩阵，这种方法在计算上有一些优越性。
Boundary Conditions
总而言之通过构造cotangent Laplacian和质量矩阵我们把离散曲面上的泊松方程转换成了一组代数方程，不过在实际求解时还需要考虑系统的边界条件。
Eigenhomers
类似于连续情况下特征函数的概念，在离散网格上Laplace算子同样有着满足eigenhomer来描述网格上函数变化的频率。
Higher-Order Elements
除了使用一阶网格外使用更高精度的单元同样可以推导出相应的FEM方法。
Point Clouds Laplacian
本节课的extra content介绍了点云上的Laplace算子。和网格相比点云缺少了点之间的连接关系，因此需要一些额外的处理。我们处理点云的主要方法是热传导方程(heat equation)，它的形式为：
\[\frac{\partial u(x, t)}{\partial t} = - \Delta u (x, t)\]利用Laplace算子的特征函数可以对\(u\)进行分解，带入热传导方程得到：
\[u(x, t) = \sum_k u^k(t) \ \phi_k(x)\] \[\sum_k (u^k)' (t) \ \phi_k(x) = - \sum_k u^k(t) \ \Delta \phi_k(x) = - \sum_k u^k(t) \ \lambda_k \phi_k(x)\]由于特征函数的正交性，等式两边的系数必须相等：
\[(u^k)' (t) = -\lambda_k u^k(t)\]这样我们就得到了\(k\)个ODE，而且每个方程的解为：
\[u^k (t) = u_0^k \cdot e^{-\lambda_k t}\]
在此基础上我们定义heat kernel为：
\[\mathcal{H}_t (x, y) = \sum_k e^{-\lambda_k t} \phi_k(x) \phi_k(y)\]利用heat kernel可以把\(u\)写成函数卷积的形式：
\[u(x, t) = \int_\mathcal{M} \mathcal{H}_t (x, y) \ u(y, 0) \ d y\]在\(\mathbb{R}^n\)上heat kernel有解析解：
\[\mathcal{H}_t (x, y) = \frac{1}{(4 \pi t)^{\frac{n}{2}}} e^{\frac{\Vert x-y \Vert_2^2}{4t}}\]不难发现它实际上就是\(n\)维正态分布的pdf。而在流形\(\mathcal{M}_m\)上，当时间\(t\)比较小时我们可以使用\(\mathbb{R}^n\)上的heat kernel进行近似：
\[\mathcal{H}_t (x, y) = \frac{1}{(4 \pi t)^{\frac{m}{2}}} e^{\frac{\Vert x-y \Vert_2^2}{4t}} \bigg( \phi(x, y) + O(t) \bigg)\]注意这里\(\phi(x, y)\)并不是特征函数，而是满足\(\phi(x, x)=1\)的光滑函数。上式说明当\(x\)和\(y\)比较接近而且时间\(t\)很小的时候，\(\mathbb{R}^n\)上的heat kernel与流形上的heat kernel是非常接近的。
回到heat equation，当时间\(t\)比较小时可以使用差商来近似导数：
\[\begin{aligned} \Delta u (x_i, t) &= -\frac{\partial u(x_i, t)}{\partial t} \approx \frac{u(x_i, 0) - u(x_i, t)}{t} \\ &= \frac{1}{t} \bigg[ u(x_i, t) - (4 \pi t)^{-\frac{m}{2}} \int_\mathcal{M} e^{\frac{\Vert x_i-y \Vert_2^2}{4t}} u(y, 0) \ d y \bigg] \end{aligned}\]对于流形上的积分则使用求和来进行近似：
\[\int_\mathcal{M} e^{\frac{\Vert x_i-y \Vert_2^2}{4t}} u(y, 0) \ d y \approx \frac{1}{K} \sum_j e^{\frac{\Vert x_i-x_j \Vert_2^2}{4t}} u(x_j, 0)\]引入系数\(W_{ij}\)就可以得到矩阵形式的Laplace算子：
\[L = D - W\] \[W_{ij} = \begin{cases} e^{\frac{\Vert x_i-x_j \Vert_2^2}{4t}} & \text{if $\Vert x_i-x_j \Vert_2 \leq \varepsilon$} \\ 0 & \text{otherwise} \\ \end{cases}\] \[D_{ii} = \sum_j W_{ji}\]
和网格上的Laplace算子相比上面推导出的算子无需样本之间的连接关系，因此特别适合点云这样的数据结构以及各种流形上的机器学习任务。
