3.3 离散微分算子

来源：Mario Botsch 等人的 Polygon Mesh Processing，第 3 章 Differential Geometry，3.3 小节 Discrete Differential Operators。

3.2 曲面 中的法向、曲率、Laplace-Beltrami 算子都假设曲面足够光滑，可以求一阶、二阶导数。但 polygon mesh 是分片线性的：每个三角形内部是平面，跨三角形边界导数不连续。

所以 3.3 的核心问题是：

如何直接从三角网格数据中，近似光滑曲面的微分几何量？

本节核心

三角网格可以看成光滑曲面的 piecewise linear approximation。

因此，离散微分算子的目标不是对网格本身硬求连续导数，而是从网格邻域中估计其背后光滑曲面的微分性质。

本节主要讨论：

• local averaging region。
• face normal 和 vertex normal。
• 分片线性函数的梯度。
• uniform graph Laplacian。
• cotangent Laplace-Beltrami。
• discrete divergence。
• discrete mean curvature。
• discrete Gaussian curvature。
• discrete curvature tensor。

3.3.1 局部平均区域

连续曲面的微分量是逐点定义的。

但在离散网格上，一个顶点本身没有足够信息。通常需要在一个局部邻域上做空间平均：

$$N(\mathbf{x})$$

如果点 $\mathbf{x}$ 是网格顶点 $v_i$，常见邻域是：

• one-ring neighborhood。
• n-ring neighborhood。
• local geodesic ball。

邻域大小会影响稳定性和精度。

邻域	优点	缺点
小邻域	保留细节，局部性强	对噪声敏感
大邻域	更稳定，抗噪	平滑过强，细节被抹掉

提示

离散微分算子的很多差异，本质上来自“在哪个局部区域上平均”和“用什么权重平均”。

常见 averaging cell

对一个中心顶点 $v_i$，书中提到几种局部面积区域。

Barycentric cell

连接三角形重心和边中点，形成顶点周围的一块局部区域。

它简单、容易计算。

Voronoi cell

用三角形外心替代重心，构造局部 Voronoi 区域。

它在误差分析上更好，但如果三角形是钝角三角形，外心可能落在三角形外部。

Mixed Voronoi cell

对钝角三角形做特殊处理，用对边中点替代外心。

这样可以让局部区域更好地平铺整个网格表面。

后续公式中经常出现的：

$$A_i$$

就是顶点 $v_i$ 对应的局部平均面积。

图源说明：根据 Botsch et al., Polygon Mesh Processing, Chapter 3 的 averaging cell 概念重绘，非原书截图。

3.3.2 法向量

很多算法需要 normal vectors：

• shading。
• curvature estimation。
• smoothing。
• remeshing。
• surface reconstruction。

三角网格上常见两类法向：

• face normal。
• vertex normal。

Face normal

对三角形：

$$T=(\mathbf{x}_i,\mathbf{x}_j,\mathbf{x}_k)$$

取两条边：

$$\mathbf{x}_j-\mathbf{x}_i, \qquad \mathbf{x}_k-\mathbf{x}_i$$

三角形单位法向为：

$$\mathbf{n}(T) = \frac{ (\mathbf{x}_j-\mathbf{x}_i)\times(\mathbf{x}_k-\mathbf{x}_i) }{ \|(\mathbf{x}_j-\mathbf{x}_i)\times(\mathbf{x}_k-\mathbf{x}_i)\| }$$

叉乘方向由三角形顶点顺序决定。

如果顶点顺序反过来，法向也会反向。

Vertex normal

顶点法向通常由 incident triangles 的 face normals 加权平均：

$$\mathbf{n}(v) = \frac{ \sum_{T\in N_1(v)} \alpha_T \mathbf{n}(T) }{ \left\| \sum_{T\in N_1(v)} \alpha_T \mathbf{n}(T) \right\| }$$

其中：

$$\alpha_T$$

是每个 incident triangle 的权重。

Vertex normal 的常见权重

Constant weights

令：

$$\alpha_T=1$$

也就是所有 incident triangles 权重相同。

优点是计算简单。

缺点是不考虑三角形大小、边长和夹角，对不规则网格可能给出不直观结果。

Area weights

令：

$$\alpha_T=|T|$$

也就是三角形面积越大，权重越大。

这种方式很高效，因为未归一化叉乘本身长度就是面积的两倍：

$$\|(\mathbf{x}_j-\mathbf{x}_i)\times(\mathbf{x}_k-\mathbf{x}_i)\| = 2|T|$$

Angle weights

令：

$$\alpha_T=\theta_T$$

其中 $\theta_T$ 是该三角形在顶点 $v$ 处的夹角。

这种方式更接近在小 geodesic disk 上平均，通常效果更好，但计算角度需要三角函数，成本略高。

信息

多数情况下，angle-weighted vertex normals 在精度和效率之间是不错的折中。

3.3.3 梯度

Laplace-Beltrami 算子定义为：

$$\Delta f=\operatorname{div}\nabla f$$

所以先要知道如何在三角网格上定义函数梯度。

设有一个分片线性函数 $f$，它在三角形三个顶点上的值为：

$$f_i,\quad f_j,\quad f_k$$

在三角形内部，用重心基函数线性插值：

$$f(\mathbf{u}) = f_i B_i(\mathbf{u}) +f_j B_j(\mathbf{u}) +f_k B_k(\mathbf{u})$$

其中 $B_i,B_j,B_k$ 是 barycentric basis functions。

因为：

$$B_i+B_j+B_k=1$$

所以：

$$\nabla B_i+\nabla B_j+\nabla B_k=0$$

函数梯度为：

$$\nabla f = f_i\nabla B_i +f_j\nabla B_j +f_k\nabla B_k$$

也可写成：

$$\nabla f = (f_j-f_i)\nabla B_j +(f_k-f_i)\nabla B_k$$

重心基函数的梯度

在一个三角形内，重心基函数是线性的，因此梯度是常量。

书中给出：

$$\nabla B_i = \frac{(\mathbf{x}_k-\mathbf{x}_j)^\perp}{2A_T}$$

其中：

• $A_T$ 是三角形面积。
• $\perp$ 表示在三角形平面内旋转 $90^\circ$。
• 梯度方向垂直于对应顶点的对边。

因此，三角形内部的 $\nabla f$ 是常量。

离散梯度不是逐顶点定义，而是通常逐三角形定义。

3.3.4 离散 Laplace-Beltrami 算子

本节讨论两个离散 Laplacian：

• uniform graph Laplacian。
• cotangent Laplacian。

Uniform Laplacian

Uniform Laplacian 只依赖网格连接关系。

对顶点 $v_i$：

$$\Delta f(v_i) = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} (f_j-f_i)$$

其中 $N_1(v_i)$ 是 $v_i$ 的一环邻居。

为什么这个公式可以看作离散的 div grad？它意味着什么？

这个公式是 Laplacian 的一种离散形式，常被称为 uniform graph Laplacian。

连续情形中：

$$\Delta f=\operatorname{div}\nabla f$$

也就是说，Laplacian 本质上可以理解为“梯度的散度”。uniform Laplacian 用最简单的邻域平均来近似这个思想。

先把公式稍微变形：

$$\Delta f(v_i) = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} (f_j-f_i)$$

由于 $f_i$ 对所有邻居都是同一个值：

$$\Delta f(v_i) = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} f_j - f_i$$

也就是：

$$\Delta f(v_i) = \operatorname{avg}_{v_j\in N_1(v_i)}(f_j) - f_i$$

所以它最直接的含义是：

Laplacian = 邻居平均值 - 当前顶点自己的值。

从离散的 $\operatorname{div}\nabla f$ 来看，可以拆成两步理解。

第一步，边上的函数差近似梯度。沿着边 $(v_i,v_j)$，函数变化量是：

$$f_j-f_i$$

它表示从 $v_i$ 走向 $v_j$ 时，函数值增加还是减少。

第二步，把所有边方向上的变化量加起来，近似散度：

$$\sum_{v_j\in N_1(v_i)}(f_j-f_i)$$

这表示中心点 $v_i$ 周围所有方向上的净变化趋势。

最后除以邻居数量：

$$|N_1(v_i)|$$

是为了取平均，减少顶点度数不同带来的影响。否则连接边更多的顶点，数值天然会更大。

因此，uniform Laplacian 可以理解为：在一环邻域里，用“所有邻居相对当前点的平均变化”来近似连续的 $\operatorname{div}\nabla f$。

这个值的符号也有直观意义：

• 如果 $\Delta f(v_i)>0$，说明邻居平均值大于当前值，当前点像一个局部低谷。
• 如果 $\Delta f(v_i)<0$，说明当前值大于邻居平均值，当前点像一个局部峰值。
• 如果 $\Delta f(v_i)=0$，说明当前值等于邻居平均值，局部达到一种平衡状态，也就是离散调和状态。

从扩散或热传导角度看，$\Delta f(v_i)$ 衡量的是当前点和周围环境的“不合群程度”。值越大，越倾向于从邻居处流入；值越小，越倾向于向邻居流出。

如果作用在坐标函数 $\mathbf{x}$ 上：

$$\Delta\mathbf{x}_i = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} (\mathbf{x}_j-\mathbf{x}_i)$$

这就是从中心点指向一环邻居平均位置的向量。

如果 f 换成顶点坐标，这个向量有什么几何意义？

如果把标量函数 $f$ 换成三维坐标函数 $\mathbf{x}$，公式变成：

$$\Delta\mathbf{x}_i = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} (\mathbf{x}_j-\mathbf{x}_i)$$

同样可以改写为：

$$\Delta\mathbf{x}_i = \operatorname{avg}_{v_j\in N_1(v_i)}(\mathbf{x}_j) - \mathbf{x}_i$$

也就是说，$\Delta\mathbf{x}_i$ 是一个从当前顶点指向一环邻居几何中心的向量。

在 Laplacian smoothing 中，可以让顶点沿这个方向移动一点：

$$\mathbf{x}_i^{new} = \mathbf{x}_i + \lambda\Delta\mathbf{x}_i$$

其中 $\lambda$ 是一个较小的步长。

这相当于让每个顶点向邻居平均位置靠拢。尖锐噪声和局部突起会被逐渐抹平，因此模型会变得更平滑。

但这也带来一个重要副作用：如果反复做这种平滑，网格通常会收缩。因为每个顶点都在向局部平均位置靠近，整体形状会逐渐被“拉紧”。

Uniform Laplacian 的问题

Uniform Laplacian 简单高效，但缺点明显：它不考虑几何位置，只看 connectivity。

因此，即使一组顶点都在同一个平面上，如果顶点分布不均匀，uniform Laplacian 也可能非零。

但连续理论中，平面区域的平均曲率应该为零。

所以 uniform Laplacian 对非均匀网格不是很好的 Laplace-Beltrami 离散化。

不过，它在某些任务中仍然有用，例如改善顶点分布的 isotropic remeshing。

Cotangent formula

更常用的 Laplace-Beltrami 离散化是 cotangent formula。

它来自混合 finite element / finite volume 的推导。

它和 uniform Laplacian 的核心区别是：uniform Laplacian 只看“谁和谁相连”，而 cotangent formula 还看三角形的具体形状。

最终公式为：

$$\Delta f(v_i) := \frac{1}{2A_i} \sum_{v_j\in N_1(v_i)} (\cot\alpha_{ij}+\cot\beta_{ij})(f_j-f_i)$$

图源说明：根据 Botsch et al., Polygon Mesh Processing, Chapter 3 的 cotangent Laplacian 权重概念重绘，非原书截图。

其中：

• $A_i$ 是顶点 $v_i$ 的局部平均面积。
• $v_j$ 是 $v_i$ 的一环邻居。
• $\alpha_{ij},\beta_{ij}$ 是边 $(i,j)$ 对面的两个角。
• $(f_j-f_i)$ 是沿边 $(i,j)$ 的函数值差。
• $\cot\alpha_{ij}+\cot\beta_{ij}$ 决定这条边对 $v_i$ 的 Laplacian 贡献有多大。

也就是说，这个公式仍然是在做“邻居值减当前值”的加权平均，只是每个邻居的权重不再相同，而是由边两侧三角形的几何形状决定。

权重为：

$$w_{ij} = \cot\alpha_{ij}+\cot\beta_{ij}$$

所以 cotangent Laplacian 可以写成：

$$\Delta f(v_i) = \frac{1}{2A_i} \sum_j w_{ij}(f_j-f_i)$$

这个 cotangent formula 应该怎么逐项理解？

可以把公式拆成三层：

$$\Delta f(v_i) = \frac{1}{2A_i} \sum_{v_j\in N_1(v_i)} w_{ij}(f_j-f_i)$$

第一层是函数差：

$$f_j-f_i$$

它表示从中心顶点 $v_i$ 沿边走到邻居 $v_j$ 时，函数值是升高还是降低。这个部分和 uniform Laplacian 一样，都是用边上的差分近似局部变化。

第二层是 cotangent 权重：

$$w_{ij}=\cot\alpha_{ij}+\cot\beta_{ij}$$

这里的 $\alpha_{ij}$ 和 $\beta_{ij}$ 不是边 $(i,j)$ 两端的角，而是这条边在左右两个三角形中的对角。

如果边 $(i,j)$ 被两个三角形共享：

$$(i,j,k), \qquad (j,i,l)$$

那么：

• $\alpha_{ij}$ 是三角形 $(i,j,k)$ 中顶点 $k$ 处的角。
• $\beta_{ij}$ 是三角形 $(j,i,l)$ 中顶点 $l$ 处的角。

直观上，cotangent 权重描述的是：边 $(i,j)$ 在局部三角形形状中有多“重要”。规则、均匀的三角形会给出比较稳定的权重；很瘦、很扁或有钝角的三角形会给出异常大甚至为负的权重。

第三层是面积归一化：

$$\frac{1}{2A_i}$$

$A_i$ 是顶点 $v_i$ 周围分配到的局部面积。除以面积的作用是把“邻域上的总变化量”变成“单位面积上的变化密度”。

那么 $A_i$ 具体怎么算？

最简单、也很常用的一种是 barycentric area。假设所有包含顶点 $v_i$ 的三角形集合为：

$$T(v_i)$$

那么：

$$A_i = \frac{1}{3} \sum_{T\in T(v_i)} |T|$$

意思是：每个 incident triangle 的面积平均分给它的三个顶点，顶点 $v_i$ 拿到相邻每个三角形面积的三分之一。

例如，$v_i$ 周围有 4 个三角形，面积分别为：

$$0.6,\quad 0.9,\quad 1.2,\quad 0.3$$

则：

$$A_i = \frac{1}{3}(0.6+0.9+1.2+0.3) = 1.0$$

更精确的做法是 Voronoi area 或 mixed Voronoi area。它不是简单地把每个三角形三等分，而是根据边中垂线、外心和钝角情况，把顶点附近真正对应的局部区域分配给 $v_i$。

在实际实现里可以这样理解：

• barycentric area：简单稳定，容易实现。
• Voronoi area：几何意义更好，但钝角三角形会麻烦。
• mixed Voronoi area：图形学中常用，对钝角三角形做特殊处理，通常更适合 cotangent Laplacian。

所以公式里的 $A_i$ 不是一个新的未知量，而是“顶点 $v_i$ 在网格表面上代表的那一小块面积”。

所以 cotangent Laplacian 的整体含义是：

在顶点 $v_i$ 周围，用三角形几何形状加权地统计邻居相对当前点的函数变化，再除以局部面积，得到曲面上的 Laplace-Beltrami 近似。

这就是它比 uniform Laplacian 更接近连续曲面 Laplace-Beltrami 的原因：它不仅知道拓扑邻接关系，还把边长、角度、局部面积这些几何信息纳入了计算。

cotangent 权重的几何含义

边 $(i,j)$ 的权重由它两侧三角形的对角决定。

如果三角形形状规则，权重通常合理。

如果三角形很瘦或出现钝角，cotangent 可能变大或变成负值。

书中指出，当：

$$\alpha_{ij}+\beta_{ij}>\pi$$

时：

$$\cot\alpha_{ij}+\cot\beta_{ij}$$

可能为负。

负权重在某些应用中会导致问题，比如参数化时可能造成三角形翻转。

注意

Cotangent Laplacian 是图形学中最常用的离散 Laplace-Beltrami，但它不是完美的。负权重和对具体三角剖分的依赖，是需要注意的两个问题。

为什么 cotangent Laplacian 更准确

相比 uniform Laplacian，cotangent formula 考虑了：

• 顶点位置。
• 三角形角度。
• 局部面积 $A_i$。

它更接近连续 Laplace-Beltrami，并广泛用于：

• surface smoothing。
• parameterization。
• shape modeling。
• curvature estimation。
• fairing。

Discrete divergence

为保持：

$$\Delta f=\operatorname{div}\nabla f$$

书中也给出离散 divergence。

设三角形上有常量向量场 $\mathbf{w}_T$，则顶点 $v_i$ 处：

$$\operatorname{div}\mathbf{w}(v_i) = \frac{1}{A_i} \sum_{T\in N_1(v_i)} \nabla B_i|_T\cdot \mathbf{w}_T A_T$$

其中：

• $\nabla B_i|_T$ 是顶点 $v_i$ 对应基函数在三角形 $T$ 内的梯度。
• $A_T$ 是三角形面积。
• $\mathbf{w}_T$ 是该三角形内的向量场。

这与离散梯度和 cotangent Laplacian 是一致的：

$$\Delta f=\operatorname{div}\nabla f$$

在离散情况下也成立。

3.3.5 离散曲率

连续曲面上有：

$$\Delta_S\mathbf{x}=-2H\mathbf{n}$$

所以对网格坐标函数 $\mathbf{x}$ 应用离散 Laplace-Beltrami，可以近似 mean curvature normal。

书中定义顶点 $v_i$ 的绝对离散平均曲率为：

$$H(v_i) = \frac{1}{2}\|\Delta\mathbf{x}_i\|$$

这里 $\Delta\mathbf{x}_i$ 是 cotangent Laplacian 作用在顶点坐标上的结果。

离散高斯曲率：角缺陷

Gaussian curvature 可以由 angle deficit 计算。

对内部顶点 $v_i$：

$$K(v_i) = \frac{1}{A_i} \left( 2\pi - \sum_{j} \theta_j \right)$$

其中：

• $\theta_j$ 是所有 incident triangles 在顶点 $v_i$ 处的角。
• $A_i$ 是顶点局部面积。

直觉是：

• 平面上，围绕一个内部点的角度和为 $2\pi$，所以 $K=0$。
• 如果角度和小于 $2\pi$，局部像凸起，$K>0$。
• 如果角度和大于 $2\pi$，局部像鞍面，$K<0$。

这来自 Gauss-Bonnet theorem。

主曲率计算

有了平均曲率：

$$H$$

和高斯曲率：

$$K$$

根据：

$$H=\frac{\kappa_1+\kappa_2}{2}$$ $$K=\kappa_1\kappa_2$$

可以解出两个主曲率：

$$\kappa_{1,2} = H\pm\sqrt{H^2-K}$$

这给出主曲率大小，但不直接给出主方向。

数值例子：角缺陷

假设内部顶点 $v_i$ 周围有 6 个三角形，每个角都是：

$$60^\circ=\frac{\pi}{3}$$

角度和：

$$6\cdot\frac{\pi}{3}=2\pi$$

因此：

$$K(v_i)=\frac{1}{A_i}(2\pi-2\pi)=0$$

这对应平坦区域。

如果周围角度和为：

$$\frac{5\pi}{3}$$

则角缺陷为：

$$2\pi-\frac{5\pi}{3} = \frac{\pi}{3}$$

如果：

$$A_i=0.5$$

则：

$$K(v_i)=\frac{\pi/3}{0.5} = \frac{2\pi}{3}$$

为正，表示局部凸起。

3.3.6 离散曲率张量

除了通过 $H,K$ 间接得到主曲率，也可以直接估计 curvature tensor。

书中介绍 Cohen-Steiner 和 Morvan 的思路：给每条 edge 定义曲率贡献。

一条 edge 的弯曲程度由它两侧三角形法向之间的 dihedral angle 决定。

设：

• $\beta(e)$ 是边 $e$ 两侧面法向的 signed dihedral angle。
• $\bar{\mathbf{e}}$ 是单位边方向。
• $\|e\cap A(v)\|$ 是边在局部区域 $A(v)$ 内的长度。

则顶点附近的曲率张量可写成：

$$C(v) = \frac{1}{A(v)} \sum_{e\in A(v)} \beta(e)\|e\cap A(v)\| \bar{\mathbf{e}}\bar{\mathbf{e}}^T$$

这个公式把边上的折角贡献累加到局部区域中。

然后可以从 $C(v)$ 的特征值和特征向量估计主曲率和主方向。

曲率张量估计的注意点

曲率张量估计依赖局部邻域大小。

常用：

• one-ring。
• two-ring。
• local geodesic disk。

如果网格不均匀，固定 n-ring 的实际几何大小可能变化很大。

此时 local geodesic disk 可能更合适。

书中也提醒，对低 valence 顶点或过小邻域，tensor averaging 可能不准确。

离散算子的整体关系

可以把本节关系整理成：

vertex values
  -> barycentric interpolation
  -> per-triangle gradient
  -> divergence of gradient
  -> cotangent Laplacian
  -> mean curvature normal

另一路：

incident triangle angles
  -> angle deficit
  -> Gaussian curvature

以及：

edge dihedral angles
  -> curvature tensor
  -> principal curvature directions

Uniform vs Cotangent Laplacian

算子	权重来源	优点	缺点
Uniform Laplacian	一环邻接数量	简单、快、只看拓扑	不适合非均匀几何
Cotangent Laplacian	对边角 cotangent + 局部面积	更接近连续 Laplace-Beltrami	可能出现负权重，依赖三角剖分

和 3.2 的关系

3.2 的连续公式：

$$\Delta_S\mathbf{x}=-2H\mathbf{n}$$

在 3.3 中变成可计算的离散流程：

$$\Delta \mathbf{x}_i = \frac{1}{2A_i} \sum_{v_j\in N_1(v_i)} (\cot\alpha_{ij}+\cot\beta_{ij}) (\mathbf{x}_j-\mathbf{x}_i)$$

然后：

$$H(v_i) = \frac{1}{2}\|\Delta\mathbf{x}_i\|$$

这就是从连续微分几何走向网格处理算法的关键桥梁。

本节记忆点

• Polygon mesh 是 piecewise linear surface，不能直接使用连续二阶导数。
• 离散微分量通常通过局部邻域平均估计。
• 顶点法向是 incident face normals 的加权平均，常用 constant、area、angle weights。
• 分片线性函数在每个三角形内部的梯度是常量。
• Uniform Laplacian：

$$\Delta f(v_i) = \frac{1}{|N_1(v_i)|} \sum_{v_j\in N_1(v_i)} (f_j-f_i)$$

• Cotangent Laplacian：

$$\Delta f(v_i) = \frac{1}{2A_i} \sum_{v_j\in N_1(v_i)} (\cot\alpha_{ij}+\cot\beta_{ij})(f_j-f_i)$$

• Discrete mean curvature：

$$H(v_i)=\frac{1}{2}\|\Delta\mathbf{x}_i\|$$

• Discrete Gaussian curvature：

$$K(v_i) = \frac{1}{A_i} \left( 2\pi-\sum_j\theta_j \right)$$

• 主曲率可由 $H,K$ 得到：

$$\kappa_{1,2}=H\pm\sqrt{H^2-K}$$

• 曲率张量可以由边的 dihedral angle 在局部区域内累加估计。

后续问题

进入第 4 章 smoothing 时可以重点关注：

1. Laplacian 如何作为平滑算子？
2. 为什么直接 Laplacian smoothing 会导致网格收缩？
3. Mean curvature flow 和 $\Delta\mathbf{x}$ 有什么关系？
4. Fairing 能量如何用 Laplacian 或 bi-Laplacian 表示？