本页只解决一个问题

怎样从单三角形 Cauchy-Riemann 残差组装全局 LSCM？

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

开始前

• 模块前置入口：先修页面
• 本页所属模块：Python 网格参数化实验
• 集中复习版：原长篇参考讲义

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()

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))),
    }

读完必须能回答

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

前后课程

上一课：边界映射与调和参数化 ｜ 下一课：Jacobian 与畸变热图