diff --git a/Documentation/doc/biblio/cgal_manual.bib b/Documentation/doc/biblio/cgal_manual.bib index e95c3acf18a..8b45267c361 100644 --- a/Documentation/doc/biblio/cgal_manual.bib +++ b/Documentation/doc/biblio/cgal_manual.bib @@ -725,6 +725,17 @@ note = {\url{ttp://hal.inria.fr/inria-00090522}} ,pages = "325--338" } +@article{cgal:gcsa-nasr-13, + journal = {Computer Graphics Forum}, + title = {{Noise-Adaptive Shape Reconstruction from Raw Point Sets}}, + author = {Simon Giraudot and David Cohen-Steiner and Pierre Alliez }, + pages = {229-238}, + volume= {32}, + number= {5}, + year = {2013}, + DOI = {10.1111/cgf.12189}, +} + @manual{ cgal:g-gmpal-96 ,author = "T. Granlund" ,title = "{GNU MP}, The {GNU} Multiple Precision Arithmetic Library, diff --git a/Installation/changes.html b/Installation/changes.html index 136e66e9479..4741500449a 100644 --- a/Installation/changes.html +++ b/Installation/changes.html @@ -266,6 +266,10 @@ and src/ directories).
  • Function CGAL::remove_outliers() has an additional parameter based on a distance threshold to make it easier and more intuitive to use.
    • +
  • New functions for automatic scale estimations: either a global + scale or multiple local scales can be estimated for both 2D and 3D + point sets based on the assumption that they sample a curve in 2D + or a surface in 3D.
    • 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 5e55112ba46..6d08de81d6f 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 @@ -36,6 +36,10 @@ - `CGAL::read_ply_custom_points()` - `CGAL::read_xyz_points()` - `CGAL::compute_average_spacing()` +- `CGAL::estimate_global_k_neighbor_scale()` +- `CGAL::estimate_global_range_scale()` +- `CGAL::estimate_local_k_neighbor_scales()` +- `CGAL::estimate_local_range_scales()` - `CGAL::remove_outliers()` - `CGAL::grid_simplify_point_set()` - `CGAL::random_simplify_point_set()` 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 83e7a140978..9c09f824f3f 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 @@ -124,7 +124,7 @@ these attributes in user-defined containers. \cgalExample{Point_set_processing_3/read_ply_points_with_colors_example.cpp} -\section Point_set_processing_3Analysis Analysis +\section Point_set_processing_3Spacing Average Spacing Function `compute_average_spacing()` computes the average spacing of all input points to their `k` nearest neighbor points, @@ -147,6 +147,58 @@ found in other \cgal components: - `bounding_box()` - `Min_sphere_of_spheres_d` +\section Point_set_processing_3Scale Automatic Scale Estimation + +Point sets are often used to sample objects with a higher dimension, +typically a curve in 2D or a surface in 3D. In such cases, finding the +scale of the objet is crucial, that is to say finding the minimal +number of points (or the minimal local range) such that the subset of +points has the appearance of a curve in 2D or a surface in 3D +\cgalCite{cgal:gcsa-nasr-13}. + +\cgal provides 2 functions that automatically estimate the scale of a +2D point set sampling a curve or a 3D point set sampling a surface: + +- `estimate_global_k_neighbor_scale()` +- `estimate_global_range_scale()` + +Functions such as `jet_estimate_normals()` or +`remove_outliers()` require a K neighbor scale while others such as +`grid_simplify_point_set()` require a range +scale. `vcm_estimate_normals()` is an example of a function that +accepts both. + +In some specific cases, the scale of a point set might not be +homogeneous (for example if the point set contains variable +noise). \cgal also provides 2 functions that automatically estimate +the scales of a point set at a set of user-defined query points: + +- `estimate_local_k_neighbor_scales()` +- `estimate_local_range_scales()` + +The 4 functions presented here work both with 2D points and 3D points +and they shouldn't be used if the point sets do not sample a curve in +2D or a surface in 3D. + +\subsection Point_set_processing_3Example_scale_estimation_global Global Scale Example + +The following example reads a 3D point set in the `xyz` format and: + +- Estimates the K neighbor global scale +- Uses it to smooth the point set +- Estimates the range global scale +- Uses it to simplify the point set + +\cgalExample{Point_set_processing_3/scale_estimation_example.cpp} + +\subsection Point_set_processing_3Example_scale_estimation_local Local Scales Example + +This second example generates a 2D point set sampling a circle with +variable noise. It then estimates the scale at 3 different query +points in the domain. + +\cgalExample{Point_set_processing_3/scale_estimation_2d_example.cpp} + \section Point_set_processing_3OutlierRemoval Outlier Removal 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 ff119f543c5..25aeee34542 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 @@ -2,6 +2,8 @@ \example Point_set_processing_3/read_write_xyz_point_set_example.cpp \example Point_set_processing_3/read_ply_points_with_colors_example.cpp \example Point_set_processing_3/average_spacing_example.cpp +\example Point_set_processing_3/scale_estimation_example.cpp +\example Point_set_processing_3/scale_estimation_2d_example.cpp \example Point_set_processing_3/remove_outliers_example.cpp \example Point_set_processing_3/grid_simplification_example.cpp \example Point_set_processing_3/grid_simplify_indices.cpp 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 92dc32c4504..3954511d182 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 @@ -53,6 +53,8 @@ if ( CGAL_FOUND ) create_single_source_cgal_program( "read_write_xyz_point_set_example.cpp" ) create_single_source_cgal_program( "read_ply_points_with_colors_example.cpp" ) create_single_source_cgal_program( "remove_outliers_example.cpp" ) + create_single_source_cgal_program( "scale_estimation_example.cpp" ) + create_single_source_cgal_program( "scale_estimation_2d_example.cpp" ) create_single_source_cgal_program( "wlop_simplify_and_regularize_point_set_example.cpp" ) create_single_source_cgal_program( "edge_aware_upsample_point_set_example.cpp" ) create_single_source_cgal_program( "structuring_example.cpp" ) diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp new file mode 100644 index 00000000000..9f1599eb2a7 --- /dev/null +++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_2d_example.cpp @@ -0,0 +1,52 @@ +#include +#include +#include +#include + +#include +#include +#include + +// Types +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_2 Point_2; + +int main (int, char**) +{ + std::vector samples; + samples.reserve (100000); + + // Generate circle with gradually variable noise + // - noise-free for points with x close to (-1) + // - noisy for points with x close to (+1) + for (std::size_t i = 0; i < 100000; ++ i) + { + FT theta = CGAL::get_default_random().get_double(0, 2. * CGAL_PI); + FT noise = 0.5 * (std::cos(theta) + 1.) * CGAL::get_default_random().get_double(0., 0.2); + int mult = (CGAL::get_default_random().get_bool() ? 1 : -1); + samples.push_back (Point_2 (std::cos(theta) * (1. + mult * noise * noise), + std::sin(theta) * (1. + mult * noise * noise))); + } + + // Search for local scales on 3 different locations + std::vector queries; + queries.reserve (3); + queries.push_back (Point_2 (-1., 0.)); + queries.push_back (Point_2 (0., 1.)); + queries.push_back (Point_2 (1., 0.)); + + std::vector k_scales; + CGAL::estimate_local_k_neighbor_scales (samples.begin(), samples.end(), + queries.begin(), queries.end(), + std::back_inserter (k_scales)); + + // Display results + std::cerr << "K-Scales found:" << std::endl + << " - On noise-free region: " << k_scales[0] << std::endl // Should be small + << " - On moderately noisy region: " << k_scales[1] << std::endl // Should be higher + << " - On very noisy region: " << k_scales[2] << std::endl; // Should be even higher + + return EXIT_SUCCESS; +} + diff --git a/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp new file mode 100644 index 00000000000..bac9900fecf --- /dev/null +++ b/Point_set_processing_3/examples/Point_set_processing_3/scale_estimation_example.cpp @@ -0,0 +1,77 @@ +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +// Concurrency +#ifdef CGAL_LINKED_WITH_TBB +typedef CGAL::Parallel_tag Concurrency_tag; +#else +typedef CGAL::Sequential_tag Concurrency_tag; +#endif + +// Types +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_3 Point_3; + +int main (int argc, char** argv) +{ + + const char* fname = (argc>1)?argv[1]:"data/sphere_20k.xyz"; + + CGAL::Timer task_timer; + + std::vector points; + std::ifstream stream(fname); + + // read input + if (!(stream + && CGAL::read_xyz_points(stream, std::back_inserter(points)))) + { + std::cerr << "Error: can't read input file" << std::endl; + return EXIT_FAILURE; + } + + // estimate k scale + task_timer.start(); + std::size_t k_scale = CGAL::estimate_global_k_neighbor_scale (points.begin(), points.end()); + task_timer.stop(); + + // Example: use estimated k as scale for jet smoothing + CGAL::jet_smooth_point_set + (points.begin(), points.end(), + static_cast(k_scale)); + + // estimate range scale + task_timer.start(); + FT range_scale = CGAL::estimate_global_range_scale (points.begin(), points.end()); + task_timer.stop(); + + // Example: use estimated range for grid simplification + points.erase (CGAL::grid_simplify_point_set (points.begin(), points.end(), range_scale), + points.end()); + + + // print some informations on runtime + std::size_t memory = CGAL::Memory_sizer().virtual_size(); + double time = task_timer.time(); + + std::cout << "Scales computed in " << time << " second(s) using " + << (memory>>20) << " MiB of memory:" << std::endl; + std::cout << " * Global K scale: " << k_scale << std::endl; + std::cout << " * Global range scale: " << range_scale << std::endl; + + + return EXIT_SUCCESS; +} + diff --git a/Point_set_processing_3/include/CGAL/estimate_scale.h b/Point_set_processing_3/include/CGAL/estimate_scale.h new file mode 100644 index 00000000000..1ac7897e0f5 --- /dev/null +++ b/Point_set_processing_3/include/CGAL/estimate_scale.h @@ -0,0 +1,765 @@ +// Copyright (c) 2013 INRIA Sophia-Antipolis (France). +// Copyright (c) 2016 GeometryFactory Sarl (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// Author(s) : Simon Giraudot + +#ifndef CGAL_ESTIMATE_SCALE_H +#define CGAL_ESTIMATE_SCALE_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + + +namespace CGAL { + + +// ---------------------------------------------------------------------------- +// Private section +// ---------------------------------------------------------------------------- +/// \cond SKIP_IN_MANUAL +namespace internal { + +template +class Quick_multiscale_approximate_knn_distance +{ + +}; + + +template +class Quick_multiscale_approximate_knn_distance +{ + typedef typename Kernel::FT FT; + typedef Search_traits_3 Tree_traits; + typedef Orthogonal_k_neighbor_search Neighbor_search; + typedef typename Neighbor_search::Tree Tree; + typedef typename Neighbor_search::iterator Iterator; + + template + struct Pmap_unary_function : public std::unary_function + { + PointPMap point_pmap; + Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { } + const typename Kernel::Point_3& operator() (const ValueType& v) const { return get(point_pmap, v); } + }; + + std::size_t m_cluster_size; + std::vector m_trees; + std::vector m_weights; + std::vector m_precomputed_factor; + +public: + + template + Quick_multiscale_approximate_knn_distance (InputIterator first, + InputIterator beyond, + PointPMap point_pmap, + std::size_t cluster_size = 25) + : m_cluster_size (cluster_size) + { + typedef Pmap_unary_function::value_type, + PointPMap> Unary_f; + + m_trees.push_back (new Tree (boost::make_transform_iterator (first, Unary_f(point_pmap)), + boost::make_transform_iterator (beyond, Unary_f(point_pmap)))); + m_weights.push_back (1.); + std::size_t nb_pts = m_trees[0]->size(); + + std::size_t nb_trees = 0; + while (nb_pts > m_cluster_size) + { + nb_trees ++; + nb_pts /= m_cluster_size; + } + + m_trees.reserve (nb_trees); + m_weights.reserve (nb_trees); + + InputIterator first_unused = beyond; + + nb_pts = m_trees[0]->size(); + for (std::size_t i = 1; i < nb_trees; ++ i) + { + first_unused + = CGAL::hierarchy_simplify_point_set (first, first_unused, point_pmap, + static_cast(m_cluster_size), 1./3.); + + m_trees.push_back (new Tree(boost::make_transform_iterator (first, Unary_f(point_pmap)), + boost::make_transform_iterator (first_unused, Unary_f(point_pmap)))); + + m_weights.push_back (m_trees[0]->size() / (FT)(m_trees.back()->size())); + } + } + + ~Quick_multiscale_approximate_knn_distance() + { + for (std::size_t i = 0; i < m_trees.size(); ++ i) + delete m_trees[i]; + } + + template + std::size_t compute_k_scale (InputIterator query, PointPMap point_pmap) + { + std::size_t out; + FT dummy; + compute_scale (query, point_pmap, out, dummy); + return out; + } + + template + FT compute_range_scale (InputIterator query, PointPMap point_pmap) + { + std::size_t dummy; + FT out; + compute_scale (query, point_pmap, dummy, out); + return out; + } + + void precompute_factors () + { + FT nb = 0.; + for (std::size_t t = 0; t < m_trees.size(); ++ t) + { + std::size_t size = (t == (m_trees.size() - 1) + ? m_trees[t]->size() + : static_cast(m_weights[t+1] / m_weights[t])); + for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i) + { + nb += m_weights[t]; + if (nb < 6.) // do not consider values under 6 + continue; + m_precomputed_factor.push_back (0.91666666 * std::log (nb)); + } + } + } + + + template + void compute_scale (InputIterator query, PointPMap point_pmap, + std::size_t& k, FT& d) + { + if (m_precomputed_factor.empty()) + precompute_factors(); + + k = 0; + d = 0.; + + FT dist_min = (std::numeric_limits::max)(); + FT sum_sq_distances = 0.; + FT nb = 0.; + std::size_t index = 0; + for (std::size_t t = 0; t < m_trees.size(); ++ t) + { + Neighbor_search search (*(m_trees[t]), get(point_pmap, *query), + static_cast((t == (m_trees.size() - 1) + ? m_trees[t]->size() + : m_weights[t+1] / m_weights[t]))); + Iterator it = search.begin(); + + if (t != 0) // Skip first point except on first scale + ++ it; + + for (; it != search.end(); ++ it) + { + sum_sq_distances += m_weights[t] * it->second; + nb += m_weights[t]; + + if (nb < 6.) // do not consider values under 6 + continue; + + // sqrt(sum_sq_distances / nb) / nb^(5/12) + // Computed in log space with precomputed factor for time optimization + FT dist = 0.5 * std::log (sum_sq_distances) - m_precomputed_factor[index ++]; + + if (dist < dist_min) + { + dist_min = dist; + k = (std::size_t)nb; + d = it->second; + } + } + } + } + +}; + + +template +class Quick_multiscale_approximate_knn_distance +{ + typedef typename Kernel::FT FT; + typedef CGAL::Point_set_2 Point_set; + typedef typename Point_set::Vertex_handle Vertex_handle; + + template + struct Pmap_unary_function : public std::unary_function + { + PointPMap point_pmap; + Pmap_unary_function (PointPMap point_pmap) : point_pmap (point_pmap) { } + const typename Kernel::Point_2& operator() (const ValueType& v) const { return get(point_pmap, v); } + }; + + template + struct Pmap_to_3d + { + PointPMap point_pmap; + typedef typename Kernel::Point_3 value_type; + typedef const value_type& reference; + typedef typename boost::property_traits::key_type key_type; + typedef boost::lvalue_property_map_tag category; + Pmap_to_3d () { } + Pmap_to_3d (PointPMap point_pmap) + : point_pmap (point_pmap) { } + friend inline value_type get (const Pmap_to_3d& ppmap, key_type i) + { + typename Kernel::Point_2 p2 = get(ppmap.point_pmap, i); + return value_type (p2.x(), p2.y(), 0.); + } + + }; + + struct Sort_by_distance_to_point + { + const typename Kernel::Point_2& ref; + Sort_by_distance_to_point (const typename Kernel::Point_2& ref) : ref (ref) { } + bool operator() (const Vertex_handle& a, const Vertex_handle& b) + { + return (CGAL::squared_distance (a->point(), ref) + < CGAL::squared_distance (b->point(), ref)); + } + }; + + + std::size_t m_cluster_size; + std::vector m_point_sets; + std::vector m_weights; + std::vector m_precomputed_factor; + +public: + + template + Quick_multiscale_approximate_knn_distance (InputIterator first, + InputIterator beyond, + PointPMap point_pmap, + std::size_t cluster_size = 25) + : m_cluster_size (cluster_size) + { + typedef Pmap_unary_function::value_type, + PointPMap> Unary_f; + + m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)), + boost::make_transform_iterator (beyond, Unary_f(point_pmap)))); + m_weights.push_back (1.); + + std::size_t nb_pts = m_point_sets[0]->number_of_vertices(); + std::size_t nb_trees = 0; + while (nb_pts > m_cluster_size) + { + nb_trees ++; + nb_pts /= m_cluster_size; + } + + m_point_sets.reserve (nb_trees); + m_weights.reserve (nb_trees); + + InputIterator first_unused = beyond; + nb_pts = m_point_sets[0]->number_of_vertices(); + + for (std::size_t i = 1; i < nb_trees; ++ i) + { + first_unused + = CGAL::hierarchy_simplify_point_set (first, first_unused, Pmap_to_3d (point_pmap), + static_cast(m_cluster_size), 1./3.); + + m_point_sets.push_back (new Point_set (boost::make_transform_iterator (first, Unary_f(point_pmap)), + boost::make_transform_iterator (first_unused, Unary_f(point_pmap)))); + + m_weights.push_back (nb_pts / (FT)(m_point_sets.back()->number_of_vertices())); + } + + m_cluster_size = cluster_size; + } + + ~Quick_multiscale_approximate_knn_distance() + { + for (std::size_t i = 0; i < m_point_sets.size(); ++ i) + delete m_point_sets[i]; + } + + template + std::size_t compute_k_scale (InputIterator query, PointPMap point_pmap) + { + std::size_t out; + FT dummy; + compute_scale (query, point_pmap, out, dummy); + return out; + } + + template + FT compute_range_scale (InputIterator query, PointPMap point_pmap) + { + std::size_t dummy; + FT out; + compute_scale (query, point_pmap, dummy, out); + return out; + } + + void precompute_factors () + { + FT nb = 0.; + for (std::size_t t = 0; t < m_point_sets.size(); ++ t) + { + std::size_t size = (t == m_point_sets.size() - 1 + ? m_point_sets[t]->number_of_vertices() + : static_cast(m_weights[t+1] / m_weights[t])); + for (std::size_t i = (t == 0 ? 0 : 1); i < size; ++ i) + { + nb += m_weights[t]; + if (nb < 6.) // do not consider values under 6 + continue; + m_precomputed_factor.push_back (1.25 * std::log (nb)); + } + } + } + + template + void compute_scale (InputIterator query, PointPMap point_pmap, + std::size_t& k, FT& d) + { + if (m_precomputed_factor.empty()) + precompute_factors(); + + k = 0; + d = 0.; + + FT dist_min = (std::numeric_limits::max)(); + FT sum_sq_distances = 0.; + FT nb = 0.; + std::size_t index = 0; + + const typename Kernel::Point_2& pquery = get(point_pmap, *query); + for (std::size_t t = 0; t < m_point_sets.size(); ++ t) + { + std::size_t size = ((t == m_point_sets.size() - 1) + ? m_point_sets[t]->number_of_vertices() + : static_cast(m_weights[t+1] / m_weights[t])); + std::vector neighbors; + neighbors.reserve (size); + m_point_sets[t]->nearest_neighbors (pquery, size, std::back_inserter (neighbors)); + + std::sort (neighbors.begin(), neighbors.end(), + Sort_by_distance_to_point (pquery)); + for (std::size_t n = (t == 0 ? 0 : 1); n < neighbors.size(); ++ n) + { + FT sq_dist = CGAL::squared_distance (pquery, neighbors[n]->point()); + + sum_sq_distances += m_weights[t] * sq_dist; + nb += m_weights[t]; + + if (nb < 6.) // do not consider values under 6 + continue; + + // sqrt(sum_sq_distances / nb) / nb^(3/4) + // Computed in log space with precomputed factor for time optimization + FT dist = 0.5 * std::log (sum_sq_distances) - m_precomputed_factor[index ++]; + + if (dist < dist_min) + { + dist_min = dist; + k = (std::size_t)nb; + d = sq_dist; + } + } + } + } + +}; + +} /* namespace internal */ +/// \endcond + + + +// ---------------------------------------------------------------------------- +// Public section +// ---------------------------------------------------------------------------- + +/// \ingroup PkgPointSetProcessing + +/// Estimates the local scale in a K nearest neighbors sense on a set +/// of user-defined query points. The computed scales correspond to +/// the smallest scales such that the K subsets of points have the +/// appearance of a surface in 3D or the appearance of a curve in 2D +/// (see \ref Point_set_processing_3Scale). +/// +/// +/// @tparam SamplesInputIterator iterator over input sample points. +/// @tparam SamplesPointPMap is a model of `ReadablePropertyMap` with +/// value type `Point_3` or `Point_2`. It can +/// be omitted if the value type of `SamplesInputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam QueriesInputIterator iterator over points where scale +/// should be computed. +/// @tparam QueriesInputIterator is a model of `ReadablePropertyMap` +/// with value type `Point_3` or `Point_2`. It +/// can be omitted if the value type of `QueriesInputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam OutputIterator is used to store the computed scales. It accepts +/// values of type `std::size_t`. +/// @tparam Kernel Geometric traits class. It can be omitted and +/// deduced automatically from the value type of `SamplesPointPMap`. +/// +/// @note This function accepts both 2D and 3D points, but sample +/// points and query must have the same dimension. + +// This variant requires all parameters. +template +OutputIterator +estimate_local_k_neighbor_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + OutputIterator output, ///< output iterator to store the computed scales + const Kernel& /*kernel*/) ///< geometric traits. +{ + typedef typename boost::property_traits::value_type Point_d; + + // Build multi-scale KD-tree + internal::Quick_multiscale_approximate_knn_distance kdtree (first, beyond, samples_pmap); + + // Compute local scales everywhere + for (QueriesInputIterator it = first_query; it != beyond_query; ++ it) + *(output ++) = kdtree.compute_k_scale (it, queries_pmap); + + return output; +} + + +/// \ingroup PkgPointSetProcessing + +/// Estimates the global scale in a K nearest neighbors sense. The +/// computed scale corresponds to the smallest scale such that the K +/// subsets of points have the appearance of a surface in 3D or the +/// appearance of a curve in 2D (see \ref Point_set_processing_3Scale). +/// +/// +/// @tparam InputIterator iterator over input points. +/// @tparam PointPMap is a model of `ReadablePropertyMap` with +/// value type `Point_3` or `Point_2`. It can +/// be omitted if the value type of `InputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam Kernel Geometric traits class. It can be omitted and +/// deduced automatically from the value type of `PointPMap`. +/// +/// @note This function accepts both 2D and 3D points. +/// +/// @return The estimated scale in the K nearest neighbors sense. +// This variant requires all parameters. +template +std::size_t +estimate_global_k_neighbor_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond, ///< past-the-end iterator over the input points. + PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + const Kernel& kernel) ///< geometric traits. +{ + std::vector scales; + estimate_local_k_neighbor_scales (first, beyond, point_pmap, + first, beyond, point_pmap, + std::back_inserter (scales), + kernel); + std::sort (scales.begin(), scales.end()); + return scales[scales.size() / 2]; +} + + +/// \ingroup PkgPointSetProcessing + +/// Estimates the local scale in a range sense on a set of +/// user-defined query points. The computed scales correspond to the +/// smallest scales such that the subsets of points included in the +/// sphere range have the appearance of a surface in 3D or the +/// appearance of a curve in 2D (see \ref Point_set_processing_3Scale). +/// +/// +/// @tparam SamplesInputIterator iterator over input sample points. +/// @tparam SamplesPointPMap is a model of `ReadablePropertyMap` with +/// value type `Point_3` or `Point_2`. It can +/// be omitted if the value type of `SamplesInputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam QueriesInputIterator iterator over points where scale +/// should be computed. +/// @tparam QueriesInputIterator is a model of `ReadablePropertyMap` +/// with value type `Point_3` or `Point_2`. It +/// can be omitted if the value type of `QueriesInputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam OutputIterator is used to store the computed scales. It accepts +/// values of type `Kernel::FT`. +/// @tparam Kernel Geometric traits class. It can be omitted and +/// deduced automatically from the value type of `SamplesPointPMap`. +/// +/// @note This function accepts both 2D and 3D points, but sample +/// points and query must have the same dimension. + +// This variant requires all parameters. +template +OutputIterator +estimate_local_range_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + OutputIterator output, ///< output iterator to store the computed scales + const Kernel& /*kernel*/) ///< geometric traits. +{ + typedef typename boost::property_traits::value_type Point_d; + + // Build multi-scale KD-tree + internal::Quick_multiscale_approximate_knn_distance kdtree (first, beyond, samples_pmap); + + // Compute local scales everywhere + for (QueriesInputIterator it = first_query; it != beyond_query; ++ it) + *(output ++) = kdtree.compute_range_scale (it, queries_pmap); + + return output; +} + + +/// \ingroup PkgPointSetProcessing + +/// Estimates the global scale in a range sense. The computed scale +/// corresponds to the smallest scale such that the subsets of points +/// inside the sphere range have the appearance of a surface in 3D or +/// the appearance of a curve in 2D (see \ref Point_set_processing_3Scale). +/// +/// +/// @tparam InputIterator iterator over input points. +/// @tparam PointPMap is a model of `ReadablePropertyMap` with +/// value type `Point_3` or `Point_2`. It can +/// be omitted if the value type of `InputIterator` is +/// convertible to `Point_3` or to `Point_2`. +/// @tparam Kernel Geometric traits class. It can be omitted and +/// deduced automatically from the value type of `PointPMap`. +/// +/// @note This function accepts both 2D and 3D points. +/// +/// @return The estimated scale in the range sense. +// This variant requires all parameters. +template +typename Kernel::FT +estimate_global_range_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond, ///< past-the-end iterator over the input points. + PointPMap point_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_3 + const Kernel& kernel) ///< geometric traits. +{ + std::vector scales; + estimate_local_range_scales (first, beyond, point_pmap, + first, beyond, point_pmap, + std::back_inserter (scales), + kernel); + std::sort (scales.begin(), scales.end()); + return std::sqrt (scales[scales.size() / 2]); +} + + +// ---------------------------------------------------------------------------- +// Useful overloads +// ---------------------------------------------------------------------------- +/// \cond SKIP_IN_MANUAL + +template +OutputIterator +estimate_local_k_neighbor_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + OutputIterator output) ///< output iterator to store the computed scales +{ + typedef typename boost::property_traits::value_type Point; + typedef typename Kernel_traits::Kernel Kernel; + return estimate_local_k_neighbor_scales (first, beyond, samples_pmap, first_query, beyond_query, + queries_pmap, output, Kernel()); +} + +template +OutputIterator +estimate_local_k_neighbor_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + OutputIterator output) ///< output iterator to store the computed scales +{ + return estimate_local_k_neighbor_scales + (first, beyond, + make_identity_property_map (typename std::iterator_traits::value_type()), + first_query, beyond_query, + make_identity_property_map (typename std::iterator_traits::value_type()), + output); + } + + +template +std::size_t +estimate_global_k_neighbor_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond, ///< past-the-end iterator over the input points. + PointPMap point_pmap) ///< property map: value_type of InputIterator -> Point_3 or Point_2 +{ + typedef typename boost::property_traits::value_type Point; + typedef typename Kernel_traits::Kernel Kernel; + return estimate_global_k_neighbor_scale (first, beyond, point_pmap, Kernel()); +} + +template +std::size_t +estimate_global_k_neighbor_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond) ///< past-the-end iterator over the input points. +{ + return estimate_global_k_neighbor_scale + (first, beyond, make_identity_property_map (typename std::iterator_traits::value_type())); +} + + +template +OutputIterator +estimate_local_range_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + SamplesPointPMap samples_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + QueriesPointPMap queries_pmap, ///< property map: value_type of InputIterator -> Point_3 or Point_2 + OutputIterator output) ///< output iterator to store the computed scales +{ + typedef typename boost::property_traits::value_type Point; + typedef typename Kernel_traits::Kernel Kernel; + return estimate_local_range_scales(first, beyond, samples_pmap, first_query, beyond_query, + queries_pmap, output, Kernel()); +} + + +template +OutputIterator +estimate_local_range_scales( + SamplesInputIterator first, ///< iterator over the first input sample. + SamplesInputIterator beyond, ///< past-the-end iterator over the input samples. + QueriesInputIterator first_query, ///< iterator over the first point where scale must be estimated + QueriesInputIterator beyond_query, ///< past-the-end iterator over the points where scale must be estimated + OutputIterator output) ///< output iterator to store the computed scales +{ + return estimate_local_range_scales + (first, beyond, + make_identity_property_map (typename std::iterator_traits::value_type()), + first_query, beyond_query, + make_identity_property_map (typename std::iterator_traits::value_type()), + output); +} + + + +template +typename Kernel_traits::value_type>::Kernel::FT +estimate_global_range_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond, ///< past-the-end iterator over the input points. + PointPMap point_pmap) ///< property map: value_type of InputIterator -> Point_3 or Point_3 +{ + typedef typename boost::property_traits::value_type Point; + typedef typename Kernel_traits::Kernel Kernel; + return estimate_global_range_scale (first, beyond, point_pmap, Kernel()); +} + + + +template +typename Kernel_traits::value_type>::Kernel::FT +estimate_global_range_scale( + InputIterator first, ///< iterator over the first input point. + InputIterator beyond) ///< past-the-end iterator over the input points. +{ + return estimate_global_range_scale + (first, beyond, make_identity_property_map (typename std::iterator_traits::value_type())); + +} +/// \endcond + +} //namespace CGAL + +#endif // CGAL_ESTIMATE_SCALE_3_H