吳光榮， 徐昌榮， 許寧， 張展鵬

（江西理工大學(xué)建筑與測繪工程學(xué)院，江西 贛州341000）

0 引 言

特征是幾何模型的重要組成部分，其對于幾何模型的外觀及其準(zhǔn)確表達(dá)具有重要作用[1].Ohtake等[2]認(rèn)為骨脊線是點(diǎn)云數(shù)據(jù)中彎曲變化程度最凸顯的區(qū)域，特征點(diǎn)與骨脊線有著公共的定義.

點(diǎn)云特征提取是目標(biāo)識別、場景重建、3D打印等各項(xiàng)任務(wù)的基礎(chǔ)，同時在數(shù)字城市建設(shè)、文物保護(hù)、三維人體模型等鄰域具有非常重要的使用價值.特征提取是一項(xiàng)重要的數(shù)據(jù)處理過程，它一直是逆向工程、模型構(gòu)建方面的研究熱點(diǎn).

Woo等[3]通過對比表面或曲線擬合提取邊緣點(diǎn)的方法，提出法矢或曲率突變的點(diǎn)為特征點(diǎn).Alrashdan等[4]將法矢向量的突變點(diǎn)作為邊界點(diǎn)，并利用神經(jīng)網(wǎng)絡(luò)來提取特征點(diǎn).Merigot等[5]學(xué)者使用Voronoi單元格的卷積協(xié)方差矩陣來描述點(diǎn)云的法矢信息，此方法雖具有一定的魯棒性，但在處理大量點(diǎn)云數(shù)據(jù)時，由于算法的運(yùn)算量呈指數(shù)增加使得運(yùn)算速度十分緩慢.黃承亮等[6]將法矢信息作為數(shù)據(jù)分割的依據(jù)，首先將點(diǎn)云數(shù)據(jù)按一定的規(guī)則連接，然后建立TIN網(wǎng)格，通過計算三角網(wǎng)格平面的法向量代替點(diǎn)的法矢，此方法較依賴數(shù)據(jù)預(yù)處理，若是將噪點(diǎn)帶入運(yùn)算，則對法矢的計算會造成較大影響.曹詩卉等[7]通過最小二乘擬合點(diǎn)的微切平面，然后估算出該點(diǎn)的法矢，此方法在一定程度上簡化了平面擬合的復(fù)雜度，但法矢的計算同樣受噪點(diǎn)的影響較大.張量等[8]提出一種基于法矢信息的點(diǎn)云數(shù)據(jù)特征分離方法，利用散亂點(diǎn)數(shù)據(jù)點(diǎn)的領(lǐng)域通過最小二乘平面擬合，對局部點(diǎn)的法矢進(jìn)行估算，該算法分割的可行性取決于領(lǐng)域點(diǎn)的判定以及平面擬合的準(zhǔn)確性.經(jīng)典的利用法矢提取特征點(diǎn)算法受噪聲點(diǎn)的影響較大，抗噪性較弱，本文提出一種改進(jìn)的利用法矢提取特征點(diǎn)方法，該方法通過細(xì)化鄰域選取，以及在曲面擬合時多次剔除誤差，從而達(dá)到準(zhǔn)確提取特征點(diǎn)的效果.

1 鄰域搜索

對于給定的點(diǎn)云數(shù)據(jù)，鄰域搜索就是為了尋求該點(diǎn)周圍的近鄰.現(xiàn)有的許多方法是利用分割來實(shí)現(xiàn)鄰域搜索，比如構(gòu)建Voronoi多邊形，但這種方法在近鄰計算時，需要遍歷所有點(diǎn)云數(shù)據(jù)，并且算法運(yùn)算量隨著分割次數(shù)的增加呈指數(shù)增加.為了避免傳統(tǒng)的對所有點(diǎn)進(jìn)行鄰域搜索，文中先引入八叉樹算法快速進(jìn)行點(diǎn)云劃分，搜索時以劃分后的格網(wǎng)為單位，該方法減少搜索點(diǎn)的數(shù)量、加快搜索速度并為之后的鄰域構(gòu)建創(chuàng)建基礎(chǔ).

1.1 八叉樹分割

八叉樹是依據(jù)點(diǎn)云的三維空間分布特點(diǎn)，分別在立方體垂直中面遞歸的八個象限切開立方體，得到八個體積相同的立方體.因此每個節(jié)點(diǎn)有八個子節(jié)點(diǎn)，被分割的節(jié)點(diǎn)稱為父節(jié)點(diǎn)，得到節(jié)點(diǎn)的稱為子節(jié)點(diǎn).分割完后判斷立方體內(nèi)點(diǎn)集的個數(shù)，可設(shè)置一個閾值，若小于則停止分割.最終得到的子節(jié)點(diǎn)稱為葉節(jié)點(diǎn)，每一個葉節(jié)點(diǎn)對應(yīng)一個葉節(jié)點(diǎn)編碼，若在分割時發(fā)現(xiàn)葉節(jié)點(diǎn)為空集，則停止分割.單元格在未達(dá)到最小尺寸或者最低復(fù)雜度時，需要對其進(jìn)行不斷的分割.八叉樹分割及節(jié)點(diǎn)編碼如圖1所示.

圖1 八叉樹分割層次結(jié)構(gòu)

1）計算初始包圍盒.算法首先根據(jù)各點(diǎn)坐標(biāo)系的大?。╔min，Xmin，Zmin）、（Xmax，Xmax，Zmax）確定初始包圍盒在 X，Y，Z 上的棱長 （Lx，Ly，Lz）. 在實(shí)際應(yīng)用中，由于點(diǎn)坐標(biāo)存在誤差，需要對現(xiàn)有坐標(biāo)進(jìn)行一定的調(diào)整， 這里設(shè)實(shí)際最小坐標(biāo)（X1，X1，Z1）小于給出的最小坐標(biāo)，加入坐標(biāo)調(diào)整系數(shù)τ.經(jīng)過調(diào)整后棱長的計算公式[9]如下：

由此可得六面體的體積計算公式V以及每個小立方體內(nèi)的點(diǎn)云密度n（設(shè)點(diǎn)云總數(shù)為N）：

確定初始包圍盒后，再用線性八叉樹對其進(jìn)行細(xì)分.分割第一步是以初始包圍盒為根節(jié)點(diǎn)，等份的分為八個小立方體，然后將八個小立方體分別作為下一步分割的根節(jié)點(diǎn)，依次分割下去，直到達(dá)到分割的停止條件.用線性八叉樹來記錄分割后每個數(shù)據(jù)的坐標(biāo)，并對每個葉節(jié)點(diǎn)進(jìn)行編號，由于八叉樹每次運(yùn)算是一分為八，在編碼每個葉節(jié)點(diǎn)時可用八進(jìn)制0～7中的一個數(shù)來記錄.每個葉節(jié)點(diǎn)對應(yīng)的八進(jìn)制數(shù)[10]如下：

式（4）中：qλ為八進(jìn)制碼 qλ∈[0，7]，λ∈[0，j-1].

立方體的空間位置序號（α，β，λ）可表示為：

式（5）中：j為八叉樹分割的次數(shù).

為了方便計算，可將上式改寫成二進(jìn)制：

則二進(jìn)制表示包圍盒空間位置編碼計算公式[11]為：

2）分割終止判斷.運(yùn)用八叉樹對根節(jié)點(diǎn)進(jìn)行不斷的分割計算量也不斷增大，k鄰域是一個確定極小范圍內(nèi)的區(qū)域數(shù)據(jù)，它的計算需依據(jù)點(diǎn)云模型、點(diǎn)云密度以及均勻度來判斷.算法的終止條件根據(jù)分割后每個小立方體內(nèi)數(shù)據(jù)點(diǎn)的個數(shù)n確定的.若n小于設(shè)定的閾值m，則算法停止.在確定終止條件時，m值是與k值正相關(guān)的，并且k值得確定是在一個大范圍點(diǎn)云數(shù)據(jù)中求得的，即得：m＝ωk，ω為相關(guān)系數(shù)，充分考慮算法的計算時間以及k鄰域判定效果，文獻(xiàn)[12]通過實(shí)驗(yàn)確定ω的最佳取值范圍為1.5～3.0.假設(shè)計算在分割了i次后停止，可得正六面體個數(shù)為8i及最小立方體的棱長：

1.2 鄰域構(gòu)建

八叉樹分割后，k鄰域的搜索理論上是在每個小立方體中搜索，但若該點(diǎn)在立方體邊緣，則其鄰域可能在其他立方體中.節(jié)點(diǎn)之間是根據(jù)六面體相連，由六面體的幾何知識可知，鄰域點(diǎn)可能在6個面、12條邊線及8個頂角上，因此在搜索時需要遍歷27個立方體.搜索鄰域點(diǎn)的一般思想是沿八叉樹向上搜索，找出當(dāng)前塊與鄰塊的共同祖先，這個祖先是向上搜索的第一個共同祖先，然后再向下搜索鄰塊.搜索過程與八叉樹的具體結(jié)構(gòu)有關(guān)[13].

假設(shè)需要搜索一點(diǎn)P的鄰域點(diǎn)，首先在該點(diǎn)所在的立方體內(nèi)進(jìn)行搜索，搜索最近鄰點(diǎn)數(shù)目達(dá)到時結(jié)束.然后依次搜索與該點(diǎn)面鄰、邊鄰、角鄰的26個相鄰的立方體，直至搜索的鄰點(diǎn)數(shù)目達(dá)到k時結(jié)束.經(jīng)過以上處理，依次得到散亂點(diǎn)云中所有數(shù)據(jù)的近鄰點(diǎn).

2 曲面擬合

2.1 最小二乘擬合曲面

最小二乘法是一種逼近理論，也是采樣數(shù)據(jù)進(jìn)行擬合時最常用的一種方法[14].在曲面擬合時，最小二乘法充分考慮周圍k鄰域信息，并且可控制在不同的應(yīng)用中k的大小，使散亂點(diǎn)云法矢估計更加標(biāo)準(zhǔn)和穩(wěn)健[15].其主要思想是找出一個最佳匹配函數(shù)使之與各坐標(biāo)點(diǎn)的距離平方和最小，也就是找出最小殘差S得到的參數(shù)值，所求表達(dá)式為：

最小二乘擬合曲面與擬合曲線原理相似，假設(shè)點(diǎn)云數(shù)據(jù)（xi，yi，zi）（i=1，2，…，n），首先通過 xi與 yi擬合曲面方程，再取擬合值zi與實(shí)際點(diǎn)的zi之差的平方和的最小值.設(shè)曲面方程為：f（x，y）=PT（x，y）A，A為方程系數(shù)，PT（x，y）為方程的一組基.為了求殘差的最小值，其殘差方程式可寫成：

P為一個n×m的矩陣，Z為實(shí)際點(diǎn)的zi坐標(biāo)，于是上式又可寫成：

根據(jù)式（13）可求的 A 值，再將 A 帶入 f（x，y）=P（x，y）A 即得擬合曲面方程.

2.2 正交最小二乘擬合曲面

以上最小二乘在曲面擬合時，將值作為一個因變量來處理，而x，y是以精確值帶入計算.實(shí)際情況x，y并不可能是一個精確值，在計算時，若是自變量的誤差超過一定范圍，需要考慮自變量存在的誤差.此時引入正交最小二乘法.

假設(shè)給出的點(diǎn)（xi，yi，zi）（i=1，2，…，n）中 xi，yi，zi的隨機(jī)誤差分別為 vxi，vyi，vzi，其協(xié)方差陣為：

則考慮自變量誤差的可z表示為：

點(diǎn)到曲面的距離殘差定義為：

則擬合方程為：

以上方程式是求殘差值ri2的最小值，而實(shí)質(zhì)上ri2表示的是點(diǎn)到曲面的垂直距離，即是求一個曲面方程，使各點(diǎn)到曲面的正交距離的平方和最小，此方法稱為正交最小二乘.

曲面觀測方程為：

可將觀測方程寫成：

綜合兩式，可得誤差方程式：

將觀測方程帶入誤差方程式可得：

激光點(diǎn)云數(shù)據(jù)一般認(rèn)為是等精度測量，根據(jù)正交最小二乘的計算準(zhǔn)則，vz2i）在等精度觀測下，按照誤差方程式VTV=min，采用間接平差法，可解出所求的曲面方程.

在前面鄰域構(gòu)建時，如果將點(diǎn)云中的異常值也帶入計算，則會降低k領(lǐng)域構(gòu)的可靠性，為保證擬合的效果，就必須剔除異常值.根據(jù)正交最小二乘的原理，通過距離對比來剔除異常值.具體步驟如下：

1）根據(jù)步驟求出曲面方程；

3）若 di＜2δ則保留該點(diǎn)，否則剔除；

4）利用剔除之后的數(shù)據(jù)重新計算曲面系數(shù)；

5）重復(fù)上述步驟 2）至 4），直到 di＜2δ時停止.

為了判定二次曲面擬合的效果，提取部分實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合實(shí)驗(yàn)：

正交最小二乘二次曲面擬合效果如圖2所示.

圖2 二次曲面擬合

3 提取特征點(diǎn)

在計算點(diǎn)云的法矢時，一般是計算在該點(diǎn)處切平面的法向量，通過求出該點(diǎn)處的擬合平面，再對等式兩邊求偏導(dǎo)可得出結(jié)果.文中采用法矢提取特征點(diǎn)，先通過八叉樹求出該點(diǎn)的近鄰，然后用正交最小二乘擬合曲面，將曲面的法矢代替該點(diǎn)切平面法矢[16].

經(jīng)過最小二乘處理后，將擬合的二次曲面作為估算面，由此得到一個曲面方程.特殊情況下若該點(diǎn)處的鄰域是一個平面，則該平面的一般表達(dá)形式為：

方程的系數(shù)可由Ax=0求出，其中：

a、b、c、d 可通過計算 ATA 的最小特征值 λmin以及特征向量（x1，x2，x3，x4）得出.

根據(jù)平面的表達(dá)式，得到平面法向量n（a，b，c），為避免各參數(shù)的相互獨(dú)立，可將法向量單位化：

若該點(diǎn)的鄰域不是一個平面，經(jīng)過正交最小二乘擬合得到一個二次曲面，可假定曲面方程為S（u，v）=S（u，v，h（u，v）），其中 h（u，v）=au2+buv+cv2+du+ev+f，可通計算曲面在該點(diǎn)出的偏導(dǎo)來求取曲面的法矢.

曲面在該點(diǎn)處的法矢可表示為：

曲面的法矢計算可表示為：

曲面的法矢可根據(jù)上式求出，但可能由于方向不同會出現(xiàn)兩個正、反的計算結(jié)果.此時需要進(jìn)行方向調(diào)整，使之一致.

文中在糾正法矢方向時，由于數(shù)據(jù)量較少，使用逐點(diǎn)比較法[17]，逐點(diǎn)比較法是建立在一定的假設(shè)條件下的，其原理較簡單，比較容易實(shí)現(xiàn).假設(shè)、分別為兩個相鄰曲面的法矢，方向?yàn)檎ńy(tǒng)一規(guī)定由里向外為正），在密集的點(diǎn)云中可假設(shè)物體不會出現(xiàn)90°的轉(zhuǎn)折，則1、之間的夾角必定小于 90°.可根據(jù)這一理論，在判定法矢方向時，先選定的方向?yàn)檎?，若、之間的夾角大于 90°，則認(rèn)為方向?yàn)樨?fù)，將取反可統(tǒng)一法矢方向.判定、的方向之后，再將與其周圍其他法矢對比，只要是夾角大于90°，則將改法矢方向取反，直至遍歷整個點(diǎn)云.通過此方法使點(diǎn)云法矢方向全部統(tǒng)一.

在不同的模型中，平坦區(qū)域內(nèi)的法矢變化一般不明顯，而在起伏變化較大的區(qū)域，比如角域、邊域，法矢變化較劇烈.在判斷該點(diǎn)是否為特征點(diǎn)時，首先設(shè)置一個閾值，若該點(diǎn)與鄰點(diǎn)的法矢夾角大于閾值，則提取判斷為特征點(diǎn)，若小于該閾值，則舍去.

4 實(shí)例分析

本實(shí)驗(yàn)的硬件環(huán)境為 Intel（R）Core（TM）i5-3230M CPU@2.60 GHz處理器、安裝內(nèi)存8.00 GB，軟件環(huán)境為MATLAB R2016a.實(shí)驗(yàn)源數(shù)據(jù)為電話聽筒模型，對點(diǎn)云數(shù)據(jù)先進(jìn)行去噪、拼接后，點(diǎn)云總數(shù)為281148個.點(diǎn)云模型如圖3所示.

圖3 電話聽筒模型圖

特征實(shí)驗(yàn)分別在相關(guān)系數(shù)為 1.5、1.8、2.0、2.2、2.4、2.6、2.8、3.0、3.2 提取點(diǎn)云數(shù)據(jù).部分提取結(jié)果如圖4所示.

圖4 不同相關(guān)系數(shù)提取效果

實(shí)驗(yàn)時，本文以相關(guān)系數(shù)為1.5～3.0為最佳的閾值選取范圍進(jìn)行實(shí)驗(yàn).通過多次的實(shí)驗(yàn)對比，分析不同的效果圖，在相關(guān)系數(shù)為2.6時提取出的模型輪廓完整、清晰，特征點(diǎn)提取層次分明，表明此時提取效果已接近最佳，并且精度也滿足實(shí)驗(yàn)要求.相關(guān)系數(shù)大于2.6時，鄰域的范圍過大，可能導(dǎo)致很多不太明顯又不是噪點(diǎn)的特征點(diǎn)被忽略，從而導(dǎo)致特征提取效果不太好.相關(guān)系數(shù)小于2.6時，鄰域范圍太小，很多法矢變化微弱的非特征點(diǎn)以及噪點(diǎn)被誤認(rèn)為特征點(diǎn)，同時加大了算法的運(yùn)行時間.

表1為在選取不同相關(guān)系數(shù)情況下，提取的點(diǎn)云數(shù)量以及簡化率.從表中可看出，相關(guān)系數(shù)較大時，提取點(diǎn)數(shù)較多，在相關(guān)系數(shù)為1.8時簡化率僅為63.59%，對比效果圖，模型左邊環(huán)形圈體征點(diǎn)由于過度提取而使提取效果不佳.隨著相關(guān)系數(shù)的逐漸增大，提取的特征點(diǎn)數(shù)呈下降趨勢，但由于點(diǎn)數(shù)的急劇減少，比如相關(guān)系數(shù)在3.0、3.2時，模型兩邊的圓提取不完整，并且中間處特征不明顯區(qū)域未能準(zhǔn)確提取.

表1 相關(guān)系數(shù)ω與簡化率的關(guān)系

為了證明本文算法的優(yōu)越性，分別用曲率、傳統(tǒng)法矢以及本文方法進(jìn)行提取實(shí)驗(yàn)對比.三種方法進(jìn)行提取的效果對比如圖5所示.

圖5 不同方法提取效果對比

從提取效果圖可看出，曲率提取以及傳統(tǒng)法矢提取受噪點(diǎn)的影響較大，在區(qū)域1不能準(zhǔn)確判別非特征點(diǎn)的敏感區(qū)域，提取較混亂，在中間一塊特征不是很明顯區(qū)域2本文方法較其他兩種方法提取效果較好，區(qū)域3由于特征點(diǎn)較少，三種方法都能準(zhǔn)確提取特征點(diǎn).本文算法在提取時由于鄰域的選擇以及對異常值的多次剔除提取的特征區(qū)域準(zhǔn)確、清晰，較曲率提取以及傳統(tǒng)法矢提取具有優(yōu)越性.

5 結(jié) 論

論文將改進(jìn)的法矢提取特征點(diǎn)算法應(yīng)用于電話聽筒點(diǎn)云數(shù)據(jù)進(jìn)行特征提取，得出以下幾點(diǎn)結(jié)論：

1）通過八叉樹分割后的點(diǎn)云，鄰域構(gòu)建時以格網(wǎng)為單位向所有可能的方位進(jìn)行鄰域搜索，加快鄰域搜索的速率同時增加了鄰域搜索的準(zhǔn)確性及可靠性.

2）正交最小二乘在曲面擬合時充分考慮自變量的誤差并引用距離對比剔除異常值，實(shí)現(xiàn)多次降噪的目的.

3）分別用本文方法及傳統(tǒng)方法進(jìn)行實(shí)驗(yàn)對比，通過效果圖可看出本文方法較傳統(tǒng)方法具有明顯的優(yōu)越性，提取效果更加清晰、準(zhǔn)確.

4）由于本文算法只是在MATLAB R2016a進(jìn)行少量數(shù)據(jù)的提取實(shí)驗(yàn)，對于大型建筑物的特征提取還有待研究.