王　雯，吳　蔚，蘇天赟

(1.中國海洋大學(xué) 信息科學(xué)與工程學(xué)院，山東 青島 266100；2.國家海洋局第一海洋研究所，山東 青島 266061)

?

一種快速二維Delaunay三角網(wǎng)點(diǎn)定位算法

王雯1，吳蔚1，蘇天赟2

(1.中國海洋大學(xué) 信息科學(xué)與工程學(xué)院，山東 青島 266100；2.國家海洋局第一海洋研究所，山東 青島 266061)

摘要：在構(gòu)建二維Delaunay三角網(wǎng)的逐點(diǎn)插入法中，定位待插點(diǎn)所在三角形的快慢是影響整個算法構(gòu)網(wǎng)速度的關(guān)鍵因素。針對目前已有算法存在的搜索路徑長、搜索路徑求解計算量大等問題，結(jié)合三角形重心的幾何性質(zhì)，對點(diǎn)定位算法進(jìn)行改進(jìn)，避免求三角形重心和相交邊的過程。實(shí)驗(yàn)結(jié)果表明，文中算法較目前其他點(diǎn)定位算法能夠有效地縮短搜索路徑，減少點(diǎn)定位的計算時間，提高Delaunay三角網(wǎng)構(gòu)網(wǎng)過程中點(diǎn)定位的效率。

關(guān)鍵詞：Delaunay三角網(wǎng)；逐點(diǎn)插入法；點(diǎn)定位算法；三角形重心

三角剖分是計算幾何領(lǐng)域的重要研究內(nèi)容之一，而其中Delaunay三角剖分由于具有許多優(yōu)良性質(zhì)，在地理信息系統(tǒng)、數(shù)字高程插值、有限元分析、幾何重構(gòu)等領(lǐng)域有著廣泛的應(yīng)用[1-2]。目前構(gòu)建Delaunay三角網(wǎng)應(yīng)用最多的方法是逐點(diǎn)插入法，而確定包含待插點(diǎn)的三角形(目標(biāo)三角形)是逐點(diǎn)插入法中的一個重要步驟，一般將其稱為點(diǎn)定位問題。通過改進(jìn)點(diǎn)定位算法，可以有效提高Delaunay三角網(wǎng)構(gòu)建效率。劉學(xué)軍等提出了一種基于面積坐標(biāo)的點(diǎn)定位算法，利用建立的三角形間的拓?fù)潢P(guān)系和三角形面積坐標(biāo)判斷定位，這種方法的計算效率高，但是由于搜索路徑不唯一，造成了搜索路徑中很多不必要的“繞彎”現(xiàn)象，增加了判斷三角形的個數(shù)[3]。劉少華等利用點(diǎn)和有向線段的正負(fù)劃分來定位,該方法的本質(zhì)與面積坐標(biāo)法相同，計算效率高,但搜索路徑同樣不唯一[4]。宋占峰等提出最速方向定位法，通過將第一個要搜索的三角形重心指向待插點(diǎn)的連線作為方向線，如果點(diǎn)不在三角形內(nèi)，則與方向線相交的邊相鄰的三角形就為下一搜索三角形，該方法計算效率較低，但搜索路徑唯一[5]。蒲浩等提出重心方向搜索法，通過判斷當(dāng)前三角形的重心點(diǎn)和待插點(diǎn)相對于三角形三條邊的關(guān)系來定位，雖然搜索路徑是目前幾種算法中最短的，但計算效率也最低[6]。張詠等提出點(diǎn)定位融和算法，將文獻(xiàn)[3]與文獻(xiàn)[6]的方法相結(jié)合，先用面積坐標(biāo)法判斷搜索路徑，當(dāng)遇到結(jié)果不唯一的情況時，將當(dāng)前判斷的三角形重心與待插點(diǎn)相連，連線經(jīng)過的與當(dāng)前判斷三角形相鄰的三角形即為下一個要判斷的三角形[7]。這種方法的搜索路徑與文獻(xiàn)[6]相同，但用面積坐標(biāo)法判斷結(jié)果不唯一時，還是要計算三角形重心和判斷連線經(jīng)過的三角形，計算效率仍有待提升。

本文針對以上算法的不足，對點(diǎn)定位算法進(jìn)行了改進(jìn)，當(dāng)遇到不確定的情形時，利用三角形重心的性質(zhì)和點(diǎn)與三角形三邊位置關(guān)系的計算結(jié)果，可以快速判斷出重心與待插點(diǎn)連線經(jīng)過的相鄰三角形，計算效率得到了提高，搜索路徑與文獻(xiàn)[6-7]相同且為目前方法中的最短路徑。

1算法概述

在構(gòu)建Delaunay三角網(wǎng)的3種主要方法中，與分治算法、三角網(wǎng)生長法相比，逐點(diǎn)插入法具有算法簡單、占用空間小、便于動態(tài)更新等優(yōu)點(diǎn)。其主要流程如下：

1)先構(gòu)建包圍所有點(diǎn)集的初始矩形凸殼并連接一條對角線形成初始三角網(wǎng)；

2)將點(diǎn)集中的點(diǎn)按照一定的規(guī)律進(jìn)行排序，依次插入到三角網(wǎng)中；

3)對于每一個待插點(diǎn)，首先要找到待插點(diǎn)所在的三角形，然后從這個三角形開始對其每一條邊的鄰接三角形進(jìn)行外接圓檢測，進(jìn)而得到包含當(dāng)前待插點(diǎn)所有影響域的空腔，再將待插點(diǎn)與空腔頂點(diǎn)相連并更新三角形的鄰接關(guān)系即完成了一個待插點(diǎn)的插入；

4)將所有點(diǎn)插入完畢后，刪除與初始矩形凸殼4個頂點(diǎn)相連的輔助三角形，即完成了Delaunay三角網(wǎng)的構(gòu)建。

逐點(diǎn)插入法的流程如圖1所示。

圖1　逐點(diǎn)插入法流程

在逐點(diǎn)插入法的Delaunay構(gòu)網(wǎng)過程中，對于每一個待插點(diǎn)都首先要進(jìn)行點(diǎn)定位過程，即定位該點(diǎn)所在三角形，只有正確找到待插點(diǎn)所在三角形才能進(jìn)行建立空腔和三角網(wǎng)更新。點(diǎn)定位算法的主要思想是先確定第一個要搜索的三角形(首三角形)，然后利用數(shù)學(xué)方法判斷點(diǎn)是否在三角形中。如果不在，則利用待插點(diǎn)和這個三角形的位置關(guān)系，找到下一個要判斷的三角形。在這個過程中，要判斷的三角形會逐漸與待插點(diǎn)靠近，直到找到目標(biāo)三角形。這些判斷過的三角形會形成一條路徑，稱為搜索路徑。點(diǎn)定位算法的效率主要受以下3個因素影響：

Factor1：首三角形與待插點(diǎn)所在三角形的距離；

Factor2：搜索路徑中三角形的數(shù)量；

Factor3：判斷點(diǎn)與三角形位置關(guān)系時的計算公式的復(fù)雜度。

對于Factor1可以采用分塊建立局部三角形區(qū)域索引，或者把前一個待插點(diǎn)最后一個更新的三角形作為首三角形的方法。由于選取首三角形的方法主要取決于逐點(diǎn)插入法中點(diǎn)的插入順序，本文對此不作重點(diǎn)闡述。對于Factor2，如果選擇的搜索路徑中三角形數(shù)量較多會增加點(diǎn)與三角形關(guān)系的判斷次數(shù)；對于Factor3，如果判斷點(diǎn)和三角形位置關(guān)系時采用復(fù)雜的計算方法，會使待插點(diǎn)與每個三角形的判斷的時間增加。這些都會對逐點(diǎn)插入法的整體效率產(chǎn)生很大影響。因此，簡化點(diǎn)定位過程，減少點(diǎn)定位過程中對點(diǎn)和邊關(guān)系的判斷次數(shù)和計算復(fù)雜度，對提高整體Delaunay的構(gòu)網(wǎng)效率有著重要作用。本文重點(diǎn)針對Factor2和Factor3兩個因子，對點(diǎn)定位算法進(jìn)行改進(jìn)，以提高Delaunay構(gòu)網(wǎng)效率。

2算法改進(jìn)

2.1分區(qū)原則

從定位角度考慮，可以把三角形所在平面分成三大區(qū)域。Ⅰ區(qū)為三角形所在區(qū)域；Ⅱ區(qū)為三角形任意一邊和另兩邊延長線組成區(qū)域；Ⅲ區(qū)為任意兩條邊的延長線交叉組成區(qū)域。待插點(diǎn)與當(dāng)前判斷三角形的位置關(guān)系可以通過待插點(diǎn)與三角形任意兩頂點(diǎn)形成新三角形的面積正負(fù)來確定[8]。如圖2所示，設(shè)平面任意△ABC，其3個頂點(diǎn)按照逆時針順序編號，坐標(biāo)分別為A(xa,ya)，B(xb,yb)，C(xc,yc)；對任意點(diǎn)P(xp,yp)，P與三角形三邊AB，BC，CA構(gòu)成的三角形面積分別為SA，SB，SC，稱邊AB，BC，CA分別為SA，SB，SC的對應(yīng)邊，則

(1)

(2)

(3)

圖2　點(diǎn)P位于不同區(qū)域的示意圖

當(dāng)SA，SB，SC均大于等于0時，P位于△ABC內(nèi)或其邊上，即為Ⅰ區(qū)，如圖2(a)所示；當(dāng)SA，SB，SC有一個小于0時，P位于相應(yīng)的小于0的S對應(yīng)邊的外側(cè)，即為Ⅱ區(qū)，如圖2(b)所示；當(dāng)SA，SB，SC有兩個小于0時，P點(diǎn)位于相應(yīng)的小于0的S對應(yīng)的兩條邊的外側(cè)，即為Ⅲ區(qū)，如圖2(c)所示。其中將點(diǎn)在三角形邊上的情況歸于Ⅰ區(qū)是為了避免在三角形邊上的點(diǎn)無法屬于任何三角形的情況；將待插點(diǎn)在三角形三邊延長線的情況歸于Ⅱ區(qū)是因?yàn)辄c(diǎn)在Ⅲ區(qū)的處理方式比Ⅱ區(qū)復(fù)雜，這樣可以簡化計算。

2.2路徑判斷

2.2.1判斷方法

對于待插點(diǎn)位于Ⅰ區(qū)和Ⅱ區(qū)的情況，面積坐標(biāo)法可以很快找到最短路徑中下一個判斷的三角形。而對于待插點(diǎn)位于Ⅲ區(qū)的情況，目前的點(diǎn)定位算法文獻(xiàn)[3-7]的判斷處理方法各有不同，這些算法有些會增加Factor2中的三角形數(shù)量，有些需要求三角形重心坐標(biāo)和相交邊，從而增加Factor3中的計算公式復(fù)雜度。本文通過對點(diǎn)定位算法的影響因素的分析，對點(diǎn)位于Ⅲ區(qū)時的處理過程進(jìn)行了優(yōu)化。

針對Factor2，文獻(xiàn)[7]已證明其算法的搜索路徑是目前已有算法中的最短路徑。根據(jù)文獻(xiàn)[7]中所述，當(dāng)待插點(diǎn)位于當(dāng)前判斷三角形的Ⅲ區(qū)時，將當(dāng)前判斷三角形的重心與待插點(diǎn)相連，判斷與連線相交的是三角形的哪條邊，則下一個要判斷的三角形就是該邊的臨邊三角形。本文算法采用與文獻(xiàn)[7]相同的搜索路徑，即目前算法中的最短路徑。

針對Factor3，當(dāng)待插點(diǎn)位于Ⅲ區(qū)時，文獻(xiàn)[5-7]都無法避免求三角形重心坐標(biāo)和相交邊的過程，每判斷一個三角形就要進(jìn)行重心坐標(biāo)計算和多次直線方程求交計算，影響點(diǎn)定位算法的計算效率。本文提出了一種簡化的判別準(zhǔn)則，即：

當(dāng)待插點(diǎn)位于Ⅲ區(qū)時，SA，SB，SC中，最小的S所對應(yīng)的邊即為當(dāng)前判斷三角形的重心與待插點(diǎn)連線穿過的邊。

采用這種判別準(zhǔn)則，可以有效避免求三角形重心和相交邊的過程，從而可以提高待插點(diǎn)位于Ⅲ區(qū)時定位三角形的計算效率。

2.2.2方法證明

算法示意圖如圖3所示。

圖3(a)所示，設(shè)△ABC為當(dāng)前判斷的三角形，P′是位于Ⅲ區(qū)的待插點(diǎn)。分別將△ABC的3個頂點(diǎn)C，A，B與對邊中點(diǎn)M1，M2，M3相連，由三角形重心定義可知CM1，AM2，BM3相交于一點(diǎn)即為重心G。將線GB沿GB方向延長，延長線將邊AB，CB延長線交叉形成的Ⅲ區(qū)分成a，b兩個區(qū)域。首先假設(shè)P′位于a區(qū)，過P′分別向AB和CB延長線做垂線，垂足分別為D和E′。將直線DP′沿DP′方向延長交GB延長線于P。過P向CB延長線作垂線，垂足為E。過G分別向AB和BC作垂線，垂足分別為F和H。則由三角形中線的性質(zhì)可知,△ABC和△APC被各自中線BM3和PM3分得的三角形面積相等，進(jìn)而可以得到S△ABP=S△CBP，即

圖3　算法示意圖

(4)

由于點(diǎn)P′位于a區(qū)，DP′長度一定小于DP。將圖3(a)中的四邊形BDPE放大，放大后如圖3(b)所示，過P′作P′D′平行PB交BD與D′，過D′向P′E′作垂線，垂足為I?？芍鰾EP∽△D′IP′，△BDP∽△D′DP′，進(jìn)而得到

由圖中可知P′I長度一定小于P′E′，則

(5)

(6)

由于待插點(diǎn)在由AB，CB延長線組成的Ⅲ區(qū)時必有SA<0，SB<0，SC>0，最小值一定在SA，SB中。進(jìn)而可以由式(6)得到SB為最小值。所以當(dāng)待插點(diǎn)位于a區(qū)時，下一個要判斷的三角形為最小值SB對應(yīng)的BC邊的臨邊三角形，此時GP′一定與邊BC相交。

同理可證，當(dāng)待插點(diǎn)位于b區(qū)時最小值SA，則下一個要判斷的三角形為SA對應(yīng)的AB邊的臨邊三角形，此時GP′一定與邊AB相交，與求重心坐標(biāo)和求相交邊的方法得到的結(jié)果一致。

因此，本文提出的判別準(zhǔn)則可以取代求重心坐標(biāo)和相交邊的過程，用于待插點(diǎn)位于Ⅲ區(qū)時定位三角形的計算。這樣，本文算法既滿足了最短搜索路徑，又簡化了定位三角形的計算和判斷過程，使得整個點(diǎn)定位算法的計算效率得到了有效的提升。

2.3實(shí)現(xiàn)步驟

設(shè)待插點(diǎn)為P(xp,yp)，*pTri為當(dāng)前判斷的三角形指針，*pTri當(dāng)前指向△ABC，三頂點(diǎn)ABC逆時針排列，坐標(biāo)分別為A(xa,ya)，B(xb,yb)，C(xc,yc)，則算法的具體步驟如下：

Step1：將首三角形作為當(dāng)前判斷三角形，即將*pTri指向首三角形；

Step2：創(chuàng)建雙精度型數(shù)組s[3]、整型變量negativenum并賦值為0；

Step3：用式(1)、(2)、(3)求出s[0]=SA，s[1]=SB,s[2]=SC的值，且每當(dāng)s[i]<0時(i=0,1,2)，negativenum自加1；

Step4：若negativenum=0，則當(dāng)前判斷的三角形即為點(diǎn)P所在三角形，算法結(jié)束，否則進(jìn)行Step5；

Step5：若negativenum=1，則將相應(yīng)小于0的s[i]所對應(yīng)的邊(s[0]對應(yīng)AB邊，s[1]對應(yīng)BC邊，s[2]對應(yīng)CA邊)的臨邊三角形作為下一個要判斷的三角形，并將*pTri指向該三角形，negativenum清零，執(zhí)行Step3，否則執(zhí)行Step6；

Step6：若negativenum=2，則比較s[0]、s[1]、s[2]的大小，將其中最小的s[i]所對應(yīng)的邊的臨邊三角形作為下一個要判斷的三角形，并將*pTri指向該三角形，negativenum清零，執(zhí)行Step3。

3實(shí)驗(yàn)結(jié)果及分析

為了測試本文算法的效率，基于點(diǎn)插入順序?yàn)殡S機(jī)排列的逐點(diǎn)插入構(gòu)網(wǎng)方法對面積坐標(biāo)法[3]、重心方向搜索法[5]、最速方向定位法[6]、點(diǎn)定位融和法[7]以及本文算法進(jìn)行了測試，測試環(huán)境為Intel(R) Core(TM)2 Duo CPU T6670 @ 2.20GHz，4GB內(nèi)存，win7，64位操作系統(tǒng)，測試點(diǎn)集為點(diǎn)數(shù)為1 000萬的均勻分布的點(diǎn)集。經(jīng)測試，這5種算法的時間復(fù)雜度均為O(N)。為使測試結(jié)果更加直觀，引入平均每個待插點(diǎn)搜索路徑中的三角形個數(shù)(NP)，NP值越小搜索路徑越短。各種方法中每個待插點(diǎn)的首三角形都是相同的，選取規(guī)則都是把前一個插入點(diǎn)最后一個更新的三角形作為首三角形。測試結(jié)果如表1所示。

表1　不同算法效率對比

從表中可以看出，對于搜索路徑，重心方向搜索法、點(diǎn)定位融和法和本文算法的NP值相同且為最小，因?yàn)檫@3種方法的搜索路徑相同且是這幾種方法中最短的；相比這3種算法，最速方向定位法的搜索路徑略長，所以NP值也略高；面積坐標(biāo)法由于搜索路徑中很多繞彎現(xiàn)象的存在，NP值最高。對于算法效率，面積坐標(biāo)法由于無需求重心坐標(biāo)和相交邊，雖然判斷的三角形最多(NP值最大)但算法效率依舊很高；而重心方向搜索法對于每個判斷三角形都要求重心坐標(biāo)和相交邊，雖然搜索路徑最短但算法效率最低，搜索時間約為面積坐標(biāo)法的2倍；最速方向定位法將重心方向搜索法中對每個三角形重心坐標(biāo)的計算簡化為只求首三角形的重心坐標(biāo)，所以搜索時間比重心方向搜索法略少，但依舊比面積坐標(biāo)法多出約60%的搜索時間；點(diǎn)定位融和法將面積坐標(biāo)法和重心方向搜索法結(jié)合，搜索時間比重心方向搜索法降低了約10%；而本文算法由于搜索路徑最短，計算效率與面積坐標(biāo)法相當(dāng)，所以最終算法效率是5種方法中最高的。

4結(jié)束語

本文在現(xiàn)有點(diǎn)定位算法基礎(chǔ)上，結(jié)合三角形重心的幾何性質(zhì)，提出了一種新的點(diǎn)定位算法，對點(diǎn)定位算法中待插點(diǎn)位于三角形兩條邊外側(cè)時的計算效率進(jìn)行了優(yōu)化，避免了求三角形重心和相交邊的過程。測試結(jié)果表明，本文算法不僅搜索路徑短，而且算法效率高，可以有效地提升Delaunay構(gòu)網(wǎng)算法的效率，具有一定的應(yīng)用價值。

參考文獻(xiàn)：

[1]謝增廣.平面點(diǎn)集Delaunay三角剖分的分治算法[J].計算機(jī)工程與設(shè)計，2012，33(1)：2652-2658.

[2]李小秋，許民獻(xiàn)，尹志永.Delaunay三角網(wǎng)關(guān)鍵技術(shù)探討[J].測繪工程，2011，20(6)：61-67.

[3]劉學(xué)軍，符鋅砂，趙建三.三角網(wǎng)數(shù)字地面模型快速構(gòu)建算法研究[J].中國公路學(xué)報，2000，13(2)：31-36.

[4]劉少華，吳東勝，羅小龍，等.Delaunay三角網(wǎng)中點(diǎn)目標(biāo)快速定位算法研究[J].測繪科學(xué)，2007，32(2)：69-71.

[5]宋占峰，蒲浩，詹振炎.基于三角網(wǎng)數(shù)字地面模型快速定位算法的研究[J].中國鐵道科學(xué)，2002，23(1)：63-66.

[6]蒲浩，宋占峰，詹振炎.快速構(gòu)建三角網(wǎng)數(shù)字地面模型方法的研究[J].中國鐵道科學(xué)，2001，22(6)：100-105.

[7]張詠，劉長星，楊瑜華，等.基于融合算法的二維Delaunay三角網(wǎng)任意點(diǎn)定位研究[J].測繪科學(xué)，2010，35(2)：84-87.

[8]JOHN B M.Transformation of trilinear and quadriplanar coordinates to and from cartesian coordinates[J].American Mineralogist,1964,49(7,8)：926-936.

[責(zé)任編輯：劉文霞]

A rapid algorithm for point positioning in 2D delaunay triangulation

WANG Wen1,WU Wei1,SU Tianyun2

(1.School of Information Science and Engineering,Ocean University of China,Qingdao 266100,China；2.State Oceanic Administration,First Institute of Oceanography,Qingdao 266061,China)

Abstract：In incremental insertion algorithms of 2D Delaunay triangulation,seeking out the triangle which the inserting point locates in is the key factor influencing the efficiency.In this paper,point positioning algorithm is improved by avoiding the calculation of gravity center and intersecting edge,which makes the use of geometric properties for triangle barycenter to solve the problem that searching path is too long and complicated to calculate in present algorithms.The experimental results show the algorithm in this paper can shorten the searching path and reduce the time of point positioning process,ultimately improve the efficiency of delaunay triangulation compared with other present point positioning algorithms.

Key words：delaunay triangulation；incremental insertion algorithm；point positioning algorithm；triangle barycenter

中圖分類號：TP391

文獻(xiàn)標(biāo)識碼：A

文章編號：1006-7949(2016)03-0025-05

作者簡介：王雯(1989-)，女，碩士研究生.

基金項(xiàng)目：國家科技重大專項(xiàng)資助項(xiàng)目(2011ZX05056-001-01);海洋公益性行業(yè)科研專項(xiàng)資助項(xiàng)目(201205001)

收稿日期：2014-11-24；修回日期:2014-12-06