1. 程式人生 > >點雲PCL估計一個點雲的表面法線

點雲PCL估計一個點雲的表面法線

表面法線是幾何體表面的重要屬性,在很多領域都有大量應用,例如:在進行光照渲染時產生符合可視習慣的效果時需要表面法線資訊才能正常進行,對於一個已知的幾何體表面,根據垂直於點表面的向量,因此推斷表面某一點的法線方向通常比較簡單。然而,由於我們獲取的點雲資料集在真實物體的表面表現為一組定點樣本,這樣就會有兩種解決方法:

使用曲面重建技術,從獲取的點雲資料集中得到取樣點對應的曲面,然後從曲面模型中計算表面法線;

直接從點雲資料集中近似推斷表面法線。

本小節將針對後一種情況進行講解,已知一個點雲資料集,在其中的每個點處直接近似計算表面法線。

理論基礎

儘管有許多不同的法線估計方法,本教程中著重講解的是其中最簡單的一個,表述如下,確定表面一點法線的問題近似於估計表面的一個相切面法線的問題,因此轉換過來以後就變成一個最小二乘法平面擬合估計問題。

注意:更多資訊,包含最小二乘法問題的數學方程式,請見[RusuDissertation]

因此估計表面法線的解決方案就變成了分析一個協方差矩陣的特徵向量和特徵值(或者PCA—主成分分析),這個協方差矩陣從查詢點的近鄰元素中建立。更具體地說,對於每一個點Pi,對應的協方差矩陣C,如下:

(1)

此處,k是點Pi鄰近點的數目,表示最近鄰元素的三維質心。這裡應該是三維的,即包括x,y,z座標,建立座標系,假設第i個點的座標是ri,質心的座標就是Σmi*ri/(Σmi)[1],姑且認為他是正確的,是協方差矩陣的第j個特徵值,是第j個特徵向量。這裡,特徵向量的意義是經過過這種特定的變換後保持方向不變。只是進行長度上的伸縮而已。而特徵值的意義是一個變換(矩陣)可由它的所有特徵向量完全表示,而每一個向量所對應的特徵值,就代表了矩陣在這一向量上的貢獻率——說的通俗一點就是能量(power)[2]。

在PCL內估計一點集對應的協方差矩陣,可以使用以下函式呼叫實現:

//定義每個表面小塊的3x3協方差矩陣的儲存物件
Eigen::Matrix3fcovariance_matrix;
//定義一個表面小塊的質心座標16-位元組對齊儲存物件
Eigen::Vector4fxyz_centroid;
//估計質心座標
compute3DCentroid(cloud,xyz_centroid);
//計算3x3協方差矩陣
computeCovarianceMatrix(cloud,xyz_centroid,covariance_matrix);

matlab中求協方差矩陣

cov(X) 求矩陣X的協方差矩陣。
diag(cov(X))得到每一個列向量的方差。sqrt(diag(cov(X)))得到每一個列的標準差。這裡與(1)還有些不同,(1)式是由給定一點,協方差矩陣是由其周圍的k鄰近個點決定。

通常,沒有數學方法能解決法線的正負向問題,如上所示,通過主成分分析法(PCA)來計算它的方向也具有二義性,無法對整個點雲資料集的法線方向進行一致性定向。圖1顯示出對一個更大資料集的兩部分產生的影響,此資料集來自於廚房環境的一部分,很明顯估計的法線方向並非完全一致,圖2展現了其對應擴充套件的高斯影象(EGI),也稱為法線球體(normal sphere),(我的理解是所有點雲的方向在一個球體中的表示。)它描述了點雲中所有法線的方向。由於資料集是2.5維,其只從一個單一的視角獲得,因此法線應該僅呈現出一半球體的擴充套件高斯影象(EGI)。然而,由於定向的不一致性,它們遍佈整個球體,如圖2所示。

                  圖1 估計廚房環境的表面法線                           圖2 法線球體

如果實際知道視點Vp,那麼這個問題的解決是非常簡單的。對所有法線定向只需要使它們一致朝向視點方向,滿足下面的方程式:

兩個向量點積大於0

下面的圖3展現了上面圖1中的資料集的所有法線被一致定向到視點後的結果演示。

圖3 將所有法線一致定向到後視點的結果

在PCL中對一個已知點的法線進行手動重新定向,你可以使用:

flipNormalTowardsViewpoint (const PointT &point, float vp_x, float vp_y, float vp_z,
Eigen
::Vector4f &normal);

注意:如果資料集是從多個捕獲視點中配準後集成的,那麼上述法線的一致性定向方法就不適用了。需要使用更復雜的演算法,更多資訊,請見[RusuDissertation]

選擇合適的尺度

如之前介紹的,在估計一個點的表面法線時,我們需要從周圍支援這個點的鄰近點著手(也稱作k鄰域)。最近鄰估計問題的具體內容又提出了另一個問題“合適的尺度”:已知一個取樣點雲資料集,k的正確取值是多少(k通過pcl::Feature::setKSearch給出)或者確定一個點r為半徑的圓內的最近鄰元素集時使用的半徑r應該取什麼值(r通過pcl::Feature::setRadiusSearch給出)。這個問題非常重要,並且在一個點特徵運算元的自動估計時(例如:使用者沒有給定閾值)是一個限制因素。為了更好地說明這個問題,以下圖示表現了選擇更小尺度(如:r值或k取相對小)與選擇更大尺度(如:r值或k值比較大)時的兩種不同效果。圖4和圖5分別為近檢視和遠檢視,兩圖中左邊部分展示選擇了一個合理的比例因子,估計的表面法線近似垂直於兩個平面,即使在互相垂直的邊沿部分,可明顯看到邊沿。如果這個尺度取的太大(右邊部分,k選擇的過大或者r選擇的過大),這樣鄰近點集將更大範圍地覆蓋鄰近表面的點,估計的點特徵表現就會扭曲失真,在兩個平面邊緣處出現旋轉表面法線,以及模糊不清的邊界,這樣就隱藏了一些細節資訊。

                    圖4選擇合理的比例因子                            圖5選擇較大的比例因子

無法深入探究更多討論,現在可粗略假設,以應用程式所需的細節需求為參考,選擇確定點的鄰域所用的尺度。簡言之,如果杯子手柄和圓柱體部分之間邊緣的曲率是重要的,那麼需要足夠小的尺度來捕獲這些細節資訊,而在其他不需要細節資訊的應用中可選擇`的尺度。

估計法線例項詳解

雖然法線估計的例子已經在上小節中給過了,這裡我們還是複習一下其中的一部分,以便更好地解釋說明本節的後續部分,在PCLPoint
Cloud Learning
)中國協助發行的書[1]
本書提供光碟的第12章例1資料夾中,開啟名為normal_estimation.cpp的程式碼檔案,同文件夾下可以找到相關的測試點雲檔案table_scene_lms400.pcd。下面詳細介紹一些原始碼中相關的函式和類。

法線估計類NormalEstimation的實際計算呼叫程式內部執行以下操作:

對點雲P中的每個點p

  1.得到p點的最近鄰元素

  2.計算p點的表面法線n

  3.檢查n的方向是否一致指向視點,如果不是則翻轉

視點座標預設為(0,0,0),可以使用以下程式碼進行更換:

setViewPoint (float vpx, float vpy, float vpz);

計算單個點的法線,使用:

computePointNormal (const pcl::PointCloud<PointInT>&cloud, const std::vector<int>&indices,
Eigen
::Vector4f &plane_parameters, float&curvature);

此處,cloud是包含點的輸入點雲,indices是點的k-最近鄰元素集索引,plane_parameters和curvature是法線估計的輸出,plane_parameters前三個座標中,以(nx, ny, nz)來表示法線,第四個座標D = nc . p_plane (centroid here) + p(這個座標值D不是很理解,可以參考開原始碼的說明)。輸出表面曲率curvature通過協方差矩陣的特徵值之間的運算估計得到,如:

圖6法線估計前原始點雲

最後利用光碟提供的CMakeLists.txt檔案,在cmake中建立工程檔案,並生成相應的可執行檔案,生成執行檔案後,將測試點雲檔案拷貝到與可執行檔案相同的路徑下或者改變原始碼中點雲的路徑,就可以運行了。如6所示,為估計所用的點雲視覺化效果,法線估計後執行結果如圖7所示,圖中的小線段表示法線。

圖7 法線估計結果視覺化

使用OpenMP加速法線估計

對於對運算速度有要求的使用者,PCL點雲庫提供了一個表面法線的附加實現程式,它使用多核/多執行緒開發規範,利用OpenMP來提高計算速度。它的類命名為pcl::NormaleEstimationOMP,並且它的應用程式介面(API100%相容單執行緒pcl::NormalEstimation,這使它適合作為一個可選提速方法。在8核系統中,可以輕鬆提速6-8倍。

參考資料:

[1] http://zhidao.baidu.com/link?url=e2XzUad5EULxhNCW7_7zHh9UKmQobW809H8VXixboioojGDma-T0MYWBG9YTo4eoUGj6KHePCGkpIXJ3sSIDUP5okdzKY3nMOrlRDp0X9Cy

[2]
http://wenku.baidu.com/view/5587010a581b6bd97f19ea74.html

[3]
朱德海、郭浩、蘇偉.點雲庫PCL學習教程(ISBN 978-7-5124-0954-5)北京航空航天出版社2012-10

轉自:http://www.xuebuyuan.com/2228903.html