統(tǒng)計(jì)學(xué)濾波和引導(dǎo)濾波相結(jié)合的點(diǎn)云數(shù)據(jù)降噪

霍佳欣 楊家志

(桂林理工大學(xué)信息科學(xué)與工程學(xué)院 廣西 桂林 541006)

0 引 言

近年來很多學(xué)者在三維點(diǎn)云的數(shù)據(jù)處理領(lǐng)域做出了很多的貢獻(xiàn)。三維點(diǎn)云數(shù)據(jù)在很多方面有廣泛的應(yīng)用,如3D重建、逆向工程、自動(dòng)駕駛、醫(yī)療等[1-4]。激光雷達(dá)掃描是目前獲取點(diǎn)云數(shù)據(jù)的主要方式。但由于一些人為或環(huán)境因素,以及激光雷達(dá)掃描裝置本身的缺陷,通常會(huì)導(dǎo)致獲取的三維數(shù)據(jù)帶有一定的噪聲,而這些噪聲數(shù)據(jù)會(huì)給后續(xù)建模和測(cè)量中的相關(guān)處理帶來嚴(yán)重的問題[5-7]。因此,三維點(diǎn)云濾波是建模前的關(guān)鍵步驟。

噪聲大多數(shù)表現(xiàn)為高頻率信息,常用的大多數(shù)去噪方法基本上都是過濾掉高頻率的信息,保存低頻率的信息。但是,三維點(diǎn)云模型中所具有的幾何特征也是屬于高頻率信息,因此過濾除去高頻率的信息有可能會(huì)錯(cuò)誤地刪除點(diǎn)云模型的尖銳幾何特征,從而丟失特征相關(guān)的信息。而通過視覺觀察物體時(shí),對(duì)物體上所具有的特征比較敏感,所以三維點(diǎn)云模型數(shù)據(jù)的處理不但要去除噪聲,并且也要保持三維點(diǎn)云數(shù)據(jù)模型特征。在分析噪聲產(chǎn)生原因的基礎(chǔ)上,可以將噪聲進(jìn)行分類,分為大尺度噪聲也就是距離點(diǎn)云模型較遠(yuǎn)的噪聲、小尺度噪聲也就是距離模型較近的小尺度噪聲兩類,并分別進(jìn)行三維點(diǎn)云數(shù)據(jù)處理。

點(diǎn)云數(shù)據(jù)處理常用的基于數(shù)學(xué)形態(tài)學(xué)濾波的算法[8]、斜率濾波[9]、TIN漸進(jìn)加密算法[10]等雖然在去除點(diǎn)云數(shù)據(jù)模型主體上的小振幅噪聲方面有很好的效果,但對(duì)于一些離群的點(diǎn)大多數(shù)不能通過自適應(yīng)的方法去除,只能手工去除,這給點(diǎn)云數(shù)據(jù)去除噪聲帶來很多的困難?；谝苿?dòng)最小二乘法(MLS)的[11-12]和基于局部最優(yōu)投影(LOP)的方法在去噪效果上有很好的表現(xiàn),但是通常噪聲點(diǎn)云模型的過渡平滑,從而丟失幾何特征。雖然雙邊濾波算法[13]可以很好地保持邊緣特征,但是該算法在時(shí)間上消耗較大。針對(duì)以上算法處理點(diǎn)云時(shí)存在的問題,本文提出一種新的去噪算法,可以把噪聲進(jìn)行合理的分類,分為大尺度的噪聲類和小尺度的噪聲類。大尺度噪聲具有的特點(diǎn)是小而密集的點(diǎn)云,具有高頻、大幅值等。存在于點(diǎn)云主體內(nèi)部或者邊界的散亂點(diǎn),可以稱之為小尺度噪聲點(diǎn),分別使用統(tǒng)計(jì)學(xué)濾波和引導(dǎo)濾波進(jìn)行去噪。

1 算法描述

現(xiàn)實(shí)世界中的三維點(diǎn)云數(shù)據(jù)模型是包含很多特征的包括尖銳特征,如邊界和角點(diǎn)等,很多三維點(diǎn)云去噪的算法對(duì)于點(diǎn)云模型的邊界信息以及角等尖銳特征往往會(huì)被忽略,這樣就會(huì)對(duì)點(diǎn)云的數(shù)據(jù)處理產(chǎn)生很大的影響。保持點(diǎn)云的尖銳特征的目的就是為了能夠重新構(gòu)造出精確的三維點(diǎn)云特征模型。濾波的目的是有效地平滑三維點(diǎn)云噪聲和消除模型中的噪聲,并保留物體表面原有的細(xì)節(jié)特征。算法的主要思想是先通過統(tǒng)計(jì)學(xué)濾波算法去除點(diǎn)云的大尺度噪聲,保留尖銳特征,之后再通過引導(dǎo)濾波算法進(jìn)行小尺度噪聲的去除。本文算法去噪流程如圖1所示。

圖1 點(diǎn)云去噪流程

1.1 算法實(shí)現(xiàn)關(guān)鍵代碼說明

1.1.1統(tǒng)計(jì)學(xué)濾波核心代碼

pcl:: StatisticalOutlierRemoval sor;

//創(chuàng)建去噪對(duì)象

sor.setInputCloud(cloud); //設(shè)置需要進(jìn)行去噪的點(diǎn)云對(duì)象

sor.setMeanK(m); //設(shè)置 m 為在進(jìn)行統(tǒng)計(jì)時(shí)考慮查詢點(diǎn)

//鄰近點(diǎn)數(shù)

sor.setStddevMulThresh(1.0);

//設(shè)置距離閾值,

sor.filter(*cloud_filtered);

//執(zhí)行去噪計(jì)算并保存點(diǎn)到

//cloud_filtered

1.1.2引導(dǎo)濾波核心代碼

pcl:: KdTreeFLANN kdtree;

//創(chuàng)建kd樹對(duì)象

kdtree.setInputCloud(cloud);

//輸入帶濾波的點(diǎn)云數(shù)據(jù)

kdtree.setEpsilon(epsilon);

//設(shè)置濾波半徑

Eigen:: Vector3d mean;

mean=neighbors_as_matrix.rowwise().mean();

//計(jì)算矩陣每行的均值

neighbors_as_matrix.transposeInPlace();

//臨近點(diǎn)矩陣轉(zhuǎn)置

Eigen:: MatrixXd centered=neighbors_as_

matrix.rowwise()-neighbors_as_matrix.colwise().mean();

//計(jì)算矩陣的中心

Eigen:: MatrixXd cov=(centered.adjoint() * centered) / double(neighbors_as_matrix.rows()-1);

//計(jì)算矩陣的協(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;

//通過公式計(jì)算線性系數(shù)

searchPointEigenType=A*searchPointEigenType+b;

//計(jì)算濾波后的點(diǎn)云數(shù)據(jù)

1.2 基于統(tǒng)計(jì)學(xué)濾波去除大尺度噪聲

激光掃描通常生成的點(diǎn)云數(shù)據(jù)集具有不同的密度,此外,激光掃描得到的點(diǎn)云數(shù)據(jù)會(huì)有稀疏的異常值產(chǎn)生,這是因?yàn)闇y(cè)量誤差,從而進(jìn)一步破壞點(diǎn)云的表達(dá)準(zhǔn)確性,使得局部點(diǎn)云特征的估計(jì)復(fù)雜化,從而導(dǎo)致錯(cuò)誤的估計(jì)結(jié)果,進(jìn)而導(dǎo)致點(diǎn)云的高層應(yīng)用表現(xiàn)不佳。明顯離群點(diǎn)的特征是在空間中分布稀疏,基于離群點(diǎn)的特征可以認(rèn)為:三維點(diǎn)云數(shù)據(jù)集中每個(gè)點(diǎn)是包括一定信息量的,越密集地分布在某個(gè)區(qū)域,則可能表示的信息量越大。噪聲信息對(duì)于三維點(diǎn)云數(shù)據(jù)模型是屬于無用信息,且包含的信息量相對(duì)來說較小,所以可以忽略不計(jì)?？紤]到離群點(diǎn)的特征,則可以定義一個(gè)密度,某處點(diǎn)云密度和定義密度作比較,如果小于定義的密度,即可以稱為點(diǎn)云無效也就是大尺度噪聲。根據(jù)大尺度噪聲的特點(diǎn)采用統(tǒng)計(jì)學(xué)濾波算法對(duì)此類噪聲進(jìn)行去噪處理。統(tǒng)計(jì)學(xué)濾波可以通過構(gòu)建K-D樹數(shù)據(jù)結(jié)構(gòu)來構(gòu)建KNN(k-nearest neighbor)[14],進(jìn)而有效對(duì)離散點(diǎn)云數(shù)據(jù)進(jìn)行劃分和管理,加快了計(jì)算速度。對(duì)于大尺度點(diǎn)云統(tǒng)計(jì)濾波算法流程如下所示。

Step1讀入存儲(chǔ)點(diǎn)云數(shù)據(jù)的PCD文件,讀取三維點(diǎn)云數(shù)據(jù)為P(p1,p2,…,pn)。

Step2對(duì)點(diǎn)云的數(shù)據(jù)點(diǎn)建立K-D樹的數(shù)據(jù)結(jié)構(gòu),通過K近鄰算法進(jìn)行K近鄰搜索。

Step3對(duì)于點(diǎn)云的每個(gè)點(diǎn)pi,定義參數(shù)K為近鄰點(diǎn)的個(gè)數(shù),根據(jù)K的值通過K近鄰算法建立鄰域,并計(jì)算點(diǎn)云中每個(gè)點(diǎn)與其最近的K個(gè)點(diǎn)的距離的平均值。

Step4在Step3的基礎(chǔ)上。計(jì)算點(diǎn)云集中所有點(diǎn)的K近鄰點(diǎn)平均距離μ和其標(biāo)準(zhǔn)差σ。

Step5式(1)中的α指的是準(zhǔn)差倍數(shù)閾值,為確定標(biāo)準(zhǔn)范圍D,設(shè)置參數(shù)μ表示為三維點(diǎn)云數(shù)據(jù)模型中所有點(diǎn)的K近鄰點(diǎn)距離的平均值,設(shè)置參數(shù)σ表示為三維點(diǎn)云數(shù)據(jù)模型中所有點(diǎn)的K近鄰點(diǎn)距離的標(biāo)準(zhǔn)差。一個(gè)點(diǎn)的距離的平均值的差值和準(zhǔn)差倍數(shù)閾值α作比較,如果低于α,則這一點(diǎn)就會(huì)被劃分為非噪聲點(diǎn)進(jìn)而保留,否則認(rèn)為該點(diǎn)是噪聲點(diǎn)進(jìn)而把這個(gè)點(diǎn)刪除。

D=μ+α×σ

(1)

Step6遍歷點(diǎn)云數(shù)據(jù)集中的所有點(diǎn),刪除大尺度噪聲的點(diǎn)云數(shù)據(jù),保存沒有刪除的點(diǎn)云數(shù)據(jù)。

1.3 基于引導(dǎo)濾波去除小尺度噪聲

導(dǎo)引圖像濾波常用于圖像,是由He等[21]提出的一種時(shí)間效率高的保持特征的平滑算子。算法的基本思想是將輸出圖像作為窗口中制導(dǎo)圖像的線性變換。受引導(dǎo)圖像濾波的啟發(fā),Han等[15]提出了引導(dǎo)點(diǎn)云濾波。由于考慮到三維點(diǎn)云數(shù)據(jù)與圖像不同,圖像包含強(qiáng)度信息,而三維點(diǎn)云數(shù)據(jù)不包含強(qiáng)度信息,直接將引導(dǎo)圖像濾波技術(shù)的方法用到三維點(diǎn)云濾波不容易實(shí)現(xiàn)。所以從三維點(diǎn)云的位置信息考慮,將引導(dǎo)圖像濾波技術(shù)擴(kuò)展到點(diǎn)云濾波技術(shù)中形成一種通用的濾波器,為引導(dǎo)點(diǎn)云濾波。由于統(tǒng)計(jì)學(xué)濾波后的點(diǎn)云數(shù)據(jù)還有一部分小尺度噪聲需要去除,引導(dǎo)濾波不僅能去除小尺度噪聲,而且能保持點(diǎn)云邊緣。通過點(diǎn)云數(shù)據(jù)建立K-D樹結(jié)構(gòu)對(duì)點(diǎn)Pi的臨近點(diǎn)搜索可以采用K近鄰搜索(KNN)或者半徑搜索。

(2)

式中:N(Pi)為點(diǎn)pi的臨近點(diǎn)。

Step2假設(shè)點(diǎn)云的局部線性模型為:

(3)

Step3通過對(duì)式(4)中點(diǎn)在KNN鄰域的極小值進(jìn)行求解,因需要對(duì)平滑效果控制進(jìn)而設(shè)置參數(shù)ε為平滑參數(shù):

(4)

Step4對(duì)函數(shù)J(ai,bi)取極小值時(shí),對(duì)ai和bi的偏導(dǎo)為零。解得:

(5)

(6)

Step5由Step4中的ai、bi,可以得到pi臨近點(diǎn)的濾波點(diǎn)云。

(7)

遍歷所有點(diǎn)云,得到濾波后的點(diǎn)云數(shù)據(jù)并保存,平滑參數(shù)取值通常是ε<1。

該方法的主要思想是利用導(dǎo)引點(diǎn)云的鄰域,將濾波后的輸出點(diǎn)推導(dǎo)為導(dǎo)引點(diǎn)云對(duì)應(yīng)點(diǎn)的線性模型。

整體去噪算法偽代碼如算法1所示。

算法1整體去噪算法

輸入:P(p1,p2,…,pn)。

輸出:P(p1,p2,…,pn)。

//輸入需要去噪的點(diǎn)云數(shù)據(jù)

Create KNN tree

//通過KNN來查詢近鄰點(diǎn)

For each piont i=1 to N

//對(duì)數(shù)據(jù)集中每個(gè)點(diǎn)遍歷

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

//計(jì)算每個(gè)點(diǎn)臨近點(diǎn)的平均距離

Caculate point KNN standard deviation

//計(jì)算每個(gè)點(diǎn)的方差

End for

Caculate all point KNN average distance μ

//計(jì)算所有點(diǎn)的臨近點(diǎn)的平均距離

Caculate all point KNN standard deviation σ

//計(jì)算所有點(diǎn)臨近點(diǎn)的方差

Determine the point is a noise point according to formula(1)

//通過公式(1)判斷點(diǎn)是否為噪聲點(diǎn)

Guided filter

Denoise the point according to section 1.2

According formula(2),(3),(4),(5),(6),(7)

OutputP(p1,p2,…,pn)

End

2 實(shí) 驗(yàn)

2.1 實(shí)驗(yàn)環(huán)境

實(shí)驗(yàn)的主要數(shù)據(jù)來自斯坦福大學(xué)點(diǎn)云模型庫(kù)[18]。選用兔子模型,對(duì)點(diǎn)云數(shù)據(jù)添加隨機(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。

2.2 實(shí)驗(yàn)參數(shù)設(shè)置

帶噪聲的點(diǎn)云數(shù)據(jù),通過統(tǒng)計(jì)學(xué)濾波算法去噪,此算法去除離群噪聲需要選擇合適的參數(shù),圖2是不同參數(shù)下的去噪結(jié)果對(duì)比,α表示標(biāo)準(zhǔn)差倍數(shù)閾值,點(diǎn)云數(shù)據(jù)的不同臨近點(diǎn)個(gè)數(shù)通過設(shè)置參數(shù)K表示,取不同值消除點(diǎn)的個(gè)數(shù)。先確定參數(shù)α的值,因?yàn)槿绻恋闹颠^小會(huì)導(dǎo)致部分離群點(diǎn)去除,如果α值過大又會(huì)使數(shù)據(jù)模型原本的數(shù)據(jù)去除?？梢钥闯霎?dāng)α=1、K=20時(shí)較α=2、K=20和α=3、K=20去噪效果好,所以α=1是合適的參數(shù)。接下來確定參數(shù)K的值,α=1、K=20較α=1、K=30,α=1、K=40,α=1、K=50的去噪效果好,所以在統(tǒng)計(jì)學(xué)濾波去除大尺度噪聲參數(shù)時(shí),合適的參數(shù)是α=1、K=20。在統(tǒng)計(jì)學(xué)濾波的基礎(chǔ)上選擇合適的引導(dǎo)濾波參數(shù),通過實(shí)驗(yàn)可以得出在進(jìn)行引導(dǎo)濾波時(shí),選擇半徑r=0.01、ε=0.1的參數(shù)值去噪效果表現(xiàn)最好,因?yàn)槿绻交瑓?shù)過大會(huì)導(dǎo)致數(shù)據(jù)模型過度平滑,失去模型原有的幾何特征。

圖2 不同參數(shù)去噪對(duì)比

表1為帶噪聲的點(diǎn)云數(shù)據(jù)通過統(tǒng)計(jì)學(xué)濾波選擇不同的參數(shù)得出的結(jié)果,由噪點(diǎn)數(shù)的多少?gòu)亩袛嗪线m的α與K的值,通過觀察當(dāng)α=1、K=20時(shí)去除噪聲點(diǎn)數(shù)更多。綜上所述,α=1、K=20是統(tǒng)計(jì)學(xué)濾波最適合參數(shù)取值。

表1 不同參數(shù)去噪點(diǎn)的個(gè)數(shù)

2.3 實(shí)驗(yàn)結(jié)果及分析

兔子的點(diǎn)云模型數(shù)據(jù)點(diǎn)個(gè)數(shù)為35 947,原始點(diǎn)云模型如圖3(a)所示,可以看出模型周圍是沒有噪聲點(diǎn)的,在原始點(diǎn)云模型上加上3 000個(gè)隨機(jī)噪聲點(diǎn),加噪聲之后點(diǎn)云數(shù)據(jù)點(diǎn)的個(gè)數(shù)為38 947,如圖3(b)所示。采用統(tǒng)計(jì)學(xué)方法之后可以明顯看出點(diǎn)云的離群點(diǎn)基本上被去除,還存在部分小尺度噪聲在兔子模型內(nèi)部,如圖3(c)所示;通過體素化網(wǎng)格濾波,兔子模型周圍還存在較多大尺度噪聲,如圖3(d)所示;半徑濾波算法雖然去除了點(diǎn)云模型的大部分大尺度噪聲,但是兔子模型周圍還有小部分大尺度噪聲,如圖3(e)所示。本文算法兔子周圍大尺度噪聲基本上被去除,小尺度噪聲也被平滑,如圖3(e)所示,可以看出本文算法去噪效果較好,兔子模型的噪聲點(diǎn)基本上被去除。

為了驗(yàn)證本文算法的去噪效果,在Bunny模型上與體素網(wǎng)格濾波算法[16]、半徑濾波算法[17]進(jìn)行了比較,其中體素網(wǎng)格算法和半徑濾波算法在點(diǎn)云庫(kù)(PCL)版本1.8.0中提供,通過C++實(shí)現(xiàn)。

通過表2算法的運(yùn)行結(jié)果可知,本文算法噪聲濾除率要比體素網(wǎng)格濾波算法提升45.5百分點(diǎn),比半徑濾波算法提升了6.7百分點(diǎn),表明本文算法具有較好的去噪效果。

表2 算法性能比較

3 結(jié) 語

本文將統(tǒng)計(jì)學(xué)濾波算法和引導(dǎo)濾波方法相結(jié)合,應(yīng)用于三維點(diǎn)云數(shù)據(jù)的去噪。對(duì)于兩類噪聲的去除,實(shí)驗(yàn)結(jié)果表明統(tǒng)計(jì)學(xué)濾波去除大尺度噪聲有很好的表現(xiàn),引導(dǎo)點(diǎn)云濾波在去除小尺度噪聲上面也有很好的表現(xiàn)。所以本文算法在去除大尺度噪聲的同時(shí)也去除了小尺度噪聲,本文算法可以有效提高噪聲濾除率,并能較好地保留有效原始點(diǎn)云數(shù)據(jù),可以應(yīng)用于激光雷達(dá)裝置采集的帶噪聲的點(diǎn)云數(shù)據(jù)通過噪聲分類進(jìn)行去噪。