Contribute NeighborClassifierFilter (#1803)

* Contribute KNNAssignFilter.

* hook up KNNAssign filter to StageFactory

* rename processOne --> doOne and make non-virtual

* remove dimension, hard-code Classification as what we update

* change name knnassignfilter --> neighborclassifierfilter

* finish rename KNNAssignFilter --> NeighborClassifierFilter

doc/stages/filters.neighborclassifier.rst

@@ -0,0 +1,78 @@
+.. _filters.neighborclassifier:
+
+filters.neighborclassifier
+===================
+
+The neighborclassifier filter allows you update the value of the classification
+for specific points to a value determined by a K-nearest neighbors vote. 
+For each point, the k nearest neighbors are queried and if more than half of
+them have the same value, the filter updates the selected point accordingly
+
+For example, if an automated classification procedure put/left erroneous 
+vegetation points near the edges of buildings which were largely classified 
+correctly, you could try using this filter to fix that problem.
+
+Similiarly, some automated classification processes result in prediction for
+only a subset of the original point cloud.  This filter could be used to
+extrapolate those predictions to the original.
+
+.. embed::
+
+Example 1
+---------
+
+This pipeline updates the Classification of all points with classification
+1 (unclassified) based on the consensus (majority) of its nearest 10 neighbors.
+
+.. code-block:: json
+    {
+      "pipeline":[
+        "autzen_class.las",
+        {
+          "type" : "filters.neighborclassifier",
+          "domain" : "Classification[1:1]",
+          "k" : 10
+        },
+        {
+          "filename":"autzen_class_refined.las"
+        }
+      ]
+    } 
+
+Example 2
+---------
+
+This pipeline moves all the classifications from "pred.txt"
+to src.las.  Any points in src.las that are not in pred.txt will be
+assigned based on the closest point in pred.txt.
+
+.. code-block:: json
+    {
+      "pipeline":[
+        "src.las",
+        {
+          "type" : "filters.neighborclassifier",
+          "k" : 1,
+          "candidate" : "pred.txt"
+        },
+        {
+          "filename":"dest.las"
+        }
+      ]
+    } 
+
+Options
+-------
+
+candidate
+  A filename which points to the point cloud containing the points which
+  will do the voting.  If not specified, defaults to the input of the filter.
+
+domain
+  A :ref:`range <ranges>` which selects points to be processed by the filter.
+  Can be specified multiple times.  Points satisfying any range will be
+  processed
+
+k
+  An integer which specifies the number of neighbors which vote on each
+  selected point.  

filters/NeighborClassifierFilter.cpp

@@ -0,0 +1,199 @@
+/******************************************************************************
+* Copyright (c) 2017, 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 "NeighborClassifierFilter.hpp"
+
+#include <pdal/PipelineManager.hpp>
+#include <pdal/StageFactory.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include "private/DimRange.hpp"
+
+#include <iostream>
+#include <utility>
+namespace pdal
+{
+
+static PluginInfo const s_info = PluginInfo(
+    "filters.neighborclassifier",
+    "Re-assign some point attributes based KNN voting",
+    "http://pdal.io/stages/filters.neighborclassifier.html" );
+
+CREATE_STATIC_PLUGIN(1, 0, NeighborClassifierFilter, Filter, s_info)
+
+NeighborClassifierFilter::NeighborClassifierFilter() : m_dim(Dimension::Id::Classification)
+{}
+
+
+NeighborClassifierFilter::~NeighborClassifierFilter()
+{}
+
+
+void NeighborClassifierFilter::addArgs(ProgramArgs& args)
+{
+    args.add("domain", "Selects which points will be subject to KNN-based assignmenassignment",
+        m_domainSpec);
+    args.add("k", "Number of nearest neighbors to consult",
+        m_k).setPositional();
+    //args.add("dimension", "Dimension on to be updated", m_dimName).setPositional();
+    Arg& candidate = args.add("candidate", "candidate file name",
+        m_candidateFile);
+}
+
+void NeighborClassifierFilter::initialize()
+{
+    for (auto const& r : m_domainSpec)
+    {
+        try
+        {
+            DimRange range;
+            range.parse(r);
+            m_domain.push_back(range);
+        }
+        catch (const DimRange::error& err)
+        {
+            throwError("Invalid 'domain' option: '" + r + "': " + err.what());
+        }
+    }
+    if (m_k < 1)
+        throwError("Invalid 'k' option: " + std::to_string(m_k) +  ", must be > 0");
+}
+void NeighborClassifierFilter::prepared(PointTableRef table)
+{
+    PointLayoutPtr layout(table.layout());
+
+    for (auto& r : m_domain)
+    {
+        r.m_id = layout->findDim(r.m_name);
+        if (r.m_id == Dimension::Id::Unknown)
+            throwError("Invalid dimension name in 'domain' option: '" +
+                r.m_name + "'.");
+    }
+    std::sort(m_domain.begin(), m_domain.end());
+    //m_dim = layout->findDim(m_dimName);
+
+    //if (m_dim == Dimension::Id::Unknown)
+    //    throwError("Dimension '" + m_dimName + "' not found.");
+}
+
+void NeighborClassifierFilter::doOneNoDomain(PointRef &point, PointRef &temp, KD3Index &kdi)
+{
+    std::vector<PointId> iSrc = kdi.neighbors(point, m_k);
+    double thresh = iSrc.size()/2.0;
+    //std::cout << "iSrc.size() " << iSrc.size() << " thresh " << thresh << std::endl;
+
+    // vote NNs
+    std::map<double, unsigned int> counts;
+    for (PointId id : iSrc)
+    {
+        temp.setPointId(id);
+        double votefor = temp.getFieldAs<double>(m_dim);
+        counts[votefor]++;
+    }
+
+    // pick winner of the vote
+    auto pr = *std::max_element(counts.begin(), counts.end(),
+        [](const std::pair<int, int>& p1, const std::pair<int, int>& p2) {
+        return p1.second < p2.second; });
+
+    // update point
+    auto oldclass = point.getFieldAs<double>(m_dim);
+    auto newclass = pr.first;
+    //std::cout << oldclass << " --> " << newclass << " count " << pr.second << std::endl;
+    if (pr.second > thresh && oldclass != newclass)
+    {    
+        point.setField(m_dim, newclass); 
+    }
+}
+
+bool NeighborClassifierFilter::doOne(PointRef& point, PointRef &temp, KD3Index &kdi)
+{   // update point.  kdi and temp both reference the NN point cloud
+
+    if (m_domain.empty())  // No domain, process all points
+        doOneNoDomain(point, temp, kdi);
+
+    for (DimRange& r : m_domain)
+    {   // process only points that satisfy a domain condition
+        if (r.valuePasses(point.getFieldAs<double>(r.m_id)))
+        {
+            doOneNoDomain(point, temp, kdi);
+            break;
+        }
+    }
+    return true;
+}
+
+PointViewPtr NeighborClassifierFilter::loadSet(const std::string& filename,
+    PointTable& table)
+{
+    PipelineManager mgr;
+
+    Stage& reader = mgr.makeReader(filename, "");
+    reader.prepare(table);
+    PointViewSet viewSet = reader.execute(table);
+    assert(viewSet.size() == 1);
+    return *viewSet.begin();
+}
+
+void NeighborClassifierFilter::filter(PointView& view)
+{
+    PointRef point_src(view, 0);
+    if (m_candidateFile.empty())
+    {   // No candidate file so NN comes from src file
+        KD3Index kdiSrc(view);
+        kdiSrc.build();
+        PointRef point_nn(view, 0);
+        for (PointId id = 0; id < view.size(); ++id)
+        {
+            point_src.setPointId(id);
+            doOne(point_src, point_nn, kdiSrc);
+        }
+    }
+    else
+    {   // NN comes from candidate file
+        PointTable candTable;
+        PointViewPtr candView = loadSet(m_candidateFile, candTable);
+        KD3Index kdiCand(*candView);
+        kdiCand.build();
+        PointRef point_nn(*candView, 0);
+        for (PointId id = 0; id < view.size(); ++id)
+        {
+            point_src.setPointId(id);
+            doOne(point_src, point_nn, kdiCand);
+        }
+    }
+}
+
+} // namespace pdal
+

filters/NeighborClassifierFilter.hpp

@@ -0,0 +1,76 @@
+/******************************************************************************
+* Copyright (c) 2017, Hobu Inc. <hobu.inc@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>
+#include <pdal/KDIndex.hpp>
+
+extern "C" int32_t NeighborClassifierFilter_ExitFunc();
+extern "C" PF_ExitFunc NeighborClassifierFilter_InitPlugin();
+
+namespace pdal
+{
+
+struct DimRange;
+
+class PDAL_DLL NeighborClassifierFilter : public Filter
+{
+public:
+    NeighborClassifierFilter();
+    ~NeighborClassifierFilter();
+
+    static void * create();
+    static int32_t destroy(void *);
+    std::string getName() const { return "filters.neighborclassifier"; }
+
+private:
+    virtual void addArgs(ProgramArgs& args);
+    virtual void prepared(PointTableRef table);
+    bool doOne(PointRef& point, PointRef& temp, KD3Index &kdi);
+    virtual void filter(PointView& view);
+    virtual void initialize();
+    void doOneNoDomain(PointRef &point, PointRef& temp, KD3Index &kdi);
+    PointViewPtr loadSet(const std::string &candFileName, PointTable &table);
+    NeighborClassifierFilter& operator=(const NeighborClassifierFilter&) = delete;
+    NeighborClassifierFilter(const NeighborClassifierFilter&) = delete;
+    StringList m_domainSpec;
+    std::vector<DimRange> m_domain;
+    int m_k;
+    Dimension::Id m_dim;
+    std::string m_dimName;
+    std::string m_candidateFile;
+};
+
+} // namespace pdal

pdal/StageFactory.cpp

@@ -55,6 +55,7 @@
 #include <filters/HAGFilter.hpp>
 #include <filters/HeadFilter.hpp>
 #include <filters/IQRFilter.hpp>
+#include <filters/NeighborClassifierFilter.hpp>
 #include <filters/KDistanceFilter.hpp>
 #include <filters/LocateFilter.hpp>
 #include <filters/LOFFilter.hpp>
@@ -291,6 +292,7 @@ StageFactory::StageFactory(bool no_plugins)
     PluginManager<Stage>::initializePlugin(HeadFilter_InitPlugin);
     PluginManager<Stage>::initializePlugin(IQRFilter_InitPlugin);
     PluginManager<Stage>::initializePlugin(KDistanceFilter_InitPlugin);
+    PluginManager<Stage>::initializePlugin(NeighborClassifierFilter_InitPlugin);
     PluginManager<Stage>::initializePlugin(LocateFilter_InitPlugin);
     PluginManager<Stage>::initializePlugin(LOFFilter_InitPlugin);
     PluginManager<Stage>::initializePlugin(MADFilter_InitPlugin);

test/unit/CMakeLists.txt

@@ -114,6 +114,7 @@ PDAL_ADD_TEST(pdal_filters_decimation_test FILES
 PDAL_ADD_TEST(pdal_filters_divider_test FILES filters/DividerFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_ferry_test FILES filters/FerryFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_groupby_test FILES filters/GroupByFilterTest.cpp)
+PDAL_ADD_TEST(pdal_filters_neighborclassifier_test FILES filters/NeighborClassifierFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_locate_test FILES filters/LocateFilterTest.cpp)
 PDAL_ADD_TEST(pdal_filters_merge_test FILES filters/MergeTest.cpp)
 PDAL_ADD_TEST(pdal_morton_order_test FILES filters/MortonOrderTest.cpp)