面向并行的動(dòng)態(tài)增量式Delaunay 三角剖分算法*

楊昊禹，劉 利，張 誠(chéng)，于 灝

清華大學(xué) 地球系統(tǒng)科學(xué)系 地球系統(tǒng)數(shù)值模擬教育部重點(diǎn)實(shí)驗(yàn)室，北京100084

1 引言

三角剖分[1]是計(jì)算幾何學(xué)領(lǐng)域中基礎(chǔ)而又重要的研究?jī)?nèi)容，其可以將平面或球面等區(qū)域中的散點(diǎn)轉(zhuǎn)化為以這些散點(diǎn)為頂點(diǎn)的三角形網(wǎng)格。三角剖分技術(shù)可應(yīng)用于眾多領(lǐng)域，例如逆向工程、計(jì)算機(jī)可視化、地理信息系統(tǒng)、有限元分析、地球系統(tǒng)模式等。作為一個(gè)基礎(chǔ)算法，三角剖分的計(jì)算效率可直接影響到上層應(yīng)用的整體效率，如何快速高效地完成三角剖分一直是業(yè)界重要的討論話(huà)題[2-4]。

三角剖分算法相關(guān)研究已經(jīng)有較久遠(yuǎn)的歷史。如今最常討論的Delaunay 三角剖分由Boris Delaunay最先提出。Delaunay 三角剖分要求每個(gè)三角形均滿(mǎn)足空?qǐng)A特性，即任意一個(gè)三角形的外接圓內(nèi)均不包含任何其他點(diǎn)。

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

地球系統(tǒng)模式是氣候演變規(guī)律研究、未來(lái)氣候預(yù)測(cè)和無(wú)縫隙數(shù)值預(yù)報(bào)所不可或缺的科學(xué)工具，由分別模擬大氣、陸面、海洋和海冰等地球系統(tǒng)圈層的分量模式通過(guò)耦合器耦合而成[5-7]。由于不同分量模式所使用的網(wǎng)格可能不同，因此需要使用插值技術(shù)來(lái)實(shí)現(xiàn)耦合變量在不同網(wǎng)格間的插值。插值技術(shù)的實(shí)現(xiàn)依賴(lài)于臨近點(diǎn)搜索技術(shù)，而三角剖分的實(shí)現(xiàn)正好可以給臨近點(diǎn)搜索提供具有通用性的基礎(chǔ)支持。

雖然三角剖分算法已具有較低時(shí)間復(fù)雜度，但隨著網(wǎng)格分辨率的不斷提高，三角剖分的開(kāi)銷(xiāo)仍會(huì)越來(lái)越大，這樣的開(kāi)銷(xiāo)仍不可忽略。為了尋求進(jìn)一步的加速，大家開(kāi)始關(guān)注并行三角剖分的研究。

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

重疊式分塊法現(xiàn)在也可分為兩個(gè)子類(lèi)，靜態(tài)重疊法[8]和動(dòng)態(tài)重疊法[9]。靜態(tài)重疊法在分塊之時(shí)就已經(jīng)確定下各分塊的大小及互相重疊的程度。動(dòng)態(tài)重疊法則允許算法在執(zhí)行中動(dòng)態(tài)改變分塊重疊區(qū)的大小。因?yàn)橹丿B區(qū)域過(guò)小會(huì)造成子塊間的公共三角形不一致，最終導(dǎo)致并行計(jì)算結(jié)果有誤（串行與并行計(jì)算結(jié)果不同），所以靜態(tài)重疊法往往使用過(guò)飽和的重疊區(qū)域以保證結(jié)果正確，但這樣會(huì)引入過(guò)多額外計(jì)算開(kāi)銷(xiāo)，拖慢整體速度。動(dòng)態(tài)重疊法因?yàn)橹С址謮K重疊區(qū)的變動(dòng)，所以可以降低上述額外開(kāi)銷(xiāo)，并使最終結(jié)果的正確性更可靠。

動(dòng)態(tài)重疊法帶來(lái)了一個(gè)新的問(wèn)題，對(duì)于單一分塊的局地三角剖分，在完成對(duì)已有點(diǎn)的三角剖分后，若重疊區(qū)擴(kuò)大，增加了若干新的點(diǎn)，如何求取考慮了新增點(diǎn)的三角剖分結(jié)果？一般說(shuō)來(lái)，有兩種實(shí)現(xiàn)方式：第一種方式是刪除已有三角剖分結(jié)果后，重新求解所有點(diǎn)（包含新增點(diǎn)）的三角剖分結(jié)果；第二種方式是在保留已有三角剖分結(jié)果基礎(chǔ)上，增量求解新增點(diǎn)對(duì)三角剖分結(jié)果的調(diào)整。毫無(wú)疑問(wèn)，后者將比前者更佳高效，本文將后者稱(chēng)為增量三角剖分算法?；趧?dòng)態(tài)重疊法設(shè)計(jì)的并行三角化軟件PatCC1[9]（parallel triangulation algorithm with commonality and parallel consistency,version 1）在實(shí)現(xiàn)增量三角剖分時(shí)采用了上述第一種方式，若能設(shè)計(jì)出第二種方式的增量實(shí)現(xiàn)，則會(huì)給PatCC1 帶來(lái)更高的性能提升。

當(dāng)前已有的三角剖分算法都無(wú)法完全滿(mǎn)足增量三角剖分需求，為此本文提出了TID（truly incremental Delaunay）算法，該算法可以高效實(shí)現(xiàn)增量三角剖分。

2 相關(guān)工作

2.1 插入法

插入法是最為常用的三角剖分算法。該類(lèi)算法首先創(chuàng)建一個(gè)假想的極大三角形以包圍住所有點(diǎn)，隨后逐步插入這些點(diǎn)，并對(duì)每個(gè)點(diǎn)進(jìn)行如下操作：定位該點(diǎn)所在的三角形，將該三角形替換為3 個(gè)或兩個(gè)子三角形（后稱(chēng)分裂操作），然后判斷相應(yīng)邊的Delaunay 合法性并進(jìn)行適當(dāng)?shù)姆D(zhuǎn)操作[10]（又稱(chēng)局部?jī)?yōu)化）。待所有點(diǎn)插入完成，最終三角剖分結(jié)果也隨之完成。

插入法的關(guān)鍵步驟即為點(diǎn)的定位,即判斷一個(gè)點(diǎn)被當(dāng)前哪個(gè)三角形所包含?，F(xiàn)有定位方式可分為以下四類(lèi)：Walking 方法[11]、Delaunay 樹(shù)方法[12]、索引矩陣法[13]以及分解定位法[14-15]。本文稱(chēng)使用前三類(lèi)定位方式的插入法為在線(xiàn)插入法，使用分解定位法的插入法為離線(xiàn)插入法。

Walking 方法的核心操作即為在當(dāng)前邊的臨近邊中選取離插入點(diǎn)最近的邊作為新的當(dāng)前邊，該方法從隨機(jī)某邊開(kāi)始，重復(fù)執(zhí)行上述操作來(lái)逐步逼近插入點(diǎn)所在的三角形。Delaunay 樹(shù)方法，會(huì)在插點(diǎn)過(guò)程中以三角形為節(jié)點(diǎn)建立一個(gè)搜索樹(shù)，點(diǎn)的定位即可通過(guò)遍歷搜索樹(shù)來(lái)完成。索引矩陣法，將整個(gè)區(qū)域分為若干矩形子域，用以粗略記錄三角形的位置信息。該方法中點(diǎn)的定位均通過(guò)初步定位、局地精確遍歷兩個(gè)步驟來(lái)完成。在分解定位法中，一個(gè)三角形除包含三個(gè)頂點(diǎn)外，還包含若干位于該三角形中的尚未插入的點(diǎn)（稱(chēng)為待插點(diǎn)）。在生成若干新三角形后，待插點(diǎn)會(huì)被分配到這些三角形中。

Walking 方法中臨邊快速選取的實(shí)現(xiàn)方法及數(shù)據(jù)結(jié)構(gòu)較為復(fù)雜；Delaunay 方法的缺點(diǎn)是內(nèi)存占用大，維護(hù)三角形間的樹(shù)形關(guān)系較為困難；索引矩陣法的缺點(diǎn)則是索引效率受總點(diǎn)數(shù)及矩陣大小的影響，輸入點(diǎn)規(guī)模的變化給定位效率帶來(lái)的波動(dòng)較大。此外，以上三種定位方法的逐點(diǎn)定位方法無(wú)法高效利用中央處理器（central processing unit,CPU）高速緩存（cache）帶來(lái)的加速效果。分解定位法對(duì)點(diǎn)進(jìn)行統(tǒng)一管理，更易通過(guò)連續(xù)訪存實(shí)現(xiàn)定位，具有較高的計(jì)算效率，但因其為離線(xiàn)算法，必須事先準(zhǔn)備好所有點(diǎn)。

2.2 分治法

分治法三角剖分主要由劃分與合并兩部分操作組成。劃分即為按照一定的方法將平面或球面的一個(gè)區(qū)域分解為若干子區(qū)域，然后分別為各子區(qū)域單獨(dú)求解三角剖分；合并則是將各子區(qū)域的三角剖分結(jié)果整合為一個(gè)整體，這過(guò)程的局部?jī)?yōu)化操作會(huì)給三角剖分的結(jié)果帶來(lái)調(diào)整。常見(jiàn)的劃分方法主要包括一維二分法（沿水平線(xiàn)或垂線(xiàn)二分）和二維二分法（沿水平線(xiàn)和垂線(xiàn)交替二分）。

一方面，當(dāng)區(qū)域中的點(diǎn)很少時(shí)，分治法的效率會(huì)遠(yuǎn)低于插入法；另一方面，分治法合并過(guò)程實(shí)現(xiàn)起來(lái)比較困難，特別是當(dāng)需要保證Delaunay 三角剖分的空?qǐng)A特性時(shí)。

3 增量算法設(shè)計(jì)

增量三角剖分的需求來(lái)源于動(dòng)態(tài)重疊式并行三角剖分算法，該類(lèi)并行算法分塊范圍的動(dòng)態(tài)變化需要在已有三角剖分的基礎(chǔ)上加入新增點(diǎn)，并求解更新后的三角剖分?？梢钥偨Y(jié)出該增量三角剖分場(chǎng)景具有以下特點(diǎn)：

（1）新增點(diǎn)大多位于已有三角形外，少量位于已有三角形內(nèi)（后文分別稱(chēng)為外增點(diǎn)和內(nèi)增點(diǎn)）。

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

（3）現(xiàn)有的動(dòng)態(tài)重疊法并行算法PatCC1 使用了分解定位插入法來(lái)完成局地三角化。

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

離線(xiàn)插入法雖然支持內(nèi)增點(diǎn)，但每次增加點(diǎn)均需要通過(guò)全局搜索對(duì)點(diǎn)進(jìn)行定位，效率低下，此外離線(xiàn)插入法也無(wú)法支持外增點(diǎn)。

分治法既不支持內(nèi)增點(diǎn)，也不支持外增點(diǎn)。這是因?yàn)榉种畏ú皇侵瘘c(diǎn)插入型算法，同時(shí)分治法的合并操作缺乏通用性，無(wú)法實(shí)現(xiàn)對(duì)已有區(qū)塊和新增區(qū)塊在任意位置下的合并。

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

（1）外增點(diǎn)的支持；

（2）高效實(shí)現(xiàn)所有新增點(diǎn)到三角形的定位。

為此，本文引入原有剖分的擴(kuò)展操作來(lái)將所有外增點(diǎn)均轉(zhuǎn)化為內(nèi)增點(diǎn)，并使用索引技術(shù)實(shí)現(xiàn)新增點(diǎn)的高效定位，設(shè)計(jì)出了增量三角剖分算法TID。

3.1 相關(guān)名詞

假想矩形：算法為了便于執(zhí)行插入法三角剖分所引入包含所有點(diǎn)的矩形。

假想點(diǎn)：假想矩形的4 個(gè)頂點(diǎn)。

原有點(diǎn)：在增加新點(diǎn)之前已有的點(diǎn)。

新增點(diǎn)：需要被插入到原有剖分中的點(diǎn)。

原有剖分：在增加新點(diǎn)之前已有的三角剖分。

外擴(kuò)剖分：原有剖分經(jīng)過(guò)擴(kuò)展后得到的新剖分。

外層三角形：頂點(diǎn)中包含假想點(diǎn)的三角形。

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

3.2 整體設(shè)計(jì)

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

（1）初始狀態(tài)的設(shè)定，主要涉及假想點(diǎn)的生成。

（2）判斷是否存在新增點(diǎn)，若存在，則執(zhí)行下一步，否則算法結(jié)束。

（3）根據(jù)新增點(diǎn)對(duì)原有剖分進(jìn)行適當(dāng)外擴(kuò)。

（4）定位新增點(diǎn)至原有剖分中。

（5）基于原有三角剖分及新增點(diǎn)，動(dòng)態(tài)更新三角剖分，返回（2）。

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

該算法可以滿(mǎn)足前述兩點(diǎn)需求，其中（3）實(shí)現(xiàn)了外增點(diǎn)的支持，（4）可以完成新增點(diǎn)的高效定位。

3.3 算法實(shí)現(xiàn)

本節(jié)將詳細(xì)討論TID 算法各步驟的具體實(shí)現(xiàn)以及算法對(duì)四點(diǎn)共圓情況的特殊處理。

3.3.1 初始化階段

與其他插入法類(lèi)似，TID 算法使用假想矩形包圍住所有點(diǎn)。在初始化階段，算法根據(jù)第一輪新增點(diǎn)的位置分布來(lái)創(chuàng)建假想矩形，保證假想矩形能夠包圍住所有第一輪新增點(diǎn)。

3.3.2 原有剖分外擴(kuò)

若原有剖分對(duì)應(yīng)的假想矩形不能包含全部已有新增點(diǎn)，則算法會(huì)首先擴(kuò)大該假想矩形，然后更新所有外層三角形的頂點(diǎn)坐標(biāo)，從而完成原有剖分外擴(kuò)。

為確定假想矩形的擴(kuò)大范圍，算法先掃描新增點(diǎn)并確定坐標(biāo)的邊界，然后分別以邊界邊長(zhǎng)的預(yù)定比例對(duì)邊界最小值進(jìn)行減小，對(duì)最大值進(jìn)行增加，進(jìn)而得到預(yù)期外擴(kuò)范圍。本文將上述預(yù)定比例稱(chēng)為外擴(kuò)參數(shù)，與該參數(shù)相關(guān)的測(cè)試結(jié)果及討論見(jiàn)4.4 節(jié)。

原有剖分外擴(kuò)過(guò)程如圖2，圖中的5 個(gè)實(shí)心圓點(diǎn)為原有點(diǎn)，4 個(gè)空心圓點(diǎn)為假想點(diǎn)，3 個(gè)空心菱形為新增點(diǎn)，圖2（a）和圖2（b）分別代表擴(kuò)展前后三角剖分的狀態(tài)。通過(guò)原有剖分的外擴(kuò)，所有外增點(diǎn)均已被轉(zhuǎn)化為內(nèi)增點(diǎn)。在外擴(kuò)后，外層三角形由于形狀的變化，其Delaunay 合法性可能會(huì)受到改變，但因?yàn)檫@些外層三角形僅為假想三角形，所以不會(huì)對(duì)最終結(jié)果產(chǎn)生影響。

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

3.3.3 新增點(diǎn)定位

新增點(diǎn)定位即將各新增點(diǎn)根據(jù)位置關(guān)系定位到外擴(kuò)剖分的某個(gè)三角形中。TID 算法通過(guò)索引矩陣法來(lái)完成高效定位。

算法首先把當(dāng)前外擴(kuò)剖分覆蓋范圍均勻劃分為若干索引塊并建立索引矩陣，再把各新增點(diǎn)放入相應(yīng)索引塊內(nèi)（將含有新增點(diǎn)的索引塊標(biāo)記為非空索引塊）；然后掃描外擴(kuò)剖分中的所有三角形，若三角形與非空索引塊的區(qū)域有重疊，則將其記錄在相應(yīng)索引塊中；最后依次在各非空索引塊內(nèi)，確定新增點(diǎn)所在的三角形，從而實(shí)現(xiàn)新增點(diǎn)的快速定位。

索引塊的數(shù)量是影響新增點(diǎn)定位效率的重要參數(shù)。索引塊過(guò)少會(huì)導(dǎo)致一個(gè)索引塊包含有很多新增點(diǎn)和三角形，從而導(dǎo)致索引效率的降低；索引塊過(guò)多雖然會(huì)降低復(fù)雜度，但會(huì)導(dǎo)致Cachemiss 增多，無(wú)法較好利用CPU 高速緩存。因此，應(yīng)該在保證索引過(guò)程不增加TID 算法復(fù)雜度的前提下，選擇最小的索引塊數(shù)量。對(duì)于不同的數(shù)據(jù)規(guī)模，TID 算法會(huì)根據(jù)公式動(dòng)態(tài)推測(cè)索引塊的最優(yōu)數(shù)量。該公式的推導(dǎo)詳見(jiàn)3.3.6 小節(jié)，相關(guān)測(cè)試結(jié)果詳見(jiàn)4.3 節(jié)。當(dāng)已有三角形數(shù)量很少時(shí)，使用索引定位法的意義不大，此時(shí)將直接使用全局搜索。

為了加快點(diǎn)與三角形位置關(guān)系的判斷，TID 算法會(huì)先使用三角形的外邊框進(jìn)行預(yù)判，只有當(dāng)點(diǎn)位于外邊框內(nèi)時(shí)，才會(huì)精確判斷該點(diǎn)是否位于三角形內(nèi)。

3.3.4 三角剖分動(dòng)態(tài)更新

在經(jīng)過(guò)新增點(diǎn)定位后，所有新增點(diǎn)均已被成組地分配到若干三角形中，稱(chēng)這些三角形為非空三角形。TID 算法通過(guò)單獨(dú)對(duì)每一個(gè)非空三角形執(zhí)行分裂操作來(lái)完成整體三角剖分的更新。分裂操作由以下三步完成：

（1）子三角形的生成；

（2）局部?jī)?yōu)化；

（3）父三角形待插點(diǎn)到子三角形的分配。

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

其次，對(duì)所有新生成的子三角形進(jìn)行局部?jī)?yōu)化。即判斷所有子三角形邊的Delaunay 合法性，若不合法則通過(guò)翻轉(zhuǎn)操作將其轉(zhuǎn)化為合法邊。翻轉(zhuǎn)操作會(huì)引入新的父三角形和子三角形。

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

最后，將父三角形剩余待插點(diǎn)按照位置關(guān)系分配至相應(yīng)子三角形內(nèi)。需注意由于翻轉(zhuǎn)操作的存在，父三角形可能存在多個(gè)，子三角形個(gè)數(shù)也可能大于4 個(gè)。

TID 算法會(huì)遞歸地分裂子三角形，當(dāng)所有子三角形都不包含任何點(diǎn)時(shí)，本次三角剖分完成。

為減少內(nèi)存使用，避免頻繁的內(nèi)存分配帶來(lái)的時(shí)間開(kāi)銷(xiāo)，本文采用數(shù)組加鏈表的方式來(lái)存儲(chǔ)待插點(diǎn)。該方法只使用O(n)的空間，同時(shí)其順序訪存可以有效利用起CPU 高速緩存。

3.3.5 四點(diǎn)共圓特例

嚴(yán)格的Delaunay 三角剖分不允許存在4 點(diǎn)共圓的情況，但共圓情況在地球系統(tǒng)模式網(wǎng)格中卻無(wú)法避免，甚至是經(jīng)常出現(xiàn)的。例如，在地球系統(tǒng)模式中最常用的經(jīng)緯網(wǎng)格中，經(jīng)度或緯度值相鄰的4 個(gè)格點(diǎn)均共圓。該情況下三角剖分會(huì)存在多種合法解。這會(huì)給算法行為帶來(lái)不確定性，若用于并行計(jì)算，則會(huì)給其帶來(lái)串并行不一致的風(fēng)險(xiǎn)。本文針對(duì)這種情況，使用了與PatCC1[9]相同的特殊判斷規(guī)則：

四點(diǎn)共圓合法性給定兩個(gè)相鄰三角形及其二者的頂點(diǎn)，若該4 個(gè)頂點(diǎn)共圓，則當(dāng)且僅當(dāng)該4 點(diǎn)中的最左點(diǎn)（若存在多個(gè)最左點(diǎn)則為最左下點(diǎn)）不被這兩個(gè)三角形共享時(shí)，這兩個(gè)三角形是合法三角形。

3.3.6 索引矩陣參數(shù)公式

給定n個(gè)原有點(diǎn)和m個(gè)新增點(diǎn)，對(duì)于TID 算法，三角剖分動(dòng)態(tài)更新階段的時(shí)間復(fù)雜度為：

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

為保證效率，本文要求索引定位法的時(shí)間復(fù)雜度既不高于構(gòu)建索引矩陣的時(shí)間復(fù)雜度，也不高于三角剖分動(dòng)態(tài)更新的時(shí)間復(fù)雜度，即：

本文選取滿(mǎn)足上述條件的最小k值作為最優(yōu)索引矩陣參數(shù)K，可以表示為：

其中，a與b為系數(shù)，可通過(guò)實(shí)際測(cè)試得來(lái)。

4 性能評(píng)估

本文的性能評(píng)估主要包括：核心三角剖分性能評(píng)估和增量場(chǎng)景下的性能評(píng)估。對(duì)于前者，本文以Github 上的三角剖分軟件作為對(duì)照組，在非增量場(chǎng)景進(jìn)行了TID 算法和同類(lèi)算法的性能對(duì)比。對(duì)于后者，本文首先對(duì)比了增量場(chǎng)景下TID 算法和PatCC1 原始算法的性能差異，然后通過(guò)比較TID 算法本身在增量場(chǎng)景和非增量場(chǎng)景的性能差異，來(lái)探究增量功能帶來(lái)的額外開(kāi)銷(xiāo)。此外，本文還分別針對(duì)索引矩陣參數(shù)和原有剖分外擴(kuò)參數(shù)進(jìn)行了相應(yīng)的實(shí)驗(yàn)探究。

4.1 實(shí)驗(yàn)配置

本文全部實(shí)驗(yàn)均在個(gè)人電腦上完成，其裝有Intel Core i9-9900K 8 核3.6 GHz 處理器，32 GB 內(nèi)存，可執(zhí)行程序均使用GNU 6.3.0 編譯器在O3 優(yōu)化等級(jí)下進(jìn)行編譯。

本文的實(shí)驗(yàn)網(wǎng)格采用了3 種不同密度，粗、中、細(xì)，以及3 種不同類(lèi)型，隨機(jī)生成網(wǎng)格、均勻網(wǎng)格、真實(shí)網(wǎng)格。其中真實(shí)網(wǎng)格選取了地球系統(tǒng)模式中的球面立方網(wǎng)格作為實(shí)驗(yàn)數(shù)據(jù)。

實(shí)驗(yàn)共使用了3 種不同插點(diǎn)方法：

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

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

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

4.2 核心三角剖分性能

核心三角剖分即TID 算法整體流程的第5 步。本實(shí)驗(yàn)的對(duì)照組為Juannavascalvente 的三角剖分軟件（https://github.com/juannavascalvente/Delaunay，后稱(chēng)Juan 軟件），其使用基于Delaunay 樹(shù)定位的插入式算法。本實(shí)驗(yàn)使用非增量式插點(diǎn)方法，分別統(tǒng)計(jì)了兩者在不同類(lèi)型、不同分辨率的網(wǎng)格下的三角剖分耗時(shí)。該耗時(shí)只包括核心三角剖分過(guò)程，去除了文件輸入輸出等無(wú)關(guān)過(guò)程的影響。

實(shí)驗(yàn)結(jié)果如表1，表中Timeout 代表實(shí)際執(zhí)行時(shí)間超過(guò)1 h，OOM（out of memory）代表程序執(zhí)行時(shí)耗盡系統(tǒng)內(nèi)存?？梢钥闯鯰ID 算法在每個(gè)網(wǎng)格下的計(jì)算效率均優(yōu)于Juan 軟件，同時(shí)網(wǎng)格規(guī)模越大，TID 算法的優(yōu)勢(shì)越明顯。因此，TID算法所采用的分解定位法在運(yùn)行速度和內(nèi)存占用上均優(yōu)于Delaunay樹(shù)定位法。

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

4.3 索引矩陣參數(shù)實(shí)驗(yàn)

索引矩陣參數(shù)決定著新增點(diǎn)定位的效率。為確定3.3.6 小節(jié)中最優(yōu)索引矩陣參數(shù)公式的系數(shù)a、b，本節(jié)使用單次增量插點(diǎn)法基于多種分辨率的隨機(jī)生成網(wǎng)格進(jìn)行了探究。

本文分別在粗網(wǎng)格、中網(wǎng)格、細(xì)網(wǎng)格下使用不同的索引矩陣參數(shù)測(cè)試了新增點(diǎn)定位的整體耗時(shí)。測(cè)試結(jié)果如圖4 所示，其橫坐標(biāo)為索引矩陣參數(shù)，縱坐標(biāo)為新增點(diǎn)定位耗時(shí)進(jìn)行歸一化后的結(jié)果，即實(shí)際耗時(shí)除以該網(wǎng)格下的最低耗時(shí)。

Fig.4 Tendency of locating time with indexing map parameter using various grids圖4 多種網(wǎng)格下定位耗時(shí)隨索引矩陣參數(shù)變化趨勢(shì)

從圖4 中數(shù)據(jù)可以看出，隨著索引矩陣參數(shù)逐漸增大，各網(wǎng)格新增點(diǎn)定位耗時(shí)均呈現(xiàn)先減小，后增大的趨勢(shì)。這主要是因?yàn)樗饕龎K較少時(shí)，每個(gè)索引塊內(nèi)點(diǎn)數(shù)及三角形數(shù)量過(guò)多，塊內(nèi)定位的開(kāi)銷(xiāo)仍比較大，隨著索引矩陣參數(shù)增大，這樣的開(kāi)銷(xiāo)越來(lái)越小。但當(dāng)索引矩陣參數(shù)過(guò)大時(shí)，過(guò)于密集的索引矩陣導(dǎo)致單索引塊中所含點(diǎn)數(shù)及三角形數(shù)量很少，在進(jìn)行塊內(nèi)定位時(shí)會(huì)產(chǎn)生較多Cachemiss 現(xiàn)象，反而拖慢定位速度。從圖中可以看出，粗、中、細(xì)網(wǎng)格的最優(yōu)矩陣參數(shù)分別為30、70、100。本文規(guī)定定位耗時(shí)不超過(guò)最低耗時(shí)0.1 倍的部分為最優(yōu)區(qū)間，則上述網(wǎng)格的最優(yōu)區(qū)間分別為[21,37]、[40,140]、[70,300]。

綜合粗、中、細(xì)網(wǎng)格的原有點(diǎn)數(shù)、新增點(diǎn)數(shù)及索引矩陣參數(shù)最優(yōu)值，本文嘗試不同的系數(shù)配置，并最終確定了合適的系數(shù)配置，a為0.2，b為0.3。該系數(shù)配置可以保證3 個(gè)網(wǎng)格下的公式推測(cè)值均落入最優(yōu)區(qū)間中，如表2 所示。為進(jìn)一步驗(yàn)證公式準(zhǔn)確性，本文構(gòu)造了點(diǎn)數(shù)更多的隨機(jī)生成網(wǎng)格作為驗(yàn)證網(wǎng)格，并測(cè)試了其不同索引矩陣參數(shù)下的定位耗時(shí)，結(jié)果如圖5。該網(wǎng)格的公式推測(cè)值也在表2 中一并列出，可見(jiàn)推測(cè)值仍位于圖中的最優(yōu)區(qū)間內(nèi)。

4.4 原有剖分外擴(kuò)參數(shù)實(shí)驗(yàn)

本節(jié)探究了3.3.2 小節(jié)中外擴(kuò)參數(shù)對(duì)新增點(diǎn)定位耗時(shí)和三角剖分動(dòng)態(tài)更新耗時(shí)的影響。該實(shí)驗(yàn)使用單次增量插點(diǎn)法，基于細(xì)密度的隨機(jī)生成網(wǎng)格，索引矩陣參數(shù)采用4.3 節(jié)中測(cè)試出的最優(yōu)值100。實(shí)驗(yàn)分別統(tǒng)計(jì)了多種外擴(kuò)參數(shù)下兩個(gè)主要階段的耗時(shí)，并記錄了非空索引塊數(shù)量占總索引塊數(shù)量的比重。

Table 2 Prediction of optimal indexing map parameters using grids with different resolution表2 不同分辨率網(wǎng)格下最優(yōu)索引矩陣參數(shù)的推測(cè)

Fig.5 Tendency of locating time with indexing map parameter using very fine grid圖5 極細(xì)網(wǎng)格下定位耗時(shí)隨索引矩陣參數(shù)變化趨勢(shì)

測(cè)試結(jié)果見(jiàn)表3，隨著外擴(kuò)系數(shù)的增大，新增點(diǎn)定位耗時(shí)的增加非常明顯，這是因?yàn)橥鈹U(kuò)系數(shù)增大會(huì)導(dǎo)致索引矩陣范圍的增大，進(jìn)而導(dǎo)致新增點(diǎn)位于越來(lái)越少的索引矩陣塊中，索引效率降低，這可從非空索引塊占比的逐漸降低中看出。另外，三角剖分動(dòng)態(tài)更新耗時(shí)隨著外擴(kuò)系數(shù)的增加而輕微增加，這主要是因?yàn)槠史址秶脑龃髮?dǎo)致外層三角形被“拉長(zhǎng)”，新增點(diǎn)在單個(gè)三角形中的分布不均勻，進(jìn)而導(dǎo)致三角形的分裂次數(shù)增多，耗時(shí)增加。由該實(shí)驗(yàn)可知，為保證程序運(yùn)行的高效性，應(yīng)使用盡量小的外擴(kuò)參數(shù)。

Table 3 Running time of two main steps using different expansion parameters表3 不同外擴(kuò)參數(shù)下的兩個(gè)主要階段耗時(shí)

4.5 增量場(chǎng)景整體性能

增量場(chǎng)景的性能實(shí)驗(yàn)由兩部分組成：TID 算法與PatCC1 原始算法對(duì)比實(shí)驗(yàn)、增量開(kāi)銷(xiāo)實(shí)驗(yàn)。這兩個(gè)性能實(shí)驗(yàn)均將輸入點(diǎn)分為10 組，每組的點(diǎn)數(shù)均已在結(jié)果表格中給出。

PatCC1 在更新三角剖分時(shí)采用了刪除原有剖分并重新求解的方法，本文以該方法作為對(duì)照組，使用多次增量式插點(diǎn)方法進(jìn)行了性能對(duì)比實(shí)驗(yàn)。結(jié)果如表4，TID 算法相比PatCC1 原始算法具有很明顯的速度優(yōu)勢(shì)?？梢?jiàn)若將TID 算法用于PatCC1 中，將會(huì)給其帶來(lái)較大的性能提升。

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

增量插點(diǎn)功能的引入也會(huì)帶來(lái)額外開(kāi)銷(xiāo)，其主要來(lái)源于索引矩陣的構(gòu)建及新增點(diǎn)的定位。為探究該開(kāi)銷(xiāo)大小，本文分別統(tǒng)計(jì)了TID 算法對(duì)相同網(wǎng)格采用不同插點(diǎn)方式（非增量式和多次增量式）時(shí)的耗時(shí)，實(shí)驗(yàn)結(jié)果見(jiàn)表5。表中序號(hào)n對(duì)應(yīng)的增量單次耗時(shí)即為增量插入第n組新點(diǎn)并更新三角剖分所消耗的時(shí)間，非增量耗時(shí)即為將第0 組到第n組的所有點(diǎn)一次性插入并進(jìn)行三角剖分的耗時(shí)。假設(shè)當(dāng)前增量耗時(shí)為T(mén)inc，上一次非增量耗時(shí)T0,當(dāng)前非增量耗時(shí)Tc，則額外開(kāi)銷(xiāo)占比。從結(jié)果中可以看出，增量功能并不會(huì)引入過(guò)多額外開(kāi)銷(xiāo)，其在較好情況下與非增量式的耗時(shí)基本相近，在較壞情況下額外開(kāi)銷(xiāo)仍低于10%。

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

5 總結(jié)

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

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