diff --git a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt index 4a6c6ec2996..a6d0001f17b 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/PackageDescription.txt @@ -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)} diff --git a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt index 7d62380935d..e17cf4bc5af 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/Point_set_processing_3.txt @@ -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 +Sample Elimination for Generating Poisson Disk Sample Sets \cgalCite{cgal:y-sefpdss}. + \subsection Point_set_processing_3Example_grid_simplification Grid Simplification Example @@ -755,6 +759,37 @@ for more details. We provide below a speed-up chart generated using the parallel 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: +- *dimensions*: The dimensions parameter specifies 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 *weight_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. r_max has to be provided if a custom *weight_function* is used. Only points within a distance of \f$r_{max}\f$ is 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 *dimensions* parameter and the number of input points \f$N\f$: + +\f$dimensions = 2:\qquad\f$ \f$r_{max} = \sqrt{\frac{V}{2\sqrt{3}N}}\f$ + +\f$dimensions = 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 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} +
+ +
+\cgalFigureCaptionBegin{Point_set_processing_3figPoisson_elimination} +Poisson elimination on point cloud with 21k points and results with \f$\frac{1}{2}\f$, \f$\frac{1}{4}\f$ and \f$\frac{1}{8}\f$ of the input size. +\cgalFigureCaptionEnd + +\subsection 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} \section Point_set_processing_3Smoothing Smoothing diff --git a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt index ea670af1b62..2d2c4413b30 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/examples.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/examples.txt @@ -23,4 +23,5 @@ \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 */ diff --git a/Point_set_processing_3/doc/Point_set_processing_3/fig/poisson_elimination.png b/Point_set_processing_3/doc/Point_set_processing_3/fig/poisson_elimination.png new file mode 100644 index 00000000000..9d023bd3850 Binary files /dev/null and b/Point_set_processing_3/doc/Point_set_processing_3/fig/poisson_elimination.png differ diff --git a/Point_set_processing_3/doc/Point_set_processing_3/fig/simplification_comparison.jpg b/Point_set_processing_3/doc/Point_set_processing_3/fig/simplification_comparison.jpg index feb1b40643f..caebff22e58 100644 Binary files a/Point_set_processing_3/doc/Point_set_processing_3/fig/simplification_comparison.jpg and b/Point_set_processing_3/doc/Point_set_processing_3/fig/simplification_comparison.jpg differ diff --git a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt index fc7a6b237f5..04d2e33493c 100644 --- a/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt +++ b/Point_set_processing_3/examples/Point_set_processing_3/CMakeLists.txt @@ -141,13 +141,3 @@ if(TARGET CGAL::Eigen3_support) else() message(STATUS "NOTICE: Some of the executables in this directory require Eigen 3.1 (or greater), and will not be compiled.") endif() - -find_path(CYCODEBASE_INCLUDE "cySampleElim.h" DOC "Path to the cyCodeBase directory (https://github.com/cemyuksel/cyCodeBase)") - -if(NOT CYCODEBASE_INCLUDE) - message( STATUS "poisson_eliminate_example will be compiled without the cyCodeBase library (https://github.com/cemyuksel/cyCodeBase)") -else() - message(CYCODEBASE_INCLUDE="${CYCODEBASE_INCLUDE}") - target_compile_definitions( poisson_eliminate_example PRIVATE CGAL_USE_CY ) - target_include_directories( poisson_eliminate_example AFTER PRIVATE ${CYCODEBASE_INCLUDE} ) -endif() diff --git a/Point_set_processing_3/examples/Point_set_processing_3/poisson_eliminate_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/poisson_eliminate_example.cpp index 1dca7de62e7..9dfc5ee7adf 100644 --- a/Point_set_processing_3/examples/Point_set_processing_3/poisson_eliminate_example.cpp +++ b/Point_set_processing_3/examples/Point_set_processing_3/poisson_eliminate_example.cpp @@ -8,12 +8,15 @@ #include #include #include -#include typedef CGAL::Simple_cartesian K; typedef K::Point_3 Point_3; -void sampling(const std::string& filename) { +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 points; if (!CGAL::IO::read_points( @@ -24,54 +27,13 @@ void sampling(const std::string& filename) { return; } - std::size_t target_size = points.size() / 5; - - bool progressive = true; - bool weight_limiting = false; - bool tiling = false; + std::size_t target_size = std::size_t(points.size() * size_factor); std::vector output; - -#ifdef CGAL_USE_CY - std::vector inputPoints, outputPoints; - std::vector out; - - for (int i = 0; i < points.size(); ++i) - inputPoints.push_back(cy::Vec3d(CGAL::to_double(points[i].x()), CGAL::to_double(points[i].y()), CGAL::to_double(points[i].z()))); - - outputPoints.resize(target_size); - - CGAL::Bbox_3 bb = CGAL::bbox_3(points.begin(), points.end()); - cy::Vec3d bl(bb.xmin(), bb.ymin(), bb.zmin()); - cy::Vec3d tr(bb.xmax(), bb.ymax(), bb.zmax()); - - cy::WeightedSampleElimination< cy::Vec3d, double, 3, int > wse; - wse.SetBoundsMin(bl); - wse.SetBoundsMax(tr); - wse.SetTiling(tiling); -#endif - output.reserve(target_size); - CGAL::Real_timer timer; - timer.start(); - CGAL::poisson_eliminate(points, target_size, std::back_inserter(output), CGAL::parameters::progressive(progressive).weight_limiting(weight_limiting).tiling(tiling)); - timer.stop(); - std::cout << std::endl << timer.time() << std::endl; + + CGAL::poisson_eliminate(points, target_size, std::back_inserter(output)); + CGAL::IO::write_points("out.xyz", output, CGAL::parameters::stream_precision(17)); - -#ifdef CGAL_USE_CY - timer.reset(); - timer.start(); - wse.SetWeightLimiting(weight_limiting); - wse.Eliminate(inputPoints.data(), int(inputPoints.size()), - outputPoints.data(), int(outputPoints.size()), progressive); - timer.stop(); - std::cout << timer.time() << std::endl; - - for (const cy::Vec3d& p : outputPoints) { - out.push_back(Point_3(p.x, p.y, p.z)); - } - CGAL::IO::write_points("out_ref.xyz", out, CGAL::parameters::stream_precision(17)); -#endif } @@ -79,8 +41,10 @@ int main(int argc, char* argv[]) { if (argc < 2) sampling(CGAL::data_file_path("points_3/radar.xyz")); - else + else if (argc < 3) sampling(argv[1]); + else if (argc < 4) + sampling(argv[1], std::atof(argv[2])); return 0; } diff --git a/Point_set_processing_3/include/CGAL/poisson_eliminate.h b/Point_set_processing_3/include/CGAL/poisson_eliminate.h index d5fffa9f7bc..1cfbd780745 100644 --- a/Point_set_processing_3/include/CGAL/poisson_eliminate.h +++ b/Point_set_processing_3/include/CGAL/poisson_eliminate.h @@ -118,11 +118,12 @@ public: } }; +template 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()(double d2, double r_max) { + 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); @@ -197,6 +198,12 @@ void pop_heap(std::vector& heap, std::vector& heap_pos \cgalParamDefault{`CGAL::Identity_property_map`} \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} @@ -206,11 +213,13 @@ void pop_heap(std::vector& heap, std::vector& heap_pos \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 \ref Point_set_processing_3PoissonElimination.} \cgalParamNEnd - \cgalParamNBegin{weight_functor} - \cgalParamDescription{a weight functor that calculates the weight of a neighbor point based on its squared distance and the `maximum_radius`.} - \cgalParamType{an instance of `std::function`.} + \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`.} + \cgalParamDefault{See \ref Point_set_processing_3PoissonElimination.} \cgalParamNEnd \cgalParamNBegin{geom_traits} @@ -227,8 +236,8 @@ void poisson_eliminate(PointRange points, std::size_t number_of_points, OutputIt using NP_helper = Point_set_processing_3_np_helper; using PointMap = typename NP_helper::Point_map; - using Point = typename boost::property_traits::value_type; - using GeomTraits = typename Kernel_traits::Kernel; + using Point = typename NP_helper::Point; + using GeomTraits = typename NP_helper::Geom_traits; using FT = typename GeomTraits::FT; using IPM = internal::Indexed_extended_point_map; PointMap point_map = NP_helper::get_point_map(points, np); @@ -255,13 +264,13 @@ void poisson_eliminate(PointRange points, std::size_t number_of_points, OutputIt // named parameter for tiling const bool tiling = parameters::choose_parameter(parameters::get_parameter(np, internal_np::tiling), false); const unsigned int ambient_dimension = CGAL::Ambient_dimension::value; - const unsigned int dimension = parameters::choose_parameter(parameters::get_parameter(np, internal_np::dimension), ambient_dimension); + const unsigned int dimension = parameters::choose_parameter(parameters::get_parameter(np, internal_np::dimension), 2); // named parameter for r_max double r_max = CGAL::to_double(parameters::choose_parameter(parameters::get_parameter(np, internal_np::maximum_radius), 2 * internal::get_maximum_radius(dimension, number_of_points, domain_size))); double r_min = CGAL::to_double(weight_limiting ? internal::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_functor), internal::Weight_functor(r_min, alpha)); + auto weight_functor = parameters::choose_parameter(parameters::get_parameter(np, internal_np::weight_function), internal::Weight_functor(r_min, alpha)); std::size_t heap_size = points.size(); std::vector tiling_points; @@ -325,7 +334,7 @@ void poisson_eliminate(PointRange points, std::size_t number_of_points, OutputIt const Point p2 = get(ipm, res[n]); double d2 = CGAL::to_double((p - p2).squared_length()); - weights[i] += weight_functor(d2, r_max); + weights[i] += weight_functor(p, p2, d2, r_max); } } @@ -356,7 +365,7 @@ void poisson_eliminate(PointRange points, std::size_t number_of_points, OutputIt const Point p2 = get(point_map, points[res[n]]); double d2 = CGAL::to_double((p - p2).squared_length()); - weights[res[n]] -= weight_functor(d2, r_max); + weights[res[n]] -= weight_functor(p2, p, d2, r_max); internal::move_down(heap, heap_pos, heap_size, weights, heap_pos[res[n]]); } diff --git a/Point_set_processing_3/test/Point_set_processing_3/test_poisson_eliminate.cpp b/Point_set_processing_3/test/Point_set_processing_3/test_poisson_eliminate.cpp index d2951b02003..caf56855830 100644 --- a/Point_set_processing_3/test/Point_set_processing_3/test_poisson_eliminate.cpp +++ b/Point_set_processing_3/test/Point_set_processing_3/test_poisson_eliminate.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include template @@ -49,17 +50,24 @@ void sampling(const std::string& filename, Check_distance_functor check_average_ std::cout << " " << points.size() << " input points" << std::endl; - std::size_t target_size = points.size() / 5; + std::size_t target_size = points.size() / 3; bool progressive = true; bool weight_limiting = false; bool tiling = false; std::vector output; + CGAL::Real_timer timer; + output.reserve(target_size); - CGAL::poisson_eliminate(points, target_size, std::back_inserter(output), CGAL::parameters::progressive(progressive).weight_limiting(weight_limiting).tiling(tiling)); + 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); @@ -73,11 +81,11 @@ int main(int argc, char* argv[]) std::cout << "testing Simple_cartesian" << std::endl; sampling>(CGAL::data_file_path("points_3/radar.xyz"), check< CGAL::Simple_cartesian>); - std::cout << std::endl << "testing Exact_predicates_inexact_constructions_kernel" << std::endl; - sampling(CGAL::data_file_path("points_3/radar.xyz"), check< CGAL::Exact_predicates_inexact_constructions_kernel>); + std::cout << std::endl << "testing Exact_predicates_inexact_constructions_kernel" << std::endl; + sampling(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::data_file_path("points_3/radar.xyz"), no_check); + std::cout << std::endl << "testing Exact_predicates_exact_constructions_kernel" << std::endl; + sampling(CGAL::data_file_path("points_3/radar.xyz"), no_check); return 0; } diff --git a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h index 4c3eeddf956..02dde581778 100644 --- a/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h +++ b/STL_Extension/include/CGAL/STL_Extension/internal/parameters_interface.h @@ -164,7 +164,7 @@ 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_functor_t, weight_functor, weight_functor) +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)