Add chamfer distance to PDALUtils and expose via kernel

doc/apps/chamfer.rst

@@ -0,0 +1,50 @@
+.. _chamfer_command:
+
+********************************************************************************
+chamfer
+********************************************************************************
+
+The ``chamfer`` command is used to compute the Chamfer distance between two
+point clouds. The Chamfer distance is computed by summing the squared distances
+between nearest neighbor correspondences of two point clouds.
+
+More formally, for two non-empty subsets :math:`X` and :math:`Y`, the Chamfer
+distance :math:`d_{CD}(X,Y)` is
+
+.. math::
+
+  d_{CD}(X,Y) = \sum_{x \in X} \operatorname*{min}_{y \in Y} ||x-y||^2_2 + \sum_{y \in Y} \operatorname*{min}_{x \in X} ||x-y||^2_2
+  
+::
+
+    $ pdal chamfer <source> <candidate>
+
+::
+
+    --source arg     Non-positional option for specifying filename of source file.
+    --candidate arg  Non-positional option for specifying filename to test against source.
+
+The algorithm makes no distinction between source and candidate files (i.e.,
+they can be transposed with no affect on the computed distance).
+
+The command returns 0 along with a JSON-formatted message summarizing the PDAL
+version, source and candidate filenames, and the Chamfer distance. Identical
+point clouds will return a Chamfer distance of 0.
+
+::
+
+    $ pdal chamfer source.las candidate.las
+    {
+      "filenames":
+      [
+        "\/path\/to\/source.las",
+        "\/path\/to\/candidate.las"
+      ],
+      "chamfer": 1.303648726,
+      "pdal_version": "1.3.0 (git-version: 191301)"
+    }
+
+.. note::
+
+  The Chamfer distance is computed for XYZ coordinates only and as such says
+  nothing about differences in other dimensions or metadata.

kernels/ChamferKernel.cpp

@@ -0,0 +1,95 @@
+/******************************************************************************
+ * Copyright (c) 2020, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following
+ * conditions are met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided
+ *       with the distribution.
+ *     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+ *       names of its contributors may be used to endorse or promote
+ *       products derived from this software without specific prior
+ *       written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+ * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+ * OF SUCH DAMAGE.
+ ****************************************************************************/
+
+#include "ChamferKernel.hpp"
+
+#include <memory>
+
+#include <pdal/PDALUtils.hpp>
+#include <pdal/PointView.hpp>
+#include <pdal/pdal_config.hpp>
+
+namespace pdal
+{
+
+static StaticPluginInfo const s_info{"kernels.chamfer", "Chamfer Kernel",
+                                     "http://pdal.io/apps/chamfer.html"};
+
+CREATE_STATIC_KERNEL(ChamferKernel, s_info)
+
+std::string ChamferKernel::getName() const
+{
+    return s_info.name;
+}
+
+void ChamferKernel::addSwitches(ProgramArgs& args)
+{
+    Arg& source = args.add("source", "Source filename", m_sourceFile);
+    source.setPositional();
+    Arg& candidate =
+        args.add("candidate", "Candidate filename", m_candidateFile);
+    candidate.setPositional();
+}
+
+PointViewPtr ChamferKernel::loadSet(const std::string& filename,
+                                    PointTableRef table)
+{
+    Stage& reader = makeReader(filename, "");
+    reader.prepare(table);
+    PointViewSet viewSet = reader.execute(table);
+    assert(viewSet.size() == 1);
+    return *viewSet.begin();
+}
+
+int ChamferKernel::execute()
+{
+    ColumnPointTable srcTable;
+    PointViewPtr srcView = loadSet(m_sourceFile, srcTable);
+
+    ColumnPointTable candTable;
+    PointViewPtr candView = loadSet(m_candidateFile, candTable);
+
+    double chamfer = Utils::computeChamfer(srcView, candView);
+
+    MetadataNode root;
+    root.add("filenames", m_sourceFile);
+    root.add("filenames", m_candidateFile);
+    root.add("chamfer", chamfer);
+    root.add("pdal_version", Config::fullVersionString());
+    Utils::toJSON(root, std::cout);
+
+    return 0;
+}
+
+} // namespace pdal

kernels/ChamferKernel.hpp

@@ -0,0 +1,60 @@
+/******************************************************************************
+ * Copyright (c) 2020, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following
+ * conditions are met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided
+ *       with the distribution.
+ *     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+ *       names of its contributors may be used to endorse or promote
+ *       products derived from this software without specific prior
+ *       written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+ * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+ * OF SUCH DAMAGE.
+ ****************************************************************************/
+
+#pragma once
+
+#include <pdal/Kernel.hpp>
+#include <pdal/Stage.hpp>
+#include <pdal/util/FileUtils.hpp>
+
+namespace pdal
+{
+
+class PointView;
+
+class PDAL_DLL ChamferKernel : public Kernel
+{
+public:
+    std::string getName() const;
+    int execute(); // overrride
+
+private:
+    virtual void addSwitches(ProgramArgs& args);
+    PointViewPtr loadSet(const std::string& filename, PointTableRef table);
+
+    std::string m_sourceFile;
+    std::string m_candidateFile;
+};
+
+} // namespace pdal

pdal/PDALUtils.cpp

@@ -364,32 +364,27 @@ double computeHausdorff(PointViewPtr srcView, PointViewPtr candView)
 {
     using namespace Dimension;
 
-    KD3Index srcIndex(*srcView);
-    srcIndex.build();
-
-    KD3Index candIndex(*candView);
-    candIndex.build();
+    KD3Index &srcIndex = srcView->build3dIndex();
+    KD3Index &candIndex = candView->build3dIndex();
 
     double maxDistSrcToCand = std::numeric_limits<double>::lowest();
     double maxDistCandToSrc = std::numeric_limits<double>::lowest();
 
-    for (PointId i = 0; i < srcView->size(); ++i)
+    for (PointRef p : *srcView)
     {
         PointIdList indices(1);
         std::vector<double> sqr_dists(1);
-        PointRef srcPoint = srcView->point(i);
-        candIndex.knnSearch(srcPoint, 1, &indices, &sqr_dists);
+        candIndex.knnSearch(p, 1, &indices, &sqr_dists);
 
         if (sqr_dists[0] > maxDistSrcToCand)
             maxDistSrcToCand = sqr_dists[0];
     }
 
-    for (PointId i = 0; i < candView->size(); ++i)
+    for (PointRef q : *candView)
     {
         PointIdList indices(1);
         std::vector<double> sqr_dists(1);
-        PointRef candPoint = candView->point(i);
-        srcIndex.knnSearch(candPoint, 1, &indices, &sqr_dists);
+        srcIndex.knnSearch(q, 1, &indices, &sqr_dists);
 
         if (sqr_dists[0] > maxDistCandToSrc)
             maxDistCandToSrc = sqr_dists[0];
@@ -401,5 +396,33 @@ double computeHausdorff(PointViewPtr srcView, PointViewPtr candView)
     return (std::max)(maxDistSrcToCand, maxDistCandToSrc);
 }
 
+double computeChamfer(PointViewPtr srcView, PointViewPtr candView)
+{
+    using namespace Dimension;
+
+    KD3Index &srcIndex = srcView->build3dIndex();
+    KD3Index &candIndex = candView->build3dIndex();
+
+    double sum1(0.0);
+    for (PointRef p : *srcView)
+    {
+        PointIdList indices(1);
+        std::vector<double> sqr_dists(1);
+        candIndex.knnSearch(p, 1, &indices, &sqr_dists);
+        sum1 += sqr_dists[0];
+    }
+
+    double sum2(0.0);
+    for (PointRef q : *candView)
+    {
+        PointIdList indices(1);
+        std::vector<double> sqr_dists(1);
+        srcIndex.knnSearch(q, 1, &indices, &sqr_dists);
+        sum2 += sqr_dists[0];
+    }
+
+    return sum1 + sum2;
+}
+
 } // namespace Utils
 } // namespace pdal

pdal/PDALUtils.hpp

@@ -273,6 +273,7 @@ bool PDAL_DLL isRemote(const std::string& path);
 bool PDAL_DLL fileExists(const std::string& path);
 std::vector<std::string> PDAL_DLL maybeGlob(const std::string& path);
 double PDAL_DLL computeHausdorff(PointViewPtr srcView, PointViewPtr candView);
+double PDAL_DLL computeChamfer(PointViewPtr srcView, PointViewPtr candView);
 
 } // namespace Utils
 } // namespace pdal

test/unit/CMakeLists.txt

@@ -338,6 +338,7 @@ PDAL_ADD_TEST(pc2pc_test FILES apps/pc2pcTest.cpp)
 if (BUILD_PIPELINE_TESTS)
     PDAL_ADD_TEST(pcpipeline_test_json FILES apps/pcpipelineTestJSON.cpp)
 endif()
+PDAL_ADD_TEST(chamfer_test FILES apps/ChamferTest.cpp)
 PDAL_ADD_TEST(hausdorff_test FILES apps/HausdorffTest.cpp)
 PDAL_ADD_TEST(random_test FILES apps/RandomTest.cpp)
 PDAL_ADD_TEST(translate_test FILES apps/TranslateTest.cpp)

test/unit/apps/ChamferTest.cpp

@@ -0,0 +1,95 @@
+/******************************************************************************
+ * Copyright (c) 2020, Bradley J Chambers (brad.chambers@gmail.com)
+ *
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following
+ * conditions are met:
+ *
+ *     * Redistributions of source code must retain the above copyright
+ *       notice, this list of conditions and the following disclaimer.
+ *     * Redistributions in binary form must reproduce the above copyright
+ *       notice, this list of conditions and the following disclaimer in
+ *       the documentation and/or other materials provided
+ *       with the distribution.
+ *     * Neither the name of Hobu, Inc. or Flaxen Geo Consulting nor the
+ *       names of its contributors may be used to endorse or promote
+ *       products derived from this software without specific prior
+ *       written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
+ * OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
+ * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+ * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
+ * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
+ * OF SUCH DAMAGE.
+ ****************************************************************************/
+
+#include <string>
+
+#include <pdal/PDALUtils.hpp>
+#include <pdal/PointView.hpp>
+#include <pdal/pdal_test_main.hpp>
+
+#include "Support.hpp"
+
+using namespace pdal;
+
+TEST(Chamfer, kernel)
+{
+    std::string A = Support::datapath("autzen/autzen-thin.las");
+    std::string B = Support::datapath("las/autzen_trim.las");
+    std::string output;
+
+    const std::string cmd =
+        Support::binpath(Support::exename("pdal")) + " chamfer " + A + " " + B;
+
+    EXPECT_EQ(Utils::run_shell_command(cmd, output), 0);
+    EXPECT_TRUE(output.find("\"chamfer\": 5.907628766e+10") !=
+                std::string::npos);
+}
+
+TEST(Chamfer, distance)
+{
+    PointTable table;
+    PointLayoutPtr layout(table.layout());
+
+    layout->registerDim(Dimension::Id::X);
+    layout->registerDim(Dimension::Id::Y);
+    layout->registerDim(Dimension::Id::Z);
+
+    PointViewPtr src(new PointView(table));
+    src->setField(Dimension::Id::X, 0, 0.0);
+    src->setField(Dimension::Id::Y, 0, 0.0);
+    src->setField(Dimension::Id::Z, 0, 0.0);
+
+    PointViewPtr cand(new PointView(table));
+    cand->setField(Dimension::Id::X, 0, 1.0);
+    cand->setField(Dimension::Id::Y, 0, 0.0);
+    cand->setField(Dimension::Id::Z, 0, 0.0);
+
+    cand->setField(Dimension::Id::X, 1, 0.0);
+    cand->setField(Dimension::Id::Y, 1, 2.0);
+    cand->setField(Dimension::Id::Z, 1, 0.0);
+
+    EXPECT_EQ(6.0, Utils::computeChamfer(src, cand));
+
+    cand->setField(Dimension::Id::X, 1, 0.0);
+    cand->setField(Dimension::Id::Y, 1, 0.0);
+    cand->setField(Dimension::Id::Z, 1, 3.0);
+
+    EXPECT_EQ(11.0, Utils::computeChamfer(src, cand));
+
+    src->setField(Dimension::Id::X, 0, 1.0);
+    src->setField(Dimension::Id::Y, 0, 1.0);
+    src->setField(Dimension::Id::Z, 0, 1.0);
+
+    EXPECT_EQ(10.0, Utils::computeChamfer(src, cand));
+}