CGAL 使用PCA计算点云法向量

一、算法原理

见：PCL 计算点云法向量并显示

1、主要函数

头文件

include <CGAL/pca_estimate_normals.h>

函数

void CGAL::pca_estimate_normals  ( PointRange &  points,  
  unsigned int  k,  
  const NamedParameters &  np = parameters::default_values()  
 ) 

  通过最近邻平面上的线性最小二乘拟合估计点范围的法方向。输出法线的方向是随机的。

• points：点云
• k： 邻域点数。
• np：下面列出的命名参数中的一个可选序列。
  

二、代码实现

#include <utility> // std::pair
#include <list>
#include <fstream>
#include <CGAL/property_map.h>      
#include <CGAL/IO/read_points.h>           // 读取点云
#include <CGAL/IO/write_points.h>          // 保存点云
#include <CGAL/pca_estimate_normals.h>    // pca计算法向量
#include <CGAL/compute_average_spacing.h> // 计算点云平均密度
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
// Types
typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;

// 定义存储点和法线的容器
typedef std::pair<Kernel::Point_3, Kernel::Vector_3> PointVectorPair;

int main(int argc, char* argv[])
{
	const std::string fname = CGAL::data_file_path("cgal//sphere_1k.xyz");
	const char* output_filename = "cgal//pca_normal.xyz";
	// -----------------------------------读取点云------------------------------------
	std::list<PointVectorPair> points;
	if (!CGAL::IO::read_points(fname, std::back_inserter(points),
		CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>())))
	{
		std::cerr << "点云读取失败！！！ " << fname << std::endl;
		return -1;
	}
	bool knn_search = true;
	const int nb_neighbors = 18; // K近邻搜索的邻域点数
	if (knn_search == true)
	{
		// ---------------------使用固定近邻点来计算法线-----------------------------
		///注意:pca_estimate_normals()需要一个范围的点以及属性映射来访问每个点的位置和法线。

		CGAL::pca_estimate_normals<CGAL::Parallel_if_available_tag>(points, nb_neighbors,
			CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>())
			.normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>()));
	}
	else
	{
		double spacing = CGAL::compute_average_spacing<CGAL::Parallel_if_available_tag>(points, nb_neighbors,
			CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>()));
		// ---------------------使用固定半径进行法线计算-----------------------------
		CGAL::pca_estimate_normals<CGAL::Parallel_if_available_tag>
			(points,
				0, // 当使用邻域半径时，K=0表示对返回的邻域数量没有限制
				CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>())
				.normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>())
				.neighbor_radius(2. * spacing)); // 使用平均间距的2倍作为搜索半径

	}
	// ---------------------------------结果保存--------------------------------------
	if (!CGAL::IO::write_XYZ(output_filename, points,
		CGAL::parameters::point_map(CGAL::First_of_pair_property_map<PointVectorPair>())
		.normal_map(CGAL::Second_of_pair_property_map<PointVectorPair>())
		.stream_precision(17)))
		return -1;

	return 0;
}