Merge pull request #1298 from PDAL/mad-and-iqr-outliers

Add MAD and IQR filters with user-selectable dimension and multiplier

doc/stages/filters.iqr.rst

@@ -0,0 +1,51 @@
+.. _filters.iqr:
+
+filters.iqr
+===============================================================================
+
+The ``filters.iqr`` filter automatically crops the input point cloud based on
+the distribution of points in the specified dimension. Specifically, we choose
+the method of Interquartile Range (IQR). The IQR is defined as the range between
+the first and third quartile (25th and 75th percentile). Upper and lower bounds
+are determined by adding 1.5 times the IQR to the third quartile or subtracting
+1.5 times the IQR from the first quartile. The multiplier, which defaults to
+1.5, can be adjusted by the user.
+
+.. note::
+
+  This method can remove real data, especially ridges and valleys in rugged
+  terrain, or tall features such as towers and rooftops in flat terrain. While
+  the number of deviations can be adjusted to account for such content-specific
+  considerations, it must be used with care.
+
+Example
+-------
+
+The sample pipeline below uses ``filters.iqr`` to automatically crop the Z
+dimension and remove possible outliers. The multiplier to determine high/low
+thresholds has been adjusted to be less agressive and to only crop those
+outliers that are greater than the third quartile plus 3 times the IQR or are
+less than the first quartile minus 3 times the IQR.
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "input.las",
+        {
+          "type":"filters.iqr",
+          "dimension":"Z",
+          "k":3.0
+        },
+        "output.laz"
+      ]
+    }
+
+Options
+-------------------------------------------------------------------------------
+
+k
+  The IQR multiplier used to determine upper/lower bounds. [Default: **1.5**]
+
+dimension
+  The name of the dimension to filter.

doc/stages/filters.mad.rst

@@ -0,0 +1,47 @@
+.. _filters.mad:
+
+filters.mad
+===============================================================================
+
+The ``filters.mad`` filter automatically crops the input point cloud based on
+the distribution of points in the specified dimension. Specifically, we choose
+the method of median absolute deviation from the median (commonly referred to as
+MAD), which is robust to outliers (as opposed to mean and standard deviation).
+
+.. note::
+
+  This method can remove real data, especially ridges and valleys in rugged
+  terrain, or tall features such as towers and rooftops in flat terrain. While
+  the number of deviations can be adjusted to account for such content-specific
+  considerations, it must be used with care.
+
+Example
+-------
+
+The sample pipeline below uses ``filters.mad`` to automatically crop the Z
+dimension and remove possible outliers. The number of deviations from the median
+has been adjusted to be less agressive and to only crop those outliers that are
+greater than four deviations from the median.
+
+.. code-block:: json
+
+    {
+      "pipeline":[
+        "input.las",
+        {
+          "type":"filters.mad",
+          "dimension":"Z",
+          "k":4.0
+        },
+        "output.laz"
+      ]
+    }
+
+Options
+-------------------------------------------------------------------------------
+
+k
+  The number of deviations from the median. [Default: **2.0**]
+
+dimension
+  The name of the dimension to filter.

filters/CMakeLists.txt

@@ -9,6 +9,8 @@ add_subdirectory(eigenvalues)
 add_subdirectory(estimaterank)
 add_subdirectory(ferry)
 add_subdirectory(hag)
+add_subdirectory(iqr)
+add_subdirectory(mad)
 add_subdirectory(merge)
 add_subdirectory(mongus)
 add_subdirectory(mortonorder)

filters/iqr/CMakeLists.txt

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

filters/iqr/IQRFilter.cpp

@@ -0,0 +1,123 @@
+/******************************************************************************
+* 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 "IQRFilter.hpp"
+
+#include <pdal/pdal_macros.hpp>
+
+#include <string>
+#include <vector>
+
+namespace pdal
+{
+
+static PluginInfo const s_info =
+    PluginInfo("filters.iqr", "Interquartile Range Filter",
+               "http://pdal.io/stages/filters.iqr.html");
+
+CREATE_STATIC_PLUGIN(1, 0, IQRFilter, Filter, s_info)
+
+std::string IQRFilter::getName() const
+{
+    return s_info.name;
+}
+
+void IQRFilter::addArgs(ProgramArgs& args)
+{
+    args.add("k", "Number of deviations", m_multiplier, 1.5);
+    args.add("dimension", "Dimension on which to calculate statistics",
+        m_dimName);
+}
+
+void IQRFilter::prepared(PointTableRef table)
+{
+    PointLayoutPtr layout(table.layout());
+    m_dimId = layout->findDim(m_dimName);
+    if (m_dimId == Dimension::Id::Unknown)
+    {
+        std::ostringstream oss;
+        oss << "Invalid dimension name in filters.iqr 'dimension' "
+            "option: '" << m_dimName << "'.";
+        throw pdal_error(oss.str());
+    }
+}
+
+PointViewSet IQRFilter::run(PointViewPtr view)
+{
+    using namespace Dimension;
+
+    PointViewSet viewSet;
+    PointViewPtr output = view->makeNew();
+
+    auto quartile = [](std::vector<double> vals, double percent)
+    {
+        std::nth_element(vals.begin(), vals.begin()+int(vals.size()*percent), vals.end());
+
+        return *(vals.begin()+int(vals.size()*percent));
+    };
+
+    std::vector<double> z(view->size());
+    for (PointId j = 0; j < view->size(); ++j)
+        z[j] = view->getFieldAs<double>(m_dimId, j);
+
+
+    double pc25 = quartile(z, 0.25);
+    log()->get(LogLevel::Debug) << "25th percentile: " << pc25 << std::endl;
+
+    double pc75 = quartile(z, 0.75);
+    log()->get(LogLevel::Debug) << "75th percentile: " << pc75 << std::endl;
+
+    double iqr = pc75-pc25;
+    log()->get(LogLevel::Debug) << "IQR: " << iqr << std::endl;
+
+    double low_fence = pc25 - m_multiplier * iqr;
+    double hi_fence = pc75 + m_multiplier * iqr;
+
+    for (PointId j = 0; j < view->size(); ++j)
+    {
+        double val = view->getFieldAs<double>(m_dimId, j);
+        if (val > low_fence && val < hi_fence)
+            output->appendPoint(*view, j);
+    }
+    log()->get(LogLevel::Debug) << "Cropping " << m_dimName
+                                << " in the range (" << low_fence
+                                << "," << hi_fence << ")" << std::endl;
+
+    viewSet.erase(view);
+    viewSet.insert(output);
+
+    return viewSet;
+}
+
+} // namespace pdal

filters/iqr/IQRFilter.hpp

@@ -0,0 +1,75 @@
+/******************************************************************************
+* 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 <string>
+
+extern "C" int32_t IQRFilter_ExitFunc();
+extern "C" PF_ExitFunc IQRFilter_InitPlugin();
+
+namespace pdal
+{
+
+class ProgramArgs;
+class PointTable;
+class PointView;
+
+class PDAL_DLL IQRFilter : public Filter
+{
+public:
+    IQRFilter() : Filter()
+    {}
+
+    static void * create();
+    static int32_t destroy(void *);
+    std::string getName() const;
+
+private:
+    double m_multiplier;
+    std::string m_dimName;
+    Dimension::Id m_dimId;
+
+    virtual void addArgs(ProgramArgs& args);
+    virtual void prepared(PointTableRef table);
+    virtual PointViewSet run(PointViewPtr view);
+
+    IQRFilter& operator=(const IQRFilter&); // not implemented
+    IQRFilter(const IQRFilter&); // not implemented
+};
+
+} // namespace pdal

filters/mad/CMakeLists.txt

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