簡介基于C實現(xiàn)的增量Delaunay三角剖分算法工程包面向計算幾何學(xué)習(xí)者與需要處理三角網(wǎng)格的C開發(fā)者。算法逐點插入并局部翻邊使剖分中所有三角形的最小角盡量最大化項目采用石溪大學(xué)顧賢峰教授提出的“半邊”數(shù)據(jù)結(jié)構(gòu)組織網(wǎng)格適合想深入理解幾何拓撲實現(xiàn)的人群。壓縮包共89個文件、4.17MB包含13個h頭文件、1個cpp主實現(xiàn)位于GaussCurvature_Sample/GaussCurvature/main.cpp、Visual Studio解決方案與工程文件、已編譯的exe及構(gòu)建日志可直接打開工程運行或按需拆解調(diào)用。已有1003人瀏覽/學(xué)習(xí)。讀者可拿到完整可編譯的Delaunay三角剖分示例厘清增量插入、局部邊翻轉(zhuǎn)與半邊結(jié)構(gòu)維護的配合關(guān)系還可參考其網(wǎng)格類封裝方式遷移到其他幾何算法項目中。1. Delaunay三角剖分為什么“最不尖銳”的網(wǎng)格要靠局部操作來維護給一堆散點Delaunay三角剖分能生成一張“盡可能勻稱”的三角形網(wǎng)格任何三角形的外接圓內(nèi)都不允許出現(xiàn)第四個點。這個空圓準則聽起來是全局約束——每一圈網(wǎng)格都要滿足它。但真正讓增量算法工程可用的是一條反直覺的性質(zhì)往一張幾十萬個三角形的網(wǎng)格里插入一個新點需要修改的只是和它相鄰的一小片區(qū)域。把其中不合法的那條邊翻一下整個網(wǎng)格的空圓條件便自動恢復(fù)。于是插入一個點的過程可以被拆成定位、拆分、翻邊三步每一步都只圍繞新點的局部鄰域做文章。搞懂這三步背后的幾何機制和C實現(xiàn)取舍才算真正明白這個算法為什么能work又為什么在某些輸入上會翻車。2. 增量插入機制一點擾動為什么不會波及整張網(wǎng)格2.1 局部 Delaunay 邊與空圓條件的等價關(guān)系增量算法的正確性基石是“局部 Delaunay 邊”這個概念。取一條共享邊 AB兩側(cè)分別是三角形 ABC 和 ABD。如果 D 落在三角形 ABC 的外接圓內(nèi)部邊 AB 就是非法的反之如果四邊形 ACBD 是凸的且 D 在 ABC 外接圓外部AB 就是合法的。這個判定與全局空圓條件之間有一個重要等價當網(wǎng)格中所有內(nèi)部邊都滿足局部 Delaunay 條件時整張網(wǎng)格就是全局 Delaunay 的。這意味著修復(fù)可以完全局部化。插入一個點之后受影響的非法邊一定與包含新點的三角形相鄰不存在隔空被破壞的邊。因為空圓條件和局部判定等價所以只要逐條處理非法邊就能恢復(fù)全局性質(zhì)不需要掃描全圖。增量算法里那句口頭禪可以這么記誰碰了新點就查誰的外接圓。這個結(jié)論也解釋了算法偽代碼為什么那么短。初始化超級三角形然后逐個插入點每插一個點只做局部修復(fù)最后刪掉超級三角形——主循環(huán)就這么幾行。真正的工作量全部藏在“如何快速定位”“如何判斷一條邊是否合法”和“如何翻邊而不弄錯鄰居關(guān)系”里。2.2 步行法定位點與三角形的包含關(guān)系插入新點后的第一件事是找到它落在哪個現(xiàn)有三角形內(nèi)部。最樸素的做法是遍歷全部三角形做包含測試每次 O(n)在幾十萬點的規(guī)模下完全不可接受。工程實現(xiàn)里更常用的是步行法walking algorithm從某個已知三角形出發(fā)逐步向目標點移動。每一步檢查目標點 p 相對當前三角形三條邊的位置。三角形按逆時針存儲頂點如果 p 在某條有向邊的右側(cè)就說明應(yīng)該從這條邊跨到對面的鄰居三角形去。直到 p 位于三條邊的同一側(cè)也就是落在當前三角形內(nèi)部定位結(jié)束。int locatePoint(const vectorPoint pts, const vectorTri tris, const Point p, int start) { int t start; for (int guard 0; guard (int)tris.size() * 4; guard) { bool inside true; for (int i 0; i 3; i) { int a tris[t].v[i]; int b tris[t].v[(i 1) % 3]; double o orient2d(p, pts[a], pts[b]); if (o 0) { // p 在邊 (a,b) 的右側(cè) int next tris[t].n[i]; if (next 0) return t; // 理論上不會走到邊界 t next; inside false; break; } } if (inside) return t; } return t; // 保護性退出正常不會執(zhí)行到這里 }這里的orient2d(p, a, b)是計算幾何里的方向測試返回正值表示 p 在從 a 到 b 的有向邊左側(cè)負值表示右側(cè)。三角形頂點按逆時針存儲時如果 p 在三條邊的左側(cè)就說明它在三角形內(nèi)部。每跨一步離目標三角形都更近一點只要網(wǎng)格拓撲正確這個過程必然收斂。代碼里的 guard 循環(huán)上限是必須的。步行法理論上最壞走 O(n) 步但浮點誤差可能讓判斷在兩條邊的臨界處來回振蕩不加保護會死循環(huán)。生產(chǎn)實現(xiàn)里還有一種常見優(yōu)化把 forbreak 展開成三次順序判斷省掉分支跳轉(zhuǎn)。在幾十萬點的規(guī)模下這個細節(jié)能帶來可感知的定位速度差異。2.3 三角形拆分與級聯(lián)翻邊定位完成后新點 p 落在三角形 ABC 內(nèi)。先把 ABC 從網(wǎng)格中摘除再用 p 與 A、B、C 生成三個新三角形 ABp、BCp、CAp。幾何上這一步等于把舊三角形劈成三塊開銷很小。但劈完之后新三角形可能與它的鄰居構(gòu)成共邊四邊形比如三角形 ABp 與鄰居 ABD 共享邊 AB如果 D 落在三角形 ABp 的外接圓內(nèi)AB 就不再是合法邊。這時需要翻邊f(xié)lip把共享邊 AB 換成 Dp兩個三角形變成 ADp 和 DBp。翻邊之后新的 ADp 邊或 DBp 邊可能又讓更外圈的三角形不合法于是翻邊必須遞歸進行。這個遞歸過程像漣漪一樣向外擴散但因為非法邊只能沿著含新點的區(qū)域傳播擴散范圍被限制在鄰域內(nèi)。翻邊必然會終止原因在于它等價于“選擇正確的對角線”。對一個凸四邊形兩條對角線中恰好有一條滿足局部 Delaunay 條件。每翻一次邊四邊形內(nèi)的最小角嚴格變大而網(wǎng)格三角形總數(shù)不變所以遞歸深度存在上界。隨機順序插入時這個深度的期望值接近常數(shù)這也正是增量算法平均效率的關(guān)鍵所在。2.4 一個完整插入序列的動線把定位、拆分、翻邊串起來一次插入完整流程可以概括為int t locatePoint(pts, tris, newPt, lastTriangle); int a tris[t].v[0], b tris[t].v[1], c tris[t].v[2]; // 創(chuàng)建三個新三角形手工維護鄰居關(guān)系 // 刪除舊三角形 t // 對三個新三角形分別執(zhí)行 legalize legalizeTriangle(ta, 1); legalizeTriangle(tb, 1); legalizeTriangle(tc, 1);每個legalizeTriangle處理一塊共享邊如果是內(nèi)部邊且不滿足局部 Delaunay 條件就翻邊然后對新形成的兩條邊繼續(xù)遞歸。每次插入的“動線”都完全一致先定位到包含新點的三角形把它拆開再把新點帶來的非法邊一條條修好。整個增量算法的骨架到這里就清楚了剩下的難點全在數(shù)據(jù)結(jié)構(gòu)——鄰居關(guān)系怎么存、翻邊時如何精準地改寫四個三角形的連接信息。這部分放到下一章展開。3. C 實現(xiàn)數(shù)據(jù)結(jié)構(gòu)、幾何謂詞與翻邊的源碼級拆解3.1 三角形的數(shù)據(jù)表示頂點數(shù)組加鄰居數(shù)組增量算法里最常用的三角形結(jié)構(gòu)是頂點數(shù)組加鄰居數(shù)組struct Tri { int v[3]; // 頂點索引按逆時針排列 int n[3]; // 鄰居索引n[i] 是邊 (v[i], v[(i1)%3]) 對面的三角形 };例如n[0]表示與邊(v[0], v[1])共享的鄰居三角形。如果某條邊沒有鄰居說明它是當前網(wǎng)格的邊界邊用 -1 表示。為什么用 int 索引而不是指針第一個原因是 vector 擴容時元素地址會移動裸指針全部失效索引不會。第二個原因是內(nèi)存局部性一個 Tri 占 24 字節(jié)存放在連續(xù)數(shù)組里幾十萬個三角形遍歷時緩存命中率遠高于用 new 逐個分配的散落對象。第三個原因是翻邊時要頻繁改寫鄰居指向索引賦值就是一次 int 拷貝比指針更簡單也更安全。這個結(jié)構(gòu)有一個隱藏約定v 數(shù)組的排列順序必須始終是逆時針。所有幾何謂詞和鄰居關(guān)系都依賴這個約定。一旦某個三角形變成順時針后面的翻邊條件會整體取反產(chǎn)生難以排查的鏡像網(wǎng)格。3.2 orient2d 與 inCircle兩個決定正確性的幾何謂詞方向測試orient2d是全部幾何判斷的地基。它計算的是向量 AB 和 AC 的叉積double orient2d(const Point a, const Point b, const Point c) { return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x); }返回值大于零表示 c 在 b 的左側(cè)即三角形 (a,b,c) 為逆時針。定位、翻邊前的凸四邊形檢查、超三角形頂點選擇都依賴這個函數(shù)的符號。外接圓判斷inCircle是另一個核心謂詞。它判斷第四個點 d 落在三角形 abc 的外接圓內(nèi)還是圓外。直接計算圓心半徑再比較距離當然可行但更高效的做法是構(gòu)造一個四階行列式double inCircle(const Point a, const Point b, const Point c, const Point d) { double adx a.x - d.x, ady a.y - d.y; double bdx b.x - d.x, bdy b.y - d.y; double cdx c.x - d.x, cdy c.y - d.y; double abd adx * bdy - ady * bdx; double bcd bdx * cdy - bdy * cdx; double cad cdx * ady - cdy * adx; double alift adx * adx ady * ady; double blift bdx * bdx bdy * bdy; double clift cdx * cdx cdy * cdy; return alift * bcd blift * cad clift * abd; }先做坐標平移把 d 作為原點能有效減小坐標值的絕對大小從而降低浮點誤差。符號約定是a、b、c 逆時針排列時返回正表示 d 在圓外返回負表示在圓內(nèi)。這兩個謂詞的實現(xiàn)都不復(fù)雜但它們共同決定了算法會不會在某個角落靜默地生成非法三角形。3.3 拆分三角形時的鄰居接線插入新點 p 到三角形 ABC 內(nèi)部時舊三角形被替換成三個新三角形。最清晰的實現(xiàn)方式是先算出舊三角形三個頂點的索引再按固定模式創(chuàng)建新三角形void insertPoint(int pi) { int t0 locatePoint(pts, tris, pts[pi], lastTriangle); const Tri old tris[t0]; int a old.v[0], b old.v[1], c old.v[2]; // 新三角形 tA 的頂點順序為 (a, b, pi) // 新三角形 tB 的頂點順序為 (b, c, pi) // 新三角形 tC 的頂點順序為 (c, a, pi) // 每個新三角形的鄰居需要同時連接舊鄰居和另外兩個新三角形 }創(chuàng)建時最需要注意的是內(nèi)部接線三個新三角形之間互連形成一個小環(huán)。例如 tA 的邊 (b, pi) 對面是 tB邊 (pi, a) 對面是 tC。手工寫這三個三角形的鄰居數(shù)組極易出錯一個常見的做法是先統(tǒng)一置 -1再用一個輔助函數(shù)linkTri把兩個三角形在指定邊上互相指向?qū)Ψ絭oid linkTri(int t1, int e1, int t2, int e2) { tris[t1].n[e1] t2; tris[t2].n[e2] t1; }然后把拆分時的所有鄰居關(guān)系寫成一行行的linkTri調(diào)用。代碼雖然長但每行只做一件事出錯了也容易定位。這類“接線錯誤”不會立刻報錯而是會在后續(xù)翻邊或輸出網(wǎng)格時以拓撲斷裂的形式出現(xiàn)排查起來非常耗時。3.4 legalize 與 flip 的源碼實現(xiàn)legalize負責(zé)檢查和修復(fù)一條共享邊。它的輸入是一個三角形索引和一條邊的索引int thirdVertex(const Tri tri, int u, int v) { for (int i 0; i 3; i) { if (tri.v[i] ! u tri.v[i] ! v) return tri.v[i]; } return -1; } void legalize(int t, int e) { int nb tris[t].n[e]; if (nb 0) return; // 邊界邊直接返回 int u tris[t].v[e]; int v tris[t].v[(e 1) % 3]; int w tris[t].v[(e 2) % 3]; int z thirdVertex(tris[nb], u, v); if (inCircle(pts[u], pts[v], pts[w], pts[z]) 0) { flip(t, e, nb, u, v, w, z); // 需要翻邊 } }flip把共享邊為 uv 的兩個三角形第三頂點分別是 w 和 z替換成共享邊為 wz 的兩個新三角形。實現(xiàn)時要特別注意新三角形的頂點必須重新排成逆時針并且四個鄰居關(guān)系要整體重接。另一種更穩(wěn)妥的寫法是“先斷后連”先把四個相關(guān)三角形的鄰居索引清零再按新拓撲重新linkTri。少做一步殘留的舊鄰居索引就會讓下一次legalize讀到錯誤的三角形。翻邊之后新形成的兩條邊可能又和它們的鄰居發(fā)生沖突所以flip內(nèi)部還需要對兩條新邊遞歸調(diào)用legalize。這與 2.3 節(jié)說的級聯(lián)翻邊是對應(yīng)的。遞歸深度由最小角增大性質(zhì)保證有界不會棧溢出。3.5 超級三角形與收尾清理增量算法需要一塊初始地盤作為定位起點這就是超級三角形。它要覆蓋所有輸入點最常用的構(gòu)造方式是先算輸入點集的包圍盒把中心作為三角形重心再向外擴到對角線長度的兩倍左右。Point minP pts[0], maxP pts[0]; for (auto p : pts) { minP.x std::min(minP.x, p.x); minP.y std::min(minP.y, p.y); maxP.x std::max(maxP.x, p.x); maxP.y std::max(maxP.y, p.y); } double d std::sqrt((maxP.x - minP.x) * (maxP.x - minP.x) (maxP.y - minP.y) * (maxP.y - minP.y)); Point c { (minP.x maxP.x) / 2.0, (minP.y maxP.y) / 2.0 }; // 三個頂點圍繞 c半徑取 2d 左右超級三角形的頂點索引一般從 n 開始保存n 是原始輸入點數(shù)。當所有點都插入完畢后任何包含超級三角形頂點的三角形都必須刪除剩下的就是最終 Delaunay 網(wǎng)格。int n pts.size(); for (auto it tris.begin(); it ! tris.end(); ) { if (it-v[0] n || it-v[1] n || it-v[2] n) { it tris.erase(it); } else { it; } }這一步看似簡單但它隱藏了一個細節(jié)刪除超級三角形后原來與超級三角形相接的邊界邊會暴露出來成為凸包邊界。如果插入順序或翻邊過程有過任何輕微錯誤這里往往會出現(xiàn)外掛的碎片三角形一眼就能看出網(wǎng)格不對。4. 避坑指南增量 Delaunay 最容易翻車的五個細節(jié)4.1 inCircle 符號約定與頂點繞序不一致現(xiàn)象網(wǎng)格生成過程幾乎正常但輸出結(jié)果像是被鏡像了部分三角形重疊或者某些內(nèi)部邊缺失。原因inCircle返回值的符號依賴三角形頂點的繞序。三角形按逆時針存儲時返回負表示第四點在圓內(nèi)如果某個三角形的頂點被存成順時針符號恰好互換翻邊條件整體取反等于一直在維護“局部非 Delaunay”網(wǎng)格。解決在每次創(chuàng)建三角形后加一個 debug 斷言檢查orient2d(pts[v0], pts[v1], pts[v2]) 0。只要有一個三角形違反逆時針約定立刻終止程序。這種斷言只在排查時開release 里可以關(guān)閉但它能攔截掉一大類隱蔽的接線錯誤。4.2 共圓點inCircle 恰好等于零現(xiàn)象某個局部區(qū)域在兩種對角線之間反復(fù)橫跳legalize遞歸不停止或者同一份輸入在不同優(yōu)化等級下得到不同結(jié)果。原因四個點精確共圓時inCircle 行列式值為零翻邊邏輯沒有定義。真實數(shù)據(jù)中完全共圓很少見但對稱的柵格點、人工生成的測試點很容易觸發(fā)。零值本身不可怕可怕的是它在不同編譯器或優(yōu)化設(shè)置下的浮點表現(xiàn)不一致。解決最常見的是引入符號擾動symbolic perturbation給每個點一個與坐標相關(guān)的微小虛擬偏移打破共圓對稱性。工程上更省事的做法是對輸入做輕量去對稱化并在inCircle返回絕對值小于某個極小閾值時默認不翻邊。對正確性要求高的場景應(yīng)該換用精確自適應(yīng)精度謂詞而不是依賴 double 的運氣。4.3 重復(fù)點與共線點沒預(yù)處理現(xiàn)象剖分過程中出現(xiàn) NaN定位死循環(huán)甚至數(shù)組越界訪問。原因重復(fù)點會讓兩個頂點在幾何上完全重合inCircle 距離項退化行列式信息丟失共線點則讓“點在三角形內(nèi)部”的判定退化為“點在線上”orient 值反復(fù)為零步行法在相鄰三角形之間不停擺動無法收斂。解決插入前先按坐標排序再去重。排序可以用std::sort配合自定義比較器去重后用唯一序列重新建點集。共線點屬于另一種退化不能硬塞進同一個剖分流程需要根據(jù)業(yè)務(wù)規(guī)則做降維處理或單獨分段。數(shù)據(jù)清洗是這一步的硬要求僥幸心理基本都會以神秘崩潰收場。4.4 超級三角形選太大現(xiàn)象靠近超級三角形的三角形面積異常大整體質(zhì)量統(tǒng)計被拉低某些翻邊操作產(chǎn)生極細長三角形。原因超級三角形頂點坐標比輸入點大幾個數(shù)量級時大三角形上的外接圓半徑也成倍放大。inCircle 里坐標差的平方項變成“大數(shù)減大數(shù)”浮點有效位數(shù)被大量消耗翻邊判斷變得不可靠。解決超三角形的半徑取包圍盒對角線長度的兩倍左右而不是拍腦袋取一個 1000 或 10000。另一個有效的補充手段是在插入前把所有點平移到以包圍盒中心為原點降低絕對坐標數(shù)值。幾何變換不改變拓撲結(jié)構(gòu)平移對整個剖分結(jié)果零影響卻能讓浮點誤差明顯變小。4.5 步行法在邊界點上的振蕩現(xiàn)象定位函數(shù)觸發(fā) guard 退出返回一個錯誤的三角形索引后續(xù)翻邊全部跑偏。原因目標點恰好在兩個三角形的公共邊上時orient 結(jié)果在兩個微小正負之間搖擺兩次判斷分別選中不同的跨越方向于是步行法在相鄰三角形之間來回跳動無法收斂。解決對 orient 的結(jié)果做鉗制當絕對值小于極小閾值時直接置為零并把“點在邊上”按“點在三角形內(nèi)”處理。閾值不能定得太大否則會把本應(yīng)通過跨越離開的鄰居點誤判成內(nèi)部點。一個實用的經(jīng)驗值是按輸入坐標量級的 1e-12 到 1e-10 來做鉗制。5. 設(shè)計取舍增量翻邊、Bowyer-Watson 與分治法的路線選擇增量翻邊并不是唯一的 Delaunay 構(gòu)建路線。另一種思路 Bowyer-Watson 算法會先找出所有外接圓包含新點的三角形把它們批量刪除再把留下的空洞與新點連成新三角形。它比增量翻邊更直觀但每次插入都需要維護一個候選集合刪除和重建的數(shù)量更多常數(shù)更大。分治法能在最壞情況下穩(wěn)定做到 O(n log n)但合并兩個子剖分時要處理復(fù)雜的凸包邊界縫合實現(xiàn)量遠大于增量法而且完全無法支持動態(tài)加點。增量法最大的工程優(yōu)勢在于局部性和動態(tài)性。局部性體現(xiàn)在每次插入只需處理與新點相鄰的少量三角形動態(tài)性體現(xiàn)在它天然支持流式添加點后續(xù)任意時刻都能繼續(xù)插入。代價是復(fù)雜度對輸入順序敏感隨機順序插入時步行法定位和翻邊次數(shù)的期望都很友好整體期望 O(n log n)如果輸入點按某種特定順序排列定位可能要穿越大半個網(wǎng)格退回 O(n^2)。所以很多實現(xiàn)會強制對點集做隨機打亂犧牲輸出順序的確定性來換取穩(wěn)定的性能。有一次我處理某公司提供的逐行掃描地形數(shù)據(jù)數(shù)十萬個點按柵格順序排列直接喂進增量實現(xiàn)生成時間從預(yù)期的幾秒膨脹到兩分多鐘。排查了很久才意識到不是算法寫錯了而是輸入順序太規(guī)整規(guī)整到每次插入的定位都要從上一個點所在的三角形跨過大半個網(wǎng)格。從那以后我每次跑增量 Delaunay 之前都會先檢查輸入點的排列習(xí)慣必要時加一次 shuffle同時把 inCircle 和 orient2d 的符號約定用斷言鎖死杜絕繞序問題在半夜偷襲。三角剖分本身不難難的是讓每一次局部修復(fù)始終聰明地工作。希望這篇拆解能幫到你把這套增量機制的底層邏輯理順。本文還有配套的精品資源點擊獲取