----------------法線求解原理--------------

表面法線是幾何體表面的重要屬性，在很多領(lǐng)域都有大量應(yīng)用，例如：在進(jìn)行光照渲染時(shí)產(chǎn)生符合可視習(xí)慣的效果時(shí)需要表面法線信息才能正常進(jìn)行，對于一個(gè)已知的幾何體表面，根據(jù)垂直于點(diǎn)表面的矢量，因此推斷表面某一點(diǎn)的法線方向通常比較簡單。然而，由于我們獲取的點(diǎn)云數(shù)據(jù)集在真實(shí)物體的表面表現(xiàn)為一組定點(diǎn)樣本，這樣就會(huì)有兩種解決方法：

（1）使用曲面重建技術(shù)，從獲取的點(diǎn)云數(shù)據(jù)集中得到采樣點(diǎn)對應(yīng)的曲面，然后從曲面模型中計(jì)算表面法線；曲面方程f(x,y,z)=0的一個(gè)法向量可以表示為n={e(df/dx), e(df/dy), e(df/dz)}.）

（2）直接從點(diǎn)云數(shù)據(jù)集中近似推斷表面法線。

下面以已知一個(gè)點(diǎn)云數(shù)據(jù)集，在其中的每個(gè)點(diǎn)處直接近似計(jì)算表面法線。為例進(jìn)行展開：

確定表面一點(diǎn)法線的問題近似于估計(jì)表面的一個(gè)相切面法線的問題，因此轉(zhuǎn)換過來以后就變成一個(gè)最小二乘法平面擬合估計(jì)問題。

最小二乘法（又稱最小平方法）是一種數(shù)學(xué)優(yōu)化技術(shù)。它通過最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。利用最小二乘法可以簡便地求得未知的數(shù)據(jù)，并使得這些求得的數(shù)據(jù)與實(shí)際數(shù)據(jù)之間誤差的平方和為最小。最小二乘法還可用于曲線擬合。其他一些優(yōu)化問題也可通過最小化能量或最大化熵用最小二乘法來表達(dá)。

對空間中的一系列散點(diǎn)，尋求一個(gè)近似平面：（按最小原則選擇的擬合平面為最小二乘擬合平面）

上述使用最小二乘我們已經(jīng)可以根據(jù)某點(diǎn)的臨近點(diǎn)閾，可以擬合出一個(gè)平面。接下來就是求平面的法線。表面法線是幾何體面的重要屬性。而點(diǎn)云數(shù)據(jù)集在真實(shí)物體的表面表現(xiàn)為一組定點(diǎn)樣本。對點(diǎn)云數(shù)據(jù)集的每個(gè)點(diǎn)的法線估計(jì)，可以看作是對表面法線的近似推斷。在PCL庫中有專門針對法向量計(jì)算的庫，但是必須了解計(jì)算原理才能記得更加深刻。

確定表面一點(diǎn)法線的問題近似于估計(jì)表面的一個(gè)相切面法線的問題，因此轉(zhuǎn)換過來以后就變成一個(gè)最小二乘法平面擬合估計(jì)問題。

下一張圖說明怎么是由求法線轉(zhuǎn)換到求最小特征值對應(yīng)的特征向量就是代表法線。

類似于上圖中的平面方程求解，只是換了一種平面表示方法。（由一般式轉(zhuǎn)換為法線式方程）

??

由上述推導(dǎo)過程可知，求法線就是求最小特征值對應(yīng)的特征向量。

----

還可以觀察上述右邊圖片中的A矩陣其實(shí)是個(gè)協(xié)方差矩陣；

在PCL內(nèi)估計(jì)一點(diǎn)集對應(yīng)的協(xié)方差矩陣，可以使用以下函數(shù)調(diào)用實(shí)現(xiàn)：

//定義每個(gè)表面小塊的3x3協(xié)方差矩陣的存儲(chǔ)對象
Eigen::Matrix3fcovariance_matrix;
//定義一個(gè)表面小塊的質(zhì)心坐標(biāo)16-字節(jié)對齊存儲(chǔ)對象
Eigen::Vector4fxyz_centroid;
//估計(jì)質(zhì)心坐標(biāo)
compute3DCentroid(cloud,xyz_centroid);
//計(jì)算3x3協(xié)方差矩陣
computeCovarianceMatrix(cloud,xyz_centroid,covariance_matrix);

通常，沒有數(shù)學(xué)方法能解決法線的正負(fù)向問題，如上所示，通過主成分分析法（PCA）來計(jì)算它的方向也具有二義性，無法對整個(gè)點(diǎn)云數(shù)據(jù)集的法線方向進(jìn)行一致性定向。圖1顯示出對一個(gè)更大數(shù)據(jù)集的兩部分產(chǎn)生的影響，此數(shù)據(jù)集來自于廚房環(huán)境的一部分，很明顯估計(jì)的法線方向并非完全一致，圖2展現(xiàn)了其對應(yīng)擴(kuò)展的高斯圖像（EGI），也稱為法線球體（normal sphere），它描述了點(diǎn)云中所有法線的方向。由于數(shù)據(jù)集是2.5維，其只從一個(gè)單一的視角獲得，因此法線應(yīng)該僅呈現(xiàn)出一半球體的擴(kuò)展高斯圖像（EGI）。然而，由于定向的不一致性，它們遍布整個(gè)球體，如圖2所示。

如果實(shí)際知道視點(diǎn)，那么這個(gè)問題的解決是非常簡單的。對所有法線定向只需要使它們一致朝向視點(diǎn)方向，滿足下面的方程式：

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?

下面的圖3展現(xiàn)了上面圖1中的數(shù)據(jù)集的所有法線被一致定向到視點(diǎn)后的結(jié)果演示。

在PCL中對一個(gè)已知點(diǎn)的法線進(jìn)行手動(dòng)重新定向，你可以使用：

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

注意：如果數(shù)據(jù)集是從多個(gè)捕獲視點(diǎn)中配準(zhǔn)后集成的，那么上述法線的一致性定向方法就不適用了。需要使用更復(fù)雜的算法。更詳細(xì)的信息見：http://pointclouds.org/documentation/tutorials/how_features_work.php#id1

選擇合適的尺度

如之前介紹的，在估計(jì)一個(gè)點(diǎn)的表面法線時(shí)，我們需要從周圍支持這個(gè)點(diǎn)的鄰近點(diǎn)著手（也稱作k鄰域）。最近鄰估計(jì)問題的具體內(nèi)容又提出了另一個(gè)問題“合適的尺度”：已知一個(gè)取樣點(diǎn)云數(shù)據(jù)集，k的正確取值是多少（k通過pcl::Feature::setKSearch給出）或者確定一個(gè)點(diǎn)r為半徑的圓內(nèi)的最近鄰元素集時(shí)使用的半徑r應(yīng)該取什么值（r通過pcl::Feature::setRadiusSearch給出）。這個(gè)問題非常重要，并且在一個(gè)點(diǎn)特征算子的自動(dòng)估計(jì)時(shí)（例如：用戶沒有給定閾值）是一個(gè)限制因素。為了更好地說明這個(gè)問題，以下圖示表現(xiàn)了選擇更小尺度（如：r值或k取相對小）與選擇更大尺度（如：r值或k值比較大）時(shí)的兩種不同效果。圖4和圖5分別為近視圖和遠(yuǎn)視圖，兩圖中左邊部分展示選擇了一個(gè)合理的比例因子，估計(jì)的表面法線近似垂直于兩個(gè)平面，即使在互相垂直的邊沿部分，可明顯看到邊沿。如果這個(gè)尺度取的太大（右邊部分），這樣鄰近點(diǎn)集將更大范圍地覆蓋鄰近表面的點(diǎn)，估計(jì)的點(diǎn)特征表現(xiàn)就會(huì)扭曲失真，在兩個(gè)平面邊緣處出現(xiàn)旋轉(zhuǎn)表面法線，以及模糊不清的邊界，這樣就隱藏了一些細(xì)節(jié)信息。

無法深入探究更多討論，現(xiàn)在可粗略假設(shè)，以應(yīng)用程序所需的細(xì)節(jié)需求為參考，選擇確定點(diǎn)的鄰域所用的尺度。簡言之，如果杯子手柄和圓柱體部分之間邊緣的曲率是重要的，那么需要足夠小的尺度來捕獲這些細(xì)節(jié)信息，而在其他不需要細(xì)節(jié)信息的應(yīng)用中可選擇大的尺度。

估計(jì)法線實(shí)例詳解：http://pointclouds.org/documentation/tutorials/index.php#features-tutorial

法線估計(jì)類NormalEstimation的實(shí)際計(jì)算調(diào)用程序內(nèi)部執(zhí)行以下操作：

對點(diǎn)云P中的每個(gè)點(diǎn)p1.得到p點(diǎn)的最近鄰元素2.計(jì)算p點(diǎn)的表面法線n3.檢查n的方向是否一致指向視點(diǎn)，如果不是則翻轉(zhuǎn)

視點(diǎn)坐標(biāo)默認(rèn)為（0,0,0），可以使用以下代碼進(jìn)行更換：

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

計(jì)算單個(gè)點(diǎn)的法線，使用：

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

此處，cloud是包含點(diǎn)的輸入點(diǎn)云，indices是點(diǎn)的k-最近鄰元素集索引，plane_parameters和curvature是法線估計(jì)的輸出，plane_parameters前三個(gè)坐標(biāo)中，以(nx, ny, nz)來表示法線，第四個(gè)坐標(biāo)D = nc . p_plane (centroid here) + p。輸出表面曲率curvature通過協(xié)方差矩陣的特征值之間的運(yùn)算估計(jì)得到，如：

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

使用OpenMP加速法線估計(jì)

對于對運(yùn)算速度有要求的用戶，PCL點(diǎn)云庫提供了一個(gè)表面法線的附加實(shí)現(xiàn)程序，它使用多核/多線程開發(fā)規(guī)范，利用OpenMP來提高計(jì)算速度。它的類命名為pcl::NormaleEstimationOMP，并且它的應(yīng)用程序接口（API）100%兼容單線程pcl::NormalEstimation，這使它適合作為一個(gè)可選提速方法。在8核系統(tǒng)中，可以輕松提速6-8倍。

-------------------曲率求解-----------------

曲率大小反映模型表面的凹凸程度。

利用上面法向量估計(jì)的PCA方法，在法向量估算的基礎(chǔ)上進(jìn)行數(shù)據(jù)點(diǎn)的曲率估算，曲率的估算方法如上式：

入描述了曲面沿法向量的變化，而和表示數(shù)據(jù)點(diǎn)在切平面上的分布情況。定義下式為數(shù)據(jù)點(diǎn)，在k鄰域內(nèi)的曲面變分。

? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??

點(diǎn)云模型在數(shù)據(jù)點(diǎn)的曲率可近似為在該點(diǎn)的曲面變分,即。

該式的基本思想實(shí)際上是使用曲面變分來近似曲率信息[1]。

曲率基礎(chǔ)知識：
平均曲率、主曲率和高斯曲率是曲率的三個(gè)基本要素。
法曲率：曲面在一點(diǎn)沿著不同方向的彎曲程度不同。或者說曲面離開切平面的速度不同。這個(gè)彎曲屬性可以用這一點(diǎn)的沿著這個(gè)方法的法曲率刻畫
主曲率：過曲面上某個(gè)點(diǎn)上具有無窮個(gè)正交曲率，其中存在一條曲線使得該曲線的曲率為極大，這個(gè)曲率為極大值Kmax，垂直于極大曲率面的曲率為極小值Kmin。這兩個(gè)曲率屬性為主曲率。他們代表著法曲率的極值。
平均曲率：是空間上曲面上某一點(diǎn)任意兩個(gè)相互垂直的正交曲率的平均值。如果一組相互垂直的正交曲率可表示為K1,K2,那么平均曲率則為：K = (K1 +K2 ) / 2。
高斯曲率：兩個(gè)主曲率的乘積即為高斯曲率，又稱總曲率，反映某點(diǎn)上總的完全程度。

#include <iostream>
#include <vector>
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
#include <pcl/features/normal_3d.h>
#include <pcl/features/principal_curvatures.h>
using namespace std;int main (int argc, char** argv){pcl::PointCloud<pcl::PointXYZ>::Ptr cloud (new pcl::PointCloud<pcl::PointXYZ>);pcl::io::loadPCDFile<pcl::PointXYZ>("A3  - Cloud.pcd", *cloud); //讀取點(diǎn)云cout << "Loaded " << cloud->points.size() << " points." << std::endl;//顯示讀取點(diǎn)云的點(diǎn)數(shù)// 計(jì)算點(diǎn)云的法線pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> normalEstimation;normalEstimation.setInputCloud (cloud);//設(shè)置鄰域點(diǎn)搜索方式pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>);normalEstimation.setSearchMethod (tree);//定義一個(gè)新的點(diǎn)云儲(chǔ)存含有法線的值pcl::PointCloud<pcl::Normal>::Ptr cloudWithNormals (new pcl::PointCloud<pcl::Normal>);//設(shè)置KD樹搜索半徑// normalEstimation.setRadiusSearch (0.03);normalEstimation.setKSearch(10);//計(jì)算出來法線的值normalEstimation.compute (*cloudWithNormals);// 建立主曲率計(jì)算pcl::PrincipalCurvaturesEstimation<pcl::PointXYZ, pcl::Normal, pcl::PrincipalCurvatures> principalCurvaturesEstimation;// 提供原始點(diǎn)云(沒有法線)principalCurvaturesEstimation.setInputCloud (cloud);// 為點(diǎn)云提供法線principalCurvaturesEstimation.setInputNormals(cloudWithNormals);// 使用與法線估算相同的KdTreeprincipalCurvaturesEstimation.setSearchMethod (tree);//principalCurvaturesEstimation.setRadiusSearch(1.0);principalCurvaturesEstimation.setKSearch(10);// 計(jì)算主曲率pcl::PointCloud<pcl::PrincipalCurvatures>::Ptr principalCurvatures (new pcl::PointCloud<pcl::PrincipalCurvatures> ());principalCurvaturesEstimation.compute (*principalCurvatures);std::cout << "output points.size: " << principalCurvatures->points.size () << std::endl;// 顯示和檢索第0點(diǎn)的主曲率。cout << "最大曲率;"<< principalCurvatures->points[0].pc1 << endl;//輸出最大曲率cout << "最小曲率:"<< principalCurvatures->points[0].pc2 << endl;//輸出最小曲率//輸出主曲率方向（最大特征值對應(yīng)的特征向量）cout << "主曲率方向;" << endl;cout << principalCurvatures->points[0].principal_curvature_x << endl;cout << principalCurvatures->points[0].principal_curvature_y << endl;cout << principalCurvatures->points[0].principal_curvature_z << endl;return 0;
}

[1]
?

總結(jié)

以上是生活随笔為你收集整理的PCL：从法线计算到曲率计算并可视化的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。