李世林 李紅軍

（北京林業(yè)大學理學院，北京，100083）

引 言

三維物體表面重建是獲取三維空間信息建立該物體的三維模型，在3D城市、3D游戲等場景建模有著重要應(yīng)用[1]；在林業(yè)信息化測量中，利用重建的樹冠形狀，準確求解樹冠體積、表面積和占地面積等信息量將會為綠地規(guī)劃提供精確的數(shù)據(jù)支持[2]；在現(xiàn)代臨床醫(yī)學中，通過重建出病人的骨骼或者腫瘤塊，將會為醫(yī)學診斷提供輔助手段[3]；由于人類活動的日益影響，許多重要的古建筑古文物都遭受到了不同程度的破壞，對這些古建筑古文物進行數(shù)字化三維重建是文物保護的有效手段之一[4]。三維激光掃描技術(shù)的快速發(fā)展，使得人們能夠快速地獲取高精度、高密度的物體表面的三維點集。所以如何準確、快速地進行三維點集的形狀分析和表面重建成為當前急需解決的重要問題。

國內(nèi)外的專家學者已經(jīng)針對二維和三維的點集形狀進行了大量的工作。1977年，文獻[5]首次以二維平面點集凸包的推廣為基礎(chǔ)，來刻畫點集的形狀；文獻[6]給出了點集形狀嚴格的數(shù)學定義；文獻[7]以Delaunary三角化為工具，采用搜索算法來“雕刻”三維點集的單連通形狀，稱為雕刻算法；文獻[8]改進“雕刻”規(guī)則，提出了三維點集表面重建的Alpha-shape算法。這種方法的優(yōu)勢在于根據(jù)設(shè)定的參數(shù)的α值大小的不同，體現(xiàn)出點集所具有的不同空洞信息，但是處理密度不均勻的數(shù)據(jù)時，一般需要對其進行很多的人工干預。1998年文獻[9-10]先后提出了一種利用點密度估計來確定α參數(shù)的方法，避免了繁瑣的人工干預，并獲得更準確的Alpha形狀；文獻[11]提出了采用雙參數(shù)Alpha-shape算法，解決了單參數(shù)在表面重建過程中精度和完整性難以兼顧的問題。

本文在前人研究的基礎(chǔ)上，考慮到原始Alpha-shape算法其參數(shù)值的選取需要過多的人工干預和在處理不均勻點云數(shù)據(jù)的重建效果不理想這兩個問題，對參數(shù)值α的選取策略進行了改進，設(shè)計出了一種依賴于近鄰點平均距離的自適應(yīng)步長的Alpha-shape表面重建算法。

1 三維點集的Alpha形狀

文獻[6]中給出了二維平面上Alpha形狀的嚴格的數(shù)學定義，本文將二維平面上的圓推廣到三維空間里的球，從而得到三維空間的Alpha形狀的定義。

定義1對于任意充分小的正實數(shù)Alpha，所有包含點集S且半徑為1/Alpha的閉合球的交稱為點集S的Alpha包，記為positive Alpha-hull，如圖1(a)所示。

定義2對于任意的負實數(shù)Alpha，所有包含點集S且半徑為-1/Alpha的閉合球的補集的交集稱為點集S的 Alpha包，記為 negative Alpha-hull，如圖1(b)所示。

定義3對于點集S中的一點p，存在一個半徑為1/Alpha閉合球，該球包含點集S中所有的點，并且點p位于球邊界上，那么點p稱為點集S的Alpha-極端。如果3個點p，q，r位于同一個閉合球的邊界上，那么這3個點稱之為Alpha-近鄰。

定義4給定一個三維空間點集S和一個任意實數(shù)Alpha，那么會得到一個這樣的空間多面體，該多面體的頂點是Alpha-極端，邊界面是分別連接Alpha-近鄰所得到三角形，這樣的空間多面體稱為三維空間點集S的Alpha形狀。

圖1 Alpha包示意圖Fig.1 Alpha-hull sketch map

為了便于理解后續(xù)的算法改進，本文提出了k近鄰平均距離的概念。

定義5對于三維空間點集S中的任意點p，點集S中距離點p最近的k個點的距離平均值稱為點p的k近鄰平均距離。

2 Alpha-shape算法與分析

1994年，Edelsbrunner等提出了三維表面重建的Alpha-shape算法[8]。算法思想是以一個半徑固定的小球在點集上“滾動”，半徑由參數(shù)α確定，若點集不能使小球通過，那么不能使小球通過的3個點所構(gòu)成的三角形就作為點集邊界面，當小球“滾動”結(jié)束后可確定整個點集的表面。

2.1 算法流程

Alpha-shape算法的具體算法流程如下，其輸入為n×3的三維空間點集S矩陣，每一行的3維向量代表空間中的一個點。輸出為m×3的邊界三角形序號集合矩陣M，其中m代表最終邊界面所包含的邊界三角形數(shù)量，每一行的3維向量代表一個邊界三角形。設(shè)三維空間點集S={p1,p2,…,pn}。

第1步從點集S中任意選擇一點p1，將與之距離小于或等于2α的點組成新的點集S1，從S1中任意選取一組點p2p3，求出過點p1p2和p3且半徑為α的球的球心o和o′。

第2步遍歷點集S1，依次求出其他點到球心o和o′的距離集合d和d′。如果d和d′中有1個集合的距離均大于等于半徑α，則可以判斷出點p1p2和p3是邊緣輪廓點，三角形p1p2p3是邊界三角形；如果d和d′中都有小于半徑α的值時，則可以判斷不是邊緣輪廓點，停止遍歷，轉(zhuǎn)到第3步。

第3步選擇S1中的下一組點按步驟1，2進行判斷，直到S1中的全部點判斷結(jié)束。

第4步選擇S中下一個點按步驟1～3進行判斷，直到S中的全部點判斷結(jié)束，輸出邊界三角形矩陣M。

2.2 算法效果分析

以一個密度不均勻的花籃模型為例對Alpha-shape算法進行效果分析，其中花籃底部密度較大，花籃頂部“把手”部分密度較小，點集模型如圖2(a)所示。當α=1時，實驗效果圖如圖2(b)所示，由于α值過小,該算法只能計算出部分表面；加大α值，當α=2.5時，底部表面重建效果較好，但是頂部“把手”部分尚未連通，如圖2(c)所示；繼續(xù)修改α值直到頂部“把手”部分連通，此時α=3.1，最后的花籃模型表面重建效果圖如圖2(d)所示。

圖2 Alpha-shape算法效果圖Fig.2 Effect diagram of Alpha-shape algorithm

觀察表面重建效果圖后發(fā)現(xiàn)，在中間“把手”和底部交接的部分，由于過大的α值，使得表面三角形偏大，其重建效果過于粗糙，不能很好地體現(xiàn)物體原有的表面凹凸信息，如圖2(e)所示；同時，過大的α值導致在處理密度較大的底部區(qū)域時需要遍歷更多的點，從而增加了大量的程序運行時間；再者，為了得到完整的物體表面重建結(jié)果，該算法在實際重建過程中需要頻繁的人工修改α參數(shù)值。綜合以上幾點，該算法在處理如林業(yè)測量等不均勻的數(shù)據(jù)時，有適應(yīng)性不強、重建效果不理想等缺點。為了完善該算法，加大該算法的適應(yīng)面，本文從α參數(shù)值的選取規(guī)則入手，提出了一種改進的Alpha-shape算法。

3 自適應(yīng)步長的Alpha-shape算法

本文在文獻[8]提出的Alpha-shape算法的基礎(chǔ)上改進了α參數(shù)值的選取規(guī)則，核心思想是利用kd-tree計算每個點的k近鄰平均距離，以該距離作為α參數(shù)值，采用動態(tài)的α值來進行后續(xù)的操作，稱之為自適應(yīng)步長的Alpha-shape算法(Variable Alpha-shape step algorithm，VAS)。

3.1 算法流程

第1步對點集S建立kd樹，從點集S中任意選擇一點p1，計算點集S中距離點p1最近的k個點的平均距離，記為α；將距離點p1小于2α的點組成新的點集S1，從S1中任意選取一組點p2p3，求出過點p1p2和p3且半徑為α的球的球心o和o′。

第2步遍歷點集S1，依次求出其他點到球心o和o′的距離集合d和d′。如果d和d′中有一個集合的距離均大于等于半徑α，則可以判斷出點p1p2和p3是邊緣輪廓點，三角形p1p2p3是邊界三角形；如果d和d′中都有小于半徑α的值時，則可以判斷不是邊緣輪廓點，停止遍歷，轉(zhuǎn)到第3步。

第3步選擇S1中的下一組點按步驟1，2進行判斷，直到S1中的全部點判斷結(jié)束。

第4步取S中下一個點按步驟1～3進行判斷，直到S中的全部點判斷結(jié)束。

需要說明的是：在第1步中，過點p1p2和p3且半徑為α的球的球心有以下3種情況：(1)點p1p2和p3所構(gòu)成的外接圓半徑大于α，此時不存在過該3點的球，但是該情況不在算法流程中出現(xiàn)，見3.2節(jié)；(2)點p1p2和p3所構(gòu)成的外接圓半徑等于α，此時球心o和o′重合；(3)點p1p2和p3所構(gòu)成的外接圓半徑小于α，此時球心o和o′關(guān)于面p1p2p3對稱。

算法中涉及點集異常處理，主要包括如下3種情況：(1)對于3個及以上的點共線的情況，所形成的邊只記距離最遠的2個點所形成的邊；(2)對于4個及以上的點共圓的情況，所形成的邊界面以多邊形三角剖分為基礎(chǔ)，避免出現(xiàn)邊界線交叉的情況；(3)對于4個及以上的點共球的情況，所形成的邊界面以多面體三角剖分為基礎(chǔ)，避免出現(xiàn)邊界面交叉的情況。

3.2 兩種算法的效果比較

VAS算法第1步中kd-tree的k近鄰查詢的平均時間復雜度為O(logn)，第2步求距離集合的時間復雜度是O(n)，第3步遍歷的時間復雜度是O(n2)，第4步遍歷的時間復雜度是O(n3)，所以該算法總的時間復雜度為O(n3)。

對于點集密度不均勻的情況下，原始的Alpha-shape算法由于其固定的α值，在處理點集的稀疏部分時，過小的α值將不能反應(yīng)點集的形狀信息，過大的α值將導致算法第1步在處理點集密集區(qū)域時，所需要遍歷的點數(shù)很多，直接加大了程序運行負擔。而VAS算法，通過動態(tài)改變α值，在保證形狀信息的前提下程序運行時間更短。同時，由于其改變α值，所以針對不均勻三維點集的表面重建效果，VAS算法更能反應(yīng)原始數(shù)據(jù)的表面信息。圖3是原始算法和VAS算法針對密度不均勻的花瓶點集的表面重建效果圖，原始算法的表面三角形偏大，其重建效果過于“粗糙”，如圖3(a)所示。而VAS算法的表面三角形較小，其表面重建效果更“細膩”，能更準確地重建出點集表面的凹凸信息，如圖3(b)所示。通過圖4的對比，改進后的VAS算法的表面重建效果要優(yōu)于原始的Alpha-shape算法。

圖3 兩種算法效果對比Fig.3 Comparison of two algorithms

4 實驗驗證

實驗針對Alpha-shape算法和VAS算法的時間效率進行對比，實驗數(shù)據(jù)分為兩類：計算機隨機生成點集(簡稱為隨機點集)和現(xiàn)實對象采集的三維點集(簡稱為采樣點集)。

4.1 隨機點集的生成

實驗所用的隨機點集分為球形點集和矩形點集，每類又分為均勻輪廓點集、非均勻輪廓點集、均勻內(nèi)部點集和非均勻內(nèi)部點集，共8組數(shù)據(jù)，每組的點數(shù)為1 000個。

球形隨機點集的解析式為

輪廓點不均勻的點集分布情況為：上下距離軸線π/4區(qū)域內(nèi)的點數(shù)是中間π/2區(qū)域內(nèi)點數(shù)的3倍，如圖4(a)所示；包含內(nèi)部點的不均勻點集分布情況為：下半球區(qū)域內(nèi)的點數(shù)是上半球區(qū)域內(nèi)的點數(shù)的3倍，如圖4(b)所示；而均勻點集是服從均勻分布的隨機生成點。

矩形隨機點集的解析式為

輪廓點不均勻的點集分布情況為：上、后、左三面的點數(shù)是下、前、右三面的點數(shù)的3倍，如圖5(a)所示；包含內(nèi)部點的不均勻點集分布情況為：上半?yún)^(qū)域內(nèi)的點數(shù)是下半?yún)^(qū)域內(nèi)的點數(shù)的3倍，如圖5(b)所示；而均勻點集是服從均勻分布的隨機生成點。

圖4 球形區(qū)域不均勻隨機點集Fig.4 Non-uniform random points in spherical region

圖5 矩形區(qū)域不均勻隨機點集Fig.5 Non-uniform random points in rectangle region

4.2 參數(shù)k的選取

關(guān)于k-近鄰平均距離中參數(shù)k的確定可以參照文獻[12]的確定方法。而文獻[13]認為近鄰點k通常在9～25之間選擇一個數(shù)，即可滿足實驗基本要求。

由于k值的大小將直接影響局部參數(shù)α的大小，如果k值過小，相對應(yīng)的局部α值將過小，此時重建并不完整，如圖2(c)所示；如果k值過大，相應(yīng)的局部α值將過大，此時程序的運行時間代價將會加大，同時重建效果并不理想，如圖2(e)所示。本文為了確定合適的參數(shù)k，對文獻[13]中提到的k值的選取范圍一一進行了重建效果實驗，實驗數(shù)據(jù)是4.1節(jié)的球形非均勻輪廓點集、矩形非均勻輪廓點集和2.2節(jié)的花籃模型。具體程序運行時間以及重建效果對比如圖6所示。從圖6可以看出，直到當k≥20時，該算法針對3種不同點集的重建效果才算完整，同時由于隨著k的取值變大，其對應(yīng)的程序運行時間將會增長，并且其增長速度在不斷增大。綜上所述，本文最終選取k=20，在實際操作過程中取得了比較好的實驗效果。

圖6 不同k值的程序運行時間及重建情況Fig.6 Operation time and reconstruction of different k values

4.3 實驗結(jié)果

上述算法的實驗是在臺式電腦上進行的，計算機配置是Intel(R)Core(TM)i7-4790 CPU@3.60 GHz處理器和4 GB內(nèi)存，程序運行環(huán)境MATLAB R2016b。表1是兩種算法對于不同點集的平均運行時間以及相應(yīng)的α參數(shù)值對比，運行時間是重復實驗20次所求的平均時間。表1中，R/S為矩形/球形;U/N為均勻/不均勻;I/O為內(nèi)部點/輪廓點。 例如：RUIpts代表矩形均勻包含內(nèi)部的點集。Alpha-shape算法的α參數(shù)值是滿足表面重建效果的最小值，其中α參數(shù)值的改變步長為0.05；而VSA算法的α參數(shù)值是k近鄰平均距離求出來的最小α值和最大α值。

表1 兩種算法針對不同的隨機點集的平均運行時間Tab.1 Average operation time of two algorithms for different random point sets

算法的實驗運行效果圖如圖7所示。

圖7 隨機點集的實驗效果圖Fig.7 Experimental effect diagram for random point sets

4.4 結(jié)果分析

從表1中的平均運行時間可以看出，本文提出的VAS算法在處理包含內(nèi)部點的點集時平均運行時間能提高30%以上。特別地，在處理密度不均勻的點集時，平均運行時間能提高60%以上，這是因為在處理密度較大的區(qū)域時，采用較小的α值，減少了需要遍歷點，從而提高了程序運行時間。但是在處理只有輪廓點的點集時，特別是均勻的點集，本文提出的算法在運行時間上多花費了近10%，這是因為輪廓點全是邊界點，減少遍歷點并不能有效地減少程序運行時間，而本文的算法在處理α值時需要消耗一定的時間，所以整體上程序的運行時間要長一些。但是對于點集密度不均勻的情況，本文提出的算法在時間效率上能提高20%以上。

4.5 激光掃描數(shù)據(jù)

本文所用的采樣點集是兩組三維激光掃描的數(shù)據(jù)，分別是包含樹枝的不均勻的樹體數(shù)據(jù)和密度比較均勻的樹冠數(shù)據(jù)，兩種算法的平均運行時間對比如表2所示，實驗運行效果圖如圖8所示。實驗表明，本文提出的VAS算法在實際采樣點集時，針對均勻和不均勻的情況下都能提高50%以上的時間運行效率。

表2 兩種算法針對采樣點集的運行時間Tab.2 Operation time of two algorithms for sampling points

5 結(jié)束語

對于三維空間離散點集的表面重建問題，本文先介紹了三維空間離散點集的Alpha形狀的相關(guān)概念。在分析表面重建的Alpha-shape算法的基礎(chǔ)上，本文提出一種改進算法——VAS算法。該算法改進了原始算法針對點集密度不均勻的情況下需要過多的人工干預α參數(shù)值的缺點，采用kd-tree和k近鄰平均距離來動態(tài)更新α值，這就使得該算法在處理點集密度較大的區(qū)域時也能以較少的遍歷次數(shù)進行表面重建，從而為三維空間離散點集的表面重建節(jié)省了大量的時間。通過隨機數(shù)據(jù)和現(xiàn)實三維模型采樣數(shù)據(jù)等多組實驗的反復驗證，實驗結(jié)果表明本文提出的改進算法與原始算法相比，除了輪廓點隨機數(shù)據(jù)的運行時間比原始算法的時間長，其余數(shù)據(jù)的時間效率都有大幅度的提高。

圖8 采樣點集的實驗效果圖Fig.8 Experimental effect diagram for sampling points

未來的工作之一是利用本文提出的算法把2015年家里加州大學圣地亞哥分校(University of California)Ery Arias-Castro教授等的周長計算工作[14]推廣到三維形狀的表面積和體積計算；把新算法的體積計算結(jié)果與文獻[15-16]的結(jié)果進行比較；另一項工作是進一步改進表面重建算法效率，采用文獻[17]所提出的基于密度的快速聚類方法，將點云數(shù)據(jù)進行聚類，對每一類固定其參數(shù)α值，從而縮短確定局部參數(shù)α值的時間，這樣在保證重建效果的基礎(chǔ)上可以進一步提高算法的時間效率。