Merge pull request #8554 from afabri/PMP-PoissonDiskSampling-GF

Add Weighted Sample Elimination
This commit is contained in:
Sébastien Loriot 2025-03-28 18:03:36 +01:00
commit f3ff704eea
18 changed files with 748 additions and 12 deletions

View File

@ -19,6 +19,7 @@
#include <CGAL/Default.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/property_map.h>
#include <CGAL/Dimension.h>
#include <boost/mpl/has_xxx.hpp>
@ -356,6 +357,22 @@ struct GetPolygonSoupGeomTraits
> ::type type;
};
namespace internal {
template <class Point, int d = CGAL::Ambient_dimension<Point>::value>
struct Vector_matching_point {
typedef typename Kernel_traits<Point>::Kernel::Vector_d type;
};
template <class Point>
struct Vector_matching_point<Point, 2> {
typedef typename Kernel_traits<Point>::Kernel::Vector_2 type;
};
template <class Point>
struct Vector_matching_point<Point, 3> {
typedef typename Kernel_traits<Point>::Kernel::Vector_3 type;
};
}
template <class PointRange, class NamedParameters, class PointMap = Default, class NormalMap = Default>
struct Point_set_processing_3_np_helper
@ -380,7 +397,7 @@ struct Point_set_processing_3_np_helper
typedef typename Geom_traits::FT FT; // public
typedef Constant_property_map<Value_type, typename Geom_traits::Vector_3> DummyNormalMap;
typedef Constant_property_map<Value_type, typename internal::Vector_matching_point<Point>::type> DummyNormalMap;
typedef typename Default::Get<NormalMap, DummyNormalMap>::type DefaultNMap;
typedef typename internal_np::Lookup_named_param_def<

View File

@ -3326,6 +3326,18 @@ pages = "207--221"
year={1997}
}
@article{cgal:y-sefpdss,
author = {Yuksel, Cem},
title = {Sample Elimination for Generating Poisson Disk Sample Sets},
journal = {Computer Graphics Forum},
volume = {34},
number = {2},
pages = {25-32},
doi = {https://doi.org/10.1111/cgf.12538},
url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/cgf.12538},
year = {2015}
}
% ----------------------------------------------------------------------------
% END OF BIBFILE
% ----------------------------------------------------------------------------

View File

@ -127,7 +127,8 @@ public:
bool
operator==(const Self& x) const {
return (var == x.var) && (index == x.index);
CGAL_kernel_assertion(var == x.var);
return index == x.index;
}
bool

View File

@ -127,7 +127,8 @@ public:
bool
operator==(const Self& x) const {
return (var == x.var) && (index == x.index);
CGAL_kernel_assertion(var == x.var);
return (index == x.index);
}
bool

View File

@ -12,8 +12,11 @@
- Added a function in the [visitor of the corefinement based methods](https://doc.cgal.org/6.1/Polygon_mesh_processing/classPMPCorefinementVisitor.html)
to know faces in the output meshes that are corresponding to input coplanar faces.
### [Algebraic Kernel](https://doc.cgal.org/6.1/Manual/packages.html#PkgAlgebraicKernelD)
### [Point Set Processing](https://doc.cgal.org/6.1/Manual/packages.html#PkgPointSetProcessing3)
- Added `poisson_eliminate()` to downsample a point cloud to a target size while providing Poisson disk property, i.e., a larger minimal distance between points.
### [Algebraic Kernel](https://doc.cgal.org/6.1/Manual/packages.html#PkgAlgebraicKernelD)
- **Breaking change**: Classes based on the RS Library are no longer provided.
### [BGL](https://doc.cgal.org/6.1/Manual/packages.html#PkgBGL)

View File

@ -32,7 +32,7 @@ format.
\cgalPkgDescriptionBegin{Point Set Processing,PkgPointSetProcessing3}
\cgalPkgPicture{point_set_processing_detail.png}
\cgalPkgSummaryBegin
\cgalPkgAuthors{Pierre Alliez, Simon Giraudot, Clément Jamin, Florent Lafarge, Quentin Mérigot, Jocelyn Meyron, Laurent Saboret, Nader Salman, Shihao Wu, and Necip Fazil Yildiran}
\cgalPkgAuthors{Pierre Alliez, Simon Giraudot, Clément Jamin, Florent Lafarge, Quentin Mérigot, Jocelyn Meyron, Sven Oesau, Laurent Saboret, Nader Salman, Shihao Wu, and Necip Fazil Yildiran}
\cgalPkgDesc{This \cgal component implements methods to analyze and process unorganized point sets. The input is an unorganized point set, possibly with normal attributes (unoriented or oriented). The point set can be analyzed to measure its average spacing, and processed through functions devoted to the simplification, outlier removal, smoothing, normal estimation, normal orientation, feature edges estimation and registration.}
\cgalPkgManuals{Chapter_Point_Set_Processing,PkgPointSetProcessing3Ref}
\cgalPkgSummaryEnd
@ -75,6 +75,7 @@ format.
- `CGAL::vcm_estimate_normals()`
- `CGAL::vcm_is_on_feature_edge()`
- `CGAL::structure_point_set()`
- `CGAL::poisson_eliminate()`
\cgalCRPSection{I/O (All Formats)}

View File

@ -650,7 +650,7 @@ functions in this component.)
\section Point_set_processing_3Simplification Simplification
Four simplification functions are devised to reduce an input point set.
Five simplification functions are devised to reduce an input point set.
Function `random_simplify_point_set()` randomly deletes a
user-specified fraction of points from the input point set. This
@ -674,6 +674,10 @@ Function `wlop_simplify_and_regularize_point_set()` not only simplifies,
but also regularizes downsampled points. This is an implementation of
the Weighted Locally Optimal Projection (WLOP) algorithm \cgalCite{wlop-2009}.
Function `poisson_eliminate()` is a greedy down sampling method that generates
a subset of the input points with Poisson disk property. This is an implementation of the
<em>Sample Elimination for Generating Poisson Disk Sample Sets</em> method \cgalCite{cgal:y-sefpdss}.
\subsection Point_set_processing_3Example_grid_simplification Grid Simplification Example
@ -749,12 +753,49 @@ A parallel version of WLOP is provided and requires the executable to be linked
<a href="https://github.com/oneapi-src/oneTBB">Intel TBB library</a>.
To control the number of threads used, the user may use the tbb::task_scheduler_init class.
See the <a href="https://software.intel.com/content/www/us/en/develop/documentation/onetbb-documentation/top.html">TBB documentation</a>
for more details. We provide below a speed-up chart generated using the parallel version of the WLOP algorithm. The machine used is a PC running Windows 7 64-bits with a 4-core i7-4700HQ@2.40GHz CPU with 8GB of RAM.
for more details. We provide below a speed-up chart generated using the parallel version of the WLOP algorithm. The machine used is a PC running Windows 7 64-bits with a 4-core i7-4700HQ\@2.40GHz CPU with 8GB of RAM.
\cgalFigureBegin{Point_set_processing_3figWLOP_parallel_performance, parallel_WLOP_performance.jpg}
Parallel WLOP speed-up, compared to the sequential version of the algorithm.
\cgalFigureEnd
\subsection Point_set_processing_3PoissonElimination Poisson Sample Elimination
The Poisson sample elimination is a greedy downsampling method that calculates a weight for each input point based on the density of its local neighborhood. Subsequently, the point with the highest weight is removed and the weights of the remaining points around it are updated until the target size is reached. A custom function to calculate the weight of a point can be provided.
The Poisson sample elimination has the following parameters:
- *dimension*: The dimension of the sampling domain of the point cloud, e.g., 2 if the point cloud is sampled from a surface, while the ambient dimension is typically 3. The default value is 2.
- *weight_function*: A custom function can be provided to calculate the weight between two points. The type of the functor is as follows:
\code{.cpp}
double(*func)(const Point &p, const Point &n, double squared_distance, double r_max)
\endcode
The default weight is \f$\left(1 - \frac{d_{p,n}}{2r_{max}}\right)^8\f$ with \f$d_{p,n}\f$ being the distance between the point p and its neighbor n.
- \f$r_{max}\f$: The \f$r_{max}\f$ parameter specifies the radius of the neighborhood, i.e., the neighboring points that are used to calculate the weight of a point. \f$r_{max}\f$ has to be provided if a custom *weight_function* is used. Only points within a distance of \f$r_{max}\f$ are used to calculate the weight of a point. A large value can thus cause a large running time. The default is calculated based in the bounding volume \f$V\f$ of the input points, the *dimension* parameter and the number of input points \f$N\f$:
\f$\mathrm{dimension} = 2:\qquad\f$ \f$r_{max} = \sqrt{\frac{V}{2\sqrt{3}N}}\f$
\f$\mathrm{dimension} = 3:\qquad\f$ \f$r_{max} = \sqrt{\frac{V}{4\sqrt{2}N}}\f$
- *progressive*: The output points of the function can be reordered to be progressive. A progressive ordering will increase the running time by a factor of at most 2 as the function is internally applied several times on increasingly smaller subsets. The default value is false.
\cgalFigureAnchor{Point_set_processing_3figPoisson_elimination}
<center>
<img src="poisson_elimination.png" style="max-width:80%;"/>
</center>
\cgalFigureCaptionBegin{Point_set_processing_3figPoisson_elimination}
Poisson elimination on point cloud randomly sampled from a mesh with 20k points showing results with \f$\frac{1}{2}\f$, \f$\frac{1}{4}\f$ and \f$\frac{1}{8}\f$ of the input size.
\cgalFigureCaptionEnd
\subsubsection Point_set_processing_3Example_poisson_elimination Poisson Sample Elimination Example
The following example reads a point cloud, applies Poisson elimination to reduce the point cloud to \f$\frac{1}{5}\f$ of its size and saves it into a file.
\cgalExample{Point_set_processing_3/poisson_eliminate_example.cpp}
\subsubsection Point_set_processing_3Example_poisson_elimination_on_mesh Poisson Sample Elimination from Mesh Example
The following example first samples a point cloud from a surface mesh and then uses Poisson elimination to provide an evenly sampled point cloud.
\cgalExample{Point_set_processing_3/poisson_eliminate_from_mesh_example.cpp}
\section Point_set_processing_3Smoothing Smoothing
@ -787,13 +828,14 @@ The following example reads a set of points with normals and smooths them via bi
Comparison for two smoothing methods: Left: Input, 250K points, normal-color mapping. Middle: Jet smoothing result, 197 seconds. Right: Bilateral smoothing result, 110 seconds.
\cgalFigureEnd
\subsubsection Point_set_processing_3Bilateral_smoothing_parallel_performance Parallel
\subsubsection Point_set_processing_3Bilateral_smoothing_parallel_performance Parallel Performance
Performance:
A parallel version of bilateral smoothing is provided and requires the executable to be linked against the
<a href="https://github.com/oneapi-src/oneTBB">Intel TBB library</a>.
The number of threads used is controlled through the tbb::task_scheduler_init class.
See the <a href="https://software.intel.com/content/www/us/en/develop/documentation/onetbb-documentation/top.html">TBB documentation</a> for more details. We provide below a speed-up chart generated using the parallel version of the bilateral smoothing algorithm. The machine used is a PC running Windows 7 64-bits with a 4-core i7-4700HQ@2.40GHz CPU with 8GB of RAM.
See the <a href="https://software.intel.com/content/www/us/en/develop/documentation/onetbb-documentation/top.html">TBB documentation</a> for more details.
We provide below a speed-up chart generated using the parallel version of the bilateral smoothing algorithm.
The machine used is a PC running Windows 7 64-bits with a 4-core i7-4700HQ\@2.40GHz CPU with 8GB of RAM.
\cgalFigureBegin{Point_set_processing_3Bilateral_smoothing_parallel_performance, parallel_bilateral_smooth_point_set_performance.jpg}
Parallel bilateral smoothing speed-up, compared to the sequential version of the algorithm.
@ -1009,7 +1051,7 @@ as well as the normal and feature edge estimation functions based on it.
Florent Lafarge with the help of Simon Giraudot contributed the point set structuring algorithm.
Started from GSoC'2019, Necip Fazil Yildiran with the help of Nicolas Mellado and Simon Giraudot introduced the wrappers for OpenGR and PointMatcher
libraries that perform registration on two point sets.
Poisson sample elimination is a reimplementation by Sven Oesau based on \cgalCite{cgal:y-sefpdss}. The original source code <a href="https://github.com/cemyuksel/cyCodeBase/">cyCodeBase</a> was used as a reference.
*/
} /* namespace CGAL */

View File

@ -5,6 +5,7 @@ Algebraic_foundations
Circulator
Stream_support
Poisson_surface_reconstruction_3
Polygon_mesh_processing
Property_map
Bounding_volumes
Principal_component_analysis

View File

@ -23,4 +23,6 @@
\example Point_set_processing_3/edges_example.cpp
\example Point_set_processing_3/structuring_example.cpp
\example Point_set_processing_3/callback_example.cpp
\example Point_set_processing_3/poisson_eliminate_example.cpp
\example Point_set_processing_3/poisson_eliminate_from_mesh_example.cpp
*/

Binary file not shown.

After

Width:  |  Height:  |  Size: 100 KiB

Binary file not shown.

Before

Width:  |  Height:  |  Size: 68 KiB

After

Width:  |  Height:  |  Size: 38 KiB

View File

@ -42,7 +42,9 @@ foreach(
edge_aware_upsample_point_set_example
structuring_example
read_ply_points_with_colors_example
write_ply_points_example)
write_ply_points_example
poisson_eliminate_example
poisson_eliminate_from_mesh_example)
create_single_source_cgal_program("${target}.cpp")
target_link_libraries(${target} PRIVATE CGAL::CGAL)
endforeach()

View File

@ -0,0 +1,46 @@
#include <vector>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/poisson_eliminate.h>
#include <CGAL/IO/write_points.h>
#include <CGAL/IO/read_points.h>
typedef CGAL::Simple_cartesian<double> K;
typedef K::Point_3 Point_3;
void sampling(const std::string& filename, double size_factor = 0.2) {
if (size_factor >= 1.0) {
std::cout << "usage poisson_eliminate_example filename size_factor" << std::endl
<< "0 < size_factor < 1" << std::endl;
return;
}
std::vector<Point_3> points;
if (!CGAL::IO::read_points(
filename,
std::back_inserter(points))) {
std::cerr << "Error: cannot read file!" << std::endl;
return;
}
std::size_t target_size = std::size_t(points.size() * size_factor);
std::vector<Point_3> output;
output.reserve(target_size);
CGAL::poisson_eliminate(points, target_size, std::back_inserter(output));
CGAL::IO::write_points("out.xyz", output, CGAL::parameters::stream_precision(17));
}
int main(int argc, char* argv[])
{
if (argc < 2)
sampling(CGAL::data_file_path("points_3/radar.xyz"));
else if (argc < 3)
sampling(argv[1]);
else if (argc < 4)
sampling(argv[1], std::atof(argv[2]));
return 0;
}

View File

@ -0,0 +1,62 @@
#include <vector>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/poisson_eliminate.h>
#include <CGAL/Polygon_mesh_processing/distance.h>
#include <CGAL/IO/write_points.h>
using Point_3 = CGAL::Simple_cartesian<double>::Point_3;
using Mesh = CGAL::Surface_mesh<Point_3>;
int main(int argc, char* argv[])
{
std::vector<Point_3> points;
std::size_t sampling_size = 20000;
std::size_t target_size = 10000;
Mesh mesh;
if (argc == 4) {
target_size = std::atoi(argv[2]);
sampling_size = std::atoi(argv[3]);
if (sampling_size < target_size) {
std::cout << "usage: poisson_eliminate_from_mesh_example input_mesh target_size sampling_size"
<< std::endl << "with sampling_size > target_size, e.g., by a factor of 2" << std::endl;
return -1;
}
}
if (argc == 3) {
target_size = std::atoi(argv[2]);
sampling_size = 2 * target_size;
}
if (argc < 2)
CGAL::IO::read_polygon_mesh(CGAL::data_file_path("meshes/elephant.off"), mesh);
else
CGAL::IO::read_polygon_mesh(argv[1], mesh);
if (mesh.number_of_faces() == 0) {
std::cout << "mesh could not be loaded" << std::endl;
std::cout << "usage: poisson_eliminate_from_mesh_example input_mesh target_size sampling_size"
<< std::endl << "with sampling_size > target_size, e.g., by a factor of 2" << std::endl;
return -1;
}
points.reserve(sampling_size);
CGAL::Polygon_mesh_processing::sample_triangle_mesh(mesh,
std::back_inserter(points),
CGAL::parameters::number_of_points_on_faces(sampling_size)
.do_sample_vertices(false)
.do_sample_edges(false));
CGAL::IO::write_points("out_sampled.xyz", points, CGAL::parameters::stream_precision(17));
std::vector<Point_3> output;
output.reserve(target_size);
CGAL::poisson_eliminate(points, target_size, std::back_inserter(output));
CGAL::IO::write_points("out_poisson_eliminated.xyz", output, CGAL::parameters::stream_precision(17));
return 0;
}

View File

@ -0,0 +1,452 @@
// Copyright (c) 2025 GeometryFactory (France).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Sven Oesau
#ifndef CGAL_POISSON_ELIMINATE_H
#define CGAL_POISSON_ELIMINATE_H
#include <CGAL/license/Point_set_processing_3.h>
#include <CGAL/Kd_tree.h>
#include <CGAL/Splitters.h>
#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/Search_traits_2.h>
#include <CGAL/Search_traits_3.h>
#include <CGAL/Search_traits_d.h>
#include <CGAL/Search_traits_adapter.h>
#include <CGAL/Named_function_parameters.h>
#include <CGAL/boost/graph/named_params_helper.h>
namespace CGAL {
namespace poisson_eliminate_impl {
double get_maximum_radius(std::size_t dimension, std::size_t sample_size, double domain_size) {
double sampleArea = domain_size / double(sample_size);
double r_max;
switch (dimension) {
case 2: r_max = CGAL::sqrt(sampleArea / (2.0 * CGAL::sqrt(3.0))); break;
case 3: r_max = std::pow(sampleArea / (4.0 * CGAL::sqrt(2.0)), 1.0 / 3.0); break;
default:
double c;
std::size_t d_start;
if ((dimension & 1)) {
c = 2.0;
d_start = 3;
}
else {
c = CGAL_PI;
d_start = 4;
}
for (std::size_t d = d_start; d <= dimension; d += 2)
c *= 2.0 * CGAL_PI / double(d);
r_max = std::pow(sampleArea / c, 1.0 / double(dimension));
break;
}
return r_max;
}
double get_minimum_radius(std::size_t input_size, std::size_t output_size, double beta, double gamma, double r_max) {
double ratio = output_size / double(input_size);
return r_max * (1 - std::pow(ratio, gamma)) * beta;
}
template<class Point, class FT = typename CGAL::Kernel_traits<Point>::FT, unsigned int dim = CGAL::Ambient_dimension<Point>::value>
Point construct(const std::array<FT, dim>& coords) {
if constexpr (dim == 2) { return Point(coords[0], coords[1]); }
else if constexpr (dim == 3) { return Point(coords[0], coords[1], coords[2]); }
else { return Point(coords); }
}
template<class Point, class FT = typename CGAL::Kernel_traits<Point>::FT>
void copy_and_replace(const Point& in, Point& out, std::size_t replace, FT value) {
std::array<FT, CGAL::Ambient_dimension<Point>::value> tmp;
for (std::size_t i = 0; i < CGAL::Ambient_dimension<Point>::value; i++)
if (i != replace)
tmp[i] = in[int(i)];
else
tmp[replace] = value;
out = construct<Point, FT, CGAL::Ambient_dimension<Point>::value>(tmp);
}
template<class PointRange, class PointMap>
class Indexed_extended_point_map {
public:
using value_type = typename boost::property_traits<PointMap>::value_type;
using reference = std::conditional_t< std::is_reference_v<typename boost::property_traits<PointMap>::reference>,
const value_type&, value_type>;
using key_type = std::size_t;
using category = boost::readable_property_map_tag;
private:
const PointRange& range;
const std::vector<value_type>& tiling_points;
const PointMap& point_map;
public:
Indexed_extended_point_map(const PointRange& input, const std::vector<value_type>& tiling, const PointMap& point_map)
: range(input), tiling_points(tiling), point_map(point_map) {}
reference operator[](key_type k) const {
CGAL_assertion(k < (range.size() + tiling_points.size()));
if (range.size() > k)
return get(point_map, *(range.begin() + k));
else
return tiling_points.at(k - range.size());
}
friend reference get(const Indexed_extended_point_map& ppmap, key_type k) {
CGAL_assertion(k < (ppmap.range.size() + ppmap.tiling_points.size()));
if ((k < ppmap.range.size()))
return get(ppmap.point_map, *(ppmap.range.begin() + k));
else
return ppmap.tiling_points.at(k - ppmap.range.size());
}
};
template<class Point>
class Weight_functor {
public:
Weight_functor(double r_min = 0, double alpha = 8) : r_min(CGAL::to_double(r_min)), alpha(CGAL::to_double(alpha)) {}
double operator()(const Point &, const Point &, double d2, double r_max) {
double d = CGAL::sqrt(d2);
if (d < r_min) d = r_min;
return std::pow(double(1) - d / r_max, alpha);
}
private:
double r_min;
double alpha;
};
void move_down(std::vector<std::size_t>& heap, std::vector<std::size_t>& heap_pos, std::size_t heap_size, std::vector<double>& weights, std::size_t idx) {
CGAL_assertion(idx <= heap.size());
std::size_t child = idx * 2 + 1;
while (child + 1 < heap_size) {
if (weights[heap[child]] < weights[heap[child + 1]]) child++;
if (weights[heap[idx]] >= weights[heap[child]])
break;
std::swap(heap[idx], heap[child]);
heap_pos[heap[idx]] = idx;
heap_pos[heap[child]] = child;
idx = child;
child = idx * 2 + 1;
}
if (child < heap_size) {
if (weights[heap[idx]] < weights[heap[child]]) {
std::swap(heap[idx], heap[child]);
heap_pos[heap[idx]] = idx;
heap_pos[heap[child]] = child;
}
}
}
void pop_heap(std::vector<std::size_t>& heap, std::vector<std::size_t>& heap_pos, std::size_t &heap_size, std::vector<double>& weights) {
std::swap(heap.front(), heap[heap_size - 1]);
heap_pos[heap.front()] = 0;
heap_pos[heap[heap_size - 1]] = heap_size - 1;
heap_size--;
move_down(heap, heap_pos, heap_size, weights, 0);
}
template<class Point, int d = CGAL::Ambient_dimension<Point>::value>
struct Compute_squared_distance {
using GeomTraits = typename CGAL::Kernel_traits<Point>::Kernel;
using FT = typename GeomTraits::FT;
FT operator()(const Point& a, const Point& b) {
return GeomTraits().squared_distance_d_object()(a, b);
}
};
template<class Point>
struct Compute_squared_distance<Point, 2> {
using GeomTraits = typename CGAL::Kernel_traits<Point>::Kernel;
using FT = typename GeomTraits::FT;
FT operator()(const Point& a, const Point& b) {
return GeomTraits().compute_squared_distance_2_object()(a, b);
}
};
template<class Point>
struct Compute_squared_distance<Point, 3> {
using GeomTraits = typename CGAL::Kernel_traits<Point>::Kernel;
using FT = typename GeomTraits::FT;
FT operator()(const Point& a, const Point& b) {
return GeomTraits().compute_squared_distance_3_object()(a, b);
}
};
template<class Point, int d = CGAL::Ambient_dimension<Point>::value>
struct Search_traits : public CGAL::Search_traits_d<typename CGAL::Kernel_traits<Point>::Kernel>{
};
template<class Point>
struct Search_traits<Point, 2> : public CGAL::Search_traits_2<typename CGAL::Kernel_traits<Point>::Kernel> {
};
template<class Point>
struct Search_traits<Point, 3> : public CGAL::Search_traits_3<typename CGAL::Kernel_traits<Point>::Kernel> {
};
}
/**
\ingroup PkgPointSetProcessing3Algorithms
performs Poisson disk elimination with a desired output size. A greedy method that calculates a weight based on the
neighborhood of each point and eliminates points until the output size is reached.
For more details, please refer to \cgalCite{cgal:y-sefpdss}.
\tparam PointRange is a model of `RandomAccessRange`. The value type of
its iterator is the key type of the named parameter `point_map`.
\tparam OutputIterator Type of the output iterator. Must accept input of the same type as the iterator of `PointRange`.
\param points input point range
\param number_of_points target output size, must be smaller that the number of points in the input range
\param output where output points are put
\param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
\cgalNamedParamsBegin
\cgalParamNBegin{point_map}
\cgalParamDescription{a property map associating points to the elements of the point set `points`}
\cgalParamType{a model of `ReadablePropertyMap` whose key type is the value type
of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
\cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
\cgalParamNEnd
\cgalParamNBegin{dimension}
\cgalParamDescription{The sampling domain of `points`, e.g., 2 if the points have been sampled from a 2d surface.}
\cgalParamType{unsigned integer}
\cgalParamDefault{2}
\cgalParamNEnd
\cgalParamNBegin{progressive}
\cgalParamDescription{reorders the points in `output` in a progressive way, i.e., the first n points in `output` with n < `number_of_points` have a Poisson disk distribution with a larger radius. }
\cgalParamType{Boolean}
\cgalParamDefault{`false`}
\cgalParamNEnd
\cgalParamNBegin{maximum_radius}
\cgalParamDescription{radius of the Poisson disk in which the neighboring points are taken into account for elimination.}
\cgalParamType{double}
\cgalParamDefault{the default is calculated from the `dimension`, the volume of the bounding box and the `number_of_points`. For more details, see parameter section in \ref Point_set_processing_3PoissonElimination.}
\cgalParamNEnd
\cgalParamNBegin{weight_function}
\cgalParamDescription{a weight function that calculates the weight of a neighbor point based on its squared distance and the `maximum_radius`.}
\cgalParamType{an instance of `std::function<double(const Point&, const Point&, double,double)>`.}
\cgalParamDefault{See parameter section in \ref Point_set_processing_3PoissonElimination.}
\cgalParamNEnd
\cgalParamNBegin{geom_traits}
\cgalParamDescription{an instance of a geometric traits class}
\cgalParamType{a model of `Kernel`}
\cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
\cgalParamNEnd
\cgalNamedParamsEnd
*/
template<class PointRange, class OutputIterator, class NamedParameters = parameters::Default_named_parameters>
void poisson_eliminate(const PointRange &points, std::size_t number_of_points, OutputIterator output, const NamedParameters& np = parameters::default_values()) {
using parameters::choose_parameter;
using parameters::get_parameter;
if (points.size() == 0)
return;
using NP_helper = Point_set_processing_3_np_helper<PointRange, NamedParameters>;
using PointMap = typename NP_helper::Const_point_map;
using Point = typename boost::property_traits<PointMap>::value_type;
using GeomTraits = typename NP_helper::Geom_traits;
using FT = typename GeomTraits::FT;
using IPM = poisson_eliminate_impl::Indexed_extended_point_map<PointRange, PointMap>;
PointMap point_map = NP_helper::get_const_point_map(points, np);
const unsigned int ambient_dimension = CGAL::Ambient_dimension<Point>::value;
using Search_base = poisson_eliminate_impl::Search_traits<Point>;
using Search_traits = CGAL::Search_traits_adapter<std::size_t, IPM, Search_base>;
using Splitter = CGAL::Sliding_midpoint<Search_traits>;
using Fuzzy_sphere = CGAL::Fuzzy_sphere<Search_traits>;
using Tree = CGAL::Kd_tree<Search_traits, Splitter, CGAL::Tag_true, CGAL::Tag_true>;
poisson_eliminate_impl::Compute_squared_distance<Point> squared_distance;
// hard coded parameters alpha, beta and gamma with values proposed in publication
const double alpha = 8;
const double beta = 0.65; // not used currently
const double gamma = 1.5; // not used currently
// named parameters for weight limiting, currently not used as the performance gain seems small and there is a high risk to get a defective output
const bool weight_limiting = parameters::choose_parameter(parameters::get_parameter(np, internal_np::weight_limiting), false);
// named parameter for progressive
const bool progressive = parameters::choose_parameter(parameters::get_parameter(np, internal_np::progressive), false);
// named parameter for tiling
const bool tiling = parameters::choose_parameter(parameters::get_parameter(np, internal_np::tiling), false);
const unsigned int dimension = parameters::choose_parameter(parameters::get_parameter(np, internal_np::dimension), 2);
CGAL_assertion(ambient_dimension >= dimension);
// Needs a multi-dimension solution
std::array<FT, ambient_dimension> lower, upper;
const Point &first = get(point_map, *(points.begin()));
for (unsigned int i = 0; i < ambient_dimension; i++)
lower[i] = upper[i] = first[i];
for (const auto &pt : points) {
const Point& p = get(point_map, pt);
for (unsigned int i = 0; i < ambient_dimension; i++) {
lower[i] = (p[i] < lower[i]) ? p[i] : lower[i];
upper[i] = (p[i] > upper[i]) ? p[i] : upper[i];
}
}
double domain_size = 1;
for (std::size_t i = 0; i < ambient_dimension; i++)
domain_size *= CGAL::to_double(upper[i] - lower[i]);
// named parameter for r_max
double r_max = CGAL::to_double(parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_radius), 2 * poisson_eliminate_impl::get_maximum_radius(dimension, number_of_points, domain_size)));
double r_min = CGAL::to_double(weight_limiting ? poisson_eliminate_impl::get_minimum_radius(points.size(), number_of_points, beta, gamma, r_max) : 0);
auto weight_functor = parameters::choose_parameter(parameters::get_parameter(np, internal_np::weight_function), poisson_eliminate_impl::Weight_functor<Point>(r_min, alpha));
std::size_t heap_size = points.size();
std::vector<Point> tiling_points;
IPM ipm(points, tiling_points, point_map);
auto tile_point = [&tiling_points, &lower, &upper, &r_max](const Point& p, std::size_t dim = 0) {
auto do_tiling = [&tiling_points, &lower, &upper, &r_max](const auto& self, const Point& p, std::size_t dim) -> void {
auto it = p.cartesian_begin();
if (lower[int(dim)] > (*(it + dim) - r_max)) {
FT v = 2 * lower[int(dim)] - (*(it + dim));
Point p2;
poisson_eliminate_impl::copy_and_replace(p, p2, dim, v);
tiling_points.emplace_back(p2);
if (dim + 1 < CGAL::Ambient_dimension<Point>::value)
self(self, tiling_points.back(), dim + 1);
}
if (upper[int(dim)] < (*(it + dim) + r_max)) {
FT v = 2 * upper[int(dim)]- (*(it + dim));
Point p2;
poisson_eliminate_impl::copy_and_replace(p, p2, dim, v);
tiling_points.emplace_back(p2);
if (dim + 1 < CGAL::Ambient_dimension<Point>::value)
self(self, tiling_points.back(), dim + 1);
}
if (dim + 1 < CGAL::Ambient_dimension<Point>::value)
self(self, p, dim + 1);
};
do_tiling(do_tiling, p, dim);
};
if (tiling) {
for (std::size_t i = 0;i<points.size();i++)
tile_point(get(point_map, *(points.begin() + i)));
}
// Tiling requires either to copy the point range or to use a special Index_property_map with a second array of points to keep the input range const.
Tree tree(boost::counting_iterator<std::size_t>(0),
boost::counting_iterator<std::size_t>(points.size() + tiling_points.size()),
Splitter(), Search_traits(ipm));
tree.build();
std::vector<std::size_t> heap(boost::counting_iterator<std::size_t>(0), boost::counting_iterator<std::size_t>(points.size()));
std::vector<std::size_t> heap_pos(points.size());
std::size_t target_points = number_of_points;
do {
// Computing weights
std::vector<double> weights(points.size(), 0);
std::vector<std::size_t> res;
for (std::size_t i = 0; i < heap_size; ++i) {
const Point& p = get(ipm, heap[i]);
res.clear();
tree.search(std::back_inserter(res), Fuzzy_sphere(p, r_max, 0, Search_traits(ipm)));
for (std::size_t n = 0; n < res.size(); n++) {
if (i == res[n] || res[n] >= heap_pos.size() || heap_pos[res[n]] >= heap_size)
continue;
const Point &p2 = get(ipm, res[n]);
double d2 = CGAL::to_double(squared_distance(p, p2));
weights[i] += weight_functor(p, p2, d2, r_max);
}
}
auto weight_order = [&](std::size_t a, std::size_t b) {
return weights[a] < weights[b];
};
std::make_heap(heap.begin(), heap.begin() + heap_size, weight_order);
// inverse mapping
for (std::size_t i = 0; i < heap.size(); i++)
heap_pos[heap[i]] = i;
while (heap_size > target_points) {
std::size_t i = heap[0];
poisson_eliminate_impl::pop_heap(heap, heap_pos, heap_size, weights);
// Adjust weights of neighbors
res.clear();
const Point& p = get(point_map, *(points.begin() + i));
tree.search(std::back_inserter(res), Fuzzy_sphere(p, r_max, 0, Search_traits(ipm)));
for (std::size_t n = 0; n < res.size(); n++) {
if (i == res[n] || res[n] >= points.size() || heap_pos[res[n]] >= heap_size)
continue;
const Point &p2 = get(point_map, *(points.begin() + res[n]));
double d2 = CGAL::to_double(squared_distance(p, p2));
weights[res[n]] -= weight_functor(p2, p, d2, r_max);
poisson_eliminate_impl::move_down(heap, heap_pos, heap_size, weights, heap_pos[res[n]]);
}
}
CGAL_assertion(heap_size == target_points);
if (progressive) {
target_points = target_points>>1;
r_max = r_max * std::pow(2.0, 1.0 / 3.0);
}
} while (progressive && target_points >= 3);
for (std::size_t i = 0; i < number_of_points; i++)
output++ = get(ipm, heap[i]);
}
} // namespace CGAL
#endif // CGAL_POISSON_ELIMINATE_H

View File

@ -23,6 +23,7 @@ endif()
create_single_source_cgal_program( "read_test.cpp" )
create_single_source_cgal_program( "test_read_write_point_set.cpp" )
create_single_source_cgal_program( "test_deprecated_io_point_set.cpp" )
create_single_source_cgal_program( "test_poisson_eliminate.cpp" )
create_single_source_cgal_program( "read_test_with_different_pmaps.cpp" )
create_single_source_cgal_program( "analysis_test.cpp" )
create_single_source_cgal_program( "remove_outliers_test.cpp" )

View File

@ -0,0 +1,88 @@
#include <CGAL/poisson_eliminate.h>
#include <vector>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/compute_average_spacing.h>
#include <CGAL/IO/read_points.h>
#include <CGAL/IO/write_points.h>
#include <CGAL/Real_timer.h>
template<class K>
void check(const std::vector<typename K::Point_3> &input, const std::vector<typename K::Point_3>& output, bool progressive) {
using FT = typename K::FT;
FT average_distance_input = CGAL::compute_average_spacing<CGAL::Parallel_if_available_tag>(input, 3);
std::cout << " " << average_distance_input << " average point distance" << std::endl;
if (output.size() < input.size()) {
FT average_distance_after = CGAL::compute_average_spacing<CGAL::Parallel_if_available_tag>(output, 3);
std::cout << " " << average_distance_after << " average point distance after elimination" << std::endl;
assert(average_distance_after > average_distance_input);
if (progressive) {
FT average_distance_progressive = CGAL::compute_average_spacing<CGAL::Parallel_if_available_tag>(CGAL::make_range(output.begin(), output.begin() + output.size() / 2), 3);
std::cout << " " << average_distance_progressive << " average point distance in the first half after progressive elimination" << std::endl;
assert(average_distance_progressive > average_distance_after);
}
}
}
void no_check(const std::vector<CGAL::Exact_predicates_exact_constructions_kernel::Point_3>&, const std::vector<CGAL::Exact_predicates_exact_constructions_kernel::Point_3>&, bool) {
}
template<class K, class Check_distance_functor>
void sampling(const std::string& filename, Check_distance_functor check_average_distance) {
using Point_3 = typename K::Point_3;
std::vector<Point_3> points;
if (!CGAL::IO::read_points(
filename,
std::back_inserter(points))) {
std::cerr << "Error: cannot read file!" << std::endl;
return;
}
std::cout << " " << points.size() << " input points" << std::endl;
std::size_t target_size = points.size() / 3;
bool progressive = true;
bool weight_limiting = false;
bool tiling = false;
std::vector<Point_3> output;
CGAL::Real_timer timer;
output.reserve(target_size);
timer.start();
CGAL::poisson_eliminate(points, target_size, std::back_inserter(output), CGAL::parameters::dimension(2).progressive(progressive).weight_limiting(weight_limiting).tiling(tiling));
timer.stop();
std::cout << timer.time() << std::endl;
std::cout << " " << output.size() << " points after elimination" << std::endl;
CGAL::IO::write_points("radar" + std::to_string(target_size) + "_2d.xyz", output, CGAL::parameters::stream_precision(17));
assert(output.size() == target_size);
check_average_distance(points, output, progressive);
std::cout << " done" << std::endl;
}
int main()
{
std::cout << "testing Simple_cartesian<double>" << std::endl;
sampling<CGAL::Simple_cartesian<double>>(CGAL::data_file_path("points_3/radar.xyz"), check< CGAL::Simple_cartesian<double>>);
std::cout << std::endl << "testing Exact_predicates_inexact_constructions_kernel" << std::endl;
sampling<CGAL::Exact_predicates_inexact_constructions_kernel>(CGAL::data_file_path("points_3/radar.xyz"), check< CGAL::Exact_predicates_inexact_constructions_kernel>);
std::cout << std::endl << "testing Exact_predicates_exact_constructions_kernel" << std::endl;
sampling<CGAL::Exact_predicates_exact_constructions_kernel>(CGAL::data_file_path("points_3/radar.xyz"), no_check);
return 0;
}

View File

@ -164,6 +164,11 @@ CGAL_add_named_parameter(region_primitive_map_t, region_primitive_map, region_pr
CGAL_add_named_parameter(postprocess_regions_t, postprocess_regions, postprocess_regions)
CGAL_add_named_parameter(sizing_function_t, sizing_function, sizing_function)
CGAL_add_named_parameter(bbox_scaling_t, bbox_scaling, bbox_scaling)
CGAL_add_named_parameter(weight_function_t, weight_function, weight_function)
CGAL_add_named_parameter(weight_limiting_t, weight_limiting, weight_limiting)
CGAL_add_named_parameter(progressive_t, progressive, progressive)
CGAL_add_named_parameter(tiling_t, tiling, tiling)
CGAL_add_named_parameter(dimension_t, dimension, dimension)
// List of named parameters that we use in the package 'Surface Mesh Simplification'
CGAL_add_named_parameter(get_cost_policy_t, get_cost_policy, get_cost)