点云：点云格网化 （师弟师妹）简单教学版

一、原理与方法

点云格网化是点云数据处理的常用方法，是一种降维或压缩的方法。（论文里烂大街了）

基本原理是将点放入不同的格网，进而以格网为基本操作单元，对点云进行分块、抽稀、语义分割等。二维格网化是投影，类似像素；三维格网是体素。

基本方法是根据范围，划分成规则或不规则的格网，将不同的点分别放入不同的格网。下面是规则格网的划分。

式中，Xmax、Xmin、Ymax、Ymin、Zmax、Zmin分别代表点云区域X、Y、Z坐标的最大值和最小值； GSD为格网分辨率。

式中，col、row、lay分别为当前点P(X, Y, Z)的格网列号、行号和层号。

二、分享给有需要的人，代码质量勿喷

下面是点云格网化与体素抽稀的C++（老）代码了（有用到Qt）。

2.1 格网化：用PDAL读点，存入 QMultiHash 中

//格网化
QMultiHash<uint, xjPoint> PointCloudGrid::GriddingPC(const QString & lasPath, const double & GSD, const bool & bVoxel)
{
	QMultiHash<uint, xjPoint> mhGrid;
 
#pragma region read las
	using namespace pdal;
	using namespace pdal::Dimension;
	pdal::Option las_opt("filename", lasPath.toStdString());//参数1："filename"（键）
	pdal::Options las_opts;
	las_opts.add(las_opt);
 
	pdal::PointTable rTable;
	pdal::LasReader las_reader;
	las_reader.setOptions(las_opts);
	las_reader.prepare(rTable);
 
	pdal::PointViewSet point_view_set = las_reader.execute(rTable);
	pdal::PointViewPtr point_view = *point_view_set.begin();
	pdal::Dimension::IdList dims = point_view->dims();
	pdal::LasHeader las_header = las_reader.header();
#pragma endregion
 
	//X-cols  Y-rows  Z-lays
	double xMax = las_header.maxX();
	double xMin = las_header.minX();
	double yMax = las_header.maxY();
	double yMin = las_header.minY();
	double zMax = las_header.maxZ();
	double zMin = las_header.minZ();
	uint cols = uint((xMax - xMin) / GSD) + 1;
	uint rows = uint((yMax - yMin) / GSD) + 1;
	uint lays = uint((zMax - zMin) / GSD) + 1;
 
	//traversal
	uint col = 0, row = 0, lay = 0, gridNumber = 0;
	for (int idx = 0; idx < point_view->size(); ++idx)
	{
		xjPoint p;
		p.x = point_view->getFieldAs<double>(Id::X, idx);
		p.y = point_view->getFieldAs<double>(Id::Y, idx);
		p.z = point_view->getFieldAs<double>(Id::Z, idx);
 
		if (las_header.hasColor())
		{
			p.red   = point_view->getFieldAs<int>(Id::Red,   idx);
			p.green = point_view->getFieldAs<int>(Id::Green, idx);
			p.blue  = point_view->getFieldAs<int>(Id::Blue,  idx);
		}
		if (las_header.hasTime())
		{
			p.GPStime = point_view->getFieldAs<double>(Id::GpsTime, idx);
		}
 
		p.intensity = point_view->getFieldAs<int>(Id::Intensity, idx);
		p.pointSourceID = point_view->getFieldAs<int>(Id::PointSourceId, idx);
		p.classification = point_view->getFieldAs<int>(Id::Classification, idx);
		p.userData = point_view->getFieldAs<int>(Id::UserData, idx);
 
		//compute r c l
		row=uint((p.y - yMin) / GSD) + 1;
		col=uint((p.x - xMin) / GSD) + 1;
		lay=uint((p.z - zMin) / GSD) + 1;
		gridNumber = (row - 1)*cols + col;
		if (bVoxel)
			gridNumber = (lay - 1)*rows*cols + (row - 1)*cols + col;
 
		//store
		mhGrid.insert(gridNumber, p);
	}
 
	return mhGrid;
}

2.2 体素抽稀：格网内点坐标取均值，用PDAL写点

//降采样-抽稀：均值
void PointCloudGrid::VoxelDownSample(const QMultiHash<uint, xjPoint> &mhGrid, const QString &resultPath)
{
#pragma region write 1as
	using namespace pdal;
	using namespace pdal::Dimension;
 
	double xoffset = 0, yoffset = 0, zoffset = 0;
 
	PointTable table;
	table.layout()->registerDim(Dimension::Id::X);
	table.layout()->registerDim(Dimension::Id::Y);
	table.layout()->registerDim(Dimension::Id::Z);
	table.layout()->registerDim(Dimension::Id::Red);
	table.layout()->registerDim(Dimension::Id::Green);
	table.layout()->registerDim(Dimension::Id::Blue);
 
	PointViewPtr view(new PointView(table));
#pragma endregion
 
	//traversal
	int idx = 0;
	QList<xjPoint> listP;
	QMultiHash<uint, int> mhKey;
	for (QMultiHash<uint, xjPoint>::const_iterator it = mhGrid.constBegin(); it != mhGrid.constEnd(); ++it)
	{
		uint gn = it.key();
		if (mhKey.contains(gn)) { continue; }
		mhKey.insert(gn, gn);
 
		//average
		double nx = 0, ny = 0, nz = 0;
		int nr = 0, ng = 0, nb = 0;
		listP = mhGrid.values(it.key());
		for each (xjPoint p in listP)
		{
			nx += p.x;
			ny += p.y;
			nz += p.z;
			nr += p.red;
			ng += p.green;
			nb += p.blue;
		}
		nx /= listP.size();
		ny /= listP.size();
		nz /= listP.size();
		nr /= listP.size();
		ng /= listP.size();
		nb /= listP.size();
 
		//store
		view->setField(Id::X, idx, nx);
		view->setField(Id::Y, idx, ny);
		view->setField(Id::Z, idx, nz);
		view->setField(Id::Red,   idx, static_cast<uint16_t>(nr));
		view->setField(Id::Green, idx, static_cast<uint16_t>(ng));
		view->setField(Id::Blue,  idx, static_cast<uint16_t>(nb));
 
		idx++;
		xoffset = nx;
		yoffset = ny;
		zoffset = nz;
	}
 
#pragma region write las
	Options xjOptions;
	xjOptions.add("filename", resultPath.toStdString());
	xjOptions.add("offset_x", xoffset);
	xjOptions.add("offset_y", yoffset);
	xjOptions.add("offset_z", zoffset);
	xjOptions.add("scale_x", 0.0001);
	xjOptions.add("scale_y", 0.0001);
	xjOptions.add("scale_z", 0.0001);
 
	BufferReader xjBufferReader;
	xjBufferReader.addView(view);
 
	StageFactory factory;
	Stage *writer = factory.createStage("writers.las");
	writer->setInput(xjBufferReader);
	writer->setOptions(xjOptions);
	writer->prepare(table);
	writer->execute(table);
#pragma endregion
}

三、试验结果