简介：面向计算共形几何课程作业与几何算法实践，这份项目框架适合高校学生和开发者搭建可编译运行的C++实验环境，重点覆盖割图、调和映射等算法模块，并配有示例模型、资源文件以及MeshLib、freeglut等第三方依赖，可在Windows、Ubuntu和macOS上使用。压缩包共373个文件、约7.69MB，以295个h头文件、17个m文件、4个cpp源文件为主，另含png/bmp/jpg预览图、txt/md说明文档、CMake配置及库文件，目录按算法、数据与资源分置，便于按模块检索。已有234人学习下载。读者可据此复现共形几何基础流程，理解工程组织、跨平台编译配置与观察器、Hodge分解等源码入口，适合作为课程作业和算法实现的参考框架。
1. 共形几何不是玄学：一张展开就撕烂的三角网格说明了什么
把一张人脸扫描网格摊平到平面上，十有八九会在鼻翼两侧翻出几个三角形，或者整块贴图被拉成一条细带。这不是代码写错了，而是"共形"两个字本身就带着硬约束：保角就意味着面积必然被牺牲，曲率集中处必然出问题。Conformal-Geometry 这套 MATLAB 资源干的事情，就是把曲面到平面、曲面到球面的保角映射拆成能算的稀疏线性系统——余切拉普拉斯、最小二乘共形映射、离散 Ricci 流，一条一条摆出来给你改。它适合已经会用 MATLAB、但被网格参数化卡住的人：做贴图展开、做曲面展平、做几何处理方向的课程设计，都能直接换掉输入数据跑自己的模型。下面从数据结构讲到调参，中间翻过的车一并写出来。
2. 半边结构与余切权重：为什么你的第一版展开总会翻几个面
2.1 半边结构不是炫技，是为了让边界能被"走"出来
三角网格最朴素的存法是两个数组：顶点坐标V（nv×3）和面片索引F（nf×3）。这个存法足够画图，但不够做参数化。共形映射要解的是顶点上的方程，而方程的边界条件来自"边界环"，边界环要靠"这条边的对面是谁"一路串下去。没有邻接关系，你只能每次用ismember现搜，几千个面还能忍，几十万个面直接卡死。
半边结构把每条有向边当成一个对象，记录起点、终点、孪生半边（反向的那条）、所属面、同面内的下一条半边。有了孪生指针，"这条边是不是边界"就退化成一次判空，边界环就是顺着"边界半边的终点出发的下一条边界半边"一路走。下面这段是我常用的构建方式，用排序代替哈希表，避免containers.Map在大网格上的性能塌方：
function [half, rings] = build_halfedge(V, F) % 输入 % V : nv x 3 顶点坐标 % F : nf x 3 面片索引，1-based，逆时针（法向朝外） % 输出 % half : (3*nf) x 5，列 = [起点, 终点, 孪生, 面号, 同面下一条] % rings : cell，每个元素是一条有序边界顶点环 nf = size(F,1); nv = size(V,1); nh = 3*nf; half = zeros(nh, 5); for f = 1:nf for k = 1:3 h = (f-1)*3 + k; k2 = mod(k,3) + 1; half(h,1) = F(f,k); % 起点 half(h,2) = F(f,k2); % 终点 half(h,4) = f; % 所属面 half(h,5) = (f-1)*3 + k2; % 同面下一条半边 end end % 用 (min,max) 顶点对排序找孪生，比哈希稳，内存也可控 key = sort([half(:,1), half(:,2)], 2); [~, ord] = sortrows(key); for ii = 1:2:nh-1 h1 = ord(ii); h2 = ord(ii+1); if isequal(key(h1,:), key(h2,:)) half(h1,3) = h2; half(h2,3) = h1; end end % 边界环提取：构造"起点 -> 边界半边"的映射 isBnd = (half(:,3) == 0); mapStart = zeros(nv,1); for h = find(isBnd)' mapStart(half(h,1)) = h; % 一个顶点若两条出边会覆盖，见 5.3 end rings = {}; used = false(nh,1); for h0 = find(isBnd)' if used(h0), continue; end ring = []; h = h0; while h > 0 && ~used(h) used(h) = true; ring(end+1) = half(h,1); h = mapStart(half(h,2)); end if numel(ring) > 2, rings{end+1} = ring; end end end这段代码里最要紧的是F的绕序。余切权重、梯度、面积全部依赖法向方向，一旦面片是顺时针，所有计算出来的角度符号都反过来，结果是"能出图，但整张图是镜像的"，而且你不会立刻发现——因为拓扑看起来完全正常。我一般会在读入之后立刻加一道自检：把所有三角面的有向面积加起来，sum(dot(cross(V(F(:,2),:)-V(F(:,1),:), V(F(:,3),:)-V(F(:,1),:)), repmat([0 0 1],nf,1), 2))，如果是负的，就执行F(:,[2 3]) = F(:,[3 2])统一翻面，别拖到后面再查。
另一个容易忽略的点是rings返回的顺序。harmonic_disk这类函数要求传入的边界是沿环有序的，顺序错了出来的不是圆盘而是自交的多边形，看上去像"参数化爆炸"。rings里每条环的顺序由半边方向决定，和面片绕序一致，所以只要F的绕序是对的，这里的顺序就天然可用。
2.2 余切权重：三个角度决定了整个矩阵的脾气
离散共形映射的核心是离散拉普拉斯算子。在三角网格上，它有一个非常干净的来源：分片线性有限元对 Dirichlet 能量的离散。对一条边(i,j)，它的权重等于两个相邻三角形在对面顶点处余切之和的一半：
$$W_{ij} = \frac{1}{2}\left(\cot\alpha_{ij} + \cot\beta_{ij}\right)$$
这个权重的好处是"只依赖角度，不依赖尺度"——把整个模型放大十倍，权重完全不变，这正是保角性质在离散层面的体现。坏处也很直白：余切在钝角处为负数。一个三角形只要有一个角超过 90 度，对边权重就可能变成负的，拉普拉斯矩阵的非对角元出现正值，半正定性被破坏。这时候线性系统照样能解，但解出来的 uv 坐标可能直接翻面。
function [L, W] = cot_laplacian(V, F) % L : 余切拉普拉斯，L(i,i)=sum(W(i,:))，L(i,j)=-W(i,j) % W : 对称权重矩阵，W(i,j) = (cot α + cot β)/2 nv = size(V,1); nf = size(F,1); m = 6*nf; % 每个面 3 条边，每条边 2 个非零 I = zeros(m,1); J = I; S = I; p = 0; for f = 1:nf idx = F(f,:); P = V(idx,:); for k = 1:3 i = idx(k); % 边起点 j = idx(mod(k,3)+1); % 边终点 o = idx(mod(k+1,3)+1); % 对角顶点 u = P(mod(k,3)+1,:) - P(mod(k+1,3)+1,:); v = P(k,:) - P(mod(k+1,3)+1,:); cotA = dot(u,v) / norm(cross(u,v)); % cot = cos/sin p = p+1; I(p)=i; J(p)=j; S(p)=0.5*cotA; p = p+1; I(p)=j; J(p)=i; S(p)=0.5*cotA; end end W = sparse(I(1:p), J(1:p), S(1:p), nv, nv); W = (W + W')/2; % 抹掉浮点不对称 W = W - spdiags(diag(W),0,nv,nv); % 清掉可能出现的自环 L = spdiags(sum(W,2), 0, nv, nv) - W; end几个参数和写法上的说明。norm(cross(u,v))等于三角形面积的两倍，把它放在分母上是正确的，因为余切本来就是cos/sin，而sin正比于面积。W = (W + W')/2不是为了美观：浮点误差会让W(i,j)和W(j,i)差一个 1e-17，sparse在后续\求解时会因为结构不对称走更慢的路径，对称化一次，后面省很多时间。W - spdiags(diag(W),...)处理的是退化情况：如果网格里出现重复顶点组成的"针状"三角形，i == j会写进对角，不清掉会让sum(W,2)重复计数。
如果模型来自扫描或布尔运算，钝角三角形会非常多，直接用原始余切权重解出来的结果常常带一片翻转。常见做法是两条路：一是翻边，把网格重连成内在 Delaunay（instrinsic Delaunay），保证所有余切非负；二是不动拓扑，在解线性系统时对负权重做截断，或者给对角加一个小正则项。前者几何更干净但要写翻边逻辑，后者一行代码但会引入轻微的角度偏差。我一般先跑截断版本确认流程通，再决定要不要上翻边。
2.3 三条路线的选型：先看拓扑，再看你要保什么
同一份L，能撑起三种完全不同的实现路线。选错了不是跑不动，是跑出来的东西不对路。
| 路线 | 适用拓扑 | 边界条件 | 保角程度 | 计算量 | 典型失败模式 |
|---|---|---|---|---|---|
| 调和映射 | 圆盘（单连通带边） | 边界必须预设 | 内部共形，边界拉伸大 | 一次稀疏求解 | 边界不圆润时内部被挤扁 |
| LSCM | 圆盘，可带多个固定点 | 固定 2 个点消退化 | 全局最小二乘意义下最优 | 一次超定最小二乘 | 大面积区域被压成窄条 |
| 离散 Ricci 流 | 圆盘 / 球面 / 环面 | 目标曲率 | 严格离散共形（保圆 Packing 度量） | 迭代数十到数百次 | 步长过大震荡不收敛 |
调和映射最省事但边界是"死"的：你给的边界形状直接决定了内部被拉成什么样，边界上弧长分布不均，内部就会出现一条被压扁的带子。LSCM 把"共形"写成最小二乘，允许整体自由变形，只在两个点上钉住，所以对边界形状不敏感，代价是它在最小二乘意义下最优，局部仍然可能有大畸变。Ricci 流走的是第三条路：不动坐标，动度量——迭代每个顶点的共形因子，让离散高斯曲率逼近目标值，收敛后再把度量铺回平面。它最慢，但它是三者里唯一"结果由曲率唯一决定"的，做共形模、做叶状结构必须走这条。
选型的判断顺序我一般是：先看亏格和边界数，圆盘且只要一张贴图就上调和或 LSCM；要球面参数化、要环面上的模数、要做保角等价类的比较，直接上 Ricci 流，别在另外两条路上绕。
3. 三条实现路线落地：调和映射、LSCM 与 Ricci 流的可运行代码
3.1 调和映射：把边界摁在圆上，内部解 Laplace
调和映射的求解思路非常直接：给边界顶点指定平面坐标，内部顶点的 uv 满足离散 Laplace 方程L_ii * uv_i + L_ib * uv_b = 0，移项就是一次线性求解。关键在于边界怎么摆。均匀角度分布在边界顶点数不均时会出问题——一条边上挤了 200 个点，另一条边只有 20 个，均匀分角度会把前者压成一团。按弧长重参数化能把这个畸变压下去一大截。
function uv = harmonic_disk(V, F, bnd) % bnd : 有序边界顶点索引（沿环顺序，与面片绕序一致） nv = size(V,1); [L, ~] = cot_laplacian(V, F); % 边界按累计弧长参数化到单位圆上 P = V(bnd,:); seg = sqrt(sum((P([2:end 1],:) - P).^2, 2)); s = [0; cumsum(seg(1:end-1))] / sum(seg); uv = zeros(nv,2); uv(bnd,1) = cos(2*pi*s); uv(bnd,2) = sin(2*pi*s); inner = setdiff((1:nv)', bnd(:)); A = L(inner, inner); B = -L(inner, bnd) * uv(bnd,:); % 边界项搬到右端 uv(inner,:) = A \ B; % 稀疏 Cholesky endA是内部顶点之间的子矩阵，规模等于内部点数，稀疏且对角占优（在没有负权重的前提下），MATLAB 的\会自动选 Cholesky，几万个点通常在一秒内。B的构造顺序不能反：L(inner,bnd)是"内部行 × 边界列"，乘上边界坐标才得到每个内部顶点受到的边界牵引。
这段代码有一个隐含前提：输入必须是圆盘拓扑。如果模型带多个洞（比如一个带两个孔的板），bnd会有多条环，直接塞进uv(bnd,:)会得到一个自交的边界，解出来的内部网格大面积翻转。正确做法是先沿最短路径把洞切开，把多连通区域变成单连通，切开的缝在两侧各复制一份顶点。这套切开逻辑不复杂但很长，实践中我一般先用genus和numel(rings)判断，多洞的直接走 Ricci 流，别在调和映射上硬撑。
另一个参数是圆的半径。半径只影响尺度不影响形状，因为 Laplace 方程是线性的，边界等比缩放，内部解也等比缩放。所以定成 1 就行，不要在这上面调参。
3.2 LSCM：把共形条件写成超定方程组
LSCM 的出发点和调和映射不同：它不预设边界，而是要求每个三角形的局部展开尽量接近"相似变换"。做法是把每个三角形先旋转到它自己的局部平面，得到三个二维点，然后用这三点构造共形条件的系数矩阵，把所有三角形堆起来形成一个超定系统A * x = 0，再固定两个点消掉平移旋转自由度，解最小二乘。
function uv = lscm(V, F, pin) % pin : [p1 p2] 两个固定顶点，用来消掉整体的平移旋转 nv = size(V,1); nf = size(F,1); A = sparse(2*nf, 2*nv); for f = 1:nf idx = F(f,:); P = V(idx,:); e1 = P(2,:) - P(1,:); e1 = e1 / norm(e1); nrm = cross(P(2,:)-P(1,:), P(3,:)-P(1,:)); area2 = norm(nrm); % 2 * 三角形面积 e2 = cross(nrm/area2, e1); Q = P - P(1,:); x = Q * e1'; % 局部 x 坐标 y = Q * e2'; % 局部 y 坐标 % 共形条件的复系数：w_k = dx + i*dy，按对边差分构造 wk = [ (x(3)-x(2)) + 1i*(y(3)-y(2)), ... (x(1)-x(3)) + 1i*(y(1)-y(3)), ... (x(2)-x(1)) + 1i*(y(2)-y(1)) ] / sqrt(area2/2); for k = 1:3 re = real(wk(k)); im = imag(wk(k)); c = 2*(idx(k)-1); A(2*f-1, c+1) = re; A(2*f-1, c+2) = -im; % 实部方程 A(2*f, c+1) = im; A(2*f, c+2) = re; % 虚部方程 end end % 固定两个点，剩余自由度解最小二乘 free = setdiff((1:nv)', pin(:)); colsFree = [2*free-1, 2*free]'; colsFree = colsFree(:); colsPin = [2*pin(:)-1, 2*pin(:)]'; colsPin = colsPin(:); d = norm(V(pin(2),:) - V(pin(1),:)); % 保尺度：固定两点距离 uvPin = [0 0; d 0]; b = -A(:, colsPin) * uvPin(:); x = zeros(2*nv,1); x(colsFree) = A(:, colsFree) \ b; % 超定，MATLAB 走 QR x(colsPin) = uvPin(:); uv = reshape(x, 2, nv)'; endsqrt(area2/2)这个缩放因子来自 LSCM 的标准推导，作用是把每个三角形的能量按面积加权，面积大的三角形话语权大，避免小三角形主导解。pin选哪两个点有讲究：选得离得太近，uvPin的d很小，整个解的尺度被压得看不出来；选在模型两端比较好。固定点的方向决定了最终结果是正着放还是歪着放，这不影响畸变，但如果要和别的结果做逐点对比，最好固定成同一个方向。
LSCM 的一个常见误解是它"没有翻转"。它在最小二乘意义下最优，但最优解仍然可以在局部翻转，尤其是细长的三角形带。判翻转的办法很简单：算每个面在 uv 平面上的有向面积，和三维面上相对于某个方向的有向面积比，符号为负的就是翻面。翻转面积占比超过 1%，基本就该考虑换路线或者先做一次重网格了。
3.3 离散 Ricci 流：不动坐标，动度量
Ricci 流的逻辑更接近几何本身：它不是去找一个映射，而是去改每个顶点处的尺度，让曲率变成想要的值。对每个顶点i有一个共形因子u_i，边长按l_ij ← l_ij * exp((u_i+u_j)/2)更新。每轮由边长用余弦定理算出三角形内角，角亏K_i = 2π - Σθ_i（内部点）或π - Σθ_i（边界点），然后u_i ← u_i + step * (K_target - K_i)，直到角亏收敛。
function [uv, u] = ricci_flow_disk(V, F, bnd, iters, step) % 圆盘拓扑，目标曲率：内部 0，边界按总曲率 2π 均分 nv = size(V,1); nf = size(F,1); E = [F(:,[1 2]); F(:,[2 3]); F(:,[3 1])]; E = unique(sort(E,2), 'rows'); % 唯一边表 ne = size(E,1); % 边 -> 列号索引 eIdx = containers.Map('KeyType','double','ValueType','double'); for e = 1:ne eIdx(E(e,1)*1e7 + E(e,2)) = e; end l = sqrt(sum((V(E(:,1),:) - V(E(:,2),:)).^2, 2)); % 初始欧氏边长 u = zeros(nv,1); isB = false(nv,1); isB(bnd) = true; Ktar = zeros(nv,1); Ktar(~isB) = 0; % 内部目标曲率为 0（平坦） for it = 1:iters % 由当前边长算每个顶点角亏 ang = zeros(nv,1); for f = 1:nf idx = F(f,:); for k = 1:3 a = idx(k); b = idx(mod(k,3)+1); c = idx(mod(k+1,3)+1); lab = l(eIdx(min(a,b)*1e7 + max(a,b))); lac = l(eIdx(min(a,c)*1e7 + max(a,c))); lbc = l(eIdx(min(b,c)*1e7 + max(b,c))); % 余弦定理求顶点 a 处的角 ang(a) = ang(a) + acos((lab^2 + lac^2 - lbc^2) / (2*lab*lac)); end end K = 2*pi - ang; K(isB) = pi - ang(isB); % 边界点用 π 而不是 2π err = norm(K - Ktar, inf); if err < 1e-6, break; end u = u + step * (Ktar - K); % 边长按共形因子更新，并限制单步幅度防震荡 u = max(min(u, 5), -5); l = l .* exp(0.5*(u(E(:,1)) + u(E(:,2)))); end uv = layout_from_lengths(E, l, bnd); % 由边长铺回平面，见下 endstep的建议取值范围是 0.1 到 0.5。取 1.0 在大部分网格上会直接震荡不收敛，因为角亏对u的响应不是线性的，一次走满会过冲。u的上下界截断是后悔药：某几个顶点的曲率在某一轮突然爆掉，exp会把边长拉到 1e10，后面再也拉不回来，一截断至少还能继续迭代。
layout_from_lengths是收尾的一步：从一个三角形开始，用 BFS 逐个把相邻顶点按已知边长放到平面上（两个圆求交），遇到已经放好的顶点就用最小二乘微调位置。这一步的误差主要来自三角不等式不严格成立——Ricci 流收敛后，绝大多数三角形满足，但个别退化三角形仍然会给出无解的圆交。我在实现里通常加一个容差：如果两圆相交超过1 + 1e-8，就截断成 1。这个容差设得太大，铺出来的网格会在这些区域明显起皱；设得太小，BFS 会卡死。
4. 结果验真：角度畸变、纹理回贴与跨工具比对
4.1 角度畸变不是看一眼图就行，要出数字
肉眼看 uv 图只能判断有没有大面积翻转，判断不了"保角质量"。真正的指标是逐三角形的角度畸变。做法是：对每个面，算出三维三角形三个内角，再算出 uv 三角形三个内角，取对应角的比值，最大比值就是这块面的拟共形畸变。
function qc = quasi_conformal_error(V, F, uv) % qc : nf x 1，每个面的拟共形畸变，1 表示完美保角 nf = size(F,1); qc = zeros(nf,1); for f = 1:nf idx = F(f,:); t3 = tri_angles(V(idx,:)); % 三维三个角 t2 = tri_angles([uv(idx,:), zeros(3,1)]); % uv 三个角 r = t3 ./ max(t2, 1e-12); qc(f) = max(max(r), max(1./r)); end end function t = tri_angles(P) % P : 3 x 3，返回三个内角（弧度），顺序与顶点顺序一致 t = zeros(3,1); for k = 1:3 a = P(mod(k,3)+1,:) - P(k,:); b = P(mod(k+1,3)+1) - P(k,:); t(k) = acos(dot(a,b)/(norm(a)*norm(b))); end end统计的时候别看均值，看 95 分位和最大值。均值会被大片平坦区域拉低，真正出问题的是那 0.1% 的面。我的经验阈值是：95 分位低于 1.15 算能用，最大值超过 3 的区域基本会出视觉瑕疵，要单独标出来回到模型上定位——通常是那儿有一片钝角三角形，或者边界上有一处剧烈转折。
4.2 纹理坐标回贴：导出这一步最容易前功尽弃
uv 算对了，导不出去等于白干。OBJ 里纹理坐标靠vt行和f v/vt索引对齐，格式错一个位置，进三维软件就是一团黑或者贴图乱飞。
function write_obj_uv(path, V, F, uv) fid = fopen(path, 'w'); if fid < 0, error('open failed: %s', path); end fprintf(fid, 'v %.6f %.6f %.6f\n', V'); % 按列循环，每行一个顶点 fprintf(fid, 'vt %.6f %.6f\n', uv'); % 每行一个纹理坐标 fprintf(fid, 'f %d/%d %d/%d %d/%d\n', [F, F]'); fclose(fid); endfprintf按列消费参数，所以V'是 3×nv，正好每列三个数一行。[F, F]'是 6×nf，每列六个数字，对应f a/a b/b c/c的格式。vt的索引和v索引在这里恰好相同，所以直接复用F，这也是为什么 uv 数组必须和 V 顶点一一对应、不能只存边界。
坐标范围也要注意。uv 落在 [0,1] 之外在 OBJ 里是合法的，纹理默认重复采样，但如果材质里设了 clamp，边界那一圈会被拉成条纹。稳妥做法是导出前统一归一化到 [0,1]，同时保留一个缩放系数，回溯的时候乘回去。
4.3 和另一套实现做交叉验证
只跑一套实现，出问题的时候你不知道是算法理解错了还是数据有问题。我习惯用两个指标做交叉核对：一是总角度畸变的分布形状，如果两套实现差了一个数量级，多半是面积加权项写错了；二是边界环的定向，把 uv 边界画出来，看它是不是和输入网格的边界环一一对应。
| 核对项 | 计算方式 | 合理范围 | 异常指向 |
|---|---|---|---|
| 95 分位角度畸变 | 4.1 节的 qc | < 1.15 | 权重或面积因子错 |
| 翻转面占比 | uv 面有向面积符号为负的比例 | < 0.1% | 钝角三角形过多 |
| uv 包围盒长宽比 | max-min 的 x/y 比 | 与模型本身接近 | 边界参数化塌缩 |
| 边界自交数 | 边界折线段两两求交 | 0 | 边界顺序或切缝错 |
这套核对不需要额外的参考实现，只靠输入模型的几何统计就能做。它的价值在于：当四个指标里只有一个异常，你大概能立刻定位到哪一步出了问题，而不是对着一张 uv 图反复猜。
5. 避坑与排查：五次翻车现场的现象、原因与解法
5.1 展开图看着挺好，贴图上文字左右反了
现象：uv 图整体形状正常，没有翻转面，但把带文字的贴图贴上去，字是镜像的。
原因：面片绕序是顺时针。所有基于有向面积的量（法向、梯度、有向面积符号）全都反了，导致 uv 平面的定向和三维曲面的定向相反，几何上等价于做了一次镜面反射。因为保角映射对镜像不敏感，畸变指标完全正常，所以很难在数值上发现。
解决：读入网格后立刻算全局有向面积和，为负就F(:,[2 3]) = F(:,[3 2])。同时在harmonic_disk里保证bnd的顺序和F的绕序一致，两者不能各用各的。
5.2 解出来一片 NaN，第一行就崩
现象：uv(inner,:) = A \ B之后，uv 里出现大量 NaN 或者 1e15 量级的值。
原因：网格里存在退化三角形——三个顶点共线，或者两个顶点坐标完全相同。这时norm(cross(u,v))等于 0，cotA变成 Inf，写进稀疏矩阵后整个A被污染。扫描数据和 CAD 导出的模型里这种面非常常见。
解决：求解前做一次清理，把面积小于1e-12 * bboxDiag^2的面删掉，把距离小于1e-9 * bboxDiag的顶点合并。清理后重新走一遍边界环提取，因为删面会改变边界。这一步不要省，我吃过一次亏——一个 4 万面的模型里只有一个退化面，查了两小时。
5.3 边界环只走出一小段就断了
现象：rings里返回的环只有几十个顶点，但模型边界明显有上千个顶点。
原因：mapStart是一个nv×1的数组，如果一个顶点有两条出边（非流形顶点，两个面片只在一点相接），后写入的会覆盖前面的，走环的时候从那一点直接跳走，剩下的边界永远走不到。
解决：把mapStart改成 cell 数组，每个顶点存一个出边列表，走环时优先选未被使用的。更进一步，先做一次非流形检测：统计每个顶点的邻接面数，如果某个顶点的邻接面构成的图不连通，就是非流形顶点，直接切开或者报错让用户处理。
5.4 矩阵解出来了，但角畸变特别大，还伴随局部翻转
现象：没有 NaN，解也算完了，但 95 分位畸变到了 2 以上，翻转面占比 5%。
原因：钝角三角形太多，余切权重出现负值，L不再半正定，解出来的调和函数在负权重区域出现"振荡"，表现为局部翻面。这个现象在低质量网格上极其普遍，但在指标上往往表现为"整体还行，局部爆炸"。
解决：三条路按成本排序。最省事的是在解之前把W里的负值截断成 0，代价是引入轻微的角度偏差；中等成本是给对角加1e-8 * trace(L)/nv的正则项，能压住振荡但不改拓扑；最彻底的是做 intrinsic Delaunay 翻边，把所有余切变正，这一步要写翻边逻辑，但做完之后所有下游算法都受益。
5.5 导出的 OBJ 在三维软件里是一团黑
现象：文件能打开，模型也在，但贴图完全不显示，或者显示成纯色。
原因：三种可能。一是vt行写在v行之后但f里没带/vt索引；二是f行的索引越界（MATLAB 是 1-based，OBJ 也是 1-based，这一步通常没问题，但如果做过顶点去重而忘了重映射F，就会越界）；三是没有法线，某些软件对无法线的模型默认关掉光照。
解决：导出后写一个自检函数，读回文件，检查max(F(:)) <= nv、numel(uv)/2 == nv、每个f行的三段索引一致，再补一段vn行（用面法向平均到顶点）。这个自检不到二十行，但能省掉大量来回。
6. 进阶技巧：用共形模和批处理脚本把结果钉死
单张模型跑通只能说明流程没崩，要判断这套实现靠不靠谱，得看它在"有已知答案"的问题上表现如何。共形模就是个很好的标尺。一个平环（annulus）在共形意义下由内外半径之比唯一确定，这个比值是共形不变量，跟你怎么参数化无关。用离散 Ricci 流把环面参数化到平面，然后量出内外边界的周长比，取对数就是模数。同一个模型换不同的网格密度、不同的初始度量，算出来的模数应该收敛到同一个值——如果差得远，说明要么曲率目标设错了，要么layout_from_lengths的铺开误差太大。
function m = annulus_modulus(V, F, bndOuter, bndInner) % 用 Ricci 流把环面铺平，返回共形模 log(R/r) [uv, ~] = ricci_flow_annulus(V, F, bndOuter, bndInner, 300, 0.3); % 分别取内外边界的 uv 坐标，算平均半径 pO = uv(bndOuter,:); pI = uv(bndInner,:); rO = mean(sqrt(sum((pO - mean(pO,1)).^2, 2))); rI = mean(sqrt(sum((pI - mean(pI,1)).^2, 2))); m = log(rO / rI); end参数上，iters给 300 是因为环面的曲率分布比圆盘更难收敛，边界曲率要按内外分别分配，总曲率得对得上欧拉示性数。step给 0.3 是折中：再大就震荡，再小则 300 轮不够。判断收敛不能只看err，还要看模数在最后 50 轮里的变化幅度，稳定在 1e-3 以内才算数。
另一类进阶用法是把整套流程包成一个可扫描的批处理。我通常写一个 config 结构体，把网格路径、路线选择、权重模式、正则项、目标曲率全塞进去，外层用parfor扫参数组合，结果落成 table：
cfgList = struct2table( ... combvec(1:3, [0 1e-8 1e-6], [0.1 0.3 0.5])'); % 路线 × 正则 × 步长 cfgList.Properties.VariableNames = {'route','reg','step'}; res = cell(height(cfgList), 4); parfor i = 1:height(cfgList) cfg = cfgList(i,:); [uv, info] = run_pipeline(meshPath, cfg); res(i,:) = {cfg, info.qc95, info.flipRatio, info.iters}; endrun_pipeline里统一做网格清理、算子构建、按route分支调用三条路线之一，最后返回畸变指标和迭代次数。这样一张表跑下来，你能直接看出"哪种网格该走哪条路线"，而不是靠感觉猜。注意parfor里不要共享稀疏矩阵，每个 worker 各自建各自的，内存换时间。
最后一个习惯：每次改动权重或者边界参数之后，别只看最终的 uv 图，先跑一遍翻转面占比和 95 分位畸变这两个数，和上一版做差值。有几次我改了一行正则项，图上看不出差别，但畸变分位数涨了 0.3，回头一查是正则加得太大，把保角性磨掉了。从那以后我每次调参都强制走一遍这两个指标，不看图先看数，省下来的返工时间比写代码多得多。希望帮到你。
本文还有配套的精品资源，点击获取
