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