基于網(wǎng)格化曲率聚類的點云分割算法

基于網(wǎng)格化曲率聚類的點云分割算法

現(xiàn)在，在3d工程中，點云數(shù)據(jù)通常被劃分為具有獨特幾何特征的拓?fù)浣Y(jié)構(gòu)區(qū)域并進(jìn)行區(qū)域重建。因此，點云數(shù)據(jù)的分割是整個剖面重建過程的重要因素之一。多年研究形成了三種方法。首先，在基于邊緣的方法、邊緣的方法和混合方法（hybdi）的方法中檢測到分散邊界點，并將邊界點連接到邊界線作為整個塊過程的基礎(chǔ)。主曲線的方法通常也被廣泛使用。然而，由于噪聲和云的不規(guī)則性，它們對噪聲和點云的不規(guī)則性非常敏感。在基于面的方法中，首先選擇合適的種子數(shù)據(jù)，然后評估云旁邊點云的局部幾何特性（如法針、曲線等），并將具有類似幾何屬性的點云添加到當(dāng)前區(qū)域，以將數(shù)據(jù)點劃分為不同的區(qū)域。由于區(qū)域生長受到閾值的影響，因此很難選擇合適的生長模式。為了克服兩種方法的限制，yokyaa和leven提出了一種混合區(qū)域分割方法，該方法結(jié)合了基于邊緣和區(qū)域生長的方法?？梢钥闯?以上的區(qū)域分割方法的本質(zhì)在于找到點鄰域的微分信息,然后利用所求出的微分信息來識別邊界點或邊界曲線.在求取曲率等微分幾何信息之前,一般首先要進(jìn)行三角化,然后進(jìn)行微分信息的求解,計算量比較大.為了克服計算量大的問題,柯映林等提出對點云數(shù)據(jù)進(jìn)行網(wǎng)格劃分建立散亂點之間的拓?fù)潢P(guān)系,并將距空間柵格中心最近點的曲率作為柵格曲率,根據(jù)空間柵格與其26近鄰的曲率差識別邊特征柵格,最后基于連通性準(zhǔn)則實現(xiàn)空間散亂數(shù)據(jù)的區(qū)域分割.本文在此基礎(chǔ)上采用一種優(yōu)化的動態(tài)網(wǎng)格劃分方法,下面具體討論這種算法.1點云聚類分析本算法的基本思想是:首先對點云數(shù)據(jù)進(jìn)行動態(tài)網(wǎng)格劃分,建立點云的拓?fù)潢P(guān)系,求出網(wǎng)格內(nèi)距離網(wǎng)格中心點最近的點Pi,擬合局部拋物面求出該點Pi的曲率等微分信息,然后用K-means方法利用特征相似性進(jìn)行點云數(shù)據(jù)的聚類,分割特征區(qū)域.1.1ymin-e,zmin-e針對散亂點云分布的不規(guī)則性,以及直接網(wǎng)格劃分時Pi點的微分特性無法代表所有點的微分特性,本文采用動態(tài)網(wǎng)格劃分的方法,對點云數(shù)據(jù)進(jìn)行網(wǎng)格化建立點云之間的拓?fù)潢P(guān)系.采用動態(tài)層法和局部重劃法相結(jié)合的動態(tài)網(wǎng)格劃分方法對點云特征突變區(qū)等部分網(wǎng)格進(jìn)行動態(tài)調(diào)整.設(shè)點云數(shù)據(jù)中最大與最小x、y、z坐標(biāo)值分別是{Xmax,Ymax,Zmax,Xmin,Ymin,Zmin},且給定容差e,經(jīng)多次實驗,容差設(shè)定為點間距的3倍.并隨機抽取3對鄰近點云,估算點云數(shù)據(jù)的密度.動態(tài)網(wǎng)格劃分的具體步驟如下:1)首先定義以點(Xmax+e,Ymax+e,Zmax+e)和(Xmin-e,Ymin-e,Zmin-e)為對角點,且表面平行于坐標(biāo)平面的空間六面體為點云的空間包圍盒.2)根據(jù)點云密度設(shè)定閾值t,將點云包圍盒沿坐標(biāo)軸方向按等間隔t劃分成空間六面體網(wǎng)格,則在包圍盒內(nèi),沿X,Y,Z方向,以t為步長分段.可以得到,在三個方向上劃分的段數(shù)分別為m,n,l:{m=ΙΝΤ((Xmax+e)-(Xmin-e)t)+1n=ΙΝΤ(Ymax+e)-(Ymin-e)t)+1l=ΙΝΤ((Ζmax+e)-(Ζmin-e)t)+1(1)那么對于空間內(nèi)任意一個點P,它都會屬于某一個位置的基本正方體(i,j,k)內(nèi),它的坐標(biāo)(xp,yp,zp)滿足:{xp∈((Xmin-e+i×t)?(Xmin-e+(i+1)×t))yp∈((Ymin-e+j×t)?(Ymin-e+(j+1)×t))zp∈((Ζmin-e+k×t)?(Ζmin-e+(k+1)×t))(2)3)統(tǒng)計每一個網(wǎng)格內(nèi)點云的個數(shù)并計算點云密度,將空網(wǎng)格刪掉.求出每一個網(wǎng)格的中心點Oi,在網(wǎng)格數(shù)據(jù)集中找到離中心點Oi最近的點Pi.對網(wǎng)格內(nèi)數(shù)據(jù)在XYZ方向進(jìn)行排序,判斷數(shù)據(jù)的分布狀況.a.如果網(wǎng)格內(nèi)的數(shù)據(jù)分布比較密,點云個數(shù)超過了拋物面擬合的需要,但是點云數(shù)據(jù)之間在X,Y,Z方向的角度偏差在一定的允許范圍內(nèi),則不必對網(wǎng)格進(jìn)行調(diào)整.b.如果網(wǎng)格內(nèi)的數(shù)據(jù)分布比較密,點云個數(shù)超過了拋物面擬合的需要,但是點云數(shù)據(jù)之間在X,Y,Z方向的角度偏差超出了允許范圍,則根據(jù)點云偏差的突變,將原網(wǎng)格在一定方向上進(jìn)行細(xì)分,添加網(wǎng)格.c.如果網(wǎng)格內(nèi)數(shù)據(jù)分布均勻,點云個數(shù)滿足拋物面擬合的需要,并且點云數(shù)據(jù)之間在X,Y,Z方向的角度偏差在一定的允許范圍內(nèi),則不必對網(wǎng)格進(jìn)行調(diào)整.d.如果網(wǎng)格內(nèi)的數(shù)據(jù)分布比較稀疏,點云個數(shù)不滿足拋物面擬合的需要,但是點云數(shù)據(jù)之間在X,Y,Z方向的角度偏差在一定的允許范圍內(nèi),則判斷鄰近網(wǎng)格內(nèi)數(shù)據(jù)的分布情況,以原網(wǎng)格為中心,邊長2t擴(kuò)大網(wǎng)格,或沿某一方向?qū)⒕W(wǎng)格與鄰近網(wǎng)格合并.然后重新判斷,直到滿足需要.1.2基于變異值分解的解析法法矢和曲率是曲面的基本特性,也是曲面特征識別的重要依據(jù)之一.空間數(shù)據(jù)點曲面曲率計算方法一般有五種,分別是坐標(biāo)轉(zhuǎn)化法、正交步進(jìn)法、交叉曲面片法、曲面三角化法和TW法.1992年stokely和wu對這五種曲率計算方法進(jìn)行了試驗比較,證明坐標(biāo)轉(zhuǎn)化法是這幾種方法中比較實用的方法,并且計算步驟亦比較簡單.本文亦采用坐標(biāo)轉(zhuǎn)化法,在每個網(wǎng)格內(nèi)進(jìn)行局部曲面擬合.其基本思想是:每一個網(wǎng)格內(nèi),求出距離網(wǎng)格中心點最近的點集中的點Pi,利用Pi的k近鄰擬合拋物面S,通過計算S的曲率特性來確定該數(shù)據(jù)點的曲率特性.每一個網(wǎng)格內(nèi)點Pi的k近鄰,記為Nd(Pi)擬合平面估算測量點法矢,目標(biāo)函數(shù)為min∑p∈Νd(pi)((p-pi)ni)2(3)經(jīng)簡單計算,可得到對稱的半正定矩陣(3階矩陣):Cv=∑p∈Νd(pi)(p-pi)(p-pi)T(4)式中:Cv最小特征值對應(yīng)的特征向量單位化,即為曲面在點Pi的法矢Ni.多次實驗表明,當(dāng)k取為24~30,曲率計算一般都能獲得較好效果.在點Pi處建立局部坐標(biāo)系(u,v,h),且坐標(biāo)原點為Pi,坐標(biāo)軸h設(shè)為法矢Ni,坐標(biāo)軸u和v可以在過點Pi垂直于Ni的平面內(nèi)任意確定.將點Pi的k近鄰由全局坐標(biāo)(xi,yi,zi)轉(zhuǎn)化為局部坐標(biāo)系(ui,vi,hi),并用轉(zhuǎn)換后數(shù)據(jù)擬合如下面的拋物面.S(u,v)=a0+a1u+a2v+a3uv+a4u2+a5v2(5)拋物面擬合通過求解目標(biāo)函數(shù)mink∑i=1[hi-(a0+a1u+a2v+a3uivi+a4u2i+a5v2i)]2實現(xiàn).應(yīng)用奇異值分解法可獲得擬合拋物面的最小二乘解.具體實現(xiàn)過程如下:Step1將網(wǎng)格鄰域內(nèi)點的坐標(biāo)代入式拋物面S方程中,設(shè)每個網(wǎng)格內(nèi)有k個點,則得到方程組AX=Z(6)式中有A=|1,u0,v0,u0,v0,u20,v201,u1,v1,u1v1,u21,v211,uk,vk,ukvk,u2k,v2k|X={a0,a1,a2,a3,a4,a5}Τ,Ζ={h0,h1,?,hk}ΤStep2采用奇異值分解法求解線性方程組(5),由X=A′Z求解曲面方程的系a,b,c.Step3通過對擬合拋物面S(u,v,h)=S(u,v,au2+buv+cv2)求一、二階偏導(dǎo)數(shù),得到Pi點(局部坐標(biāo)系的原點)的平均曲率與高斯曲率:Su=(1,0,2au+bv)Suu=(0,0,2a)Sv=(0,1,bu+2cv)Svv=(0,0,2c)Suv=(0,0,b)Pi點的法矢為:n=Su×Sv|Su×Sv|由拋物面的第一、第二基本形式,可以得到:E=SuSu,F=SuSv,G=SvSv,L=Suu·n,M=Suv·n,N=Svv·n則得到Pi點的平均曲率(H)和高斯曲率(K):Κ=LΝ-Μ2EG-F2Η=EΝ-2FΜ+GL2(EG-F2)(7)1.3網(wǎng)格點云產(chǎn)業(yè)分割的實現(xiàn)通過計算得到每一個空間網(wǎng)格的平均曲率、高斯曲率和法矢等微分幾何信息后,通過曲率、法矢與曲面的關(guān)系,采用K-means的方法,對點云數(shù)據(jù)進(jìn)行區(qū)域分割.首先為了更快更準(zhǔn)的進(jìn)行分割,根據(jù)如下圖1,2所示的曲率分布圖,選擇比較合適的K個(K=2)初始聚類中心(初始聚類中心也可以隨機選取的),并把每個對象分配給離它最近的中心,從而得到兩個初始聚類.然后,計算出當(dāng)前每個聚類的中心作為新的聚類中心,并把每個對象重新分配到最近的中心.如果新的聚類的質(zhì)量優(yōu)于原先的聚類,則用新聚類代替原聚類.循環(huán)執(zhí)行這一過程直至聚類質(zhì)量不再提高為止或設(shè)定一定的迭代次數(shù).用K1,K2表示初始分割后的兩個區(qū)域,經(jīng)過了初始分割后,每一個空間網(wǎng)格的點,都帶有一個屬性值Kind,其值等于K1,K2中的一個.利用網(wǎng)格的拓?fù)潢P(guān)系,比較每一個網(wǎng)格和其鄰近的26個網(wǎng)格之間的Kind屬性以及曲率的變化趨勢,來檢查和糾正初始分割結(jié)果,并對點云數(shù)據(jù)進(jìn)行重新分割.具體實現(xiàn)過程如下:按照網(wǎng)格的索引號(i,j,k),搜索每一個網(wǎng)格和其鄰近網(wǎng)格的分類情況和曲率的變化趨勢.可以分為下面三種情況:1)如果該網(wǎng)格內(nèi)的Kind值與其他鄰近網(wǎng)格的Kind值全不相同,則該點為分類錯誤的點,重新對其進(jìn)行分類.2)如果該網(wǎng)格內(nèi)的Kind值與其他鄰近網(wǎng)格的Kind值不完全相同,則利用曲率的變化趨勢,判斷其分類是否正確.3)如果該網(wǎng)格內(nèi)的Kind值與其他鄰近網(wǎng)格的Kind值全相同,并且曲率無明顯變化趨勢,則證明該點是分類正確的點.依照上面的方法,完成所有網(wǎng)格點云數(shù)據(jù)的檢查和糾正,在利用網(wǎng)格的拓?fù)潢P(guān)系,將曲率值沒有明顯的區(qū)別,但是,物理上屬于不用曲面的點云數(shù)據(jù)進(jìn)行分割,得到比較滿足精度的分割后的數(shù)據(jù).2點云數(shù)據(jù)區(qū)域分割作者提出的算法已在C#平臺上,以directX為圖形顯示工具編程實現(xiàn).本實例使用由萊卡公司的近距離掃描儀HDS4500掃描獲得的點云數(shù)據(jù),如下圖3(a)所示,所掃描的是古建筑天花上梁的一部分,包含13917個散亂點.在進(jìn)行空間網(wǎng)格劃分時,考慮到點與點之間的距離在6?mm左右,所以網(wǎng)格尺寸是50?mm.古建筑的構(gòu)件結(jié)構(gòu)一般是由平面和圓柱或是圓弧等幾何特征組成,從下圖3(a)可以看出該構(gòu)件是由平面和圓柱面組成.采用本文方法對點云數(shù)據(jù)進(jìn)行區(qū)域分割,由曲率和二次曲面的關(guān)系,可以得知經(jīng)過該算法初始分割后,只能將平面和圓柱面初始分割開,如下圖3(b)所示是點云數(shù)據(jù),是經(jīng)過曲率的分割后的數(shù)據(jù)結(jié)果顯示,黑色為平面幾何特征,黃色為圓柱面幾何特征.對初始分割后的點云數(shù)據(jù),通過空間網(wǎng)格拓?fù)潢P(guān)系的檢查、糾正和二次分割后,得到如下圖3(c)所顯示的比較滿意的點云數(shù)據(jù)區(qū)域分割效果.用不同的顏色表示不同的區(qū)域.3采用三角化運算和動態(tài)網(wǎng)格本文通過點云數(shù)據(jù)動態(tài)網(wǎng)格的空間劃分建立測