尤 磊, 晏成名, 宋新宇

(信陽師范學院 a.計算機與信息技術學院; b.河南省物聯(lián)網(wǎng)與智能安防工程研究中心;c. 數(shù)學與統(tǒng)計學院, 河南 信陽 464000)

0 引言

Delaunay三角剖分具有最小角最大化與空外接圓的特征，這使得平面點集的三角剖分具有局部性、唯一性與合理性的優(yōu)良性質[1-2]。這些性質使得Delaunay三角剖分在表面重建、數(shù)字地形構建、有限元分析與虛擬現(xiàn)實等領域得到廣泛應用[3-4]。

傳統(tǒng)上Delaunay三角剖分算法分為3類：分治合并法、逐點插入法與三角網(wǎng)生長算法。分治合并法將點集劃分為多個子集，將各子集上構建的三角網(wǎng)進行合并[5-6]，以完成點集的Delaunay三角剖分；逐點插入法通過在三角網(wǎng)中定位新加入點的位置，并由此優(yōu)化三角網(wǎng)，該方法需要頻繁遍歷三角網(wǎng)中的三角形；三角網(wǎng)生長算法是在點集中搜索基邊的第3點以形成一個新的三角形[7]，可利用Delaunay三角網(wǎng)的局部性逐步減少點集規(guī)模以提升構網(wǎng)效率。概括而言，分治合并算法效率最高[8]。

相對于串行逐點構建三角形，并行Delaunay三角剖分可進一步提升構建效率。BLELLOCH等[9]通過點集在坐標軸上的中軸線構建Delaunay邊劃分子集，再分別構建子集的三角網(wǎng)；李堅等[10]通過對點集四叉樹剖分得到點集的若干子集，分別構建各子集的三角網(wǎng)后再融合得到全局三角網(wǎng)。張春亢等[11]通過構建三角網(wǎng)墻對點集縱橫向切開得到若干子集，各子集并行構網(wǎng)后再將子集構建的三角網(wǎng)直接合并，然而構建網(wǎng)墻需要使用額外的控制點，并需要刪除子集邊界處錯誤的三角形。FUNKE等[12]通過構建子集在邊界處的邊界三角網(wǎng)以融合子集的三角網(wǎng)，并分別采用共享內(nèi)存與分布式內(nèi)存的方式實現(xiàn)并行算法。上述并行算法均取得較好效果，但也存在如下問題：(1)點集子集劃分時，未考慮點集中點的分布情況；(2)在子網(wǎng)合并時都需要一些額外的輔助工作，如刪除邊界處的三角形或構建邊界三角網(wǎng)等。

以優(yōu)先點為中心的三角網(wǎng)生長算法[8]邊構網(wǎng)邊刪除已完成構建的點，通過逐漸減少點集的規(guī)模以提升構網(wǎng)效率。并行Delaunay三角剖分可提升計算效率，分治合并算法效率最高且適合構建并行算法?；诖?，本文結合上述3種算法優(yōu)勢，提出一種平面點集的自適應二分的并行Delaunay三角剖分算法。其主要過程包括：根據(jù)平面點集中的點的幾何分布，構建點集二分的引導線；基于引導線采用優(yōu)先點的生長算法，構建點集二分的優(yōu)先點，然后根據(jù)優(yōu)先點構建的三角形對點集二分；點集逐級并行二分，直至子集中點的數(shù)量小于設定分割閾值；然后分別對各個子集并行執(zhí)行Delaunay三角剖分算法。本文算法優(yōu)勢主要體現(xiàn)在：(1)自適應確定引導線以盡可能地均勻二分點集；(2)點集二分與各子集的Delaunay三角剖分均可并行運行；(3)點集二分過程中生成的三角網(wǎng)與子集構建的三角網(wǎng)均可直接合并，不需要額外處理。

1 自適應二分

分治或并行算法通常采用坐標軸分割的方式將點集劃分為不同子集[9-10，12]。其將點集在X(或Y)軸的起始區(qū)間均分為等距的若干段，X(或Y)坐標位于同一段的點屬于同一個子集。然而在實際應用中，點集中的點并不均勻分布。如果按照坐標軸等距劃分的方式，容易導致各子集間點的數(shù)量相差很大，難以發(fā)揮分治或并行算法的優(yōu)勢。為確保各子集間點的數(shù)量接近，提出自適應二分的點集劃分策略，其根據(jù)點集中點的分布計算二分的引導線，再根據(jù)引導線將點集劃分為兩個數(shù)量相近的子集。

1.1 引導線構建

點集P的協(xié)方差矩陣的特征值對應的特征向量可反映點集在不同方向上的離散程度。圖1中實線與虛線所在方向分別是點集協(xié)方差矩陣的最大與第二大特征值對應特征向量的方向，分別稱之為第一主方向和第二主方向。

圖1 點集的第一與第二主方向示意圖Fig. 1 Schematic diagram of the first and second main directions of the point set

圖1中實線與虛線相交于點集的重心點po，即po坐標是點集中所有點坐標的平均值。不難發(fā)現(xiàn)，無論點集如何分布，過點p且沿著第二主方向的直線l可將點集劃分為2個數(shù)量相近的子集。因此本文以第二主方向作為引導線方向對點集進行二分。

(1)

(2)

1.2 優(yōu)先點選定與拓展

得到引導起點ps與終點pe后，使用以優(yōu)先點為中心的Delaunay三角網(wǎng)生長算法[8]沿著引導線構建點集P的局部三角網(wǎng)，構建優(yōu)先點選擇策略確保優(yōu)先點從ps開始并沿著引導線搜索，直至引導終點pe結束。

優(yōu)先點選擇策略為：設當前優(yōu)先點為pc，pcn0與pcn1分別為pc構建三角形頂點中的另外2個頂點，若pcn0到引導終點pe的距離小于pcn1到pe的距離，且pcn0到引導線的距離小于pcn1到引導線的距離，則點pcn0為下一個優(yōu)先點。

Delaunay三角網(wǎng)的空外接圓性質可確保一個點必與其最近的一個點組成一條Delaunay邊。因此對于引導起點ps，可計算得到以ps為起點的一條Delaunay邊es。此時，以點ps為第一個優(yōu)先點，以邊es為第一條拓展邊(若恰巧邊es的兩個頂點都是凸包點，則有可能拓展失敗，此時以邊es的逆向邊為拓展邊)得到一個三角形，然后拓展新生成三角形中以點ps為頂點的另外一條邊，直至點ps完成拓展。此時，根據(jù)優(yōu)先點選擇策略得到下一個優(yōu)先點pcn，然后繼續(xù)逐優(yōu)先點地拓展直至下一個優(yōu)先點為引導終點pe，并將點pe可構成的Delaunay三角形全部生成。在拓展過程中，每增加一個三角形即更新該三角形三個頂點的點夾角和以提升拓展效率。

圖1所示點集的第一次拓展過程得到的三角形如圖2中藍色三角形所示，優(yōu)先點的引導線如圖2中紅色虛線所示。

圖2 優(yōu)先點拓展與點集二分示意圖Fig. 2 Schematic diagram of priority point expansion and point set dichotomy

1.3 點集二分

完成所有優(yōu)先點的拓展后，得到若干圍繞優(yōu)先點且相互鄰接的三角形網(wǎng)格。這些三角形將點集P分為3部分：(1)優(yōu)先點；(2)三角形左側未完成拓展的點；(3)三角形右側未完成拓展的點。由于優(yōu)先點已經(jīng)在優(yōu)先點拓展過程中得到，且這部分點屬于已完成拓展的點，其相應的三角形網(wǎng)格可直接保存至最終三角形網(wǎng)格中，后續(xù)構網(wǎng)也不需要考慮這些點。從而，點集二分是根據(jù)優(yōu)先點拓展后得到的三角形的幾何空間位置，將點集劃分為左右兩個未完成拓展的點的集合PL和PR。

根據(jù)圖2，未完成拓展的點有2種類型：未生成三角形的點與已生成部分三角形的點。對于未生成三角形的一個點pi，可直接根據(jù)該點與引導線的位置判定：如果三點按次序(ps、pe與pi)構成三角形的面積為正值，上述三點呈逆時針次序排列，則點pi屬于引導線右側未完成的拓展點集PR；若三角形面積為負值，則pi屬于引導線左側未完成的拓展點集PL。對于已生成部分三角形的點而言，首先判斷該點所在三角形的另外兩個點是否是優(yōu)先點，如果都是優(yōu)先點，則按照優(yōu)先點的生成次序和該點按序構成三角形的面積判定，面積為正則屬于PR，為負則屬于PL。

上述過程可將未生成三角形的點劃分到子集PR與PL中。對于已生成部分三角形的點，由于其生成三角形的頂點中有一個優(yōu)先點，此時可采用類別傳播策略：設未確定歸屬的點為pi，優(yōu)先點為pcn，另一點為pj，則若三點形成夾角∠pipcnpj角度不大于90°，則點pi與pj屬于同一類別；若pj已判定類別，則將pj的類別傳播給pi；若pj尚未判定類別，則設置pi=pj，繼續(xù)采取類別傳播策略計算pi的類別，直至可以計算出pi的類別。

經(jīng)過上述處理過程，絕大部分的點都可以準確劃分類別(但也存在一些特殊情況，將在1.6節(jié)給出一個后續(xù)拓展過程)。即點集P可根據(jù)優(yōu)先點拓展得到的三角形二分得到PL與PR兩個需要后續(xù)待生成三角網(wǎng)的子集。

1.4 邊界邊計算

根據(jù)圖2，子集PL與PR中點已生成的邊形成子集PL與PR的邊界。這些邊界上的邊稱之為邊界邊，同時也是優(yōu)先點所生成三角形的邊界邊的逆向邊。因此在后續(xù)對子集PL與PR構建三角形網(wǎng)格時，保留上述邊界邊可確保子集在三角形網(wǎng)格生成時受邊界邊的約束，進而可以直接將子集得到的三角網(wǎng)融入最終三角形網(wǎng)格中。

邊界邊的計算方法如下所述：對于優(yōu)先點生成三角形網(wǎng)格中的邊，若該邊沒有逆向邊，則該邊是優(yōu)先點生成三角形網(wǎng)格的邊界邊；若已經(jīng)有逆向邊，則該邊不是優(yōu)先點生成三角形網(wǎng)格的邊界邊；若該邊的兩個端點都屬于PL，則該邊的逆向邊是PL的邊界邊；同理得到PR的邊界邊。

在對子集進行自適應二分時，子集中已有點生成邊界邊。此時可根據(jù)邊界邊的兩個端點隸屬于PL與PR的規(guī)則進行處理。即，若邊界邊的兩個端點都屬于PL(PR)，則該邊界邊是PL(PR)的邊界邊。

根據(jù)圖2，在引導起點ps與引導終點pe附近可能存在著一條邊界邊的兩個端點分別屬于子集PL與PR的情況。這種特殊情況將在1.6節(jié)進行處理。

1.5 子集再二分

在得到子集PL與PR之后，若PL或PR中點的數(shù)量小于設定的子集分割閾值Ns，則無須對該子集繼續(xù)劃分，否則繼續(xù)對該子集進行劃分，直至所有子集中點的數(shù)量小于等于Ns。圖1所示點集(700個點)的一個最終子集劃分(Ns為50)，如圖3所示。形狀小且紅色的點為優(yōu)先點，形狀大且顏色各異的點為各個子集的點，不同子集采用不同的顏色表示?？梢钥闯?，點密度小的區(qū)域的子集較少，點密度大的區(qū)域的子集較多。

經(jīng)過上述步驟可得到點集P的若干個子集，每個子集中的點的數(shù)量不超過Ns，且包含相應的邊界邊。此時可采用并行算法對每個子集構建Delaunay三角網(wǎng)。本文將以優(yōu)先點為中心的Delaunay三角網(wǎng)生成算法[8]進行改進，使其支持邊界邊約束的平面點集Delaunay三角網(wǎng)生成。

圖3 所示點集的一個子集劃分結果Fig. 3 The result of a subset partition of point set

1.6 特殊情況處理

點集的幾何拓撲結構復雜，上述計算過程難以適用于所有情況，下面列出幾種特殊情況下的處理策略。

1.6.1 點類別難以判定

由于引導起點ps與引導終點pe更靠近點集的邊界，經(jīng)常出現(xiàn)其他點與點ps或pe夾角是鈍角的情況，如圖4(a)所示，點A與C是優(yōu)先點，點B可根據(jù)點A與C判定，屬于左側PL點集，E點可通過類別傳播判定屬于PL。但由于∠ECD>90°，因此E的類別不傳播給D，從而達到正確劃分類別的目的。對于圖4(b)，此時∠ECB和∠ECD均大于90°，從而B和D均不能將其類別傳播給E，這就導致E的類別難以判定。

圖4 點類別判定的2種情況Fig. 4 Two cases of class determination

1.6.2 跨類別邊界邊

在引導起點ps與引導終點pe位置處，有時出現(xiàn)邊界邊的一個端點屬于PL，另一個端點屬于PR的情況。圖4(a)中邊ED就是一條跨界邊界邊。此時，該邊界邊既不屬于PL，也不屬于PR。這導致后續(xù)網(wǎng)格生成不完整。

1.6.3 后續(xù)拓展策略

上述兩種情況中，一種情況是由于點的歸屬不好判定而導致點難以劃分給任意一個類別，另一種情況是邊不屬于任意一個子集的邊界邊。若將上述點和邊隨意劃分為任意一個類別，都將導致三角形生成錯誤。當出現(xiàn)上述情況時，由于點或邊并不歸屬于任一個子集，因此后續(xù)子集構網(wǎng)時也將不涉及該點或該邊，這將使最終構建的三角形出現(xiàn)空洞現(xiàn)象。為避免空洞出現(xiàn)，可根據(jù)點的角度和構建后續(xù)拓展策略。在算法初始化時，初始化點集每個點的角度和，凸包點的角度和為360°減去該點形成兩條凸包邊的外夾角角度，其他各點均為0°。在點集二分時，保存當前處理點集中各點形成的半邊及各點所形成夾角的夾角和，然后在所有子集網(wǎng)格生成之后，再對夾角和不是360°的點進行拓展處理：若以該點為起點的半邊的逆向邊不存在，則構建此逆向邊并拓展。

2 并行算法設計

2.1 算法基本過程

根據(jù)上節(jié)所述，采用點集自適應二分策略將點集劃分為兩個子集，子集再二分時與其他子集無直接聯(lián)系；在邊界邊的控制下，子集構建Delaunay三角網(wǎng)的過程是獨立的，從而子集在劃分過程中產(chǎn)生的三角形與子集構建的三角形可直接匯總至全局Delaunay三角網(wǎng)中?；谏鲜?，可構建一個并行的Delaunay三角網(wǎng)生成算法，其并行計算過程包括子集二分與子集Delaunay三角網(wǎng)構建。

并行算法的主要步驟為：初始化所有點的角度和；并行對點集進行自適應二分以得到若干個子集，使得每個子集中點的數(shù)量不超過Ns；并行計算每個子集的Delaunay三角網(wǎng)；再對點夾角和未達到360°的點做拓展。

下面給出主要算法的偽代碼。

算法1: 點集自適應二分

輸入.點集Q與Q的邊界邊;

輸出.點集Q的左右子集QL和QR，QL與QR的邊界邊，優(yōu)先點的三角網(wǎng)格與點集Q的點夾角和。

Step1. 角度和初始化；

Step2. 引導線構建；

Step3. 引導起始點構建；

Step4. 優(yōu)先點選定與拓展；

Step5. 點集二分得到左右子集QL和QR；

Step6. 獲取子集QL和QR的邊界邊。

為使點集二分過程和子集Delaunay三角網(wǎng)構建便于并行計算，可將子集及其邊界邊形成一個結構體存入容器中。

算法2: 并行點集自適應二分

輸入. 點集Q，其邊界邊結構體實例的容器A與子集分割閾值Ns;

輸出.不超過子集分割閾值Ns個點的子集與其邊界邊結構體實例的容器B。

Step1. 并行訪問容器A中每一個結構體實例，若其點集中點個數(shù)不大于Ns，則將該結構體插入B的尾部，否則轉入Step2；

Step2. 對點集與其邊界邊的結構體執(zhí)行點集自適應二分，并將其二分結果按照點集中點的數(shù)量分別插入A和B的尾部；

Step3. 刪除當前訪問的A中每一個結構體實例；

Step4. 轉入Step1直至A中元素為空。

在得到點集與其邊界邊結構體實例容器的基礎上，采用文獻[8]的改進并行算法對容器中B的子集并行構建Delaunay三角網(wǎng)。

2.2 算法復雜度分析

2.2.1 時間復雜度

以點的訪問作為算法的基本運算,算法的核心步驟為點集自適應二分和Delaunay三角網(wǎng)生成。

點集自適應二分時，計算凸包點的最優(yōu)時間復雜度為O(nlgn)[13]；引導線構建需要計算點集協(xié)方差矩陣的特征向量。盡管矩陣運算的時間復雜度為O(n3)，但平面點集的協(xié)方差矩陣的階數(shù)為2，相對于點集規(guī)模而言，協(xié)方差矩陣的計算量小于遍歷點集的計算量O(n)。從而，引導線構建時計算點集在第二主方向上的投影以及引導起始點選定的時間復雜度均可認為是O(n)。

優(yōu)先點擴展和子集Delaunay三角網(wǎng)生成均建立在以優(yōu)先點為中心的Delaunay三角網(wǎng)生成算法[8]的基礎上，因此，此部分的時間復雜度為：在最差情況下時間復雜度為O(n2)，最優(yōu)情況下時間復雜度為O(nlgn)。

本文的點集二分是把問題規(guī)模為n的問題分解2個規(guī)模為n/2的問題，本質上雖未提升時間復雜度，但使用并行計算可提升計算效率。

綜合上述，本文所提算法的時間復雜度為：在最差情況下時間復雜度為O(n2)，最優(yōu)情況下時間復雜度為O(nlgn)。

2.2.2 空間復雜度

在算法運行過程中，需要記錄點夾角和、點集二分后得到的邊界邊與運行Delaunay三角網(wǎng)生成算法。點夾角和是針對每一個點的，因此其空間復雜度為O(n)；邊界邊是以三角形的邊為基礎得到的，根據(jù)歐拉公式可知頂點和三角形的個數(shù)是線性關系，因此邊界邊的存儲的空間復雜度為O(n)；以優(yōu)先點為中心的Delaunay三角網(wǎng)生成算法的空間復雜度也是O(n)。因此本文算法的空間復雜度是O(n)。

3 實驗

采用OpenMP多線程并發(fā)編程API在VC++2017開發(fā)環(huán)境下實現(xiàn)本文所提算法，采用PCL點云庫[14]對點集進行訪問與存儲，在CPU型號為i7-8700K，內(nèi)存16 GB的臺式機上開展實驗。

為驗證算法有效性，將文獻[8]算法(算法A)與本文算法(算法B)進行比較以驗證本文算法的并行效果。分別以1.0萬、1.5萬、2.0萬、2.5萬與3.0萬個點的規(guī)模、在不同子集分割閾值下開展實驗。分別記錄本文算法點集二分最終結果中子集點數(shù)量的最大值、最小值與平均值、點集二分與子集三角網(wǎng)構建的時間。實驗結果數(shù)據(jù)如表1所示。算法B構建的一個三角網(wǎng)如圖5所示，為清晰顯示，此處顯示圖1所示點集構建的三角網(wǎng)為例，圖中同一顏色且相鄰的點屬于劃分后的同一個子集。

表1 算法運行時間比較Tab. 1 Comparison of algorithm running time

從表1中可以看出，本文的并行算法有較好表現(xiàn)，構網(wǎng)時間減少明顯。在不同的子集分割閾值下，點集劃分時間與子集構網(wǎng)時間有一定的偏差，但總體上子集構網(wǎng)時間略長，這是因為點集劃分時生成的三角形數(shù)量小于子集構網(wǎng)生成的三角形數(shù)量；子集分割閾值越大，點集劃分使用時間越少，這是因為閾值越大，子集包括的點數(shù)量越多，劃分的子集數(shù)量越少，劃分次數(shù)越少，因此耗時較少；子集構網(wǎng)時間并未跟隨子集分割閾值變化呈現(xiàn)規(guī)律性變化，這可能是因為子集構網(wǎng)時生成三角形較多，且更易受點幾何拓撲結構的影響；值得指出的是，算法B的總用時受子集分割閾值的影響不大，對于同一點集，不同閾值下的總用時相差不大，但總用時與點集規(guī)模直接相關。

圖5 算法B構建的圖1所示點集的三角網(wǎng)Fig. 5 Triangulation of point set shown in Fig. 1 constructed by algorithm B

表1中列出點集二分后子集中點數(shù)量的最大值、最小值與平均值。點數(shù)量最大值取決于分割閾值；點數(shù)量最小值并無明顯規(guī)律，也有一些數(shù)量極少的子集出現(xiàn)。當子集中點的數(shù)量大于分割閾值時，需要對該子集繼續(xù)二分，二分時產(chǎn)生的優(yōu)先點并不屬于后續(xù)的兩個二分集合；同時考慮到點的幾何分布復雜，這些都可能導致左右兩個子集點數(shù)量分布不均勻，但從平均值來看，平均值均略大于子集分割閾值的二分之一。

為直觀顯示算法A與算法B在運行時間上的差異，圖6給出算法A與算法B在子集分割閾值為900時的運行時間比較圖。

圖6 算法A與算法B的運行時間比較圖Fig. 6 Running time comparison of algorithm A and algorithm B

從圖6可以看出，與算法A相比，算法B隨著點集規(guī)模的增長趨勢更趨近于線性增長，這也充分說明本文所提算法的有效性。

4 結論

為提升Delaunay三角網(wǎng)的構建效率，提出了一種并行計算算法，該算法結合分治算法和生長算法的優(yōu)勢，通過并行算法將點集劃分為若干子集，然后并行地構建子集的三角網(wǎng)，計算過程中產(chǎn)生的三角形就是最終的三角形，不需要額外的合并操作，這充分發(fā)揮了Delaunay三角網(wǎng)的局部性與全局性，即只要確保每個點構建的三角形是Delaunay三角形，那么構建的結果就是最終輸出結果，從而確保本文算法可有效提升構網(wǎng)效率。

實際應用中，存在只需要構建點集的局部Delaunay三角網(wǎng)的情況。此種情況下，只需要根據(jù)使用情況構建局部Delaunay三角網(wǎng)即可。本文提出的基于引導線構建Delaunay三角網(wǎng)的策略可為構建局部Delaunay三角網(wǎng)提供一種思路：即根據(jù)需求構建一條或多條引導線，然后沿著引導線構建三角形，即可得到Delaunay局部的三角網(wǎng)。

在平面Delaunay三角網(wǎng)構建過程中，若一點的三角形已經(jīng)構建完成，則該點不參與后續(xù)的構網(wǎng)過程。這對于非凸包點來說，即是該點的點夾角和達到360°，而對于凸包點來說，需要結合其相鄰凸包點的凸包邊來確定其點夾角和。基于文獻[8]，使用點夾角和來移除優(yōu)先點和執(zhí)行后續(xù)拓展過程，這也導致本文算法時間復雜度受凸包計算復雜度的限制。嘗試構建不需要凸包計算的改進算法可進一步提升構網(wǎng)性能。

以OpenMP多線程并發(fā)編程API構建并行算法，然而基于GPU的多線程計算技術具有更好的計算效率[15-16]。將本文算法改進為基于GPU的并行算法也將進一步提升構網(wǎng)效率。

與二維點集不同，對于三維點集，點的點夾角和達到360°并不意味著該點的三角網(wǎng)已經(jīng)構建完成。這意味著難以將本文算法直接應用于三維點集的Delaunay網(wǎng)構建。但點集的自適應二分思想依然適用于三維點集，將本文算法進行改進使其適用于三維點集是下一步的研究內(nèi)容。