Add Transformation algorithm

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

src/alg.hpp

@@ -125,6 +125,27 @@ struct Translate : public Alg
     virtual void finalize(std::vector<std::unique_ptr<PipelineManager>>& 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<std::string> tileOutputFiles;
+
+    // impl
+    virtual void addArgs() override;
+    virtual bool checkArgs() override;
+    virtual void preparePipelines(std::vector<std::unique_ptr<PipelineManager>>& pipelines) override;
+    virtual void finalize(std::vector<std::unique_ptr<PipelineManager>>& pipelines) override;
+};
+
 struct Density : public Alg
 {
     // parameters from the user

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);

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 <iostream>
+#include <filesystem>
+#include <thread>
+
+#include <pdal/PipelineManager.hpp>
+#include <pdal/Stage.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+#include <pdal/pdal_types.hpp>
+#include <pdal/Polygon.hpp>
+
+#include <gdal_utils.h>
+
+#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<PipelineManager> pipeline(ParallelJobInfo *tile, bool invert, std::string matrix)
+{
+    std::unique_ptr<PipelineManager> 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<std::unique_ptr<PipelineManager>>& 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<std::unique_ptr<PipelineManager>>&)
+{
+    if (tileOutputFiles.empty())
+        return;
+
+    // now build a new output VPC
+    std::vector<std::string> args;
+    args.push_back("--output=" + outputFile);
+    for (std::string f : tileOutputFiles)
+        args.push_back(f);
+    buildVpc(args);
+}