張玉存， 李亞彬， 付獻(xiàn)斌

(1.燕山大學(xué) 電氣工程學(xué)院，河北 秦皇島 066004；2.河北環(huán)境工程學(xué)院 信息工程系，河北 秦皇島 066102)

1 引 言

隨著三維掃描技術(shù)的發(fā)展與廣泛應(yīng)用，三維掃描點(diǎn)云數(shù)據(jù)處理技術(shù)的研究也日益興起[1]。由于掃描過(guò)程中測(cè)量設(shè)備的限制和環(huán)境干擾，導(dǎo)致測(cè)得的三維點(diǎn)云數(shù)據(jù)不可避免地存在噪聲，點(diǎn)云消噪已成為三維幾何數(shù)據(jù)處理領(lǐng)域的重要研究方向[2]。

目前，眾多學(xué)者提出了多種點(diǎn)云模型去噪算法，主要有通過(guò)網(wǎng)格濾波算法、局部突出特征信息等去噪方法。其中根據(jù)網(wǎng)格濾波算法去噪的方法中以雙邊濾波光順去噪最為普遍，F(xiàn)leishman S等[3]介紹了一種各向異性的網(wǎng)格去噪算法，通過(guò)對(duì)圖像去噪的雙邊過(guò)濾算法采用它來(lái)進(jìn)行3D網(wǎng)格化，消除了網(wǎng)格中的噪聲。Zheng等[4]提出一種雙邊標(biāo)準(zhǔn)濾波器，通過(guò)空間距離和信號(hào)距離，該特征保護(hù)方案在視覺(jué)上和數(shù)字上都比以前的方法更有效果，特別是在具有尖銳特征或不規(guī)則采樣的區(qū)域；江亮亮等[5]提出一種用于網(wǎng)格去噪的自適應(yīng)雙邊濾波器，根據(jù)特征強(qiáng)度自適應(yīng)地調(diào)整雙邊濾波器的濾波參數(shù)進(jìn)一步提高了去噪性能；Lu 等[6]引入有效光滑面法線的規(guī)范雙邊濾波器和一個(gè)與過(guò)濾面法線估計(jì)頂點(diǎn)法線的有效方案，在頂點(diǎn)更新中，提出了一種新的迭代頂點(diǎn)更新算法來(lái)有效地更新頂點(diǎn)位置；崔鑫等[7]利用主元分析法估算點(diǎn)云的曲率和法向量，根據(jù)曲率特征標(biāo)識(shí)點(diǎn)云數(shù)據(jù)的特征區(qū)域，使用雙邊濾波器對(duì)點(diǎn)云模型進(jìn)行平滑去噪；Zhang等[8]建立了基于擴(kuò)散濾波一般原理的三維擴(kuò)散濾波方程，推導(dǎo)了尺度參數(shù)和最大迭代次數(shù)的理論公式，實(shí)現(xiàn)了自適應(yīng)降噪。但是由于雙邊濾波的平滑濾波因子和特征保持因子在不同特征區(qū)域的權(quán)重不能自適應(yīng)調(diào)整，所以在濾波過(guò)程中可能會(huì)出現(xiàn)過(guò)光順，因此研究者根據(jù)點(diǎn)云數(shù)據(jù)局部突出特征信息，如將模型的尖銳部分與其他部分分割，將點(diǎn)云數(shù)據(jù)分為平坦區(qū)與非平坦區(qū)等特征進(jìn)行去噪。吳祿慎等[9]根據(jù)點(diǎn)云平均曲率的局部特征權(quán)值，針對(duì)不同特征區(qū)域分別采用鄰域距離平均濾波算法和自適應(yīng)雙邊濾波算法進(jìn)行去噪濾波；Zheng等[10]根據(jù)特征類(lèi)型為多個(gè)法線分配特征點(diǎn)提取尖銳的特征結(jié)構(gòu)，根據(jù)過(guò)濾的法線更新點(diǎn)位置，可以穩(wěn)健地濾除高水平的噪聲，同時(shí)保持重要的幾何特征完整；Mattei等[11]提出了一種移動(dòng)魯棒主成分分析(mrpca)，點(diǎn)的估計(jì)位置通過(guò)局部平均來(lái)計(jì)算，通過(guò)加權(quán)最小化來(lái)保留尖銳特征，利用權(quán)重測(cè)量本地鄰域中的法向矢量之間的相似性進(jìn)行消噪；Wang等[12]提出了一種有效的方法來(lái)處理基于新型滾動(dòng)制導(dǎo)正常濾波器不同尺度的幾何特征，平滑小規(guī)模的幾何特征，同時(shí)保留大規(guī)模的功能；樊宇等[13]運(yùn)用激光掃描光刀法掃描得到點(diǎn)云數(shù)據(jù)，提出了一種基于幾何關(guān)系的三角形濾波法則，能夠較好地進(jìn)行去噪處理。

為避免對(duì)曲率較大的數(shù)據(jù)點(diǎn)的影響以保證邊界特征得以保持，本文提出基于曲率約束的點(diǎn)云分割去噪方法。利用最小二乘擬合曲面與過(guò)鄰域點(diǎn)作擬合曲面的法截面，結(jié)合相交曲線的旋轉(zhuǎn)角估算鄰域內(nèi)點(diǎn)的曲率，以鄰域點(diǎn)為根節(jié)點(diǎn)搜索下一鄰域點(diǎn)，根據(jù)曲率的約束對(duì)點(diǎn)云數(shù)據(jù)分割；對(duì)分割后點(diǎn)集合利用一種映射關(guān)系對(duì)噪聲點(diǎn)進(jìn)行線性回歸，移動(dòng)噪聲點(diǎn)到新位置，實(shí)現(xiàn)點(diǎn)模型去噪。

2 曲率約束分割算法

對(duì)于三維點(diǎn)云數(shù)據(jù)集P={Pi=(xi,yi,zi)T,i=1,2…N}，Pi為點(diǎn)云中一個(gè)采樣點(diǎn)，N為采樣點(diǎn)中心，建立坐標(biāo)軸，以數(shù)據(jù)的第一列為x軸，第二列為y軸，第三列為z軸，進(jìn)行曲面分割。

(1)

選取Pi點(diǎn)及其鄰域內(nèi)點(diǎn)，圖1為同面內(nèi)的鄰域點(diǎn)，圖2為不同曲面內(nèi)的鄰域點(diǎn)。根據(jù)Pi點(diǎn)與鄰域點(diǎn)所形成三角面片結(jié)構(gòu)，構(gòu)建不同三角面片的平面擬合，即：

y=AKx+BKy+CKz+DK

圖1 選取點(diǎn)P及其鄰域點(diǎn)Fig.1 Point P and its neighborhood point

圖2 鄰域點(diǎn)在不同的曲面Fig.2 Neighborhood points are on different surfaces

設(shè)三角面片的法向量a=(Aa，Ba，Ca)，b=(Ab，Bb，Cb)，AK、BK、CK、DK為每一個(gè)三角面片所生成平面的系數(shù)，其中K為鄰域點(diǎn)的個(gè)數(shù)。一個(gè)固定三角面片的法向量與其它法向量的夾角為θ，即：

(2)

當(dāng)在θn的矩陣中所有的角度差值在一定的閾值內(nèi)，則Pi的鄰域內(nèi)沒(méi)有兩曲面相交的點(diǎn)或兩曲面夾角相對(duì)較?。划?dāng)角度差值超過(guò)閾值時(shí)，則該鄰域內(nèi)存在曲面相交角度較大的交界點(diǎn)。

對(duì)于存在曲面相交點(diǎn)的鄰域點(diǎn)進(jìn)行曲面擬合，如圖3、圖4所示。

圖3 采樣點(diǎn)鄰域數(shù)據(jù)點(diǎn)Fig.3 Sampling point neighborhood data point

圖4 鄰域數(shù)據(jù)點(diǎn)擬合曲面Fig.4 Neighborhood data point fitting surface

設(shè)擬合曲面方程為：

f(x,y)=a1x2+a2y2+b1xy+b2yz+c1x+c2y+d

將各點(diǎn)投影至曲面上，通過(guò)微分幾何計(jì)算曲面上各點(diǎn)的高斯曲率K、平均曲率H：

(3)

(4)

(5)

曲率的差值的數(shù)據(jù)為：

(6)

取兩組數(shù)據(jù)μj=(X1，X2)做為初始聚類(lèi)中心，使得差值數(shù)據(jù)分為兩類(lèi)，則目標(biāo)函數(shù)為：

(7)

μj為分類(lèi)的簇中心，更新迭代使得簇成員到簇中心差值和為最小。

提取兩類(lèi)中最大值與最小值：

(8)

(9)

則定義采樣點(diǎn)的特性為：

(10)

當(dāng)一曲面中不存在相交點(diǎn)但在某段的曲率偏高時(shí)，不利于后期對(duì)數(shù)據(jù)的擬合處理，在曲率較高的情況下運(yùn)用區(qū)域生長(zhǎng)的方法對(duì)曲率較高的曲面進(jìn)行分割，可減小數(shù)據(jù)后期處理時(shí)的擬合偏差。

對(duì)于除去相交點(diǎn)后的點(diǎn)云數(shù)據(jù)，對(duì)Pi點(diǎn)鄰域內(nèi)的點(diǎn)進(jìn)行曲面擬合得曲面S0，將Pi點(diǎn)的鄰域點(diǎn)投影至擬合曲面，擬合曲面方程為：

f(x,y,z)=a1x2+a2y2+a3z2+b1xy+b2yz+
b3xz+c1x+c2y+c3z+d

圖5 曲面的法截面Fig.5 The normal section of a surface

(11)

記法截面的方程為g(x,y,z)=0，法截面S1與曲面S0的相交曲線為l0，則l0的方程表示為：

(12)

(13)

(14)

(15)

(16)

(17)

(18)

則法截面與曲面相交的弧線的旋轉(zhuǎn)角為：

(19)

旋轉(zhuǎn)θ值越大則曲線l0的平均曲率越小，旋轉(zhuǎn)θ值越小則曲線l0的平均曲率越大。

圖6 同一曲率范圍內(nèi)的點(diǎn)Fig.6 Points in the range of the same curvature

圖7 不同曲率范圍內(nèi)的點(diǎn)Fig.7 Points in the range of different curvature

3 曲率約束消噪方法

以點(diǎn)云數(shù)據(jù)點(diǎn)為球心，半徑為r的球體范圍內(nèi)的點(diǎn)云數(shù)據(jù)點(diǎn)數(shù)為該點(diǎn)的密度值[14]。在以點(diǎn)云數(shù)據(jù)點(diǎn)為球心，半徑為r的球體范圍內(nèi)尋找點(diǎn)云個(gè)數(shù)。當(dāng)該數(shù)據(jù)點(diǎn)為大尺度噪聲點(diǎn)或離群點(diǎn)時(shí)，以該點(diǎn)為球心半徑為r的球體內(nèi)的點(diǎn)云個(gè)數(shù)較少，即該點(diǎn)的密度值偏小；當(dāng)該數(shù)據(jù)點(diǎn)為正常值時(shí)，以該點(diǎn)為球心半徑為r的球體內(nèi)的點(diǎn)云個(gè)數(shù)較多，即密度值較大。通過(guò)設(shè)定球體內(nèi)的點(diǎn)云數(shù)的閾值范圍，可將大尺度噪聲及離群點(diǎn)去除。

對(duì)于集合MN{N=1,2,3,…}中含有噪聲的點(diǎn)云數(shù)據(jù)，其中不乏含有噪聲點(diǎn)距離真實(shí)點(diǎn)較遠(yuǎn)，偏離主體點(diǎn)云。對(duì)集合MN中的點(diǎn)Pi，i為集合中點(diǎn)云的個(gè)數(shù)，計(jì)算其到所有臨近點(diǎn)的距離值，設(shè)定閾值r，得到在閾值r的范圍的點(diǎn)云數(shù)據(jù)個(gè)數(shù)成為該點(diǎn)基于距離r的密度，記為density(Pi,r)。根據(jù)點(diǎn)云密度設(shè)定在距離為r的范圍內(nèi)含有的數(shù)據(jù)點(diǎn)個(gè)數(shù)的閾值ξ。當(dāng)Pi點(diǎn)的密度值density(Pi,r)<ξ，則判定點(diǎn)Pi為噪聲點(diǎn)，應(yīng)予剔除；當(dāng)Pi點(diǎn)的密度值density(Pi,r)≥ξ，則判定點(diǎn)Pi為合格點(diǎn)。

在集合MN中含有噪聲的點(diǎn)云數(shù)據(jù)中，建立一個(gè)三次曲面來(lái)逼近數(shù)據(jù)點(diǎn)，三次曲面方程為：

f(x,y)=a9x3+a8y3+a7x2y+a6xy2+a5x2+
a4y2+a3xy+a2x+a1y+a0

對(duì)比圖8、圖9可知鄰域點(diǎn)在經(jīng)過(guò)三次曲面擬合后，曲面在噪聲點(diǎn)處的曲率變化優(yōu)于二次曲面擬合，根據(jù)擬合的曲面，求取各點(diǎn)在曲面的投影點(diǎn)，即每一點(diǎn)距離曲面最近的點(diǎn)為該點(diǎn)的投影點(diǎn)，(x0，y0，z0)為數(shù)據(jù)點(diǎn)，(x，y，f(x，y))為曲面上的點(diǎn)，則距離d滿足下式：

d2=(x-x0)2+(y-y0)2+(f(x,y)-z0)2

圖8 二次曲面集合內(nèi)的點(diǎn)Fig.8 Points in a set of two surfaces

圖9 三次曲面集合內(nèi)的點(diǎn)Fig.9 Points in a set of three surfaces

求d的最小值，則：

(20)

解得投影點(diǎn)坐標(biāo)為：

(21)

取曲率較大的數(shù)據(jù)點(diǎn)構(gòu)成協(xié)方差矩陣：

(22)

(23)

式中f()為一種映射。對(duì)噪聲點(diǎn)光順平移，設(shè)噪聲點(diǎn)移動(dòng)距離為(Δxi，Δyi，Δzi)，則：

(24)

(25)

式中：ω為各分量方向的權(quán)重系數(shù)。

移動(dòng)值與真實(shí)值之間距離插值最小時(shí)，移動(dòng)后的點(diǎn)最接近真實(shí)點(diǎn)，設(shè)t*=(txi，tyi，tzi)，u*=(Δxi，Δyi，Δzi)，ω*=(ωxi，ωyi，ωzi)表示不同系數(shù)權(quán)值，則其目標(biāo)函數(shù)為：

(26)

求得系數(shù)：

(27)

將ω*代入式(25)得到移動(dòng)距離(Δxi，Δyi，Δzi)，光順去除噪聲點(diǎn)。當(dāng)移動(dòng)完成，對(duì)移動(dòng)后的點(diǎn)進(jìn)行二次曲面擬合，將各點(diǎn)投影至曲面，做各點(diǎn)的切線，以某一點(diǎn)為基準(zhǔn)點(diǎn)根據(jù)式(18)判斷角度θ。根據(jù)閾值η，當(dāng)所有角度滿足θ≤η，對(duì)下一數(shù)據(jù)點(diǎn)集合進(jìn)行光順處理；當(dāng)不滿足θ>η時(shí)，則對(duì)移動(dòng)后的點(diǎn)進(jìn)行三次曲面擬合，重復(fù)式(19)～(27)步驟，至所有點(diǎn)滿足條件，完成點(diǎn)云數(shù)據(jù)的光順去噪。

4 實(shí)驗(yàn)與分析

本文利用曲率的約束對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行分割處理，對(duì)每個(gè)區(qū)域進(jìn)行去噪，在去噪的同時(shí)可以保持曲率變化較大區(qū)域的細(xì)節(jié)特征；設(shè)計(jì)了一種映射關(guān)系，使得噪聲點(diǎn)利用曲率的變化回歸至真實(shí)點(diǎn)區(qū)域。將用雙邊濾波的處理結(jié)果與本文處理結(jié)果進(jìn)行了比較。

圖10中展示點(diǎn)云模型去除大尺度噪聲結(jié)果，其中，圖10(a)為含有大尺度噪聲點(diǎn)云模型，圖10(b)為去除大尺度噪聲點(diǎn)云模型的結(jié)果。通過(guò)對(duì)比圖10(a)實(shí)線矩形框和虛線矩形框與圖10(b)實(shí)線矩形框和虛線矩形框中的點(diǎn)云數(shù)據(jù)，可知通過(guò)設(shè)定點(diǎn)云密度的閾值可以將大尺度點(diǎn)云數(shù)據(jù)噪聲進(jìn)行去除，并且點(diǎn)云模型保持原來(lái)特征。

圖10 點(diǎn)云模型去除大尺度噪聲結(jié)果Fig.10 Point cloud model to remove large scale noise results

表1顯示點(diǎn)云模型經(jīng)過(guò)點(diǎn)云密度閾值去除大尺度噪聲的結(jié)果，設(shè)定密度球的半徑r為0.09 mm，原始點(diǎn)云數(shù)據(jù)為12 950個(gè)，經(jīng)點(diǎn)云密度閾值去除大尺度噪聲后的點(diǎn)云數(shù)為10 132個(gè)。去除大量大尺度噪聲與離群點(diǎn)，并使得模型的頂點(diǎn)與邊緣特征得到良好的保持。

表1 大尺度噪聲去噪結(jié)果Tab.1 The result of large scale noise de-noising

Fandisk噪聲模型在經(jīng)過(guò)雙邊濾波處理與本文算法處理后的比較，如圖11所示。在圖11(d)與圖11(e)中即圖11(b)與圖11(c)中的區(qū)域1，通過(guò)對(duì)比圖11(d)與圖11(e)可以看出在曲率變化較大的區(qū)域經(jīng)過(guò)雙邊濾波處理后仍有明顯噪聲影響，表面凸起，有明顯噪聲點(diǎn)殘留；本文算法對(duì)曲率變化較大的區(qū)域的特征保持效果較好，噪聲去除效果較好。在圖11(f)與圖11(g)中即圖11(b)與圖11(c)中的區(qū)域2，通過(guò)對(duì)比圖11(f)與圖11(g)可以看出在兩個(gè)曲面相交的鄰域范圍內(nèi)，雙邊濾波處理后的曲面有明顯噪聲存在，模型存在較多的噪聲棱角；本文處理后的曲面光順效果較好，且模型無(wú)明顯凸起棱角。在圖11(h)與圖11(i)中即圖11(b)與圖11(c)中的區(qū)域3，通過(guò)對(duì)比圖11(h)與圖11(i)可以看出，本文算法在邊緣特征保持與曲面光順?lè)矫嫘Ч^好。

圖11 Fandisk噪聲模型去噪結(jié)果比較Fig.11 Comparison of de-noising results of Fandisk noise model

圖12為Fandisk噪聲模型邊緣及頂點(diǎn)去噪結(jié)果。通過(guò)對(duì)比圖12(a)、圖12(b)中的區(qū)域1和區(qū)域2中的邊緣信息，可知雙邊濾波進(jìn)行消噪時(shí)對(duì)模型的邊緣信息保留不夠完整，存在較多噪聲棱角，邊緣信息處理粗糙。對(duì)比圖12(a)、圖12(b)中的的區(qū)域3，可知在尖銳的頂點(diǎn)處，雙邊濾波處理將頂點(diǎn)數(shù)據(jù)平順，頂點(diǎn)位置及邊緣發(fā)生變化；本文算法的去噪結(jié)果對(duì)模型的邊緣及頂點(diǎn)特征保存完整。

由圖11、12及表2可以看出：在運(yùn)用雙邊濾波進(jìn)行消噪時(shí)，選取最佳參數(shù)σc=10，σs=20；運(yùn)用本文的算法時(shí)選取密度球的半徑r為0.09 mm。

表2 雙邊濾波與本文算法的去噪結(jié)果對(duì)比Tab.2 Comparison of denoising results between bilateral filtering and our algorithm

經(jīng)過(guò)雙邊濾波的算法處理，能夠去除部分噪聲數(shù)據(jù)，但在模型的邊緣特征及曲率較大區(qū)域的特征保持方面小規(guī)模稍差；本文算法處理后，能夠有效地去除數(shù)據(jù)的冗余，同時(shí)能夠保持模型的邊緣特征及曲率較大區(qū)域的特征。

圖12 Fandisk噪聲模型邊緣及頂點(diǎn)去噪結(jié)果Fig.12 Edge and vertex denoising results of Fandisk noise models

5 結(jié) 論

本文提出了一種基于曲率約束的點(diǎn)云分割去噪方法。根據(jù)采樣點(diǎn)與其鄰域點(diǎn)在局部擬合曲面上的投影點(diǎn)之間的關(guān)系構(gòu)建曲率權(quán)重約束條件，下一采樣點(diǎn)的與其鄰域點(diǎn)在局部擬合曲面上投影點(diǎn)之間同樣采用曲率權(quán)重約束條件，判斷兩鄰域是否添加至同一點(diǎn)云數(shù)據(jù)集合中，對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行曲率分割；運(yùn)用協(xié)方差分析法對(duì)每一點(diǎn)云數(shù)據(jù)集合得出噪聲點(diǎn)的特征向量，構(gòu)建一種映射關(guān)系，實(shí)現(xiàn)點(diǎn)云數(shù)據(jù)去噪。實(shí)驗(yàn)結(jié)果表明，本文算法在確保消除點(diǎn)云噪聲的同時(shí)，能有效地保存點(diǎn)云曲面變化較大時(shí)的信息特征，提高去噪效果，防止消噪過(guò)程中點(diǎn)云模型的扭曲變形。