一、需求

将las格式的点云数据转换为cesium能加载渲染的3dtiles格式(pnts)。

二、具体实现步骤

2.1 依赖开源库

gdal库：点云数据多为投影坐标系，cesium渲染需要的是地理坐标系，所以需要gdal库进行坐标转换。

LASzip库：读写las的库。

Qt 、osg等库不是必须，可以不用。我的代码其他地方需要这两个库，所以就用上了。

2.2点云读取

使用LASzip库读取点云到内存。这里在具体程序实现的时候，需要注意分块策略，因为有些点云数据动则几十个G，如果不进行分块处理，对硬件压力太大了。具体的分块思路应该是现对点云数据构建索引，我代码中用的四叉树构建的点云索引；然后根据点云数量及范围对点云进行分块，每次处理的时候针对某一块进行处理。而且分块后方便后续并行加速处理。

2.3 lod构建

第一步：根据2.2中的分块结果，读取其中的某一块点云数据，统计其点个数和外界包围盒。

第二步：判断统计的点云数据个数是否大于指定的阈值(该值是3dtiles中每个pnts文件的最大点数量)。如果小于指定阈值，则该块不需要进行分块，将该块中的全部点写入pnts。如果小于指定阈值，则以该块外界包围盒(矩形)中心为分界点，创建四个子节点块，该块则被称为根节点。

第三步：针对第二步中得到的根节点块和子节点块，继续按照步骤二中的规则对四个子节点进行分块。具体是否分块还是按照节点中的点云数与指定阈值的大小来决定。以此循环，直到所有的节点都不能再分块为止。这样就完成了整个点云的lod构建。

上述方法描述起来比较拗口，其实用程序实现的时候，只需要进行递归调用即可实现，部分实现代码如下：

QSharedPointer<scially::OSGIndexNode> Build3dtiles(const std::vector<PointCI>* pointSet, std::vector<unsigned int>& pointIndex, osg::BoundingBox boundingBox, osg::BoundingBox boundingBoxLevel0, const std::string& saveFilePath, const std::string& strBlock, unsigned int level, unsigned int childNo, ExportMode exportMode, const scially::SpatialTransform& transform, const scially::TileStorage& storage) { // filename of self, left, right std::string saveFileName, pageName1, pageName2, pageName3, pageName4; if (level == 0) { saveFileName = strBlock; } else { char tmpSaveFileName[100]; sprintf(tmpSaveFileName, "%s%s%d%s%d", strBlock.c_str(), "_L", level, "_", childNo); saveFileName.assign(tmpSaveFileName); saveFileName = /*saveFilePath + "/" +*/ saveFileName; } //********** 1 2 ********** //********** 3 4 ********** char tmpPageName1[100], tmpPageName2[100], tmpPageName3[100], tmpPageName4[100]; sprintf(tmpPageName1, "%s%s%d%s%d", strBlock.c_str(), "_L", level + 1, "_", childNo * 4); pageName1.assign(tmpPageName1); sprintf(tmpPageName2, "%s%s%d%s%d", strBlock.c_str(), "_L", level + 1, "_", childNo * 4 + 1); pageName2.assign(tmpPageName2); sprintf(tmpPageName3, "%s%s%d%s%d", strBlock.c_str(), "_L", level + 1, "_", childNo * 4 + 2); pageName3.assign(tmpPageName3); sprintf(tmpPageName4, "%s%s%d%s%d", strBlock.c_str(), "_L", level + 1, "_", childNo * 4 + 3); pageName4.assign(tmpPageName4); // handle leaf case if (pointIndex.size() < _maxPointNumPerOneNode || level >= _maxTreeLevel) { //double rangeValue = boundingBox.radius() / scially::SPLIT_PIXEL / 32; double geometricError = scially::osgBoundingSize(boundingBoxLevel0) / scially::SPLIT_PIXEL / 16 / (pow(2.0, level)); auto pntsNode = QSharedPointer<scially::PntsTile>::create(); pntsNode->tileName() = QString::fromStdString(saveFileName); pntsNode->fileName() = QString::fromStdString(strBlock); pntsNode->geometricError() = geometricError; pntsNode->tileFolder() = QString::fromStdString(saveFilePath); pntsNode->m_suffix = ".pnts"; pntsNode->m_refine = "ADD"; pntsNode->boundingBox() = boundingBox; QString outFileName = pntsNode->relativeNodePath(pntsNode->m_suffix); QByteArray pntsBuffer; osg::Vec3d tileCenter = transform.transform(boundingBox.center()); if (write2Pnts(pointSet, pointIndex, osg::Vec3d(0, 0, 0), pntsBuffer)) { if (!storage.saveFile(outFileName, pntsBuffer)) { qWarning("saveFile pnts failed!"); return nullptr; } return pntsNode; } else { qWarning("write2Pnts failed!"); return nullptr; } } // prepare box AxisInfo midAxisInfo; osg::BoundingBox boundingBox1; osg::BoundingBox boundingBox2; osg::BoundingBox boundingBox3; osg::BoundingBox boundingBox4; midAxisInfo = FindMidAxis(boundingBox, boundingBox1, boundingBox2, boundingBox3, boundingBox4); // split self, child1, child2, child3, child4 float interval = (float)pointIndex.size() / (float)_maxPointNumPerOneNode; int count = -1; std::vector<unsigned int> selfPointSetIndex; std::vector<unsigned int> child1PointSetIndex; std::vector<unsigned int> child2PointSetIndex; std::vector<unsigned int> child3PointSetIndex; std::vector<unsigned int> child4PointSetIndex; for (int i = 0; i < pointIndex.size(); i++) { int tmp = int((float)i / interval); if (tmp > count && selfPointSetIndex.size() < _maxPointNumPerOneNode) { count = tmp; selfPointSetIndex.push_back(pointIndex[i]); } else { PointCI tmpPoint = pointSet->at(pointIndex[i]); if (tmpPoint.P[0]<midAxisInfo.midx && tmpPoint.P[1]>midAxisInfo.midy) { child1PointSetIndex.push_back(pointIndex[i]); } else if (tmpPoint.P[0] > midAxisInfo.midx && tmpPoint.P[1] > midAxisInfo.midy) { child2PointSetIndex.push_back(pointIndex[i]); } else if (tmpPoint.P[0] < midAxisInfo.midx && tmpPoint.P[1] < midAxisInfo.midy) { child3PointSetIndex.push_back(pointIndex[i]); } else { child4PointSetIndex.push_back(pointIndex[i]); } } } // export //double rangeValue = boundingBox.radius() / scially::SPLIT_PIXEL / 32; double geometricError = scially::osgBoundingSize(boundingBoxLevel0) / scially::SPLIT_PIXEL / 16 / (pow(2.0, level)); auto pntsNode = QSharedPointer<scially::PntsTile>::create(); pntsNode->tileName() = QString::fromStdString(saveFileName); pntsNode->fileName() = QString::fromStdString(strBlock); pntsNode->geometricError() = geometricError; pntsNode->tileFolder() = QString::fromLocal8Bit(saveFilePath.c_str()); pntsNode->m_suffix = ".pnts"; pntsNode->m_refine = "ADD"; pntsNode->boundingBox() = boundingBox; QString outFileName = pntsNode->relativeNodePath(pntsNode->m_suffix); QByteArray pntsBuffer; osg::Vec3d tileCenter = transform.transform(boundingBox.center()); if (write2Pnts(pointSet, selfPointSetIndex, osg::Vec3d(0, 0, 0), pntsBuffer)) { if (!storage.saveFile(outFileName, pntsBuffer)) { qWarning("saveFile pnts failed!"); return nullptr; } } else { qWarning("write2Pnts failed!"); return false; } // recursive left if (child1PointSetIndex.size()) { auto n = Build3dtiles(pointSet, child1PointSetIndex, boundingBox1, boundingBoxLevel0, saveFilePath, strBlock, level + 1, childNo * 4, exportMode, transform, storage); if (n) { pntsNode->append(n); } child1PointSetIndex.swap(std::vector<unsigned int>()); } if (child2PointSetIndex.size()) { auto n = Build3dtiles(pointSet, child2PointSetIndex, boundingBox2, boundingBoxLevel0, saveFilePath, strBlock, level + 1, childNo * 4 + 1, exportMode, transform, storage); if (n) { pntsNode->append(n); } child2PointSetIndex.swap(std::vector<unsigned int>()); } if (child3PointSetIndex.size()) { auto n = Build3dtiles(pointSet, child3PointSetIndex, boundingBox3, boundingBoxLevel0, saveFilePath, strBlock, level + 1, childNo * 4 + 2, exportMode, transform, storage); if (n) { pntsNode->append(n); } child3PointSetIndex.swap(std::vector<unsigned int>()); } if (child4PointSetIndex.size()) { auto n = Build3dtiles(pointSet, child4PointSetIndex, boundingBox4, boundingBoxLevel0, saveFilePath, strBlock, level + 1, childNo * 4 + 3, exportMode, transform, storage); if (n) { pntsNode->append(n); } child4PointSetIndex.swap(std::vector<unsigned int>()); } return pntsNode; }

2.4 写pnts

按照2.3构建起点云的lod后，我们得到的是一个包含多个节点的四叉树，四叉树中每个节点代表一块点数量不超过指定阈值的点云数据。我们要将所有的节点写入pnts文件。例如根节点点云数量为s，指定的阈值是t，则我们在写根节点的pnts文件时，需要从s中均匀的抽取出t个点。只需要计算出间隔取样距离即可：float interval = s/t;按interval间隔取点，写入pnts。下面的代码是单个pnts文件的写入实现，完全按照cesium官方给出的pnts格式进行写入的。写颜色的时候，有强度渲染和RGB渲染之分。下面是写pnts代码截取。

//写pnts代码 bool write2Pnts(const std::vector<PointCI>* pointSet, std::vector<unsigned int>& pointIndex, const osg::Vec3d& center, QByteArray& buffer) { if (pointIndex.size() <= 0) { return false; } osg::ref_ptr<osg::Vec3Array> pointArray = new osg::Vec3Array; osg::ref_ptr<osg::Vec4Array> colorArray = new osg::Vec4Array; for (std::vector<unsigned int>::iterator i = pointIndex.begin(); i != pointIndex.end(); i++) { PointCI tmpPoint = pointSet->at(*i); pointArray->push_back(tmpPoint.P); if (_colorMode == ColorMode::Debug) { colorArray->push_back(_colorBar[0]); } else if (_colorMode == ColorMode::RGB) { colorArray->push_back( osg::Vec4(Color8BitsToFloat(tmpPoint.C[0]), Color8BitsToFloat(tmpPoint.C[1]), Color8BitsToFloat(tmpPoint.C[2]), 1.f)); } else if (_colorMode == ColorMode::IntensityGrey) { colorArray->push_back(_colorBar[tmpPoint.I]); } else if (_colorMode == ColorMode::IntensityBlueWhiteRed) { colorArray->push_back(_colorBar[tmpPoint.I]); } else if (_colorMode == ColorMode::IntensityHeightBlend) { float x = (tmpPoint.P.z() - this->_boundingBoxGlobal.zMin()) / (this->_boundingBoxGlobal.zMax() - this->_boundingBoxGlobal.zMin()); //// sigmoid //x = (x - 0.5) * 4; //x = 1. / (1. + exp(-5 * x)); int index = x * 255; index = std::max(0, std::min(255, index)); osg::Vec4 color = _colorBar[index]; color = color * 0.9 + color * (tmpPoint.I / 255.) * 0.1; color.w() = 1.0; colorArray->push_back(color); } } QDataStream dataStream(&buffer, QIODevice::WriteOnly); dataStream.setByteOrder(QDataStream::LittleEndian); // TODO: batch table // feature table QByteArray featureTableJsonByte; { QJsonObject featureTableJson; featureTableJson["POINTS_LENGTH"] = (int)pointArray->size(); if (center != osg::Vec3d(0, 0, 0)) { featureTableJson["RTC_CENTER"] = QJsonArray({ center.x(), center.y(), center.z() }); } QJsonObject posjson; posjson.insert("byteOffset", 0); featureTableJson["POSITION"] = posjson; if (colorArray->size() != 0) { QJsonObject colorjson; long long length = pointArray->size() * 3 * 4; colorjson.insert("byteOffset", length); featureTableJson["RGB"] = colorjson; } featureTableJsonByte = QJsonDocument(featureTableJson) .toJson(QJsonDocument::Compact); while (featureTableJsonByte.size() % 8 != 0) { featureTableJsonByte.append(' '); } } QByteArray featureTableBinary; { std::ostringstream oss; for (int i = 0; i < pointArray->size(); ++i) { osg::Vec3f pt = pointArray->at(i); oss.write((char*)&(pt[0]), sizeof(float)); oss.write((char*)&(pt[1]), sizeof(float)); oss.write((char*)&(pt[2]), sizeof(float)); } for (int i = 0; i < colorArray->size(); ++i) { osg::Vec4f color = colorArray->at(i); uint8_t r = color[0] * 255; uint8_t g = color[1] * 255; uint8_t b = color[2] * 255; oss.write((char*)&(r), sizeof(uint8_t)); oss.write((char*)&(g), sizeof(uint8_t)); oss.write((char*)&(b), sizeof(uint8_t)); //oss.write((char*)&(colorn[3]), sizeof(int)); } oss.flush(); std::string str = oss.str(); featureTableBinary.append(str.data(), str.size()); while (featureTableJsonByte.size() % 8 != 0) { featureTableJsonByte.append(' '); } } uint32_t version = 1; uint32_t featureTableJSONByteLength = featureTableJsonByte.size(); uint32_t featureTableBinaryByteLength = featureTableBinary.size(); uint32_t batchTableJSONByteLength = 0; uint32_t batchTableBinaryByteLength = 0; uint32_t headerLength = 28 + featureTableJSONByteLength + featureTableBinaryByteLength + batchTableJSONByteLength + batchTableBinaryByteLength; //header { dataStream.writeRawData("pnts", 4); dataStream << version; dataStream << headerLength; dataStream << featureTableJSONByteLength; dataStream << featureTableBinaryByteLength; dataStream << batchTableJSONByteLength; dataStream << batchTableBinaryByteLength; } dataStream.writeRawData(featureTableJsonByte.data(), featureTableJsonByte.size()); dataStream.writeRawData(featureTableBinary.data(), featureTableBinary.size()); return true; }

2.5写json

按照3dtiles数据格式要求，需要有一个总的tileset.json文件，然后每个tile中需要有根节点的json文件。