簡(jiǎn)介本資源是一套基于Bowyer-Watson算法實(shí)現(xiàn)Delaunay三角剖分與Voronoi圖生成的Matlab完整仿真方案面向本科及碩士階段的科研學(xué)習(xí)者適用于計(jì)算幾何、空間分析、地理信息系統(tǒng)、路徑規(guī)劃及圖像處理等方向的基礎(chǔ)建模需求。壓縮包共含4個(gè)文件3張結(jié)果可視化PNG圖1個(gè)核心m腳本總大小僅92KB輕量易部署適配Matlab 2014a/2019a環(huán)境內(nèi)含可直接運(yùn)行的代碼及對(duì)應(yīng)效果圖便于理解算法流程與幾何結(jié)構(gòu)關(guān)系。目前已有132人下載學(xué)習(xí)資源聚焦經(jīng)典計(jì)算幾何問(wèn)題提供從點(diǎn)集輸入、增量構(gòu)網(wǎng)、外接圓判據(jù)到對(duì)偶圖轉(zhuǎn)換的全流程實(shí)現(xiàn)代碼注釋清晰、邏輯分層明確特別適合初學(xué)者掌握Delaunay/Voronoi內(nèi)在關(guān)聯(lián)與Matlab向量化編程技巧。 做網(wǎng)格生成和空間分析這幾年Delaunay三角網(wǎng)絡(luò)和Voronoi泰森多邊形一直是我繞不開(kāi)的兩個(gè)東西。最近整理舊項(xiàng)目時(shí)翻到一套用MATLAB寫(xiě)的代碼基于Bowyer-Watson算法從零構(gòu)建這兩個(gè)結(jié)構(gòu)代碼雖然不算長(zhǎng)但把計(jì)算幾何里很多經(jīng)典細(xì)節(jié)都踩了一遍。如果你也在做點(diǎn)集剖分、GIS鄰近分析、區(qū)域劃分或者只是想在MATLAB里徹底搞懂Delaunay/Voronoi的底層邏輯這篇應(yīng)該能幫你省下不少折騰時(shí)間。1. 項(xiàng)目背景與核心思路Delaunay和Voronoi到底在算什么1.1 這兩個(gè)幾何結(jié)構(gòu)為什么總是一起出現(xiàn)Delaunay三角網(wǎng)絡(luò)和Voronoi泰森多邊形本質(zhì)上是一對(duì)對(duì)偶結(jié)構(gòu)。Delaunay三角網(wǎng)滿(mǎn)足一個(gè)很關(guān)鍵的性質(zhì)每個(gè)三角形的外接圓內(nèi)不包含其他點(diǎn)稱(chēng)為空外接圓性質(zhì)。而Voronoi多邊形是把平面劃分成若干區(qū)域使得每個(gè)區(qū)域內(nèi)的點(diǎn)到該區(qū)域控制點(diǎn)的距離最近。這兩個(gè)結(jié)構(gòu)的對(duì)偶關(guān)系體現(xiàn)在Voronoi圖的每條邊恰好對(duì)應(yīng)Delaunay三角網(wǎng)的一條邊而且是垂直平分的關(guān)系Voronoi圖的每個(gè)頂點(diǎn)恰好是某個(gè)Delaunay三角形的外心。這個(gè)對(duì)偶關(guān)系在實(shí)際項(xiàng)目里非常有用。比如做無(wú)線(xiàn)基站覆蓋分析你想知道每個(gè)基站管轄哪片區(qū)域直接算Voronoi做地形建模、有限元網(wǎng)格劃分你需要一個(gè)不產(chǎn)生細(xì)長(zhǎng)三角形的剖分直接上Delaunay。很多場(chǎng)景下算出了Delaunay基本就等于拿到了Voronoi反過(guò)來(lái)也一樣。這也是為什么我在項(xiàng)目里把這兩個(gè)東西放在一起實(shí)現(xiàn)一次構(gòu)建兩邊受益。1.2 Bowyer-Watson算法的三步核心Bowyer-Watson算法是逐點(diǎn)插入法里最經(jīng)典的一種核心思路可以拆成三步先構(gòu)造一個(gè)足夠大的超三角形把點(diǎn)集中所有點(diǎn)都包含進(jìn)去保證初始狀態(tài)下存在一個(gè)可以合法插入點(diǎn)的基礎(chǔ)三角網(wǎng)。逐個(gè)插入點(diǎn)P找到所有外接圓包含P的三角形這些三角形被稱(chēng)為壞三角形。把它們?nèi)縿h除留下一個(gè)多邊形的空洞這個(gè)空洞叫影響空腔。用P和空腔邊界上的所有邊逐條連接生成新的三角形恢復(fù)三角網(wǎng)結(jié)構(gòu)。循環(huán)執(zhí)行第2、3步直到所有點(diǎn)插入完畢最后再把所有包含超三角形頂點(diǎn)的三角形刪除。整個(gè)過(guò)程看起來(lái)不復(fù)雜但每一步都有隱含的邊界情況超三角形取多大才不會(huì)影響最終結(jié)果壞三角形的判斷用什么精度空腔邊界如何提取才算干凈去重和退化怎么處理。第四章里我會(huì)逐段拆MATLAB代碼把這些問(wèn)題都攤開(kāi)來(lái)說(shuō)。1.3 從Delaunay到Voronoi的映射關(guān)系因?yàn)镈elaunay和Voronoi對(duì)偶所以生成Voronoi其實(shí)不需要重新跑算法。我采用的做法是先完整構(gòu)建Delaunay三角網(wǎng)然后對(duì)每個(gè)三角形求外心。這樣一來(lái)每個(gè)Delaunay三角形對(duì)應(yīng)一個(gè)Voronoi頂點(diǎn)每條Delaunay內(nèi)部邊對(duì)應(yīng)一條Voronoi邊每個(gè)Delaunay點(diǎn)對(duì)應(yīng)一個(gè)Voronoi胞元。實(shí)際操作時(shí)只需要遍歷三角形列表計(jì)算外心再根據(jù)三角形之間的鄰接關(guān)系把共享同一條邊的兩個(gè)三角形的外心連起來(lái)就是一條Voronoi邊。構(gòu)建每個(gè)點(diǎn)的胞元時(shí)找出所有以該點(diǎn)為頂點(diǎn)的三角形把這些三角形的外心按角度排序并連線(xiàn)就得到了該點(diǎn)對(duì)應(yīng)的泰森多邊形。這個(gè)映射關(guān)系看起來(lái)很直接但真正讓代碼“跑得穩(wěn)”的細(xì)節(jié)全在外心和鄰接關(guān)系的計(jì)算上下文會(huì)詳細(xì)講。2. 動(dòng)手寫(xiě)代碼前必須先做的三個(gè)設(shè)計(jì)決策2.1 點(diǎn)集和三角形用什么數(shù)據(jù)結(jié)構(gòu)最順手寫(xiě)MATLAB代碼前首先要把數(shù)據(jù)結(jié)構(gòu)定下來(lái)否則代碼一長(zhǎng)就會(huì)亂。我的做法是點(diǎn)集用N×2的矩陣P存儲(chǔ)每行是一個(gè)點(diǎn)的x、y坐標(biāo)坐標(biāo)編號(hào)從1到N。三角形用一個(gè)Ntri×3的矩陣triList存儲(chǔ)每一行存三個(gè)點(diǎn)的索引代表一個(gè)三角形。注意這里存的是索引不是坐標(biāo)因?yàn)樗饕鋈ブ亍⑴判?、鄰接查找時(shí)都很快。為什么不用cell數(shù)組存點(diǎn)坐標(biāo)因?yàn)樗饕僮骺梢员苊饷看蝿?chuàng)建新數(shù)組時(shí)的拷貝開(kāi)銷(xiāo)尤其在逐點(diǎn)插入過(guò)程中壞三角形的刪除和新三角形的添加會(huì)頻繁發(fā)生。我用triList矩陣配合邏輯索引能夠用一行代碼過(guò)濾掉壞三角形非常符合MATLAB的向量化風(fēng)格。另外還有一個(gè)用途最終刪除超三角形時(shí)只需要檢查三角形三個(gè)頂點(diǎn)索引是否大于N就能快速篩掉比比較坐標(biāo)靠譜得多。2.2 超三角形選多大才算“足夠大”超三角形是整個(gè)Bowyer-Watson算法的啟動(dòng)條件。如果它不夠大會(huì)有兩種情況要么是某些點(diǎn)落在超三角形外部導(dǎo)致初始剖分不合法要么是三角形外接圓覆蓋范圍不夠插入點(diǎn)時(shí)會(huì)漏掉部分壞三角形。這兩種情況都會(huì)讓最終三角網(wǎng)出現(xiàn)撕裂或重復(fù)三角形。我習(xí)慣這么取超三角形頂點(diǎn)先計(jì)算點(diǎn)集的包圍盒取中心點(diǎn)c和包圍盒對(duì)角線(xiàn)長(zhǎng)度d然后以c為中心構(gòu)造一個(gè)邊長(zhǎng)約為10倍d的等邊三角形。把三個(gè)頂點(diǎn)放在離點(diǎn)群足夠遠(yuǎn)的位置確保所有點(diǎn)包括它們的極端外接圓都被包含在內(nèi)。需要提醒的是等邊三角形優(yōu)于直角或鈍角三角形因?yàn)樗耐饨訄A半徑和邊長(zhǎng)比值較小能減少后續(xù)計(jì)算外接圓時(shí)的浮點(diǎn)誤差。超三角形也不是越大越好。太大了會(huì)讓初始三角形外接圓半徑變得非常大判斷壞三角形時(shí)就會(huì)出現(xiàn)大數(shù)吃小數(shù)的問(wèn)題導(dǎo)致外心坐標(biāo)誤差明顯增大。我實(shí)測(cè)下來(lái)10倍對(duì)角線(xiàn)是個(gè)比較穩(wěn)的經(jīng)驗(yàn)值既不干擾最終結(jié)果也不至于讓數(shù)值精度失控。2.3 外接圓判斷的幾何式與浮點(diǎn)精度Bowyer-Watson的壞三角形判斷核心是判斷一個(gè)點(diǎn)P是否落在三角形ABC的外接圓內(nèi)部。求外心坐標(biāo)最直接的方法是解線(xiàn)性方程組也可以用下列公式function [cx, cy, r2] circumcircle(A, B, C) % A、B、C是兩行或三行坐標(biāo)點(diǎn)這里按2D點(diǎn)處理 d 2 * (A(1)*(B(2)-C(2)) B(1)*(C(2)-A(2)) C(1)*(A(2)-B(2))); if abs(d) 1e-12 cx inf; cy inf; r2 inf; return; end ax A(1); ay A(2); bx B(1); by B(2); cx_ C(1); cy_ C(2); a2 ax*ax ay*ay; b2 bx*bx by*by; c2 cx_*cx_ cy_*cy_; ux (a2*(by-cy_) b2*(cy_-ay) c2*(ay-by)) / d; uy (a2*(cx_-bx) b2*(ax-cx_) c2*(bx-ax)) / d; cx ux; cy uy; r2 (ax-ux)^2 (ay-uy)^2; end判斷P是否在圓內(nèi)時(shí)不能直接用dist(P, center) r因?yàn)楦↑c(diǎn)運(yùn)算中邊界情況極多。我統(tǒng)一用平方距離比較dist2 r2 tol其中tol取一個(gè)與點(diǎn)集尺度相關(guān)的值比如點(diǎn)集包圍盒對(duì)角線(xiàn)長(zhǎng)度的1e-9倍。引入容差后位于圓上的點(diǎn)也能穩(wěn)定地進(jìn)入刪除流程避免因?yàn)?.0000001的誤差導(dǎo)致插入失敗這是我在實(shí)際調(diào)試中踩過(guò)的最典型的坑。3. MATLAB核心實(shí)現(xiàn)Bowyer-Watson算法逐段拆解3.1 初始化讀入點(diǎn)集、去重、構(gòu)造超三角形寫(xiě)了一個(gè)主函數(shù)delaunay_voronoi_demo.m整個(gè)流程從隨機(jī)生成測(cè)試點(diǎn)開(kāi)始。實(shí)際項(xiàng)目里可以直接替換成自己的坐標(biāo)點(diǎn)。% 生成隨機(jī)測(cè)試點(diǎn)稍微加點(diǎn)簇狀分布更接近真實(shí)場(chǎng)景 rng(42); N 60; P [randn(N,1)*0.3 0.5, randn(N,1)*0.3 0.5]; P [P; rand(N,1)*3 1, rand(N,1)*3 1]; P unique(round(P, 10), rows, stable); N size(P, 1);這里unique去重很關(guān)鍵。如果你的數(shù)據(jù)里存在重復(fù)點(diǎn)Bowyer-Watson會(huì)在插入重復(fù)點(diǎn)時(shí)產(chǎn)生零面積三角形進(jìn)而讓外接圓計(jì)算出現(xiàn)分母接近0的情況后面所有判斷都會(huì)失效。round(P,10)是為了避免浮點(diǎn)運(yùn)算導(dǎo)致“同一坐標(biāo)但差一位小數(shù)”的點(diǎn)沒(méi)被當(dāng)成重復(fù)點(diǎn)。構(gòu)造超三角形我單獨(dú)抽了一個(gè)子函數(shù)保持主流程清晰function T superTriangle(P) xmin min(P(:,1)); xmax max(P(:,1)); ymin min(P(:,2)); ymax max(P(:,2)); cx (xmin xmax) / 2; cy (ymin ymax) / 2; d sqrt((xmax - xmin)^2 (ymax - ymin)^2); R 10 * d; T [cx, cy 2*R; cx - sqrt(3)*R, cy - R; cx sqrt(3)*R, cy - R]; end用中心點(diǎn)和固定半徑構(gòu)造等邊三角形是為了讓三個(gè)頂點(diǎn)關(guān)于點(diǎn)集中心對(duì)稱(chēng)分布從源頭上減小外心計(jì)算的病態(tài)程度。這一步不需要太精確但要記好最終結(jié)果里包含這三個(gè)頂點(diǎn)的三角形都要?jiǎng)h掉。3.2 主循環(huán)逐點(diǎn)插入、找壞三角形、重建局部三角網(wǎng)主循環(huán)是整個(gè)算法的核心我把它拆成了幾個(gè)清晰的步驟避免在一大坨代碼里繞暈。% 初始化把超三角形頂點(diǎn)也加入坐標(biāo)矩陣統(tǒng)一編號(hào) T superTriangle(P); P_all [P; T]; triList [N1, N2, N3]; % 初始三角形 for i 1:N p P(i,:); % 1. 找出所有外接圓包含p的三角形 [cx, cy, r2] arrayfun((t) triangleCircum(P_all, triList(t,:)), 1:size(triList,1)); dist2 (cx - p(1)).^2 (cy - p(2)).^2; tol 1e-9 * (max(P_all(:,1)) - min(P_all(:,1))); badIdx find(dist2 r2 tol); if isempty(badIdx) error(當(dāng)前點(diǎn)未落在任何三角形外接圓內(nèi)超三角形可能不夠大); end badTri triList(badIdx, :); triList(badIdx, :) []; % 2. 提取影響空腔邊界統(tǒng)計(jì)每條邊出現(xiàn)的次數(shù) edges [badTri(:,1) badTri(:,2); badTri(:,2) badTri(:,3); badTri(:,3) badTri(:,1)]; edges sort(edges, 2); [uniqueEdges, ~, ic] unique(edges, rows, stable); edgeCount accumarray(ic, 1); boundaryEdges uniqueEdges(edgeCount 1, :); % 3. 連接新點(diǎn)與邊界邊生成新三角形 newTris [boundaryEdges(:,1), boundaryEdges(:,2), repmat(i, size(boundaryEdges,1), 1)]; triList [triList; newTris]; end這段代碼里最值得講的是邊界提取。刪除壞三角形后留下的空腔邊界本質(zhì)上是這些壞三角形里出現(xiàn)了奇數(shù)次的邊。出現(xiàn)兩次的邊就是兩個(gè)壞三角形共享的內(nèi)部邊刪除三角形后它不應(yīng)該再存在只出現(xiàn)一次的邊才是空腔的真正外邊界。我用sort把邊的頂點(diǎn)排序再用unique配合accumarray統(tǒng)計(jì)每條邊的出現(xiàn)次數(shù)一行代碼就把內(nèi)部邊和邊界邊分離開(kāi)了。這個(gè)方法比逐個(gè)鄰接遍歷高效得多而且在MATLAB里屬于標(biāo)準(zhǔn)操作。循環(huán)結(jié)束后還差一步刪除所有包含超三角形頂點(diǎn)的三角形。triList(any(triList N, 2), :) [];超過(guò)N的索引全部是超三角形的三個(gè)頂點(diǎn)之一直接刪掉即可。到這一步一個(gè)不依賴(lài)MATLAB任何內(nèi)置剖分函數(shù)的Delaunay三角網(wǎng)就構(gòu)建完成了。3.3 去重和邊界檢查別讓重復(fù)三角形毀了結(jié)果很多人寫(xiě)B(tài)owyer-Watson時(shí)會(huì)遇到最后結(jié)果里出現(xiàn)重復(fù)三角形或者三角形三條邊共線(xiàn)的情況這通常不是算法邏輯問(wèn)題而是點(diǎn)集去重不徹底、容差設(shè)置不合適導(dǎo)致的。我在循環(huán)前已經(jīng)對(duì)原始點(diǎn)做了unique去重但插入過(guò)程中因?yàn)槿莶钤蛉匀挥锌赡苌擅娣e接近0的三角形。所以在三角形生成后我加了一個(gè)兜底檢查過(guò)濾掉面積小于某閾值的三角形。A P_all(triList(:,1), :); B P_all(triList(:,2), :); C P_all(triList(:,3), :); area2 abs((B(:,1)-A(:,1)).*(C(:,2)-A(:,2)) - (B(:,2)-A(:,2)).*(C(:,1)-A(:,1))); triList(area2 1e-12, :) [];這個(gè)過(guò)濾不能隨便刪因?yàn)槊娣e過(guò)小的三角形雖然不影響后續(xù)Voronoi的生成但會(huì)在可視化時(shí)畫(huà)出一堆肉眼幾乎看不見(jiàn)的碎邊影響圖形輸出的美觀(guān)度。更重要的是它會(huì)降低后續(xù)鄰接查找的效率。3.4 生成Voronoi從三角形外心到泰森多邊形Delaunay三角網(wǎng)輸出得到后Voronoi就省力了。先給每個(gè)三角形算外心再聚合到每個(gè)原始點(diǎn)上。% 計(jì)算每個(gè)Delaunay三角形外心 m size(triList, 1); vorVertices zeros(m, 2); for t 1:m [vorVertices(t,1), vorVertices(t,2)] triangleCircum(P_all, triList(t,:)); end如果只需要Voronoi邊可以直接遍歷每條內(nèi)部Delaunay邊找到共享邊的兩個(gè)三角形把兩個(gè)外心連線(xiàn)。但我在項(xiàng)目里更多是需要“每個(gè)點(diǎn)的泰森多邊形”所以我按點(diǎn)聚合所有鄰接三角形的外心vorCells cell(N, 1); for i 1:N triIdx find(any(triList i, 2)); if isempty(triIdx) continue; end cellVerts vorVertices(triIdx, :); % 按角度排序保證多邊形頂點(diǎn)順序正確 angles atan2(cellVerts(:,2) - mean(cellVerts(:,2)), cellVerts(:,1) - mean(cellVerts(:,1))); [~, order] sort(angles); vorCells{i} cellVerts(order, :); end這里有一點(diǎn)要注意Voronoi胞元在凸包邊界處是無(wú)限區(qū)域外心可能跑到無(wú)窮遠(yuǎn)。我在可視化時(shí)會(huì)把邊界外的點(diǎn)裁剪到畫(huà)布范圍內(nèi)MATLAB的axis命令雖然可以顯示但裁剪操作最好自己做否則連線(xiàn)會(huì)畫(huà)出很長(zhǎng)的斜線(xiàn)。第四章可視化部分我會(huì)演示如何既保留邊界射線(xiàn)的效果又不破壞畫(huà)面。4. 實(shí)操記錄從零跑通整個(gè)流程4.1 完整的函數(shù)組織方式寫(xiě)代碼的時(shí)候我把功能拆成三個(gè)文件主腳本、核心算法函數(shù)、工具函數(shù)這樣調(diào)試和維護(hù)都方便。主腳本負(fù)責(zé)生成點(diǎn)、調(diào)用算法、畫(huà)圖核心算法函數(shù)delaunayByBowyerWatson負(fù)責(zé)構(gòu)建三角網(wǎng)工具函數(shù)circumcircle和superTriangle放在同一個(gè)文件末尾作為嵌套函數(shù)或者單獨(dú)拆出來(lái)都行。如果你的點(diǎn)集有幾個(gè)萬(wàn)級(jí)別逐點(diǎn)插入配for循環(huán)在MATLAB里會(huì)明顯變慢。這時(shí)候有兩個(gè)方向一是把壞三角形查找向量化像我在第三章里那樣一次性計(jì)算所有三角形外接圓而不是再套一層for二是考慮加入網(wǎng)格空間索引減少每次插入時(shí)掃描的三角形數(shù)量。對(duì)于幾千點(diǎn)的規(guī)模來(lái)說(shuō)向量化已經(jīng)足夠我的實(shí)測(cè)結(jié)果是在普通筆記本上60個(gè)點(diǎn)可以瞬間完成2000個(gè)點(diǎn)大約需要0.3到1秒完全可以用在項(xiàng)目調(diào)試中。4.2 可視化把Delaunay和Voronoi畫(huà)在一張圖上畫(huà)圖是檢驗(yàn)算法正確性最直觀(guān)的手段。我習(xí)慣把Delaunay三角網(wǎng)用淺色線(xiàn)畫(huà)出來(lái)Voronoi多邊形用另一條顏色畫(huà)在同一張圖上能很清楚地看到對(duì)偶關(guān)系。figure; hold on; % 畫(huà)Delaunay三角網(wǎng) triplot(triList, P_all(:,1), P_all(:,2), Color, [0.6 0.6 0.6], LineWidth, 0.5); % 畫(huà)Voronoi for i 1:N if ~isempty(vorCells{i}) patch(vorCells{i}(:,1), vorCells{i}(:,2), w, EdgeColor, [0.8 0.2 0.2], LineWidth, 1.2); end end plot(P(:,1), P(:,2), ko, MarkerFaceColor, k, MarkerSize, 4); axis equal;一個(gè)非常容易踩的坑是axis equal。如果沒(méi)有等比例坐標(biāo)軸Voronoi的多邊形形狀會(huì)被拉伸看起來(lái)像是算法錯(cuò)誤實(shí)際上只是顯示比例問(wèn)題。我第一次排查了半天最后發(fā)現(xiàn)是坐標(biāo)軸沒(méi)等比例白白浪費(fèi)了一個(gè)小時(shí)。4.3 對(duì)比驗(yàn)證用MATLAB內(nèi)置函數(shù)當(dāng)“裁判”自己寫(xiě)的算法對(duì)不對(duì)最直接的辦法是拿MATLAB自帶的delaunayTriangulation做交叉驗(yàn)證。雖然項(xiàng)目目的是練習(xí)底層實(shí)現(xiàn)但驗(yàn)證環(huán)節(jié)不能省。DT delaunayTriangulation(P); triBuiltin DT.ConnectivityList; % 比較三角形數(shù)量是否一致 fprintf(自定義 Delaunay 三角形數(shù)量: %d\n, size(triList, 1)); fprintf(內(nèi)置 Delaunay 三角形數(shù)量: %d\n, size(triBuiltin, 1));注意三角形數(shù)量一致不代表兩個(gè)剖分完全一樣頂點(diǎn)索引順序和三角形集合可能有微小差異這是Delaunay剖分不唯一導(dǎo)致的正常情況。更實(shí)際的驗(yàn)證方式是檢查空外接圓性質(zhì)和邊界數(shù)量所有三角形外接圓內(nèi)不含其它點(diǎn)凸包邊界上的邊數(shù)量等于凸包頂點(diǎn)數(shù)。我寫(xiě)了一個(gè)自動(dòng)化檢查腳本每個(gè)點(diǎn)插入后都斷言一下當(dāng)前三角網(wǎng)的三角形數(shù)量滿(mǎn)足歐拉公式這樣能盡早暴露問(wèn)題而不是等最后畫(huà)圖才發(fā)現(xiàn)。5. 典型問(wèn)題與排查實(shí)錄5.1 外接圓誤判浮點(diǎn)容差到底怎么設(shè)置才合理這是整個(gè)項(xiàng)目里我踩得最深的坑。一開(kāi)始我把容差寫(xiě)成了固定值1e-10結(jié)果點(diǎn)集坐標(biāo)范圍比較大的時(shí)候比如坐標(biāo)在10000量級(jí)1e-10完全被浮點(diǎn)精度淹沒(méi)導(dǎo)致一部分恰好在圓上的點(diǎn)沒(méi)有被判定為壞三角形最后生成Delaunay網(wǎng)格時(shí)出現(xiàn)明顯的孔洞。后來(lái)我把容差改成相對(duì)值與點(diǎn)集包圍盒對(duì)角線(xiàn)長(zhǎng)度掛鉤tol 1e-9 * diagLen。diagLen是點(diǎn)集范圍的量級(jí)這樣做的好處是無(wú)論你的數(shù)據(jù)是單位正方形內(nèi)的小規(guī)模坐標(biāo)還是幾十萬(wàn)坐標(biāo)值的地理坐標(biāo)容差都能跟著縮放。建議讀者在實(shí)現(xiàn)時(shí)一定不要用固定絕對(duì)值容差這算是一條通用經(jīng)驗(yàn)。5.2 四點(diǎn)共圓退化情況下會(huì)出現(xiàn)非唯一剖分當(dāng)四個(gè)點(diǎn)恰好共圓時(shí)Delaunay剖分并不唯一。Bowyer-Watson在這種情況下會(huì)隨機(jī)選擇兩條對(duì)角線(xiàn)中的一條不能說(shuō)算法錯(cuò)誤但可能和你預(yù)期的結(jié)果不一致而且有時(shí)會(huì)在相鄰位置產(chǎn)生一個(gè)極其扁平的三角形看起來(lái)像是bug。我處理這類(lèi)問(wèn)題的方式是加入微小擾動(dòng)或者接受非唯一剖分。如果點(diǎn)集來(lái)自實(shí)際數(shù)據(jù)精確保存了原始坐標(biāo)那就接受算法給出的結(jié)果如果點(diǎn)集是人工生成的比如規(guī)則網(wǎng)格點(diǎn)我會(huì)在預(yù)處理階段給每個(gè)點(diǎn)加一個(gè)約1e-8量級(jí)的隨機(jī)擾動(dòng)打破共圓狀態(tài)保證剖分結(jié)果的穩(wěn)定性。這個(gè)技巧在有限元網(wǎng)格生成里也常用不只是為了“好看”。5.3 凸包邊界上的Voronoi射線(xiàn)怎么處理理論上凸包邊界邊的Voronoi邊是一條射線(xiàn)延伸到無(wú)窮遠(yuǎn)。但在屏幕上畫(huà)圖或者做區(qū)域統(tǒng)計(jì)時(shí)不能真的畫(huà)到無(wú)窮遠(yuǎn)需要做一個(gè)裁剪。我的做法是直接取所有Voronoi胞元頂點(diǎn)如果某個(gè)頂點(diǎn)坐標(biāo)超過(guò)畫(huà)布范圍的1.5倍就把它拉回到畫(huà)布邊緣然后再連線(xiàn)。這樣做雖然犧牲了一點(diǎn)點(diǎn)幾何精確度但對(duì)視覺(jué)呈現(xiàn)和后續(xù)的區(qū)域面積統(tǒng)計(jì)影響很小。如果你需要精確的無(wú)限區(qū)域Voronoi建議用polyshape和intersect做裁剪或者直接參考CGAL等成熟庫(kù)的思路。5.4 算法速度慢當(dāng)點(diǎn)集規(guī)模達(dá)到5000以上時(shí)怎么辦逐點(diǎn)插入的Bowyer-Watson算法如果不加任何優(yōu)化復(fù)雜度在點(diǎn)集規(guī)模增大后會(huì)明顯變差。我在測(cè)試2000個(gè)點(diǎn)的時(shí)候還能接受到5000個(gè)點(diǎn)時(shí)每次插入都要掃描全部現(xiàn)有三角形耗時(shí)開(kāi)始變得明顯。如果遇到大數(shù)據(jù)集我的建議不是在MATLAB里硬堆代碼優(yōu)化而是換一個(gè)思路先用delaunayTriangulation內(nèi)置函數(shù)跑通業(yè)務(wù)流程再用自己寫(xiě)的版本做小規(guī)模驗(yàn)證和教學(xué)?；蛘邔?duì)算法加入空間索引比如把平面分成格網(wǎng)插入點(diǎn)P時(shí)只檢查P所在格網(wǎng)及其鄰接格網(wǎng)內(nèi)的三角形這樣能把掃描范圍縮小一個(gè)量級(jí)。這兩種方法我都試過(guò)后者的編碼量確實(shí)比較大所以除非是性能敏感場(chǎng)景否則沒(méi)必要一開(kāi)始就做。一路實(shí)現(xiàn)下來(lái)我最大的感受是Bowyer-Watson本質(zhì)上不復(fù)雜真正考驗(yàn)人的地方全在邊界處理和數(shù)據(jù)結(jié)構(gòu)的組織上。超三角形大小、容差設(shè)置、邊界邊提取方式、Voronoi頂點(diǎn)裁剪每一個(gè)細(xì)節(jié)都可能讓最終結(jié)果從“看起來(lái)合理”變成“嚴(yán)格正確”。如果你打算在MATLAB里復(fù)現(xiàn)這套流程建議把本文第三章的核心代碼完整跑一遍再結(jié)合第五章的踩坑清單做排查。等你自己親手調(diào)通一遍這兩個(gè)計(jì)算幾何里的經(jīng)典結(jié)構(gòu)就會(huì)真正變成你工具箱里的東西而不是只會(huì)調(diào)包的黑盒。本文還有配套的精品資源點(diǎn)擊獲取