面向并行的動態增量式Delaunay 三角剖分算法*

楊昊禹，劉 利，張 誠，于 灝

清華大學 地球系統科學系 地球系統數值模擬教育部重點實驗室，北京100084

1 引言

三角剖分[1]是計算幾何學領域中基礎而又重要的研究內容，其可以將平面或球面等區域中的散點轉化為以這些散點為頂點的三角形網格。三角剖分技術可應用于眾多領域，例如逆向工程、計算機可視化、地理信息系統、有限元分析、地球系統模式等。作為一個基礎算法，三角剖分的計算效率可直接影響到上層應用的整體效率，如何快速高效地完成三角剖分一直是業界重要的討論話題[2-4]。

三角剖分算法相關研究已經有較久遠的歷史。如今最常討論的Delaunay 三角剖分由Boris Delaunay最先提出。Delaunay 三角剖分要求每個三角形均滿足空圓特性，即任意一個三角形的外接圓內均不包含任何其他點。

目前二維平面Delaunay 三角剖分算法（簡稱為三角剖分算法）可分為多種不同類型，包括插入法、構造法、分治法、高維鑲嵌法等，其中最常用的算法為插入法和分治法。插入法邏輯簡單、易于實現，分治法通常擁有最低的時間復雜度，可達到O(nlgn)，其中n表示點數。這兩種算法討論詳見第2 章。

地球系統模式是氣候演變規律研究、未來氣候預測和無縫隙數值預報所不可或缺的科學工具，由分別模擬大氣、陸面、海洋和海冰等地球系統圈層的分量模式通過耦合器耦合而成[5-7]。由于不同分量模式所使用的網格可能不同，因此需要使用插值技術來實現耦合變量在不同網格間的插值。插值技術的實現依賴于臨近點搜索技術，而三角剖分的實現正好可以給臨近點搜索提供具有通用性的基礎支持。

雖然三角剖分算法已具有較低時間復雜度，但隨著網格分辨率的不斷提高，三角剖分的開銷仍會越來越大，這樣的開銷仍不可忽略。為了尋求進一步的加速，大家開始關注并行三角剖分的研究。

并行三角剖分算法按照并行分塊策略的不同，可以分為非重疊分塊法和重疊式分塊法兩類。非重疊分塊法將全球分為若干互不重疊的子塊，并行地對各塊進行局地三角剖分，最后通過串行縫合操作將各子塊合并為最終結果。重疊式分塊法會在分塊時保持各子塊間的部分重疊，避免了合并時的串行縫合操作。雖然重疊會引入少量額外計算，但因其可以避免串行縫合，因此具有更高的并行可擴展性。在并行規模越來越大的今天，重疊式分塊法比非重疊分塊法有更大優勢。

重疊式分塊法現在也可分為兩個子類，靜態重疊法[8]和動態重疊法[9]。靜態重疊法在分塊之時就已經確定下各分塊的大小及互相重疊的程度。動態重疊法則允許算法在執行中動態改變分塊重疊區的大小。因為重疊區域過小會造成子塊間的公共三角形不一致，最終導致并行計算結果有誤（串行與并行計算結果不同），所以靜態重疊法往往使用過飽和的重疊區域以保證結果正確，但這樣會引入過多額外計算開銷，拖慢整體速度。動態重疊法因為支持分塊重疊區的變動，所以可以降低上述額外開銷，并使最終結果的正確性更可靠。

動態重疊法帶來了一個新的問題，對于單一分塊的局地三角剖分，在完成對已有點的三角剖分后，若重疊區擴大，增加了若干新的點，如何求取考慮了新增點的三角剖分結果？一般說來，有兩種實現方式：第一種方式是刪除已有三角剖分結果后，重新求解所有點（包含新增點）的三角剖分結果；第二種方式是在保留已有三角剖分結果基礎上，增量求解新增點對三角剖分結果的調整。毫無疑問，后者將比前者更佳高效，本文將后者稱為增量三角剖分算法。基于動態重疊法設計的并行三角化軟件PatCC1[9]（parallel triangulation algorithm with commonality and parallel consistency,version 1）在實現增量三角剖分時采用了上述第一種方式，若能設計出第二種方式的增量實現，則會給PatCC1 帶來更高的性能提升。

當前已有的三角剖分算法都無法完全滿足增量三角剖分需求，為此本文提出了TID（truly incremental Delaunay）算法，該算法可以高效實現增量三角剖分。

2 相關工作

2.1 插入法

插入法是最為常用的三角剖分算法。該類算法首先創建一個假想的極大三角形以包圍住所有點，隨后逐步插入這些點，并對每個點進行如下操作：定位該點所在的三角形，將該三角形替換為3 個或兩個子三角形（后稱分裂操作），然后判斷相應邊的Delaunay 合法性并進行適當的翻轉操作[10]（又稱局部優化）。待所有點插入完成，最終三角剖分結果也隨之完成。

插入法的關鍵步驟即為點的定位,即判斷一個點被當前哪個三角形所包含。現有定位方式可分為以下四類：Walking 方法[11]、Delaunay 樹方法[12]、索引矩陣法[13]以及分解定位法[14-15]。本文稱使用前三類定位方式的插入法為在線插入法，使用分解定位法的插入法為離線插入法。

Walking 方法的核心操作即為在當前邊的臨近邊中選取離插入點最近的邊作為新的當前邊，該方法從隨機某邊開始，重復執行上述操作來逐步逼近插入點所在的三角形。Delaunay 樹方法，會在插點過程中以三角形為節點建立一個搜索樹，點的定位即可通過遍歷搜索樹來完成。索引矩陣法，將整個區域分為若干矩形子域，用以粗略記錄三角形的位置信息。該方法中點的定位均通過初步定位、局地精確遍歷兩個步驟來完成。在分解定位法中，一個三角形除包含三個頂點外，還包含若干位于該三角形中的尚未插入的點（稱為待插點）。在生成若干新三角形后，待插點會被分配到這些三角形中。

Walking 方法中臨邊快速選取的實現方法及數據結構較為復雜；Delaunay 方法的缺點是內存占用大，維護三角形間的樹形關系較為困難；索引矩陣法的缺點則是索引效率受總點數及矩陣大小的影響，輸入點規模的變化給定位效率帶來的波動較大。此外，以上三種定位方法的逐點定位方法無法高效利用中央處理器（central processing unit,CPU）高速緩存（cache）帶來的加速效果。分解定位法對點進行統一管理，更易通過連續訪存實現定位，具有較高的計算效率，但因其為離線算法，必須事先準備好所有點。

2.2 分治法

分治法三角剖分主要由劃分與合并兩部分操作組成。劃分即為按照一定的方法將平面或球面的一個區域分解為若干子區域，然后分別為各子區域單獨求解三角剖分；合并則是將各子區域的三角剖分結果整合為一個整體，這過程的局部優化操作會給三角剖分的結果帶來調整。常見的劃分方法主要包括一維二分法（沿水平線或垂線二分）和二維二分法（沿水平線和垂線交替二分）。

一方面，當區域中的點很少時，分治法的效率會遠低于插入法；另一方面，分治法合并過程實現起來比較困難，特別是當需要保證Delaunay 三角剖分的空圓特性時。

3 增量算法設計

增量三角剖分的需求來源于動態重疊式并行三角剖分算法，該類并行算法分塊范圍的動態變化需要在已有三角剖分的基礎上加入新增點，并求解更新后的三角剖分。可以總結出該增量三角剖分場景具有以下特點：

（1）新增點大多位于已有三角形外，少量位于已有三角形內（后文分別稱為外增點和內增點）。

（2）性能需求：新增點并非逐一加入，而是批量式加入。

（3）現有的動態重疊法并行算法PatCC1 使用了分解定位插入法來完成局地三角化。

對于增量三角化的實現，在線插入法可以支持內增點，但無法高效支持外增點。雖然可以采取折中策略避免出現外增點，即設置一個非常大的初始三角形，但卻會帶來較大性能損失。例如，極大初始三角形會導致索引矩陣覆蓋范圍很大，但常用的卻只為中間少量索引塊；極大初始三角形也會使Delaunay 樹不平衡，增加樹深，增加單次搜索的耗時。

離線插入法雖然支持內增點，但每次增加點均需要通過全局搜索對點進行定位，效率低下，此外離線插入法也無法支持外增點。

分治法既不支持內增點，也不支持外增點。這是因為分治法不是逐點插入型算法，同時分治法的合并操作缺乏通用性，無法實現對已有區塊和新增區塊在任意位置下的合并。

綜上所述，無法基于分治法而只能基于插入法來實現增量三角剖分。由于分解定位插入法具有很好的數據訪問局部性，且PatCC1 也采用該算法來實現局地三角剖分，因此基于分解定位插入法來實現增量三角剖分，既有利于取得高效率，也有利于重用PatCC1 中分解定位插入法的實現，但必須解決以下兩個問題：

（1）外增點的支持；

（2）高效實現所有新增點到三角形的定位。

為此，本文引入原有剖分的擴展操作來將所有外增點均轉化為內增點，并使用索引技術實現新增點的高效定位，設計出了增量三角剖分算法TID。

3.1 相關名詞

假想矩形：算法為了便于執行插入法三角剖分所引入包含所有點的矩形。

假想點：假想矩形的4 個頂點。

原有點：在增加新點之前已有的點。

新增點：需要被插入到原有剖分中的點。

原有剖分：在增加新點之前已有的三角剖分。

外擴剖分：原有剖分經過擴展后得到的新剖分。

外層三角形：頂點中包含假想點的三角形。

三角形的外邊框：能包圍住該三角形，且各邊均垂直或平行于坐標軸橫軸的最小矩形。

3.2 整體設計

增量算法整體流程如圖1，主要由以下步驟組成：

（1）初始狀態的設定，主要涉及假想點的生成。

（2）判斷是否存在新增點，若存在，則執行下一步，否則算法結束。

（3）根據新增點對原有剖分進行適當外擴。

（4）定位新增點至原有剖分中。

（5）基于原有三角剖分及新增點，動態更新三角剖分，返回（2）。

Fig.1 Main flowchart of TID圖1 TID 算法流程圖

該算法可以滿足前述兩點需求，其中（3）實現了外增點的支持，（4）可以完成新增點的高效定位。

3.3 算法實現

本節將詳細討論TID 算法各步驟的具體實現以及算法對四點共圓情況的特殊處理。

3.3.1 初始化階段

與其他插入法類似，TID 算法使用假想矩形包圍住所有點。在初始化階段，算法根據第一輪新增點的位置分布來創建假想矩形，保證假想矩形能夠包圍住所有第一輪新增點。

3.3.2 原有剖分外擴

若原有剖分對應的假想矩形不能包含全部已有新增點，則算法會首先擴大該假想矩形，然后更新所有外層三角形的頂點坐標，從而完成原有剖分外擴。

為確定假想矩形的擴大范圍，算法先掃描新增點并確定坐標的邊界，然后分別以邊界邊長的預定比例對邊界最小值進行減小，對最大值進行增加，進而得到預期外擴范圍。本文將上述預定比例稱為外擴參數，與該參數相關的測試結果及討論見4.4 節。

原有剖分外擴過程如圖2，圖中的5 個實心圓點為原有點，4 個空心圓點為假想點，3 個空心菱形為新增點，圖2（a）和圖2（b）分別代表擴展前后三角剖分的狀態。通過原有剖分的外擴，所有外增點均已被轉化為內增點。在外擴后，外層三角形由于形狀的變化，其Delaunay 合法性可能會受到改變，但因為這些外層三角形僅為假想三角形，所以不會對最終結果產生影響。

Fig.2 Expansion of existing triangulation圖2 原有剖分外擴示意圖

3.3.3 新增點定位

新增點定位即將各新增點根據位置關系定位到外擴剖分的某個三角形中。TID 算法通過索引矩陣法來完成高效定位。

算法首先把當前外擴剖分覆蓋范圍均勻劃分為若干索引塊并建立索引矩陣，再把各新增點放入相應索引塊內（將含有新增點的索引塊標記為非空索引塊）；然后掃描外擴剖分中的所有三角形，若三角形與非空索引塊的區域有重疊，則將其記錄在相應索引塊中；最后依次在各非空索引塊內，確定新增點所在的三角形，從而實現新增點的快速定位。

索引塊的數量是影響新增點定位效率的重要參數。索引塊過少會導致一個索引塊包含有很多新增點和三角形，從而導致索引效率的降低；索引塊過多雖然會降低復雜度，但會導致Cachemiss 增多，無法較好利用CPU 高速緩存。因此，應該在保證索引過程不增加TID 算法復雜度的前提下，選擇最小的索引塊數量。對于不同的數據規模，TID 算法會根據公式動態推測索引塊的最優數量。該公式的推導詳見3.3.6 小節，相關測試結果詳見4.3 節。當已有三角形數量很少時，使用索引定位法的意義不大，此時將直接使用全局搜索。

為了加快點與三角形位置關系的判斷，TID 算法會先使用三角形的外邊框進行預判，只有當點位于外邊框內時，才會精確判斷該點是否位于三角形內。

3.3.4 三角剖分動態更新

在經過新增點定位后，所有新增點均已被成組地分配到若干三角形中，稱這些三角形為非空三角形。TID 算法通過單獨對每一個非空三角形執行分裂操作來完成整體三角剖分的更新。分裂操作由以下三步完成：

（1）子三角形的生成；

（2）局部優化；

（3）父三角形待插點到子三角形的分配。

首先，TID 算法從父三角形待插點中取出離父三角形重心最近的點，以該點作為插入點，連接該點與父三角形各頂點形成多個子三角形。如圖3 所示，當插入點位于父三角形內時，會形成3 個子三角形；當插入點位于父三角形某一邊上時，則在形成兩個子三角形同時，將共享該邊的另一個三角形也分裂為兩個子三角形，共生成4 個子三角形。

其次，對所有新生成的子三角形進行局部優化。即判斷所有子三角形邊的Delaunay 合法性，若不合法則通過翻轉操作將其轉化為合法邊。翻轉操作會引入新的父三角形和子三角形。

Fig.3 Generation of new triangles圖3 子三角形生成示例

最后，將父三角形剩余待插點按照位置關系分配至相應子三角形內。需注意由于翻轉操作的存在，父三角形可能存在多個，子三角形個數也可能大于4 個。

TID 算法會遞歸地分裂子三角形，當所有子三角形都不包含任何點時，本次三角剖分完成。

為減少內存使用，避免頻繁的內存分配帶來的時間開銷，本文采用數組加鏈表的方式來存儲待插點。該方法只使用O(n)的空間，同時其順序訪存可以有效利用起CPU 高速緩存。

3.3.5 四點共圓特例

嚴格的Delaunay 三角剖分不允許存在4 點共圓的情況，但共圓情況在地球系統模式網格中卻無法避免，甚至是經常出現的。例如，在地球系統模式中最常用的經緯網格中，經度或緯度值相鄰的4 個格點均共圓。該情況下三角剖分會存在多種合法解。這會給算法行為帶來不確定性，若用于并行計算，則會給其帶來串并行不一致的風險。本文針對這種情況，使用了與PatCC1[9]相同的特殊判斷規則：

四點共圓合法性給定兩個相鄰三角形及其二者的頂點，若該4 個頂點共圓，則當且僅當該4 點中的最左點（若存在多個最左點則為最左下點）不被這兩個三角形共享時，這兩個三角形是合法三角形。

3.3.6 索引矩陣參數公式

給定n個原有點和m個新增點，對于TID 算法，三角剖分動態更新階段的時間復雜度為：

給定k×k的索引矩陣，本文稱k為索引矩陣參數。構建該索引矩陣的時間復雜度為O(m+n)，即O(n)。假設新增點在外圈均勻分布，則非空索引塊位于索引矩陣的外層，非空索引塊數為O(k)，各索引塊平均包含三角形數量為，各非空索引塊平均包含新增點數為，因此索引矩陣法的定位復雜度為上述三者之積，即。

為保證效率，本文要求索引定位法的時間復雜度既不高于構建索引矩陣的時間復雜度，也不高于三角剖分動態更新的時間復雜度，即：

本文選取滿足上述條件的最小k值作為最優索引矩陣參數K，可以表示為：

其中，a與b為系數，可通過實際測試得來。

4 性能評估

本文的性能評估主要包括：核心三角剖分性能評估和增量場景下的性能評估。對于前者，本文以Github 上的三角剖分軟件作為對照組，在非增量場景進行了TID 算法和同類算法的性能對比。對于后者，本文首先對比了增量場景下TID 算法和PatCC1 原始算法的性能差異，然后通過比較TID 算法本身在增量場景和非增量場景的性能差異，來探究增量功能帶來的額外開銷。此外，本文還分別針對索引矩陣參數和原有剖分外擴參數進行了相應的實驗探究。

4.1 實驗配置

本文全部實驗均在個人電腦上完成，其裝有Intel Core i9-9900K 8 核3.6 GHz 處理器，32 GB 內存，可執行程序均使用GNU 6.3.0 編譯器在O3 優化等級下進行編譯。

本文的實驗網格采用了3 種不同密度，粗、中、細，以及3 種不同類型，隨機生成網格、均勻網格、真實網格。其中真實網格選取了地球系統模式中的球面立方網格作為實驗數據。

實驗共使用了3 種不同插點方法：

（1）非增量式：一次性加入所有點并生成三角剖分。

（2）單次增量式：將點分為兩組，依次加入并生成三角剖分。

（3）多次增量式：將點分為多組，依次加入并生成三角剖分。

4.2 核心三角剖分性能

核心三角剖分即TID 算法整體流程的第5 步。本實驗的對照組為Juannavascalvente 的三角剖分軟件（https://github.com/juannavascalvente/Delaunay，后稱Juan 軟件），其使用基于Delaunay 樹定位的插入式算法。本實驗使用非增量式插點方法，分別統計了兩者在不同類型、不同分辨率的網格下的三角剖分耗時。該耗時只包括核心三角剖分過程，去除了文件輸入輸出等無關過程的影響。

實驗結果如表1，表中Timeout 代表實際執行時間超過1 h，OOM（out of memory）代表程序執行時耗盡系統內存。可以看出TID 算法在每個網格下的計算效率均優于Juan 軟件，同時網格規模越大，TID 算法的優勢越明顯。因此，TID算法所采用的分解定位法在運行速度和內存占用上均優于Delaunay樹定位法。

Table 1 Comparison of kernel triangulation time between Juan algorithm and TID algorithm表1 Juan 和TID 算法核心三角剖分耗時

4.3 索引矩陣參數實驗

索引矩陣參數決定著新增點定位的效率。為確定3.3.6 小節中最優索引矩陣參數公式的系數a、b，本節使用單次增量插點法基于多種分辨率的隨機生成網格進行了探究。

本文分別在粗網格、中網格、細網格下使用不同的索引矩陣參數測試了新增點定位的整體耗時。測試結果如圖4 所示，其橫坐標為索引矩陣參數，縱坐標為新增點定位耗時進行歸一化后的結果，即實際耗時除以該網格下的最低耗時。

Fig.4 Tendency of locating time with indexing map parameter using various grids圖4 多種網格下定位耗時隨索引矩陣參數變化趨勢

從圖4 中數據可以看出，隨著索引矩陣參數逐漸增大，各網格新增點定位耗時均呈現先減小，后增大的趨勢。這主要是因為索引塊較少時，每個索引塊內點數及三角形數量過多，塊內定位的開銷仍比較大，隨著索引矩陣參數增大，這樣的開銷越來越小。但當索引矩陣參數過大時，過于密集的索引矩陣導致單索引塊中所含點數及三角形數量很少，在進行塊內定位時會產生較多Cachemiss 現象，反而拖慢定位速度。從圖中可以看出，粗、中、細網格的最優矩陣參數分別為30、70、100。本文規定定位耗時不超過最低耗時0.1 倍的部分為最優區間，則上述網格的最優區間分別為[21,37]、[40,140]、[70,300]。

綜合粗、中、細網格的原有點數、新增點數及索引矩陣參數最優值，本文嘗試不同的系數配置，并最終確定了合適的系數配置，a為0.2，b為0.3。該系數配置可以保證3 個網格下的公式推測值均落入最優區間中，如表2 所示。為進一步驗證公式準確性，本文構造了點數更多的隨機生成網格作為驗證網格，并測試了其不同索引矩陣參數下的定位耗時，結果如圖5。該網格的公式推測值也在表2 中一并列出，可見推測值仍位于圖中的最優區間內。

4.4 原有剖分外擴參數實驗

本節探究了3.3.2 小節中外擴參數對新增點定位耗時和三角剖分動態更新耗時的影響。該實驗使用單次增量插點法，基于細密度的隨機生成網格，索引矩陣參數采用4.3 節中測試出的最優值100。實驗分別統計了多種外擴參數下兩個主要階段的耗時，并記錄了非空索引塊數量占總索引塊數量的比重。

Table 2 Prediction of optimal indexing map parameters using grids with different resolution表2 不同分辨率網格下最優索引矩陣參數的推測

Fig.5 Tendency of locating time with indexing map parameter using very fine grid圖5 極細網格下定位耗時隨索引矩陣參數變化趨勢

測試結果見表3，隨著外擴系數的增大，新增點定位耗時的增加非常明顯，這是因為外擴系數增大會導致索引矩陣范圍的增大，進而導致新增點位于越來越少的索引矩陣塊中，索引效率降低，這可從非空索引塊占比的逐漸降低中看出。另外，三角剖分動態更新耗時隨著外擴系數的增加而輕微增加，這主要是因為剖分范圍的增大導致外層三角形被“拉長”，新增點在單個三角形中的分布不均勻，進而導致三角形的分裂次數增多，耗時增加。由該實驗可知，為保證程序運行的高效性，應使用盡量小的外擴參數。

Table 3 Running time of two main steps using different expansion parameters表3 不同外擴參數下的兩個主要階段耗時

4.5 增量場景整體性能

增量場景的性能實驗由兩部分組成：TID 算法與PatCC1 原始算法對比實驗、增量開銷實驗。這兩個性能實驗均將輸入點分為10 組，每組的點數均已在結果表格中給出。

PatCC1 在更新三角剖分時采用了刪除原有剖分并重新求解的方法，本文以該方法作為對照組，使用多次增量式插點方法進行了性能對比實驗。結果如表4，TID 算法相比PatCC1 原始算法具有很明顯的速度優勢。可見若將TID 算法用于PatCC1 中，將會給其帶來較大的性能提升。

Table 4 Comparison of triangulation time between TID algorithm and PatCC1 algorithm表4 TID 算法與PatCC1 原始算法增量三角剖分耗時

增量插點功能的引入也會帶來額外開銷，其主要來源于索引矩陣的構建及新增點的定位。為探究該開銷大小，本文分別統計了TID 算法對相同網格采用不同插點方式（非增量式和多次增量式）時的耗時，實驗結果見表5。表中序號n對應的增量單次耗時即為增量插入第n組新點并更新三角剖分所消耗的時間，非增量耗時即為將第0 組到第n組的所有點一次性插入并進行三角剖分的耗時。假設當前增量耗時為Tinc，上一次非增量耗時T0,當前非增量耗時Tc，則額外開銷占比。從結果中可以看出，增量功能并不會引入過多額外開銷，其在較好情況下與非增量式的耗時基本相近，在較壞情況下額外開銷仍低于10%。

Table 5 Comparison of triangulation time between incremental TID algorithm and non-incremental TID algorithm表5 TID 算法增量式與非增量式三角剖分耗時

5 總結

本文基于分解定位插入法設計了一種高效的增量三角剖分算法TID，該算法支持在已有三角剖分基礎上增加任意位置的點并更新三角剖分。該算法不僅具有增量插點特性，其與同類非增量軟件相比也具有更高的計算效率。此外，合法化特殊規則的引入確保該算法能夠生成唯一的三角剖分結果。

TID算法現已成功用于并行三角剖分軟件PatCC1。TID 算法源代碼見https://github.com/WireFisher/TID。