点云投影至xoy平面生成强度图像

一、备注：

        重新更改过一次的版本，由于新电脑还没来得及装VS 配环境，所以没有测过，后续环境配好了会测一下，做个效果图。先贴上代码吧。按道理应该是没什么问题的，有报错欢迎底下留言！

后续会补上 地面点提取  提取之后 离群点剔除 相关代码

一、算法原理

        就很简单 将点云正射投影至xoy平面，设置固定分辨率，构建格网，求格网内的强度极大值，将强度极大值映射为0~255赋给像素值

 二、代码

// opencv
#include <opencv2/opencv.hpp>   
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
// pcl
#include <pcl/io/pcd_io.h>
#include <pcl/point_types.h>
//lablas
#include <liblas/liblas.hpp>

using namespace std;
using namespace cv;

typedef  pcl::PointCloud<pcl::PointXYZI>::Ptr      pcXYZIPtr;
typedef  pcl::PointCloud<pcl::PointXYZI>           pcXYZI;

/*!
* @brief 读取las格式点云数据函数
* param[in] file_name:las文件  
* param[in] point_cloud:点云指针   
*/
bool readLasFile(const string &file_name, const pcXYZIPtr &point_cloud)
{
	std::ifstream ifs;
	ifs.open(file_name, std::ios::in | std::ios::binary);
	if (ifs.bad())
	{
		cout << "Failed to load las File" << endl;
	}
	liblas::ReaderFactory f;
	liblas::Reader reader = f.CreateWithStream(ifs);
	int number_points = reader.GetHeader().GetPointRecordsCount();//获取las数据点的个数
	point_cloud->width = number_points;
	point_cloud->height = 1;
	point_cloud->is_dense = false;
	point_cloud->resize(point_cloud->width * point_cloud->height);

	liblas::Header const &header = reader.GetHeader();
	int i = 0;
	while (reader.ReadNextPoint())
	{
		point_cloud->points[i].x = reader.GetPoint().GetX();
		point_cloud->points[i].y = reader.GetPoint().GetY();
		point_cloud->points[i].z = reader.GetPoint().GetZ();
		point_cloud->points[i].intensity = reader.GetPoint().GetIntensity();
		
		i++;
	}
	// 检查 强度信息是否存在
	if (point_cloud->points[0].intensity == 0) //check if the intensity exsit
	{
		cout << "Warning! Point cloud intensity may not be imported properly, check the scalar field's name.\n" << endl;
	}
	return 1;
}

/*!
* @brief 点云转强度图像函数
* param[in] cloud：点云指针  
* param[in] resolution：点云转图像分辨率  
* param[out] img: 点云生成的强度图像
*/
cv::Mat pointCloud2imgI(const pcXYZIPtr &cloud, double resolution)
{
	float minx, miny, maxx, maxy, mini, maxi;
	minx = miny = mini = FLT_MAX;
	maxx = maxy = maxi = -FLT_MAX;

	// 获取点云的最大最小 x、y点的坐标
	for (int i = 0; i < cloud->points.size(); i++)
	{
		minx = min(minx, cloud->points[i].x);
		miny = min(miny, cloud->points[i].y);
		maxx = max(maxx, cloud->points[i].x);
		maxy = max(maxy, cloud->points[i].y);
        maxi = max(manxi, cloud->points[i].intensity);
	}

	double lx = maxx - minx;					//点云长度
	double ly = maxy - miny;					//点云宽度
	int rows = round(ly / resolution);			//图像高度
	int clos = round(lx / resolution);			//图像宽度

	cv::Mat img = cv::Mat::zeros(rows, clos, CV_8UC3);

	//强度格网矩阵
	vector<vector<float>> vec_grid_intensity;   

	//格网分配空间 及初始化
	vec_grid_intensity.resize(rows);
	//初始化格网
	for (int i = 0; i < rows; i++)
	{
		vec_grid_intensity[i].resize(clos);
	}

	//依次将点压入所在格网
	for (int i = 0; i < cloud->points.size(); i++)
	{
		int m = (maxy - cloud->points[i].y) / resolution;
		int n = (cloud->points[i].x - minx) / resolution;

		if (m > 0 && m < rows && n > 0 && n < clos)
		{
		    // 将格网中的点云的最大强度值作为格网的强度值
			vec_grid_intensity[m][n] = max(cloud->points[i].intensity, vec_grid_intensity[m][n]);	
		}
	}

	//强度拉伸到0-255之间赋值
	for (int i = 0; i < rows; i++)
	{
		uchar * data = img.ptr<uchar>(i);
		for (int j = 0; j < clos; j++)
		{
			if (vec_grid_intensity[i][j] > 0)
			{
				int pixel = (int)(vec_grid_intensity[i][j]  / maxi * 255);
				//各通道像素赋值
				data[3 * j] = pixel;
				data[3 * j + 1] = pixel;
				data[3 * j + 2] = pixel;
			}
		}
	}
	return img;
}

int main()
{
	const string las_path = "D:/实验/点云/光谷大道.las";
	const double resolution = 0.05;

	pcXYZIPtr cloud(new pcXYZI);

	readLasFile(las_path, cloud);
	Mat cloud_imgI = pointCloud2imgI(cloud, resolution);

	imwrite("cloud_img.jpg", cloud_imgI);
}

三、结果

原始点云

强度图像