1 // Copyright (c) 2014 INRIA Sophia-Antipolis (France).
2 // All rights reserved.
3 //
4 // This file is part of CGAL (www.cgal.org).
5 //
6 // $URL: https://github.com/CGAL/cgal/blob/v5.3/Point_set_processing_3/include/CGAL/vcm_estimate_normals.h $
7 // $Id: vcm_estimate_normals.h 9939011 2020-06-21T15:54:21+02:00 Mael Rouxel-Labbé
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 // Author(s) : Jocelyn Meyron and Quentin Mérigot
11 //
12
13 #ifndef CGAL_VCM_ESTIMATE_NORMALS_H
14 #define CGAL_VCM_ESTIMATE_NORMALS_H
15
16 #include <CGAL/license/Point_set_processing_3.h>
17
18 #include <CGAL/disable_warnings.h>
19
20 #include <CGAL/Point_set_processing_3/internal/Voronoi_covariance_3/voronoi_covariance_3.h>
21
22 #include <CGAL/property_map.h>
23 #include <CGAL/point_set_processing_assertions.h>
24 #include <CGAL/Delaunay_triangulation_3.h>
25 #include <CGAL/Kd_tree.h>
26 #include <CGAL/Search_traits_3.h>
27 #include <CGAL/Search_traits_adapter.h>
28 #include <CGAL/property_map.h>
29 #include <CGAL/Orthogonal_k_neighbor_search.h>
30 #include <CGAL/Fuzzy_sphere.h>
31
32 #include <CGAL/boost/graph/Named_function_parameters.h>
33 #include <CGAL/boost/graph/named_params_helper.h>
34
35 #include <CGAL/Default_diagonalize_traits.h>
36
37 #include <iterator>
38 #include <vector>
39
40 namespace CGAL {
41
42
43 // ----------------------------------------------------------------------------
44 // Private section
45 // ----------------------------------------------------------------------------
46 namespace internal {
47
48 /// @cond SKIP_IN_MANUAL
49 /// Computes the VCM for each point in the property map.
50 /// The matrix is computed by intersecting the Voronoi cell
51 /// of a point and a sphere whose radius is `offset_radius` and discretized
52 /// by `N` planes.
53 ///
54 /// @tparam ForwardIterator iterator over input points.
55 /// @tparam PointMap is a model of `ReadablePropertyMap` with a value_type = `Kernel::Point_3`.
56 /// @tparam K Geometric traits class.
57 /// @tparam Covariance Covariance matrix type. It is similar to an array with a length of 6.
58 template < typename ForwardIterator,
59 typename PointMap,
60 class K,
61 class Covariance
62 >
63 void
vcm_offset(ForwardIterator first,ForwardIterator beyond,PointMap point_map,std::vector<Covariance> & cov,double offset_radius,std::size_t N,const K &)64 vcm_offset (ForwardIterator first, ///< iterator over the first input point.
65 ForwardIterator beyond, ///< past-the-end iterator over the input points.
66 PointMap point_map, ///< property map: value_type of ForwardIterator -> Point_3.
67 std::vector<Covariance> &cov, ///< vector of covariance matrices.
68 double offset_radius, ///< radius of the sphere.
69 std::size_t N, ///< number of planes used to discretize the sphere.
70 const K & /*kernel*/) ///< geometric traits.
71 {
72 // Sphere discretization
73 typename CGAL::Voronoi_covariance_3::Sphere_discretization<K> sphere(offset_radius, N);
74
75 // Compute the Delaunay Triangulation
76 std::vector<typename K::Point_3> points;
77 points.reserve(std::distance(first, beyond));
78 for (ForwardIterator it = first; it != beyond; ++it)
79 points.push_back(get(point_map, *it));
80
81 typedef Delaunay_triangulation_3<K> DT;
82 DT dt(points.begin(), points.end());
83
84 cov.clear();
85 cov.reserve(points.size());
86 // Compute the VCM
87 for (typename std::vector<typename K::Point_3>::iterator
88 it = points.begin(); it != points.end(); ++it)
89 {
90 typename DT::Vertex_handle vh = dt.nearest_vertex(*it);
91 cov.push_back(
92 Voronoi_covariance_3::voronoi_covariance_3(dt, vh, sphere)
93 );
94 }
95 }
96 /// @endcond
97
98 /// @cond SKIP_IN_MANUAL
99 // Convolve using a radius.
100 template < class ForwardIterator,
101 class PointMap,
102 class K,
103 class Covariance
104 >
105 void
vcm_convolve(ForwardIterator first,ForwardIterator beyond,PointMap point_map,const std::vector<Covariance> & cov,std::vector<Covariance> & ncov,double convolution_radius,const K &)106 vcm_convolve (ForwardIterator first,
107 ForwardIterator beyond,
108 PointMap point_map,
109 const std::vector<Covariance> &cov,
110 std::vector<Covariance> &ncov,
111 double convolution_radius,
112 const K &)
113 {
114 typedef std::pair<typename K::Point_3, std::size_t> Tree_point;
115 typedef First_of_pair_property_map< Tree_point > Tree_map;
116 typedef Search_traits_3<K> Traits_base;
117 typedef Search_traits_adapter<Tree_point, Tree_map, Traits_base> Traits;
118 typedef Kd_tree<Traits> Tree;
119 typedef Fuzzy_sphere<Traits> Fuzzy_sphere;
120
121 // Kd tree
122 Tree tree;
123 tree.reserve(cov.size());
124 std::size_t i=0;
125 for (ForwardIterator it = first; it != beyond; ++it, ++i)
126 tree.insert( Tree_point(get(point_map, *it), i) );
127
128 // Convolving
129 ncov.clear();
130 ncov.reserve(cov.size());
131 for (ForwardIterator it = first; it != beyond; ++it) {
132 std::vector<Tree_point> nn;
133 tree.search(std::back_inserter(nn),
134 Fuzzy_sphere (get(point_map, *it), convolution_radius));
135
136 Covariance m;
137 std::fill(m.begin(), m.end(), typename K::FT(0));
138 for (std::size_t k = 0; k < nn.size(); ++k)
139 {
140 std::size_t index = nn[k].second;
141 for (int i=0; i<6; ++i)
142 m[i] += cov[index][i];
143 }
144 ncov.push_back(m);
145 }
146 }
147 /// @endcond
148
149 /// @cond SKIP_IN_MANUAL
150 // Convolve using neighbors.
151 template < class ForwardIterator,
152 class PointMap,
153 class K,
154 class Covariance
155 >
156 void
vcm_convolve(ForwardIterator first,ForwardIterator beyond,PointMap point_map,const std::vector<Covariance> & cov,std::vector<Covariance> & ncov,unsigned int nb_neighbors_convolve,const K &)157 vcm_convolve (ForwardIterator first,
158 ForwardIterator beyond,
159 PointMap point_map,
160 const std::vector<Covariance> &cov,
161 std::vector<Covariance> &ncov,
162 unsigned int nb_neighbors_convolve,
163 const K &)
164 {
165 typedef std::pair<typename K::Point_3, std::size_t> Tree_point;
166 typedef First_of_pair_property_map< Tree_point > Tree_map;
167 typedef Search_traits_3<K> Traits_base;
168 typedef Search_traits_adapter<Tree_point, Tree_map, Traits_base> Traits;
169 typedef Orthogonal_k_neighbor_search<Traits> Neighbor_search;
170 typedef typename Neighbor_search::Tree Tree;
171
172 // Search tree
173 Tree tree;
174 tree.reserve(cov.size());
175 std::size_t i=0;
176 for (ForwardIterator it = first; it != beyond; ++it, ++i)
177 tree.insert( Tree_point(get(point_map, *it), i) );
178
179 // Convolving
180 ncov.clear();
181 ncov.reserve(cov.size());
182 for (ForwardIterator it = first; it != beyond; ++it) {
183 Neighbor_search search(tree, get(point_map, *it), nb_neighbors_convolve);
184
185 Covariance m;
186 for (typename Neighbor_search::iterator nit = search.begin();
187 nit != search.end();
188 ++nit)
189 {
190 std::size_t index = nit->first.second;
191 for (int i=0; i<6; ++i)
192 m[i] += cov[index][i];
193 }
194
195 ncov.push_back(m);
196 }
197 }
198 /// @endcond
199
200 } // namespace internal
201
202 // ----------------------------------------------------------------------------
203 // Public section
204 // ----------------------------------------------------------------------------
205
206 /**
207 \ingroup PkgPointSetProcessing3Algorithms
208 computes the Voronoi Covariance Measure (VCM) of a point cloud,
209 a construction that can be used for normal estimation and sharp feature detection.
210
211 The VCM associates to each point the covariance matrix of its Voronoi
212 cell intersected with the ball of radius `offset_radius`.
213 In addition, if the second radius `convolution_radius` is positive, the covariance matrices are smoothed
214 via a convolution process. More specifically, each covariance matrix is replaced by
215 the average of the matrices of the points located at a distance at most `convolution_radius`.
216 The choice for parameter `offset_radius` should refer to the geometry of the underlying surface while
217 the choice for parameter `convolution_radius` should refer to the noise level in the point cloud.
218 For example, if the point cloud is a uniform and noise-free sampling of a smooth surface,
219 `offset_radius` should be set to the minimum local feature size of the surface, while `convolution_radius` can be set to zero.
220
221 The Voronoi covariance matrix of each vertex is stored in an array `a` of length 6 and is as follow:
222
223 <center>
224 \f$ \begin{bmatrix}
225 a[0] & a[1] & a[2] \\
226 a[1] & a[3] & a[4] \\
227 a[2] & a[4] & a[5] \\
228 \end{bmatrix}\f$
229 </center>
230
231 \tparam PointRange is a model of `Range`. The value type of
232 its iterator is the key type of the named parameter `point_map`.
233
234 \param points input point range
235 \param ccov output range of covariance matrices.
236 \param offset_radius offset_radius.
237 \param convolution_radius convolution_radius.
238 \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
239
240 \cgalNamedParamsBegin
241 \cgalParamNBegin{point_map}
242 \cgalParamDescription{a property map associating points to the elements of the point set `points`}
243 \cgalParamType{a model of `ReadWritePropertyMap` whose key type is the value type
244 of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
245 \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
246 \cgalParamNEnd
247
248 \cgalParamNBegin{geom_traits}
249 \cgalParamDescription{an instance of a geometric traits class}
250 \cgalParamType{a model of `Kernel`}
251 \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
252 \cgalParamNEnd
253 \cgalNamedParamsEnd
254
255 \sa `CGAL::vcm_is_on_feature_edge()`
256 \sa `CGAL::vcm_estimate_normals()`
257
258 */
259 template <typename PointRange,
260 typename NamedParameters>
261 void
compute_vcm(const PointRange & points,std::vector<std::array<double,6>> & ccov,double offset_radius,double convolution_radius,const NamedParameters & np)262 compute_vcm (const PointRange& points,
263 std::vector< std::array<double, 6> > &ccov,
264 double offset_radius,
265 double convolution_radius,
266 const NamedParameters& np)
267 {
268 using parameters::choose_parameter;
269 using parameters::get_parameter;
270
271 // basic geometric types
272 typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
273 typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
274
275 PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
276 Kernel kernel;
277
278 // First, compute the VCM for each point
279 std::vector< std::array<double, 6> > cov;
280 std::size_t N = 20;
281 internal::vcm_offset (points.begin(), points.end(),
282 point_map,
283 cov,
284 offset_radius,
285 N,
286 kernel);
287 // Then, convolve it (only when convolution_radius != 0)
288 if (convolution_radius == 0) {
289 ccov.reserve(cov.size());
290 std::copy(cov.begin(), cov.end(), std::back_inserter(ccov));
291 } else {
292 internal::vcm_convolve(points.begin(), points.end(),
293 point_map,
294 cov,
295 ccov,
296 convolution_radius,
297 kernel);
298 }
299 }
300
301 /// \cond SKIP_IN_MANUAL
302 // variant with default NP
303 template <typename PointRange>
304 void
compute_vcm(const PointRange & points,std::vector<std::array<double,6>> & ccov,double offset_radius,double convolution_radius)305 compute_vcm (const PointRange& points,
306 std::vector< std::array<double, 6> > &ccov,
307 double offset_radius,
308 double convolution_radius)
309 {
310 compute_vcm (points, ccov, offset_radius, convolution_radius,
311 CGAL::Point_set_processing_3::parameters::all_default (points));
312 }
313
314 /// \endcond
315
316 /// \cond SKIP_IN_MANUAL
317 template <typename PointRange,
318 typename NamedParameters
319 >
320 void
321 vcm_estimate_normals_internal (PointRange& points,
322 double offset_radius, ///< offset radius.
323 double convolution_radius, ///< convolution radius.
324 const NamedParameters& np,
325 int nb_neighbors_convolve = -1 ///< number of neighbors used during the convolution.
326 )
327 {
328 using parameters::choose_parameter;
329 using parameters::get_parameter;
330
331 // basic geometric types
332 typedef typename CGAL::GetPointMap<PointRange, NamedParameters>::type PointMap;
333 typedef typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::type NormalMap;
334 typedef typename Point_set_processing_3::GetK<PointRange, NamedParameters>::Kernel Kernel;
335 typedef typename GetDiagonalizeTraits<NamedParameters, double, 3>::type DiagonalizeTraits;
336
337 CGAL_static_assertion_msg(!(boost::is_same<NormalMap,
338 typename Point_set_processing_3::GetNormalMap<PointRange, NamedParameters>::NoMap>::value),
339 "Error: no normal map");
340
341 PointMap point_map = choose_parameter<PointMap>(get_parameter(np, internal_np::point_map));
342 NormalMap normal_map = choose_parameter<NormalMap>(get_parameter(np, internal_np::normal_map));
343
344 typedef std::array<double, 6> Covariance;
345
346 // Compute the VCM and convolve it
347 std::vector<Covariance> cov;
348 if (nb_neighbors_convolve == -1) {
349 compute_vcm(points,
350 cov,
351 offset_radius,
352 convolution_radius,
353 np);
354 } else {
355 internal::vcm_offset(points.begin(), points.end(),
356 point_map,
357 cov,
358 offset_radius,
359 20,
360 Kernel());
361
362 if (nb_neighbors_convolve > 0)
363 {
364 std::vector<Covariance> ccov;
365 ccov.reserve(cov.size());
366 internal::vcm_convolve(points.begin(), points.end(),
367 point_map,
368 cov,
369 ccov,
370 (unsigned int) nb_neighbors_convolve,
371 Kernel());
372
373 cov.clear();
374 std::copy(ccov.begin(), ccov.end(), std::back_inserter(cov));
375 }
376 }
377
378 // And finally, compute the normals
379 int i = 0;
380 for (typename PointRange::iterator it = points.begin(); it != points.end(); ++it) {
381 std::array<double, 3> enormal = {{ 0,0,0 }};
382 DiagonalizeTraits::extract_largest_eigenvector_of_covariance_matrix
383 (cov[i], enormal);
384
385 typename Kernel::Vector_3 normal(enormal[0],
386 enormal[1],
387 enormal[2]);
388 put(normal_map, *it, normal);
389 i++;
390 }
391 }
392 /// @endcond
393
394
395 /**
396 \ingroup PkgPointSetProcessing3Algorithms
397 Estimates normal directions of the range of `points`
398 using the Voronoi Covariance Measure with a radius for the convolution.
399 The output normals are randomly oriented.
400
401 See `compute_vcm()` for a detailed description of the parameters `offset_radius` and `convolution_radius`
402 and of the Voronoi Covariance Measure.
403
404 \tparam PointRange is a model of `Range`. The value type of
405 its iterator is the key type of the named parameter `point_map`.
406
407 \param points input point range
408 \param offset_radius offset_radius.
409 \param convolution_radius convolution_radius.
410 \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
411
412 \cgalNamedParamsBegin
413 \cgalParamNBegin{point_map}
414 \cgalParamDescription{a property map associating points to the elements of the point set `points`}
415 \cgalParamType{a model of `ReadWritePropertyMap` whose key type is the value type
416 of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
417 \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
418 \cgalParamNEnd
419
420 \cgalParamNBegin{normal_map}
421 \cgalParamDescription{a property map associating normals to the elements of the point set `points`}
422 \cgalParamType{a model of `ReadWritePropertyMap` whose key type is the value type
423 of the iterator of `PointRange` and whose value type is `geom_traits::Vector_3`}
424 \cgalParamNEnd
425
426 \cgalParamNBegin{diagonalize_traits}
427 \cgalParamDescription{the solver used for diagonalizing covariance matrices}
428 \cgalParamType{a class model of `DiagonalizeTraits`}
429 \cgalParamDefault{If Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined
430 then an overload using `Eigen_diagonalize_traits` is provided.
431 Otherwise, the internal implementation `CGAL::Diagonalize_traits` is used}
432 \cgalParamNEnd
433
434 \cgalParamNBegin{geom_traits}
435 \cgalParamDescription{an instance of a geometric traits class}
436 \cgalParamType{a model of `Kernel`}
437 \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
438 \cgalParamNEnd
439 \cgalNamedParamsEnd
440 */
441 template <typename PointRange,
442 typename NamedParameters
443 >
444 void
vcm_estimate_normals(PointRange & points,double offset_radius,double convolution_radius,const NamedParameters & np)445 vcm_estimate_normals (PointRange& points,
446 double offset_radius,
447 double convolution_radius,
448 const NamedParameters& np
449 )
450 {
451 vcm_estimate_normals_internal(points, offset_radius, convolution_radius, np);
452 }
453
454 /// \cond SKIP_IN_MANUAL
455 // variant with default NP
456 template <typename PointRange>
457 void
vcm_estimate_normals(PointRange & points,double offset_radius,double convolution_radius)458 vcm_estimate_normals (PointRange& points,
459 double offset_radius, ///< offset radius.
460 double convolution_radius) ///< convolution radius.
461 {
462 return vcm_estimate_normals
463 (points, offset_radius, convolution_radius,
464 CGAL::Point_set_processing_3::parameters::all_default(points));
465 }
466
467 /// \endcond
468
469
470 /**
471 \ingroup PkgPointSetProcessing3Algorithms
472 Estimates normal directions of the range of `points`
473 using the Voronoi Covariance Measure with a number of neighbors for the convolution.
474 The output normals are randomly oriented.
475
476 See `compute_vcm()` for a detailed description of the parameter `offset_radius`
477 and of the Voronoi Covariance Measure.
478
479 \tparam PointRange is a model of `Range`. The value type of
480 its iterator is the key type of the named parameter `point_map`.
481
482 \param points input point range
483 \param offset_radius offset_radius.
484 \param k number of neighbor points used for convolution.
485 \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
486
487 \cgalNamedParamsBegin
488 \cgalParamNBegin{point_map}
489 \cgalParamDescription{a property map associating points to the elements of the point set `points`}
490 \cgalParamType{a model of `ReadWritePropertyMap` whose key type is the value type
491 of the iterator of `PointRange` and whose value type is `geom_traits::Point_3`}
492 \cgalParamDefault{`CGAL::Identity_property_map<geom_traits::Point_3>`}
493 \cgalParamNEnd
494
495 \cgalParamNBegin{normal_map}
496 \cgalParamDescription{a property map associating normals to the elements of the point set `points`}
497 \cgalParamType{a model of `ReadWritePropertyMap` whose key type is the value type
498 of the iterator of `PointRange` and whose value type is `geom_traits::Vector_3`}
499 \cgalParamNEnd
500
501 \cgalParamNBegin{diagonalize_traits}
502 \cgalParamDescription{the solver used for diagonalizing covariance matrices}
503 \cgalParamType{a class model of `DiagonalizeTraits`}
504 \cgalParamDefault{If Eigen 3 (or greater) is available and `CGAL_EIGEN3_ENABLED` is defined
505 then an overload using `Eigen_diagonalize_traits` is provided.
506 Otherwise, the internal implementation `CGAL::Diagonalize_traits` is used}
507 \cgalParamNEnd
508
509 \cgalParamNBegin{geom_traits}
510 \cgalParamDescription{an instance of a geometric traits class}
511 \cgalParamType{a model of `Kernel`}
512 \cgalParamDefault{a \cgal Kernel deduced from the point type, using `CGAL::Kernel_traits`}
513 \cgalParamNEnd
514 \cgalNamedParamsEnd
515 */
516 template < typename PointRange,
517 typename NamedParameters
518 >
519 void
vcm_estimate_normals(PointRange & points,double offset_radius,unsigned int k,const NamedParameters & np)520 vcm_estimate_normals (PointRange& points,
521 double offset_radius,
522 unsigned int k,
523 const NamedParameters& np
524 )
525 {
526 vcm_estimate_normals_internal(points, offset_radius, 0, np, k);
527 }
528
529 /// \cond SKIP_IN_MANUAL
530 // variant with default NP
531 template <typename PointRange>
532 void
vcm_estimate_normals(PointRange & points,double offset_radius,unsigned int k)533 vcm_estimate_normals (PointRange& points,
534 double offset_radius, ///< offset radius.
535 unsigned int k)
536 {
537 return vcm_estimate_normals
538 (points, offset_radius, k,
539 CGAL::Point_set_processing_3::parameters::all_default(points));
540 }
541
542
543 /// \endcond
544
545 } // namespace CGAL
546
547 #include <CGAL/enable_warnings.h>
548
549 #endif // CGAL_VCM_ESTIMATE_NORMALS_H
550