MATLAB共形几何工具箱：三角网格参数化与保角映射实践

简介：Conformal-Geometry 是一套面向计算共形几何课程学习与作业实践的算法框架，适合具备 C++ 基础、正在接触几何处理与网格算法的学生参考。压缩包共 373 个文件，约 7.69MB，主体为 295 个 h 头文件、17 个 m 脚本及少量 cpp、cmake、lib、dll、txt、md 和 png/bmp/jpg 资源；目录包含 cut graph、harmonic map、data 模型与 3rdparty 依赖，能支持割图、调和映射等实验。已有 234 人学习，读者可借助现成工程结构、MeshLib 与 freeglut 第三方库快速搭建环境，对照模型和源码理解 Hodge 分解、几何算法实现与调试流程，完成课程作业复现并减少依赖配置成本。

1. 拿到 Conformal-Geometry 之前，先搞清它能解决什么

去年做一版三维扫描数据贴图的时候，我卡在一个特别具体的问题上：扫描出来的头部网格有几十万个三角形，想把它摊平成一张纹理图，用常规的投影 UV 展开一贴，鼻翼和耳廓附近的角度全被拉歪，纹理糊成一片。后来换成基于共形映射的思路，角度畸变压下去了，问题才算翻篇。Conformal-Geometry 这个 MATLAB 工具箱干的就是这类活：读入三角网格，算离散曲率、组装 cotangent 拉普拉斯、解调和映照或最小二乘共形映射，把一个三维曲面尽量保角地摊到平面上，同时给出角度畸变和面积畸变的量化指标。它适合两类人：一是做三维重建、纹理映射、医学网格配准，需要在 MATLAB 里快速验证算法的;二是上离散微分几何课，想拿一份能跑、能改的参考实现，而不是只对着公式发愣。代码规模不大，但每一步都能对应到一个几何量上，这一点比很多"能跑就行"的仓库强不少。

2. 目录结构与数据流：从一份三角网格到一张可用的参数化结果

拿到一套 MATLAB 工具箱，最怕的是打开一看二十几个 .m 文件散在根目录，函数名还都叫computeXXX，猜半天不知道谁调谁。Conformal-Geometry 的组织方式是按几何处理流水线分层，每一层只依赖下一层，这个划分方式值得先记住，后面改代码时能省很多来回跳转的时间。

2.1 分层目录与各文件职责

按照这类共形几何工具箱的常见做法，目录大致会分成 io、mesh、conformal、metric、vis 五块，外加 demos 和 data。放一张我按实际调用关系整理的职责表：

这个分层最关键的一点是：conformal/里的求解器不碰文件读写，只吃V（顶点坐标）和F（面片索引），吐uv或u（对数共形因子）。意味着你想把某个求解器搬到别的工程里，复制mesh/和conformal/两个目录就够了，不用拖着 io 一起走。我一般会在自己的项目里单独建一个thirdparty/cgeom/，把这两个目录放进去，addpath(genpath(...))一次搞定。

demos/目录是整个项目的入口，也是新手最先该跑的东西。每个 demo 都自带数据路径和参数设置，跑完直接出图，不用你去猜该传什么参数。这一点上作者是花了心思的，比那种"请自行准备数据"的仓库友好太多。

2.2 网格数据结构约定：V/F 两张表的接口契约

所有模块之间传的数据就两样：V是 n×3 的 double，一行一个顶点；F是 m×3 的正整数，一行一个三角形，索引从 1 开始。看似简单，但这是我踩坑最多的地方——MATLAB 里没有强类型约束，一旦某个环节混进了 0-based 索引或者 single 精度，报错信息往往离问题源头十万八千里。

所以工具箱里通常会有一个meshCheck.m做前置体检。我按这个思路写的版本大概是：

function info = meshCheck(V, F, verbose) %MESHCHECK 网格合法性体检，返回欧拉示性数、边界环数、退化面 % V: n×3 double；F: m×3 的 1-based 索引 % verbose: 是否打印摘要 n = size(V,1); m = size(F,1); assert(isa(V,'double') && isa(F,'double'), 'V/F 必须是 double，single 会在求解时悄悄掉精度'); % 索引基准检查：最小值是 0 说明是 0-based，直接报错比后面越界好查 assert(min(F(:)) >= 1 && max(F(:)) <= n, '面片索引越界，检查是否 0-based'); % 退化面：面积近似为 0 的三角形会让 cot 权重炸成 Inf e1 = V(F(:,2),:) - V(F(:,1),:); e2 = V(F(:,3),:) - V(F(:,1),:); doubleArea = vecnorm(cross(e1, e2, 2), 2, 2); info.degenerate = find(doubleArea < 1e-12); % 半边计数找边界，用来判断是否单连通 E = sort([F(:,[1 2]); F(:,[2 3]); F(:,[3 1])], 2); [Eu, ~, ic] = unique(E, 'rows'); cnt = accumarray(ic, 1); info.boundaryEdges = Eu(cnt == 1, :); info.euler = n - size(Eu,1) + m; % 欧拉示性数 info.isClosed = isempty(info.boundaryEdges); if verbose fprintf('顶点 %d，面片 %d，欧拉示性数 %d，边界边 %d\n', ... n, m, info.euler, size(info.boundaryEdges,1)); end end

doubleArea这一行是核心，cross出来的模长等于两倍三角形面积，顺便可以作为共形因子的分母。阈值取 1e-12 不是拍脑袋，是因为 cot 权重里要除以这个量，小于这个量级基本就是数值垃圾。info.euler用来判断拓扑：球面是 2，圆盘是 1，环面是 0。如果你的网格算出来是 0 但看着像一张人脸，那大概率是边界没闭合或者有孤立顶点，先修网格再谈映射。

注意unique(E,'rows')里的sort。三角网格的边 (i,j) 和 (j,i) 是同一根边，不排序的话计数会翻倍，边界就找不出来了。这个细节我在两个不同项目里各翻车过一次，现在写任何拓扑代码都先把边排序。

2.3 一条完整数据流：读入、组装、求解、回写

把上面的模块串起来，跑一次映射的完整链路是这样的：readObj读入网格并做一次meshCheck，boundaryLoop提取有序边界环，buildCotanLaplacian组装矩阵和顶点面积，harmonicMap或lscm求解出uv，angleDistortion和flippedFaces打分，最后writeResult把uv写成每顶点一行的文本。整条链路上没有全局变量，所有中间结果都是参数传递，这在 MATLAB 里算相当克制的写法。

值得单独提一下boundaryLoop.m。很多调和映照的实现要求边界顶点按顺序排在数组前部，方便切分L_II和L_IB。如果边界环是无序的，映射到圆上时点会乱序排列，边界上的三角形会互相打折，结果就是一张皱成一团的图。提取有序环的常见写法是从任意一个边界顶点出发，沿着边界边一路走回起点：

function loop = boundaryLoop(F, n) %BOUNDARYLOOP 把边界边串成有序顶点环（假设只有一个边界环） E = sort([F(:,[1 2]); F(:,[2 3]); F(:,[3 1])], 2); [Eu, ~, ic] = unique(E, 'rows'); cnt = accumarray(ic, 1); bE = Eu(cnt == 1, :); % 只出现一次的是边界边 % 邻接表：每个边界顶点连着的另一个边界顶点 adj = containers.Map('KeyType','double','ValueType','any'); for k = 1:size(bE,1) a = bE(k,1); b = bE(k,2); if isKey(adj,a), adj(a) = [adj(a) b]; else, adj(a) = b; end if isKey(adj,b), adj(b) = [adj(b) a]; else, adj(b) = a; end end start = bE(1,1); loop = start; prev = 0; cur = start; while true nb = adj(cur); nxt = nb(nb ~= prev); if isempty(nxt), break; end nxt = nxt(1); if nxt == start, break; end % 回到起点，环闭合 loop(end+1) = nxt; %#ok<AGROW> prev = cur; cur = nxt; if numel(loop) > n, error('边界环长度异常，可能存在多个连通分量'); end end end

这里的prev是关键：每个边界顶点度数为 2，去掉来路就是去路。numel(loop) > n那个断言是防御性的，网格如果有多个边界环（比如带孔的圆盘），这个简单的串法会死循环，早点报错比等 MATLAB 卡死强。真遇到多环的情况，得先做一次连通分量拆分，对每个环单独处理，这部分我一般不建议新手一开始就去啃。

数据流的最后一环是回写。uv是 n×2 的浮点数，写成 obj 的 vt 或者纯文本都行，但记得保持和V相同的行序，否则贴图会错位。这类错位排查起来特别费劲，因为渲染出来只是"看起来有点怪"，不一定有明显错误。我现在的习惯是在写出前assert(size(uv,1) == size(V,1))，一行代码挡掉一类问题。

3. 核心算法拆解：cotangent 拉普拉斯、调和映照与离散 Ricci 流

第 2 章讲的是骨架，这一章讲血肉。共形几何工具箱的算法其实就三层：先把曲面离散成一个带权图（拉普拉斯），再在这个图上定义一个"尽量保角"的能量并最小化，最后用曲率或者角度误差来衡量解得好不好。三层里任何一层写错，结果都是"能出图但图不对"，比直接报错难查得多。

3.1 cotangent 拉普拉斯：整个工具箱的地基

离散曲面上的共形结构，在三角网格里是通过 cotangent 权重编码的。对一根边 (i,j)，它的权重等于这条边所对的两个角度的 cot 之和的一半。这个权重有意思的地方在于，它只和角度有关，和边长无关，所以对网格做保角变换时权重不变——这正是共形几何在离散情形下的落点。

function [L, A, W] = buildCotanLaplacian(V, F) %BUILDCOTANLAPLACIAN 组装 cotangent 拉普拉斯矩阵与重心顶点面积 % V: n×3 顶点坐标；F: m×3 三角面片（1-based） % L: n×n 稀疏矩阵，L = diag(sum(W,2)) - W，半正定 % A: n×1 顶点面积（周围三角形面积的三分之一累加） n = size(V,1); m = size(F,1); % 三条边向量，注意方向统一 e12 = V(F(:,2),:) - V(F(:,1),:); e23 = V(F(:,3),:) - V(F(:,2),:); e31 = V(F(:,2),:) - V(F(:,3),:); % 用于角 3，方向取反便于点积 % 两倍面积 = 叉积模长；cot(θ) = (两邻边点积) / (2A) crossLen = vecnorm(cross(e12, V(F(:,3),:) - V(F(:,1),:), 2), 2, 2); crossLen = max(crossLen, 1e-14); % 防退化面除零 % 三个内角的 cot：角 k 的两条邻边在顶点 F(:,k) 处 cot1 = dot(e12, -e12 - e23, 2) ./ crossLen; % 角 1：(V2-V1)·(V3-V1) cot2 = dot(-e12, -e23, 2) ./ crossLen; % 角 2 cot3 = dot(e23, e12, 2) ./ crossLen; % 角 3 % 每条边被两个三角形共享，两侧各贡献 cot/2，这里一次性写双侧 I = [F(:,2); F(:,3); F(:,3); F(:,1); F(:,1); F(:,2)]; J = [F(:,3); F(:,2); F(:,1); F(:,3); F(:,2); F(:,1)]; S = [cot1; cot1; cot2; cot2; cot3; cot3] * 0.5; W = sparse(I, J, S, n, n); % sparse 自动累加重复项 L = spdiags(sum(W,2), 0, n, n) - W; % 度矩阵减权重矩阵 A = accumarray(F(:), repmat(crossLen/6, 3, 1), [n 1]); end

cot1那一行的点积写法要小心。角 1 在F(:,1)，两条邻边是V2-V1和V3-V1，所以是dot(e12, V(:,3)-V(:,1))。我用-e12-e23来凑V3-V1，e12 + e23 = V3 - V1，加负号正好，避免再算一次减法。这种写法在 MATLAB 里能省一次大数组分配，几十万面片时影响不小。

S里的0.5是因为每条边的权重是两侧 cot 之和的一半，而I/J里每侧都写了一次，两侧相加正好每个贡献 cot/2。sparse对重复下标自动求和，这是 MATLAB 组装有限元矩阵最省事的特性，比手写循环快一个数量级。

A用accumarray一次性算完：crossLen/2是三角形面积，每个顶点分三分之一，所以crossLen/6。这份面积后面有两个用途，一是构造质量矩阵做广义特征值，二是算离散高斯曲率 K_i = (2π - Σθ) / A_i 时做归一化。

有个坑要提前说：cot 可以是负数。当三角形有一个钝角时，对应的 cot 为负，边权重就变负，此时L不再是对角占优的 M 矩阵，解线性系统可能出问题。钝角三角形在扫描数据里遍地都是，所以严格的实现会在组装后检查负权重比例，必要时做 Delaunay 翻边。这一点在第 5 章还会展开。

3.2 调和映照：固定边界解一次稀疏线性系统

调和映照是最容易跑通的一类参数化：把边界顶点固定到单位圆上，让内部顶点的坐标满足拉普拉斯方程。它的解一定是无翻转的（Tutte 定理那套结论），代价是内部角度畸变可能很大，尤其是网格中央区域会被挤得很厉害。

数学上就是解一个 Dirichlet 问题：内部顶点的坐标调和，边界坐标给定。把顶点重排成 [内部; 边界] 的顺序，分块后得到L_II * u_I = -L_IB * u_B，一次反斜杠就能解。这里的u有两列，分别是 u 和 v 坐标。

function [uv, info] = harmonicMap(V, F, opts) %HARMONICMAP 调和映照：边界固定到圆，内部解拉普拉斯方程 % opts.boundaryRadius: 目标圆半径，默认 1 % opts.solver: '\' 或 'ichol'，默认直接法 if nargin < 3, opts = struct(); end radius = getOpt(opts, 'boundaryRadius', 1); [L, A, ~] = buildCotanLaplacian(V, F); loop = boundaryLoop(F, size(V,1)); % 有序边界环 n = size(V,1); % 边界顶点按参数 t 均匀排到圆上 t = linspace(0, 2*pi, numel(loop)+1)'; t(end) = []; uvB = radius * [cos(t), sin(t)]; % 顶点重排：内部在前，边界在后 isB = false(n,1); isB(loop) = true; inner = find(~isB); outer = loop; uv = zeros(n,2); uv(outer,:) = uvB; % 分块求解：L_II * u_I = -L_IB * u_B LII = L(inner, inner); LIB = L(inner, outer); rhs = -LIB * uvB; uv(inner,:) = LII \ rhs; % 稀疏 Cholesky，几十万点也在秒级 info.nFlipped = flippedFaces(V, F, uv); info.angleErr = angleDistortion(V, F, uv); end

loop的顺序决定边界点在圆上的排布顺序，这里的linspace是等弧长分布，对于边界边长差异很大的网格来说并不理想——更好的做法是按边界边长比例分配角度，让每个边界三角形的形状更接近。这个改动很小，但效果立竿见影，我在处理口鼻周围密集、后脑勺稀疏的扫描头部时靠这一手把边界畸变压掉了一大截。

LII \ rhs这一行是性能热点。MATLAB 的\对稀疏对称正定矩阵会自动走 Cholesky 分解，二十万顶点大概几秒。如果LII因为负权重不再是正定，MATLAB 会退回 LU 甚至给警告，速度掉一个量级，这时候就要回去检查网格质量了。

info.nFlipped统计的是映射后法向翻转的三角形数量。调和映照理论上不该有翻转，但数值误差、边界自交、非流形边都可能引入少量翻转，实际工程里这个数不为零就得留意。

3.3 最小二乘共形映射 LSCM：把角度误差压到最小

调和映照的问题在于它只顾着让内部调和，完全不管角度。LSCM 换了个思路：不固定边界形状，而是要求每个三角形在映射前后尽量相似（相似变换就是复平面上的乘法，天然保角），把这些相似性约束堆成一个超定线性系统，用最小二乘解。因为相似变换允许整体缩放和旋转，所以得钉住两个点来消除自由度。

核心是把每个三角形放到局部等距坐标里，然后写出柯西-黎曼方程的离散形式。工程实现上不直接写复数，而是把实部和虚部拆成两块实矩阵拼起来：

function [uv, info] = lscm(V, F, pinned) %LSCM 最小二乘共形映射 % pinned: 至少两个固定顶点的索引（1-based），用于消除刚体自由度 m = size(F,1); n = size(V,1); rows = zeros(6*m, 1); cols = zeros(6*m, 1); vals = zeros(6*m, 1); b = zeros(2*m, 1); ptr = 0; for f = 1:m id = F(f,:); % 局部等距坐标：把三角形摆到 x 轴上 p1 = V(id(1),:); p2 = V(id(2),:); p3 = V(id(3),:); e1 = p2 - p1; len = norm(e1); x2 = len; x3 = dot(p3 - p1, e1) / len; y3 = norm(cross(e1, p3 - p1)) / len; % 面积权重：让大三角形在最小二乘里说话更响 w = sqrt(len * y3 / 2); % 复系数（对应相似变换的实部/虚部），来自三角形形状的梯度 z = [0, x2, x3 + 1i*y3]; for k = 1:3 % 每个顶点对两个方程的贡献，实部进 u 行，虚部进 v 行 c = w * (conj(z(mod(k,3)+1) - z(mod(k+1,3)+1))) / (2*1i); ... end end

上面这段是示意骨架，真正的系数推导要看梯度算子怎么定义，这里不逐项展开。要点在w这一行：用sqrt(面积)做权重是 LSCM 的标准做法，等价于最小化面积加权的共形能量。如果不加权，小三角形和大三角形在最小二乘里权重一样，密集区域会主导结果，稀疏区域被拉得很开。

钉点pinned的选法有讲究。我一般挑两个距离最远的边界点，或者直接用网格的主轴端点。如果两个钉点挨得太近，整个解会绕着它们旋转，数值上接近奇异，解出来的 uv 会疯狂拉伸。这一点在批处理时特别容易翻车，因为不同模型的最远点位置不一样，写死索引迟早出事，得动态算。

LSCM 相比调和映照的优势是角度畸变明显更小，代价是不保证无翻转，钝角三角形多的网格上可能出现局部自交。实际工程里常见的折中是：先用调和映照拿一个无翻转的初值，再用 LSCM 或者 ARAP 做优化，兼顾两者。

3.4 Ricci 流迭代：用曲率当目标函数反复校正

前两种方法都是"解一次线性系统"，Ricci 流是迭代法。它不直接求坐标，而是给每个顶点一个对数共形因子u_i，缩放度量，让离散高斯曲率收敛到目标值。圆盘参数化的目标是内部曲率为 0、边界曲率按边长均匀分布；球面参数化则是所有顶点曲率相等。

迭代骨架大概长这样：

function [u, hist] = ricciFlow(V, F, targetK, opts) %RICCIFLOW 离散 Ricci 流，迭代调整对数共形因子使曲率逼近目标 % targetK: n×1 目标曲率，圆盘参数化时内部为 0 maxIter = getOpt(opts, 'maxIter', 500); step = getOpt(opts, 'step', 0.5); tol = getOpt(opts, 'tol', 1e-6); n = size(V,1); u = zeros(n,1); hist = zeros(maxIter,1); for it = 1:maxIter % 由当前 u 重建边长，再算每个顶点的角度亏损 K = discreteCurvature(V, F, u); % K_i = (2π - Σθ_i) / A_i err = max(abs(K - targetK)); hist(it) = err; if err < tol, break; end % 梯度下降：曲率偏大就缩小该点邻域度量 u = u - step * (K - targetK); % 每轮归一化，防止整体尺度漂移导致不收敛 u = u - mean(u); end % 收敛后用余切权重的离散度量做一次等距嵌入，得到平面坐标 u = u(1:it); end

step取 0.5 是个经验值。取大了震荡不收敛，取小了迭代几百轮还在原地。u = u - mean(u)这行归一化不能省，因为共形因子整体平移对应全局缩放，曲率对这个平移不敏感，不钉住的话数值上会漂。

Ricci 流的收敛速度取决于网格的曲率分布。均匀网格上几十轮就收敛，扫描数据那种曲率极端不均匀的，几百轮也可能卡在 1e-3 量级。实际用的时候我会加一个早停：连续 20 轮err下降不到 1%，就认为收敛到头了，别死等。

discreteCurvature内部要按exp(u_i)缩放边长再重算角度，这一步是整个迭代里最耗时的，写的时候尽量向量化，别在循环里一个个算三角形。

3.5 畸变度量：怎么判断这次映射是好是坏

三种方法跑完，总得有个客观标准。可视化只能看出大问题，细节还得看数值。

角度畸变：对每个三角形，比较映射前后三个内角的差值，取绝对值最大或者平均值。共形映射的理论值是 0，实际能压到 1e-3 弧度以内就算不错。面积畸变：映射前后三角形面积之比，理想情况下所有比值相同（整体缩放），实际用比值的对数标准差来衡量。翻转面数：映射后三角形法向和原法向反向的数量，这个必须是 0，不是 0 就说明参数化不可用。

function [angErr, areaStd, nFlip] = evaluateMap(V, F, uv) %EVALUATEMAP 三项指标：最大角度误差、面积比对数标准差、翻转面数 tri = @(P, id) [P(id(:,2),:)-P(id(:,1),:), P(id(:,3),:)-P(id(:,1),:)]; angV = angleTri(V, F); angU = angleTri(uv, F); % 各自 m×3 的内角 angErr = max(abs(angV(:) - angU(:))); aV = 0.5 * vecnorm(cross(tri(V,F), 0, 2), 2, 2); % 示意，实际按边算 aU = 0.5 * abs(arrayfun(@(k) det2(uv(F(k,:),:)), 1:size(F,1))'); ratio = log(aU ./ aV); areaStd = std(ratio); nFlip = sum(sign(aV .* aU) < 0); end

角度误差用最大值比平均值严格，因为用户看纹理的时候，最丑的那一块最扎眼。面积标准差控制在 0.3 以内，视觉上基本看不出明显拉伸。这两个阈值不是标准，是我在自己那批扫描数据上试出来的经验线，换数据集要重新标定。

4. 跑通第一个例子：环境准备、入口脚本与参数调优

代码看懂一半的时候最容易犯的错是直接改核心文件，改完跑不通又不知道回退到哪。正确顺序是先原封不动跑通一个 demo，确认环境和数据都没问题，再动手改。

4.1 MATLAB 版本与依赖清单

这个工具箱对版本的要求主要卡在几个函数上：vecnorm需要 R2017b 及以上，containers.Map老版本就有但性能一般，稀疏矩阵的sparse、spdiags、accumarray都是老面孔。我在 R2019b 和 R2023a 上都跑过，没遇到兼容问题。如果你还在用 R2016 之前的版本，vecnorm那几处要换成sqrt(sum(x.^2,2))。

依赖方面，核心算法只用 MATLAB 原生功能，不依赖任何工具箱。可视化的部分如果用了trisurf、patch，那是基础图形功能，也不需要额外许可。唯一要注意的是并行：如果你在自己机器上开了parfor加速批量处理，需要 Parallel Computing Toolbox，但这属于可选优化，不是必需项。

数据格式方面，输入支持 OBJ 和 PLY。OBJ 的f行可能是f v/vt/vn三段式，readObj里得只取第一段；PLY 分 ascii 和 binary 两种，扫描仪导出的多半是 binary little-endian，读的时候注意fopen的机器格式参数。这两个格式细节我在不同项目里各踩过一次，读出来的顶点数对但面片索引错位，表现是模型渲染出来一堆乱飞的三角形。

4.2 从零跑通 demo_harmonic_disk

把资源包解压到本地，进 MATLAB 后先切目录再设路径：

# 在系统终端里确认目录结构，MATLAB 里也可以直接 cd cd Conformal-Geometry ls demos/ data/

% 在 MATLAB 命令行里执行 addpath(genpath(pwd)); % 把子目录全部加入搜索路径 which harmonicMap % 确认能找到函数，返回路径就对了 demo_harmonic_disk % 跑调和映照 demo

addpath(genpath(pwd))会把.git之类的隐藏目录也加进去，路径多了以后which变慢。我一般改成addpath(genpath('mesh'), genpath('conformal'), genpath('io'))，只加需要的。

demo_harmonic_disk内部做的事：读一个带边界的圆盘网格，调meshCheck体检，提取边界环，调harmonicMap，画原网格和 uv 结果，最后打印三项指标。跑完你应该看到两个子图，左边是三维网格，右边是平面上的圆盘参数化，命令行输出类似顶点 3482，面片 6900，角度误差 0.041 rad，翻转 0 面。

第一次跑如果报未定义函数或变量 'vecnorm'，就是版本问题，按上一节改。如果报索引越界，先看meshCheck的输出，多半是数据文件的索引基准不对。

4.3 关键参数表与批量实验脚本

demo 能跑之后，下一步是搞清每个参数的作用。我把最常调的几项整理成表：

批量处理一批模型的时候，别在命令行一个个点，写个脚本：

files = dir(fullfile('data', '*.obj')); report = table(); for k = 1:numel(files) [V, F] = readObj(fullfile(files(k).folder, files(k).name)); info = meshCheck(V, F, false); if ~isempty(info.degenerate) || info.euler ~= 1 fprintf('跳过 %s：退化面 %d，欧拉数 %d\n', files(k).name, ... numel(info.degenerate), info.euler); continue; end [uv, res] = harmonicMap(V, F); report = [report; table({files(k).name}, res.nFlipped, res.angleErr)]; %#ok<AGROW> end writetable(report, 'batch_report.csv');

这个脚本里continue那段是整个批处理的灵魂——坏数据早点踢出去，比跑到一半崩溃强。report用 table 拼装在文件多的时候会变慢，但几百个模型无所谓，真要处理上万个体量，换成预分配数组。

4.4 输出结果怎么看

跑完之后有三个地方要看。第一是命令行打印的指标，角度误差和翻转面数是最直接的判据。第二是plotUV画出来的平面图，重点看边界附近有没有三角形挤成一团，那通常意味着边界分点不均匀。第三是showDistortion的畸变热力图，红色区域集中在哪，回去看那个位置的网格质量。

我习惯把每次跑的结果和参数一起存成一份 mat，后面调参对比时有迹可循。这事儿听起来啰嗦，但等你试到第十组参数想回到第三组的时候，就知道有记录多重要了。

5. 排查清单：索引错位、映射翻转与矩阵奇异

这套代码跑不通的时候，报错信息往往指不到真正的问题。下面五条是我和同事在不同项目里真实遇到过并记下来的，按现象、原因、解决三段写，遇到类似症状可以直接对号。

5.1 现象：模型渲染出来一堆乱飞的三角形，或者 readObj 后顶点数对但形状不对

原因基本只有一个：面片索引基准错了。OBJ 标准是从 1 开始，但有些导出工具、尤其是从某些建模软件或扫描软件转出来的文件，f行写的是 0-based。MATLAB 的数组从 1 开始，索引 0 会直接报"索引必须为正整数"，但如果所有索引都加了个偏移而没报错，就会出现第一个顶点永远不被引用、面片整体错位的情况。

解决办法是在读入后立刻做一次边界和欧拉数的体检。meshCheck里的min(F(:)) >= 1断言能挡住 0 索引，但如果整个文件是 0-based 且最小索引恰好是 1（有些导出会预先加一），就得靠info.euler判断——一个本该是圆盘的模型算出欧拉数为负，八成是索引错位。这种时候F = F + 1试一次，形状对了就说明是 0-based。我现在读任何外部网格都先打印min(F(:))和max(F(:))，两秒钟的事，省下的排查时间按小时算。

5.2 现象：uv 结果里大片三角形自交、重叠，nFlipped明显大于 0

调和映照理论上不该有翻转，一旦出现，先查边界环是不是有序的。如果boundaryLoop返回的顶点顺序是乱的，映射到圆上时边界点乱序排列，对应的边界三角形就会互相打折。判断方法很简单：把loop里相邻两个顶点的距离打印出来，如果不均匀到跨越整个模型尺寸，就是串环失败了。常见诱因是网格有多个边界环（带孔），而简单的串法只处理单环，走到分叉处就断了。

另一种可能是网格非流形：某条边被三个以上三角形共享，或者有孤立顶点。这类网格在数学上根本不是流形，任何参数化方法都会出问题。用meshCheck里的边计数扩展一下，统计每条边被几个面共享，大于 2 的直接列出来，先修网格再谈映射。我一般用 MeshLab 之类的工具做一次非流形清理，比在 MATLAB 里写修复逻辑快得多。

5.3 现象：解线性系统时 MATLAB 警告"矩阵接近奇异或缩放错误"，或者解出来的 uv 数值大到离谱

直接原因是LII不是正定的。cotangent 权重在钝角三角形上是负的，负权重多了以后矩阵失去对角占优性质，Cholesky 分解失败。诊断方法是打印full(LII)的最小特征值，负的或者接近零就中招了。

三个处理方向，按代价从低到高：一是给LII的对角线加一个很小的正则项LII + 1e-8 * speye(size(LII))，代价是解略微偏离精确解，但通常看不出来；二是做一次 Delaunay 翻边，把钝角三角形消掉，这需要额外的拓扑操作，代码量不小；三是换求解器，用minres或者gmres这类迭代法，配合不完全 Cholesky 预条件。我一般先用正则项顶着，等结果质量真的不够了再上翻边。

还有一个容易忽略的点：如果网格里有孤立顶点（不属于任何三角形），它在L里对应的整行整列都是 0，LII必然奇异。meshCheck里加一句find(sum(W,2) == 0)就能提前发现。

5.4 现象：边界附近的三角形畸变特别严重，中央区域反而还好

这不是 bug，是边界处理策略的问题。harmonicMap里用linspace等弧长分配边界角度，对边界边长差异大的网格来说，短边对应的角度和长边一样，映射后短边被拉长、长边被压短，边界附近就炸了。改进是按边界边长累积长度来分配角度参数，让每个边界三角形的形状尽量保持。

% 按边长比例分配边界角度，替代等弧长 seg = vecnorm(V(loop,:) - V(loop([2:end 1]),:), 2, 2); t = [0; cumsum(seg)]; t = t(1:end-1) / sum(seg) * 2*pi; uvB = radius * [cos(t), sin(t)];

这段改动只有三行，效果却很明显。我在一个头部模型上对比过，边界最大角度误差从 0.31 弧度降到 0.08 弧度。如果你的边界是特征边（比如嘴唇轮廓），还应该把这些特征点固定到圆上的特定位置，而不是任由边长决定。

5.5 现象：同一份数据跑两次，结果数值不一样

MATLAB 的稀疏矩阵求解在某些情况下不是严格确定性的，尤其是开了多线程的 BLAS 之后，浮点累加顺序会变，最后一位有差异。正常情况下这点差异在 1e-12 量级，不影响结果。但如果你看到的是肉眼可见的差异，那多半是迭代法（比如 Ricci 流）的初始值或者随机性导致的。

ricciFlow的初值u = zeros(n,1)是确定的，所以问题一般出在早停条件上——如果err在某轮刚好落在tol边界附近，两次运行的迭代次数可能差一轮，结果就有细微差别。要完全可复现，把tol调紧一点，或者干脆固定maxIter不用早停。另外rng的设置也检查一下，如果你的代码里有用到随机采样的部分（比如随机选钉点），记得在脚本开头固定种子。

6. 二次开发：把共形映射接进自己的处理管线

真正用起来的时候，你很少会单独跑这个工具箱，更多是把它当成一个环节塞进已有的管线里。我踩过的坑集中在两处：数据接口和验证方式。

数据接口上，建议在工具箱和你自己的代码之间加一层薄适配，别让V/F的约定渗透到整个工程。比如你自己的数据结构是struct或者pointCloud对象，就写一个toCgeom.m和一个fromCgeom.m，两边互转。我吃过亏：一开始图省事，在业务代码里到处直接传F矩阵，后来想换另一套参数化实现，发现要改三十几个文件。现在无论接什么第三方库，适配层都是第一步。

function m = toCgeom(mesh) %TOCGEOM 你自己的网格结构 -> 工具箱的 V/F 约定 m.V = double(mesh.vertices); % 强制 double m.F = double(mesh.faces); % 索引从 1 开始 assert(min(m.F(:)) >= 1, '内部结构必须用 1-based 索引'); end function mesh = fromCgeom(m, uv) %FROMCGEOM 把 uv 写回你自己的结构，保持行序一致 assert(size(uv,1) == size(m.V,1), 'uv 与顶点数不匹配'); mesh = m; mesh.uv = uv; end

验证方式也得提前定好。我现在的做法是每次改完求解器，先跑一组固定的回归模型，比对三项指标：翻转面数必须为 0，角度误差不能比改动前差，面积标准差在阈值内。这三项过了才算改对。光看图是不够的，图上很小的畸变在纹理上看不出来，等到批量处理几百个模型才会集体爆发。

一个具体的优化技巧：如果你的管线只需处理近似圆盘拓扑的网格，别上 Ricci 流，调和映照加边界按边长分点就够了，速度快一个数量级。Ricci 流的价值在于处理任意拓扑和精确控制曲率，代价是几百轮迭代。选型先看你的拓扑复杂度，再谈精度。

从那以后我每次往管线里接新的几何算法，都强制先跑一遍小样本回归，把指标记在文件里，再上全量数据。这个习惯救过我好几次，尤其是那种"改了个看起来无关紧要的参数，结果整体畸变悄悄变大"的情况。希望帮到你。

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