自动关联目录：离散曲面与共形算法

本章目标

连续理论给出度量、曲率、Laplacian 和共形结构；算法必须把它们转成有限维数据。本章从三角网格出发，逐步推导：

$$\text{边长} \to \text{三角形角度与面积} \to \text{角度亏损} \to \text{余切能量} \to \text{稀疏 Laplace 矩阵} \to \text{调和/共形参数化}.$$

1. 三角剖分曲面

一个定向三角网格记为

$$M=(V,E,F),$$

其中$V$ 是顶点，$E$ 是边，$F$ 是有向三角形。几何嵌入给每个顶点位置$p_i\in\mathbb R^3$。

1.1 组合结构与几何结构要分开

• 组合结构：哪些顶点、边、面相邻；
• 分片线性几何：每条边的长度$l_{ij}$；
• 外在嵌入：顶点在$\mathbb R^3$ 中的位置。

离散内在度量只需边长，并要求每个三角形满足严格三角不等式：

$$l_{ij}+l_{jk}>l_{ki}.$$

同一组边长可以有不同三维折叠，但内在角度、面积和离散曲率相同。

1.2 半边结构

每条无向边拆成两个定向半边$h=(i\to j)$ 与$h^*=(j\to i)$。每个半边至少记录：

origin(h) : 起点
next(h)   : 同一面中的下一条半边
twin(h)   : 反向半边；边界处可为空
face(h)   : 所属面
edge(h)   : 所属无向边

常数时间可完成：

• 沿面循环：h -> next(h)；
• 穿过边：h -> twin(h)；
• 绕顶点：h -> twin(next(h))（具体公式随约定变化）。

2. 从边长恢复单个三角形

三角形$(i,j,k)$ 的边长记为

$$a=l_{jk},\quad b=l_{ki},\quad c=l_{ij}.$$

余弦定理：

$$\cos\theta_i=\frac{b^2+c^2-a^2}{2bc}.$$

Heron 公式：

$$s=\frac{a+b+c}{2}, \qquad A=\sqrt{s(s-a)(s-b)(s-c)}.$$

也可在局部二维坐标中放置

$$q_i=(0,0),\quad q_j=(c,0),$$

$$q_k=\left( \frac{b^2+c^2-a^2}{2c}, \sqrt{b^2-x_k^2} \right).$$

所有三角形内的梯度、Jacobian 和能量都可在这个局部坐标中计算。

3. 分片平坦度量与锥奇点

每个三角形内部是欧氏平面，沿边按等长方式粘合。曲率只集中在顶点，顶点邻域像一个欧氏圆锥。

内部顶点$i$ 的总角为

$$\Theta_i=\sum_{t\ni i}\theta_i^t.$$

离散高斯曲率（角度亏损）为

$$\boxed{K_i=2\pi-\Theta_i.}$$

边界顶点的边界曲率为

$$\boxed{k_i=\pi-\Theta_i.}$$

这里$K_i$ 是积分曲率，不是除以面积后的曲率密度。若需近似点值曲率，可除以 barycentric 或 Voronoi dual area。

4. 离散 Gauss–Bonnet 的推导

设网格有边界。将所有顶点曲率求和：

$$\sum_{i\in V^\circ}(2\pi-\Theta_i) +\sum_{i\in\partial V}(\pi-\Theta_i).$$

每个三角形内角和为$\pi$，所以

$$\sum_i\Theta_i=\pi F.$$

记内部顶点数$V_i$、边界顶点数$V_b$：

$$\sum K_i=2\pi V_i+\pi V_b-\pi F.$$

对三角网格，数面-边关联：

$$3F=2E_i+E_b,$$

每条边界分支是闭合多边形，所以$E_b=V_b$。代数化简得到

$$\boxed{ \sum_{i\in V^\circ}K_i+ \sum_{i\in\partial V}k_i =2\pi(V-E+F)=2\pi\chi(M). }$$

离散 Gauss–Bonnet 对任意三角剖分严格成立，不是收敛意义下的近似。

5. 分片线性有限元基函数

对每个顶点$i$，定义帽函数$\phi_i$：

$$\phi_i(p_j)=\delta_{ij},$$

并在每个三角形内线性。任意分片线性函数

$$u_h=\sum_i u_i\phi_i.$$

在每个三角形内$\nabla\phi_i$ 为常向量。

5.1 单个三角形上的梯度

在二维有向三角形$(q_i,q_j,q_k)$ 中，令$R_{90}(x,y)=(-y,x)$，有

$$\nabla\phi_i =\frac{R_{90}(q_k-q_j)}{2A},$$

以及循环置换。因为$\phi_i$ 在对边$(j,k)$ 上为 0，梯度垂直于该边；大小为“高度的倒数”。

6. 余切权重的有限元推导

连续 Dirichlet 能量的双线性型为

$$a(u,v)=\int_M\langle\nabla u,\nabla v\rangle dA.$$

离散 stiffness matrix

$$L_{ij}=\int_M\langle\nabla\phi_i,\nabla\phi_j\rangle dA.$$

在包含边$(i,j)$ 的单个三角形$(i,j,k)$ 上：

$$\int_t\langle\nabla\phi_i,\nabla\phi_j\rangle dA =-\frac12\cot\theta_k.$$

6.1 计算过程

由梯度公式：

$$\langle\nabla\phi_i,\nabla\phi_j\rangle =\frac{(q_k-q_j)\cdot(q_i-q_k)}{4A^2}.$$

分子对应角$k$ 两条边的负点积：

$$(q_k-q_j)\cdot(q_i-q_k) =-\|q_j-q_k\|\|q_i-q_k\|\cos\theta_k.$$

又

$$2A=\|q_j-q_k\|\|q_i-q_k\|\sin\theta_k.$$

乘以三角形面积$A$ 后得到$-\frac12\cot\theta_k$。

对于内部边$(i,j)$ 两侧对角$\alpha_{ij},\beta_{ij}$：

$$w_{ij}=\frac12(\cot\alpha_{ij}+\cot\beta_{ij}),$$

$$L_{ij}=-w_{ij},\quad i\ne j,$$

$$L_{ii}=\sum_{j\in N(i)}w_{ij}.$$

7. 离散 Dirichlet 能量

$$E_D(u)=\frac12u^TLu =\frac12\sum_{(i,j)\in E}w_{ij}(u_i-u_j)^2.$$

由于它来自积分$\int\|\nabla u_h\|^2$，即使某些余切权重为负，总矩阵仍为半正定（对无退化网格）。但负权重会破坏离散最大值原理。

7.1 零空间

每行和为 0：

$$L\mathbf1=0.$$

对有$c$ 个连通分支的网格，$\dim\ker L=c$。求解 Poisson 方程$Lu=b$ 时必须满足每个分支上

$$\mathbf1^Tb=0$$

或施加 Dirichlet/锚点约束。

8. 质量矩阵与强形式 Laplacian

stiffness matrix$L$ 表示弱形式。要近似点值$\Delta u$，还需质量矩阵

$$M_{ij}=\int_M\phi_i\phi_jdA.$$

单个三角形的 consistent mass matrix：

$$M_t=\frac{A_t}{12} \begin{bmatrix} 2&1&1\\1&2&1\\1&1&2 \end{bmatrix}.$$

lumped mass 把每个三角形$A_t/3$ 分给三个顶点，使$M$ 对角。

强形式近似为

$$\Delta_hu\approx -M^{-1}Lu$$

（符号取决于约定）。谱问题是广义特征问题

$$L\phi=\lambda M\phi.$$

9. 离散 Dirichlet 问题

将顶点分成内部$I$ 和边界$B$：

$$\begin{bmatrix} L_{II}&L_{IB}\\ L_{BI}&L_{BB} \end{bmatrix} \begin{bmatrix}u_I\\u_B\end{bmatrix}=0.$$

给定边界$u_B$，内部满足

$$\boxed{L_{II}u_I=-L_{IB}u_B.}$$

若权重非负，内部值是邻居的凸组合：

$$u_i=\frac{\sum_jw_{ij}u_j}{\sum_jw_{ij}}.$$

因此离散最大值原理成立。

10. Tutte 嵌入定理

设平面图三连通，把外边界固定为严格凸多边形，每个内部顶点放为邻居的严格凸组合：

$$u_i=\sum_{j\in N(i)}\lambda_{ij}u_j, \qquad \lambda_{ij}>0,quad\sum_j\lambda_{ij}=1.$$

则所得直线嵌入没有边交叉，内部面不翻转。

10.1 证明思想

1. 每个内部顶点位于邻居凸包内；
2. 任意方向投影后，内部最大值只能沿连到边界的“最大值平台”传播；
3. 若发生折叠或交叉，可构造内部极值/分离结构，与三连通及凸边界矛盾。

10.2 为什么余切权重不自动满足定理

钝角的余切为负，因此$w_{ij}$ 可能为负。Tutte 的凸组合假设失效。常见选择：

• uniform weight$w_{ij}=1$：正，但畸变较大；
• mean value coordinates：对合适网格为正；
• intrinsic Delaunay triangulation：使内部边余切和非负；
• 直接用带 barrier 的无翻转优化。

11. 调和参数化

给边界顶点二维位置$z_B=(u_B,v_B)$，分别解

$$L_{II}u_I=-L_{IB}u_B,$$

$$L_{II}v_I=-L_{IB}v_B.$$

这最小化

$$E_D(z)=E_D(u)+E_D(v)$$

且固定边界。它是调和映射，但一般不严格共形，因为固定边界形状可能迫使内部发生角度畸变。

12. 单个三角形的参数化 Jacobian

设原三角形局部二维矩阵

$$X= \begin{bmatrix} x_j-x_i&x_k-x_i\\ y_j-y_i&y_k-y_i \end{bmatrix},$$

UV 三角形矩阵

$$U= \begin{bmatrix} u_j-u_i&u_k-u_i\\ v_j-v_i&v_k-v_i \end{bmatrix}.$$

仿射映射 Jacobian 为

$$\boxed{J_t=UX^{-1}.}$$

用 SVD$J_t=R_1\operatorname{diag}(\sigma_1,\sigma_2)R_2^T$ 可计算：

• 共形畸变：$K_t=\sigma_1/\sigma_2$；
• 面积缩放：$J_t^A=\det J_t=\sigma_1\sigma_2$；
• symmetric Dirichlet：$\sigma_1^2+\sigma_2^2+\sigma_1^{-2}+\sigma_2^{-2}$；
• 翻转：$\det J_t<0$。

13. 共形能量与 Cauchy–Riemann 残差

平面映射$f=(u,v)$ 的共形条件：

$$u_x-v_y=0, \qquad u_y+v_x=0.$$

定义 least-squares conformal energy

$$E_C(f)=\frac12\int_M \left[(u_x-v_y)^2+(u_y+v_x)^2\right]dA.$$

展开：

$$\begin{aligned} 2E_C &=\int_M (u_x^2+u_y^2+v_x^2+v_y^2),dA\\ &\quad-2\int_M(u_xv_y-u_yv_x),dA. \end{aligned}$$

所以

$$\boxed{E_C=E_D-A(f),}$$

其中$A(f)$ 是有向映射面积（具体$1/2$ 系数随定义约定）。这再次给出：Dirichlet 能量减面积正是非共形误差。

14. LSCM 的单三角形组装

在局部复坐标$z_i=x_i+iy_i$ 中，线性 UV 函数$w=u+iv$ 的反全纯导数在三角形内为常数。可将每个三角形贡献写成

$$E_t=\frac1{4A_t} \left| (w_j-w_i)(z_k-z_i) -(w_k-w_i)(z_j-z_i) \right|_{\text{anti}}^2,$$

更方便的实现形式是直接离散 Cauchy–Riemann 残差。设

$$\nabla u=\sum_{a\in\{i,j,k\}}u_a\nabla\phi_a, \qquad \nabla v=\sum_av_a\nabla\phi_a.$$

三角形残差：

$$r_t= \begin{bmatrix} u_x-v_y\\u_y+v_x \end{bmatrix} =C_tx_t,$$

$$x_t=(u_i,u_j,u_k,v_i,v_j,v_k)^T.$$

局部能量

$$E_t=A_t\|C_tx_t\|^2,$$

局部 Hessian

$$H_t=2A_tC_t^TC_t.$$

把$H_t$ scatter-add 到全局矩阵$H$。因为平移、旋转和缩放不改变共形性，系统有相似变换自由度；通常固定两个不同顶点的 UV 坐标，或添加规范化约束。

15. LSCM 约束系统

将所有未知量排列为

$$x=(u_1,\ldots,u_n,v_1,\ldots,v_n)^T.$$

问题：

$$\min_x\frac12x^THx \quad\text{s.t.}\quad Cx=d.$$

KKT 系统：

$$\boxed{ \begin{bmatrix} H&C^T\\C&0 \end{bmatrix} \begin{bmatrix}x\\\lambda\end{bmatrix} = \begin{bmatrix}0\\d\end{bmatrix}. }$$

或者消去固定变量，只解自由变量：

$$H_{FF}x_F=-H_{FC}x_C.$$

锚点应尽量相隔较远，否则尺度规范化病态。

16. LSCM 不保证无翻转

LSCM 最小化二次能量，但可行域没有包含$\det J_t>0$ 约束。特别是：

• 网格不是圆盘；
• 锚点选择差；
• 三角形退化；
• 边界高度非凸；
• 能量极小值经过退化状态；

都可能出现翻转。需要显式 barrier，例如对$\det J_t$ 加

$$-\mu\log\det J_t$$

或使用本身在$\det J\to0^+$ 时发散的 symmetric Dirichlet 能量并配合 line search。

17. 离散共形等价的另一种定义：顶点缩放

参数化算法把网格映到平面；离散共形几何还研究“哪些分片平坦度量属于同一共形类”。常用 vertex scaling 定义：

$$\boxed{ \tilde l_{ij}=e^{(u_i+u_j)/2}l_{ij}. }$$

$u_i$ 是离散共形因子。改变$u$ 会改变角度和顶点曲率，但保留一个离散意义下的共形结构。

当三角不等式失效或需要跨越不同 Delaunay 区域时，执行 intrinsic edge flip。四边形翻边的边长更新由 Ptolemy 关系联系，这是离散 Möbius 几何的重要结构。

18. 曲率对共形因子的导数

离散曲率

$$K_i(u)=2\pi-\sum_{t\ni i}\theta_i^t(u).$$

在一个固定 Delaunay 区域内，其 Jacobian 与余切 Laplacian 密切相关：

$$\frac{\partial K_i}{\partial u_j} \sim -w_{ij},\qquad i\ne j.$$

因此处方目标曲率$\bar K$ 可定义势能

$$\nabla E(u)=K(u)-\bar K,$$

Hessian 是半正定 Laplace 型矩阵。固定整体尺度后，可用 Newton 法求解。这是离散 Ricci flow 章节的基础。

19. 误差评估必须分开

全局统计不要只报告平均值。应至少给：最大值、面积加权中位数、95% 分位数、翻转数和退化数。

20. Python：组装余切 Laplacian

import numpy as np
from scipy.sparse import coo_matrix


def cotangent(u, v, eps=1e-15):
    cross = np.linalg.norm(np.cross(u, v))
    if cross <= eps:
        raise ValueError("degenerate triangle")
    return float(np.dot(u, v) / cross)


def cotangent_laplacian(vertices, faces):
    n = len(vertices)
    weights = {}

    for i, j, k in np.asarray(faces, dtype=int):
        pi, pj, pk = vertices[[i, j, k]]
        cot_i = cotangent(pj - pi, pk - pi)
        cot_j = cotangent(pk - pj, pi - pj)
        cot_k = cotangent(pi - pk, pj - pk)

        for (a, b), value in [
            ((j, k), 0.5 * cot_i),
            ((k, i), 0.5 * cot_j),
            ((i, j), 0.5 * cot_k),
        ]:
            edge = tuple(sorted((a, b)))
            weights[edge] = weights.get(edge, 0.0) + value

    rows, cols, data = [], [], []
    diagonal = np.zeros(n)
    for (i, j), w in weights.items():
        rows += [i, j]
        cols += [j, i]
        data += [-w, -w]
        diagonal[i] += w
        diagonal[j] += w

    rows.extend(range(n))
    cols.extend(range(n))
    data.extend(diagonal)
    return coo_matrix((data, (rows, cols)), shape=(n, n)).tocsr()


def validate_laplacian(L, atol=1e-10):
    symmetry = np.max(np.abs((L - L.T).data), initial=0.0)
    row_sum = np.max(np.abs(np.asarray(L.sum(axis=1)).ravel()))
    if symmetry > atol or row_sum > atol:
        raise ValueError((symmetry, row_sum))

21. Python：单三角形 Jacobian 与畸变

def local_triangle_2d(p0, p1, p2):
    e1 = p1 - p0
    length = np.linalg.norm(e1)
    if length <= 1e-15:
        raise ValueError("zero edge")
    x_axis = e1 / length
    normal = np.cross(e1, p2 - p0)
    area2 = np.linalg.norm(normal)
    if area2 <= 1e-15:
        raise ValueError("degenerate triangle")
    normal /= area2
    y_axis = np.cross(normal, x_axis)
    return np.array([
        [0.0, 0.0],
        [length, 0.0],
        [np.dot(p2 - p0, x_axis), np.dot(p2 - p0, y_axis)],
    ])


def triangle_distortion(points_3d, uv):
    x = local_triangle_2d(*points_3d)
    X = np.column_stack((x[1] - x[0], x[2] - x[0]))
    U = np.column_stack((uv[1] - uv[0], uv[2] - uv[0]))
    J = U @ np.linalg.inv(X)
    sigma = np.linalg.svd(J, compute_uv=False)
    det = np.linalg.det(J)
    K = np.inf if sigma[1] <= 1e-15 else sigma[0] / sigma[1]
    return {"J": J, "det": det, "sigma": sigma, "K": K}

22. 选算法的决策表

23. 习题

计算题

1. 对一个等边三角形手算局部 stiffness matrix。
2. 对两个三角形组成的正方形组装全局$L$，验证$L\mathbf1=0$。
3. 构造一个钝角三角形，观察哪条边的余切贡献为负。
4. 对$J=\operatorname{diag}(2,1/2)$ 计算所有表中畸变指标。
5. 从顶点角度亏损推导四面体表面的总曲率。

证明题

1. 完整推导$\int_t\nabla\phi_i\cdot\nabla\phi_j=-\frac12\cot\theta_k$。
2. 证明离散 Dirichlet 能量对常数平移不变。
3. 在所有权重为正时证明离散最大值原理。
4. 推导$E_C=E_D-A$。
5. 说明为何$K=1$ 不能排除反射，必须检查行列式。

实现题

1. 实现质量矩阵、stiffness matrix，并通过线性函数 patch test。
2. 实现 uniform Tutte 与 cotangent harmonic 两种参数化，对比翻转和$K_t$。
3. 实现 LSCM 局部矩阵组装，用两个锚点求解。
4. 输出每个三角形的$\log K_t$、$\log|\det J_t|$ 热图。
5. 对同一网格执行 intrinsic Delaunay flips，比较负权重数量。

24. 本章结论

余切 Laplacian、调和参数化和 LSCM 不是三个孤立配方：

$$\boxed{ \text{PL 度量} \to \text{有限元梯度} \to \text{Dirichlet/共形能量} \to \text{稀疏二次优化} }$$

每个算法保证都附带假设：正权重、凸边界、圆盘拓扑、非退化三角形或方向保持初值。实际代码必须显式检查这些假设。

下一篇用完整代码把拓扑检查、矩阵组装、边界条件、求解与畸变统计串起来：Python 实验：三角网格共形参数化。