利用PCL点云下采样实现数据体素化

PCL

PCL(Point Cloud Library) 库集成了针对大体量级别的空间点数据处理所需要的算法和操作，降低了处理相关需求的复杂度，对快速建立点云数据文档和渲染有着很好的作用。

体素化Voxelization

体素化是通过用空间均匀大小的体素网格(voxel grid)来模拟模型或者点云的几何形态的过程，实现模型体素化的方式有很多，比如基于八叉树的三模网格模型体素化，基于GPU并利用渲染管线中fragment shader部分实现的栅格化插值。本章主要讨论用PCL的下采样方式来实现模型点云的体素化，并且实现可视化。

算法

1. 读入点云数据。
2. 获取AABB盒子。
3. 对点云数据进行下采样过滤。（实现体素化）
4. 渲染可视化。

【建立结构体】

typedef struct point_XYZ{
	
		double point_X;
		double point_Y;
		double point_Z;
		double point_R;
};

【读入点云数据】

	//point clouds
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud(new pcl::PointCloud<pcl::PointXYZ>());
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_txt(new pcl::PointCloud<pcl::PointXYZ>());
	pcl::PointCloud<pcl::PointXYZ>::Ptr cloud_filter(new pcl::PointCloud<pcl::PointXYZ>());

	FILE *file;
	int point_Sum; 
	point_XYZ point;
	vector<point_XYZ> vec_Point;
	
	file = fopen(argv[1], "r");
	if(file){
	
		while(fscanf(file, "%lf %lf %lf", &point.point_X, &point.point_Y, &point.point_Z) != EOF){
		
			vec_Point.push_back(point);
		}
	}
	else{
	
		cout << "Could not load data from " << argv[1] <<" file..." <<endl;
		return 0;
	}
	
	fclose(file);
	point_Sum = vec_Point.size();
	
	for(size_t i = 0; i< cloud -> points.size(); ++i){
	
		cloud -> points[i].x = vec_Point[i].point_X;
		cloud -> points[i].y = vec_Point[i].point_Y;
		cloud -> points[i].z = vec_Point[i].point_Z;
	}

或者直接从PCD文件中读入点云信息

	if(pcl::io::loadPCDFile(pcdName, *cloud) == -1){
	
		cout << "Could not open the file: " << pcdName << "..." <<endl;
		return -1;
	}

【获取AABB】

	pcl::MomentOfInertiaEstimation <pcl::PointXYZ> feature_extractor;
	feature_extractor.setInputCloud(cloud);
	feature_extractor.compute();
	
	pcl::PointXYZ min_point_AABB;
	pcl::PointXYZ max_point_AABB;
	
	feature_extractor.getAABB(min_point_AABB, max_point_AABB);

【下采样体素化】

	//VoxelGrid filtering ***PCL***
	pcl::VoxelGrid<pcl::PointXYZ> fil;
	fil.setInputCloud(cloud);
	fil.setLeafSize(atof(argv[2]), atof(argv[2]), atof(argv[2]));
	fil.filter(*cloud_filter);
	
	int filterPointSize = cloud_filter -> points.size();

下采样的过程就是把点云放到规定大小的单位体素盒中去，如果一个局部的点群都落在了一个单位体素盒子里，那么这个盒子中所有的点都会用一个重心点来表示。通过这个过程，点云数据会大大被过滤，但是模型的体态还保持着完整性。

如果不用PCL库也可以实现，下面的核心带面可供参考：

	int ijk0 = static_cast<int> (ceil(vec_Point[i].point_X * inverse_voxel_size.index_X) - (bound_min.point_X * inverse_voxel_size.index_X));
	int ijk1 = static_cast<int> (floor(vec_Point[i].point_Y * inverse_voxel_size.index_Y) - (bound_min.point_Y * inverse_voxel_size.index_Y));
	int ijk2 = static_cast<int> (floor(vec_Point[i].point_Z * inverse_voxel_size.index_Z) - (bound_min.point_Z * inverse_voxel_size.index_Z));

如上求出每个点对应的x，y，z后，就可以对每个单位体素求其中点的重心点，并且把重心点插入到每个单位体素中。

【可视化渲染】

	pcl::visualization::PCLVisualizer::Ptr viewer (new pcl::visualization::PCLVisualizer("Voxelization"));
	viewer -> setBackgroundColor(0,0,0);
	
	//Add AABB box
	viewer -> addCube (min_point_AABB.x, max_point_AABB.x, min_point_AABB.y, max_point_AABB.y, min_point_AABB.z, max_point_AABB.z,
		1.0, 1.0, 1.0, "AABB");
	viewer -> setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION, pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, "AABB");

两种体素化渲染方式：

1. Solid cube

	for(int i=0; i<filterPointSize; i++){
	
			double x = cloud_filter->points[i].x;
			double y = cloud_filter->points[i].y;
			double z = cloud_filter->points[i].z;
		
			Eigen::Vector3f center(floor(x / voxelSize)*voxelSize + voxelSize/2, floor(y / voxelSize)*voxelSize + voxelSize/2, floor(z / voxelSize)*voxelSize + voxelSize/2);
		
			Eigen::Quaternionf rotation(1,0,0,0);
			string cube = "AABB"+to_string(i);
			viewer -> addCube(center,rotation,voxelSize, voxelSize, voxelSize, cube);
			
	}

1. Wireframe cube

	//rendering by z axis
	pcl::visualization::PointCloudColorHandlerGenericField<pcl::PointXYZ> v_color(cloud,"z"); 
	viewer -> addPointCloud<pcl::PointXYZ>(cloud, v_color, "vertices color");
			
	for(int i=0; i<filterPointSize; i++){
	
			double x = cloud_filter->points[i].x;
			double y = cloud_filter->points[i].y;
			double z = cloud_filter->points[i].z;
		
			Eigen::Vector3f center(floor(x / voxelSize)*voxelSize + voxelSize/2, floor(y / voxelSize)*voxelSize + voxelSize/2, floor(z / voxelSize)*voxelSize + voxelSize/2);
		
			Eigen::Quaternionf rotation(1,0,0,0);
			string cube = "AABB"+to_string(i);
			viewer -> addCube(center,rotation,voxelSize, voxelSize, voxelSize, cube);
			
			//shape rendering
			viewer -> setShapeRenderingProperties(pcl::visualization::PCL_VISUALIZER_REPRESENTATION, pcl::visualization::PCL_VISUALIZER_REPRESENTATION_WIREFRAME, cube);
	}

【示例图】