---

---

前言

Delaunay三角剖分最初應用于2維領域，而與Greedy三角化算法的結合，使之成為目前在三維重建領域最為基礎的算法原理之一，很多學者針對其原理進行改進用以三維點云模型的構建。

---

一、PCL點云三角化

1.1 Delaunay三角剖分

定義：假設點集中的一條邊e（兩個端點為a,b），e若滿足下列條件，則稱之為Delaunay邊：存在一個圓經過a,b兩點，圓內(圓上最多三點共圓)不含點集中任何其他的點。而Delaunay三角化就是指三角網格均是由Delaunay邊組成，并滿足最小角最大原則（在點集可能形成的三角剖分中，Delaunay三角剖分所形成的三角形的最小角最大）。

針對以上定義，目前已提出了很多經典的剖分算法，如Lawson算法、Bowyer-Watson算法。以上算法都很有意思，通過點插法實現，具體原理可以查看以下鏈接。
技術分享：Delaunay三角剖分算法介紹

1.2 貪婪三角化

PCL中采用將三維點云投影到二維平面的方法來實現三角剖分， 具體采用貪婪三角化算法。
其過程為：
1：計算點云中點的法線，再將點云通過法線投影到二維坐標平面。
2：使用基于Delaunay三角剖分的空間區域增長算法完成平面點集的三角化。
3：根據投影點云的連接關系確定原始三維點云間的拓撲關系，最終得到曲面模型。

PCL中的NormalEstimation和GreedyProjectionTriangulation類實現該計算過程。
源代碼：

FFN和SFN是指兩個不同方向的邊緣鄰域集，在connectPoint方法里完成計算。

/** \brief Index of the current query point **/
int R_;std::vector<int> ffn_;
std::vector<int> sfn_;

// Locating FFN and SFN to adapt distance thresholddouble sqr_source_dist = (coords_[R_] - coords_[source_[R_]]).squaredNorm ();double sqr_ffn_dist = (coords_[R_] - coords_[ffn_[R_]]).squaredNorm ();double sqr_sfn_dist = (coords_[R_] - coords_[sfn_[R_]]).squaredNorm ();double max_sqr_fn_dist = (std::max)(sqr_ffn_dist, sqr_sfn_dist);double sqr_dist_threshold = (std::min)(sqr_max_edge, sqr_mu * sqrDists[1]); //sqr_mu * sqr_avg_conn_dist);if (max_sqr_fn_dist > sqrDists[nnn_-1]){if (0 == increase_nnn4fn)PCL_WARN("Not enough neighbors are considered: ffn or sfn out of range! Consider increasing nnn_... Setting R=%d to be BOUNDARY!\n", R_);increase_nnn4fn++;state_[R_] = BOUNDARY;continue;}double max_sqr_fns_dist = (std::max)(sqr_source_dist, max_sqr_fn_dist);if (max_sqr_fns_dist > sqrDists[nnn_-1]){if (0 == increase_nnn4s)PCL_WARN("Not enough neighbors are considered: source of R=%d is out of range! Consider increasing nnn_...\n", R_);increase_nnn4s++;}

計算法線：

// Get the normal estimate at the current point 
const Eigen::Vector3f nc = (*input_)[(*indices_)[R_]].getNormalVector3fMap ();

三角化：

// Triangulatingif (angles_[2].visible == false){if ( !( (angles_[0].index == ffn_[R_] && angles_[1].index == sfn_[R_]) || (angles_[0].index == sfn_[R_] && angles_[1].index == ffn_[R_]) ) ){state_[R_] = BOUNDARY;}else{if ((source_[R_] == angles_[0].index) || (source_[R_] == angles_[1].index))state_[R_] = BOUNDARY;else{if (sqr_max_edge < (coords_[ffn_[R_]] - coords_[sfn_[R_]]).squaredNorm ()){state_[R_] = BOUNDARY;}else{tmp_ = coords_[source_[R_]] - proj_qp_;uvn_s[0] = tmp_.dot(u_);uvn_s[1] = tmp_.dot(v_);double angleS = std::atan2(uvn_s[1], uvn_s[0]);double dif = angles_[1].angle - angles_[0].angle;if ((angles_[0].angle < angleS) && (angleS < angles_[1].angle)){if (dif < 2*M_PI - maximum_angle_)state_[R_] = BOUNDARY;elsecloseTriangle (polygons);}else{if (dif >= maximum_angle_)state_[R_] = BOUNDARY;elsecloseTriangle (polygons);}}}}continue;}

源碼中大量代碼關注于三角形的連接問題。
最后調用MeshConstruction類的reconstruct方法進行表面重建。

template <typename PointInT> void
MeshConstruction<PointInT>::reconstruct (std::vector<pcl::Vertices> &polygons)
{if (!initCompute ()){polygons.clear ();return;}// Check if a space search locator was givenif (check_tree_){if (!tree_){if (input_->isOrganized ())tree_.reset (new pcl::search::OrganizedNeighbor<PointInT> ());elsetree_.reset (new pcl::search::KdTree<PointInT> (false));}// Send the surface dataset to the spatial locatortree_->setInputCloud (input_, indices_);}// Set up the output dataset//polygons.clear ();//polygons.reserve (2 * indices_->size ()); /// NOTE: usually the number of triangles is around twice the number of vertices// Perform the actual surface reconstructionperformReconstruction (polygons);deinitCompute ();
}

二、程序示例

//----------------------------------法線計算-----------------------------------pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> n;pcl::PointCloud<pcl::Normal>::Ptr normals(new pcl::PointCloud<pcl::Normal>);pcl::search::KdTree<pcl::PointXYZ>::Ptr tree1(new pcl::search::KdTree<pcl::PointXYZ>);tree->setInputCloud(cloud);n.setInputCloud(cloud);n.setSearchMethod(tree1);n.setKSearch(20);n.compute(*normals);pcl::PointCloud<pcl::PointNormal>::Ptr cloud_with_normals(new pcl::PointCloud<pcl::PointNormal>);pcl::concatenateFields(*cloud, *normals, *cloud_with_normals);//連接點云和法線pcl::search::KdTree<pcl::PointNormal>::Ptr tree2(new pcl::search::KdTree<pcl::PointNormal>);tree2->setInputCloud(cloud_with_normals);pcl::GreedyProjectionTriangulation<pcl::PointNormal> gp;pcl::PolygonMesh triangles;gp.setSearchRadius(0.025);//設置搜索半徑，即連接點的最大距離gp.setMu(2.5); //加權因子，對于每個樣本點，其映射所選球的半徑由mu與離樣本點最近點的距離乘積決定，用以解決點云密度不均勻的問題，mu一般取值2.5-3gp.setMaximumNearestNeighbors(600); //最大領域點個數gp.setMaximumSurfaceAngle(M_PI / 4);//臨近點的法線和樣本點法線的最大偏離角度gp.setMinimumAngle(M_PI / 18); //三角形最小角gp.setMaximumAngle(2 * M_PI / 3);//三角形最大角gp.setNormalConsistency(false); //保證法線朝向一致gp.setInputCloud(cloud_with_normals);gp.setSearchMethod(tree2);   gp.reconstruct(triangles);

總結

徹底理解三角化的源代碼比較困難，缺少相關學習資料，歡迎共同研究，提出意見。