自动关联目录：Python 实验：从 Tutte、调和映射到 LSCM

本页保留完整可运行程序，适合最后串联。第一次实验请进入模块 6：Python 网格参数化实验，分 6 步完成。

实验目标

本实验不调用现成参数化函数，而是直接实现三种方法：

Uniform Tutte：所有邻边权重为 1，凸圆边界；满足正权重条件，用来验证无翻转保证。
Cotangent harmonic：用余切权重最小化离散 Dirichlet 能量；通常角度畸变更小，但负权重时不再自动保证嵌入。
LSCM：最小化每个三角形的 Cauchy–Riemann 残差，边界自由，只固定两个锚点。

三者使用同一套输入和误差评估：

\text{拓扑检查} \to \text{局部度量} \to \text{稀疏系统} \to \text{UV} \to \text{逐三角形 Jacobian/SVD/Beltrami}.

1. 实验假设

输入网格必须满足：

单连通圆盘拓扑；
恰好一条边界环；
所有面朝向一致；
没有重复面、孤立顶点、非流形边；
每个三角形面积严格为正。

如果输入是球面或环面，先切开；如果有多个洞，需要先设计切图，不能把所有边界随意塞到同一个圆上。

2. 依赖与运行方式

python -m pip install numpy scipy matplotlib

本文代码按出现顺序放在同一个 parameterization.py 中即可运行。

3. 构造一个弯曲圆盘网格

我们从规则矩形三角剖分开始，再将内部顶点弯到三维。原始网格拓扑仍是圆盘，但三维内在边长不再等于平面矩形的边长。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from collections import Counter, defaultdict, deque
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve


def curved_rectangle_mesh(nx=24, ny=16):
    vertices = []
    for j in range(ny + 1):
        y = j / ny
        for i in range(nx + 1):
            x = i / nx
            z = 0.28 * np.sin(np.pi * x) * np.sin(np.pi * y)
            vertices.append((x, y, z))

    def vid(i, j):
        return j * (nx + 1) + i

    faces = []
    for j in range(ny):
        for i in range(nx):
            a, b = vid(i, j), vid(i + 1, j)
            c, d = vid(i, j + 1), vid(i + 1, j + 1)
            # 交错对角线，两个三角形都保持一致方向。
            if (i + j) % 2 == 0:
                faces += [(a, b, d), (a, d, c)]
            else:
                faces += [(a, b, c), (b, d, c)]

    return np.asarray(vertices, float), np.asarray(faces, int)

4. 拓扑预检与边界环

4.1 从有向面恢复边

对一致定向曲面，内部无向边出现两次且方向相反；边界边只出现一次。记录这唯一方向，就得到边界的诱导方向。

def mesh_edges_and_boundary_loops(vertex_count, faces):
    use_count = Counter()
    directed_occurrence = defaultdict(list)

    for face in faces:
        i, j, k = map(int, face)
        for a, b in [(i, j), (j, k), (k, i)]:
            edge = tuple(sorted((a, b)))
            use_count[edge] += 1
            directed_occurrence[edge].append((a, b))

    non_manifold = [e for e, count in use_count.items() if count > 2]
    if non_manifold:
        raise ValueError(f"non-manifold edges: {non_manifold[:8]}")

    inconsistent = []
    for edge, occurrences in directed_occurrence.items():
        if len(occurrences) == 2:
            (a, b), (c, d) = occurrences
            if (a, b) != (d, c):
                inconsistent.append(edge)
    if inconsistent:
        raise ValueError(f"inconsistent face orientation: {inconsistent[:8]}")

    boundary_halfedges = [
        directed_occurrence[e][0]
        for e, count in use_count.items()
        if count == 1
    ]

    outgoing = defaultdict(list)
    incoming = defaultdict(list)
    for a, b in boundary_halfedges:
        outgoing[a].append(b)
        incoming[b].append(a)

    boundary_vertices = set(outgoing) | set(incoming)
    for v in boundary_vertices:
        if len(outgoing[v]) != 1 or len(incoming[v]) != 1:
            raise ValueError(f"boundary is not a 1-manifold at vertex {v}")

    loops = []
    unused = set(boundary_halfedges)
    while unused:
        start_edge = min(unused)
        start, current = start_edge
        loop = [start]
        unused.remove(start_edge)

        while current != start:
            loop.append(current)
            nxt = outgoing[current][0]
            edge = (current, nxt)
            if edge not in unused:
                raise ValueError("boundary loop closes incorrectly")
            unused.remove(edge)
            current = nxt
        loops.append(np.asarray(loop, dtype=int))

    V, E, F = vertex_count, len(use_count), len(faces)
    report = {
        "V": V,
        "E": E,
        "F": F,
        "chi": V - E + F,
        "boundary_components": len(loops),
    }
    return sorted(use_count), loops, report


def validate_disk(vertices, faces):
    edges, loops, report = mesh_edges_and_boundary_loops(len(vertices), faces)
    if report["chi"] != 1 or report["boundary_components"] != 1:
        raise ValueError(f"expected disk topology, got {report}")

    adjacency = defaultdict(list)
    for i, j in edges:
        adjacency[i].append(j)
        adjacency[j].append(i)
    visited = {0}
    queue = deque([0])
    while queue:
        i = queue.popleft()
        for j in adjacency[i]:
            if j not in visited:
                visited.add(j)
                queue.append(j)
    if len(visited) != len(vertices):
        raise ValueError("mesh is disconnected or contains isolated vertices")

    double_areas = []
    for i, j, k in faces:
        area2 = np.linalg.norm(
            np.cross(vertices[j] - vertices[i], vertices[k] - vertices[i])
        )
        double_areas.append(area2)
    if min(double_areas) <= 1e-14:
        raise ValueError("mesh contains a degenerate triangle")
    return edges, loops[0], report

4.2 为什么 $\chi=1,b=1$ 还不完全够

对连通、可定向流形，

\chi=2-2g-b.

$\chi=1,b=1$ 推出 $g=0$ 。但如果网格不是流形或不连通，这个结论不能直接使用。因此实际预检还应做连通分量、顶点邻域和面朝向检查；上面代码展示的是本实验需要的最低版本。

5. 将边界按弧长放到单位圆

若按顶点编号等间距摆放，原网格边界采样不均匀时会制造额外切向拉伸。按三维边界边长分配圆周角：

\theta_k =2\pi \frac{\sum_{r<k}\|p_{r+1}-p_r\|} {\sum_r\|p_{r+1}-p_r\|}.

def boundary_to_circle(vertices, boundary):
    closed = np.r_[boundary, boundary[0]]
    lengths = np.linalg.norm(
        vertices[closed[1:]] - vertices[closed[:-1]], axis=1
    )
    perimeter = lengths.sum()
    if perimeter <= 1e-14:
        raise ValueError("zero boundary perimeter")

    cumulative = np.r_[0.0, np.cumsum(lengths[:-1])]
    theta = 2.0 * np.pi * cumulative / perimeter
    return np.column_stack((np.cos(theta), np.sin(theta)))

6. 两类边权

6.1 Uniform 权重

w_{ij}=1.

它不利用几何，但严格为正。配合凸边界时符合 Tutte 定理的核心假设。

6.2 Cotangent 权重

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

它离散 Dirichlet 能量，但钝角可能产生负权重。

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


def edge_weights(vertices, faces, mode):
    weights = {}

    if mode == "uniform":
        for i, j, k in faces:
            for edge in [(i, j), (j, k), (k, i)]:
                weights[tuple(sorted(map(int, edge)))] = 1.0
        return weights

    if mode != "cotangent":
        raise ValueError(f"unknown weight mode: {mode}")

    for i, j, k in faces:
        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 edge, contribution in [
            ((j, k), 0.5 * cot_i),
            ((k, i), 0.5 * cot_j),
            ((i, j), 0.5 * cot_k),
        ]:
            edge = tuple(sorted(map(int, edge)))
            weights[edge] = weights.get(edge, 0.0) + contribution
    return weights

7. 组装 Laplace 矩阵

采用正半定约定

(Lu)_i=\sum_{j\in N(i)}w_{ij}(u_i-u_j).

def laplacian_from_weights(vertex_count, weights):
    rows, cols, values = [], [], []
    diagonal = np.zeros(vertex_count)

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

    rows.extend(range(vertex_count))
    cols.extend(range(vertex_count))
    values.extend(diagonal)
    return coo_matrix(
        (values, (rows, cols)),
        shape=(vertex_count, vertex_count),
    ).tocsr()


def validate_laplacian(L, atol=1e-10):
    asymmetry = L - L.T
    symmetry_error = (
        0.0 if asymmetry.nnz == 0 else np.max(np.abs(asymmetry.data))
    )
    row_sum_error = np.max(np.abs(np.asarray(L.sum(axis=1)).ravel()))
    if symmetry_error > atol or row_sum_error > atol:
        raise ValueError({
            "symmetry_error": symmetry_error,
            "row_sum_error": row_sum_error,
        })

8. 求解 Tutte / harmonic 参数化

将顶点分为内部 $I$ 与边界 $B$ ：

L_{II}u_I=-L_{IB}u_B, \qquad L_{II}v_I=-L_{IB}v_B.

def harmonic_parameterization(vertices, faces, boundary, mode):
    n = len(vertices)
    boundary = np.asarray(boundary, dtype=int)
    is_boundary = np.zeros(n, dtype=bool)
    is_boundary[boundary] = True
    interior = np.flatnonzero(~is_boundary)

    uv = np.zeros((n, 2), dtype=float)
    uv[boundary] = boundary_to_circle(vertices, boundary)

    weights = edge_weights(vertices, faces, mode)
    L = laplacian_from_weights(n, weights)
    validate_laplacian(L)

    A = L[interior][:, interior].tocsr()
    rhs = -(L[interior][:, boundary] @ uv[boundary])
    uv[interior, 0] = spsolve(A, rhs[:, 0])
    uv[interior, 1] = spsolve(A, rhs[:, 1])

    residual = L[interior] @ uv
    diagnostics = {
        "negative_edge_weights": sum(w < -1e-14 for w in weights.values()),
        "linear_residual_inf": float(np.max(np.abs(residual))),
    }
    return uv, diagnostics

Uniform 版本若仍翻转，优先怀疑输入拓扑、边界方向或面朝向，而不是调参。Cotangent 版本翻转则不必然是代码错误，因为负权重时 Tutte 假设不成立。

9. LSCM 的局部矩阵

在单个三角形局部坐标 $(x,y)$ 中，分片线性基函数梯度为常数：

\nabla\phi_a=(g_{ax},g_{ay}).

对

u=\sum_au_a\phi_a, \qquad v=\sum_av_a\phi_a,

Cauchy–Riemann 残差为

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,

C_t= \begin{bmatrix} g_{ix}&g_{jx}&g_{kx}&-g_{iy}&-g_{jy}&-g_{ky}\\ g_{iy}&g_{jy}&g_{ky}&g_{ix}&g_{jx}&g_{kx} \end{bmatrix}.

局部能量和 Hessian：

E_t=A_t\|C_tx_t\|^2, \qquad H_t=A_tC_t^TC_t.

def local_triangle_2d(p0, p1, p2, eps=1e-15):
    e1 = p1 - p0
    length = np.linalg.norm(e1)
    normal = np.cross(e1, p2 - p0)
    area2 = np.linalg.norm(normal)
    if length <= eps or area2 <= eps:
        raise ValueError("degenerate triangle")

    x_axis = e1 / length
    normal = 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 basis_gradients(q):
    x0, y0 = q[0]
    x1, y1 = q[1]
    x2, y2 = q[2]
    area2 = (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
    if area2 <= 1e-15:
        raise ValueError("local triangle orientation is invalid")
    gradients = np.array([
        [y1 - y2, x2 - x1],
        [y2 - y0, x0 - x2],
        [y0 - y1, x1 - x0],
    ]) / area2
    return gradients, 0.5 * area2


def lscm_hessian(vertices, faces):
    n = len(vertices)
    rows, cols, values = [], [], []

    for face in faces:
        i, j, k = map(int, face)
        q = local_triangle_2d(vertices[i], vertices[j], vertices[k])
        grad, area = basis_gradients(q)
        gx, gy = grad[:, 0], grad[:, 1]
        C = np.block([
            [gx, -gy],
            [gy,  gx],
        ]).reshape(2, 6)
        Ht = area * (C.T @ C)
        dofs = np.array([i, j, k, n + i, n + j, n + k])

        for local_a, global_a in enumerate(dofs):
            for local_b, global_b in enumerate(dofs):
                rows.append(global_a)
                cols.append(global_b)
                values.append(Ht[local_a, local_b])

    return coo_matrix(
        (values, (rows, cols)), shape=(2 * n, 2 * n)
    ).tocsr()

np.block(...).reshape(2, 6) 容易因数组形状被误用。这里 gx,gy 都是一维长度 3 数组，最终两行分别是 [gx,-gy] 与 [gy,gx]。实现后应对单个三角形手工验证相似变换能量为 0。

10. 固定两个锚点求 LSCM

LSCM 能量对平移、旋转、统一缩放不敏感。固定两个不同顶点的完整 UV 坐标可消除四个实自由度。选择边界上三维距离最远的一对点：

def farthest_boundary_pair(vertices, boundary):
    points = vertices[boundary]
    diff = points[:, None, :] - points[None, :, :]
    squared = np.sum(diff * diff, axis=2)
    a, b = np.unravel_index(np.argmax(squared), squared.shape)
    return int(boundary[a]), int(boundary[b])


def lscm_parameterization(vertices, faces, boundary):
    n = len(vertices)
    H = lscm_hessian(vertices, faces)
    a, b = farthest_boundary_pair(vertices, boundary)

    fixed_dofs = np.array([a, b, n + a, n + b], dtype=int)
    fixed_values = np.array([0.0, 1.0, 0.0, 0.0])
    is_fixed = np.zeros(2 * n, dtype=bool)
    is_fixed[fixed_dofs] = True
    free = np.flatnonzero(~is_fixed)

    rhs = -(H[free][:, fixed_dofs] @ fixed_values)
    solution = np.zeros(2 * n)
    solution[fixed_dofs] = fixed_values
    solution[free] = spsolve(H[free][:, free].tocsr(), rhs)

    uv = np.column_stack((solution[:n], solution[n:]))
    residual = H[free] @ solution
    return uv, {
        "anchors": (a, b),
        "linear_residual_inf": float(np.max(np.abs(residual))),
    }

11. 逐三角形 Jacobian

原三角形局部边矩阵

X_t=[q_j-q_i\;q_k-q_i],

UV 边矩阵

U_t=[z_j-z_i\;z_k-z_i].

映射 Jacobian

J_t=U_tX_t^{-1}.

def triangle_jacobian(points_3d, uv_triangle):
    q = local_triangle_2d(*points_3d)
    X = np.column_stack((q[1] - q[0], q[2] - q[0]))
    U = np.column_stack((uv_triangle[1] - uv_triangle[0],
                         uv_triangle[2] - uv_triangle[0]))
    return U @ np.linalg.inv(X)


def parameterization_metrics(vertices, faces, uv):
    dets, conformal_K, mu_abs, areas_3d = [], [], [], []

    for i, j, k in faces:
        points = vertices[[i, j, k]]
        J = triangle_jacobian(points, uv[[i, j, k]])
        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]
        mu = 1.0 if not np.isfinite(K) else (K - 1.0) / (K + 1.0)
        area = 0.5 * np.linalg.norm(
            np.cross(points[1] - points[0], points[2] - points[0])
        )

        dets.append(det)
        conformal_K.append(K)
        mu_abs.append(mu)
        areas_3d.append(area)

    dets = np.asarray(dets)
    conformal_K = np.asarray(conformal_K)
    mu_abs = np.asarray(mu_abs)
    areas_3d = np.asarray(areas_3d)

    finite_K = conformal_K[np.isfinite(conformal_K)]
    return {
        "flipped": int(np.sum(dets < 0.0)),
        "degenerate": int(np.sum(np.abs(dets) <= 1e-14)),
        "K_median": float(np.median(finite_K)),
        "K_p95": float(np.percentile(finite_K, 95)),
        "K_max": float(np.max(finite_K)),
        "mu_p95": float(np.percentile(mu_abs, 95)),
        "log_area_p95_abs": float(
            np.percentile(np.abs(np.log(np.maximum(np.abs(dets), 1e-300))), 95)
        ),
        "face_K": conformal_K,
        "face_det": dets,
        "face_area": areas_3d,
    }

11.1 为什么同时看 $K$ 与 $\det J$

$K=1$ 表示两个方向缩放相同，但可能是统一放大、缩小或反射；
$\det J$ 给面积缩放和方向；
反射满足 $K=1$ 但 $\det J<0$ ；
$|\mu|$ 与 $K$ 等价，只是压缩到 $[0,1)$ 。

12. 可视化网格与误差热图

def plot_parameterizations(faces, results):
    fig, axes = plt.subplots(
        len(results), 2, figsize=(11, 4 * len(results)), constrained_layout=True
    )
    if len(results) == 1:
        axes = np.asarray([axes])

    for row, (name, uv, metrics, diagnostics) in enumerate(results):
        tri = mtri.Triangulation(uv[:, 0], uv[:, 1], faces)

        axes[row, 0].triplot(tri, color="black", linewidth=0.35)
        axes[row, 0].set_aspect("equal")
        axes[row, 0].set_title(
            f"{name}: flipped={metrics['flipped']}, "
            f"residual={diagnostics['linear_residual_inf']:.2e}"
        )

        color = np.log(np.maximum(metrics["face_K"], 1.0))
        image = axes[row, 1].tripcolor(
            tri, facecolors=color, shading="flat", cmap="magma"
        )
        axes[row, 1].triplot(tri, color="white", alpha=0.18, linewidth=0.2)
        axes[row, 1].set_aspect("equal")
        axes[row, 1].set_title("log conformal distortion log(K)")
        fig.colorbar(image, ax=axes[row, 1])

    plt.show()

13. 完整实验入口

def main():
    vertices, faces = curved_rectangle_mesh(nx=24, ny=16)
    edges, boundary, topology = validate_disk(vertices, faces)
    print("topology:", topology)

    experiments = []
    for name, mode in [
        ("uniform Tutte", "uniform"),
        ("cotangent harmonic", "cotangent"),
    ]:
        uv, diagnostics = harmonic_parameterization(
            vertices, faces, boundary, mode
        )
        metrics = parameterization_metrics(vertices, faces, uv)
        experiments.append((name, uv, metrics, diagnostics))

    uv, diagnostics = lscm_parameterization(vertices, faces, boundary)
    metrics = parameterization_metrics(vertices, faces, uv)
    experiments.append(("LSCM", uv, metrics, diagnostics))

    for name, _, metrics, diagnostics in experiments:
        summary = {
            key: value
            for key, value in metrics.items()
            if not key.startswith("face_")
        }
        print("\n", name)
        print("diagnostics:", diagnostics)
        print("distortion:", summary)

    plot_parameterizations(faces, experiments)


if __name__ == "__main__":
    main()

14. 结果应该怎样解释

不要预设某个方法在所有指标上最好。一般预期：

方法	主要优化目标	你应检查的现象
Uniform Tutte	正权重重心条件	凸边界下翻转应为 0，但角度畸变可能较大
Cotangent harmonic	固定圆边界的 Dirichlet 能量	角度通常改善；若有负权重，检查翻转
LSCM	Cauchy–Riemann 残差	边界可自由变形， $K$ 通常较低；不保证无翻转

如果三个方法给出完全相同结果，可能是测试网格太平坦/太对称，无法区分算法。可增大曲面弯曲、改变采样密度或换成带强凹边界的圆盘网格。

14.1 本文代码的实测结果

在 nx=24, ny=16 的默认网格上共有 425 个顶点、1192 条边、768 个三角形， $\chi=1$ 且只有一条边界。实际运行得到：

| 方法 | 负边权数 | 线性残差 $\|r\|_\infty$ | 翻转 | $K$ 中位数 | $K$ 95% | $K$ 最大值 | $|\mu|$ 95% | |---|---😐---😐---😐---😐---😐---😐---😐 | Uniform Tutte | 0 | $2.78\times10^{-15}$ | 0 | 1.318 | 2.475 | 6.151 | 0.424 | | Cotangent harmonic | 188 | $1.24\times10^{-15}$ | 0 | 1.251 | 2.283 | 6.190 | 0.391 | | LSCM | 不适用 | $1.78\times10^{-15}$ | 0 | 1.017 | 1.024 | 1.026 | 0.0116 |

这里能看到三点：

三种线性系统都被高精度求解，几何差异不是由“方程没解准”造成的。
Cotangent 权重存在负值，但本例仍没有翻转；负权重表示 Tutte 的理论保证失效，不表示结果必然失败。
LSCM 允许边界自由变化，因此能显著降低角度畸变；固定为圆的两种方法承担了额外边界约束。这个比较不能简单解读为“LSCM 在任何任务都更好”。

15. 单元测试

15.1 相似变换的共形误差应为 0

对单三角形局部坐标 $q$ ，令

U=sRq+t.

应得到 $\sigma_1=\sigma_2=s$ 、 $K=1$ 、无翻转。

def test_similarity_distortion():
    points = np.array([
        [0.0, 0.0, 0.0],
        [2.0, 0.0, 0.0],
        [0.4, 1.3, 0.0],
    ])
    q = local_triangle_2d(*points)
    angle, scale = 0.7, 2.3
    R = np.array([
        [np.cos(angle), -np.sin(angle)],
        [np.sin(angle),  np.cos(angle)],
    ])
    uv = scale * (q @ R.T) + np.array([3.0, -1.0])
    J = triangle_jacobian(points, uv)
    sigma = np.linalg.svd(J, compute_uv=False)
    assert np.allclose(sigma, [scale, scale])
    assert np.linalg.det(J) > 0.0


def test_boundary_of_boundary(vertices, faces):
    # 可复用前篇文章的 B1/B2 组装；这里强调算法不变量。
    from numpy.testing import assert_allclose
    edges, B1, B2 = boundary_matrices(len(vertices), faces)
    assert_allclose((B1 @ B2).toarray(), 0.0, atol=1e-12)

15.2 LSCM 局部 patch test

把一个三角形施加任意相似变换，代入 $C_tx_t$ ，残差应接近 0。若不为 0，通常是：

u/v 自由度排列错误；
局部坐标朝向错误；
基函数梯度符号错误；
CR 方程行的正负号不一致。

16. 数值问题与改进

16.1 病态矩阵

瘦长三角形使基函数梯度很大，矩阵条件数上升。检查：

最小角；
最短边/最长边比；
迭代求解残差；
锚点距离。

16.2 负余切权重

可做 intrinsic Delaunay edge flip，使

\alpha_{ij}+\beta_{ij}\le\pi

从而 $w_{ij}\ge0$ 。注意 extrinsic 三维 Delaunay 与 intrinsic Delaunay 不是同一个概念。

16.3 LSCM 翻转

可先用 Tutte 得到无翻转初值，再最小化 symmetric Dirichlet 能量，并在线搜索中保持 $\det J_t>0$ 。这将二次一次求解升级为非线性优化，但能显式保护局部双射。

16.4 大网格

使用稀疏 Cholesky 或预条件 CG；
对多个右端项复用分解；
不要形成稠密 $A^TA$ ；
以面积加权统计误差；
用相对残差而非只看绝对残差。

17. 从教学代码走向真实模型

真实 OBJ/PLY 流程至少需要：

合并重复顶点并移除零面积面；
统一面朝向；
检查非流形边/顶点；
提取连通分量和边界环；
对闭合或高亏格模型计算切图；
对多个 charts 做 packing；
在接缝两侧复制顶点并维护原始顶点对应；
输出纹理坐标索引，而不假设位置索引与 UV 索引相同。

18. 进一步实验

把圆边界换成正方形，比较边界角附近的 $K$ 。
将边界按顶点等间距与按弧长分配对比。
扰动一个内部顶点制造钝角，观察负余切权重与翻转。
改变 LSCM 锚点距离，记录条件数和残差。
用 analytic disk/hemisphere 构造已知共形映射，做收敛实验：网格加密时误差是否下降。
实现 symmetric Dirichlet + barrier，对比 LSCM。
在输出 UV 上计算离散 Beltrami 系数方向并画椭圆 glyph。

19. 本章结论

一个可信的几何实验必须同时有：

\boxed{ \text{理论假设检查} +\text{线性系统残差} +\text{几何误差统计} +\text{可视化} +\text{单元测试}. }

“图看起来不错”不能替代 $\det J$ 、 $K$ 、 $|\mu|$ 与残差。反过来，单一平均误差也不能替代误差热图和极端分位数。

理论回顾见：离散曲面与共形算法。下一阶段是闭合高亏格曲面的周期参数化：黎曼面、Abel 微分与参数域。