Add a Radial Density filter (#1373)

This filter creates a new dimension called RadialDensity that records the
density of points within a sphere of specified radius and centered at the
current point.

doc/stages/filters.radialdensity.rst

@@ -0,0 +1,47 @@
+.. _filters.radialdensity:
+
+===============================================================================
+filters.radialdensity
+===============================================================================
+
+The Radial Density filter creates a new attribute ``RadialDensity`` that
+contains the density of points in a sphere of given radius.
+
+The density at each point is computed by counting the number of points falling
+within a sphere of given radius (default is 1.0) and centered at the current
+point. The number of neighbors (including the query point) is then normalized
+by the volume of the sphere, defined as
+
+.. math::
+
+  V = \frac{4}{3} \pi r^3
+
+The radius :math:`r` can be adjusted by changing the ``radius`` option.
+
+Example
+-------------------------------------------------------------------------------
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "input.las",
+        {
+          "type":"filters.radialdensity",
+          "radius":2.0
+        },
+        {
+          "type":"writers.bpf",
+          "filename":"output.bpf",
+          "output_dims":"X,Y,Z,RadialDensity"
+        }
+      ]
+    }
+
+
+Options
+-------------------------------------------------------------------------------
+
+radius
+  Radius. [Default: **1.0**]
+

filters/CMakeLists.txt

@@ -19,6 +19,7 @@ add_subdirectory(mortonorder)
 add_subdirectory(normal)
 add_subdirectory(outlier)
 add_subdirectory(pmf)
+add_subdirectory(radialdensity)
 add_subdirectory(randomize)
 add_subdirectory(range)
 add_subdirectory(reprojection)

filters/radialdensity/CMakeLists.txt

@@ -0,0 +1,2 @@
+PDAL_ADD_DRIVER(filter radialdensity "RadialDensityFilter.cpp" "RadialDensityFilter.hpp" objects)
+set(PDAL_TARGET_OBJECTS ${PDAL_TARGET_OBJECTS} ${objects} PARENT_SCOPE)

filters/radialdensity/RadialDensityFilter.cpp

@@ -0,0 +1,92 @@
+/******************************************************************************
+* 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 "RadialDensityFilter.hpp"
+
+#include <pdal/KDIndex.hpp>
+#include <pdal/pdal_macros.hpp>
+
+#include <string>
+#include <vector>
+
+namespace pdal
+{
+
+static PluginInfo const s_info =
+    PluginInfo("filters.radialdensity", "RadialDensity Filter",
+               "http://pdal.io/stages/filters.radialdensity.html");
+
+CREATE_STATIC_PLUGIN(1, 0, RadialDensityFilter, Filter, s_info)
+
+std::string RadialDensityFilter::getName() const
+{
+    return s_info.name;
+}
+
+void RadialDensityFilter::addArgs(ProgramArgs& args)
+{
+    args.add("radius", "Radius", m_rad, 1.0);
+}
+
+void RadialDensityFilter::addDimensions(PointLayoutPtr layout)
+{
+    using namespace Dimension;
+    m_rdens = layout->registerOrAssignDim("RadialDensity", Type::Double);
+}
+
+void RadialDensityFilter::filter(PointView& view)
+{
+    using namespace Dimension;
+
+    // Build the 3D KD-tree.
+    log()->get(LogLevel::Debug) << "Building 3D KD-tree...\n";
+    KD3Index index(view);
+    index.build();
+
+    // Search for neighboring points within the specified radius. The number of
+    // neighbors (which includes the query point) is normalized by the volume
+    // of the search sphere and recorded as the density.
+    log()->get(LogLevel::Debug) << "Computing densities...\n";
+    double factor = 1.0 / ((4.0 / 3.0) * 3.14159 * (m_rad * m_rad * m_rad));
+    for (PointId i = 0; i < view.size(); ++i)
+    {
+        double x = view.getFieldAs<double>(Id::X, i);
+        double y = view.getFieldAs<double>(Id::Y, i);
+        double z = view.getFieldAs<double>(Id::Z, i);
+        std::vector<PointId> pts = index.radius(x, y, z, m_rad);
+        view.setField(m_rdens, i, pts.size() * factor);
+    } 
+}
+
+} // namespace pdal

filters/radialdensity/RadialDensityFilter.hpp

@@ -0,0 +1,74 @@
+/******************************************************************************
+* 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>
+#include <pdal/plugin.hpp>
+
+#include <memory>
+
+extern "C" int32_t RadialDensityFilter_ExitFunc();
+extern "C" PF_ExitFunc RadialDensityFilter_InitPlugin();
+
+namespace pdal
+{
+
+class PointLayout;
+class PointView;
+class ProgramArgs;
+
+class PDAL_DLL RadialDensityFilter : public Filter
+{
+public:
+    RadialDensityFilter() : Filter()
+    {}
+
+    static void * create();
+    static int32_t destroy(void *);
+    std::string getName() const;
+
+private:
+    Dimension::Id m_rdens;
+    double m_rad;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void addDimensions(PointLayoutPtr layout);
+    virtual void filter(PointView& view);
+
+    RadialDensityFilter& operator=(const RadialDensityFilter&); // not implemented
+    RadialDensityFilter(const RadialDensityFilter&); // not implemented
+};
+
+} // namespace pdal

src/StageFactory.cpp

@@ -58,6 +58,7 @@
 #include <normal/NormalFilter.hpp>
 #include <outlier/OutlierFilter.hpp>
 #include <pmf/PMFFilter.hpp>
+#include <radialdensity/RadialDensityFilter.hpp>
 #include <range/RangeFilter.hpp>
 #include <reprojection/ReprojectionFilter.hpp>
 #include <sample/SampleFilter.hpp>
@@ -251,6 +252,7 @@ StageFactory::StageFactory(bool no_plugins)
     PluginManager::initializePlugin(NormalFilter_InitPlugin);
     PluginManager::initializePlugin(OutlierFilter_InitPlugin);
     PluginManager::initializePlugin(PMFFilter_InitPlugin);
+    PluginManager::initializePlugin(RadialDensityFilter_InitPlugin);
     PluginManager::initializePlugin(RangeFilter_InitPlugin);
     PluginManager::initializePlugin(ReprojectionFilter_InitPlugin);
     PluginManager::initializePlugin(SampleFilter_InitPlugin);