DDG课程笔记2-曲面与曲率

这个系列是对CMU离散微分几何课程(CMU CS 15-458/858: Discrete Differential Geometry)的总结和回顾。有了外代数和外微分的基础以后，接下来我们开始讨论曲面和曲率的一些性质。

Surfaces

Parameterized Surface

对于三维空间中的曲面，我们可以把它看成是二维平面通过映射\(f: M \rightarrow \mathbb{R}^3\)嵌入到三维空间中，这样的曲面称为参数化曲面(Parameterized Surface)。映射\(f\)描述了二维平面如何在三维空间中进行「伸展」，对\(f\)进行微分就可以把平面向量\(X\)映射到空间中得到对应的切向量\(df(X)\)，换句话说\(df\)告诉了我们平面如何在局部进行了「伸展」。

将\(df\)写成矩阵的形式就得到了\(f\)对应的Jacobi矩阵(Jacobian matrix):

\[df(X) = J_f X\] \[J_f = \begin{bmatrix} \frac{\partial f^1}{\partial x^1} & \cdots & \frac{\partial f^1}{\partial x^n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f^n}{\partial x^1} & \cdots & \frac{\partial f^n}{\partial x^n} \\ \end{bmatrix}\]

当曲面通过\(f\)进行变换时可能会遇到一些退化的情况，比如说我们可以定义\(f\)把平面映射到球面：

\[f(u, v) = (\cos (u) \sin (v), \sin (u) \sin (v), \cos (v))\] \[\begin{aligned} df = &(-\sin (u) \sin (v), & & \cos (u) \sin (v), & &0 & &) & & du + \\ &(\cos (u) \cos (v), & &\cos (v) \sin (u), & &-\sin (v) & &) & &dv \end{aligned}\]

如果把\(f\)认为是平面地图到球面的映射，则\(u\)方向对应东西方向而\(v\)方向对应南北方向。当\(v=0\)我们会位于球面北极的位置，此时取\(X=(1, 0)\)得到\(df(X) = (0, 0, 0)\)。也就是说沿\(X\)方向前进我们始终会留在原地，换句话说我们没有办法在两极定义东西方向。

为了避免上述这种退化问题同时方便后面的分析和推导这里引入immersed surface的概念：定义\(f: U \rightarrow \mathbb{R}^n\)是一个immersion当且仅当\(df\)不会退化，即\(df(X) \vert_p = 0 \Leftrightarrow X \vert_p = 0\)。在本课程中除非特殊说明我们均假定\(f\)是一个immersion。

Induced Metric

对于参数化曲面我们总是希望可以在参数平面上研究曲面的在空间中的性质，比如说角度、长度、面积等等。以角度为例，显然参数平面上向量\(X\)和\(Y\)的夹角并不代表映射到空间后两个方向上的夹角。此外，对于同一个曲面可能存在多个不同的参数化方法，因此参数平面上向量的夹角并没有告诉我们任何关于曲面自身的信息。

尽管参数平面上向量的夹角没有什么特殊的意义，但通过微分映射\(df\)我们可以把参数平面上的向量转换到曲面的切平面上，从而得到曲面的信息。具体地我们定义切平面上\(df(X)\)和\(df(Y)\)的内积为：

\[g(X, Y) = \langle df(X), df(Y) \rangle\]

称\(g\)为\(f\)的诱导度量(induced metric)。将参数平面上的任意两个向量\(X\)和\(Y\)带入诱导度量得到：

\[g(X, Y) = \langle df(X), df(Y) \rangle = (J_f X)^T (J_f Y) = X^T J_f^T J_f Y = X^T \mathbf{I} Y\]

其中\(\mathbf{I} = J_f^T J_f\)称为曲面的第一基本形式(first fundamental form)。

第一基本形式完全描述了曲面的度量性质，我们可以以此计算曲面上的夹角、长度、面积等几何量。

Normals

除了切向量外我们往往还关心曲面的法向量。我们定义高斯映射(Gauss Map)为参数平面上的点到对应单位法向量的映射。显然单位法向量到在单位球面上，因此高斯映射也可以看做是参数平面到单位球上的映射。

我们还可以定义区域内的平均法向：

\[\frac{1}{A(\Omega)} \int_\Omega N dA\]

其中被积项\(N dA\)为大小等于面积微元的法向量，称为vector area。

我们将differential 2-form的运算法则推广到向量，用向量叉乘\(\times\)代替标量乘法可以得到：

\[\begin{aligned} df \wedge df (X, Y) &= df(X) \times df(Y) - df(Y) \times df(X) \\ &= 2 df(X) \times df(Y) \\ &= 2 N dA(X, Y) \end{aligned}\]

即\(N dA = \frac{1}{2} df \wedge df\)，再根据Stokes定理得到：

\[2 \int_\Omega N dA = \int_\Omega df \wedge df = \int_\Omega d (f df) = \int_{\partial \Omega} f df\]

上式表明对于给定的映射\(f\)，具有相同边界的的曲面会拥有相同的vector area；同时也表明对于封闭曲面其平均法向为0。

Discrete Surfaces

本节的最后来讨论一下曲面的离散表示。目前空间曲面最常用的表达形式仍然是三角网格，三角网格的连接关系表示曲面的拓扑信息，而网格顶点的坐标则表达了曲面的几何信息。

如果我们忽略网格的几何信息仅考虑拓扑关系，我们也可以把网格看做是一个流形(manifold)。那么接下来的问题是如何把流形嵌入到空间中？很显然我们可以构造顶点坐标的函数\(f\)将平面网格映射到空间中。也就是说我们只需要知道网格和映射\(f\)就足以表示空间曲面了。

离散曲面的优势在于曲面的每个面片都是平面，因此我们之前在\(\mathbb{R}^n\)推导的公式仍然是有效的。实际上，我们之前推导的所有离散微分算子都可以直接应用到离散曲面上而无需任何额外的处理；此外我们在推导时也没有考虑任何的几何关系，也就是说我们可以仅依赖于拓扑关系就可以对网格进行处理。

Smooth Curvature

Curvature of Curves

接下来介绍曲率相关的知识，首先来考虑光滑曲线的曲率。类似于参数平面，我们可以把曲线看成是直线\(L\)经过映射\(\gamma\)嵌入到平面（空间）中，此时称\(\gamma(t)\)为一条参数化曲线：

显然对于同一条给定的曲线存在不同的参数化方法，其中最重要的一种参数化是使用曲线的弧长作为参数，称为弧长参数化(arc-length parameterized)：

\[s = s(t) = \int_0^t \lVert \dot \gamma(t) \rVert dt\]

对于使用弧长参数化的空间曲线，曲线上的每一点都对应着一个局部坐标系称为Frenet标架(Frenet frame)，对应的三个坐标轴分别为切向\(T\)、法向\(N\)和从法向\(B\)：

\[T(s) = \frac{d}{ds} \gamma(s)\]
\[N(s) = \frac{d}{ds} T(s)\]
\[B(s) = T(s) \times N(s)\]

可以证明使用弧长参数化时曲线的切向量长度为1，因此三个坐标轴对应的向量均为单位向量。同时三个向量还满足Frenet-Serret公式(Frenet–Serret formulas)：

\[\frac{d}{ds} \begin{bmatrix} T \\ N \\ B \end{bmatrix} = \begin{bmatrix} 0 & \kappa & 0 \\ -\kappa & 0 & \tau \\ 0 & -\tau & 0 \end{bmatrix} \begin{bmatrix} T \\ N \\ B \end{bmatrix}\]

其中\(\kappa = \langle N, \frac{d}{ds}T \rangle\)描述了切向的变化率，称为曲率(curvature)；\(\tau = -\langle N, \frac{d}{ds}B \rangle\)描述了从法向的变化率，称为挠率(torsion)。

此外我们还可以考虑曲面上的曲线：设曲线上一点的曲面法向为\(N_M\)、从法向为\(B_M = T \times N_M\)，定义曲线的法曲率(normal curvature)\(\kappa_n\)和测地曲率(geodesic curvature)\(\kappa_g\)：

\[\kappa_n = \langle N_M, \frac{d}{ds}T \rangle\] \[\kappa_g = \langle B_M, \frac{d}{ds}T \rangle\]

法曲率代表曲线为了维持在曲面上所具有的曲率，而测地曲率则描述了曲线偏离测地线的程度。

Curvature of Surfaces

在曲线曲率的基础上开始介绍曲面的曲率。回忆Gauss map告诉了我们曲面上任意点的法向，对Gauss map进行微分得到Weingarten map \(dN\)，\(dN(X)\)表示法向如何沿方向\(X\)变化：

在此基础上我们定义曲面的法曲率(normal curvature)为法向沿切向\(df(X)\)的变化率：

\[\kappa_N(X) = \frac{\langle df(X), dN(X) \rangle}{\vert df(X) \vert^2}\]

举个简单的例子，考虑参数化柱面\(f(u, v) = (\cos(u), \sin(u), v)\)。分别计算Gauss map和Weingarten map得到:

\[df = (-\sin(u), \cos(u), 0) du + (0, 0, 1) dv\] \[N = (-\sin(u), \cos(u), 0) \times (0, 0, 1) = (\cos(u), \sin(u), 0)\] \[dN = (-\sin(u), \cos(u), 0) du\]

带入法曲率计算公式得到：

\[\kappa_N \bigg(\frac{\partial}{\partial u} \bigg) = \frac{\langle df(\frac{\partial}{\partial u}), dN(\frac{\partial}{\partial u}) \rangle}{\vert df(\frac{\partial}{\partial u}) \vert^2} = (-\sin(u), \cos(u), 0) \cdot (-\sin(u), \cos(u), 0) = 1\] \[\kappa_N \bigg(\frac{\partial}{\partial v} \bigg) = \frac{\langle df(\frac{\partial}{\partial v}), dN(\frac{\partial}{\partial v}) \rangle}{\vert df(\frac{\partial}{\partial v}) \vert^2} = (0, 0, 1) \cdot 0 = 0\]

这表示柱面在横截面内的曲率为1(圆弧)，垂直于截面的曲率为0(没有弯曲)，与我们在几何上的认知是一致的。

Principal Curvature

在曲面的切平面上存在两个相互垂直的方向\(X_1\)和\(X_2\)使得曲率分别取最小和最大值，这两个方向称为主方向(principal directions)对应的曲率称为主曲率(principal curvatures)。主曲率的重要性质如下：

\[g(X_1, X_2) = 0\] \[dN(X_i) = \kappa_i \cdot df(X_i)\]

Shape Operator

注意到\(dN\)垂直于法向\(N\)，即\(dN\)位于曲面的切平面上，因此存在唯一的线性映射\(S\)使得：

\[df(SX) = dN(X) \Leftrightarrow df \circ S = dN\]

线性映射\(S\)称为shape operator。对于主方向有：

\[S X_i = \kappa_i X_i\]

即主方向是\(S\)的特征向量，主曲率为对应的特征值。需要说明的是\(S\)未必是对称的，参数平面上\(X_1\)与\(X_2\)也不一定是相互垂直的。

此外我们还可以通过shape operator来得到曲面的第二基本形式(second fundamental form) \(\mathbf{II}\)：

\[\mathbf{II}(X, Y) = g(SX, Y) = dN(X) \cdot df(Y)\]

利用第一基本形式和第二基本形式也可以把法曲率写成：

\[\kappa_n(X) = \frac{\mathbf{II}(X, X)}{\mathbf{I}(X, X)}\]

Gaussian and Mean Curvature

通过主曲率我们可以定义曲面的高斯曲率(Gaussian curvature)\(K\)和平均曲率(mean curvature)\(H\)：

\[K = \kappa_1 \kappa_2\] \[H = \frac{1}{2} (\kappa_1 + \kappa_2)\]

其中，高斯曲率是曲面的内蕴不变量。根据Gauss-Bonnet定理(Gauss-Bonnet theorem)封闭曲面满足：

\[\int_M K dA = 2 \pi \chi\]

其中\(\chi = 2 - 2g\)为欧拉示性数(Euler characteristic)，\(g\)为曲面的亏格(genus)。对于有边界的开放曲面，Gauss-Bonnet定理仍然成立：

\[\int_M K dA + \int_{\partial M} \kappa_g ds = 2 \pi \chi\]

换言之，不管是对曲面还是对边界进行变形只要不改变曲面的拓扑性质总曲率的积分保持不变。

此外如果已知高斯曲率和平均曲率，还可以反推主曲率：

\[\kappa_1 = H - \sqrt{H^2 - K}\] \[\kappa_2 = H + \sqrt{H^2 - K}\]

Discrete Curvature

Scalar Curvatures

最后我们来讨论离散曲面上的曲率。定义顶点\(i\)的angle defect \(\Omega_i\)为\(2 \pi\)与相邻面内角和的差：

\[\Omega_i = 2 \pi - \sum_{ijk} \theta_i^{jk}\]

不难发现angle defect表示顶点的「平坦」程度，angle defect越小表示越平坦(平面上angle defect为0)。

考虑离散曲面的Gauss map，我们可以把每个相邻三角形的法向放到原点，此时：

方向的顶点都位于单位球上；
三角形的转角(dihedral angles)等于单位球法向的内角；
顶点的每一个内角等于球面上的转角；
顶点的angle defect等于球面上法向包围的面积。

对于凸多面体，对每个顶点的angle defect求和相当于计算单位球的表面积(\(4 \pi\))，也就是说angle defect的和(积分)是某种拓扑不变量。实际上angle defect就是离散曲面的高斯曲率，而且同样满足Gauss-Bonnet定理：

\[\sum_{i \in V} \Omega_i = 2 \pi \chi\]

离散曲面上的平均曲率则定义为法向夹角的一半乘以边长：

\[H_{ij} = \frac{1}{2} l_{ij} \varphi_{ij}\]

Curvature Normals

在几何处理中更为常用的曲率是曲率向量，即大小等于曲率方向与法向相同的向量。前面已经推导了面积法向(area normal)为\(\frac{1}{2} df \wedge df = N dA\)，类似地可以推导平均曲率法向(mean curvature normal)和高斯曲率法向(Gauss curvature normal)：

\[\frac{1}{2} df \wedge dN = H N dA\] \[\frac{1}{2} dN \wedge dN = K N dA\]

将面积法向\(N dA\)在顶点的dual cell上进行积分得到离散曲面的面积法向：

\[\begin{aligned} \int_C N dA &= \frac{1}{3} \int_\Omega N dA = \frac{1}{6} \int_{\partial \Omega} f \times df \\ &= \frac{1}{6} \sum_{ij \in \partial \Omega} \int_{e_{ij}} f \times df \\ &= \frac{1}{6} \sum_{ij \in \partial \Omega} \frac{f_i + f_j}{2} \times (f_j - f_i) \\ &= \frac{1}{6} \sum_{ij \in \partial \Omega} f_i \times f_j \end{aligned}\]

类似地，对平均曲率法向和高斯曲率法向进行积分可以得到离散曲面对应的曲率向量：

\[(HN)_i = \int_C H N dA = \int_C df \wedge dN = \frac{1}{2} \sum_{ij \in E} (\cot \alpha_{ij} + \cot \beta_{ij}) (f_i - f_j)\] \[(KN)_i = \int_C K N dA = \int_C dN \wedge dN = \frac{1}{2} \sum_{ij \in E} \frac{\varphi_{ij}}{l_{ij}} (f_j - f_i)\]

Steiner’s Formula

离散曲面上计算曲率的难点在于曲面不是光滑的，在很多地方曲率没有良好的定义。Steiner公式(Steiner’s Formula)提供了一种对离散曲率的近似求法：假设有一个半径为\(\varepsilon\)的球沿着凸多面体表面运动，在球运动的同时对多面体进行扩张，那么扩张后的体积可以利用半径\(\varepsilon\)进行展开

\[\text{volume}(A + B_\varepsilon) = \text{volume}(A) + \sum_{k=1}^n \Phi_k(A) \varepsilon^k\]

稍后我们会发现系数\(\Phi_k(A)\)与曲率有着深刻的联系。

首先考虑高斯曲率，显然扩张后的平面和柱面上高斯曲率为0，只有原来顶点对应的球面上存在非0的曲率。对于球面部分任一点的高斯曲率为\(K = \frac{1}{\varepsilon^2}\)，在顶点\(i\)对应的球面上进行积分得到\(K_i = (\Omega_i \varepsilon^2) \frac{1}{\varepsilon^2} =\Omega_i\)，恰好等于顶点的angle defect。同时根据Gauss-Bonnet定理，高斯曲率的积分为：

\[K(\varepsilon) = \sum_{i \in V} \Omega_i = 2 \pi \chi\]

然后考虑平均曲率，扩张后的平面部分平均曲率为0。对于柱面部分，每一点平均曲率为\(H = \frac{1}{2 \varepsilon}\)。对柱面进行积分得到每一条边的平均曲率为\(H_{ij} = \frac{1}{2 \varepsilon} l_{ij} \varphi_{ij} \varepsilon = \frac{1}{2} l_{ij} \varphi_{ij}\)。在顶点对应的球面上，平均曲率为\(H_i = (\Omega_i \varepsilon^2) \frac{1}{\varepsilon} = \Omega_i \varepsilon\)。把它们加起来得到总平均曲率：

\[H(\varepsilon) = \frac{1}{2} \sum_{ij \in E} l_{ij} \varphi_{ij} + \varepsilon \sum_{i \in V} \Omega_i\]

接下来考虑扩张后的表面积：平面部分保持不变，每个柱面面积为\(l_{ij} \varphi_{ij} \varepsilon\)，每个球面面积为\(\Omega_i \varepsilon^2\)。因此总的表面积为：

\[A(\varepsilon) = \sum_{ijk \in F} A_{ijk} + \varepsilon \sum_{ij \in E} l_{ij} \varphi_{ij} + \varepsilon^2 \sum_{i \in V} \Omega_i\]

最后来考虑扩张后的体积：

\[V(\varepsilon) = V_0 + \varepsilon \sum_{ijk \in F} A_{ijk} + \frac{1}{2} \varepsilon^2 \sum_{ij \in E} l_{ij} \varphi_{ij} + \frac{1}{3} \varepsilon^3 \sum_{i \in V} \Omega_i\]

以上的推导表明多面体的体积、表面积、曲率之间存在微分关系：

\[V(\varepsilon) \stackrel{d}{\longrightarrow} A(\varepsilon) \stackrel{d}{\longrightarrow} H(\varepsilon) \stackrel{d}{\longrightarrow} G(\varepsilon) \stackrel{d}{\longrightarrow} 0\]

同时也表明曲率描述了多面体体积扩张的速率。

Curvature Variations

假设\(f: M \rightarrow \mathbb{R}^3\)是一个光滑的映射，则曲面的体积、表面积、曲率可以表示为：

\[\text{volume}(f) = \frac{1}{3} \int_M N \cdot f dA\]
\[\text{area}(f) = \int_M dA\]
\[\text{mean}(f) = \int_M H dA\]
\[\text{Gauss}(f) = \int_M K dA = 2 \pi \chi\]

对这些量求变分可以得到很多有趣的性质。对于封闭离散曲面，其体积可以表示为：

\[\text{volume}(f) = \frac{1}{6} \sum_{ijk \in F} f_i \cdot (f_j \times f_k)\]

对于曲面上任意顶点\(i\)计算体积对于顶点位置的梯度：

\[\nabla_{f_i} \text{volume}(f) = \frac{1}{6} \nabla_{f_i} \sum_{ijk \in F} f_i \cdot (f_j \times f_k) = \frac{1}{6} \sum_{ijk \in F} f_j \times f_k = \int_{C_i} N dA\]

这说明体积对于顶点\(i\)的梯度恰好等于法向在顶点dual cell上的积分，从而得到\(\delta \text{volume}(f) = N\)。其几何意义为顶点的法向即为增大曲面体积最快的方向。

接下来考虑曲面的表面积，离散曲面的表面积为：

\[\text{area}(f) = \sum_{ijk \in F} A_{ijk}\]

在三角形上，面积对顶点的梯度为：

\[\nabla_p A = \frac{1}{2} N \times e\]

将上式带入到曲面表面积公式中，同样对顶点求梯度得到：

\[\nabla_{f_i} \text{area}(f) = \frac{1}{2} \sum_{ij \in E} (\cot \alpha_{ij} + \cot \beta_{ij}) (f_i - f_j) = (H N)_i\]

详细证明可参见Mean curvature from area。因此\(\delta \text{area}(f) = H N\)，这表示平均曲率向量的方向为增大表面积最快的方向。

对平均曲率求梯度得到：

\[\nabla_{f_i} \text{mean}(f) = \frac{1}{2} \sum_{ij \in E} \nabla_{f_i} (l_{ij} \varphi_{ij}) = \frac{1}{2} \sum_{ij \in E} (\nabla_{f_i} l_{ij}) \varphi_{ij} + l_{ij} (\nabla_{f_i} \varphi_{ij})\]

根据Schläfli公式(Schläfli formula)：

\[\sum_{ij \in E} l_{ij} (\nabla_{f_i} \varphi_{ij}) = 0\]

因此，

\[\nabla_{f_i} \text{mean}(f) = \frac{1}{2} \sum_{ij \in E} (\nabla_{f_i} l_{ij}) \varphi_{ij} = \frac{1}{2} \sum_{ij \in E} \frac{\varphi_{ij}}{l_{ij}} (f_j - f_i) = (K N)_i\]

即平均曲率的变分等于高斯曲率\(\delta \text{mean}(f) = K N\)。

综合以上可以得到变分关系：

\[\text{volume} \stackrel{\delta f}{\longrightarrow} \text{area} \stackrel{\delta f}{\longrightarrow} \text{mean} \stackrel{\delta f}{\longrightarrow} \text{Gauss} \stackrel{\delta f}{\longrightarrow} 0\]

Curvature Flow

最后我们来介绍曲率流(curvature flow)算法对网格进行处理。由于曲率和网格的体积表面积等几何量存在变分关系，因此沿曲率方向移动顶点可以改变网格的几何量，从而实现对网格的降噪（滤波）。曲率流的基本步骤为：

计算网格上每个顶点的曲率；
沿着曲率向量的方向以一定的速度移动顶点；
重复以上步骤直至收敛。

从变分的角度上看，曲率流相当于最小化网格的能量。假设网格的总能量为\(E(f)\)，通过梯度下降来最小化能量得到迭代公式：

\[\frac{d}{dt} f(t) = -\delta E(f(t))\]

在时间域和空间域上进行离散得到：

\[\frac{f_i^{k+1} - f_i^{k}}{\tau} = -\nabla_{f_i} E(f)\]

移项后得到迭代公式：

\[f_i^{k+1} = f_i^{k} - \tau \nabla_{f_i} E(f)\]

因此我们只需要选择合适的能量\(E(f)\)以及对应的梯度\(\delta E(f)\)就可以进行迭代了。曲率流算法中最出名的是平均曲率流(mean curvature flow)：取曲面总能量为曲面表面积\(E(f) = \int_M dA\)，此时梯度项为顶点的平均曲率\(\delta E = (HN)_i\)，建立迭代格式：

\[f_i^{k+1} = f_i^{k} - \frac{\tau }{2} \sum_{ij \in E} (\cot \alpha_{ij} + \cot \beta_{ij}) (f_i^k - f_j^k)\]

通过平均曲率流可以使曲面面积在每次迭代过程中不断减少，当算法收敛时得到的曲面为给定边界条件下表面积最小的曲面，称为极小曲面(minimal surface)。在日常生活中，很多建筑的流线型外观就是利用极小曲面来进行设计的。