From 389ae91de0fac67064eccaf83417263edc157c7e Mon Sep 17 00:00:00 2001 From: Bradley J Chambers Date: Thu, 1 Sep 2016 13:36:04 -0400 Subject: [PATCH 1/2] Add MAD and IQR filters with user-selectable dimension and multiplier --- doc/stages/filters.iqr.rst | 51 +++++++++++++++ doc/stages/filters.mad.rst | 47 ++++++++++++++ filters/CMakeLists.txt | 2 + filters/iqr/CMakeLists.txt | 2 + filters/iqr/IQRFilter.cpp | 124 +++++++++++++++++++++++++++++++++++++ filters/iqr/IQRFilter.hpp | 75 ++++++++++++++++++++++ filters/mad/CMakeLists.txt | 2 + filters/mad/MADFilter.cpp | 124 +++++++++++++++++++++++++++++++++++++ filters/mad/MADFilter.hpp | 75 ++++++++++++++++++++++ src/StageFactory.cpp | 4 ++ 10 files changed, 506 insertions(+) create mode 100644 doc/stages/filters.iqr.rst create mode 100644 doc/stages/filters.mad.rst create mode 100644 filters/iqr/CMakeLists.txt create mode 100644 filters/iqr/IQRFilter.cpp create mode 100644 filters/iqr/IQRFilter.hpp create mode 100644 filters/mad/CMakeLists.txt create mode 100644 filters/mad/MADFilter.cpp create mode 100644 filters/mad/MADFilter.hpp diff --git a/doc/stages/filters.iqr.rst b/doc/stages/filters.iqr.rst new file mode 100644 index 0000000000..ba34fe7c78 --- /dev/null +++ b/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. diff --git a/doc/stages/filters.mad.rst b/doc/stages/filters.mad.rst new file mode 100644 index 0000000000..e9bbcf7da6 --- /dev/null +++ b/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. diff --git a/filters/CMakeLists.txt b/filters/CMakeLists.txt index 78d0b25f6f..67473ed9b5 100644 --- a/filters/CMakeLists.txt +++ b/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) diff --git a/filters/iqr/CMakeLists.txt b/filters/iqr/CMakeLists.txt new file mode 100644 index 0000000000..5a1c3ae01a --- /dev/null +++ b/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) diff --git a/filters/iqr/IQRFilter.cpp b/filters/iqr/IQRFilter.cpp new file mode 100644 index 0000000000..a2952abd8c --- /dev/null +++ b/filters/iqr/IQRFilter.cpp @@ -0,0 +1,124 @@ +/****************************************************************************** +* 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 + +#include +#include + +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); +} + +// ready or prepared to make sure dim exists? does addArgs come first? (i think so) +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 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 z(view->size()); + for (PointId j = 0; j < view->size(); ++j) + z[j] = view->getFieldAs(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(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 diff --git a/filters/iqr/IQRFilter.hpp b/filters/iqr/IQRFilter.hpp new file mode 100644 index 0000000000..03e9c54d37 --- /dev/null +++ b/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 +#include + +#include + +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 diff --git a/filters/mad/CMakeLists.txt b/filters/mad/CMakeLists.txt new file mode 100644 index 0000000000..c695efb363 --- /dev/null +++ b/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) diff --git a/filters/mad/MADFilter.cpp b/filters/mad/MADFilter.cpp new file mode 100644 index 0000000000..2a95acacac --- /dev/null +++ b/filters/mad/MADFilter.cpp @@ -0,0 +1,124 @@ +/****************************************************************************** +* 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 "MADFilter.hpp" + +#include + +#include +#include + +namespace pdal +{ + +static PluginInfo const s_info = + PluginInfo("filters.mad", "Median Absolute Deviation Filter", + "http://pdal.io/stages/filters.mad.html"); + +CREATE_STATIC_PLUGIN(1, 0, MADFilter, Filter, s_info) + +std::string MADFilter::getName() const +{ + return s_info.name; +} + +void MADFilter::addArgs(ProgramArgs& args) +{ + args.add("k", "Number of deviations", m_multiplier, 2.0); + args.add("dimension", "Dimension on which to calculate statistics", + m_dimName); +} + +// ready or prepared to make sure dim exists? does addArgs come first? (i think so) +void MADFilter::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.mad 'dimension' " + "option: '" << m_dimName << "'."; + throw pdal_error(oss.str()); + } +} + +PointViewSet MADFilter::run(PointViewPtr view) +{ + using namespace Dimension; + + PointViewSet viewSet; + PointViewPtr output = view->makeNew(); + + auto median = [](std::vector vals) + { + std::nth_element(vals.begin(), vals.begin()+vals.size()/2, vals.end()); + + return *(vals.begin()+vals.size()/2); + }; + + std::vector z(view->size()); + for (PointId j = 0; j < view->size(); ++j) + z[j] = view->getFieldAs(m_dimId, j); + + double m = median(z); + log()->get(LogLevel::Debug) << "Median value: " << m << std::endl; + + std::vector a(view->size()); + for (PointId j = 0; j < view->size(); ++j) + a[j] = std::fabs(view->getFieldAs(Id::Z, j)-m); + + double mad = median(a)*1.4862; + log()->get(LogLevel::Debug) << "MAD: " << mad << std::endl; + + for (PointId j = 0; j < view->size(); ++j) + { + if (a[j]/mad < m_multiplier) + output->appendPoint(*view, j); + } + + double low_fence = m - m_multiplier * mad; + double hi_fence = m + m_multiplier * mad; + + 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 diff --git a/filters/mad/MADFilter.hpp b/filters/mad/MADFilter.hpp new file mode 100644 index 0000000000..eee699810a --- /dev/null +++ b/filters/mad/MADFilter.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 +#include + +#include + +extern "C" int32_t MADFilter_ExitFunc(); +extern "C" PF_ExitFunc MADFilter_InitPlugin(); + +namespace pdal +{ + +class ProgramArgs; +class PointTable; +class PointView; + +class PDAL_DLL MADFilter : public Filter +{ +public: + MADFilter() : 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); + + MADFilter& operator=(const MADFilter&); // not implemented + MADFilter(const MADFilter&); // not implemented +}; + +} // namespace pdal diff --git a/src/StageFactory.cpp b/src/StageFactory.cpp index ad4d1f2740..289feeecc6 100644 --- a/src/StageFactory.cpp +++ b/src/StageFactory.cpp @@ -48,6 +48,8 @@ #include #include #include +#include +#include #include #include #include @@ -235,6 +237,8 @@ StageFactory::StageFactory(bool no_plugins) PluginManager::initializePlugin(EstimateRankFilter_InitPlugin); PluginManager::initializePlugin(FerryFilter_InitPlugin); PluginManager::initializePlugin(HAGFilter_InitPlugin); + PluginManager::initializePlugin(IQRFilter_InitPlugin); + PluginManager::initializePlugin(MADFilter_InitPlugin); PluginManager::initializePlugin(MergeFilter_InitPlugin); PluginManager::initializePlugin(MongusFilter_InitPlugin); PluginManager::initializePlugin(MortonOrderFilter_InitPlugin); From 38b45cd569ba576acb75e99218e876636a9d4359 Mon Sep 17 00:00:00 2001 From: Bradley J Chambers Date: Thu, 1 Sep 2016 15:19:25 -0400 Subject: [PATCH 2/2] Remove old comment to self --- filters/iqr/IQRFilter.cpp | 1 - filters/mad/MADFilter.cpp | 1 - 2 files changed, 2 deletions(-) diff --git a/filters/iqr/IQRFilter.cpp b/filters/iqr/IQRFilter.cpp index a2952abd8c..de47a24ce2 100644 --- a/filters/iqr/IQRFilter.cpp +++ b/filters/iqr/IQRFilter.cpp @@ -60,7 +60,6 @@ void IQRFilter::addArgs(ProgramArgs& args) m_dimName); } -// ready or prepared to make sure dim exists? does addArgs come first? (i think so) void IQRFilter::prepared(PointTableRef table) { PointLayoutPtr layout(table.layout()); diff --git a/filters/mad/MADFilter.cpp b/filters/mad/MADFilter.cpp index 2a95acacac..59d79dfb1a 100644 --- a/filters/mad/MADFilter.cpp +++ b/filters/mad/MADFilter.cpp @@ -60,7 +60,6 @@ void MADFilter::addArgs(ProgramArgs& args) m_dimName); } -// ready or prepared to make sure dim exists? does addArgs come first? (i think so) void MADFilter::prepared(PointTableRef table) { PointLayoutPtr layout(table.layout());