Add VoxelCenterNearestNeighborFilter and VoxelCentroidNearestNeighbor…

…Filter (#1603) For each populated voxel, use the coordinates of the voxel center or centroid respectively to find the nearest neighbor in the input PointView and append it to the output PointView.
PDAL · Jun 19, 2017 · dc428f4 · dc428f4
1 parent 30c0e9a
commit dc428f4
Show file tree

Hide file tree

Showing 7 changed files with 464 additions and 0 deletions.
diff --git a/doc/stages/filters.voxelcenternearestneighbor.rst b/doc/stages/filters.voxelcenternearestneighbor.rst
@@ -0,0 +1,49 @@
+.. _filters.voxelcenternearestneighbor:
+
+===============================================================================
+filters.voxelcenternearestneighbor
+===============================================================================
+
+VoxelCenterNearestNeighbor is a voxel-based sampling filter. The input point
+cloud is divided into 3D voxels at the given cell size. For each populated
+voxel, the coordinates of the voxel center are used as the query point in a 3D
+nearest neighbor search. The nearest neighbor is then added to the output point
+cloud, along with any existing dimensions.
+
+.. note::
+
+    This is similar to the existing :ref:`filters.voxelgrid`. However, in the
+    case of the VoxelGrid, the centroid of the points falling within the voxel
+    is added to the output point cloud. The drawback with this approach is that
+    all dimensional data is lost, and the sampled cloud now consists of only XYZ
+    coordinates.
+
+Example
+-------
+
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "input.las",
+        {
+          "type":"filters.voxelcenternearestneighbor",
+          "cell":10.0
+        },
+        "output.las"
+      ]
+    }
+
+
+.. seealso::
+
+    :ref:`filters.voxelcentroidnearestneighbor` offers a similar solution, using
+    as the query point the centroid of all points falling within the voxel as
+    opposed to the voxel center coordinates.
+
+Options
+-------------------------------------------------------------------------------
+
+cell
+  Cell size in the X, Y, and Z dimension. [Default: **1.0**]
diff --git a/doc/stages/filters.voxelcentroidnearestneighbor.rst b/doc/stages/filters.voxelcentroidnearestneighbor.rst
@@ -0,0 +1,47 @@
+.. _filters.voxelcentroidnearestneighbor:
+
+===============================================================================
+filters.voxelcentroidnearestneighbor
+===============================================================================
+
+VoxelCentroidNearestNeighbor is a voxel-based sampling filter. The input point
+cloud is divided into 3D voxels at the given cell size. For each populated
+voxel, the centroid of the points within that voxel is computed. This centroid
+is used as the query point in a 3D nearest neighbor search. The nearest neighbor
+is then added to the output point cloud, along with any existing dimensions.
+
+.. note::
+
+    This is similar to the existing :ref:`filters.voxelgrid`. However, in the
+    case of the VoxelGrid, the centroid itself is added to the output point
+    cloud. The drawback with this approach is that all dimensional data is lost,
+    and the sampled cloud now consists of only XYZ coordinates.
+
+Example
+-------
+
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "input.las",
+        {
+          "type":"filters.voxelcentroidnearestneighbor",
+          "cell":10.0
+        },
+        "output.las"
+      ]
+    }
+
+
+.. seealso::
+
+    :ref:`filters.voxelcenternearestneighbor` offers a similar solution, using
+    the voxel center as opposed to the voxel centroid for the query point.
+
+Options
+-------------------------------------------------------------------------------
+
+cell
+  Cell size in the X, Y, and Z dimension. [Default: **1.0**]
diff --git a/filters/VoxelCenterNearestNeighborFilter.cpp b/filters/VoxelCenterNearestNeighborFilter.cpp
@@ -0,0 +1,107 @@
+/******************************************************************************
+ * Copyright (c) 2017, 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 "VoxelCenterNearestNeighborFilter.hpp"
+
+#include <pdal/KDIndex.hpp>
+#include <pdal/pdal_macros.hpp>
+
+#include <set>
+#include <string>
+#include <tuple>
+#include <vector>
+
+namespace pdal
+{
+
+static PluginInfo const s_info =
+    PluginInfo("filters.voxelcenternearestneighbor",
+               "Voxel Center Nearest Neighbor Filter",
+               "http://pdal.io/stages/filters.voxelcenternearestneighbor.html");
+
+CREATE_STATIC_PLUGIN(1, 0, VoxelCenterNearestNeighborFilter, Filter, s_info)
+
+std::string VoxelCenterNearestNeighborFilter::getName() const
+{
+    return s_info.name;
+}
+
+void VoxelCenterNearestNeighborFilter::addArgs(ProgramArgs& args)
+{
+    args.add("cell", "Cell size", m_cell, 1.0);
+}
+
+PointViewSet VoxelCenterNearestNeighborFilter::run(PointViewPtr view)
+{
+    BOX3D bounds;
+    view->calculateBounds(bounds);
+
+    KD3Index kdi(*view);
+    kdi.build();
+
+    // Make an initial pass through the input PointView to detect populated
+    // voxels.
+    std::set<std::tuple<size_t, size_t, size_t>> populated_voxels;
+    for (PointId id = 0; id < view->size(); ++id)
+    {
+        double y = view->getFieldAs<double>(Dimension::Id::Y, id);
+        double x = view->getFieldAs<double>(Dimension::Id::X, id);
+        double z = view->getFieldAs<double>(Dimension::Id::Z, id);
+        size_t r = static_cast<size_t>(floor(y - bounds.miny) / m_cell);
+        size_t c = static_cast<size_t>(floor(x - bounds.minx) / m_cell);
+        size_t d = static_cast<size_t>(floor(z - bounds.minz) / m_cell);
+        populated_voxels.emplace(std::make_tuple(r, c, d));
+    }
+
+    // Make a second pass through the populated voxels to find the nearest
+    // neighbor to each voxel center.
+    PointViewPtr output = view->makeNew();
+    for (auto const& t : populated_voxels)
+    {
+        auto& r = std::get<0>(t);
+        auto& c = std::get<1>(t);
+        auto& d = std::get<2>(t);
+        double y = bounds.miny + (r + 0.5) * m_cell;
+        double x = bounds.minx + (c + 0.5) * m_cell;
+        double z = bounds.minz + (d + 0.5) * m_cell;
+        std::vector<PointId> neighbors = kdi.neighbors(x, y, z, 1);
+        output->appendPoint(*view, neighbors[0]);
+    }
+
+    PointViewSet viewSet;
+    viewSet.insert(output);
+    return viewSet;
+}
+
+} // namespace pdal
diff --git a/filters/VoxelCenterNearestNeighborFilter.hpp b/filters/VoxelCenterNearestNeighborFilter.hpp
@@ -0,0 +1,75 @@
+/******************************************************************************
+ * Copyright (c) 2017, 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>
+#include <pdal/plugin.hpp>
+
+#include <cstdint>
+#include <string>
+
+extern "C" int32_t VoxelCenterNearestNeighborFilter_ExitFunc();
+extern "C" PF_ExitFunc VoxelCenterNearestNeighborFilter_InitPlugin();
+
+namespace pdal
+{
+
+class PointLayout;
+class PointView;
+
+class PDAL_DLL VoxelCenterNearestNeighborFilter : public Filter
+{
+public:
+    VoxelCenterNearestNeighborFilter() : Filter()
+    {
+    }
+
+    static void* create();
+    static int32_t destroy(void*);
+    std::string getName() const;
+
+private:
+    double m_cell;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual PointViewSet run(PointViewPtr view);
+
+    VoxelCenterNearestNeighborFilter&
+    operator=(const VoxelCenterNearestNeighborFilter&); // not implemented
+    VoxelCenterNearestNeighborFilter(
+        const VoxelCenterNearestNeighborFilter&); // not implemented
+};
+
+} // namespace pdal
diff --git a/filters/VoxelCentroidNearestNeighborFilter.cpp b/filters/VoxelCentroidNearestNeighborFilter.cpp
@@ -0,0 +1,107 @@
+/******************************************************************************
+ * Copyright (c) 2017, 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 "VoxelCentroidNearestNeighborFilter.hpp"
+
+#include <pdal/EigenUtils.hpp>
+#include <pdal/KDIndex.hpp>
+#include <pdal/pdal_macros.hpp>
+
+#include <map>
+#include <string>
+#include <tuple>
+#include <vector>
+
+#include <Eigen/Dense>
+
+namespace pdal
+{
+
+static PluginInfo const s_info = PluginInfo(
+    "filters.voxelcentroidnearestneighbor",
+    "Voxel Centroid Nearest Neighbor Filter",
+    "http://pdal.io/stages/filters.voxelcentroidnearestneighbor.html");
+
+CREATE_STATIC_PLUGIN(1, 0, VoxelCentroidNearestNeighborFilter, Filter, s_info)
+
+std::string VoxelCentroidNearestNeighborFilter::getName() const
+{
+    return s_info.name;
+}
+
+void VoxelCentroidNearestNeighborFilter::addArgs(ProgramArgs& args)
+{
+    args.add("cell", "Cell size", m_cell, 1.0);
+}
+
+PointViewSet VoxelCentroidNearestNeighborFilter::run(PointViewPtr view)
+{
+    BOX3D bounds;
+    view->calculateBounds(bounds);
+
+    KD3Index kdi(*view);
+    kdi.build();
+
+    // Make an initial pass through the input PointView to index PointIds by
+    // row, column, and depth.
+    std::map<std::tuple<size_t, size_t, size_t>, std::vector<PointId>>
+        populated_voxel_ids;
+    for (PointId id = 0; id < view->size(); ++id)
+    {
+        double y = view->getFieldAs<double>(Dimension::Id::Y, id);
+        double x = view->getFieldAs<double>(Dimension::Id::X, id);
+        double z = view->getFieldAs<double>(Dimension::Id::Z, id);
+        size_t r = static_cast<size_t>(floor(y - bounds.miny) / m_cell);
+        size_t c = static_cast<size_t>(floor(x - bounds.minx) / m_cell);
+        size_t d = static_cast<size_t>(floor(z - bounds.minz) / m_cell);
+        populated_voxel_ids[std::make_tuple(r, c, d)].push_back(id);
+    }
+
+    // Make a second pass through the populated voxels to compute the voxel
+    // centroid and to find its nearest neighbor.
+    PointViewPtr output = view->makeNew();
+    for (auto const& t : populated_voxel_ids)
+    {
+        Eigen::Vector3f centroid = eigen::computeCentroid(*view, t.second);
+        std::vector<PointId> neighbors =
+            kdi.neighbors(centroid[0], centroid[1], centroid[2], 1);
+        output->appendPoint(*view, neighbors[0]);
+    }
+
+    PointViewSet viewSet;
+    viewSet.insert(output);
+    return viewSet;
+}
+
+} // namespace pdal