李可 高清維 盧一相 孫冬 竺德

一組點的凸包是指包含這些點的最小凸多邊形[1],凸包問題是計算幾何和圖形學中最基礎的問題之一.平面凸包可以用來處理聚類分析[2?5],圍欄問題和城市規(guī)劃問題等等.在GIS 應用[6]中,凸包可以快速獲得區(qū)域地理的基本輪廓和最有效邊界信息,方便定位與數(shù)字化建模.在模式識別學科的研究中,凸包的應用對于優(yōu)化數(shù)據(jù)結構和降低數(shù)據(jù)計算規(guī)模提供了新思路.凸包算法[7?10]在人工智能、人臉識別[11?12]等前沿領域均有不同程度的應用.

常用的二維平面點集凸包算法有Graham 算法[13]、Jarvis 步進算法[14]和窮舉算法等,這些經(jīng)典算法在處理少量數(shù)據(jù)時表現(xiàn)出色,但是當平面點集數(shù)據(jù)量超過106數(shù)量級時效率較差,所以不能推廣到具有大規(guī)模數(shù)據(jù)集的工程中使用.為解決實際工程應用中具有超大規(guī)模的平面點集的凸包計算問題,出現(xiàn)了許多快速凸包算法的研究.文獻[15]中提出了一種提高排序效率的排序算法,使得在理論上凸包算法的時間復雜度可以降低為排序算法的時間復雜度.文獻[1]提出將二維快速凸包算法與廣義超越算法相結合,當輸入包含非極值點時,算法運行更快,并且占用內(nèi)存更少,但是當使用浮點算法時,可能會導致嚴重的錯誤.文獻[16?17]提出了基于GPU 加速的方式來提高凸包問題求解的效率.文獻[10]提出了一種計算平面自由曲面精確凸包的高效實時算法,該算法建立在圓弧構造的近似凸包上,通過對圓弧進行簡單的幾何檢驗,來確定近似的凸殼線段.文獻[18]描述了平面點集凸包的四種空間效率算法,這些算法的輸出與輸入位于同一位置,并且只使用少量的額外內(nèi)存.文獻[19]提出了使用主成分分析法對點集預處理來改善凸包算法效率的遞歸算法,但計算過程中進行多次坐標變換,使得該算法較為復雜.文獻[20]利用二維平面空間象限的幾何和對稱特性提出對稱凸包算法,對點集數(shù)較大的計算效果有明顯的改善.

通過對平面點集所在區(qū)域進行正交化分割,可以得到一種基于點集所在區(qū)域正交化分割的新算法,該算法使得分割得到的點集子集的凸包極點生成過程可以在時間上同步進行.對于大規(guī)模的點集數(shù)據(jù),分割的層次越深,并行化的程度就越高,運算的時間花銷就會越低.段落內(nèi)容安排如下: 第1 節(jié)詳細介紹了該算法的原理,包括點集區(qū)域正交化分割的規(guī)則和具體過程,以及單個子集中凸包極點的生成步驟.第2 節(jié)在特殊情況和平均情況下分析了凸包極點生成算法的時間復雜度.第3 節(jié)給出一些實驗數(shù)據(jù),比較了該算法與其他一些凸包算法的運行效率,還測試了大規(guī)模平面點集數(shù)據(jù)計算凸包的時間花銷.第4 節(jié)對全文進行了總結,并提出進一步的研究方向.

1 算法原理

本文算法主要包括點集數(shù)據(jù)的正交化分割和凸包極點的生成2 個部分.1)對平面點集進行正交化分割,以獲取不相干的點集子集簇.實驗中,在點集數(shù)據(jù)的數(shù)量為500 萬時,將正交化分割的層級設置為5.2)對預處理后的點集子集序列進行操作.首先使用初始凸包極點在序列中設置區(qū)間,根據(jù)一定的判斷準則,檢索序列子區(qū)間中的凸包極點,并且在檢索過程中拋棄掉對于凸包極點生成沒有貢獻的冗余點,而檢索得到的凸包極點又重新應用于區(qū)間設置,直到所有的凸包極點都被檢索到.

1.1 點集數(shù)據(jù)正交化分割

對于一個有限的二維點集,由平面幾何知識可知,點集中的部分點可以構成一個凸包(凸多邊形)將其余所有點包含在內(nèi).一般情況下,確定一個點是否為凸包極點,需要保證該點相對于點集中其余所有點,該點作為凸包極點是成立的.然而這樣確定凸包極點的過程是極其冗余的,為此提出了對點集數(shù)據(jù)正交化分割的概念: 將二維點集沿豎直方向分割成一系列互不相關的子集,每個子集中包含部分凸包極點,在驗證子集中的某個點是否為凸包極點時,僅需保證該點相對于該子集中的所有點其作為凸包極點成立.點集正交化分割后,一個大的凸包求解問題就變成了一簇互不影響的小問題求解.

對點集數(shù)據(jù)的正交化分割,可以保證在生成凸包極點的過程中點集數(shù)據(jù)不被重復使用,避免計算冗余.不同層級的分割過程中,根據(jù)已知的凸包極點,拋棄對于凸包極點生成沒有貢獻的無效點,可以減少數(shù)據(jù)冗余給算法帶來的無效操作.

平面點集P={pi(xi,yi),i=1,2,3,···,n},對P進行遍歷比較操作獲取4 個極值點:x方向上的極小值點pj和極大值點pk,y方向上的極小值點pm和極大值點pn,易知這4 個極值點均為凸包極點.這4 個極值點將點集所在平面區(qū)域分割成5 個部分,如圖1 所示,點集的4 個極值點分別位于矩形的4 條邊上,邊角上4 個直角三角形區(qū)域內(nèi)的點集數(shù)據(jù)相互獨立,而中間四邊形區(qū)域內(nèi)的點集數(shù)據(jù)是對于凸包極點生成沒有貢獻的無效點,可以在算法執(zhí)行中將其舍棄.

圖1 極值點結構圖Fig.1 Structure diagram of extreme points

判斷點集P中的點是否在由4 個極值點構成的四邊形中,將位于四邊形中的無效點從點集P中舍去,具體判斷方式是計算點pi(xi,yi) 是否滿足如下的不等式組:

式中,k1和b1是點pj和點pn所在直線的斜率和截距,k2和b2是點pn和點pk所在直線的斜率和截距,k3和b3是點pk和點pm所在直線的斜率和截距,k4和b4是點pm和點pj所在直線的斜率和截距.舍去滿足不等式組的無效點后,剩余的點通過簡單的區(qū)間極值判斷,構成新的點集子集P1、P2、P3和P4,分別對應圖1 矩形中的左上角、右上角、右下角和左下角所在直角三角形中的離散點數(shù)據(jù).

以P1為例,對其進行第2 層級的正交化分割,如圖2 所示,P1中的2 個極值點pj和pn所在直線ljn將矩形區(qū)域分成2 個部分,下半部分空白區(qū)域內(nèi)的點已經(jīng)在第1 次正交化分割的操作中舍去,而上半部分的點可以繼續(xù)分割.計算P1中的點到直線ljn的距離,獲取距離直線最遠的點pjnl,該點也是凸包極點,以pjnl的x方向上的值xjnl為分界點,將P1中的點分割成P11和P12兩個子集,P11中的點分布在圖2 左側下面的一個直角三角形區(qū)域中,其x值小于xjnl,P12中的點分布在圖2 右側上面的一個直角三角形區(qū)域中,其x值大于xjnl.由圖2 還可以看出,分割過程中點pj、pn和pjnl圍成的鈍角三角形區(qū)域內(nèi)的無效點也要被舍去,這個過程與上文中判斷點是否在四邊形中基本相同,只需修改部分不等式即可實現(xiàn).

圖2 點集子集第二層級正交化分割圖Fig.2 Diagram of the second level orthogonalization of the subset of point sets

得到P11和P12這個層級的全部點集子集后,再對這些點集子集進行下一個層級的正交化分割.P11中對應的兩極值點為pj和pjnl,P12中對應的兩極值點為pjnl和pn,其分割的步驟與對和P1同一層級上的點集的操作過程完全相同.

點集P正交化分割操作全部完成后,得到最小的點集子集簇對點集子集進行凸包極點求解前,需要對其中的點進行x方向上的排序,排序后得到子集序列{Si,i=1,2,3,···,m}.對于序列中x值相同的多個點,由點集P1或P3得到的子集按y值從小到大排序,由P2或P4得到的子集按y值從大到小排序,使得每個子集中的首尾兩點均為凸包極點.在大量排序算法的研究中,許多非基于比較的排序方式已經(jīng)實現(xiàn)了具有線性時間復雜度的排序算法,如計數(shù)排序,桶排序等,實驗中采用桶排序?qū)c集進行排序.

點集P所在平面區(qū)域正交化分割的流程如圖3所示.圖3 中僅給出前兩個層級正交化分割的過程,之后層級間的發(fā)展與第2 層的分割形式完全相同.由圖3 可以看出,子集的分割操作B 僅依賴于子集內(nèi)部的參數(shù),即子集中直接繼承得到的2 個極值點和距離這兩極值點所在直線最遠的點,與其他子集互不相關,所以點集正交化分割的過程是并行分割過程.在使用多核CPU 并行處理的實現(xiàn)中,首先在主進程中對點集P進行A 操作,分割出P1、P2、P3、P4,然后在主進程下開辟出4 個線程,將P1、P2、P3、P4同時放入到這4 個線程中.在每一個線程中,對對應的子集執(zhí)行B 操作,每個子集被分割成2 個次子集,當任意一個線程結束時,在該線程下再開辟出兩個線程并將兩個次子集加入進去.之后執(zhí)行相同的操作,直到某條進程中新開辟線程的次數(shù)達到預設值,預設值即為正交化分割的層數(shù).在得到最次的子集后,對其執(zhí)行C 操作獲取最小凸包極點,C 操作即為第2 節(jié)中凸包極點的生成操作.在圖3中,如果P1對應的線程中操作已經(jīng)完成,而P2對應的線程中操作還在進行中,P1對應線程下會立即開辟新的線程進行對P11和P12的操作,所以對于一個點集P,正交化分割過程花費的時間取決于并行化流程所有進程中花費時間最多的一條進程.

圖3 算法并行化流程圖Fig.3 Parallel flow chart of the algorithm

1.2 凸包極點的生成

算法1.單個子集序列的凸包極點生成算法

預處理.通過區(qū)域正交化分割獲取子集序列Si.

初始化.初始化序列Si.

1)遍歷起點設置為序列首點,遍歷終點設置為序列尾點,當前遍歷點設置為遍歷起點,最大距離點設置為遍歷終點;

2)最大距離值設置為0;

3)設置一個容器,存儲每次計算新得到的遍歷終點,容器最底部存入1)中的遍歷終點;

遍歷終點計算.利用循環(huán)迭代的方式依次計算點集中的凸包極點;

4)獲取遍歷起點與遍歷終點,計算兩點所在直線信息,當前遍歷點設置為遍歷起點,并將最大距離更新為0;

5)計算當前遍歷點到直線的距離.若為負,將其從序列中拋棄;若為正,比較其與最大距離的值;若其大于最大距離,則將該值作為最大距離的值,并更新該位置作為最大距離點;

6)當前遍歷點沿序列向下一點移動,若當前遍歷點等于遍歷終點,則進入步驟7);若不等于遍歷終點,則返回步驟5);

7)若最大距離值大于0,則將步驟5)最后記錄的最大距離點位置作為遍歷終點并進入步驟4);若小于或等于0,則將此時的遍歷終點作為遍歷起點,并將容器最外面的一個遍歷終點拋棄掉.若此時容器不為空,則將其最外面的點作為遍歷終點并進入步驟4);若為空,則結束.

8)輸出序列中剩余的點,即為凸包極點.

算法1 中點到直線距離的計算,由于計算歐氏距離需要涉及到開方等復雜運算,所以該算法使用點沿y軸方向到直線的距離?y.設目標直線為y=kx+b,點(x0,y0)到直線在垂直方向的距離為?y,這個距離乘以常值可以得到實際的歐氏距離 ?L. 由點集子集P1、P2、P3、P4分別得到的最小子集序列Si在計算距離時有略微的不同,其與凸包極點和直線的相對位置有關,由圖1可以看出,子集P1和P2分割得到的最小子集序列中凸包極點存在于直線的上方,所以?y=y0?(kx0+b),而子集P3和P4得到的最小子集序列中凸包極點存在于直線的下方,所以 ?y=(kx0+b)?y0.

最小點集子集序列的獲得是點集所在區(qū)域正交化分割得到的結果,其凸包極點的生成相互獨立,且在實現(xiàn)上步驟相同,僅有參數(shù)的差異,因而使用并行算法對其處理,可以極大地降低時間花銷,提高運行效率.各個子集序列內(nèi)的凸包極點均得到后,加以簡單合并操作就可以得到完整的凸包點集,而且這個點集是有序的.

2 時間復雜度分析

對于任意一個點集子集,由于不同的點集在平面上的分布均不相同,凸包極點在點集中的分布也是隨機的,所以由算法本身很難計算出具有一般意義的時間復雜度公式,但是當假定極點分布為一些極端情況時可以得到時間復雜度的上下界.如圖4點集遍歷示意圖,每一條線段表示一個遍歷區(qū)間,半圓曲線表示遍歷區(qū)間上的遍歷過程,線段上的點表示凸包極點.點集序列被完整遍歷的次數(shù)κ與點集序列中點的個數(shù)n的乘積 O (κ·n) 即為本算法的時間復雜度,下面計算κ的取值范圍: 假設凸包極點的個數(shù)為h,當凸包極點集中在遍歷終點且遍歷過程見圖4 的左圖,點集序列被完整遍歷的次數(shù)κ=h ?1最大;當凸包極點集中在遍歷起點且遍歷過程見圖4 的右圖,點集序列被完整遍歷的次數(shù)κ=2最小,所以該算法的時間復雜度上下限為O(2n)≤O(κ·n)≤O((h ?1)·n).

圖4 點集遍歷示意圖Fig.4 Diagram of point sets traversal

實際上只有很特殊的平面點集的凸包才會出現(xiàn)最好或最壞的情況,因此還需要研究具有近似平均情況的時間復雜度.如圖5 左圖,假設離散點在其分布區(qū)域內(nèi)是均勻分布的,且凸包極點在離散點中是均勻?qū)ΨQ存在的,則此時凸包極點的個數(shù)為2λ+1,即任意兩個凸包極點之間都會有或者都沒有一個新的凸包極點,λ是任意非負整數(shù).圖5 右圖是極點個數(shù)h=5 時的遍歷過程,圖上的數(shù)字對應遍歷過程的順序,設其序列長度為4 時,按照遍歷順序,其需要遍歷的區(qū)間長度依次為4211211,這種遍歷過程是最為對稱的形式,所以可以得到平均情況下的最好時間復雜度.對于一個有序的點集序列,通過歸納法可以得到時間復雜度的具體表達式:

圖5 理想凸包極點分布圖(h=2,3,5,9)和h=5 的遍歷過程Fig.5 Ideal convex hull pole map (h=2,3,5,9) and a traversal process with h=5

1)當點集中有2 個凸包極點,設其序列長度為1,則需遍歷的區(qū)間總長度為1;

2)當點集中有3 個凸包極點,設其序列長度為2,則需遍歷的區(qū)間總長度為211;

3)當點集中有5 個凸包極點,設其序列長度為4,則需遍歷的區(qū)間總長度為4211211;

4)當點集中有9 個凸包極點,設其序列長度為8,則需遍歷的區(qū)間總長度為842112114211211;

5)當點集中有2λ+1個凸包極點,設其序列長度為2λ,則需遍歷的區(qū)間長度和為=2λ+

令 2λ+1=h,可知,極點個數(shù)為h時的時間復雜度為=O((1+log2(h?1)·n).

凸包求解常用的Graham 算法,分治算法和快包法的時間復雜度為 O (nlogn),Jarvis 步進算法的時間復雜度為 O (h·n),通過分析不同的時間復雜度函數(shù)可以發(fā)現(xiàn),基于區(qū)域正交化分割算法的時間復雜度在最好和最壞的情況下均優(yōu)于上述凸包算法的時間復雜度.

3 算法實驗

為了驗證算法的正確性,使用該算法設計的程序?qū)ris 數(shù)據(jù)集樣本的任意兩個特征構成的二維點數(shù)據(jù)集進行了測試.Iris 數(shù)據(jù)集的中文名是安德森鳶尾花卉數(shù)據(jù)集,該數(shù)據(jù)集包含150 個樣本,對應數(shù)據(jù)集的每行數(shù)據(jù),每行數(shù)據(jù)包含一個樣本的5 個特征: 花萼長度<1>、花萼寬度<2>、花瓣長度<3>、花瓣寬度<4>和樣本的類別信息<5>,所以Iris 數(shù)據(jù)集是一個150 行5 列的二維表.圖6 展示了該算法獲得凸包極點構成的凸多邊形,圖6(a)和圖6(d) 展示了凸包極點分布相對均勻的情形,圖6(b)和圖6(e)展示了凸包極點分布較為特殊的情形.

圖6 Iris 數(shù)據(jù)集的部分凸包圖形Fig.6 Graphs of partial convex hull of Iris data set

為獲得實驗中不同規(guī)模點集對應的最佳分割層級數(shù),對規(guī)模在5000 到1000 萬之間的點集進行了實驗.表1 展示了不同規(guī)模點集在5 層正交化分割下最小子集對應的點數(shù),其中每個規(guī)模的點集分割數(shù)據(jù)均采用30 組不同類型點集的計算結果平均值.先前的實驗中,大量的測試顯示當點集點數(shù)在500萬左右時,正交化層級設置為5 可以得到最好的實驗結果.依照500 萬規(guī)模點集5 層分割為參考,分析其第5 層與第4 層關系可知,最后一次分割得到的子集規(guī)模需要保證在1000 以下.表1中,加粗數(shù)字對應的層級數(shù)即為對應規(guī)模下點集的正交化分割層級設置參考.縱觀表1 中的數(shù)據(jù),點集的規(guī)模每增加一個數(shù)量級,對應的層級設置參考也相應的增加1,所以任意大小點集的正交化分割層級設置可以得到如下初步結論: 對于規(guī)模小于1000 的點集,不進行正交化分割,1 萬左右的點集層級設置為1,10 萬左右的點集層級設置為3,100 萬左右的點集層級設置為4,1000 萬左右的點集層級設置為5,以此類推,之后點集規(guī)模每增加1 個數(shù)量級,正交化分割層級增加1.由于這個結論是基于大量實驗的數(shù)值結論,而非理論推導,所以在具體的實驗中,層級的選擇可以適應性的 ± 1.

表1 點集不同層級分割的子集點數(shù)Table 1 The number of points in the subset corresponding to different levels of segmentation

利用計算機隨機生成的二維數(shù)據(jù)點點集,對Graham 算法、Jarvis 算法、文獻[19]的算法和基于區(qū)域正交化分割的算法計算凸包極點的運行時間進行了對比測試.為了使不同算法得到的結果具有可比性,所以基于區(qū)域正交化分割的算法和對比所使用的Graham 算法、Jarvis 步進算法以及文獻[19]的算法均是由C++編程實現(xiàn),并在VS2017 軟件上運行獲取實驗數(shù)據(jù),而且沒有使用任何的外部加速.所有的實驗均在同一平臺上運行,使用的計算機配置為Inter(R)i5-4460 CPU,3.20 GHz,4 GB內(nèi)存.這些算法分別取點數(shù)量相同但類型不同的10 組點集測試,每組運行10 次取平均值.表2 展示的是10 組點集數(shù)據(jù)運行結果的平均值,算法1 是提出的算法只對點集進行正交化分割而沒有并行處理的情況.由表2 可以看出,當二維點集數(shù)據(jù)的數(shù)量大于20 萬時,使用Graham 算法和Jarvis 算法已經(jīng)無法在1 秒內(nèi)完成凸包問題的求解,而沒有并行化處理的算法1與文獻[19]算法相比效果略差,文獻[19]算法雖然仍能夠在1 秒內(nèi)完成求解,但是與基于區(qū)域正交化分割的算法相比有一個數(shù)量級上的差距.

表2 各算法運行時間 (s)Table 2 Runtime for different algorithms (s)

為測試基于區(qū)域正交化分割的算法對于超大規(guī)模點集的計算能力,表3 給出了平面點集數(shù)據(jù)數(shù)量級在105到106之間時,該算法和文獻[19]算法處理凸包問題的運行時間.表3 中凸包數(shù)量對應的兩列數(shù)據(jù)左側為文獻[19]的數(shù)據(jù),右側為本文算法的數(shù)據(jù).在表3 數(shù)據(jù)中,點集點數(shù)為20 萬、50 萬、100 萬、200 萬和500 萬時,該算法求解凸包問題花費的時間均小于文獻[19] 算法花費的時間;在點集點數(shù)為50 萬時,文獻[19]的算法已經(jīng)無法在1 秒內(nèi)完成,而該算法只需要0.1 秒的時間;在點數(shù)為200 萬時,文獻[19]需要花費7 秒的時間,而該算法只需要0.4 秒的時間;在點數(shù)為500 萬時,文獻[19]需要18.5 秒,而本文算法可以在1 秒內(nèi)完成.

表3 不同數(shù)量級點集的運行時間 (s)Table 3 Runtime comparison for different orders of magnitude point sets (s)

為進一步驗證正交化分割凸包算法的有效性,將其應用到CT 肺圖像的分割實驗中.在近幾年的研究中,對于如何從CT 圖像中準確分割出肺區(qū)域,已有基于分形幾何和最小凸包的肺區(qū)域分割方法[21],該方法相較于傳統(tǒng)方法更加準確且效率更高.文獻[21]方法的步驟是: 首先,使用閾值法獲得肺的初始區(qū)域,并依據(jù)CT 圖像上下層結構相似的性質(zhì)以及肋骨和骨組織的位置關系移除其中的氣管、支氣管、背景區(qū)域等非肺組織.然后,利用網(wǎng)格線將肺區(qū)域分割成互不重疊的小區(qū)域塊,計算各個含邊界區(qū)塊的分形維數(shù)并根據(jù)邊界的全局和局部性質(zhì)確定分形維數(shù)閾值,通過該閾值選定需要修復的邊界區(qū)塊.最后,利用最小凸包方法修復肺邊界.

在文獻[21]中,對于由閾值選定的待修復的邊界區(qū)塊,使用Jarvis 步進算法作為最小凸包法來對其進行肺邊界的修復.實驗時對待修復的邊界區(qū)塊圖像進行雙線性二次插值操作,來增加區(qū)塊中的肺陰影的點數(shù),并使用提出的算法作為最小凸包法修復肺邊界,而其余的操作均與文獻[21]的操作相同.圖7(a)展示了一幅肺的原始CT 圖像,圖7(b)是使用最小凸包法修復肺邊界操作前的預處理圖,此時已經(jīng)移除了非肺組織并對圖像進行了二值化操作,圖中肺內(nèi)部的空洞區(qū)域是預處理過程中圖像形態(tài)學腐蝕膨脹的結果.圖7(c)是由文獻[21]修復得到的圖像,圖7(d)是提出的算法修復得到的圖像.通過比較可以發(fā)現(xiàn),在肺圖像的兩側外輪廓上2 種方法修復的結果沒有明顯的區(qū)別,但是在肺圖像中間部位的邊界上,由文獻[21]算法得到左側肺圖像仍有部分凹陷處未能修復,在右側肺的圖像中甚至有一個嚴重的人工偽影,圖7(c)中深入到肺內(nèi)部的一個細長的空白區(qū)域.而使用基于區(qū)域正交化分割的算法設計的方法修復的肺圖像并沒有嚴重的人工痕跡,而且肺的輪廓表現(xiàn)得更加光滑流暢.在時間花費上,雖然實驗過程中對待修復的邊界區(qū)塊加入了插值處理,但是花費時間仍然小于使用Jarvis 步進算法作為最小凸包法的文獻[21]所花費的時間.

圖7 肺圖像的修復Fig.7 Restoration of lung image

4 結束語

對于大規(guī)模平面點集凸包問題的求解,本文提出了基于點集所在區(qū)域正交化分割的凸包極點生成算法.算法前期對點集數(shù)據(jù)所在區(qū)域進行正交化分割,可以避免數(shù)據(jù)冗余和無效計算,后期凸包極點生成的過程中,通過不斷地剔除掉點集中的非凸包極點,可以提高程序的運行效率.當點集數(shù)據(jù)的規(guī)模在500 萬時該算法仍然可以在1 秒內(nèi)完成凸包問題的求解,這對于工程應用有著重要的意義.由大量實驗得到的正交化分割層級設置的結論,也為該算法在實際問題中的應用提供了重要的參考.一系列的數(shù)值實驗結果表明,該算法準確、高效、實用性強,并且正交化分割的層級參數(shù)設置具備很高的魯棒性.