王　雯，吳　蔚，蘇天赟

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

?

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

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

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

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

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

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

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

1算法概述

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

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

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

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

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

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

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

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

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

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

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

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

2算法改進(jìn)

2.1分區(qū)原則

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

(1)

(2)

(3)

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

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

2.2路徑判斷

2.2.1判斷方法

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

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

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

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

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

2.2.2方法證明

算法示意圖如圖3所示。

圖3(a)所示，設(shè)△ABC為當(dāng)前判斷的三角形，P′是位于Ⅲ區(qū)的待插點(diǎn)。分別將△ABC的3個(gè)頂點(diǎn)C，A，B與對(duì)邊中點(diǎn)M1，M2，M3相連，由三角形重心定義可知CM1，AM2，BM3相交于一點(diǎn)即為重心G。將線GB沿GB方向延長(zhǎng)，延長(zhǎng)線將邊AB，CB延長(zhǎng)線交叉形成的Ⅲ區(qū)分成a，b兩個(gè)區(qū)域。首先假設(shè)P′位于a區(qū)，過P′分別向AB和CB延長(zhǎng)線做垂線，垂足分別為D和E′。將直線DP′沿DP′方向延長(zhǎng)交GB延長(zhǎng)線于P。過P向CB延長(zhǎng)線作垂線，垂足為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′長(zhǎng)度一定小于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長(zhǎng)度一定小于P′E′，則

(5)

(6)

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

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

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

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

設(shè)待插點(diǎn)為P(xp,yp)，*pTri為當(dāng)前判斷的三角形指針，*pTri當(dāng)前指向△ABC，三頂點(diǎn)ABC逆時(shí)針排列，坐標(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時(shí)(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]所對(duì)應(yīng)的邊(s[0]對(duì)應(yīng)AB邊，s[1]對(duì)應(yīng)BC邊，s[2]對(duì)應(yīng)CA邊)的臨邊三角形作為下一個(gè)要判斷的三角形，并將*pTri指向該三角形，negativenum清零，執(zhí)行Step3，否則執(zhí)行Step6；

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

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

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

表1　不同算法效率對(duì)比

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

4結(jié)束語(yǔ)

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

參考文獻(xiàn)：

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

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

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

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

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

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

[7]張?jiān)?，劉長(zhǎng)星，楊瑜華，等.基于融合算法的二維Delaunay三角網(wǎng)任意點(diǎn)定位研究[J].測(cè)繪科學(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

中圖分類號(hào)：TP391

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

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

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

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

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