肖正濤，高 健，吳東慶，張攬宇

(1. 廣東工業(yè)大學(xué)省部共建精密電子制造技術(shù)與裝備國(guó)家重點(diǎn)實(shí)驗(yàn)室，廣州 510006；2. 廣東工貿(mào)職業(yè)技術(shù)學(xué)院機(jī)電工程學(xué)院， 廣州 510510；3. 仲愷農(nóng)業(yè)工程學(xué)院計(jì)算科學(xué)學(xué)院， 廣州 510220)

0 引言

近年來，隨著深度傳感技術(shù)的發(fā)展與普及，三維點(diǎn)云的獲取變得越來越容易。由于從三維點(diǎn)云中提取特征不易受比例、旋轉(zhuǎn)和光照的影響，而且比從RGB圖像中得到的信息更準(zhǔn)確，因而深度傳感技術(shù)得到了越來越多的關(guān)注，被廣泛應(yīng)用在自動(dòng)化加工生產(chǎn)線、無人駕駛汽車等領(lǐng)域。

由于實(shí)際中得到的三維點(diǎn)云數(shù)據(jù)量龐大，且往往包含有一些噪聲信息，在使用前通常要進(jìn)行一些預(yù)處理[1]，包括去噪、降采樣(downsampling[2])、點(diǎn)云數(shù)據(jù)結(jié)構(gòu)化等，以提升點(diǎn)云的最終處理效果和效率。其中，降采樣是點(diǎn)云預(yù)處理過程中一個(gè)非常關(guān)鍵的步驟，它不僅會(huì)影響點(diǎn)云的處理速度，還會(huì)影響最終的處理效果。對(duì)點(diǎn)云進(jìn)行降采樣處理簡(jiǎn)而言之就是對(duì)點(diǎn)云進(jìn)行精簡(jiǎn)[3]，或者說就是減少點(diǎn)云中點(diǎn)的數(shù)量。點(diǎn)云降采樣的方法有很多，如體素網(wǎng)格(voxel grid[4])法、系統(tǒng)采樣(systematic sampling[5])、隨機(jī)采樣(random sampling[6])、最遠(yuǎn)點(diǎn)采樣(farthest point sampling[7])等，這些不同的降采樣方法通常有各自適合的使用場(chǎng)景。

體素網(wǎng)格法，也稱體素柵格法，是一種常用的降采樣方法。該方法先對(duì)三維點(diǎn)云建立軸向包圍盒(axis-aligned bounding box, 簡(jiǎn)稱AABB[8])，然后沿各個(gè)坐標(biāo)軸方向?qū)鼑蟹殖蒼等份，接著計(jì)算每一個(gè)體素中所有點(diǎn)的重心，最后將該重心作為該體素的采樣值。體素網(wǎng)格法降采樣的優(yōu)點(diǎn)是計(jì)算效率高、采樣點(diǎn)分布均勻，因而該方法得到廣泛使用。

點(diǎn)對(duì)特征(point pair feature，簡(jiǎn)稱PPF)三維物體識(shí)別算法由Drost B[9]在2010提出，是目前綜合性能表現(xiàn)最佳的三維物體識(shí)別算法之一[10]，被廣泛應(yīng)用在機(jī)器人自主作業(yè)系統(tǒng)。Birdal T等[11]在OpenCV[12]中實(shí)現(xiàn)了PPF三維物體識(shí)別算法模塊。在Birdal的源代碼中，點(diǎn)云預(yù)處理階段使用了體素網(wǎng)格降采樣方法。筆者查看了Birdal在2014年9月提交的最初版本和在2020年1月提交的最新版本，Birdal在將體素的三維索引轉(zhuǎn)換為一維索引時(shí)一直存在著錯(cuò)誤。

本文針對(duì)OpenCV點(diǎn)對(duì)特征三維物體識(shí)別模塊中體素網(wǎng)格降采樣算法進(jìn)行了深入研究，詳細(xì)講解了體素網(wǎng)格降采樣方法的網(wǎng)格歸屬劃分和體素的索引轉(zhuǎn)換，分析了Birdal方法存在的問題，以及Birdal方法在降采樣時(shí)可能會(huì)導(dǎo)致的后果，并給出了正確的體素三維索引到一維索引的轉(zhuǎn)換表達(dá)式。通過實(shí)驗(yàn)比較，在極端情形，Birdal的方法不能正確降采樣，而本文方法能正確完成降采樣任務(wù)。

1 體素網(wǎng)格降采樣方法

1.1 體素網(wǎng)格方法的實(shí)現(xiàn)

如圖1a所示，對(duì)點(diǎn)云建立軸向包圍盒，并沿各個(gè)坐標(biāo)軸方向?qū)鼑蟹殖蒼等份。在包圍盒最上、最右和最后添加如圖1b所示的虛擬體素，以接納位于包圍盒最上、最右和最后表面或棱線上的點(diǎn)。圖1c將體素分組，不同的組用不同的顏色表示。圖1d對(duì)包括假想體素在內(nèi)的包圍盒沿x、y、z三個(gè)坐標(biāo)軸方向分別標(biāo)記上索引，索引從0開始計(jì)數(shù)。

點(diǎn)云中任意點(diǎn)所在體素的三維索引(i,j,k)的計(jì)算式如下所示：

(1)

式中，i,j,k分別表示x、y、z軸方向從0開始的索引值，trunc表示截?cái)嗳≌?也稱舍棄小數(shù)取整)，n表示包圍盒沿x、y、z軸分成的份數(shù)(從1開始計(jì)數(shù))，dx表示點(diǎn)到包圍盒左側(cè)面的距離，Lx表示包圍盒沿x軸方向的邊長(zhǎng)，dy表示點(diǎn)到包圍盒前表面的距離，Ly表示包圍盒沿y軸方向的邊長(zhǎng)，dz表示點(diǎn)到包圍盒下表面的距離，Lz表示包圍盒沿z軸方向的邊長(zhǎng)。從式(1)可知，i,j,k的范圍是0～n。

為了程序處理上的便利，我們需要將所有體素(包括虛擬體素)的三維索引轉(zhuǎn)換為一維索引，如圖1d所示。如果按照先x軸方向，然后y軸方向，最后z軸方向的順序，轉(zhuǎn)換式如下：

(c)將體素分組，不同的組 用不同的顏色表示 (d) 體素的三維索引到 一維索引的轉(zhuǎn)換

idx=i*(n+1)(n+1)+j*(n+1)+k

(2)

式中，idx表示體素的一維索引，i,j,k∈[0,n]。這樣，每一個(gè)三維索引(i,j,k)就唯一對(duì)應(yīng)著一個(gè)一維索引idx。每一個(gè)一維索引idx也同樣唯一對(duì)應(yīng)著一個(gè)三維索引(i,j,k)，從式(2)可推知一維索引到三維索引的計(jì)算式如下：

(3)

式(2)和式(3)建立起了體素三維索引與體素一維索引之間的一一對(duì)應(yīng)關(guān)系。

(a) 對(duì)點(diǎn)云建立軸向包圍盒， 并對(duì)包圍盒進(jìn)行等分 (b) 假想在包圍盒最上、最右 和最后添加虛擬的體素

式(1)對(duì)包圍盒內(nèi)和包圍表面(包括棱線和頂點(diǎn))上的點(diǎn)進(jìn)行了歸屬劃分，但不夠直觀，尤其是位于包圍盒Top，Right和Back三個(gè)表面(包括相應(yīng)的交線和頂點(diǎn))上的點(diǎn)。將假想添加的體素分解，如圖 2所示，每一組虛擬體素分別標(biāo)記為Ψ1,…,Ψ7。標(biāo)記出包圍盒的Front，Top和Right面，與Front，Top和Right面對(duì)應(yīng)的則是Back，Bottom和Left面。包圍盒Top，Right和Back三個(gè)表面(包括相應(yīng)的交線和交點(diǎn))上的點(diǎn)的歸屬分別是：頂點(diǎn)D屬于體素Ψ4，直線AD(除D點(diǎn))屬于體素Ψ5，直線BD(除D點(diǎn))屬于體素Ψ2，直線CD(除D點(diǎn))屬于體素Ψ6；Top面(除直線AD和直線BD)屬于體素Ψ1，Right面(除直線AD和直線CD)屬于體素Ψ7，Back面(除直線BD和直線CD)屬于體素Ψ3。這種歸屬劃分與式(1)是一致的。

圖2 分解圖效果

完成了體素網(wǎng)格的建立，空間點(diǎn)的劃分和體素三維索引到一維索引的轉(zhuǎn)換后，接著計(jì)算每一個(gè)體素中所有點(diǎn)的重心，最后將該重心作為該體素的采樣值。

1.2 對(duì)Birdal算法的分析

在Birdal實(shí)現(xiàn)的點(diǎn)對(duì)特征三維物體識(shí)別模塊的源代碼中，Birdal使用如下的表達(dá)式將體素的三維索引轉(zhuǎn)換為一維索引：

idx=i·n·n+j·n+k

(4)

其中，i,j,k的范圍是0～n。

筆者認(rèn)為式(4)是錯(cuò)誤的，式(4)會(huì)直接導(dǎo)致不同的體素對(duì)應(yīng)同一個(gè)一維索引。如索引為(1,0,0)的體素和索引為(0,n,0)的體素是不同的，但按照式(4)，它們對(duì)應(yīng)的一維索引皆為n2，后果就是將不同體素中的點(diǎn)當(dāng)作屬于同一個(gè)體素的點(diǎn)。

更進(jìn)一步探討，式(4)會(huì)將哪些不同的體素對(duì)應(yīng)為同一個(gè)體素呢?這個(gè)問題的數(shù)學(xué)描述就是：當(dāng)i1=i2,j1=j2,k1=k2不同時(shí)成立時(shí)，求方程i1·n·n+j1·n+k1=i2·n·n+j2·n+k2在i1,j1,k1,i2,j2,k2∈[0,n]且n≥2時(shí)的非負(fù)整數(shù)解。這是一個(gè)求解整系數(shù)線性方程組的整數(shù)解問題[13]，特別之處是此時(shí)方程組中的方程只有1個(gè)，未知數(shù)有6個(gè)。為了便于分析，我們將這個(gè)方程變換成如下形式：

n2(i1-i2)+n(j1-j2)+(k1-k2)=0

(5)

(6)

由于i1,j1,k1,i2,j2,k2∈[0,n]，因而i1-i2,j1-j2,k1-k2∈[-n,n]，根據(jù)式(6)，可得到下面這組線性不等式[14]：

(7)

不等式組(7)的解有6組，如下：

(8)

此時(shí)，式(6)即可轉(zhuǎn)化為下面6組方程組：

(9)

分析式(9)中的6組方程組，每組方程組等號(hào)右邊皆有一項(xiàng)不是n就是-n，也就是說，對(duì)于每組方程組，i1-i2,j1-j2,k1-k2中必有一項(xiàng)等于n或-n。由于i1,j1,k1,i2,j2,k2∈[0,n]，可知每組方程組中的i1,i2或j1,j2或k1,k2必有一組只會(huì)是0和n。

在三維索引(i,j,k)中，當(dāng)其中一個(gè)索引分量為0或n時(shí)，此三維索引對(duì)應(yīng)的體素一定位于包圍盒表面，如圖1d所示。我們已經(jīng)知道，Birdal用來轉(zhuǎn)換體素三維索引到一維索引的計(jì)算式(4)會(huì)導(dǎo)致不同的體素可能對(duì)應(yīng)同一個(gè)一維索引。現(xiàn)在我們更確切地知道，式(4)會(huì)將位于包圍盒表面的不同體素對(duì)應(yīng)為同一個(gè)一維索引。因而，Birdal算法對(duì)位于包圍盒表面上的點(diǎn)有可能不能正確降采樣。

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

本文的實(shí)驗(yàn)點(diǎn)云數(shù)據(jù)有兩個(gè)：①圖3a中人工制作的長(zhǎng)方體點(diǎn)云，該點(diǎn)云的所有點(diǎn)皆位于一個(gè)長(zhǎng)方體的6個(gè)表面，且長(zhǎng)方體的各個(gè)表面與相應(yīng)的坐標(biāo)面皆平行，這個(gè)特殊點(diǎn)云的軸向包圍盒為該長(zhǎng)方體自身；②圖3c中的加機(jī)工零件點(diǎn)云，該點(diǎn)云是用CREAFORM HandySCAN 700手持式三維激光掃描儀掃描圖3b所示的一堆機(jī)加工零件獲得，該點(diǎn)云只有非常有限的幾個(gè)點(diǎn)位于軸向包圍盒表面。圖3a和圖3c中的綠色線框?yàn)楦髯渣c(diǎn)云的軸向包圍盒。本文實(shí)驗(yàn)以O(shè)penCV和PCL (Point Cloud Library)為計(jì)算平臺(tái)，計(jì)算機(jī)配置是Intel 酷睿I5-7200U，12 GB內(nèi)存。本文分別使用Birdal算法和本文方法對(duì)三維點(diǎn)云進(jìn)行降采樣，并從降采樣的直觀效果、計(jì)算時(shí)間、降采樣前后變化等多個(gè)方面進(jìn)行了比較和分析。

長(zhǎng)方體點(diǎn)云和機(jī)加工零件點(diǎn)云包含的點(diǎn)數(shù)、點(diǎn)云軸向包圍盒的尺寸大小以及體素網(wǎng)格降采樣時(shí)包圍盒每邊將要分成的等分?jǐn)?shù)如表 1所示。

表1 實(shí)驗(yàn)點(diǎn)云的一些數(shù)據(jù)

用Birdal方法和本文方法分別對(duì)圖3a所示的長(zhǎng)方體點(diǎn)云進(jìn)行降采樣處理，包圍盒沿各個(gè)坐標(biāo)軸方向分成的份數(shù)皆為n=5。Birdal方法得到的降采樣結(jié)果如圖4a所示，除了x方向包圍盒左右兩個(gè)表面有點(diǎn)外，其余所有點(diǎn)皆落在包圍盒內(nèi)部，也就是說，包圍盒其余4個(gè)表面上皆沒有點(diǎn)。本文方法得到的降采樣結(jié)果如圖4b所示，包圍盒6個(gè)表面皆有點(diǎn)，只在靠左、靠前和靠下的位置有一些采樣點(diǎn)位于包圍盒靠近表面的內(nèi)部。直觀上看，兩種方法對(duì)長(zhǎng)方體點(diǎn)云進(jìn)行降采樣得到的結(jié)果差異非常大。

(a) 長(zhǎng)方體點(diǎn)云 (b) 機(jī)加工零件 (c) 機(jī)加工零件對(duì)應(yīng)的點(diǎn)云

用Birdal方法和本文方法分別對(duì)圖3c所示的機(jī)加工零件點(diǎn)云進(jìn)行降采樣處理，包圍盒沿各個(gè)坐標(biāo)軸方向分成的份數(shù)皆為n=40。Birdal方法得到的降采樣結(jié)果如圖4c所示，本文方法得到的降采樣結(jié)果如圖4d所示。直觀上看，兩種方法對(duì)機(jī)加工零件點(diǎn)云進(jìn)行降采樣得到的結(jié)果無明顯差異。

(a) Birdal方法，包圍盒半透明顯示 (b) 本文方法

表2展示了Birdal方法和本文方法對(duì)長(zhǎng)方體點(diǎn)云和機(jī)加工零件點(diǎn)云降采樣的一些數(shù)據(jù)。無論是對(duì)長(zhǎng)方體點(diǎn)云，還是對(duì)機(jī)加工零件點(diǎn)云，Birdal方法和本文方法在計(jì)算時(shí)間上沒有顯著差異。對(duì)長(zhǎng)方體點(diǎn)云來說，無論是降采樣后的點(diǎn)數(shù)量，還是降采樣后的點(diǎn)云到原始點(diǎn)云的距離，Birdal方法和本文方法皆有非常顯著的差異；從降采樣后的點(diǎn)云到原始點(diǎn)云的距離來看，Birdal方法得到的值皆比本文方法得到的值要大，說明Birdal方法得到的點(diǎn)云與原始點(diǎn)云的差異比本文方法要大。對(duì)機(jī)加工零件點(diǎn)云來說，無論是降采樣后的點(diǎn)數(shù)量，還是降采樣后的點(diǎn)云到原始點(diǎn)云的距離，Birdal方法和本文方法得到的結(jié)果是一樣的。

表2 降采樣結(jié)果

在計(jì)算降采樣后的點(diǎn)云到原始點(diǎn)云的距離時(shí)，我們可以得到降采樣后的點(diǎn)云中每一點(diǎn)與原始點(diǎn)云的最小距離，將這些最小距離值繪成頻數(shù)分布直方圖，如圖5所示。在圖5中，橫軸表示降采樣后的點(diǎn)云中的點(diǎn)到原始點(diǎn)云的最小距離，橫軸被分成若干個(gè)間隔，以方便進(jìn)行頻數(shù)統(tǒng)計(jì)?？v軸表示降采樣后的點(diǎn)云中到原始點(diǎn)云的距離在某個(gè)范圍內(nèi)的點(diǎn)有多少。為了便于直觀對(duì)比分析，圖5a和圖5b的橫軸和縱軸的刻度范圍、距離間隔皆一致，圖5c和圖5d也是如此。對(duì)長(zhǎng)方體點(diǎn)云而言，從圖5a和圖5b可知，Birdal方法得到的點(diǎn)云中有相當(dāng)一部分點(diǎn)與原始點(diǎn)云的距離偏大，而本文方法得到的點(diǎn)云總體上與原始點(diǎn)云的距離更小，也即本文方法得到點(diǎn)云更接近原始點(diǎn)云。圖5c和圖5d的頻數(shù)分布直方圖基本相同，表明兩種方法得到點(diǎn)云也基本相同。

(a) Birdal’s (b) Ours

表3展示了不同方法得到的降采樣點(diǎn)云之間的距離值，用于比較兩個(gè)點(diǎn)云之間的差異情況。通常情況下，若A，B兩個(gè)點(diǎn)云不同，那么從A點(diǎn)云到B點(diǎn)云的距離，與從B點(diǎn)云到A點(diǎn)云的距離往往是不同的。若A到B的距離越小，則表明A與B越相似。若A到B的距離為0，則表明A是B的子集，若同時(shí)B到A的距離也為0，則表明B也是A的子集，此時(shí)可知A與B相同。表3中，對(duì)降采樣后的長(zhǎng)方體點(diǎn)云而言，除了最小值這一項(xiàng)為0外，其余皆不為0。最小值這一項(xiàng)為0，表明Birdal方法得到的點(diǎn)云與本文方法得到的點(diǎn)云有重合的部分；其余項(xiàng)不為0，表明Birdal方法得到的點(diǎn)云與本文方法得到的點(diǎn)云在整體上是不相同的。表3中，對(duì)降采樣后的機(jī)加工零件點(diǎn)云而言，Birdal方法所得的點(diǎn)云到本文方法所得的點(diǎn)云的距離值皆為0，且本文方法所得的點(diǎn)云到Birdal方法所得的點(diǎn)云的距離值也皆為0，這表明兩種方法得到的降采樣后的機(jī)加工零件點(diǎn)云是一樣的。

表3 降采樣后的點(diǎn)云比較

從上面的結(jié)果和分析我們可知：Birdal方法對(duì)特殊的長(zhǎng)方體點(diǎn)云不能正確降采樣，得到的結(jié)果與原始點(diǎn)云有較大的差異，Birdal方法對(duì)加工零件點(diǎn)云可以正確降采樣；本文方法無論是對(duì)特殊的長(zhǎng)方體點(diǎn)云還是機(jī)加工零件點(diǎn)云皆能正確降采樣。

對(duì)Birdal方法而言，若點(diǎn)云包圍盒表面有點(diǎn)，且這些點(diǎn)所在的體素索引滿足方程(5)，則點(diǎn)云將不能被正確降采樣。當(dāng)這兩個(gè)條件不能同時(shí)滿足，使用Birdal方法就可以對(duì)點(diǎn)云進(jìn)行正確降采樣。這也就是為何機(jī)加工零件點(diǎn)云存在著幾個(gè)位于包圍盒表面上的點(diǎn)，而Birdal方法仍可以對(duì)該點(diǎn)云正確降采樣的原因所在。

3 結(jié)論

針對(duì)OpenCV中體素網(wǎng)格降采樣算法存在的錯(cuò)誤，本文提出了一種體素網(wǎng)格降采樣算法，該算法的關(guān)鍵是對(duì)點(diǎn)云中的每一個(gè)點(diǎn)正確劃分網(wǎng)格歸屬，并將網(wǎng)格的三維索引轉(zhuǎn)換為與之一一對(duì)應(yīng)的一維索引。理論分析和實(shí)驗(yàn)結(jié)果均表明，當(dāng)點(diǎn)云中的點(diǎn)位于包圍盒表面時(shí)，OpenCV中的體素網(wǎng)格降采樣算法不能對(duì)點(diǎn)云進(jìn)行正確降采樣，而本文提出的體素網(wǎng)格降采樣方法則可以正確有效地對(duì)各種情形的點(diǎn)云進(jìn)行降采樣。實(shí)驗(yàn)結(jié)果同時(shí)表明，兩種降采樣方法在計(jì)算時(shí)間上無明顯差異。