简介面向计算共形几何课程作业与几何算法实践这份项目框架适合高校学生和开发者搭建可编译运行的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 半边结构不是炫技是为了让边界能被走出来三角网格最朴素的存法是两个数组顶点坐标Vnv×3和面片索引Fnf×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(ii1); 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(end1) half(h,1); h mapStart(half(h,2)); end if numel(ring) 2, rings{end1} 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(k1,3)1); % 对角顶点 u P(mod(k,3)1,:) - P(mod(k1,3)1,:); v P(k,:) - P(mod(k1,3)1,:); cotA dot(u,v) / norm(cross(u,v)); % cot cos/sin p p1; I(p)i; J(p)j; S(p)0.5*cotA; p p1; 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-17sparse在后续\求解时会因为结构不对称走更慢的路径对称化一次后面省很多时间。W - spdiags(diag(W),...)处理的是退化情况如果网格里出现重复顶点组成的针状三角形i j会写进对角不清掉会让sum(W,2)重复计数。如果模型来自扫描或布尔运算钝角三角形会非常多直接用原始余切权重解出来的结果常常带一片翻转。常见做法是两条路一是翻边把网格重连成内在 Delaunayinstrinsic 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, c1) re; A(2*f-1, c2) -im; % 实部方程 A(2*f, c1) im; A(2*f, c2) 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_iu_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(k1,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(k1,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))等于 0cotA变成 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-basedOBJ 也是 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扫参数组合结果落成 tablecfgList 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回头一查是正则加得太大把保角性磨掉了。从那以后我每次调参都强制走一遍这两个指标不看图先看数省下来的返工时间比写代码多得多。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/26 18:42:02 · 阅读 1.2k 三角网格保角参数化实战：余切拉普拉斯、LSCM与Ricci流
简介面向计算共形几何课程作业与几何算法实践这份项目框架适合高校学生和开发者搭建可编译运行的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 半边结构不是炫技是为了让边界能被走出来三角网格最朴素的存法是两个数组顶点坐标Vnv×3和面片索引Fnf×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(ii1); 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(end1) half(h,1); h mapStart(half(h,2)); end if numel(ring) 2, rings{end1} 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(k1,3)1); % 对角顶点 u P(mod(k,3)1,:) - P(mod(k1,3)1,:); v P(k,:) - P(mod(k1,3)1,:); cotA dot(u,v) / norm(cross(u,v)); % cot cos/sin p p1; I(p)i; J(p)j; S(p)0.5*cotA; p p1; 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-17sparse在后续\求解时会因为结构不对称走更慢的路径对称化一次后面省很多时间。W - spdiags(diag(W),...)处理的是退化情况如果网格里出现重复顶点组成的针状三角形i j会写进对角不清掉会让sum(W,2)重复计数。如果模型来自扫描或布尔运算钝角三角形会非常多直接用原始余切权重解出来的结果常常带一片翻转。常见做法是两条路一是翻边把网格重连成内在 Delaunayinstrinsic 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, c1) re; A(2*f-1, c2) -im; % 实部方程 A(2*f, c1) im; A(2*f, c2) 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_iu_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(k1,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(k1,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))等于 0cotA变成 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-basedOBJ 也是 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扫参数组合结果落成 tablecfgList 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回头一查是正则加得太大把保角性磨掉了。从那以后我每次调参都强制走一遍这两个指标不看图先看数省下来的返工时间比写代码多得多。希望帮到你。本文还有配套的精品资源点击获取