CCG 06-02:边界映射与调和参数化

本页只解决一个问题

怎样设置圆周边界、组装 Laplacian 并求内部 UV?

本页是原长篇讲义的单概念版本。第一次阅读只追踪本页定义的对象;看见未掌握的前置符号时,先沿链接退回,不要跳过定义硬背结论。

开始前

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

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

θk=2π∑r<k∥pr+1−pr∥∑r∥pr+1−pr∥.\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 权重

wij=1.w_{ij}=1.

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

6.2 Cotangent 权重

wij=12(cot⁡αij+cot⁡βij).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=∑j∈N(i)wij(ui−uj).(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 参数化

将顶点分为内部II 与边界BB:

LIIuI=−LIBuB,LIIvI=−LIBvB.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 假设不成立。

读完必须能回答

  1. 不看公式,用自己的话回答开头的核心问题。
  2. 指出本页最重要公式中每个符号的类型:标量、向量、矩阵、映射还是集合。
  3. 找出定理成立所需的假设,并举一个删除某条假设后可能失败的例子。
  4. 对最小例子手算一遍;若本页含代码,再比较手算与程序输出。

前后课程

上一课:实验准备、网格构造与拓扑预检 | 下一课:LSCM 局部矩阵与锚点