diff --git a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt index 18916e87b3f..e2c64e40efb 100644 --- a/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt +++ b/Point_set_processing_3/doc/Point_set_processing_3/NamedParameters.txt @@ -93,7 +93,7 @@ along sharp features. The range of possible values is `[0, 1]`. \n \cgalNPBegin{neighbor_radius} \anchor PSP_neighbor_radius is the spherical neighborhood radius.\n \b Type: floating scalar value\n -Default value: `-1`, in which case the functions that use it will compute an adapted value automatically. +Default value: `-1`, in which case it is not used. \cgalNPEnd \cgalNPBegin{number_of_output_points} \anchor PSP_number_of_output_points 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 426329c0a29..23cbd2c4c57 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 @@ -281,11 +281,10 @@ points has the appearance of a curve in 2D or a surface in 3D - `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. +Functions such as `grid_simplify_point_set()` require a range scale +while `jet_estimate_normals()`, `remove_outliers()` or +`vcm_estimate_normals()` are examples of functions that accepts both a +K neighbor scale or a range scale. In some specific cases, the scale of a point set might not be homogeneous (for example if the point set contains variable @@ -325,7 +324,9 @@ points in the domain. Function `remove_outliers()` deletes a user-specified fraction of outliers from an input point set. More specifically, it sorts the input points in increasing order of average squared distances to their -`k` nearest neighbors and deletes the points with largest value. +nearest neighbors and deletes the points with largest value. The user +can either specify a fixed number of nearest neighbors or a fixed +spherical neighborhood radius. \subsection Point_set_processing_3Example_outlier_removal Example @@ -446,13 +447,16 @@ Two smoothing functions are devised to smooth an input point set. Function `jet_smooth_point_set()` smooths the input point set by projecting each point onto a smooth parametric surface patch -(so-called jet surface) fitted over its `k` nearest neighbors. +(so-called jet surface) fitted over its nearest neighbors. Function `bilateral_smooth_point_set()` smooths the input point set by -iteratively projecting each point onto the implicit surface patch fitted over its `k` nearest neighbors. +iteratively projecting each point onto the implicit surface patch fitted over its nearest neighbors. Bilateral projection preserves sharp features according to the normal (gradient) information. Normals are thus required as input. For more details, see section 4 of \cgalCite{ear-2013}. +For both functions, the user can either specify a fixed number of +nearest neighbors or a fixed spherical neighborhood radius. + \subsection Point_set_processing_3Example_jet_smoothing Jet Smoothing Example The following example generates a set of 9 points close to the `xy` @@ -487,13 +491,13 @@ functions provide an estimate of the normal to \b S at each point. The result is an unoriented normal vector for each input point. Function `jet_estimate_normals()` estimates the normal direction -at each point from the input set by fitting a jet surface over its `k` +at each point from the input set by fitting a jet surface over its nearest neighbors. The default jet is a quadric surface. This algorithm is well suited to point sets scattered over curved surfaces. Function `pca_estimate_normals()` estimates the normal direction at each point from the set by linear least squares fitting of a plane -over its `k` nearest neighbors. This algorithm is simpler and +over its nearest neighbors. This algorithm is simpler and faster than `jet_estimate_normals()`. Function `vcm_estimate_normals()` estimates the normal direction @@ -502,13 +506,16 @@ of the point set. This algorithm is more complex and slower than the previous algorithms. It is based on the article \cgalCite{cgal:mog-vbcfe-11}. +For these three functions, the user can either specify a fixed number +of nearest neighbors or a fixed spherical neighborhood radius. + \section Point_set_processing_3NormalOrientation Normal Orientation Function `mst_orient_normals()` orients the normals of a set of points with unoriented normals using the method described by Hoppe et al. in Surface reconstruction from unorganized points \cgalCite{cgal:hddms-srup-92}. More specifically, this method constructs a -Riemannian graph over the input points (the graph of the `k` +Riemannian graph over the input points (the graph of the nearest neighbor points) and propagates a seed normal orientation within a minimum spanning tree computed over this graph. The result is an oriented normal vector for each input unoriented normal, except for @@ -521,8 +528,9 @@ Normal orientation of a sampled cube surface. Left: unoriented normals. Right: o \subsection Point_set_processing_3Example_normals Example The following example reads a point set from a file, estimates the -normals through PCA over the 6 nearest neighbors and orients the -normals: +normals through PCA (either over the 18 nearest neighbors or using a +spherical neighborhood radius of twice the average spacing) and +orients the normals: \cgalExample{Point_set_processing_3/normals_example.cpp} diff --git a/Point_set_processing_3/examples/Point_set_processing_3/normals_example.cpp b/Point_set_processing_3/examples/Point_set_processing_3/normals_example.cpp index 85df999904d..25a28429661 100644 --- a/Point_set_processing_3/examples/Point_set_processing_3/normals_example.cpp +++ b/Point_set_processing_3/examples/Point_set_processing_3/normals_example.cpp @@ -1,4 +1,5 @@ #include +#include #include #include #include @@ -42,11 +43,31 @@ int main(int argc, char*argv[]) // Note: pca_estimate_normals() requiresa range of points // as well as property maps to access each point's position and normal. const int nb_neighbors = 18; // K-nearest neighbors = 3 rings - CGAL::pca_estimate_normals - (points, nb_neighbors, - CGAL::parameters::point_map(CGAL::First_of_pair_property_map()). - normal_map(CGAL::Second_of_pair_property_map())); + if (argc > 2 && std::strcmp(argv[2], "-r") == 0) // Use a fixed neighborhood radius + { + // First compute a spacing using the K parameter + double spacing + = CGAL::compute_average_spacing + (points, nb_neighbors, + CGAL::parameters::point_map(CGAL::First_of_pair_property_map())); + + // Then, estimate normals with a fixed radius + CGAL::pca_estimate_normals + (points, + 0, // when using a neighborhood radius, K=0 means no limit on the number of neighbors returns + CGAL::parameters::point_map(CGAL::First_of_pair_property_map()). + normal_map(CGAL::Second_of_pair_property_map()). + neighbor_radius(2. * spacing)); // use 2*spacing as neighborhood radius + } + else // Use a fixed number of neighbors + { + CGAL::pca_estimate_normals + (points, nb_neighbors, + CGAL::parameters::point_map(CGAL::First_of_pair_property_map()). + normal_map(CGAL::Second_of_pair_property_map())); + } + // Orients normals. // Note: mst_orient_normals() requires a range of points // as well as property maps to access each point's position and normal. diff --git a/Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/neighbor_query.h b/Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/neighbor_query.h new file mode 100644 index 00000000000..85921d2f2eb --- /dev/null +++ b/Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/neighbor_query.h @@ -0,0 +1,113 @@ +// Copyright (c) 2019 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$ +// SPDX-License-Identifier: GPL-3.0+ +// +// Author(s) : Simon Giraudot + +#ifndef CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H +#define CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H + +#include + +#include +#include +#include +#include + +#include + +namespace CGAL { +namespace Point_set_processing_3 { +namespace internal { + +struct Maximum_points_reached_exception : public std::exception { }; + +template +void neighbor_query (const Point& query, + const CGAL::Kd_tree& tree, + unsigned int k, + FT neighbor_radius, + PointContainer& points) +{ + typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; + typedef typename Neighbor_search::iterator Search_iterator; + typedef CGAL::Fuzzy_sphere Sphere; + + if (neighbor_radius != FT(0)) + { + Sphere fs (query, neighbor_radius, 0, tree.traits()); + + // if k=0, no limit on the number of neighbors returned + if (k == 0) + k = (std::numeric_limits::max)(); + + try + { + std::function back_insert_with_limit + = [&](const Point& point) -> void + { + points.push_back (point); + if (points.size() == k) + throw Maximum_points_reached_exception(); + }; + + auto function_output_iterator + = boost::make_function_output_iterator (back_insert_with_limit); + + tree.search (function_output_iterator, fs); + } + catch (const Maximum_points_reached_exception&) + { } + + // Fallback, if less than 3 points are return, search for the 3 + // first points + if (points.size() < 3) + k = 3; + // Else, no need to search for K nearest neighbors + else + k = 0; + } + + if (k != 0) + { + // Gather set of (k+1) neighboring points. + // Perform k+1 queries (as in point set, the query point is + // output first). Search may be aborted if k is greater + // than number of input points. + points.reserve(k+1); + Neighbor_search search(tree,query,k+1); + Search_iterator search_iterator = search.begin(); + unsigned int i; + for(i=0;i<(k+1);i++) + { + if(search_iterator == search.end()) + break; // premature ending + points.push_back(search_iterator->first); + search_iterator++; + } + CGAL_point_set_processing_precondition(points.size() >= 1); + } +} + +} } } // namespace CGAL::Point_set_processing_3::internal + +#endif // CGAL_PSP_INTERNAL_NEIGHBOR_QUERY_H diff --git a/Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h b/Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h index 5b12c20d770..93836493036 100644 --- a/Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h +++ b/Point_set_processing_3/include/CGAL/bilateral_smooth_point_set.h @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -185,58 +186,6 @@ compute_denoise_projection( return Pwn(update_point, update_normal); } -/// Computes neighbors from kdtree. -/// -/// \pre `k >= 2`. -/// -/// @tparam Kernel Geometric traits class. -/// @tparam Tree KD-tree. -/// -/// @return neighbors pwn of query point. -template < typename Kernel, - typename Tree> -std::vector, - CGAL_PSP3_DEFAULT_ALLOCATOR > > -compute_kdtree_neighbors( - const CGAL::Point_with_normal_3& query, ///< 3D point - const Tree& tree, ///< KD-tree - unsigned int k ///< number of neighbors -) -{ - // basic geometric types - typedef CGAL::Point_with_normal_3 Pwn; - - // types for K nearest neighbors search - typedef bilateral_smooth_point_set_internal::Kd_tree_traits Tree_traits; - typedef CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Search_iterator; - - // performs k + 1 queries (if unique the query point is - // output first). search may be aborted when k is greater - // than number of input points - Neighbor_search search(tree, query, k+1); - Search_iterator search_iterator = search.begin(); - ++search_iterator; - unsigned int i; - std::vector - , CGAL_PSP3_DEFAULT_ALLOCATOR > - > neighbor_pwns; - - for(i = 0; i < (k+1); ++i) - { - if(search_iterator == search.end()) - break; // premature ending - - Pwn pwn = search_iterator->first; - neighbor_pwns.push_back(pwn); - ++search_iterator; - } - - // output - return neighbor_pwns; -} - - /// Computes max-spacing of one query point from K nearest neighbors. /// /// \pre `k >= 2`. @@ -302,6 +251,7 @@ class Compute_pwns_neighbors typedef typename Kernel::FT FT; unsigned int m_k; + FT m_neighbor_radius; const Tree & m_tree; const Pwns & m_pwns; Pwns_neighbors & m_pwns_neighbors; @@ -309,11 +259,12 @@ class Compute_pwns_neighbors cpp11::atomic& interrupted; public: - Compute_pwns_neighbors(unsigned int k, const Tree &tree, + Compute_pwns_neighbors(unsigned int k, FT neighbor_radius, const Tree &tree, const Pwns &pwns, Pwns_neighbors &neighbors, cpp11::atomic& advancement, cpp11::atomic& interrupted) - : m_k(k), m_tree(tree), m_pwns(pwns), m_pwns_neighbors(neighbors) + : m_k(k), m_neighbor_radius (neighbor_radius), m_tree(tree) + , m_pwns(pwns), m_pwns_neighbors(neighbors) , advancement (advancement), interrupted (interrupted) {} void operator() ( const tbb::blocked_range& r ) const @@ -322,9 +273,10 @@ public: { if (interrupted) break; + + CGAL::Point_set_processing_3::internal::neighbor_query + (m_pwns[i], m_tree, m_k, m_neighbor_radius, m_pwns_neighbors[i]); - m_pwns_neighbors[i] = bilateral_smooth_point_set_internal:: - compute_kdtree_neighbors(m_pwns[i], m_tree, m_k); ++ advancement; } } @@ -392,7 +344,7 @@ public: \ingroup PkgPointSetProcessing3Algorithms This function smooths an input point set by iteratively projecting each - point onto the implicit surface patch fitted over its k nearest neighbors. + point onto the implicit surface patch fitted over its nearest neighbors. Bilateral projection preserves sharp features according to the normal (gradient) information. Both point positions and normals will be modified. For more details, please see section 4 in \cgalCite{ear-2013}. @@ -422,6 +374,12 @@ public: If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd \cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type `geom_traits::Vector_3`.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood radius. If + provided, the neighborhood of a query point is computed with a fixed spherical + radius instead of a fixed number of neighbors. In that case, the parameter + `k` is used as a limit on the number of points returned by each spherical + query (to avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{sharpness_angle} controls the sharpness of the result.\cgalParamEnd \cgalParamBegin{callback} an instance of `std::function`. It is called regularly when the @@ -477,7 +435,8 @@ bilateral_smooth_point_set( PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap()); - + FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), FT(0)); + // copy points and normals Pwns pwns; for(typename PointRange::iterator it = points.begin(); it != points.end(); ++it) @@ -546,7 +505,7 @@ bilateral_smooth_point_set( Point_set_processing_3::internal::Parallel_callback parallel_callback (callback, 2 * nb_points); - Compute_pwns_neighbors f(k, tree, pwns, pwns_neighbors, + Compute_pwns_neighbors f(k, neighbor_radius, tree, pwns, pwns_neighbors, parallel_callback.advancement(), parallel_callback.interrupted()); tbb::parallel_for(tbb::blocked_range(0, nb_points), f); @@ -570,8 +529,9 @@ bilateral_smooth_point_set( std::size_t nb = 0; for(pwn_iter = pwns.begin(); pwn_iter != pwns.end(); ++pwn_iter, ++pwns_iter, ++ nb) { - *pwns_iter = bilateral_smooth_point_set_internal:: - compute_kdtree_neighbors(*pwn_iter, tree, k); + CGAL::Point_set_processing_3::internal::neighbor_query + (*pwn_iter, tree, k, neighbor_radius, *pwns_iter); + if (callback && !callback ((nb+1) / double(2. * nb_points))) return std::numeric_limits::quiet_NaN(); } diff --git a/Point_set_processing_3/include/CGAL/jet_estimate_normals.h b/Point_set_processing_3/include/CGAL/jet_estimate_normals.h index 9d627e9c8c1..dfcb1190306 100644 --- a/Point_set_processing_3/include/CGAL/jet_estimate_normals.h +++ b/Point_set_processing_3/include/CGAL/jet_estimate_normals.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -74,39 +75,22 @@ typename Kernel::Vector_3 jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at Tree& tree, ///< KD-tree unsigned int k, ///< number of neighbors + typename Kernel::FT neighbor_radius, unsigned int degree_fitting) { // basic geometric types typedef typename Kernel::Point_3 Point; - // types for K nearest neighbors search - typedef typename CGAL::Search_traits_3 Tree_traits; - typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Search_iterator; - // types for jet fitting typedef Monge_via_jet_fitting< Kernel, Simple_cartesian, SvdTraits> Monge_jet_fitting; typedef typename Monge_jet_fitting::Monge_form Monge_form; - // Gather set of (k+1) neighboring points. - // Perform k+1 queries (as in point set, the query point is - // output first). Search may be aborted if k is greater - // than number of input points. - std::vector points; points.reserve(k+1); - Neighbor_search search(tree,query,k+1); - Search_iterator search_iterator = search.begin(); - unsigned int i; - for(i=0;i<(k+1);i++) - { - if(search_iterator == search.end()) - break; // premature ending - points.push_back(search_iterator->first); - search_iterator++; - } - CGAL_point_set_processing_precondition(points.size() >= 1); - + std::vector points; + CGAL::Point_set_processing_3::internal::neighbor_query + (query, tree, k, neighbor_radius, points); + // performs jet fitting Monge_jet_fitting monge_fit; const unsigned int degree_monge = 1; // we seek for normal and not more. @@ -120,10 +104,12 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute #ifdef CGAL_LINKED_WITH_TBB template class Jet_estimate_normals { + typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; const Tree& tree; const unsigned int k; + const FT neighbor_radius; const unsigned int degree_fitting; const std::vector& input; std::vector& output; @@ -131,11 +117,13 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute cpp11::atomic& interrupted; public: - Jet_estimate_normals(Tree& tree, unsigned int k, std::vector& points, + Jet_estimate_normals(Tree& tree, unsigned int k, FT neighbor_radius, + std::vector& points, unsigned int degree_fitting, std::vector& output, cpp11::atomic& advancement, cpp11::atomic& interrupted) - : tree(tree), k (k), degree_fitting (degree_fitting), input (points), output (output) + : tree(tree), k (k), neighbor_radius (neighbor_radius) + , degree_fitting (degree_fitting), input (points), output (output) , advancement (advancement) , interrupted (interrupted) { } @@ -146,7 +134,7 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute { if (interrupted) break; - output[i] = CGAL::internal::jet_estimate_normal(input[i], tree, k, degree_fitting); + output[i] = CGAL::internal::jet_estimate_normal(input[i], tree, k, neighbor_radius, degree_fitting); ++ advancement; } } @@ -168,7 +156,7 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute /** \ingroup PkgPointSetProcessing3Algorithms Estimates normal directions of the range of `points` - using jet fitting on the k nearest neighbors. + using jet fitting on the nearest neighbors. The output normals are randomly oriented. \pre `k >= 2` @@ -188,6 +176,12 @@ jet_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd \cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type `geom_traits::Vector_3`.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood radius. If + provided, the neighborhood of a query point is computed with a fixed spherical + radius instead of a fixed number of neighbors. In that case, the parameter + `k` is used as a limit on the number of points returned by each spherical + query (to avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{degree_fitting} degree of jet fitting.\cgalParamEnd \cgalParamBegin{svd_traits} template parameter for the class `Monge_via_jet_fitting`. If \ref thirdpartyEigen "Eigen" 3.2 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined, @@ -220,6 +214,7 @@ jet_estimate_normals( typedef typename Point_set_processing_3::GetPointMap::type PointMap; typedef typename Point_set_processing_3::GetNormalMap::type NormalMap; typedef typename Point_set_processing_3::GetK::Kernel Kernel; + typedef typename Kernel::FT FT; typedef typename GetSvdTraits::type SvdTraits; CGAL_static_assertion_msg(!(boost::is_same& callback = choose_param(get_param(np, internal_np::callback), std::function()); @@ -251,7 +248,7 @@ jet_estimate_normals( CGAL_point_set_processing_precondition(points.begin() != points.end()); // precondition: at least 2 nearest neighbors - CGAL_point_set_processing_precondition(k >= 2); + CGAL_point_set_processing_precondition(k >= 2 || neighbor_radius > FT(0)); std::size_t memory = CGAL::Memory_sizer().virtual_size(); CGAL_TRACE(" %ld Mb allocated\n", memory>>20); CGAL_TRACE(" Creates KD-tree\n"); @@ -282,7 +279,8 @@ jet_estimate_normals( std::vector normals (kd_tree_points.size (), CGAL::NULL_VECTOR); CGAL::internal::Jet_estimate_normals - f (tree, k, kd_tree_points, degree_fitting, normals, + f (tree, k, neighbor_radius, + kd_tree_points, degree_fitting, normals, parallel_callback.advancement(), parallel_callback.interrupted()); tbb::parallel_for(tbb::blocked_range(0, kd_tree_points.size ()), f); @@ -301,7 +299,7 @@ jet_estimate_normals( { Vector normal = internal::jet_estimate_normal( get(point_map,*it), - tree, k, degree_fitting); + tree, k, neighbor_radius, degree_fitting); put(normal_map, *it, normal); // normal_map[it] = normal if (callback && !callback ((nb+1) / double(kd_tree_points.size()))) diff --git a/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h b/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h index 498097cdb3b..a2069571ffe 100644 --- a/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h +++ b/Point_set_processing_3/include/CGAL/jet_smooth_point_set.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -74,39 +75,22 @@ jet_smooth_point( const typename Kernel::Point_3& query, ///< 3D point to project Tree& tree, ///< KD-tree const unsigned int k, ///< number of neighbors. + typename Kernel::FT neighbor_radius, const unsigned int degree_fitting, const unsigned int degree_monge) { // basic geometric types typedef typename Kernel::Point_3 Point; - // types for K nearest neighbors search - typedef typename CGAL::Search_traits_3 Tree_traits; - typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Search_iterator; - // types for jet fitting typedef Monge_via_jet_fitting< Kernel, Simple_cartesian, SvdTraits> Monge_jet_fitting; typedef typename Monge_jet_fitting::Monge_form Monge_form; - - // Gather set of (k+1) neighboring points. - // Performs k + 1 queries (if unique the query point is - // output first). Search may be aborted if k is greater - // than number of input points. - std::vector points; points.reserve(k+1); - Neighbor_search search(tree,query,k+1); - Search_iterator search_iterator = search.begin(); - unsigned int i; - for(i=0;i<(k+1);i++) - { - if(search_iterator == search.end()) - break; // premature ending - points.push_back(search_iterator->first); - search_iterator++; - } - CGAL_point_set_processing_precondition(points.size() >= 1); + + std::vector points; + CGAL::Point_set_processing_3::internal::neighbor_query + (query, tree, k, neighbor_radius, points); // performs jet fitting Monge_jet_fitting monge_fit; @@ -123,6 +107,7 @@ jet_smooth_point( typedef typename Kernel::Point_3 Point; const Tree& tree; const unsigned int k; + const typename Kernel::FT neighbor_radius; unsigned int degree_fitting; unsigned int degree_monge; const std::vector& input; @@ -131,12 +116,14 @@ jet_smooth_point( cpp11::atomic& interrupted; public: - Jet_smooth_pwns (Tree& tree, unsigned int k, std::vector& points, - unsigned int degree_fitting, unsigned int degree_monge, std::vector& output, + Jet_smooth_pwns (Tree& tree, unsigned int k, typename Kernel::FT neighbor_radius, + std::vector& points, + unsigned int degree_fitting, unsigned int degree_monge, std::vector& output, cpp11::atomic& advancement, cpp11::atomic& interrupted) - : tree(tree), k (k), degree_fitting (degree_fitting), - degree_monge (degree_monge), input (points), output (output) + : tree(tree), k (k), neighbor_radius(neighbor_radius) + , degree_fitting (degree_fitting) + , degree_monge (degree_monge), input (points), output (output) , advancement (advancement) , interrupted (interrupted) { } @@ -148,6 +135,7 @@ jet_smooth_point( if (interrupted) break; output[i] = CGAL::internal::jet_smooth_point(input[i], tree, k, + neighbor_radius, degree_fitting, degree_monge); ++ advancement; @@ -170,7 +158,7 @@ jet_smooth_point( /** \ingroup PkgPointSetProcessing3Algorithms - Smoothes the range of `points` using jet fitting on the k + Smoothes the range of `points` using jet fitting on the nearest neighbors and reprojection onto the jet. As this method relocates the points, it should not be called on containers sorted w.r.t. point locations. @@ -190,6 +178,12 @@ jet_smooth_point( \cgalNamedParamsBegin \cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`. If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood radius. If + provided, the neighborhood of a query point is computed with a fixed spherical + radius instead of a fixed number of neighbors. In that case, the parameter + `k` is used as a limit on the number of points returned by each spherical + query (to avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{degree_fitting} degree of jet fitting.\cgalParamEnd \cgalParamBegin{degree_monge} Monge degree.\cgalParamEnd \cgalParamBegin{svd_traits} template parameter for the class `Monge_via_jet_fitting`. If @@ -228,6 +222,8 @@ jet_smooth_point_set( "Error: no SVD traits"); PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); + typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), + typename Kernel::FT(0)); unsigned int degree_fitting = choose_param(get_param(np, internal_np::degree_fitting), 2); unsigned int degree_monge = choose_param(get_param(np, internal_np::degree_monge), 2); const std::function& callback = choose_param(get_param(np, internal_np::callback), @@ -271,8 +267,8 @@ jet_smooth_point_set( std::vector mutated_points (kd_tree_points.size (), CGAL::ORIGIN); CGAL::internal::Jet_smooth_pwns - f (tree, k, kd_tree_points, degree_fitting, degree_monge, - mutated_points, + f (tree, k, neighbor_radius, kd_tree_points, degree_fitting, degree_monge, + mutated_points, parallel_callback.advancement(), parallel_callback.interrupted()); tbb::parallel_for(tbb::blocked_range(0, kd_tree_points.size ()), f); @@ -293,7 +289,7 @@ jet_smooth_point_set( const typename boost::property_traits::reference p = get(point_map, *it); put(point_map, *it , internal::jet_smooth_point( - p,tree,k,degree_fitting,degree_monge) ); + p,tree,k,neighbor_radius,degree_fitting,degree_monge) ); if (callback && !callback ((nb+1) / double(kd_tree_points.size()))) break; } diff --git a/Point_set_processing_3/include/CGAL/mst_orient_normals.h b/Point_set_processing_3/include/CGAL/mst_orient_normals.h index 8cc551a31bc..4ed8e68d7e7 100644 --- a/Point_set_processing_3/include/CGAL/mst_orient_normals.h +++ b/Point_set_processing_3/include/CGAL/mst_orient_normals.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -319,6 +320,7 @@ create_riemannian_graph( IndexMap index_map, ///< property map ForwardIterator -> index ConstrainedMap constrained_map, ///< property map ForwardIterator -> bool unsigned int k, ///< number of neighbors + typename Kernel::FT neighbor_radius, const Kernel& /*kernel*/) ///< geometric traits. { // Input points types @@ -331,7 +333,6 @@ create_riemannian_graph( typedef Euclidean_distance_vertex_handle_3 KDistance; typedef Orthogonal_k_neighbor_search Neighbor_search; typedef typename Neighbor_search::Tree Tree; - typedef typename Neighbor_search::iterator Search_iterator; // Riemannian_graph types typedef internal::Riemannian_graph Riemannian_graph; @@ -395,21 +396,15 @@ create_riemannian_graph( std::size_t it_index = get(index_map,it); Vector_ref it_normal_vector = get(normal_map,*it); - // Gather set of (k+1) neighboring points. - // Perform k+1 queries (as in point set, the query point is - // output first). Search may be aborted if k is greater - // than number of input points. - Point_ref point = get(point_map, *it); Point_vertex_handle_3 point_wrapper(point.x(), point.y(), point.z(), it); - Neighbor_search search(*tree, point_wrapper, k+1); - Search_iterator search_iterator = search.begin(); - for(std::size_t i=0;i<(k+1);i++) - { - if(search_iterator == search.end()) - break; // premature ending + std::vector neighbor_points; + CGAL::Point_set_processing_3::internal::neighbor_query + (point_wrapper, *tree, k, neighbor_radius, neighbor_points); - ForwardIterator neighbor = search_iterator->first; + for (std::size_t i = 0; i < neighbor_points.size(); ++ i) + { + ForwardIterator neighbor = neighbor_points[i]; std::size_t neighbor_index = get(index_map,neighbor); if (neighbor_index > it_index) // undirected graph { @@ -431,8 +426,6 @@ create_riemannian_graph( weight = 0; // safety check riemannian_graph_weight_map[e] = (float)weight; } - - search_iterator++; } // Check if point is source @@ -595,6 +588,12 @@ create_mst_graph( If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd \cgalParamBegin{normal_map} a model of `ReadWritePropertyMap` with value type `geom_traits::Vector_3`.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood radius. If + provided, the neighborhood of a query point is computed with a fixed spherical + radius instead of a fixed number of neighbors. In that case, the parameter + `k` is used as a limit on the number of points returned by each spherical + query (to avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{point_is_constrained_map} a model of `ReadablePropertyMap` with value type `bool`. Points with a `true` value will be used as seed points: their normal will be considered as already oriented, it won't be altered and it will be propagated to its neighbors. If this parameter is omitted, @@ -628,6 +627,8 @@ mst_orient_normals( PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap()); + typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), + typename Kernel::FT(0)); ConstrainedMap constrained_map = choose_param(get_param(np, internal_np::point_is_constrained), ConstrainedMap()); Kernel kernel; @@ -677,12 +678,14 @@ mst_orient_normals( point_map, normal_map, kernel)), k, + neighbor_radius, kernel); else riemannian_graph = create_riemannian_graph(points.begin(), points.end(), point_map, normal_map, index_map, constrained_map, k, + neighbor_radius, kernel); // Creates a Minimum Spanning Tree starting at source_point diff --git a/Point_set_processing_3/include/CGAL/pca_estimate_normals.h b/Point_set_processing_3/include/CGAL/pca_estimate_normals.h index 2cdabd3e57d..92d93f86b33 100644 --- a/Point_set_processing_3/include/CGAL/pca_estimate_normals.h +++ b/Point_set_processing_3/include/CGAL/pca_estimate_normals.h @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -57,7 +58,6 @@ namespace CGAL { /// \cond SKIP_IN_MANUAL namespace internal { - /// Estimates normal direction using linear least /// squares fitting of a plane on the K nearest neighbors. /// @@ -73,34 +73,17 @@ template < typename Kernel, typename Kernel::Vector_3 pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute the normal at const Tree& tree, ///< KD-tree - unsigned int k) ///< number of neighbors + unsigned int k, ///< number of neighbors + typename Kernel::FT neighbor_radius) { // basic geometric types typedef typename Kernel::Point_3 Point; typedef typename Kernel::Plane_3 Plane; - // types for K nearest neighbors search - typedef typename CGAL::Search_traits_3 Tree_traits; - typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Search_iterator; - - // Gather set of (k+1) neighboring points. - // Perform k+1 queries (as in point set, the query point is - // output first). Search may be aborted if k is greater - // than number of input points. - std::vector points; points.reserve(k+1); - Neighbor_search search(tree,query,k+1); - Search_iterator search_iterator = search.begin(); - unsigned int i; - for(i=0;i<(k+1);i++) - { - if(search_iterator == search.end()) - break; // premature ending - points.push_back(search_iterator->first); - search_iterator++; - } - CGAL_point_set_processing_precondition(points.size() >= 1); - + std::vector points; + CGAL::Point_set_processing_3::internal::neighbor_query + (query, tree, k, neighbor_radius, points); + // performs plane fitting by point-based PCA Plane plane; linear_least_squares_fitting_3(points.begin(),points.end(),plane,Dimension_tag<0>()); @@ -113,21 +96,25 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute #ifdef CGAL_LINKED_WITH_TBB template class PCA_estimate_normals { + typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; typedef typename Kernel::Vector_3 Vector; const Tree& tree; const unsigned int k; + const FT neighbor_radius; const std::vector& input; std::vector& output; cpp11::atomic& advancement; cpp11::atomic& interrupted; public: - PCA_estimate_normals(Tree& tree, unsigned int k, std::vector& points, - std::vector& output, + PCA_estimate_normals(Tree& tree, unsigned int k, FT neighbor_radius, + std::vector& points, + std::vector& output, cpp11::atomic& advancement, cpp11::atomic& interrupted) - : tree(tree), k (k), input (points), output (output) + : tree(tree), k (k), neighbor_radius (neighbor_radius) + , input (points), output (output) , advancement (advancement) , interrupted (interrupted) { } @@ -138,7 +125,7 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute { if (interrupted) break; - output[i] = CGAL::internal::pca_estimate_normal(input[i], tree, k); + output[i] = CGAL::internal::pca_estimate_normal(input[i], tree, k, neighbor_radius); ++ advancement; } } @@ -159,7 +146,7 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute /** \ingroup PkgPointSetProcessing3Algorithms Estimates normal directions of the range of `points` - by linear least squares fitting of a plane over the k nearest neighbors. + by linear least squares fitting of a plane over the nearest neighbors. The output normals are randomly oriented. \pre `k >= 2` @@ -174,20 +161,29 @@ pca_estimate_normal(const typename Kernel::Point_3& query, ///< point to compute \param k number of neighbors \param np optional sequence of \ref psp_namedparameters "Named Parameters" among the ones listed below. - \cgalNamedParamsBegin - \cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`. - If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd - \cgalParamBegin{normal_map} a model of `WritablePropertyMap` with value type + \cgalNamedParamsBegin \cgalParamBegin{point_map} a model of + `ReadablePropertyMap` with value type `geom_traits::Point_3`. If + this parameter is omitted, + `CGAL::Identity_property_map` is + used.\cgalParamEnd \cgalParamBegin{normal_map} a model of + `WritablePropertyMap` with value type `geom_traits::Vector_3`.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood + radius. If provided, the neighborhood of a query point is + computed with a fixed spherical radius instead of a fixed number + of neighbors. In that case, the parameter `k` is used as a limit + on the number of points returned by each spherical query (to + avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{callback} an instance of - `std::function`. It is called regularly when the - algorithm is running: the current advancement (between 0. and - 1.) is passed as parameter. If it returns `true`, then the - algorithm continues its execution normally; if it returns - `false`, the algorithm is stopped and the remaining normals are - left unchanged.\cgalParamEnd - \cgalParamBegin{geom_traits} an instance of a geometric traits class, model of `Kernel`\cgalParamEnd - \cgalNamedParamsEnd + `std::function`. It is called regularly when the + algorithm is running: the current advancement (between 0. and 1.) + is passed as parameter. If it returns `true`, then the algorithm + continues its execution normally; if it returns `false`, the + algorithm is stopped and the remaining normals are left + unchanged.\cgalParamEnd \cgalParamBegin{geom_traits} an instance + of a geometric traits class, model of `Kernel`\cgalParamEnd + \cgalNamedParamsEnd */ template ::type PointMap; typedef typename Point_set_processing_3::GetNormalMap::type NormalMap; typedef typename Point_set_processing_3::GetK::Kernel Kernel; + typedef typename Kernel::FT FT; CGAL_static_assertion_msg(!(boost::is_same::NoMap>::value), @@ -213,6 +210,7 @@ pca_estimate_normals( PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); NormalMap normal_map = choose_param(get_param(np, internal_np::normal_map), NormalMap()); + FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), FT(0)); const std::function& callback = choose_param(get_param(np, internal_np::callback), std::function()); @@ -263,7 +261,7 @@ pca_estimate_normals( std::vector normals (kd_tree_points.size (), CGAL::NULL_VECTOR); CGAL::internal::PCA_estimate_normals - f (tree, k, kd_tree_points, normals, + f (tree, k, neighbor_radius, kd_tree_points, normals, parallel_callback.advancement(), parallel_callback.interrupted()); tbb::parallel_for(tbb::blocked_range(0, kd_tree_points.size ()), f); @@ -283,7 +281,7 @@ pca_estimate_normals( Vector normal = internal::pca_estimate_normal( get(point_map,*it), tree, - k); + k, neighbor_radius); put(normal_map, *it, normal); // normal_map[it] = normal if (callback && !callback ((nb+1) / double(kd_tree_points.size()))) diff --git a/Point_set_processing_3/include/CGAL/remove_outliers.h b/Point_set_processing_3/include/CGAL/remove_outliers.h index 64118d54e76..5e781975b4c 100644 --- a/Point_set_processing_3/include/CGAL/remove_outliers.h +++ b/Point_set_processing_3/include/CGAL/remove_outliers.h @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -63,33 +64,16 @@ typename Kernel::FT compute_avg_knn_sq_distance_3( const typename Kernel::Point_3& query, ///< 3D point to project Tree& tree, ///< KD-tree - unsigned int k) ///< number of neighbors + unsigned int k, ///< number of neighbors + typename Kernel::FT neighbor_radius) { // geometric types typedef typename Kernel::FT FT; typedef typename Kernel::Point_3 Point; - // types for K nearest neighbors search - typedef typename CGAL::Search_traits_3 Tree_traits; - typedef typename CGAL::Orthogonal_k_neighbor_search Neighbor_search; - typedef typename Neighbor_search::iterator Search_iterator; - - // Gather set of (k+1) neighboring points. - // Perform k+1 queries (if in point set, the query point is - // output first). Search may be aborted if k is greater - // than number of input points. - std::vector points; points.reserve(k+1); - Neighbor_search search(tree,query,k+1); - Search_iterator search_iterator = search.begin(); - unsigned int i; - for(i=0;i<(k+1);i++) - { - if(search_iterator == search.end()) - break; // premature ending - points.push_back(search_iterator->first); - search_iterator++; - } - CGAL_point_set_processing_precondition(points.size() >= 1); + std::vector points; + CGAL::Point_set_processing_3::internal::neighbor_query + (query, tree, k, neighbor_radius, points); // compute average squared distance typename Kernel::Compute_squared_distance_3 sqd; @@ -111,7 +95,7 @@ compute_avg_knn_sq_distance_3( /** \ingroup PkgPointSetProcessing3Algorithms Removes outliers: - - computes average squared distance to the K nearest neighbors, + - computes average squared distance to the nearest neighbors, - and sorts the points in increasing order of average distance. This method modifies the order of input points so as to pack all remaining points first, @@ -130,6 +114,13 @@ compute_avg_knn_sq_distance_3( \cgalNamedParamsBegin \cgalParamBegin{point_map} a model of `ReadablePropertyMap` with value type `geom_traits::Point_3`. If this parameter is omitted, `CGAL::Identity_property_map` is used.\cgalParamEnd + \cgalParamBegin{neighbor_radius} spherical neighborhood + radius. If provided, the neighborhood of a query point is + computed with a fixed spherical radius instead of a fixed number + of neighbors. In that case, the parameter `k` is used as a limit + on the number of points returned by each spherical query (to + avoid overly large number of points in high density areas). If no + limit is wanted, use `k=0`.\cgalParamEnd \cgalParamBegin{threshold_percent} maximum percentage of points to remove.\cgalParamEnd \cgalParamBegin{threshold_distance} minimum distance for a point to be considered as outlier (distance here is the square root of the average squared distance to K nearest neighbors).\cgalParamEnd @@ -169,6 +160,8 @@ remove_outliers( typedef typename Point_set_processing_3::GetK::Kernel Kernel; PointMap point_map = choose_param(get_param(np, internal_np::point_map), PointMap()); + typename Kernel::FT neighbor_radius = choose_param(get_param(np, internal_np::neighbor_radius), + typename Kernel::FT(0)); double threshold_percent = choose_param(get_param(np, internal_np::threshold_percent), 10.); double threshold_distance = choose_param(get_param(np, internal_np::threshold_distance), 0.); const std::function& callback = choose_param(get_param(np, internal_np::callback), @@ -213,7 +206,7 @@ remove_outliers( { FT sq_distance = internal::compute_avg_knn_sq_distance_3( get(point_map,*it), - tree, k); + tree, k, neighbor_radius); sorted_points.insert( std::make_pair(sq_distance, *it) ); if (callback && !callback ((nb+1) / double(kd_tree_points.size()))) return points.end(); diff --git a/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h b/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h index b3b5789bab7..00dcfda9160 100644 --- a/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h +++ b/Point_set_processing_3/include/CGAL/wlop_simplify_and_regularize_point_set.h @@ -442,8 +442,8 @@ public: \cgalParamBegin{neighbor_radius} spherical neighborhood radius. This is a key parameter that needs to be finely tuned. The result will be irregular if too small, but a larger value will impact the runtime. In practice, choosing a radius such that the neighborhood of each sample point includes at least two rings - of neighboring sample points gives satisfactory result. The default value is set to 8 times the average - spacing of the point set.\cgalParamEnd + of neighboring sample points gives satisfactory result. If this parameter is not provided, it is + automatically set to 8 times the average spacing of the point set.\cgalParamEnd \cgalParamBegin{number_of_iterations} number of iterations to solve the optimsation problem. The default value is 35. More iterations give a more regular result but increase the runtime.\cgalParamEnd \cgalParamBegin{require_uniform_sampling} an optional preprocessing, which will give better result if the