diff --git a/doc/stages/filters.sample.rst b/doc/stages/filters.sample.rst index 84150c160c..5b6a34d1a6 100755 --- a/doc/stages/filters.sample.rst +++ b/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 diff --git a/filters/SampleFilter.cpp b/filters/SampleFilter.cpp index 10903c7cb1..87bbbdf465 100644 --- a/filters/SampleFilter.cpp +++ b/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 #include -#include -#include -#include -#include #include -#include 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(Id::X); + double y = point.getFieldAs(Id::Y); + double z = point.getFieldAs(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 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 diff --git a/filters/SampleFilter.hpp b/filters/SampleFilter.hpp index aa298d97f5..56aee244b5 100644 --- a/filters/SampleFilter.hpp +++ b/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 - -#include +#include namespace pdal { -class Options; - -class PDAL_DLL SampleFilter : public pdal::Filter +class PDAL_DLL SampleFilter : public Filter, public Streamable { + using Voxel = std::tuple; + using Coord = std::tuple; + using CoordList = std::vector; + 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 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 diff --git a/filters/VoxelPoissonFilter.cpp b/filters/VoxelPoissonFilter.cpp deleted file mode 100644 index 19a715ff04..0000000000 --- a/filters/VoxelPoissonFilter.cpp +++ /dev/null @@ -1,171 +0,0 @@ -/****************************************************************************** - * Copyright (c) 2020, 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 "VoxelPoissonFilter.hpp" - -namespace pdal -{ - -static StaticPluginInfo const s_info{ - "filters.voxelpoisson", "Voxelized Poisson Sampling", - "http://pdal.io/stages/filters.voxelpoisson.html"}; - -CREATE_STATIC_STAGE(VoxelPoissonFilter, s_info) - -VoxelPoissonFilter::VoxelPoissonFilter() {} - -std::string VoxelPoissonFilter::getName() const -{ - return s_info.name; -} - -void VoxelPoissonFilter::addArgs(ProgramArgs& args) -{ - args.add("cell", "Cell size", m_cell, 1.0); - // consider specifying radius instead, relationship being that the cell - // size should be radius/sqrt(2), or radius should be sqrt(2)*cell -} - -void VoxelPoissonFilter::ready(PointTableRef) -{ - m_populatedVoxels.clear(); - m_radius = std::sqrt(2.0) * m_cell * 0.5; - log()->get(LogLevel::Debug) - << "cell " << m_cell << ", radius " << m_radius << std::endl; -} - -PointViewSet VoxelPoissonFilter::run(PointViewPtr view) -{ - PointViewPtr output = view->makeNew(); - PointRef point(*view); - for (PointId id = 0; id < view->size(); ++id) - { - point.setPointId(id); - if (voxelize(point)) - output->appendPoint(*view, id); - } - - PointViewSet viewSet; - viewSet.insert(output); - return viewSet; -} - -bool VoxelPoissonFilter::voxelize(PointRef& point) -{ - /* - * Calculate the voxel coordinates for the incoming point. - * gx, gy, gz will be the global coordinates from (0, 0, 0). - */ - double x = point.getFieldAs(Dimension::Id::X); - double y = point.getFieldAs(Dimension::Id::Y); - double z = point.getFieldAs(Dimension::Id::Z); - if (m_populatedVoxels.empty()) - { - m_originX = x - (m_cell / 2); - m_originY = y - (m_cell / 2); - m_originZ = z - (m_cell / 2); - } - - // Offset by origin. - x -= m_originX; - y -= m_originY; - z -= m_originZ; - - Voxel v = std::make_tuple((int)(std::floor(x / m_cell)), - (int)(std::floor(y / m_cell)), - (int)(std::floor(z / m_cell))); - - // this is silly, works for now - double xd = point.getFieldAs(Dimension::Id::X); - double yd = point.getFieldAs(Dimension::Id::Y); - double zd = point.getFieldAs(Dimension::Id::Z); - - // get all points in m_populateVoxels +/- 1 voxel away from v in each - // dimension compute minimum distance to any of these points (there is - // probably an optimization we can make here), and only insert point if - // greater than radius away - - double mindist = std::numeric_limits::max(); - 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) - { - auto it = m_populatedVoxels.find(std::make_tuple(xi, yi, zi)); - if (it == m_populatedVoxels.end()) - continue; - std::vector> coords = - m_populatedVoxels[std::make_tuple(xi, yi, zi)]; - for (std::tuple const& coord : coords) - { - double xv = std::get<0>(coord); - double yv = std::get<1>(coord); - double zv = std::get<2>(coord); - double dist = std::sqrt((xv - xd) * (xv - xd) + - (yv - yd) * (yv - yd) + - (zv - zd) * (zv - zd)); - if (dist < mindist) - mindist = dist; - } - } - } - } - - if (mindist > m_radius) - { - std::tuple vd = std::make_tuple(xd, yd, zd); - auto it = m_populatedVoxels.find(v); - if (it != m_populatedVoxels.end()) - { - m_populatedVoxels[v].push_back(vd); - } - else - { - std::vector> coord; - coord.push_back(vd); - m_populatedVoxels.emplace(v, coord); - } - return true; - } - - return false; -} - -bool VoxelPoissonFilter::processOne(PointRef& point) -{ - return voxelize(point); -} - -} // namespace pdal diff --git a/filters/VoxelPoissonFilter.hpp b/filters/VoxelPoissonFilter.hpp deleted file mode 100644 index 5db47f432b..0000000000 --- a/filters/VoxelPoissonFilter.hpp +++ /dev/null @@ -1,74 +0,0 @@ -/****************************************************************************** - * Copyright (c) 2020, 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 - -namespace pdal -{ - -class PointLayout; -class PointView; - -class PDAL_DLL VoxelPoissonFilter : public Filter, public Streamable -{ - using Voxel = std::tuple; - -public: - VoxelPoissonFilter(); - VoxelPoissonFilter& operator=(const VoxelPoissonFilter&) = delete; - VoxelPoissonFilter(const VoxelPoissonFilter&) = delete; - - std::string getName() const override; - -private: - virtual void addArgs(ProgramArgs& args) override; - virtual PointViewSet run(PointViewPtr view) override; - virtual void ready(PointTableRef) override; - virtual bool processOne(PointRef& point) override; - - bool voxelize(PointRef& point); - - double m_cell; - double m_radius; - double m_originX; - double m_originY; - double m_originZ; - std::map>> - m_populatedVoxels; -}; - -} // namespace pdal