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

一、原理与方法

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

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

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

式中，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 lasusing 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-laysdouble 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;//traversaluint 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 lrow=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;//storemhGrid.insert(gridNumber, p);}return mhGrid;
}

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

//降采样-抽稀：均值
void PointCloudGrid::VoxelDownSample(const QMultiHash<uint, xjPoint> &mhGrid, const QString &resultPath)
{
#pragma region write 1asusing 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//traversalint 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);//averagedouble 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();//storeview->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 lasOptions 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
}

三、试验结果