Add outlier detection criterion and make minor edit to filters.lof do…

…cumentation (#2504)

* Add outlier detection criterion and make minor edit to filters.lof documentation

This commit adds three filters that are implementations of the outlier
detection criterion specified in Weyrich, T et al. “Post-Processing of Scanned
3D Surface Data.” Proceedings of Eurographics Symposium on Point-Based Graphics
2004 (2004): 85–94. Print.

* Nearest-Neighbor Reciprocity Criterion available as filters.reciprocity

* Miniball Criterion available as filters.miniball

* Plane Fit Criterion available as filters.planefit

* Removing debug header enables us to compile on Windows

* Update Eigen matrix from float to double

doc/references.rst

@@ -34,6 +34,8 @@ Reference
 
 .. [Alexa2003] Alexa, Marc, et al. "Computing and rendering point set surfaces." Visualization and Computer Graphics, IEEE Transactions on 9.1 (2003): 3-15.
 
+.. [Breunig2000] Breunig, M.M., Kriegel, H.-P., Ng, R.T., Sander, J., 2000. LOF: Identifying Density-Based Local Outliers. Proc. 2000 Acm Sigmod Int. Conf. Manag. Data 1–12.
+
 .. [Chen2012] Chen, Ziyue et al. “Upward-Fusion Urban DTM Generating Method Using Airborne Lidar Data.” ISPRS Journal of Photogrammetry and Remote Sensing 72 (2012): 121–130.
 
 .. [Cook1986] Cook, Robert L. "Stochastic sampling in computer graphics." *ACM Transactions on Graphics (TOG)* 5.1 (1986): 51-72.
@@ -42,6 +44,8 @@ Reference
 
 .. [Dippe1985] Dippé, Mark AZ, and Erling Henry Wold. "Antialiasing through stochastic sampling." *ACM Siggraph Computer Graphics* 19.3 (1985): 69-78.
 
+.. [Fischer2010] Fischer, Kaspar, Bernd Gärtner, and Martin Kutz. “Fast Smallest-Enclosing-Ball Computation in High Dimensions.” 26473 (2010): 630–641. Web.
+
 .. [Guinard2017] Guinard S., Landrieu L. "Weakly Supervised Segmented-Aided Classification of Urban Scenes From 3D LIDAR Point Clouds." Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLII-1/W1, 151-157, 2017
 
 .. [Kazhdan2006] Kazhdan, Michael, Matthew Bolitho, and Hugues Hoppe. "Poisson surface reconstruction." Proceedings of the fourth Eurographics symposium on Geometry processing. Vol. 7. 2006.
@@ -56,4 +60,6 @@ Reference
 
 .. [Rusu2008] Rusu, Radu Bogdan, et al. "Towards 3D point cloud based object maps for household environments." Robotics and Autonomous Systems 56.11 (2008): 927-941.
 
+.. [Weyrich2004] Weyrich, T et al. “Post-Processing of Scanned 3D Surface Data.” Proceedings of Eurographics Symposium on Point-Based Graphics 2004 (2004): 85–94. Print.
+
 .. [Zhang2003] Zhang, Keqi, et al. "A progressive morphological filter for removing nonground measurements from airborne LIDAR data." Geoscience and Remote Sensing, IEEE Transactions on 41.4 (2003): 872-882.

doc/stages/filters.lof.rst

@@ -37,9 +37,9 @@ users of this filter should find instructive.
 Example
 -------
 
-The sample pipeline below uses computes the LOF with a neighborhood of 20
-neighbors, followed by a range filter to crop out points whose
-``LocalOutlierFactor`` exceeds 1.2 before writing the output.
+The sample pipeline below computes the LOF with a neighborhood of 20 neighbors,
+followed by a range filter to crop out points whose ``LocalOutlierFactor``
+exceeds 1.2 before writing the output.
 
 .. code-block:: json
 
@@ -62,5 +62,3 @@ Options
 _`minpts`
   The number of k nearest neighbors. [Default: 10]
 
-.. [Breunig2000] Breunig, M.M., Kriegel, H.-P., Ng, R.T., Sander, J., 2000. LOF: Identifying Density-Based Local Outliers. Proc. 2000 Acm Sigmod Int. Conf. Manag. Data 1–12.
-

doc/stages/filters.miniball.rst

@@ -0,0 +1,53 @@
+.. _filters.miniball:
+
+filters.miniball
+===============================================================================
+
+The **Miniball Criterion** was introduced in [Weyrich2004]_ and is based on the
+assumption that points that are distant to the cluster built by their
+k-neighborhood are likely to be outliers. First, the smallest enclosing ball is
+computed for the k-neighborhood, giving a center point and radius
+[Fischer2010]_. The miniball criterion is then computed by comparing the
+distance (from the current point to the miniball center) to the radius of the
+miniball.
+
+The author suggests that the Miniball Criterion is more robust than the
+:ref:`Plane Fit Criterion <filters.planefit>` around high-frequency details,
+but demonstrates poor outlier detection for points close to a smooth surface.
+
+The filter creates a single new dimension, ``Miniball``, that records the
+Miniball criterion for the current point.
+
+.. note::
+
+  To inspect the newly created, non-standard dimensions, be sure to write to an
+  output format that can support arbitrary dimensions, such as BPF.
+
+.. embed::
+
+Example
+-------
+
+The sample pipeline below computes the Miniball criterion with a neighborhood
+of 8 neighbors. We do not apply a fixed threshold to single out outliers based
+on the Miniball criterion as the range of values can vary from one dataset to
+another. In general, higher values indicate the likelihood of a point being an
+outlier.
+
+.. code-block:: json
+
+  [
+      "input.las",
+      {
+          "type":"filters.miniball",
+          "knn":8
+      }
+      "output.laz"
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+knn
+  The number of k nearest neighbors. [Default: 8]
+

doc/stages/filters.planefit.rst

@@ -0,0 +1,58 @@
+.. _filters.planefit:
+
+filters.planefit
+===============================================================================
+
+The **Plane Fit Criterion** was introduced in [Weyrich2004]_ and computes the
+deviation of a point from a manifold approximating its neighbors.  First, a
+plane is fit to each point's k-neighborhood by performing an eigenvalue
+decomposition. Next, the mean point to plane distance is computed by
+considering all points within the neighborhood. This is compared to the point
+to plane distance of the current point giving rise to the k-neighborhood. As
+the mean distance of the k-neighborhood approaches 0, the Plane Fit criterion
+will tend toward 1. As point to plane distance of the current point approaches
+0, the Plane Fit criterion will tend toward 0.
+
+The author suggests that the Plane Fit Criterion is well suited to outlier
+detection when considering noisy reconstructions of smooth surfaces, but
+produces poor results around small features and creases.
+
+The filter creates a single new dimension, ``PlaneFit``, that records the
+Plane Fit criterion for the current point.
+
+.. note::
+
+  To inspect the newly created, non-standard dimensions, be sure to write to an
+  output format that can support arbitrary dimensions, such as BPF.
+
+.. embed::
+
+Example
+-------
+
+The sample pipeline below computes the Plane Fit criterion with a neighborhood
+of 8 neighbors. We do not apply a fixed threshold to single out outliers based
+on the Plane Fit criterion as the range of values can vary from one dataset to
+another. In general, higher values indicate the likelihood of a point being an
+outlier.
+
+.. code-block:: json
+
+  [
+      "input.las",
+      {
+          "type":"filters.planefit",
+          "knn":8
+      }
+      "output.laz"
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+knn
+  The number of k nearest neighbors. [Default: 8]
+
+threads
+  The number of threads used for computing the plane fit criterion. [Default: 1]
+

doc/stages/filters.reciprocity.rst

@@ -0,0 +1,55 @@
+.. _filters.reciprocity:
+
+filters.reciprocity
+===============================================================================
+
+The **Nearest-Neighbor Reciprocity Criterion** was introduced in [Weyrich2004]_
+and is based on a simple assumption, that valid points may be in the
+k-neighborhood of an outlier, but the outlier will most likely not be part of
+the valid point's k-neighborhood.
+
+The author suggests that the Nearest-Neighbor Reciprocity Criterion is more
+robust than both the :ref:`Plane Fit <filters.planefit>` and :ref:`Miniball
+<filters.miniball>` Criterion, being equally sensitive around smooth and
+detailed regions. The criterion does however produce invalid reslts near
+manifold borders.
+
+The filter creates a single new dimension, ``Reciprocity``, that records the
+percentage of points(in the range 0 to 100) that are considered uni-directional
+neighbors of the current point. 
+
+.. note::
+
+  To inspect the newly created, non-standard dimensions, be sure to write to an
+  output format that can support arbitrary dimensions, such as BPF.
+
+.. embed::
+
+Example
+-------
+
+The sample pipeline below computes reciprocity with a neighborhood of 8
+neighbors, followed by a range filter to crop out points whose ``Reciprocity``
+percentage is less than 98% before writing the output.
+
+.. code-block:: json
+
+  [
+      "input.las",
+      {
+          "type":"filters.reciprocity",
+          "knn":8
+      },
+      {
+          "type":"filters.range",
+          "limits":"Reciprocity[:98.0]"
+      },
+      "output.laz"
+  ]
+
+Options
+-------------------------------------------------------------------------------
+
+knn
+  The number of k nearest neighbors. [Default: 8]
+

filters/MiniballFilter.cpp

@@ -0,0 +1,151 @@
+/******************************************************************************
+ * Copyright (c) 2019, 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.
+ ****************************************************************************/
+
+// PDAL implementation of the miniball criterion presented in T. Weyrich, M.
+// Pauly, R. Keiser, S. Heinzle, S. Scandella, and M. Gross, “Post-processing
+// of Scanned 3D Surface Data,” Proc. Eurographics Symp.  Point-Based Graph.
+// 2004, pp. 85–94, 2004.
+
+#include "MiniballFilter.hpp"
+
+#include <pdal/KDIndex.hpp>
+#include <pdal/util/ProgramArgs.hpp>
+
+#include "private/miniball/Seb.h"
+
+#include <cmath>
+#include <string>
+#include <thread>
+#include <vector>
+
+namespace pdal
+{
+
+using namespace Dimension;
+
+static StaticPluginInfo const s_info
+{
+    "filters.miniball",
+    "Miniball (Kutz et al., 2003)",
+    "http://pdal.io/stages/filters.miniball.html"
+};
+
+CREATE_STATIC_STAGE(MiniballFilter, s_info)
+
+std::string MiniballFilter::getName() const
+{
+    return s_info.name;
+}
+
+void MiniballFilter::addArgs(ProgramArgs& args)
+{
+    args.add("knn", "k-Nearest neighbors", m_knn, 8);
+    args.add("threads", "Number of threads used to run this filter", m_threads,
+             1);
+}
+
+void MiniballFilter::addDimensions(PointLayoutPtr layout)
+{
+    m_miniball =
+        layout->registerOrAssignDim("Miniball", Dimension::Type::Double);
+}
+
+void MiniballFilter::filter(PointView& view)
+{
+    KD3Index& kdi = view.build3dIndex();
+
+    point_count_t nloops = view.size();
+    std::vector<std::thread> threadPool(m_threads);
+    for (int t = 0; t < m_threads; t++)
+    {
+        threadPool[t] = std::thread(std::bind(
+            [&](const PointId start, const PointId end, const PointId t) {
+                for (PointId i = start; i < end; i++)
+                    setMiniball(view, i, kdi);
+            },
+            t * nloops / m_threads,
+            (t + 1) == m_threads ? nloops : (t + 1) * nloops / m_threads, t));
+    }
+    for (auto& t : threadPool)
+        t.join();
+}
+
+void MiniballFilter::setMiniball(PointView& view, const PointId& i,
+                                 const KD3Index& kdi)
+{
+    typedef double FT;
+    typedef Seb::Point<FT> Point;
+    typedef std::vector<Point> PointVector;
+    typedef Seb::Smallest_enclosing_ball<FT> Miniball;
+
+    double X = view.getFieldAs<double>(Dimension::Id::X, i);
+    double Y = view.getFieldAs<double>(Dimension::Id::Y, i);
+    double Z = view.getFieldAs<double>(Dimension::Id::Z, i);
+
+    // Find k-nearest neighbors of i.
+    auto ni = kdi.neighbors(i, m_knn + 1);
+
+    PointVector S;
+    std::vector<double> coords(3);
+    for (auto const& j : ni)
+    {
+        if (j == i)
+            continue;
+        coords[0] = view.getFieldAs<double>(Dimension::Id::X, j);
+        coords[1] = view.getFieldAs<double>(Dimension::Id::Y, j);
+        coords[2] = view.getFieldAs<double>(Dimension::Id::Z, j);
+        S.push_back(Point(3, coords.begin()));
+    }
+
+    // add neighbors to Miniball mb(3, S)
+    Miniball mb(3, S);
+
+    // obtain radius r = mb.radius();
+    FT radius = mb.radius();
+
+    // obtain center = mb.center_begin()
+    Miniball::Coordinate_iterator center_it = mb.center_begin();
+    double x = center_it[0];
+    double y = center_it[1];
+    double z = center_it[2];
+
+    // compute distance d from p to center
+    double d =
+        std::sqrt((X - x) * (X - x) + (Y - y) * (Y - y) + (Z - z) * (Z - z));
+
+    double miniball = d / (d + 2 * radius / (std::sqrt(3)));
+    view.setField(m_miniball, i, miniball);
+}
+
+} // namespace pdal