EuclideanClusterExtraction這個名字起的很奇怪,歐式距離聚類這個該如何理解?歐式距離只是一種距離測度的方法呀!有了一個Cluster在里面,我以為是某一種聚類算法,層次聚類?k-NN聚類?K-Means?還是模糊聚類?感覺很奇怪,看下代碼吧。
找一個實例cluster_extraction.cpp的main入口函數。
找到computer函數,該方法中定義了一個pcl::EuclideanClusterExtraction<pcl::PointXYZ> ec;對象,接着就是參數賦值。
具體的算法執行在ec.extract (cluster_indices)中。因此進入EuclideanClusterExtraction查看具體的方法。
1 void compute (const pcl::PCLPointCloud2::ConstPtr &input, std::vector<pcl::PCLPointCloud2::Ptr> &output, 2 int min, int max, double tolerance) 3 { 4 // Convert data to PointCloud<T> 5 PointCloud<pcl::PointXYZ>::Ptr xyz (new PointCloud<pcl::PointXYZ>); 6 fromPCLPointCloud2 (*input, *xyz); 7 8 // Estimate 9 TicToc tt; 10 tt.tic (); 11 12 print_highlight (stderr, "Computing "); 13 14 // Creating the KdTree object for the search method of the extraction 15 pcl::search::KdTree<pcl::PointXYZ>::Ptr tree (new pcl::search::KdTree<pcl::PointXYZ>); 16 tree->setInputCloud (xyz); 17 18 std::vector<pcl::PointIndices> cluster_indices; 19 pcl::EuclideanClusterExtraction<pcl::PointXYZ> ec; 20 ec.setClusterTolerance (tolerance); 21 ec.setMinClusterSize (min); 22 ec.setMaxClusterSize (max); 23 ec.setSearchMethod (tree); 24 ec.setInputCloud (xyz); 25 ec.extract (cluster_indices); 26 27 print_info ("[done, "); print_value ("%g", tt.toc ()); print_info (" ms : "); print_value ("%d", cluster_indices.size ()); print_info (" clusters]\n"); 28 29 output.reserve (cluster_indices.size ()); 30 for (std::vector<pcl::PointIndices>::const_iterator it = cluster_indices.begin (); it != cluster_indices.end (); ++it) 31 { 32 pcl::ExtractIndices<pcl::PCLPointCloud2> extract; 33 extract.setInputCloud (input); 34 extract.setIndices (boost::make_shared<const pcl::PointIndices> (*it)); 35 pcl::PCLPointCloud2::Ptr out (new pcl::PCLPointCloud2); 36 extract.filter (*out); 37 output.push_back (out); 38 } 39 }
可以在pcl_segmentation項目下的extract_clusters.h文件中查看EuclideanClusterExtraction的定義,可知這是一個模板類。
1 template <typename PointT> 2 class EuclideanClusterExtraction: public PCLBase<PointT>
實現文件在項目下的extract_clusters.hpp中,(還有一個extract_clusters.cpp文件),查看其中的extract方法:
1 template <typename PointT> void pcl::EuclideanClusterExtraction<PointT>::extract (std::vector<PointIndices> &clusters) 2 { 3 if (!initCompute () || 4 (input_ != 0 && input_->points.empty ()) || 5 (indices_ != 0 && indices_->empty ())) 6 { 7 clusters.clear (); 8 return; 9 } 10 11 // Initialize the spatial locator 12 if (!tree_) 13 { 14 if (input_->isOrganized ()) 15 tree_.reset (new pcl::search::OrganizedNeighbor<PointT> ()); 16 else 17 tree_.reset (new pcl::search::KdTree<PointT> (false)); 18 } 19 20 // Send the input dataset to the spatial locator 21 tree_->setInputCloud (input_, indices_); 22 extractEuclideanClusters (*input_, *indices_, tree_, static_cast<float> (cluster_tolerance_), clusters, min_pts_per_cluster_, max_pts_per_cluster_); 23 24 //tree_->setInputCloud (input_); 25 //extractEuclideanClusters (*input_, tree_, cluster_tolerance_, clusters, min_pts_per_cluster_, max_pts_per_cluster_); 26 27 // Sort the clusters based on their size (largest one first) 28 std::sort (clusters.rbegin (), clusters.rend (), comparePointClusters); 29 30 deinitCompute (); 31 }
注意在extract_clusters.h中定義了四個
1 template <typename PointT> void 2 extractEuclideanClusters ( 3 const PointCloud<PointT> &cloud, const boost::shared_ptr<search::Search<PointT> > &tree, 4 float tolerance, std::vector<PointIndices> &clusters, 5 unsigned int min_pts_per_cluster = 1, unsigned int max_pts_per_cluster = (std::numeric_limits<int>::max) ());
1 template <typename PointT> void 2 extractEuclideanClusters ( 3 const PointCloud<PointT> &cloud, const std::vector<int> &indices, 4 const boost::shared_ptr<search::Search<PointT> > &tree, float tolerance, std::vector<PointIndices> &clusters, 5 unsigned int min_pts_per_cluster = 1, unsigned int max_pts_per_cluster = (std::numeric_limits<int>::max) ());
1 template <typename PointT, typename Normal> void 2 extractEuclideanClusters ( 3 const PointCloud<PointT> &cloud, const PointCloud<Normal> &normals, 4 float tolerance, const boost::shared_ptr<KdTree<PointT> > &tree, 5 std::vector<PointIndices> &clusters, double eps_angle, 6 unsigned int min_pts_per_cluster = 1, 7 unsigned int max_pts_per_cluster = (std::numeric_limits<int>::max) ())
1 template <typename PointT, typename Normal> 2 void extractEuclideanClusters ( 3 const PointCloud<PointT> &cloud, const PointCloud<Normal> &normals, 4 const std::vector<int> &indices, const boost::shared_ptr<KdTree<PointT> > &tree, 5 float tolerance, std::vector<PointIndices> &clusters, double eps_angle, 6 unsigned int min_pts_per_cluster = 1, 7 unsigned int max_pts_per_cluster = (std::numeric_limits<int>::max) ())
前兩個的方法實現在文件extract_clusters.hpp中,后兩個直接在頭文件中就以內聯函數的形式實現了,兩個大同小異。擇其中第一個加點注釋,發現其實是采用的區域生長算法實現的分割。理解錯誤了,區域生長需要種子點,這里應該叫層次聚類方法。
1 template <typename PointT, typename Normal> void 2 extractEuclideanClusters ( 3 const PointCloud<PointT> &cloud, const PointCloud<Normal> &normals, 4 float tolerance, const boost::shared_ptr<KdTree<PointT> > &tree, 5 std::vector<PointIndices> &clusters, double eps_angle, 6 unsigned int min_pts_per_cluster = 1, 7 unsigned int max_pts_per_cluster = (std::numeric_limits<int>::max) ()) 8 { 9 if (tree->getInputCloud ()->points.size () != cloud.points.size ()) 10 { 11 PCL_ERROR ("[pcl::extractEuclideanClusters] Tree built for a different point cloud dataset (%lu) than the input cloud (%lu)!\n", tree->getInputCloud ()->points.size (), cloud.points.size ()); 12 return; 13 } 14 if (cloud.points.size () != normals.points.size ()) 15 { 16 PCL_ERROR ("[pcl::extractEuclideanClusters] Number of points in the input point cloud (%lu) different than normals (%lu)!\n", cloud.points.size (), normals.points.size ()); 17 return; 18 } 19 20 // Create a bool vector of processed point indices, and initialize it to false 21 std::vector<bool> processed (cloud.points.size (), false); 22 23 std::vector<int> nn_indices; 24 std::vector<float> nn_distances; 25 // Process all points in the indices vector 26 for (size_t i = 0; i < cloud.points.size (); ++i)//遍歷點雲中的每一個點 27 { 28 if (processed[i])//如果該點已經處理則跳過 29 continue; 30 31 std::vector<unsigned int> seed_queue;//定義一個種子隊列 32 int sq_idx = 0; 33 seed_queue.push_back (static_cast<int> (i));//加入一個種子 34 35 processed[i] = true; 36 37 while (sq_idx < static_cast<int> (seed_queue.size ()))//遍歷每一個種子 38 { 39 // Search for sq_idx 40 if (!tree->radiusSearch (seed_queue[sq_idx], tolerance, nn_indices, nn_distances))//沒找到近鄰點就繼續 41 { 42 sq_idx++; 43 continue; 44 } 45 46 for (size_t j = 1; j < nn_indices.size (); ++j) // nn_indices[0] should be sq_idx 47 { 48 if (processed[nn_indices[j]]) // Has this point been processed before ?種子點的近鄰點中如果已經處理就跳出此次循環繼續 49 continue; 50 51 //processed[nn_indices[j]] = true; 52 // [-1;1] 53 double dot_p = normals.points[i].normal[0] * normals.points[nn_indices[j]].normal[0] + 54 normals.points[i].normal[1] * normals.points[nn_indices[j]].normal[1] + 55 normals.points[i].normal[2] * normals.points[nn_indices[j]].normal[2]; 56 if ( fabs (acos (dot_p)) < eps_angle ) //根據內積求夾角,法向量都是單位向量,種子點和近鄰點的法向量夾角小於閾值, 57 { 58 processed[nn_indices[j]] = true; 59 seed_queue.push_back (nn_indices[j]);//將此種子點的臨近點作為新的種子點。 60 } 61 } 62 63 sq_idx++; 64 } 65 66 // If this queue is satisfactory, add to the clusters 67 if (seed_queue.size () >= min_pts_per_cluster && seed_queue.size () <= max_pts_per_cluster) 68 { 69 pcl::PointIndices r; 70 r.indices.resize (seed_queue.size ()); 71 for (size_t j = 0; j < seed_queue.size (); ++j) 72 r.indices[j] = seed_queue[j]; 73 74 // These two lines should not be needed: (can anyone confirm?) -FF 75 std::sort (r.indices.begin (), r.indices.end ()); 76 r.indices.erase (std::unique (r.indices.begin (), r.indices.end ()), r.indices.end ()); 77 78 r.header = cloud.header; 79 clusters.push_back (r); // We could avoid a copy by working directly in the vector 80 } 81 } 82 }
