霍佳欣 楊家志
(桂林理工大學(xué)信息科學(xué)與工程學(xué)院 廣西 桂林 541006)
近年來很多學(xué)者在三維點云的數(shù)據(jù)處理領(lǐng)域做出了很多的貢獻(xiàn)。三維點云數(shù)據(jù)在很多方面有廣泛的應(yīng)用,如3D重建、逆向工程、自動駕駛、醫(yī)療等[1-4]。激光雷達(dá)掃描是目前獲取點云數(shù)據(jù)的主要方式。但由于一些人為或環(huán)境因素,以及激光雷達(dá)掃描裝置本身的缺陷,通常會導(dǎo)致獲取的三維數(shù)據(jù)帶有一定的噪聲,而這些噪聲數(shù)據(jù)會給后續(xù)建模和測量中的相關(guān)處理帶來嚴(yán)重的問題[5-7]。因此,三維點云濾波是建模前的關(guān)鍵步驟。
噪聲大多數(shù)表現(xiàn)為高頻率信息,常用的大多數(shù)去噪方法基本上都是過濾掉高頻率的信息,保存低頻率的信息。但是,三維點云模型中所具有的幾何特征也是屬于高頻率信息,因此過濾除去高頻率的信息有可能會錯誤地刪除點云模型的尖銳幾何特征,從而丟失特征相關(guān)的信息。而通過視覺觀察物體時,對物體上所具有的特征比較敏感,所以三維點云模型數(shù)據(jù)的處理不但要去除噪聲,并且也要保持三維點云數(shù)據(jù)模型特征。在分析噪聲產(chǎn)生原因的基礎(chǔ)上,可以將噪聲進(jìn)行分類,分為大尺度噪聲也就是距離點云模型較遠(yuǎn)的噪聲、小尺度噪聲也就是距離模型較近的小尺度噪聲兩類,并分別進(jìn)行三維點云數(shù)據(jù)處理。
點云數(shù)據(jù)處理常用的基于數(shù)學(xué)形態(tài)學(xué)濾波的算法[8]、斜率濾波[9]、TIN漸進(jìn)加密算法[10]等雖然在去除點云數(shù)據(jù)模型主體上的小振幅噪聲方面有很好的效果,但對于一些離群的點大多數(shù)不能通過自適應(yīng)的方法去除,只能手工去除,這給點云數(shù)據(jù)去除噪聲帶來很多的困難?;谝苿幼钚《朔?MLS)的[11-12]和基于局部最優(yōu)投影(LOP)的方法在去噪效果上有很好的表現(xiàn),但是通常噪聲點云模型的過渡平滑,從而丟失幾何特征。雖然雙邊濾波算法[13]可以很好地保持邊緣特征,但是該算法在時間上消耗較大。針對以上算法處理點云時存在的問題,本文提出一種新的去噪算法,可以把噪聲進(jìn)行合理的分類,分為大尺度的噪聲類和小尺度的噪聲類。大尺度噪聲具有的特點是小而密集的點云,具有高頻、大幅值等。存在于點云主體內(nèi)部或者邊界的散亂點,可以稱之為小尺度噪聲點,分別使用統(tǒng)計學(xué)濾波和引導(dǎo)濾波進(jìn)行去噪。
現(xiàn)實世界中的三維點云數(shù)據(jù)模型是包含很多特征的包括尖銳特征,如邊界和角點等,很多三維點云去噪的算法對于點云模型的邊界信息以及角等尖銳特征往往會被忽略,這樣就會對點云的數(shù)據(jù)處理產(chǎn)生很大的影響。保持點云的尖銳特征的目的就是為了能夠重新構(gòu)造出精確的三維點云特征模型。濾波的目的是有效地平滑三維點云噪聲和消除模型中的噪聲,并保留物體表面原有的細(xì)節(jié)特征。算法的主要思想是先通過統(tǒng)計學(xué)濾波算法去除點云的大尺度噪聲,保留尖銳特征,之后再通過引導(dǎo)濾波算法進(jìn)行小尺度噪聲的去除。本文算法去噪流程如圖1所示。
圖1 點云去噪流程
1.1.1統(tǒng)計學(xué)濾波核心代碼
pcl:: StatisticalOutlierRemoval
//創(chuàng)建去噪對象
sor.setInputCloud(cloud); //設(shè)置需要進(jìn)行去噪的點云對象
sor.setMeanK(m); //設(shè)置 m 為在進(jìn)行統(tǒng)計時考慮查詢點
//鄰近點數(shù)
sor.setStddevMulThresh(1.0);
//設(shè)置距離閾值,
sor.filter(*cloud_filtered);
//執(zhí)行去噪計算并保存點到
//cloud_filtered
1.1.2引導(dǎo)濾波核心代碼
pcl:: KdTreeFLANN
//創(chuàng)建kd樹對象
kdtree.setInputCloud(cloud);
//輸入帶濾波的點云數(shù)據(jù)
kdtree.setEpsilon(epsilon);
//設(shè)置濾波半徑
Eigen:: Vector3d mean;
mean=neighbors_as_matrix.rowwise().mean();
//計算矩陣每行的均值
neighbors_as_matrix.transposeInPlace();
//臨近點矩陣轉(zhuǎn)置
Eigen:: MatrixXd centered=neighbors_as_
matrix.rowwise()-neighbors_as_matrix.colwise().mean();
//計算矩陣的中心
Eigen:: MatrixXd cov=(centered.adjoint() * centered) / double(neighbors_as_matrix.rows()-1);
//計算矩陣的協(xié)方差
Eigen:: MatrixXd e=(cov + epsilon * Eigen:: MatrixXd:: Identity(3,3));
e=e.inverse();
Eigen:: MatrixXd A=cov*e;
Eigen:: MatrixXd b=mean-A*mean;
//通過公式計算線性系數(shù)
searchPointEigenType=A*searchPointEigenType+b;
//計算濾波后的點云數(shù)據(jù)
激光掃描通常生成的點云數(shù)據(jù)集具有不同的密度,此外,激光掃描得到的點云數(shù)據(jù)會有稀疏的異常值產(chǎn)生,這是因為測量誤差,從而進(jìn)一步破壞點云的表達(dá)準(zhǔn)確性,使得局部點云特征的估計復(fù)雜化,從而導(dǎo)致錯誤的估計結(jié)果,進(jìn)而導(dǎo)致點云的高層應(yīng)用表現(xiàn)不佳。明顯離群點的特征是在空間中分布稀疏,基于離群點的特征可以認(rèn)為:三維點云數(shù)據(jù)集中每個點是包括一定信息量的,越密集地分布在某個區(qū)域,則可能表示的信息量越大。噪聲信息對于三維點云數(shù)據(jù)模型是屬于無用信息,且包含的信息量相對來說較小,所以可以忽略不計。考慮到離群點的特征,則可以定義一個密度,某處點云密度和定義密度作比較,如果小于定義的密度,即可以稱為點云無效也就是大尺度噪聲。根據(jù)大尺度噪聲的特點采用統(tǒng)計學(xué)濾波算法對此類噪聲進(jìn)行去噪處理。統(tǒng)計學(xué)濾波可以通過構(gòu)建K-D樹數(shù)據(jù)結(jié)構(gòu)來構(gòu)建KNN(k-nearest neighbor)[14],進(jìn)而有效對離散點云數(shù)據(jù)進(jìn)行劃分和管理,加快了計算速度。對于大尺度點云統(tǒng)計濾波算法流程如下所示。
Step1讀入存儲點云數(shù)據(jù)的PCD文件,讀取三維點云數(shù)據(jù)為P(p1,p2,…,pn)。
Step2對點云的數(shù)據(jù)點建立K-D樹的數(shù)據(jù)結(jié)構(gòu),通過K近鄰算法進(jìn)行K近鄰搜索。
Step3對于點云的每個點pi,定義參數(shù)K為近鄰點的個數(shù),根據(jù)K的值通過K近鄰算法建立鄰域,并計算點云中每個點與其最近的K個點的距離的平均值。
Step4在Step3的基礎(chǔ)上。計算點云集中所有點的K近鄰點平均距離μ和其標(biāo)準(zhǔn)差σ。
Step5式(1)中的α指的是準(zhǔn)差倍數(shù)閾值,為確定標(biāo)準(zhǔn)范圍D,設(shè)置參數(shù)μ表示為三維點云數(shù)據(jù)模型中所有點的K近鄰點距離的平均值,設(shè)置參數(shù)σ表示為三維點云數(shù)據(jù)模型中所有點的K近鄰點距離的標(biāo)準(zhǔn)差。一個點的距離的平均值的差值和準(zhǔn)差倍數(shù)閾值α作比較,如果低于α,則這一點就會被劃分為非噪聲點進(jìn)而保留,否則認(rèn)為該點是噪聲點進(jìn)而把這個點刪除。
D=μ+α×σ
(1)
Step6遍歷點云數(shù)據(jù)集中的所有點,刪除大尺度噪聲的點云數(shù)據(jù),保存沒有刪除的點云數(shù)據(jù)。
導(dǎo)引圖像濾波常用于圖像,是由He等[21]提出的一種時間效率高的保持特征的平滑算子。算法的基本思想是將輸出圖像作為窗口中制導(dǎo)圖像的線性變換。受引導(dǎo)圖像濾波的啟發(fā),Han等[15]提出了引導(dǎo)點云濾波。由于考慮到三維點云數(shù)據(jù)與圖像不同,圖像包含強度信息,而三維點云數(shù)據(jù)不包含強度信息,直接將引導(dǎo)圖像濾波技術(shù)的方法用到三維點云濾波不容易實現(xiàn)。所以從三維點云的位置信息考慮,將引導(dǎo)圖像濾波技術(shù)擴(kuò)展到點云濾波技術(shù)中形成一種通用的濾波器,為引導(dǎo)點云濾波。由于統(tǒng)計學(xué)濾波后的點云數(shù)據(jù)還有一部分小尺度噪聲需要去除,引導(dǎo)濾波不僅能去除小尺度噪聲,而且能保持點云邊緣。通過點云數(shù)據(jù)建立K-D樹結(jié)構(gòu)對點Pi的臨近點搜索可以采用K近鄰搜索(KNN)或者半徑搜索。
(2)
式中:N(Pi)為點pi的臨近點。
Step2假設(shè)點云的局部線性模型為:
(3)
Step3通過對式(4)中點在KNN鄰域的極小值進(jìn)行求解,因需要對平滑效果控制進(jìn)而設(shè)置參數(shù)ε為平滑參數(shù):
(4)
Step4對函數(shù)J(ai,bi)取極小值時,對ai和bi的偏導(dǎo)為零。解得:
(5)
(6)
Step5由Step4中的ai、bi,可以得到pi臨近點的濾波點云。
(7)
遍歷所有點云,得到濾波后的點云數(shù)據(jù)并保存,平滑參數(shù)取值通常是ε<1。
該方法的主要思想是利用導(dǎo)引點云的鄰域,將濾波后的輸出點推導(dǎo)為導(dǎo)引點云對應(yīng)點的線性模型。
整體去噪算法偽代碼如算法1所示。
算法1整體去噪算法
輸入:P(p1,p2,…,pn)。
輸出:P(p1,p2,…,pn)。
//輸入需要去噪的點云數(shù)據(jù)
Create KNN tree
//通過KNN來查詢近鄰點
For each piont i=1 to N
//對數(shù)據(jù)集中每個點遍歷
Build all point K-D tree
//創(chuàng)建K-D樹
Build all point KNN
End for
Statistic filter
For each point i=1 to N
Caculate point KNN average distance
//計算每個點臨近點的平均距離
Caculate point KNN standard deviation
//計算每個點的方差
End for
Caculate all point KNN average distance μ
//計算所有點的臨近點的平均距離
Caculate all point KNN standard deviation σ
//計算所有點臨近點的方差
Determine the point is a noise point according to formula(1)
//通過公式(1)判斷點是否為噪聲點
Guided filter
Denoise the point according to section 1.2
According formula(2),(3),(4),(5),(6),(7)
OutputP(p1,p2,…,pn)
End
實驗的主要數(shù)據(jù)來自斯坦福大學(xué)點云模型庫[18]。選用兔子模型,對點云數(shù)據(jù)添加隨機(jī)噪聲。計算機(jī)CPU為Intel(R)Core(TM)i5-3317UCPU@1.70 GHz,8 GB內(nèi)存; 操作系統(tǒng)為 Windows 10 64位; 編程環(huán)境為Visual Studio 2017+Point Cloud Library(PCL) 1.8.0。
帶噪聲的點云數(shù)據(jù),通過統(tǒng)計學(xué)濾波算法去噪,此算法去除離群噪聲需要選擇合適的參數(shù),圖2是不同參數(shù)下的去噪結(jié)果對比,α表示標(biāo)準(zhǔn)差倍數(shù)閾值,點云數(shù)據(jù)的不同臨近點個數(shù)通過設(shè)置參數(shù)K表示,取不同值消除點的個數(shù)。先確定參數(shù)α的值,因為如果α的值過小會導(dǎo)致部分離群點去除,如果α值過大又會使數(shù)據(jù)模型原本的數(shù)據(jù)去除??梢钥闯霎?dāng)α=1、K=20時較α=2、K=20和α=3、K=20去噪效果好,所以α=1是合適的參數(shù)。接下來確定參數(shù)K的值,α=1、K=20較α=1、K=30,α=1、K=40,α=1、K=50的去噪效果好,所以在統(tǒng)計學(xué)濾波去除大尺度噪聲參數(shù)時,合適的參數(shù)是α=1、K=20。在統(tǒng)計學(xué)濾波的基礎(chǔ)上選擇合適的引導(dǎo)濾波參數(shù),通過實驗可以得出在進(jìn)行引導(dǎo)濾波時,選擇半徑r=0.01、ε=0.1的參數(shù)值去噪效果表現(xiàn)最好,因為如果平滑參數(shù)過大會導(dǎo)致數(shù)據(jù)模型過度平滑,失去模型原有的幾何特征。
圖2 不同參數(shù)去噪對比
表1為帶噪聲的點云數(shù)據(jù)通過統(tǒng)計學(xué)濾波選擇不同的參數(shù)得出的結(jié)果,由噪點數(shù)的多少從而判斷合適的α與K的值,通過觀察當(dāng)α=1、K=20時去除噪聲點數(shù)更多。綜上所述,α=1、K=20是統(tǒng)計學(xué)濾波最適合參數(shù)取值。
表1 不同參數(shù)去噪點的個數(shù)
兔子的點云模型數(shù)據(jù)點個數(shù)為35 947,原始點云模型如圖3(a)所示,可以看出模型周圍是沒有噪聲點的,在原始點云模型上加上3 000個隨機(jī)噪聲點,加噪聲之后點云數(shù)據(jù)點的個數(shù)為38 947,如圖3(b)所示。采用統(tǒng)計學(xué)方法之后可以明顯看出點云的離群點基本上被去除,還存在部分小尺度噪聲在兔子模型內(nèi)部,如圖3(c)所示;通過體素化網(wǎng)格濾波,兔子模型周圍還存在較多大尺度噪聲,如圖3(d)所示;半徑濾波算法雖然去除了點云模型的大部分大尺度噪聲,但是兔子模型周圍還有小部分大尺度噪聲,如圖3(e)所示。本文算法兔子周圍大尺度噪聲基本上被去除,小尺度噪聲也被平滑,如圖3(e)所示,可以看出本文算法去噪效果較好,兔子模型的噪聲點基本上被去除。
為了驗證本文算法的去噪效果,在Bunny模型上與體素網(wǎng)格濾波算法[16]、半徑濾波算法[17]進(jìn)行了比較,其中體素網(wǎng)格算法和半徑濾波算法在點云庫(PCL)版本1.8.0中提供,通過C++實現(xiàn)。
通過表2算法的運行結(jié)果可知,本文算法噪聲濾除率要比體素網(wǎng)格濾波算法提升45.5百分點,比半徑濾波算法提升了6.7百分點,表明本文算法具有較好的去噪效果。
表2 算法性能比較
本文將統(tǒng)計學(xué)濾波算法和引導(dǎo)濾波方法相結(jié)合,應(yīng)用于三維點云數(shù)據(jù)的去噪。對于兩類噪聲的去除,實驗結(jié)果表明統(tǒng)計學(xué)濾波去除大尺度噪聲有很好的表現(xiàn),引導(dǎo)點云濾波在去除小尺度噪聲上面也有很好的表現(xiàn)。所以本文算法在去除大尺度噪聲的同時也去除了小尺度噪聲,本文算法可以有效提高噪聲濾除率,并能較好地保留有效原始點云數(shù)據(jù),可以應(yīng)用于激光雷達(dá)裝置采集的帶噪聲的點云數(shù)據(jù)通過噪聲分類進(jìn)行去噪。