散亂點(diǎn)云的圓柱提取和擬合

周 穎 林 意

(江南大學(xué)人工智能與計(jì)算機(jī)學(xué)院 江蘇 無(wú)錫 214122)

0 引 言

隨著激光掃描設(shè)備在各個(gè)領(lǐng)域的應(yīng)用,人們可以用較低的成本快速獲取物體表面的三維坐標(biāo),其中對(duì)圓柱點(diǎn)云擬合的應(yīng)用是工業(yè)上比較常見(jiàn)的問(wèn)題。圓柱擬合在工業(yè)上的應(yīng)用包括特征提取[1]、對(duì)象識(shí)別[2]、自主導(dǎo)航[3]、幾何質(zhì)量檢測(cè)[4]、施工測(cè)試[5]等。目前,研究人員對(duì)圓柱擬合做了一系列的研究,解決了工業(yè)上的一些難題。常用的圓柱擬合方法歸納起來(lái)有:遺傳算法[6]、高斯圖法[7]、坐標(biāo)轉(zhuǎn)換[8]、圓度判別法[9]等。這些方法原理復(fù)雜、計(jì)算量大。對(duì)此,張益澤等[10]提出了任意初值的圓柱擬合方法,通過(guò)對(duì)傳統(tǒng)的誤差方程進(jìn)行改進(jìn),減少了各參數(shù)之間的關(guān)聯(lián),在方程中賦予任意初值,利用附有條件方程的間接平差求解圓柱參數(shù);劉支亮等[11]在文獻(xiàn)[10]的基礎(chǔ)上,提出了自適應(yīng)調(diào)整權(quán)陣的穩(wěn)健加權(quán)總體算法,以3倍準(zhǔn)差為閾值,剔除粗差,構(gòu)建穩(wěn)健的最小二乘擬合算法。迭代計(jì)算對(duì)初值的要求很高,研究人員對(duì)迭代初值做了進(jìn)一步的改進(jìn):Wu等[12]通過(guò)協(xié)方差估計(jì)點(diǎn)云的法向量,再根據(jù)最小二乘法來(lái)估計(jì)軸線方向向量;申旭等[13]提出了三點(diǎn)共線法求得直線方向向量,把直線方向向量作為初始值代入誤差方程迭代求解。此外,Paláncz等[14]利用期望最大化和最大似然估計(jì)來(lái)尋找高斯混合分布情況下的圓柱參數(shù);Al-Subaihi[15]提出了一種基于最小化正交平方距離和最小的圓柱擬合方法,即最小化點(diǎn)到圓柱的正交距離,并用梯度下降法求解;梁雪等[16]研究了一種一邊擬合一邊去噪的方法來(lái)獲取圓柱參數(shù),先隨機(jī)選取一定量的點(diǎn),根據(jù)圓柱面的方程擬合圓柱,再選取一個(gè)閾值,把每一個(gè)點(diǎn)代入方程,計(jì)算的結(jié)果與閾值相比較,若小于閾值,則為內(nèi)點(diǎn),否則為噪聲點(diǎn)。以上算法大部分都需要迭代計(jì)算,計(jì)算量大,計(jì)算時(shí)間長(zhǎng),然而迭代對(duì)初始值的要求很高,容易陷入局部最優(yōu)陷阱,對(duì)計(jì)算結(jié)果產(chǎn)生影響。

在工業(yè)測(cè)量中,受到場(chǎng)地和儀器的限制,很難采集到規(guī)則的沒(méi)有噪聲的點(diǎn)云。由于采集到的點(diǎn)云是雜亂的,不能對(duì)初始點(diǎn)云直接進(jìn)行圓柱擬合,因此,本文先對(duì)點(diǎn)云進(jìn)行預(yù)處理,采用統(tǒng)計(jì)濾波消除噪聲點(diǎn),然后用體素濾波器進(jìn)行降采樣[17],得到數(shù)目較少且沒(méi)有噪聲的點(diǎn)云。接下來(lái)通過(guò)區(qū)域生長(zhǎng)算法對(duì)點(diǎn)云進(jìn)行分割,提取分割后的圓柱點(diǎn)云,再進(jìn)行圓柱面擬合。圓柱面可以看作是一組到定直線的距離為定值的點(diǎn)的集合,這一特征表明了圓柱需要7個(gè)參數(shù)來(lái)唯一確定,這7個(gè)參數(shù)分別是圓柱中心軸線的方向向量(a,b,c)、中心軸線上一點(diǎn)(x0,y0,z0),以及圓柱面的半徑r,如圖1所示,計(jì)算點(diǎn)到直線的距離,列出誤差方程,利用奇異值分解求出軸線方向向量,再通過(guò)旋轉(zhuǎn),使得軸線與z軸平行,將點(diǎn)云平行投影到xoy平面,把三維問(wèn)題轉(zhuǎn)換為了二維平面問(wèn)題,再根據(jù)最小二乘原理擬合圓求出半徑。本文方法理論簡(jiǎn)單,易于實(shí)現(xiàn),而且避免了迭代計(jì)算。

1 點(diǎn)云預(yù)處理

1.1 統(tǒng)計(jì)濾波器

1.2 體素降采樣

密集的點(diǎn)云在進(jìn)行分割和特征提取的時(shí)候會(huì)增加計(jì)算難度和計(jì)算時(shí)間,用體素濾波器進(jìn)行降采樣能在減少點(diǎn)云數(shù)量的同時(shí)不改變點(diǎn)云的幾何結(jié)構(gòu)。體素是一個(gè)三維空間的小立方體,在進(jìn)行降采樣時(shí),體素內(nèi)所有的點(diǎn)用它們的重心來(lái)近似表示,最后,所有體素的重心組成濾波后的點(diǎn)云。如圖2所示,圖2(a)為創(chuàng)建的體積為l×l×l的包圍盒,圖2(b)為降采樣后的用體素重心表示的點(diǎn)云。體素邊長(zhǎng)l的大小控制著采樣率,當(dāng)l過(guò)小,采樣后的點(diǎn)云仍然密集,效果不明顯;當(dāng)l過(guò)大,采樣后的點(diǎn)云過(guò)于稀疏,導(dǎo)致局部特征丟失。

(a) (b)圖2 體素降采樣示例

1.3 圓柱提取

點(diǎn)云分割有聚類分割[18]、區(qū)域生長(zhǎng)[16]、RANSA方法分割[19]和語(yǔ)義分割[20]等,本節(jié)采用區(qū)域生長(zhǎng)算法對(duì)點(diǎn)云進(jìn)行分割從而提取圓柱。區(qū)域生成算法的基本思想是將具有相同特征的點(diǎn)加入到同一區(qū)域,本文的相同特征是指兩個(gè)點(diǎn)云的法向量的夾角小于設(shè)置的平滑閾值。設(shè)一組點(diǎn)云數(shù)據(jù)為P={P1,P2,…,Pn},區(qū)域生成具體步驟包括種子點(diǎn)的選取、生長(zhǎng)準(zhǔn)則、終止條件。種子點(diǎn)選取的是曲率值最小的點(diǎn),即最平滑的點(diǎn)。具體步驟如下:

(1) 對(duì)點(diǎn)Pi用kd_tree查找其鄰域。kd_tree是一種帶有約束條件的二叉查找樹(shù),每一級(jí)都在當(dāng)前維度上劃分開(kāi)所有的節(jié)點(diǎn),每個(gè)節(jié)點(diǎn)都表示一個(gè)超平面。kd_tree的示例如圖3所示。

圖3 kd_tree示例

(2) 設(shè)Ci為Pi的k個(gè)鄰域點(diǎn)構(gòu)成的協(xié)方差矩陣,則有:

假設(shè)求得式(3)的特征值滿足λi1≤λi2≤λi3,對(duì)應(yīng)的特征向量為vi1、vi2、vi3,則最小特征值λi1對(duì)應(yīng)的特征向量vi1即為Pi的法向量,鄰域曲率[21]為:

鄰域曲率反映了該點(diǎn)的鄰域所構(gòu)成的曲面的彎曲程度,σi越小表明曲面越平坦越接近平面。計(jì)算出每一點(diǎn)的鄰域曲率σi,對(duì)σi進(jìn)行排序,選擇σi最小的點(diǎn)作為種子點(diǎn)。

(3) 在進(jìn)行區(qū)域生長(zhǎng)的時(shí)候,計(jì)算鄰域點(diǎn)與種子點(diǎn)的法向量之間的夾角,若夾角小于設(shè)置的平滑閾值,則把該鄰域點(diǎn)加入到當(dāng)前區(qū)域。

(4) 當(dāng)所有滿足條件的點(diǎn)都被合并到當(dāng)前區(qū)域之后,在沒(méi)有被合并的點(diǎn)中選擇曲率最小的點(diǎn)作為種子點(diǎn),返回第(3)步。終止條件是所有點(diǎn)都被合并。

(5) 點(diǎn)云分割完成之后,直接把圓柱點(diǎn)云提取出來(lái)即可。

2 圓柱擬合

2.1 直線擬合

設(shè)圓柱中心軸線所在的直線為L(zhǎng),L的方向向量為(a,b,c),對(duì)L的方向向量進(jìn)行歸一化處理,得到a2+b2+c2=1。設(shè)L上某一點(diǎn)的坐標(biāo)為(x0,y0,z0),則L的標(biāo)準(zhǔn)式方程為:

把式(5)寫(xiě)成參數(shù)方程為:

式中:t為參數(shù)變量。則L可以表示為:

L=v0+Dt

(7)

式中:v0=(x0,y0,z0)T;D=(a,b,c)T。

圖4 點(diǎn)到直線的距離

建立目標(biāo)方程:

式中:n為測(cè)量點(diǎn)的個(gè)數(shù)。

首先,求初始值v0。對(duì)v0求導(dǎo)并令其導(dǎo)數(shù)為零:

故:

又有:

因此:

即證明了L過(guò)測(cè)量點(diǎn)的重心。

把式(9)展開(kāi)用點(diǎn)的坐標(biāo)表示:

(a(xi-x0)+b(yi-y0)+c(zi-z0))2]

(15)

對(duì)各測(cè)量點(diǎn)去重心:

式(15)變形為:

對(duì)各參數(shù)求導(dǎo)并令其為0,變換成矩陣形式為:

對(duì)矩陣進(jìn)行特征值分解,最大特征值對(duì)應(yīng)的特征向量即為L(zhǎng)的方向向量。

2.2 半徑擬合

旋轉(zhuǎn)點(diǎn)云使L與z軸平行,再將點(diǎn)云投影到xoy平面上,如圖5所示,將三維問(wèn)題轉(zhuǎn)換為二維問(wèn)題,再利用最小二乘原理擬合圓。

圖5 擬合半徑

首先將點(diǎn)云繞x軸旋轉(zhuǎn)θ角度,使L與yoz平面平行,旋轉(zhuǎn)矩陣為R1:

其中:

再將點(diǎn)云繞y軸旋轉(zhuǎn)β角度,使L與z軸平行,旋轉(zhuǎn)矩陣為R2:

其中:

設(shè)點(diǎn)云的半徑為r,圓心坐標(biāo)為(xa,ya),根據(jù)最小二乘原理寫(xiě)出誤差方程:

對(duì)各參數(shù)求導(dǎo)并等于0,令:

變換為矩陣形式為:

(26)

解方程求出圓心坐標(biāo)(xa,ya),代入求導(dǎo)公式中,求出r,再通過(guò)反坐標(biāo)變換求出在圓柱上的圓心點(diǎn)。

3 實(shí)例與分析

為驗(yàn)證本方法的有效性,采用三維掃描儀對(duì)桌面進(jìn)行掃描,對(duì)得到的數(shù)據(jù)進(jìn)行處理,可視化如圖6(a)所示。按照第2節(jié)提到的方法,首先對(duì)初始點(diǎn)云點(diǎn)云進(jìn)行統(tǒng)計(jì)濾波去噪,本文選擇的鄰域點(diǎn)個(gè)數(shù)是100,標(biāo)準(zhǔn)差倍數(shù)是0.8,去噪后的點(diǎn)云如圖6(b)所示,數(shù)目為130 251,可以看出,原始點(diǎn)云中的噪聲點(diǎn)基本已去除,靠近圓柱邊緣的點(diǎn)還存在小部分,在進(jìn)行區(qū)域生長(zhǎng)的時(shí)候會(huì)剔除掉。為了縮短后續(xù)區(qū)域分割的時(shí)間和提高圓柱擬合的效率,接下來(lái)對(duì)點(diǎn)云進(jìn)行降采樣,降低點(diǎn)云的密度。設(shè)置體素邊長(zhǎng)為0.03 m時(shí),點(diǎn)云數(shù)目為18 428,如圖7(a)所示,可以很明顯地看到點(diǎn)云密度降低了,并且保留了點(diǎn)云的形狀特征;設(shè)置體素邊長(zhǎng)為0.06 m時(shí),點(diǎn)云數(shù)目為5 237,如圖7(b)所示,點(diǎn)云更加稀疏,且形狀特征沒(méi)有發(fā)生變化。本文選用邊長(zhǎng)為0.06 m的體素進(jìn)行降采樣。

(a) 原始點(diǎn)云 (b) 去噪后的點(diǎn)云圖6 點(diǎn)云濾波去噪

(a) 體素邊長(zhǎng)0.03 m(b) 體素邊長(zhǎng)0.06 m圖7 不同邊長(zhǎng)的體素對(duì)點(diǎn)云進(jìn)行降采樣

在進(jìn)行區(qū)域生長(zhǎng)分割時(shí),根據(jù)實(shí)際情況選擇合適的曲率閾值和平滑閾值,本文中設(shè)置的平滑閾值為0.05π,曲率閾值為0.01,將點(diǎn)云分割成平面、圓柱和邊緣噪聲點(diǎn),再通過(guò)設(shè)置最小簇類數(shù)目將噪聲點(diǎn)剔除,分割完成之后,把圓柱點(diǎn)云提取出來(lái),如圖8所示。

圖8 提取出來(lái)的圓柱點(diǎn)云

提取出圓柱點(diǎn)云之后,按照本文方法對(duì)圓柱進(jìn)行擬合,求出的軸線方向向量為(0.025 31,-0.999 67,-0.005),擬合半徑r=0.762 87 m,圓柱的軸線方程可以表示為:

式中:t為參數(shù)變量。

平面圓的圓度是指點(diǎn)到圓心的距離與半徑之差,同理擴(kuò)展到三維空間,圓柱的圓度[4]是指圓柱面上的點(diǎn)到中心軸線的距離與半徑之差,計(jì)算結(jié)果如表1所示,圖9為圓度折線。對(duì)圓度進(jìn)行統(tǒng)計(jì),圓度均值為-0.01 mm,標(biāo)準(zhǔn)差為0.45 mm。

表1 部分觀測(cè)點(diǎn)的圓度值

圖9 圓度折線

為了進(jìn)一步驗(yàn)證本文方法的正確性和有效性,采用文獻(xiàn)[4]中的算法和數(shù)據(jù)進(jìn)行圓柱擬合,本文算法中的(x0,y0,z0)與文獻(xiàn)[4]中的(x0,y0,z0)都表示軸線上一點(diǎn),但是意義不同,所以沒(méi)有可比較性,只需要比較另外四個(gè)參數(shù)值即可。表2是本文方法與文獻(xiàn)[6]和文獻(xiàn)[10]中方法的比較,圖10是三種方法的圓度比較,可以看出,三種方法的差異是比較小的,且本文方法的圓度標(biāo)準(zhǔn)差略小,說(shuō)明了本文方法的正確性。

表2 與文獻(xiàn)[6]和文獻(xiàn)[10]結(jié)果對(duì)比

圖10 圓度對(duì)比

4 結(jié) 語(yǔ)

本文介紹了從散亂點(diǎn)云中提取圓柱點(diǎn)云并進(jìn)行圓柱擬合的一種方法。本文方法首先用統(tǒng)計(jì)濾波器和體素濾波器對(duì)點(diǎn)云進(jìn)行預(yù)處理,然后采用區(qū)域生長(zhǎng)算法對(duì)點(diǎn)云進(jìn)行分割,得到擬合算法需要的數(shù)據(jù),然后利用點(diǎn)到直線的距離建立方程,求出軸線,然后通過(guò)旋轉(zhuǎn)和投影,將三維問(wèn)題轉(zhuǎn)換為二維平面問(wèn)題,最后利用最小二乘法擬合出半徑。本文方法原理簡(jiǎn)單,便于實(shí)現(xiàn),而且具有很好的精度,最大的特點(diǎn)是避免了迭代。通過(guò)實(shí)驗(yàn)分析和比較,證明了本文方法的有效性和正確性。

在用濾波統(tǒng)計(jì)器進(jìn)行去噪的時(shí)候,鄰域點(diǎn)的個(gè)數(shù)和標(biāo)準(zhǔn)差倍數(shù)的選取對(duì)結(jié)果是有直接影響的,這里值得進(jìn)一步研究。