Add avg option to NNdistance filter. (#2076)

* Add nn distance dimension.

* Add general nearest neighbor distance filter.

* Migrate KDistanceFilter to more generic NNDistanceFilter.
Close #2071

* Full test.

* Add some test comments.

* Use distances, not square distances when computing averages.

* Remove unnecessary braces.

filters/NNDistanceFilter.cpp

@@ -0,0 +1,142 @@
+/******************************************************************************
+* Copyright (c) 2016, 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 "NNDistanceFilter.hpp"
+
+#include <string>
+#include <vector>
+
+#include <pdal/KDIndex.hpp>
+
+namespace pdal
+{
+
+static PluginInfo const s_info
+{
+    "filters.nndistance",
+    "NN-Distance Filter",
+    "http://pdal.io/stages/filters.nndistance.html"
+};
+
+CREATE_STATIC_STAGE(NNDistanceFilter, s_info)
+
+std::string NNDistanceFilter::getName() const
+{
+    return s_info.name;
+}
+
+
+NNDistanceFilter::NNDistanceFilter() : Filter()
+{}
+
+
+std::istream& operator>>(std::istream& in, NNDistanceFilter::Mode& mode)
+{
+    std::string s;
+    in >> s;
+
+    s = Utils::tolower(s);
+    if (s == "kth")
+        mode = NNDistanceFilter::Mode::Kth;
+    else if (s == "avg")
+        mode = NNDistanceFilter::Mode::Average;
+    else
+        in.setstate(std::ios_base::failbit);
+    return in;
+}
+
+
+std::ostream& operator<<(std::ostream& out, const NNDistanceFilter::Mode& mode)
+{
+    switch (mode)
+    {
+    case NNDistanceFilter::Mode::Kth:
+        out << "kth";
+    case NNDistanceFilter::Mode::Average:
+        out << "avg";
+    }
+    return out;
+}
+
+
+void NNDistanceFilter::addArgs(ProgramArgs& args)
+{
+    args.add("mode", "Distance computation mode (kth, avg)", m_mode, Mode::Kth);
+    args.add("k", "k neighbors", m_k, size_t(10));
+}
+
+
+void NNDistanceFilter::addDimensions(PointLayoutPtr layout)
+{
+    layout->registerDim(Dimension::Id::NNDistance);
+}
+
+
+void NNDistanceFilter::filter(PointView& view)
+{
+    using namespace Dimension;
+
+    // Build the 3D KD-tree.
+    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
+    KD3Index& index = view.build3dIndex();
+
+    // Increment the minimum number of points, as knnSearch will be returning
+    // the query point along with the neighbors.
+    size_t k = m_k + 1;
+
+    // Compute the k-distance for each point. The k-distance is the Euclidean
+    // distance to k-th nearest neighbor.
+    log()->get(LogLevel::Debug) << "Computing k-distances...\n";
+    for (PointId idx = 0; idx < view.size(); ++idx)
+    {
+        std::vector<PointId> indices(k);
+        std::vector<double> sqr_dists(k);
+        index.knnSearch(idx, k, &indices, &sqr_dists);
+        double val;
+        if (m_mode == Mode::Kth)
+            val = std::sqrt(sqr_dists[k - 1]);
+        else // m_mode == Mode::Average
+        {
+            val = 0;
+
+            // We start at 1 since index 0 is the test point.
+            for (size_t i = 1; i < k; ++i)
+                val += std::sqrt(sqr_dists[i]);
+            val /= (k - 1);
+        }
+        view.setField(Dimension::Id::NNDistance, idx, val);
+    }
+}
+
+} // namespace pdal

filters/NNDistanceFilter.hpp

@@ -0,0 +1,76 @@
+/******************************************************************************
+* Copyright (c) 2016, 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/Filter.hpp>
+
+namespace pdal
+{
+
+class PointLayout;
+class PointView;
+class ProgramArgs;
+
+class PDAL_DLL NNDistanceFilter : public Filter
+{
+public:
+    NNDistanceFilter();
+
+    NNDistanceFilter& operator=(const NNDistanceFilter&) = delete;
+    NNDistanceFilter(const NNDistanceFilter&) = delete;
+
+    std::string getName() const;
+
+private:
+    enum class Mode
+    {
+        Kth,
+        Average
+    };
+
+    size_t m_k;
+    Mode m_mode;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void addDimensions(PointLayoutPtr layout);
+    virtual void filter(PointView& view);
+
+    friend std::istream& operator>>(std::istream& in,
+        NNDistanceFilter::Mode& mode);
+    friend std::ostream& operator<<(std::ostream& in,
+        const NNDistanceFilter::Mode& mode);
+};
+
+} // namespace pdal

pdal/Dimension.json

@@ -370,6 +370,10 @@
     "name": "Omit",
     "type": "uint8_t",
     "description": "Used to shallowly mark a point as being omitted without removing it"
+    },
+    {
+    "name": "NNDistance",
+    "type": "double",
+    "description": "Distance metric related to a point's nearest neighbors."
     }
-
 ] }

test/unit/CMakeLists.txt

@@ -137,6 +137,7 @@ PDAL_ADD_TEST(pdal_io_text_writer_test FILES io/TextWriterTest.cpp)
 #
 PDAL_ADD_TEST(pdal_filters_assign_test FILES filters/AssignFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_chipper_test FILES filters/ChipperTest.cpp)
+PDAL_ADD_TEST(pdal_filters_nndistance_test FILES filters/NNDistanceTest.cpp)
 target_include_directories(pdal_filters_chipper_test PRIVATE ${PDAL_JSONCPP_INCLUDE_DIR})
 PDAL_ADD_TEST(pdal_filters_colorinterp_test
     FILES

test/unit/filters/NNDistanceTest.cpp

@@ -0,0 +1,124 @@
+/******************************************************************************
+* Copyright (c) 2018, Hobu Inc. (info@hobu.co)
+*
+* 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 <filters/NNDistanceFilter.hpp>
+#include <pdal/pdal_test_main.hpp>
+#include <pdal/StageFactory.hpp>
+
+#include "Support.hpp"
+
+using namespace pdal;
+
+TEST(NNDistanceTest, kdist)
+{
+    StageFactory f;
+    Stage *reader(f.createStage("readers.faux"));
+    Stage *filter(f.createStage("filters.nndistance"));
+
+    // Make a 10x10x10 grid of points.
+    Options rOpts;
+    rOpts.add("mode", "grid");
+    rOpts.add("bounds", "([0, 10],[0,10],[0,10])");
+    rOpts.add("count", 1000);
+    reader->setOptions(rOpts);
+
+    // Kth distance with k == 3
+    Options fOpts;
+    fOpts.add("mode", "kth");
+    fOpts.add("k", 3);
+    filter->setOptions(fOpts);
+    filter->setInput(*reader);
+
+    PointTable t;
+    filter->prepare(t);
+    PointViewSet s = filter->execute(t);
+    PointViewPtr v = *(s.begin());
+
+    for (PointId i = 0; i < v->size(); ++i)
+        EXPECT_EQ(v->getFieldAs<double>(Dimension::Id::NNDistance, i), 1);
+
+    // Kth distance with k == 4
+    Options opts2;
+    opts2.add("k", 4);
+    filter->setOptions(opts2);
+
+    PointTable t2;
+    filter->prepare(t2);
+    s = filter->execute(t2);
+    v = *(s.begin());
+
+    for (PointId i = 0; i < v->size(); ++i)
+    {
+        double d = v->getFieldAs<double>(Dimension::Id::NNDistance, i);
+        double x = v->getFieldAs<double>(Dimension::Id::X, i);
+        double y = v->getFieldAs<double>(Dimension::Id::Y, i);
+        double z = v->getFieldAs<double>(Dimension::Id::Z, i);
+        if ((x == 9 || x == 0) && (y == 9 || y == 0) && (z == 9 || z == 0))
+            EXPECT_EQ(d, sqrt(2));
+    }
+
+    // Test for avg distance.
+    Options opts3;
+    opts3.add("mode", "avg");
+    opts3.add("k", 6);
+    filter->setOptions(opts3);
+
+    PointTable t3;
+    filter->prepare(t3);
+    s = filter->execute(t3);
+    v = *(s.begin());
+
+    for (PointId i = 0; i < v->size(); ++i)
+    {
+        double d = v->getFieldAs<double>(Dimension::Id::NNDistance, i);
+        double x = v->getFieldAs<double>(Dimension::Id::X, i);
+        double y = v->getFieldAs<double>(Dimension::Id::Y, i);
+        double z = v->getFieldAs<double>(Dimension::Id::Z, i);
+        int xe = (x == 9 || x == 0) ? 1 : 0;
+        int ye = (y == 9 || y == 0) ? 1 : 0;
+        int ze = (z == 9 || z == 0) ? 1 : 0;
+        bool corner = (xe + ye + ze == 3);
+        bool edge = (xe + ye + ze == 2);
+        bool face = (xe + ye + ze == 1);
+        if (corner)
+            EXPECT_EQ(d, (1 + std::sqrt(2)) / 2.0);
+        else if (edge)
+            EXPECT_EQ(d, (2 + std::sqrt(2)) / 3.0);
+        else if (face)
+            EXPECT_EQ(d, (5 + std::sqrt(2)) / 6.0);
+        else
+            EXPECT_EQ(d, 1);
+    }
+}
+