☰
共形几何与网格参数化:LSCM原理、C++实现与工程实践
2026/10/17 3:30:57 网站建设 项目流程

简介:这份计算共形几何课程的C++项目框架,专为需要动手实现几何算法的学生设计,覆盖割图(cut graph)与调和映射(harmonic map)核心算法,附带三维模型与可交互查看器,并已实测兼容Windows 10、Ubuntu 18.04、macOS 10.15.5。压缩包共373个文件,以.h头文件与.cpp源码为主,另有17个MATLAB脚本、若干位图纹理及Markdown/TXT说明,含Eigen等依赖,整体仅7.69MB。现有234人学习该资源。项目目录划分清晰,cutgraph与harmonic_map分别对应两大算法模块,3rdparty中集成MeshLib与freeglut库,读者可边阅读源码边运行查看器,直观理解割图、调和映射及Hodge分解的代码实现,并可基于框架扩展算法实验。

1. 共形几何:网格参数化里那个"保角但不保面积"的数学工具

打开一路三维扫描出来的曲面模型,想往上贴一张细节丰富的纹理,结果UV展开之后棋盘格全部拧成麻花,怎么调边界都收不干净——这是做网格参数化最常见的一种翻车现场。共形几何(Conformal Geometry)在这个场景里扮演的角色,就是给你一套"尽量不失真"的映射:它保证贴图上的直角在贴到曲面之后仍然是直角,图案不会被斜切或压扁,代价是它不保证面积。做纹理映射、四边形网格生成、形状配准的从业者,几乎都要跟这个方向打交道。这篇文章就围绕"怎么把一个普通三角网格参数化到平面,同时把角度畸变控制到最小"展开,先讲清离散共形映射的数学落点,再给出一套可复现的C++实现,最后把pin点选择、畸变度量和五个高频踩坑点逐个拆开。

2. 离散共形映射的两条路线:COT权重与最小二乘共形映射,为什么我选后者

2.1 共形条件在网格上的落点:角度不变、边长按比例缩放

连续情形下,一个复平面上的映射 f(z) 共形,当且仅当它满足柯西-黎曼方程,几何意义是"无穷小圆映到无穷小圆"。放到三角网格上,这个条件退化成逐三角形条件:每个三角形从三维曲面映射到二维参数平面时,其局部雅可比矩阵应该近似一个"旋转 + 均匀缩放"的相似变换。

具体说,对某个三角形,三维空间里的两条切边映射到二维参数平面后得到两条边向量。如果映射共形,这两条边向量应该保持互相垂直且等长,也就是雅可比矩阵满足 J^T J = λ² I。这里 λ 是该三角形处的局部缩放因子。角度不变就是由"雅可比矩阵没有各向异性"保证的;边长按比例缩放则由 λ 描述。之所以几乎所有离散共形算法都从这两条落点出发,是因为它们把几何条件变成了可计算的代数额,后面COT和LSCM只是对同一个条件的两种离散化方式。

一个容易被忽略的事实是:离散共形映射并不是真的逐点共形。它只是在每个三角形上追求"尽可能共形",误差会集中到顶点和边上。网格越密、三角形越接近正三角形,离散解越逼近连续解。这也解释了一个现象:同一套算法在不同质量的网格上结果天差地别,网格本身的质量往往比算法选择更决定最终畸变。

2.2 COT-DCP:写起来最短,但对网格质量最敏感

COT权重的思路是先把共形能量写成二次型。对每条内部边 (i,j),定义权重 w_ij = (cot α_ij + cot β_ij) / 2,其中 α_ij 和 β_ij 是这条边两侧对角的两个角度。用这个权重构造离散Dirichlet能量:

D(u) = (1/2) Σ_edges w_ij ‖u_i - u_j‖²

这个能量衡量的是"调和程度"。再引入映射后面积项 A(u),共形能量写成 E(u) = D(u) - A(u)。关键是,面积项对参数坐标也是二次的,所以整个共形能量仍然是二次型,最小化它只需要解一个对称稀疏线性系统。这就是DCP(离散共形参数化)最吸引人的地方:一条代码路径,一次线性求解,不需要迭代。

它的问题是数值上的脆弱。当三角形出现大于90度的钝角时,cot值为负,边权重变负。负权重意味着能量不再是凸的,解完可能得到一堆翻转三角形。实践中COT-DCP对扫描模型那种细长三角形密集的网格非常敏感,一个角落的退化三角形就能毁掉整个参数化。另外这个二次型如果没有额外约束,旋转和缩放自由度会留在核里,必须固定边界或固定足够多的点,边界条件又直接影响结果质量。

2.3 LSCM:把共形条件写成线性最小二乘

LSCM(最小二乘共形映射)走的是另一条路:不追求让共形能量恰好为零,而是对每个三角形建立一个误差方程,让整体的共形误差最小。做法是先把每个三维三角形投影到自身平面,得到局部二维坐标,然后要求参数平面上的三角形映射到这个局部平面时尽可能相似。

这一步会把每个三角形分解出两条线性方程。一条表达"两个切向量等长",另一条表达"两个切向量正交"。每个三角形贡献两个方程,一个F个三角形、n个顶点的网格,最后得到一个 2F × 2n 的超定稀疏线性系统。用最小二乘解它,就是LSCM的全部内容。

相比COT-DCP,LSCM有三个实际优势:第一,它不依赖边权重正负,钝角三角形不会直接给系统引入负系数的玄学问题;第二,它天然适合自由边界,边界顶点的位置由方程自行决定,不需要先指定一个凸边界;第三,它把问题统一成"解最小二乘",数值求解器的选择面更宽,调试经验也更通用。代价是需要多花一点代码去构造每三角形的方程,而且超定系统的存储和求解比单稀疏矩阵略复杂。

2.4 路线对比:何时用LSCM,何时必须在COT上加约束

从工程选型看,这两条路线不是谁取代谁,而是适用场景不同。COT-DCP代码短、矩阵小,适合实验性质地快速验证一个想法,或者在网格质量很好、三角形接近正三角形时使用。LSCM适合正式落地:扫描模型、带有大量钝角三角形的网格、需要自由边界的结果,都建议直接上LSCM。

我一般这样决策:如果网格是我自己生成的、三角形质量可控,优先尝试COT-DCP,因为它的后处理逻辑简单;如果输入网格来自外部扫描或格式转换,三角形质量完全不可控,直接走LSCM。下表是两者在我的工程视角下的核心差异。

维度COT-DCPLSCM
系统类型单个对称线性系统超定最小二乘系统
三角形角度敏感度高,钝角引发负权重中,退化三角形仍会影响
边界处理通常需要固定边界自由边界,自动决定边界形状
代码量少,核心约40行中,核心约90行
结果稳定性网格差时容易翻转网格差时相对稳定
典型场景重网格化、实验验证纹理映射、形状配准、生产管线

如果你已经确定要走COT-DCP,一个必要动作是预处理网格:把所有过于细长的三角形合并或重三角化,至少把角度小于5度的三角形处理掉,否则负权重就会变成调试黑匣子。而LSCM虽然对网格质量的容忍度更高,也不是完全免疫,狭长三角形同样会造成局部大畸变,只是症状从"翻转一片"变成"局部拉伸一大块"。

3. 用C++实现最小LSCM解算器:从网格数据到稀疏矩阵

3.1 预处理输入:顶点、三角形与固定点

LSCM解算器的输入是三样:顶点三维坐标、三角形索引、两个固定顶点及其目标二维坐标。固定点用来消除整体平移、旋转和缩放这四个刚体自由度,少一个系统就是奇异的。下面的代码用一个结构体组织输入,局部基函数负责把三维三角形压到自身平面。

// lscm_solver.cpp:最小LSCM求解器核心装配 // 依赖:常见稀疏矩阵求解库,接口与Eigen接近 #include <Eigen/Sparse> #include <Eigen/SparseCholesky> #include <vector> #include <array> #include <cmath> using Vec3 = Eigen::Vector3d; struct MeshData { int nV; // 顶点数量 std::vector<std::array<int,3>> F; // 三角形顶点索引 std::vector<Vec3> V; // 三维顶点坐标 }; // 把三维三角形压到自身平面,得到局部二维边向量 a, b // a 对应 v0->v1,b 对应 v0->v2 void localBasis(const Vec3& p0, const Vec3& p1, const Vec3& p2, double a[2], double b[2]) { Vec3 e1 = (p1 - p0).normalized(); Vec3 t = p2 - p0; Vec3 e2 = (t - e1 * t.dot(e1)).normalized(); a[0] = (p1 - p0).norm(); a[1] = 0.0; b[0] = t.dot(e1); b[1] = t.dot(e2); }

localBasis 做的事情是建立三角形自身的二维坐标系,原点在 p0,x轴沿 v0->v1,y轴垂直于 x 轴且落在三角形平面内。这个局部二维坐标很重要,因为后续的共形方程需要把三维形状和二维参数坐标放在同一个参照系里比较。注意 t.dot(e1) 可能为负,那说明三角形在三维空间里是逆时针投影,没有关系,面积由 det 捕捉。

3.2 对每个三角形生成两条共形方程

每三角形的两条方程来自"参数空间到局部平面的雅可比矩阵近似相似变换"。设局部平面上的两条边向量为 a、b,参数平面上的两条边向量为 (u1-u0, u2-u0) 和 (v1-v0, v2-v0),构造矩阵 S = [a b] 和 U = [du dv],则雅可比 J = U S⁻¹。共形条件拆成两行:J 的对角元相等,副对角元互为相反数。这两行都是参数坐标的线性组合,可以装配进全局矩阵。

// 装配 LSCM 超定方程 A * x_free = rhs // 变量排列:前 n 个是 u 坐标,后 n 个是 v 坐标 void buildLSCM(const MeshData& m, int pin0, int pin1, const std::array<double,2>& uv0, const std::array<double,2>& uv1, Eigen::SparseMatrix<double>& A, Eigen::VectorXd& rhs) { const int n = m.nV; std::vector<bool> isPin(2 * n, false); std::vector<double> pinVal(2 * n, 0.0); isPin[pin0] = true; isPin[pin1] = true; isPin[n + pin0] = true; isPin[n + pin1] = true; pinVal[pin0] = uv0[0]; pinVal[pin1] = uv1[0]; pinVal[n + pin0] = uv0[1]; pinVal[n + pin1] = uv1[1]; std::vector<int> freeCol(2 * n, -1); int nFree = 0; for (int c = 0; c < 2 * n; ++c) if (!isPin[c]) freeCol[c] = nFree++; const int nEq = (int)m.F.size() * 2; std::vector<Eigen::Triplet<double>> trips; trips.reserve(nEq * 8); rhs.setZero(nEq); // 写系数:若目标列被固定,移动到右端;否则写入自由列 auto coef = [&](int row, int col, double w) { if (isPin[col]) rhs(row) -= w * pinVal[col]; else trips.emplace_back(row, freeCol[col], w); }; for (int f = 0; f < (int)m.F.size(); ++f) { int i0 = m.F[f][0], i1 = m.F[f][1], i2 = m.F[f][2]; double a[2], b[2]; localBasis(m.V[i0], m.V[i1], m.V[i2], a, b); double detS = a[0] * b[1] - a[1] * b[0]; if (std::abs(detS) < 1e-12) detS = 1e-12; // 退化三角形保护 double si00 = b[1] / detS, si01 = -b[0] / detS; double si10 = -a[1] / detS, si11 = a[0] / detS; double w = 1.0 / std::sqrt(detS); // 按面积加权 int r0 = f * 2, r1 = f * 2 + 1; // 方程1:J00 - J11 = 0 coef(r0, i0, (-si00 - si10) * w); coef(r0, i1, si00 * w); coef(r0, i2, si10 * w); coef(r0, n+i0, ( si01 + si11) * w); coef(r0, n+i1, -si01 * w); coef(r0, n+i2, -si11 * w); // 方程2:J01 + J10 = 0 coef(r1, i0, (-si01 - si11) * w); coef(r1, i1, si01 * w); coef(r1, i2, si11 * w); coef(r1, n+i0, (-si00 - si10) * w); coef(r1, n+i1, si00 * w); coef(r1, n+i2, si10 * w); } A.resize(nEq, nFree); A.setFromTriplets(trips.begin(), trips.end()); }

这里的系数符号要仔细核对:变量 u_i 位于列 i,v_i 位于列 n+i。每个三角形两行方程里,同一个属性坐标被展开成五个非零系数,它们在行上的分布是"当前顶点取负、两个末端顶点取正",对应边向量 (u1-u0) 的差结构。权重 w 取面积平方根的倒数,目的是让面积大的三角形不要因为绝对误差大而在最小二乘里过度主导,想让面积大者权重更高时改成 w = std::sqrt(detS) 即可。

3.3 装配稀疏矩阵并解正规方程

超定系统 A x_free = rhs 的行数通常是列数的四倍左右,不能直接求逆。常见做法是把它转成正规方程 A^T A x_free = A^T rhs,再对 A^T A 做稀疏Cholesky分解。A^T A 是对称半正定矩阵,在固定点约束下正定,SimplicialLDLT这类分解器可以直接处理。

// 主流程:求解LSCM并还原完整二维坐标 Eigen::VectorXd solveLSCM(const MeshData& m, int pin0, int pin1, const std::array<double,2>& uv0, const std::array<double,2>& uv1) { const int n = m.nV; Eigen::SparseMatrix<double> A; Eigen::VectorXd rhs; buildLSCM(m, pin0, pin1, uv0, uv1, A, rhs); // 解法1:正规方程 + 稀疏Cholesky,适合中等规模网格 Eigen::SparseMatrix<double> AtA = A.transpose() * A; Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver; solver.compute(AtA); Eigen::VectorXd xFree = solver.solve(A.transpose() * rhs); // 把自由变量和固定变量合并回完整坐标向量 std::vector<bool> isPin(2 * n, false); std::vector<double> pinVal(2 * n, 0.0); isPin[pin0] = true; isPin[pin1] = true; isPin[n + pin0] = true; isPin[n + pin1] = true; pinVal[pin0] = uv0[0]; pinVal[pin1] = uv1[0]; pinVal[n + pin0] = uv0[1]; pinVal[n + pin1] = uv1[1]; Eigen::VectorXd uv(2 * n); for (int c = 0, k = 0; c < 2 * n; ++c) { if (isPin[c]) uv(c) = pinVal[c]; else uv(c) = xFree(k++); } return uv; // 前 n 项是 u,后 n 项是 v }

正规方程有两个隐患:一是 A^T A 的条件数等于 A 条件数的平方,网格退化时求解误差会被放大;二是 A^T A 的填充量可能明显高于 A 本身,大网格内存压力翻倍。如果网格顶点数超过十万,我更倾向于用迭代最小二乘,比如LeastSquaresConjugateGradient,或者直接用稀疏QR分解,虽然慢一点但数值上稳很多。

3.4 刚算完后的三项自检

算完 uv 坐标后不要立刻进下游,先做三件事。第一,统计每个三角形映射后的有符号面积,数一下有多少面积变成负数,翻转率超过 0.1% 就要回头检查固定点和网格预处理。第二,把 uv 坐标的缩放范围拉回到模型包围盒比例,检查两个固定点的目标距离是不是让整体缩放偏离常识太多。第三,随机抽 20 个三角形,手工用局部基函数重算一次共形方程残差,确认装配逻辑没有系统性符号错误。

这三项自检每次最多花十分钟,但能拦掉大多数下游纹理错误。纹理贴图出了问题往往要花一晚上排查,而结果是模型问题还是UV问题都说不清,那时候就真成玄学调试了。

4. 边界与约束:pin点怎么选、边界怎么松、畸变怎么量

4.1 固定边界、自由边界与只固定两点

参数化边界条件决定了整体由内向外的形貌。固定边界意味着预先指定边界顶点在参数平面上的位置,常见方案是映射到圆形或矩形。优势是结果边界规整,适合后续做四边形网格或者批量处理;代价是内部畸变会被不合适的边界形状放大。自由边界则让边界顶点参与最小二乘优化,边界自动收缩成与曲面曲率分布匹配的形状,畸变更小,但边界不再规整。

LSCM默认就是自由边界,这也是它处理高曲率区域时不容易撕扯的原因。只固定两个顶点是自由边界和固定边界之间的一种折中:平移、旋转、缩放被锁死,但边界顶点仍然自由优化。对纹理映射场景,我通常先跑自由边界,看边界形状是否符合下游需求,不符合再考虑固定边界。

固定边界时有个常见误解:把边界点按弧长均匀放到圆上就是最合理的。实际上,共形映射对边长比例天然敏感,弧长均匀分配会把曲率集中区域在边界上挤成一团,结果中间区域被迫拉伸。更稳的做法是先跑一次自由边界LSCM,把得到的边界 uv 坐标按角度重新采样到圆上,再用这些坐标做固定边界约束。

4.2 pin点的数量与位置怎么定

pin点的作用只是消除刚体自由度,数量两个就够,再多就会开始扭曲结果。位置选择有三个原则:第一,两个pin点不要选在相邻顶点,距离太近等价于固定了一条极短的边,相当于把整个参数化的尺度锚定在一个噪声很大的局部量上;第二,尽量避开曲率极值点,因为曲率极值处的共形方程残差本来就大,固定它会把误差压到它周围;第三,两个pin点之间的测地线路径应该大致穿过模型中心,让整体缩放分布更均匀。

目标坐标的赋值也有讲究。把uv0设为 (0,0)、uv1设为 (1,0) 是最省事的,但会导致输出 uv 处于一个固定尺度。如果下游纹理分辨率是 4096,而你输出 uv 才 0~1,会出现不必要的浮点精度损失。我一般会根据模型包围盒对角线长度先把目标距离设成与三维尺度同量级,比如三维模型直径为 10,就把两个pin的目标距离设为 8~12,这样缩放因子在 1 附近,数值更干净。

pin点还会影响翻转的位置。固定点周围容易出现负面积三角形,原因是固定的坐标和周围自由优化的顶点之间产生了冲突。如果翻转集中在某个pin点附近,把该pin点换到对面低曲率区域往往就能解决。

4.3 用三个正交指标判断参数化质量

评估参数化质量不能只盯着纹理贴图看,需要有可量化的指标。角度畸变、面积畸变、翻转率这三个指标互相独立,覆盖了共形映射的两类误差和一类致命错误。

角度畸变用每个三角形的共形能量残差表达,具体是 |J|² / area - 1 分布,或者统计三角形三个内角在映射前后的最大偏差。面积畸变用映射后三角形面积与三维面积的比值分布衡量,理想共形映射允许这个比值在不同三角形之间变化,因为缩放因子 λ 本来就逐点不同,但面积比在相邻三角形之间应该连续。翻转率就是负面积三角形占比,正常结果应该为0。

指标计算方式阈值经验说明
角度畸变共形残差或内角偏差均值均值 < 3度衡量"直角是否保留"
面积畸变面积比分布标准差标准差 < 0.3衡量"纹理密度是否均匀"
翻转率负面积三角形占比应为 0任何非零都要排查

这三个指标配合起来还有一个额外用途:判断共形参数化是否适合当前问题。做纹理映射时,面积畸变往往比角度畸变更影响观感。共形映射在这个指标上不设防,面积畸变过大就需要考虑分块参数化或者引入面积约束,这就引出第六章的进阶方向。

5. 共形参数化的五处高频翻车点:现象、原因、解法

5.1 大面积三角形翻转:病态三角形与负权重

现象:解出来的 uv 坐标里,局部一大片三角形的有符号面积为负,纹理呈现镜面翻转和褶皱,严重时整个区域像揉皱的纸。原因通常有两个:网格里存在极端钝角或狭长三角形,破坏了共形方程的局部条件;或者 pin 点固定在一个病态区域,把周围扭曲放大。解决顺序是先查网格质量,把退化三角形做局部边翻转或重新三角化,再试把 pin 点挪到曲率平缓区域。如果翻转仍然大面积存在,检查局部基函数里的 detS 是否出现极端小值,那是退化三角形混入的信号。

5.2 角度全对但纹理拉伸:共形不保面积的代价

现象:棋盘格纹理贴上去,直角全部保持,但格子的密度分布明显不均匀,高曲率区域格子密、平坦区域格子疏。原因是共形映射只约束角度,面积畸变完全由模型的曲率分布决定。这不是bug,是数学特性。解决路径有几条:最简单的是接受角度保真,在纹理坐标里额外做一次Mipmap偏移补偿;更彻底的是换成分块参数化,把高曲率区域单独切开;追求单张UV时就只能引入面积正则项,牺牲一点角度换面积均匀,具体做法见第六章。

5.3 边界绕圈与自交:pin点与边界松弛打架

现象:自由边界情况下,边界顶点在 uv 平面上绕成一个螺旋或出现自交,三角形本身没有翻转,但边界整体不合理。原因是固定边界虽然放松了,但两个 pin 点的位置仍然给边界施加了隐性约束,边界绕圈往往是解得了一个局部最优,整体旋转序没有理顺。解决方法是调整 pin 点位置,把两个 pin 点尽量放在边界两侧;如果模型有孔洞,优先把 pin 点选在孔洞边界上,用孔洞的拓扑位置锚定整体旋转。另一个有效手段是解完后再对所有 uv 坐标做一次整体 Procrustes 对齐,把旋转残差直接消除。

5.4 求解器出NaN:非流形网格吃掉了稀疏矩阵

现象:求解器返回 NaN,或者解出来的 uv 坐标溢出到十万量级。原因几乎都是输入网格有问题:重复顶点、重复边、非流形边的三角形共享、法向不一致导致局部基函数的 e2 计算出错。这类问题在稀疏矩阵里表现为某一行完全为零,系统奇异,Cholesky 分解直接失败。解决路径是先做网格预处理:查重复顶点索引、查每条无向边是否恰好被两个三角形共享、检查 e2 计算时出现除以零的位置并人工修正。很多情况下你会发现网格里藏着几千个零面积三角形,删掉它们之后 NaN 自动消失。

5.5 大网格内存爆炸:别把超定系统直接喂给稠密求解器

现象:十万顶点以上的网格,装配完成后程序内存暴涨,或者求解时间超过十几分钟。原因通常是显式构造了稠密矩阵,或者正规方程的填充因子失控。A^T A 的填充量受网格的拓扑结构影响,三角形strip类网格填充很低,但随机连接的高复杂网格可能膨胀十倍以上。解决方法是绕开正规方程,改用迭代法。另一个常见优化是把三角形索引重新排序,让相邻三角形的顶点索引聚拢,可以显著降低矩阵的带宽。对百万面级别的网格,我还会在装配前把 u 和 v 两个变量块重排成交错排列,让同一顶点的 u、v 列相邻,稀疏性会好很多。

6. 进阶:用共形因子做形状描述子,用拟共形迭代补面积

6.1 共形因子:局部拉伸变成一个标量场

共形映射在每个三角形上有一个局部缩放因子 λ,它等于参数平面三角形面积与三维空间三角形面积的比值开根号。这个 λ 拼接到所有三角形上就得到共形因子场,它刻画了曲面的局部"膨胀程度"。共形因子对形状非常敏感,两个在欧氏空间里看起来接近的模型,只要局部曲率分布不同,共形因子场就会出现明显的差异。

一个实用的技巧是:对不同模型做共形参数化后,把共形因子场归一化到相同的总尺度下,再在模型表面做颜色映射。这样不同模型之间的共形因子分布可以直接对比,用来做形状分类、局部异常检测或者质量检查。对医学表面配准这类场景,共形因子甚至可以当作一个稳定的局部特征描述子,比单纯用曲率更抗噪声,因为它是通过全局优化得到的,局部噪声会被最小二乘分摊掉。

6.2 从纯共形到低畸变:一个两段式迭代

如果下游确实需要同时控制角度和面积,一个可靠的方案是把 LSCM 当作热身,再做一轮面积正则化修正。第一轮用 LSCM 得到共形映射 uv;第二轮把面积畸变信息反馈进去,把每个三角形的面积项加权进能量。迭代几步之后,角度稍有让步但面积分布变得均匀。这个两段式思路比一开始就用复杂能量函数好调得多,因为它把两个目标分开了。

我现在每次跑完 LSCM,都会顺手把共形因子场输出一份伪彩图存盘。这不仅能快速判断参数化质量,还能在下游纹理贴图出问题时多一个排查视角。共形几何这个东西,刚接触时觉得全是数学符号,用顺手之后会发现它其实是一把很实用的工程工具,前提是你把网格质量守住、把指标量化住,然后一步步验证。希望这些经验对你有帮助。

本文还有配套的精品资源,点击获取

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询