Apply streamable Poisson concept directly to SampleFilter

* No need for new VoxelPoisson filter

* Still not happy with a visible artifact that appears to be along voxel
  boundaries

doc/stages/filters.sample.rst

@@ -8,39 +8,55 @@ The practice of performing Poisson sampling via "Dart Throwing" was introduced
 in the mid-1980's by [Cook1986]_ and [Dippe1985]_, and has been applied to
 point clouds in other software [Mesh2009]_.
 
-The sampling can be performed in a single pass through the point cloud.
-To begin,
-each input point is assumed to be kept. As we iterate through the kept points,
-we retrieve all neighbors within a given ``radius``, and mark these neighbors as
-points to be discarded. All remaining kept points are appended to the output
-``PointView``. The full layout (i.e., the dimensions) of the input ``PointView``
-is kept in tact (the same cannot be said for :ref:`filters.voxelgrid`).
+Our implementation of Poisson sampling is made streamable by voxelizing the
+space and only adding points to the output ``PointView`` if they do not violate
+the minimum distance criterion (as specified by ``radius``). The voxelization
+allows several optimizations, first by checking for existing points within the
+same voxel as the point under consideration, which are mostly likely to
+violate the minimum distance criterion. Furthermore, we can easily visit
+neighboring voxels (limiting the search to those that are populated) without
+the need to create a KD-tree from the entire input ``PointView`` first and
+performing costly spatial searches.
 
 .. seealso::
 
-    :ref:`filters.decimation`, :ref:`filters.fps` and
-    :ref:`filters.relaxationdartthrowing` also perform decimation.
+    :ref:`filters.decimation`, :ref:`filters.fps`,
+    :ref:`filters.relaxationdartthrowing`,
+    :ref:`filters.voxelcenternearestneighbor`,
+    :ref:`filters.voxelcentroidnearestneighbor`, and :ref:`voxeldownsize` also
+    perform decimation.
 
 .. note::
 
-    The ``shuffle`` option does not reorder points in the PointView, but
-    shuffles the order in which the points are visited while processing, which
-    can improve the quality of the result.
+    Starting with PDAL v2.3, the ``filters.sample`` now supports streaming
+    mode. As a result, there is no longer an option to ``shuffle`` points (or
+    to provide a ``seed`` for the shuffle).
+
+.. note::
+
+    Starting with PDAL v2.3, a ``cell`` option has been added that works with
+    the existing ``radius``. The user must provide one or the other, but not
+    both. The provided option will be used to automatically compute the other.
+    The relationship between ``cell`` and ``radius`` is such that the
+    ``radius`` defines the radius of a sphere that circumscribes a voxel with
+    edge length defined by ``cell``.
 
 .. embed::
 
+.. streamable::
+
 Options
 -------------------------------------------------------------------------------
 
-radius
-  Minimum distance between samples. [Default: 1.0]
-
-shuffle
-  Choose whether or not to shuffle order in which points are visited. [Default:
-  true]
+cell
+  Voxel cell size. If ``radius`` is set, ``cell`` is automatically computed
+  such that the cell is circumscribed by the sphere defined by ``radius``.
 
-seed
-  Seed for random number generator, used only with shuffle.
+radius
+  Minimum distance between samples. If ``cell`` is set, ``radius`` is
+  automatically computed to defined a sphere that circumscribes the voxel cell.
+  Whether specified or derived, ``radius`` defines the minimum allowable
+  distance between points.
 
 .. include:: filter_opts.rst
 

filters/SampleFilter.cpp

@@ -1,5 +1,5 @@
 /******************************************************************************
- * Copyright (c) 2016-2017, Bradley J Chambers (brad.chambers@gmail.com)
+ * Copyright (c) 2016-2017, 2020 Bradley J Chambers (brad.chambers@gmail.com)
  *
  * All rights reserved.
  *
@@ -34,25 +34,18 @@
 
 #include "SampleFilter.hpp"
 
-#include <pdal/KDIndex.hpp>
 #include <pdal/util/ProgramArgs.hpp>
-#include <pdal/util/Utils.hpp>
 
-#include <chrono>
-#include <numeric>
-#include <random>
 #include <string>
-#include <vector>
 
 namespace pdal
 {
 
-static StaticPluginInfo const s_info
-{
-    "filters.sample",
-    "Subsampling filter",
-    "http://pdal.io/stages/filters.sample.html"
-};
+using namespace Dimension;
+
+static StaticPluginInfo const s_info{
+    "filters.sample", "Subsampling filter",
+    "http://pdal.io/stages/filters.sample.html"};
 
 CREATE_STATIC_STAGE(SampleFilter, s_info)
 
@@ -61,69 +54,153 @@ std::string SampleFilter::getName() const
     return s_info.name;
 }
 
-
 void SampleFilter::addArgs(ProgramArgs& args)
 {
-    args.add("radius", "Minimum radius", m_radius, 1.0);
-    args.add("shuffle", "Shuffle points prior to sampling?", m_shuffle, true);
-    m_seedArg = &args.add("seed", "Random number generator seed", m_seed);
+    m_cellArg = &args.add("cell", "Cell size", m_cell);
+    m_radiusArg = &args.add("radius", "Minimum radius", m_radius);
+}
+
+void SampleFilter::prepared(PointTableRef table)
+{
+    if (m_cellArg->set() && m_radiusArg->set())
+        throwError("Must set only one of 'cell' or 'radius'.");
+
+    if (!m_cellArg->set() && !m_radiusArg->set())
+        throwError("Must set 'cell' or 'radius' but not both.");
 }
 
+void SampleFilter::ready(PointTableRef)
+{
+    m_populatedVoxels.clear();
+
+    if (m_cellArg->set())
+        m_radius = m_cell / 2.0 * std::sqrt(3.0);
 
-PointViewSet SampleFilter::run(PointViewPtr inView)
+    if (m_radiusArg->set())
+        m_cell = 2.0 * m_radius / std::sqrt(3.0);
+
+    log()->get(LogLevel::Debug)
+        << "cell " << m_cell << ", radius " << m_radius << std::endl;
+}
+
+PointViewSet SampleFilter::run(PointViewPtr view)
 {
-    point_count_t np = inView->size();
+    PointViewPtr output = view->makeNew();
+    for (PointRef point : *view)
+    {
+        if (voxelize(point))
+            output->appendPoint(*view, point.pointId());
+    }
 
-    // Return empty PointViewSet if the input PointView has no points.
-    // Otherwise, make a new output PointView.
     PointViewSet viewSet;
-    if (!np)
-        return viewSet;
-    PointViewPtr outView = inView->makeNew();
+    viewSet.insert(output);
+    return viewSet;
+}
 
-    // Build the 3D KD-tree.
-    KD3Index& index = inView->build3dIndex();
+bool SampleFilter::voxelize(PointRef& point)
+{
+    double x = point.getFieldAs<double>(Id::X);
+    double y = point.getFieldAs<double>(Id::Y);
+    double z = point.getFieldAs<double>(Id::Z);
 
-    PointIdList shuffledIds(np);
-    std::iota(shuffledIds.begin(), shuffledIds.end(), 0);
-    if (m_shuffle)
+    // Use the coordinates of the first point as origin to offset data and
+    // derive integer voxel indices.
+    if (m_populatedVoxels.empty())
     {
-        if (!m_seedArg->set())
-            m_seed =
-                std::chrono::system_clock::now().time_since_epoch().count();
-        std::shuffle(shuffledIds.begin(), shuffledIds.end(),
-                     std::mt19937(m_seed));
+        m_originX = x; // - (m_cell / 2);
+        m_originY = y; // - (m_cell / 2);
+        m_originZ = z; // - (m_cell / 2);
     }
 
-    // All points are marked as kept (1) by default. As they are masked by
-    // neighbors within the user-specified radius, their value is changed to 0.
-    std::vector<int> keep(np, 1);
+    // Get voxel indices for current point.
+    Voxel v = std::make_tuple((int)(std::floor((x - m_originX) / m_cell)),
+                              (int)(std::floor((y - m_originY) / m_cell)),
+                              (int)(std::floor((z - m_originZ) / m_cell)));
 
-    // We are able to subsample in a single pass over the shuffled indices.
-    for (PointId const& i : shuffledIds)
+    // Check current voxel before any of the neighbors. We will most often have
+    // points that are too close in the point's enclosing voxel, thus saving
+    // cycles.
+    if (m_populatedVoxels.find(v) != m_populatedVoxels.end())
     {
-        // If a point is masked, it is forever masked, and cannot be part of the
-        // sampled cloud. Otherwise, the current index is appended to the output
-        // PointView.
-        if (keep[i] == 0)
-            continue;
-        outView->appendPoint(*inView, i);
-
-        // We now proceed to mask all neighbors within m_radius of the kept
-        // point.
-        PointIdList ids = index.radius(i, m_radius);
-        for (PointId const& id : ids)
-            keep[id] = 0;
+        // Get list of coordinates in candidate voxel.
+        CoordList coords = m_populatedVoxels[v];
+        for (Coord const& coord : coords)
+        {
+            // Compute Euclidean distance between current point and
+            // candidate voxel.
+            double xv = std::get<0>(coord);
+            double yv = std::get<1>(coord);
+            double zv = std::get<2>(coord);
+            double dist = std::sqrt((xv - x) * (xv - x) + (yv - y) * (yv - y) +
+                                    (zv - z) * (zv - z));
+
+            // If any point is closer than the minimum radius, we can
+            // immediately return false, as the minimum distance
+            // criterion is violated.
+            if (dist < m_radius)
+                return false;
+        }
     }
 
-    // Simply calculate the percentage of retained points.
-    double frac = (double)outView->size() / (double)inView->size();
-    log()->get(LogLevel::Debug2)
-        << "Retaining " << outView->size() << " of " << inView->size()
-        << " points (" << 100 * frac << "%)\n";
+    // Iterate over immediate neighbors of current voxel, computing minimum
+    // distance between any already added point and the current point.
+    for (int xi = std::get<0>(v) - 1; xi < std::get<0>(v) + 2; ++xi)
+    {
+        for (int yi = std::get<1>(v) - 1; yi < std::get<1>(v) + 2; ++yi)
+        {
+            for (int zi = std::get<2>(v) - 1; zi < std::get<2>(v) + 2; ++zi)
+            {
+                Voxel candidate = std::make_tuple(xi, yi, zi);
+
+                // We have already visited the center voxel, and can skip it.
+                if (v == candidate)
+                    continue;
+
+                // Check that candidate voxel is occupied.
+                if (m_populatedVoxels.find(candidate) ==
+                    m_populatedVoxels.end())
+                    continue;
+
+                // Get list of coordinates in candidate voxel.
+                CoordList coords = m_populatedVoxels[candidate];
+                for (Coord const& coord : coords)
+                {
+                    // Compute Euclidean distance between current point and
+                    // candidate voxel.
+                    double xv = std::get<0>(coord);
+                    double yv = std::get<1>(coord);
+                    double zv = std::get<2>(coord);
+                    double dist =
+                        std::sqrt((xv - x) * (xv - x) + (yv - y) * (yv - y) +
+                                  (zv - z) * (zv - z));
+
+                    // If any point is closer than the minimum radius, we can
+                    // immediately return false, as the minimum distance
+                    // criterion is violated.
+                    if (dist < m_radius)
+                        return false;
+                }
+            }
+        }
+    }
 
-    viewSet.insert(outView);
-    return viewSet;
+    Coord coord = std::make_tuple(x, y, z);
+    if (m_populatedVoxels.find(v) != m_populatedVoxels.end())
+    {
+        m_populatedVoxels[v].push_back(coord);
+    }
+    else
+    {
+        CoordList coords;
+        coords.push_back(coord);
+        m_populatedVoxels.emplace(v, coords);
+    }
+    return true;
+}
+
+bool SampleFilter::processOne(PointRef& point)
+{
+    return voxelize(point);
 }
 
 } // namespace pdal

filters/SampleFilter.hpp

@@ -1,5 +1,5 @@
 /******************************************************************************
- * Copyright (c) 2016, Bradley J Chambers (brad.chambers@gmail.com)
+ * Copyright (c) 2016, 2020 Bradley J Chambers (brad.chambers@gmail.com)
  *
  * All rights reserved.
  *
@@ -35,32 +35,41 @@
 #pragma once
 
 #include <pdal/Filter.hpp>
-
-#include <string>
+#include <pdal/Streamable.hpp>
 
 namespace pdal
 {
 
-class Options;
-
-class PDAL_DLL SampleFilter : public pdal::Filter
+class PDAL_DLL SampleFilter : public Filter, public Streamable
 {
+    using Voxel = std::tuple<int, int, int>;
+    using Coord = std::tuple<double, double, double>;
+    using CoordList = std::vector<Coord>;
+
 public:
-    SampleFilter() : Filter()
-    {}
+    SampleFilter() : Filter() {}
     SampleFilter& operator=(const SampleFilter&) = delete;
     SampleFilter(const SampleFilter&) = delete;
 
     std::string getName() const;
 
 private:
+    double m_cell;
+    Arg* m_cellArg;
     double m_radius;
-    bool m_shuffle;
-    Arg* m_seedArg;
-    unsigned m_seed;
+    Arg* m_radiusArg;
+    double m_originX;
+    double m_originY;
+    double m_originZ;
+    std::map<Voxel, CoordList> m_populatedVoxels;
 
     virtual void addArgs(ProgramArgs& args);
+    virtual void prepared(PointTableRef table);
+    virtual bool processOne(PointRef& point);
+    virtual void ready(PointTableRef);
     virtual PointViewSet run(PointViewPtr view);
+
+    bool voxelize(PointRef& point);
 };
 
 } // namespace pdal