Discrete Differential Geometry 离散微分几何笔记

前言

最近在上Keenan Crane的DDG课程，感觉过去在学CG时涉及到mesh处理很多不明就里的操作都有了理论上的支持，并深深地感受到了微分几何在CG领域的强大作用。离散微分几何（DDG）的这一研究主题的目的就是如何将微分几何中已经解决的问题转变为计算机可以处理的离散问题，并关心这样的转化会带来什么样的差异。

本系列笔记参考了Keenan Crane DDG 2019的课程内容及他的一本小书DDG: An Applied Introduction。感兴趣的朋友可自行前往查阅英文原版教材及习题。

在学习离散微分几何的过程中，最大的感触就是对连续版本的微分几何不够熟悉，基础不够扎实，以至于在推导离散版本的公式的时候遇到了很多困难。但换一个方向说，这也恰恰是离散微分几何的魅力。在微分几何的理论世界中，很多概念无疑是很难理解的，或者说至少不是那么直观的。但是当这些问题被成功转化为离散的版本之后，我们就可以直观地体会到微分几何中所学习到的性质是多么的优美和简洁。

概览

在此先列出本系列笔记打算覆盖的内容：

微分几何基础介绍。
Exterior Calculus基础介绍。
离散几何的基本形式，如单纯形、复合形、网格等。
离散几何的基本操作：基本运算、求导、积分、拉普拉斯算子（即求解Laplacian/Poisson Equation的过程）等。
离散几何的操作在CG中的应用，如平滑，参数化，优化网格等。

和Machine learning中一个核心的概念类似的是，DDG中也存在所谓的No free lunch原则，即将某一个微分几何概念离散化的过程必然导致其某些性质的缺失。换言之，对于任何一个微分几何中的概念，如斜率、曲率、微分算子、梯度算子等等，我们都可以找出许多种方法来对其离散化，而这些方法都只保留了原概念中一部分的性质。以下是Keenan总结的解决一个离散微分几何问题的通用思路：

对于要离散化的微分几何概念，先找出多个在Smooth setting中等价的定义。
找到对应于每一种定义的离散化表示。
研究这些离散化表示都保留了哪些性质，并针对不同的应用选择合适的离散化方式。

在原书中针对这个问题给出了一个很经典的关于curvature离散化的讨论，此处不做赘述。

工具

对离散微分几何的研究涉及到以下的知识点:

Exterior Calculus：这是研究离散微分几何最有用的工具之一，必须灵活掌握。
传统微分几何。
CG中基于Halfedge的网格表示。
线性方程组的求解及常见的优化问题的求解。

微分几何基础

在第一篇笔记中我们首先复习一下基础的微分几何知识，这里只讨论进行DDG学习所必要的内容，细致的微分几何请参照数学系的教材进行学习。

微分几何关心什么？

微分几何主要关心流形（Manifold）的性质。流形可以理解为欧氏空间(Euclidean Space)中曲线和曲面等概念的延伸，是一个局部近似欧氏空间的可弯曲空间。我们对欧氏空间已经非常熟悉，在欧氏空间中有曲线，有曲面，而微分几何则关心这些曲线和曲面在更广义的空间中有怎么样的表现，其中包括：

在广义空间中如何定义函数？
如何对函数求导和积分？
如何求极值或求解一个优化问题？
流形的长度？面积？体积？
如何求解微分方程？
……

我们将这样一个广义的空间称作 $M$ ，即Manifold。微分几何的核心假定就是 $M$ 在局部上应该与 $\mathbb{R}^n$ 近似，即存在一个双射（bijection）使得在 $M$ 中的任意一个点 $\alpha$ 的邻域 $U_\alpha$ ，都有 $\phi_\alpha: U_\alpha \rightarrow \mathbb{R}^n$ 。我们将这样的一个双射 $\phi_\alpha$ 称作charts，表示其局部近似性就像棋盘一样，每一格都近似于一个欧氏空间。

自然而然地，我们需要思考如何用数学来表示这样一个广义的空间，这可以通过一个将空间 $M$ 中的点映射到某个实数的函数来表示，即 $f: M \rightarrow \mathbb{R}^n$ 。同样我们可以定义从一个流形到另一个流形的函数 $f: M \rightarrow N$ 。前者我们马上会在下文中展开讨论，而后者则可以理解为两个流形之间的线性近似，这会为我们后面引入切空间（Tangent space） $T_pM$ 埋下伏笔。

曲线

定义与度量

最简单的流形就是一维的流形。我们可以将常见的三维曲线定义为从一个一维流形到一个三维实数空间的变换 $\gamma: I \rightarrow \mathbb{R}^3$ 。这里的 $I$ 表示interval，与曲面的 $M$ 表示广义的Manifold对应。类似的，我们关心在原空间中的切向量 $X$ 在变换后的空间会是什么样。首先我们知道其距离可以表示为：

$\lvert d\gamma(X) \lvert = \sqrt{g(X,X)} = \sqrt{d\gamma(X) \cdot d\gamma(X)}$ 其中 $g(\cdot)$ 就是metric。这里 $(\cdot)$ 是广义上的内积，与欧氏空间中的向量内积不同（一个简单的例子就是函数的内积即 $f(X) \cdot f(Y) = \int_{M}f(X)f(Y)dM$ ）。为了表示区分，我们将会使用 $\langle X, Y \rangle$ 表示传统的向量内积。

在三维曲线上最常见的变换是等距变换，它对于任意 $X$ 都满足 $\lvert d\gamma(X)\lvert = \lvert X\lvert$ 。

Frenet Frame

接下来我们来系统地看一下三维曲线上的切线、法线、和副法线构成的Frenet Frame。首先我们进行如下的定义：

是中曲线的单位切向量。
- 首先注意 $T$ 是在 $\mathbb{R}^3$ 上的，与 $I$ 无关，也就是这个定义是落在欧氏空间内的。
- 其次要注意 $T$ 是单位向量，即 $\lvert T\lvert = 1$ 。
接下来考虑单位切向量的变化率，这可以表示为。
- $\kappa$ 被称作曲率（Curvature），可以理解为曲线上密切圆的半径的倒数，即密切圆半径越大，曲线越近似于直线，密切圆半径越小，曲线越近似于点。
- $N$ 被称作法向量（Principal normal），是垂直于 $T$ 的一个单位向量，可以视作quarter turn。
- $T$ 之所以必定垂直与 $N$ ，是因为二者都是单位向量，任何有平行于 $T$ 的变化量都会使得 $T$ 的长度变化，因此变化的方向只可能发生在法向上。
我们将称作副法向量（Binormal），这样就构成了一个以点为原点的正交系，即Frenet frame。
- 此时可以推导出 $T,N,B$ 的变化率满足如下关系： $\begin{bmatrix} dT\\ dN\\ dB \end{bmatrix} = \begin{bmatrix} 0 & -\kappa & 0\\ \kappa & 0 & -\tau\\ 0 & \tau & 0 \end{bmatrix}\begin{bmatrix} T\\ N\\ B \end{bmatrix}$
- 上式中 $\tau$ 为挠率，即 $B$ 沿 $N$ 方向变化的速度，可以理解为扭转的速率。

曲面

定义

从曲线自然而然地可以延伸到曲面。平滑的曲面是微分几何中研究的一个重点，也是DDG最为关心的部分。从工业界的应用来说，很多时候我们并不关心实体，而更关心表面的变化，尤其在可视化中，我们并不关心一个三维物体的内部，而只关心从外面看到这个三维物体应该是什么样的。在这样的情况下，对曲面的性质及变换的研究就显得格外的重要。

最简单的定义曲面的方式就是一个将某个二维流形变换到实数空间的映射，简单说来，即 $f: M \rightarrow \mathbb{R}^3$ 。当然了，根据上文的介绍，我们也可以将曲面看做是一个流形，而 $f$ 是这个流形的数学表示。

曲面的微分（切向量）

我们关心曲面映射的微分 $df$ 。 $df$ 将 $M$ 上的某一点 $p$ 处的切向量 $X$ 变换到实数空间中的 $df(X)$ 。我们之所以关心微分（differential）而不是这个映射本身，是因为直接刻画 $f$ 可能是一件很复杂的事情，因为 $M$ 是可以弯曲的。换言之，我们更关心 $M$ 上面的某一个切向量（相当于先将M在某个点附近用欧氏空间近似，再研究该点附近的一个向量）会如何随着 $f$ 的改变而改变——原空间 $M$ 中的一个向量 $X$ ，在 $f$ 使得 $M$ 变为 $f(M)$ 的同时，随之变为 $df(X)$ 。换言之，微分 $df$ 刻画了如何将一个向量从某个空间变换到另一个空间。至于为什么是切向量，可以参照An Introduction to Manifold的第二章，该章节系统地证明了欧氏空间中的切向量与导数的等价性。换言之，因为流形在局部的性质很好，所以我们可以用线性空间的方法去研究它。表达这样一个关系最重要的式子，就是 $T_p(\mathbb{R}^n) \cong \mathcal{D}_p(\mathbb{R}^n)$ ，它表示点 $p$ 的切空间与该点处的导数空间是等价的，从而有：

$v = \sum v^i e_i = \left. \sum v^i \frac{\partial}{\partial x^i} \right\rvert_p$

其中 $e_i$ 指的是切空间中的基，而 $\frac{\partial}{\partial x^i} \rvert_p$ 指的是欧氏空间中 $p$ 点的偏微分，二者是等价的。

值得注意的是，如果我们考虑的不是单个向量，而是一个向量场，那么 $df$ 同样可以将整个向量场从某个空间变换到另一个空间。我们定义一个 $\mathbb{R}^n$ 中的一个开子集 $U$ 上的向量场为：

$X_p = \left. \sum a^i(p) \frac{\partial}{\partial x^i} \right\rvert_p, p \in U, a^i(p) \in \mathbb{R}$

注意我们重复使用了 $X$ 这个符号，在前面它表示的是一个向量，这里表示的是一个向量场函数。我们将 $X$ 所处的空间 $T_pM$ 定义为切空间，即 $M$ 的所有切向量的表示。我们同样可以用导数的方法来理解向量场，假定对于开集 $U$ 上有连续函数 $g$ ，那么我们可以定义一个新的函数 $Xg(p)$ 如下：

$(Xg)(p) = X_pg = \sum a^i(p)\frac{\partial g}{\partial x^i}(p)$ 这就是向量场的导数表示。

到这里可能都看懵了，讲的啥跟啥啊。其实讲了这么多，我们就是希望把欧氏空间中的导数和流形中的切向量（切向量场、微分）联系起来。下面这个例子可以很好地展示如何将传统意义上的导数和几何上的切向量联系起来。

Example 1: Devivative of a real-value function

对于一个实值函数而言，我们经常将它画成下图左边的样子，即 $\phi(x)$ 随着 $x$ 的变化而变化的一张函数图像。事实上，从微分几何的视角来看，我们可以认为 $x$ 所在的空间为 $M$ （这里 $M$ 也是 $\mathbb{R}^1$ ），而 $\phi(x)$ 所在的空间为 $\mathbb{R}^1$ 。从这个角度看，在M空间中的切向量 $X$ 经过 $\phi$ 变换后所对应的变化 $d\phi(X)$ 恰好就是 $\phi'$ ，即 $\phi$ 的导数。我们将这个关系写作：

$\phi' = d\phi(X)$

这就将几何上的微分算子和欧氏空间中的导数联系了起来，当然了，这里的流形 $M$ 本身也是一个欧氏空间，因此这个对应是很简单也很直接的。我们用这个例子只是为了直观地展示一下微分究竟做了一件什么事情。

derivative of a real-value function

接下来让我们再看另一个关于方向导数的例子。

Example 2: Devivative of a real-value function

假定现在我们的函数变成了二维上的实值函数，即 $\phi: \mathbb{R}^2 \rightarrow \mathbb{R}$ ，那么对应的微分 $d\phi(X)$ 是什么呢？答案是方向导数（directional derivative），即在某一点 $p$ 给定一个方向 $X$ ，函数 $\phi$ 沿该方向的斜率。要理解这一点也很简单，我们的流形空间是 $M = \mathbb{R}^2$ ，所以 $X$ 应该是 $\mathbb{R}^2$ 上的一个切向量，给定这个向量我们想了解 $\phi(x)$ 会如何变化，自然而然就会联想到方向导数。当然了，这里的 $M$ 也仍然是欧氏空间。但我们引入这一系列微分几何的讨论，都是为了使得我们讨论的性质在 $M$ 未必是欧氏空间时依然成立。

度量（Metric）

有了曲面微分的定义，我们的第一个问题就是切向量 $dX$ 变到 $df(X)$ 之后其长度发生了什么变化？通常我们采用如下的方式来表示经过 $f$ 变换后的 $X$ 的长度：

$\lvert df(X) \rvert = \sqrt{g(X,X)} = \sqrt{df(X) \cdot df(X)}$

这与前面曲线的长度的定义是一致的。有了度量，我们可以刻画一些特殊的变换，如等距变换(isometric)和保角变换(Conformal)：

等距变换：要求 $f$ 不发生拉伸，即 $\lvert df(X)\lvert = \lvert X\lvert$ 。
保角变换：要求不改变两向量之间的夹角，即，其中是一个取值恒为正数的函数。
- 保角变换允许曲面发生拉伸，但不允许发生剪切。
- 保角变换总是存在，这一点很重要。

法向量（Normal）

有了切向量，我们可以对应地定义法向量，即，如果向量 $u\in \mathbb{R}^3$ 满足 $df(X)\cdot u = 0$ ，那么 $u$ 就是 $p$ 处的法向量。特别的，我们将单位法向量记作 $N$ 。在有向曲面上， $N$ 的方向可以唯一选定，而在无向曲面上（如莫比乌斯环）， $N$ 的方向不定。由于 $N$ 的长度为1，我们可以同样可以将法向量表示为一个从流形到单位球的映射，即 $N: M \rightarrow S^2$ 。这里可以观察到一些有趣的相似点：

$f$ 将 $M$ 映射到 $\mathbb{R}^3$ ，我们可以用 $df(X)$ 来刻画切向量 $X$ 在欧氏空间中本身的变化。
$N$ 将 $M$ 映射到 $S^2$ ，我们同样可以用 $dN(X)$ 来刻画切向量 $X$ 在欧氏空间中法向上的变化。

注意，二者的输入都是切向量 $X$ ，这是因为 $M$ 的切空间与欧氏空间近似。这里我们将 $N$ 称作Gauss map，将 $dN$ 称作Weingarten map。

曲率（Curvature）

法向曲率（Normal curvature）

前面我们已经利用Frenet Frame研究了曲线上某一点处的切向量，法向量和副法向量，并由此引出了曲率和挠率的定义。这个定义如何引申到曲面上呢？答案是我们可以继续研究曲面上的曲线的曲率并将之作为曲面上某一点的曲率。

这引出了法向曲率（Normal curvature）的定义，即假定在 $p$ 点处存在一个由 $df(X)$ 和 $N$ 张成的平面 $\pi$ 与曲面相交得到一条曲线，那么该曲线的曲率就是这个曲面在该点处的法向曲率，记作：

$\kappa_n(X) = \frac{df(X) \cdot dN(X)}{\lvert df(X)\lvert ^2}$

这个式子来自于前面提到的Frenet frame性质，即 $dN = \kappa T - \tau B$ ，我们可以利用 $T$ 去点积 $dN$ 来提取其中的曲率部分，即 $T \cdot dN = T \cdot(\kappa T - \tau B) = \kappa \lvert T\lvert ^2$ ，从而得到上式。

Shape operator

注意在法向曲率的定义中我们选定了切向量 $X$ ，一个随之而来的问题是：选择什么方向 $X_1$ 和 $X_2$ 会得到最大和最小的法向曲率 $\kappa_1$ 和 $\kappa_2$ 。为了解决这个问题，我们需要引入一个新的算子——Shape operator，它是一个从切空间到切空间的映射，可表示为：

$df(SX) = dN(X), S: TM \rightarrow TM$

上式表达的含义是， $S$ 作用于切空间中的 $X$ 后再经过 $f$ 变换，等价于 $dN$ 变换，即切向量 $X$ 在欧氏空间中法向的变化。它与 $dN$ 的区别就在于 $S$ 变换后仍然在 $TM$ 空间，而 $dN$ 变换后在 $\mathbb{R}^n$ 空间。当然了，这里有一个隐藏的性质，那就是为什么 $N$ 的变化总可以用 $TM$ 空间中的一个向量来表示？这是因为单位法向量 $N$ 的长度确定，因此它的变化必然发生在切平面上，也就必然落在 $TM$ 空间中。

总的来说， $S$ 算子接受一个切向量 $X$ 作为输入，然后告诉你在变换后的空间中沿 $X$ 方向法向量的变化率，即曲率。对于曲面而言 $S$ 是一个矩阵，这个矩阵的eigenvalue和eigenvector就反映了最大曲率和最小曲率以及这些曲率对应的方向，即：

$SX_i = \kappa_i X_i$

Mean and Gaussian curvarture

由于曲面所在的切空间是二维的，知道了principal curvature，我们可以很轻松地计算出任何一个方向的曲率。当然了，principal curvature并不容易计算，因此我们往往会关心两个替代品：Mean curvature $H = \frac{\kappa_1 + \kappa_2}{2}$ 和Gaussian curvature $K = \kappa_1\kappa_2$ 。以下对二者的性质作简要介绍：

mean可以认为是logical or，即关心比较大的曲率；而gaussian可以认为是logical and，即关心比较小的曲率。
gaussian为0时，即至少有一个主曲率方向为0，换言之，该方向是平坦的（想象一个柱面）。此时该曲面称作developable，即可以不经伸缩就变换为一个平面。
mean为0时，表示最大曲率和最小曲率互为相反数，此时曲面为马鞍面，也称作minimal surface，即给定曲面边界时面积最小的曲面形式。
$\kappa_1 = \kappa_2$ 时，任意方向的曲率相等，此时该点称作umbilic point。

Fundamental forms

最后简要提一句Fundamental form，即基本形（然而并没有什么用）。第一基本形即induced metric，即：

$I(X,Y) = g(X,Y)$

第二基本形是Shape operator变换后的切向量与另一个切向量的负内积，即：

$II(X,Y) = -g(SX, Y) = -dN(X)\cdot d(Y)$

引入坐标系

到这里为止基本就是离散微分几何会用到的微分几何概念了，毕竟我们凡人只和三维欧氏空间打交道，掌握好曲线和曲面已经可以开始进一步的学习了。然而前面的讨论都是抽象的，即我们并不关心原空间 $M$ 是什么样子，也不关心变换后的欧氏空间有几维，这有利于我们学习抽象的概念而不是拘泥于具体的运算。不过呢，最终DDG还是要落到实际的运算中的，而运算的第一步就是引入坐标系。

我们可以将讨论的范围限制在 $f: \mathbb{R}^2 \rightarrow \mathbb{R}^3$ 中，即认为 $M$ 是 $\mathbb{R}^2$ 而 $f(M)$ 是 $\mathbb{R}^3$ 。这样的话，由于原空间和变换后的空间都是欧氏空间，我们就可以利用正交坐标系来表示前面学习到的算子了。

首先我们来看一下 $f$ 是什么，由于在 $\mathbb{R}^2$ 中切向量就是一个该空间中的向量本身， $f$ 无非就是一个将二维向量map到三维向量的变换，即：

$f(x_1, x_2) = (f_1(x_1, x_2). f_2(x_1, x_2), f_3(x_1, x_2))$

那 $df$ 无非就是一个雅克比矩阵，即:

$df = J = \begin{bmatrix} \frac{\partial f_1}{ \partial x_1} & \frac{\partial f_1}{ \partial x_2}\\ \frac{\partial f_2}{ \partial x_1} & \frac{\partial f_2}{ \partial x_2}\\ \frac{\partial f_3}{ \partial x_1} & \frac{\partial f_3}{ \partial x_2} \end{bmatrix}$

那么将微分作用于切向量也很简单，即 $df(\mathbf{x}) = \mathbf{Jx}$ 。

那metric作为两个向量的内积，其表示也很简单：

$g(\mathbf{x},\mathbf{y}) = (\mathbf{Jx})^T(\mathbf{Jy}) = \mathbf{x}^T\mathbf{I}\mathbf{y}$

其中 $\mathbf{I} = \mathbf{J}^T\mathbf{J}$ 就是第一基本形。类似地我们可以导出其他算子在正交坐标系中的表示，此处不作赘述。

总结

本文介绍了DDG所使用到的工具并初步介绍了相关的微分几何知识。尽管微分几何的证明细节很多，但重要的是把握以下的核心思路：

流形是曲线和曲面的高维抽象表示。
流形通常是非线性的，但是局部可以用一个线性空间去近似。
有了线性空间我们就可以使用线性代数来研究流形在局部上的性质。
流形上某一点的切向量与对应的欧氏空间中的导数存在duality。
我们当然可以通过分析切向量来研究流形的性质，但是，利用对偶性，分析函数显然是一个更好的选择。

之后我们还会补上外微积分的内容。有了这些基础知识，我们就可以开始推导离散版本的微分几何概念了，并开始研究如何将这些离散的概念应用到实际的网格处理中。在文末我们简要提几个以后会重点解决的问题：

如何对一个离散的网格求导？
如何对一个离散的网格积分？
给定网格的初始状态和一个扰动或一个动力源（如重力场、碰撞等等），如何求解网格的变化？
如果网格表示的是具体的实物，如何表示实物的物理性质，如弹性、刚度、挠度等等？

这些问题都涉及到将连续的微分几何概念化为离散的网格中的概念来进行求解，其应用也不限于计算机图形学，有限元（FEM）的计算与DDG密切相关。每当一想到DDG如此有趣，就希望我自己不要偷懒摸鱼，可以早日写完这个系列（误）。