CCG 05-09:误差指标与 Python 实现

本页只解决一个问题

如何分别验证残差、翻转、角度畸变和面积畸变?

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

开始前

19. 误差评估必须分开

指标公式测量内容
翻转det⁡Jt≤0\det J_t\le0局部不可逆
共形畸变Kt=σ1/σ2K_t=\sigma_1/\sigma_2各向异性
Beltrami$\mu_t
面积畸变log⁡(σ1σ2)\log(\sigma_1\sigma_2)局部密度变化
Dirichletσ12+σ22\sigma_1^2+\sigma_2^2总拉伸,但偏好收缩
symmetric Dirichlet∣J∣F2+∣J−1∣F2|J|_F^2+|J^{-1}|_F^2同时惩罚压缩和拉伸

全局统计不要只报告平均值。应至少给:最大值、面积加权中位数、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}

读完必须能回答

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

前后课程

上一课:顶点缩放与曲率 Jacobian | 下一课:参数化算法选择