Modified odm_georef to output laz/las

Former-commit-id: ebd6cf83a2
pull/1161/head
Piero Toffanin 2018-12-26 12:51:03 -05:00
rodzic de80696426
commit 9346eef114
5 zmienionych plików z 101 dodań i 67 usunięć

Wyświetl plik

@ -26,6 +26,14 @@ set(OpenCV_LIBS opencv_core opencv_imgproc opencv_highgui)
include_directories(${PCL_ROOT}/include/pcl-${PCL_VERSION_MAJOR}.${PCL_VERSION_MINOR})
include_directories(${EIGEN_ROOT})
# PDAL and jsoncpp
find_package(PDAL REQUIRED CONFIG)
message(${PDAL_INCLUDE_DIRS})
include_directories(${PDAL_INCLUDE_DIRS})
include_directories("${PROJECT_SOURCE_DIR}/../../SuperBuild/src/pdal/vendor/jsoncpp/dist")
link_directories(${PDAL_LIBRARY_DIRS})
add_definitions(${PDAL_DEFINITIONS})
# Add source directory
aux_source_directory("./src" SRC_LIST)
@ -33,4 +41,4 @@ aux_source_directory("./src" SRC_LIST)
add_executable(${PROJECT_NAME} ${SRC_LIST})
# Link
target_link_libraries(${PROJECT_NAME} ${PCL_COMMON_LIBRARIES} ${PCL_IO_LIBRARIES} ${PCL_SURFACE_LIBRARIES} ${PROJ4_LIBRARY} ${OpenCV_LIBS})
target_link_libraries(${PROJECT_NAME} ${PCL_COMMON_LIBRARIES} ${PCL_IO_LIBRARIES} ${PCL_SURFACE_LIBRARIES} ${PROJ4_LIBRARY} ${OpenCV_LIBS} jsoncpp ${PDAL_LIBRARIES})

Wyświetl plik

@ -13,8 +13,6 @@ using namespace std;
// This
#include "Georef.hpp"
#define IS_BIG_ENDIAN (*(uint16_t *)"\0\xff" < 0x100)
std::ostream& operator<<(std::ostream &os, const GeorefSystem &geo)
{
return os << setiosflags(ios::fixed) << setprecision(7) << geo.system_ << "\n" << geo.eastingOffset_ << " " << geo.northingOffset_;
@ -413,6 +411,16 @@ void Georef::parseArguments(int argc, char *argv[])
log_ << "Reading point cloud from: " << inputPointCloudFilename_ << "\n";
georeferencePointCloud_ = true;
}
else if(argument == "-outputPointCloudSrs" && argIndex < argc)
{
argIndex++;
if (argIndex >= argc)
{
throw GeorefException("Argument '" + argument + "' expects 1 more input following it, but no more inputs were provided.");
}
outputPointCloudSrs_ = std::string(argv[argIndex]);
log_ << "Using SRS for point cloud: " << outputPointCloudSrs_ << "\n";
}
else if(argument == "-gcpFile" && argIndex < argc)
{
argIndex++;
@ -590,10 +598,13 @@ void Georef::printHelp()
// log_ << "The resized resolution used in bundler.\n\n";
log_ << "\"-outputFile <path>\" (optional, default <inputFile>_geo)" << "\n";
log_ << "\"Output obj file that will contain the georeferenced texture mesh.\n\n";
log_ << "Output obj file that will contain the georeferenced texture mesh.\n\n";
log_ << "\"-outputPointCloudFile <path>\" (mandatory if georeferencing a point cloud)" << "\n";
log_ << "\"Output ply file that will contain the georeferenced point cloud.\n\n";
log_ << "Output las/laz file that will contain the georeferenced point cloud.\n\n";
log_ << "\"-outputPointCloudSrs <srs>\" (recommended if georeferencing a point cloud)" << "\n";
log_ << "Output spatial reference system to use for output point cloud.\n\n";
log_.setIsPrintingInCout(printInCoutPop);
}
@ -1406,72 +1417,43 @@ void Georef::performFinalTransform(Mat4 &transMat, pcl::TextureMesh &mesh, pcl::
template <typename Scalar>
void Georef::transformPointCloud(const char *inputFile, const Eigen::Transform<Scalar, 3, Eigen::Affine> &transform, const char *outputFile){
try{
pcl::PointCloud<pcl::PointXYZRGBNormal> pointCloud;
if(pcl::io::loadPLYFile<pcl::PointXYZRGBNormal> (inputFile, pointCloud) == -1) {
throw GeorefException("Error when reading point cloud:\n" + std::string(inputFile) + "\n");
}
else
{
log_ << "Successfully loaded " << pointCloud.size() << " points with corresponding normals from file.\n";
}
log_ << "Writing transformed point cloud to " << outputFile << "...\n";
log_ << "Transforming point cloud...\n";
// We don't use PCL's built-in functions
// because PCL does not support double coordinates
// precision
std::ofstream f (outputFile);
f << "ply" << std::endl;
// PDAL pipeline: ply reader --> matrix transform --> las writer.
if (IS_BIG_ENDIAN){
f << "format binary_big_endian 1.0" << std::endl;
pdal::Options inPlyOpts;
inPlyOpts.add("filename", inputFile);
pdal::PointTable table;
pdal::PlyReader plyReader;
plyReader.setOptions(inPlyOpts);
pdal::MatrixTransformFilter<double> transformFilter(transform);
transformFilter.setInput(plyReader);
pdal::Options outLasOpts;
outLasOpts.add("filename", outputFile);
std::string s(outputFile);
if (s.substr(s.find_last_of(".") + 1) == "laz"){
outLasOpts.add("compression", "laszip");
}
outLasOpts.add("offset_x", georefSystem_.eastingOffset_);
outLasOpts.add("offset_y", georefSystem_.northingOffset_);
outLasOpts.add("offset_z", 0);
if (!outputPointCloudSrs_.empty()){
outLasOpts.add("a_srs", outputPointCloudSrs_);
}else{
f << "format binary_little_endian 1.0" << std::endl;
log_ << "WARNING! A point cloud will be generated without a SRS, because -outputPointCloudSrs was not specified!\n";
}
const char *type = "double";
if (sizeof(Scalar) == sizeof(float)){
type = "float";
}
f << "element vertex " << pointCloud.size() << std::endl
<< "property " << type << " x" << std::endl
<< "property " << type << " y" << std::endl
<< "property " << type << " z" << std::endl
<< "property uchar red" << std::endl
<< "property uchar green" << std::endl
<< "property uchar blue" << std::endl
<< "end_header" << std::endl;
struct PlyPoint{
Scalar x;
Scalar y;
Scalar z;
uint8_t r;
uint8_t g;
uint8_t b;
} p;
size_t psize = sizeof(Scalar) * 3 + sizeof(uint8_t) * 3;
// Transform
for (unsigned int i = 0; i < pointCloud.size(); i++){
Scalar x = static_cast<Scalar>(pointCloud[i].x);
Scalar y = static_cast<Scalar>(pointCloud[i].y);
Scalar z = static_cast<Scalar>(pointCloud[i].z);
p.r = pointCloud[i].r;
p.g = pointCloud[i].g;
p.b = pointCloud[i].b;
p.x = static_cast<Scalar> (transform (0, 0) * x + transform (0, 1) * y + transform (0, 2) * z + transform (0, 3));
p.y = static_cast<Scalar> (transform (1, 0) * x + transform (1, 1) * y + transform (1, 2) * z + transform (1, 3));
p.z = static_cast<Scalar> (transform (2, 0) * x + transform (2, 1) * y + transform (2, 2) * z + transform (2, 3));
f.write(reinterpret_cast<char*>(&p), psize);
// TODO: normals can be computed using the inverse transpose
// https://paroj.github.io/gltut/Illumination/Tut09%20Normal%20Transformation.html
}
f.close();
pdal::LasWriter lasWriter;
lasWriter.setOptions(outLasOpts);
lasWriter.setInput(transformFilter);
lasWriter.prepare(table);
lasWriter.execute(table);
log_ << "Point cloud file saved.\n";
}

Wyświetl plik

@ -18,6 +18,9 @@
// Transformation
#include "FindTransform.hpp"
// PDAL matrix transform filter
#include "MatrixTransformFilter.hpp"
/*!
* \brief The GeorefSystem struct is used to store information about a georeference system.
*/
@ -275,7 +278,8 @@ private:
std::string outputObjFilename_; /**< The path to the output mesh obj file. **/
std::string inputPointCloudFilename_; /**< The path to the input point cloud file. **/
std::string outputPointCloudFilename_; /**< The path to the output point cloud file. **/
std::string georefFilename_; /**< The path to the output offset file. **/
std::string georefFilename_; /**< The path to the output offset file. **/
std::string outputPointCloudSrs_; /**< The spatial reference system of the point cloud file to be written. Can be an EPSG string (e.g. “EPSG:26910”) or a WKT string. **/
bool georeferencePointCloud_;
bool exportCoordinateFile_;

Wyświetl plik

@ -0,0 +1,2 @@
#include "MatrixTransformFilter.hpp"

Wyświetl plik

@ -0,0 +1,38 @@
#pragma once
#include <pcl/common/eigen.h>
#include <pdal/PointTable.hpp>
#include <pdal/PointView.hpp>
#include <pdal/io/PlyReader.hpp>
#include <pdal/io/LasWriter.hpp>
#include <pdal/Options.hpp>
#include <pdal/Filter.hpp>
namespace pdal{
template <typename Scalar>
class MatrixTransformFilter : public Filter{
Eigen::Transform<Scalar, 3, Eigen::Affine> transform;
public:
MatrixTransformFilter(const Eigen::Transform<Scalar, 3, Eigen::Affine> &transform)
: transform(transform){};
std::string getName() const { return "MatrixTransformFilter"; }
virtual void filter(PointView &view)
{
for (PointId id = 0; id < view.size(); ++id)
{
Scalar x = view.getFieldAs<Scalar>(Dimension::Id::X, id);
Scalar y = view.getFieldAs<Scalar>(Dimension::Id::Y, id);
Scalar z = view.getFieldAs<Scalar>(Dimension::Id::Z, id);
view.setField(pdal::Dimension::Id::X, id, transform (0, 0) * x + transform (0, 1) * y + transform (0, 2) * z + transform (0, 3));
view.setField(pdal::Dimension::Id::Y, id, transform (1, 0) * x + transform (1, 1) * y + transform (1, 2) * z + transform (1, 3));
view.setField(pdal::Dimension::Id::Z, id, transform (2, 0) * x + transform (2, 1) * y + transform (2, 2) * z + transform (2, 3));
}
}
// TODO: implement streaming mode
};
}