diff --git a/CMakeLists.txt b/CMakeLists.txt index db1f194..505fb00 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,6 +24,7 @@ add_executable(pdal_wrench src/to_raster.cpp src/to_raster_tin.cpp src/to_vector.cpp + src/transform.cpp src/translate.cpp src/utils.cpp src/vpc.cpp diff --git a/src/alg.hpp b/src/alg.hpp index e3e8a8b..6832b5c 100644 --- a/src/alg.hpp +++ b/src/alg.hpp @@ -125,6 +125,27 @@ struct Translate : public Alg virtual void finalize(std::vector>& pipelines) override; }; +struct Transform : public Alg +{ + // parameters from the user + std::string outputFile; + bool invert; + std::string matrix; // whitespace delimited 4x4 transformation matrix + std::string outputFormat; // las / laz / copc + + // args - initialized in addArgs() + pdal::Arg* argOutput = nullptr; + pdal::Arg* argOutputFormat = nullptr; + + std::vector tileOutputFiles; + + // impl + virtual void addArgs() override; + virtual bool checkArgs() override; + virtual void preparePipelines(std::vector>& pipelines) override; + virtual void finalize(std::vector>& pipelines) override; +}; + struct Density : public Alg { // parameters from the user diff --git a/src/main.cpp b/src/main.cpp index f4f2ff4..fb44ef1 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -123,6 +123,11 @@ int main(int argc, char* argv[]) Translate translate; runAlg(args, translate); } + else if (cmd == "transformation") + { + Transform transform; + runAlg(args, transform); + } else if (cmd == "tile") { runTile(args); diff --git a/src/transform.cpp b/src/transform.cpp new file mode 100644 index 0000000..1db2079 --- /dev/null +++ b/src/transform.cpp @@ -0,0 +1,171 @@ +/***************************************************************************** + * Copyright (c) 2023, Lutra Consulting Ltd. and Hobu, Inc. * + * * + * All rights reserved. * + * * + * This program is free software; you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation; either version 3 of the License, or * + * (at your option) any later version. * + * * + ****************************************************************************/ + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include "utils.hpp" +#include "alg.hpp" +#include "vpc.hpp" + +using namespace pdal; + +namespace fs = std::filesystem; + + +void Transform::addArgs() +{ + argOutput = &programArgs.add("output,o", "Output point cloud file", outputFile); + argOutputFormat = &programArgs.add("output-format", "Output format (las/laz/copc)", outputFormat); + programArgs.add("invert", "If set to true, applies the inverse of the provided transformation matrix", invert); + programArgs.add("matrix", "A whitespace-delimited transformation matrix. The matrix is assumed to be presented in row-major order. Only matrices with sixteen elements are allowed.", matrix); +} + +bool Transform::checkArgs() +{ + if (!argOutput->set()) + { + std::cerr << "missing output" << std::endl; + return false; + } + + // TODO: or use the same format as the reader? + if (argOutputFormat->set()) + { + if (outputFormat != "las" && outputFormat != "laz") + { + std::cerr << "unknown output format: " << outputFormat << std::endl; + return false; + } + } + else + outputFormat = "las"; // uncompressed by default + + return true; +} + + +static std::unique_ptr pipeline(ParallelJobInfo *tile, bool invert, std::string matrix) +{ + std::unique_ptr manager( new PipelineManager ); + + Stage& r = manager->makeReader( tile->inputFilenames[0], ""); + + Stage *last = &r; + + // filtering + if (!tile->filterBounds.empty()) + { + Options filter_opts; + filter_opts.add(pdal::Option("bounds", tile->filterBounds)); + + if (readerSupportsBounds(r)) + { + // Reader of the format can do the filtering - use that whenever possible! + r.addOptions(filter_opts); + } + else + { + // Reader can't do the filtering - do it with a filter + last = &manager->makeFilter( "filters.crop", *last, filter_opts); + } + } + if (!tile->filterExpression.empty()) + { + Options filter_opts; + filter_opts.add(pdal::Option("expression", tile->filterExpression)); + last = &manager->makeFilter( "filters.expression", *last, filter_opts); + } + + Options transformation_opts; + transformation_opts.add(pdal::Option("invert", invert)); + transformation_opts.add(pdal::Option("matrix", matrix)); + last = &manager->makeFilter( "filters.transformation", *last, transformation_opts ); + + pdal::Options writer_opts; + // TODO: maybe adjust scale as well - depending on the CRS + writer_opts.add(pdal::Option("forward", "header,scale,vlr")); + writer_opts.add(pdal::Option("offset_x", "auto")); + writer_opts.add(pdal::Option("offset_y", "auto")); + writer_opts.add(pdal::Option("offset_z", "auto")); + + (void)manager->makeWriter( tile->outputFilename, "", *last, writer_opts); + + return manager; +} + + +void Transform::preparePipelines(std::vector>& pipelines) +{ + if (ends_with(inputFile, ".vpc")) + { + if (!ends_with(outputFile, ".vpc")) + { + std::cerr << "If input file is a VPC, output should be VPC too." << std::endl; + return; + } + + // for /tmp/hello.vpc we will use /tmp/hello dir for all results + fs::path outputParentDir = fs::path(outputFile).parent_path(); + fs::path outputSubdir = outputParentDir / fs::path(outputFile).stem(); + fs::create_directories(outputSubdir); + + // VPC handling + VirtualPointCloud vpc; + if (!vpc.read(inputFile)) + return; + + for (const VirtualPointCloud::File& f : vpc.files) + { + ParallelJobInfo tile(ParallelJobInfo::FileBased, BOX2D(), filterExpression, filterBounds); + tile.inputFilenames.push_back(f.filename); + + // for input file /x/y/z.las that goes to /tmp/hello.vpc, + // individual output file will be called /tmp/hello/z.las + fs::path inputBasename = fs::path(f.filename).stem(); + tile.outputFilename = (outputSubdir / inputBasename).string() + "." + outputFormat; + + tileOutputFiles.push_back(tile.outputFilename); + + pipelines.push_back(pipeline(&tile, invert, matrix)); + } + } + else + { + ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression, filterBounds); + tile.inputFilenames.push_back(inputFile); + tile.outputFilename = outputFile; + pipelines.push_back(pipeline(&tile, invert, matrix)); + } +} + +void Transform::finalize(std::vector>&) +{ + if (tileOutputFiles.empty()) + return; + + // now build a new output VPC + std::vector args; + args.push_back("--output=" + outputFile); + for (std::string f : tileOutputFiles) + args.push_back(f); + buildVpc(args); +}