尤 磊, 晏成名, 宋新宇

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

0 引言

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

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

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

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

1 自適應(yīng)二分

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

1.1 引導(dǎo)線構(gòu)建

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

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

圖1中實(shí)線與虛線相交于點(diǎn)集的重心點(diǎn)po，即po坐標(biāo)是點(diǎn)集中所有點(diǎn)坐標(biāo)的平均值。不難發(fā)現(xiàn)，無(wú)論點(diǎn)集如何分布，過(guò)點(diǎn)p且沿著第二主方向的直線l可將點(diǎn)集劃分為2個(gè)數(shù)量相近的子集。因此本文以第二主方向作為引導(dǎo)線方向?qū)c(diǎn)集進(jìn)行二分。

(1)

(2)

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

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

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

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

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

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

1.3 點(diǎn)集二分

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

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

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

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

1.4 邊界邊計(jì)算

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

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

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

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

1.5 子集再二分

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

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

圖3 所示點(diǎn)集的一個(gè)子集劃分結(jié)果Fig. 3 The result of a subset partition of point set

1.6 特殊情況處理

點(diǎn)集的幾何拓?fù)浣Y(jié)構(gòu)復(fù)雜，上述計(jì)算過(guò)程難以適用于所有情況，下面列出幾種特殊情況下的處理策略。

1.6.1 點(diǎn)類(lèi)別難以判定

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

圖4 點(diǎn)類(lèi)別判定的2種情況Fig. 4 Two cases of class determination

1.6.2 跨類(lèi)別邊界邊

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

1.6.3 后續(xù)拓展策略

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

2 并行算法設(shè)計(jì)

2.1 算法基本過(guò)程

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

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

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

算法1: 點(diǎn)集自適應(yīng)二分

輸入.點(diǎn)集Q與Q的邊界邊;

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

Step1. 角度和初始化；

Step2. 引導(dǎo)線構(gòu)建；

Step3. 引導(dǎo)起始點(diǎn)構(gòu)建；

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

Step5. 點(diǎn)集二分得到左右子集QL和QR；

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

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

算法2: 并行點(diǎn)集自適應(yīng)二分

輸入. 點(diǎn)集Q，其邊界邊結(jié)構(gòu)體實(shí)例的容器A與子集分割閾值Ns;

輸出.不超過(guò)子集分割閾值Ns個(gè)點(diǎn)的子集與其邊界邊結(jié)構(gòu)體實(shí)例的容器B。

Step1. 并行訪問(wèn)容器A中每一個(gè)結(jié)構(gòu)體實(shí)例，若其點(diǎn)集中點(diǎn)個(gè)數(shù)不大于Ns，則將該結(jié)構(gòu)體插入B的尾部，否則轉(zhuǎn)入Step2；

Step2. 對(duì)點(diǎn)集與其邊界邊的結(jié)構(gòu)體執(zhí)行點(diǎn)集自適應(yīng)二分，并將其二分結(jié)果按照點(diǎn)集中點(diǎn)的數(shù)量分別插入A和B的尾部；

Step3. 刪除當(dāng)前訪問(wèn)的A中每一個(gè)結(jié)構(gòu)體實(shí)例；

Step4. 轉(zhuǎn)入Step1直至A中元素為空。

在得到點(diǎn)集與其邊界邊結(jié)構(gòu)體實(shí)例容器的基礎(chǔ)上，采用文獻(xiàn)[8]的改進(jìn)并行算法對(duì)容器中B的子集并行構(gòu)建Delaunay三角網(wǎng)。

2.2 算法復(fù)雜度分析

2.2.1 時(shí)間復(fù)雜度

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

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

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

本文的點(diǎn)集二分是把問(wèn)題規(guī)模為n的問(wèn)題分解2個(gè)規(guī)模為n/2的問(wèn)題，本質(zhì)上雖未提升時(shí)間復(fù)雜度，但使用并行計(jì)算可提升計(jì)算效率。

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

2.2.2 空間復(fù)雜度

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

3 實(shí)驗(yàn)

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

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

表1 算法運(yùn)行時(shí)間比較Tab. 1 Comparison of algorithm running time

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

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

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

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

圖6 算法A與算法B的運(yùn)行時(shí)間比較圖Fig. 6 Running time comparison of algorithm A and algorithm B

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

4 結(jié)論

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

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

在平面Delaunay三角網(wǎng)構(gòu)建過(guò)程中，若一點(diǎn)的三角形已經(jīng)構(gòu)建完成，則該點(diǎn)不參與后續(xù)的構(gòu)網(wǎng)過(guò)程。這對(duì)于非凸包點(diǎn)來(lái)說(shuō)，即是該點(diǎn)的點(diǎn)夾角和達(dá)到360°，而對(duì)于凸包點(diǎn)來(lái)說(shuō)，需要結(jié)合其相鄰?fù)拱c(diǎn)的凸包邊來(lái)確定其點(diǎn)夾角和?；谖墨I(xiàn)[8]，使用點(diǎn)夾角和來(lái)移除優(yōu)先點(diǎn)和執(zhí)行后續(xù)拓展過(guò)程，這也導(dǎo)致本文算法時(shí)間復(fù)雜度受凸包計(jì)算復(fù)雜度的限制。嘗試構(gòu)建不需要凸包計(jì)算的改進(jìn)算法可進(jìn)一步提升構(gòu)網(wǎng)性能。

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

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