Compare alg - using M3C2 (#83)

* tool to compare point clouds

* only add the tool if pdal version 2.10 or higher

* download two datasets for testing compare tool

* test for compare

* add skip condition to avoid failing tests

* fix

* remove debug leftover

* compare should never support streaming

* add step-sample instead of  sample-pct

* update pipeline

* use build directory

* values changed due to sampling change

* update year

* format

* remove

* update parameters

* can be empty, with default value 0.0

* rename parameter

* if subsamplingCellSize == 0 skip sampling

* do not store writer

* update test

* add another test without subsampling-cell-size set

* update data path

* update readme

Author: JanCaha

CMakeLists.txt

@@ -30,6 +30,7 @@ add_executable(pdal_wrench
     src/utils.cpp
     src/vpc.cpp
     src/height_above_ground.cpp
+    src/compare.cpp
 
     src/tile/tile.cpp
     src/tile/BufferCache.cpp

README.md

@@ -218,6 +218,15 @@ pdal_wrench height_above_ground --input=data.las --output=data_hag.las --algorit
 pdal_wrench height_above_ground --input=data.las --output=data_hag.las --algorithm=delaunay --replace-z=true --delaunay-count=10
 ```
 
+## compare
+
+Compares two point clouds using M3C2 algorithm and outputs a point cloud with new dimensions: m3c2_distance, m3c2_uncertainty, m3c2_significant, m3c2_std_dev1, m3c2_std_dev2, m3c2_count1 and m3c2_count2. The input data is subsampled to create set of core points, subsampling can be modified using subsampling-cell-size parameter, if it is set to 0.0, no subsampling is done and all points are used.
+
+```
+pdal_wrench compare --input=first.las --input-compare=second.las --output=changes.las --subsampling-cell-size=1.0 --normal-radius=2.0 --cyl-radius=2.0 --cyl-halflen=5.0 --reg-error=0.0 --cyl-orientation=up
+```
+
+
 # Virtual Point Clouds (VPC)
 
 This is similar to GDAL's VRT - a single file referring to other files that contain actual data. Software then may handle all data as a single dataset.
@@ -271,3 +280,4 @@ When algorithms create derived VPCs, by default they use uncompressed LAS, but `
 | height_above_ground | multi-threaded | per file |
 | filter_noise | multi-threaded | per file |
 | classify_ground | multi-threaded | per file |
+| compare | no | only available in PDAL Version > 2.10 |

src/alg.hpp

@@ -441,6 +441,34 @@ struct HeightAboveGround : public Alg
 
     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 ComparePointClouds : public Alg
+{   
+    ComparePointClouds() { isStreaming = false; }
+
+    // parameters from the user
+    std::string outputFile;
+
+    std::string comparedInputFile;
+    double subsamplingCellSize;  // cell size for Poisson sampling
+    double normalRadius = 2.0;
+    double cylRadius = 2.0;
+    double cylHalflen = 5.0;
+    double regError = 0.0;
+    std::string cylOrientation = "up";
+
+    // args - initialized in addArgs()
+    pdal::Arg* argOutput = nullptr;
+    pdal::Arg* argOutputFormat = nullptr;
+    pdal::Arg* argComparedInputFile = nullptr;
+    pdal::Arg* argOrientation = nullptr;
+
     // impl
     virtual void addArgs() override;
     virtual bool checkArgs() override;

src/compare.cpp

@@ -0,0 +1,184 @@
+/*****************************************************************************
+ *   Copyright (c) 2026, 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 <filesystem>
+#include <iostream>
+#include <thread>
+
+#include <pdal/PipelineManager.hpp>
+#include <pdal/Polygon.hpp>
+#include <pdal/Stage.hpp>
+#include <pdal/pdal_types.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include <gdal_utils.h>
+
+#include "alg.hpp"
+#include "utils.hpp"
+#include "vpc.hpp"
+
+using namespace pdal;
+
+void ComparePointClouds::addArgs() 
+{
+  argOutput = &programArgs.add("output,o", "Output point cloud file", outputFile);
+  argComparedInputFile = &programArgs.add("input-compare", "Point cloud file to compare against input", comparedInputFile);
+  
+  programArgs.add("subsampling-cell-size", "Minimum spacing between points", subsamplingCellSize, 0.0);
+
+  programArgs.add( "normal-radius", "Radius of the sphere around each core point that defines the neighbors from which normals are calculated.", normalRadius, 2.0);
+  programArgs.add("cyl-radius", "Radius of the cylinder inside of which points are searched for when calculating change", cylRadius, 2.0);
+  programArgs.add("cyl-halflen", "The half-length of the cylinder of neighbors used for calculating change", cylHalflen, 5.0);
+  programArgs.add("reg-error", "Registration error", regError, 0.0);
+  argOrientation = &programArgs.add( "cyl-orientation", "Which direction to orient the cylinder/normal vector used for comparison between the two point clouds. (up, origin, none)", cylOrientation, "up");
+}
+
+bool ComparePointClouds::checkArgs() 
+{
+
+  if (ends_with(inputFile, ".vpc")) {
+    std::cerr << "input cannot be a VPC file" << std::endl;
+    return false;
+  }
+
+  if (ends_with(comparedInputFile, ".vpc")) {
+    std::cerr << "compared input cannot be a VPC file" << std::endl;
+    return false;
+  }
+
+  if (!argOutput->set()) {
+    std::cerr << "missing output" << std::endl;
+    return false;
+  }
+
+  if (!argComparedInputFile->set()) {
+    std::cerr << "missing compared input file" << std::endl;
+    return false;
+  }
+
+  if (argOrientation->set()) {
+    if (cylOrientation != "up" && cylOrientation != "origin" &&
+        cylOrientation != "none") {
+      std::cerr << "unknown orientation: " << cylOrientation << std::endl;
+      return false;
+    }
+  }
+
+  return true;
+}
+
+static std::unique_ptr<PipelineManager> pipeline(ParallelJobInfo *tile, std::string compareFile, double subsamplingCellSize, double normalRadius, double cylRadius, double cylHalflen, double regError, std::string cylOrientation)
+{
+  std::unique_ptr<PipelineManager> manager(new PipelineManager);
+
+  Stage &reader1 = makeReader(manager.get(), tile->inputFilenames[0]);
+  Stage *last = &reader1;
+
+  // filtering
+  if (!tile->filterBounds.empty())
+  {
+      Options filter_opts;
+      filter_opts.add(pdal::Option("bounds", tile->filterBounds));
+
+      if (readerSupportsBounds(reader1))
+      {
+          // Reader of the format can do the filtering - use that whenever possible!
+          reader1.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);
+  }     
+  
+  Stage *reader1FilteredCropped = last;
+
+  Stage &reader2 = makeReader(manager.get(), compareFile);
+  last = &reader2;
+
+  // filtering
+  if (!tile->filterBounds.empty())
+  {
+      Options filter_opts;
+      filter_opts.add(pdal::Option("bounds", tile->filterBounds));
+
+      if (readerSupportsBounds(reader2))
+      {
+          // Reader of the format can do the filtering - use that whenever possible!
+          reader2.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);
+  }
+  Stage *reader2FilteredCropped = last;
+
+  Stage *corePoints = nullptr;
+  if (subsamplingCellSize != 0.0)
+  {
+    pdal::Options sample_opts;
+    sample_opts.add(pdal::Option("cell", subsamplingCellSize));
+    corePoints = &manager->makeFilter("filters.sample", *reader1FilteredCropped, sample_opts);
+  }
+  else
+  {
+    corePoints = reader1FilteredCropped;
+  }
+
+  Options compare_opts;
+  compare_opts.add(pdal::Option("normal_radius", normalRadius));
+  compare_opts.add(pdal::Option("cyl_radius", cylRadius));
+  compare_opts.add(pdal::Option("cyl_halflen", cylHalflen));
+  compare_opts.add(pdal::Option("reg_error", regError));
+  compare_opts.add(pdal::Option("orientation", cylOrientation));
+
+  Stage *filterM3c2 = &manager->makeFilter("filters.m3c2", compare_opts);
+
+  // inputs to M3C2 filter are origin file, compared file and core points (in order)
+  filterM3c2->setInput(*reader1FilteredCropped);
+  filterM3c2->setInput(*reader2FilteredCropped);
+  filterM3c2->setInput(*corePoints);
+  
+  makeWriter(manager.get(), tile->outputFilename, filterM3c2);
+
+  return manager;
+}
+
+void ComparePointClouds::preparePipelines(std::vector<std::unique_ptr<PipelineManager>> &pipelines)
+{
+  ParallelJobInfo tile(ParallelJobInfo::Single, BOX2D(), filterExpression,
+                       filterBounds);
+  tile.inputFilenames.push_back(inputFile);
+  tile.outputFilename = outputFile;
+  pipelines.push_back(pipeline(&tile, comparedInputFile, subsamplingCellSize, normalRadius, cylRadius, cylHalflen, regError, cylOrientation));
+}
+
+void ComparePointClouds::finalize( std::vector<std::unique_ptr<PipelineManager>> &)
+{
+}

src/main.cpp

@@ -21,6 +21,7 @@
 
 #include <pdal/util/FileUtils.hpp>
 #include <pdal/pdal_config.hpp>
+#include <pdal/pdal_features.hpp>
 
 #include "gdal_version.h"
 
@@ -42,6 +43,10 @@ void printUsage()
     std::cout << "   build_vpc              Creates a virtual point cloud" << std::endl;
     std::cout << "   classify_ground        Classify ground points" << std::endl;
     std::cout << "   clip                   Outputs only points that are inside of the clipping polygons" << std::endl;
+// TODO this check can be removed once PDAL 2.10 is widely used    
+#if PDAL_VERSION_MAJOR > 2 || ( PDAL_VERSION_MAJOR == 2 && PDAL_VERSION_MINOR >= 10 )
+    std::cout << "   compare                Compares two point clouds for differences" << std::endl;
+#endif    
     std::cout << "   density                Exports a raster where each cell contains number of points" << std::endl;
     std::cout << "   filter_noise           Classify noise points" << std::endl;
     std::cout << "   height_above_ground    Calculates height above ground for each point" << std::endl;
@@ -154,6 +159,14 @@ int main(int argc, char* argv[])
         ClassifyGround classifyGround;
         runAlg(args, classifyGround);
     }
+// TODO this check can be removed once PDAL 2.10 is is widely used    
+#if PDAL_VERSION_MAJOR > 2 || ( PDAL_VERSION_MAJOR == 2 && PDAL_VERSION_MINOR >= 10 )
+    else if (cmd == "compare")
+    {
+        ComparePointClouds comparison;
+        runAlg(args, comparison);
+    }
+#endif    
     else
     {
         std::cerr << "unknown command: " << cmd << std::endl;

tests/conftest.py

@@ -17,15 +17,25 @@ def _prepare_data():
 
     base_data = utils.test_data_filepath("stadium-utm.laz")
 
-    if not base_data.exists():
-        # PDAL autzen stadium dataset
-        url = "https://media.githubusercontent.com/media/PDAL/data/refs/heads/main/autzen/stadium-utm.laz"
-
-        r = requests.get(url, timeout=10 * 60)
-
-        with open(base_data, "wb") as f:
-            f.write(r.content)
-            # Run the pdal_wrench command
+    # PDAL autzen stadium dataset
+    utils.download_data(
+        "https://media.githubusercontent.com/media/PDAL/data/refs/heads/main/autzen/stadium-utm.laz", base_data
+    )
+
+    autzen_2010_data = utils.test_data_filepath("autzen-bmx-2010.las")
+    autzen_2023_data = utils.test_data_filepath("autzen-bmx-2023.las")
+
+    utils.download_data(
+        "https://raw.githubusercontent.com/PDAL/PDAL/a057739ef0c620d4668752becef25f2b196bf174/test/data/autzen/autzen-bmx-2010.las",
+        autzen_2010_data,
+    )
+    utils.download_data(
+        "https://raw.githubusercontent.com/PDAL/PDAL/a057739ef0c620d4668752becef25f2b196bf174/test/data/autzen/autzen-bmx-2023.las",
+        autzen_2023_data,
+    )
+
+    pdal.Reader(autzen_2010_data.as_posix()).pipeline().execute()
+    pdal.Reader(autzen_2023_data.as_posix()).pipeline().execute()
 
     laz_file = pdal.Reader(base_data.as_posix()).pipeline()
     number_points = laz_file.execute()
@@ -217,3 +227,15 @@ def main_copc_file() -> str:
 def main_copc_file_without_classification() -> str:
     "Return path to the main copc file without classification"
     return utils.test_data_filepath("stadium-utm-not-classified.copc.laz").as_posix()
+
+
+@pytest.fixture
+def autzen_2010_file() -> str:
+    "Return path to the autzen 2010 las file"
+    return utils.test_data_filepath("autzen-bmx-2010.las").as_posix()
+
+
+@pytest.fixture
+def autzen_2023_file() -> str:
+    "Return path to the autzen 2023 las file"
+    return utils.test_data_filepath("autzen-bmx-2023.las").as_posix()