1 // Copyright (c) 2015 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/Shape_detection/include/CGAL/Regularization/regularize_planes.h $
7 // $Id: regularize_planes.h 0779373 2020-03-26T13:31:46+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Florent Lafarge, Simon Giraudot
12 //
13
14 #ifndef CGAL_REGULARIZE_PLANES_H
15 #define CGAL_REGULARIZE_PLANES_H
16
17 #include <CGAL/license/Shape_detection.h>
18
19 #include <CGAL/centroid.h>
20 #include <CGAL/squared_distance_3.h>
21 #include <CGAL/Shape_detection/Efficient_RANSAC.h>
22
23 namespace CGAL {
24
25 // ----------------------------------------------------------------------------
26 // Private section
27 // ----------------------------------------------------------------------------
28 /// \cond SKIP_IN_MANUAL
29 namespace internal {
30 namespace PlaneRegularization {
31
32 template <typename Traits>
33 struct Plane_cluster
34 {
35 bool is_free;
36 std::vector<std::size_t> planes;
37 std::vector<std::size_t> coplanar_group;
38 std::vector<std::size_t> orthogonal_clusters;
39 typename Traits::Vector_3 normal;
40 typename Traits::FT cosangle_symmetry;
41 typename Traits::FT area;
42 typename Traits::FT cosangle_centroid;
43
Plane_clusterPlane_cluster44 Plane_cluster()
45 : is_free (true)
46 , normal (typename Traits::FT(0.),
47 typename Traits::FT(0.),
48 typename Traits::FT(1.))
49 , cosangle_symmetry (typename Traits::FT(0.))
50 , area (typename Traits::FT(0.))
51 , cosangle_centroid (typename Traits::FT(0.))
52 { }
53 };
54
55
56 template <typename Traits>
regularize_normal(const typename Traits::Vector_3 & n,const typename Traits::Vector_3 & symmetry_direction,typename Traits::FT cos_symmetry)57 typename Traits::Vector_3 regularize_normal
58 (const typename Traits::Vector_3& n,
59 const typename Traits::Vector_3& symmetry_direction,
60 typename Traits::FT cos_symmetry)
61 {
62 typedef typename Traits::FT FT;
63 typedef typename Traits::Point_3 Point;
64 typedef typename Traits::Vector_3 Vector;
65 typedef typename Traits::Line_3 Line;
66 typedef typename Traits::Plane_3 Plane;
67
68 Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction;
69
70 Plane plane_symmetry (pt_symmetry, symmetry_direction);
71 Point pt_normal = CGAL::ORIGIN + n;
72
73 if (n != symmetry_direction || n != -symmetry_direction)
74 {
75 Plane plane_cut (CGAL::ORIGIN, pt_normal, CGAL::ORIGIN + symmetry_direction);
76 Line line;
77 CGAL::Object ob_1 = CGAL::intersection(plane_cut, plane_symmetry);
78 if (!assign(line, ob_1))
79 return n;
80
81 FT delta = std::sqrt ((FT)1. - cos_symmetry * cos_symmetry);
82
83 Point projected_origin = line.projection (CGAL::ORIGIN);
84 Vector line_vector (line);
85 line_vector = line_vector / std::sqrt (line_vector * line_vector);
86 Point pt1 = projected_origin + delta * line_vector;
87 Point pt2 = projected_origin - delta * line_vector;
88
89 if (CGAL::squared_distance (pt_normal, pt1) <= CGAL::squared_distance (pt_normal, pt2))
90 return Vector (CGAL::ORIGIN, pt1);
91 else
92 return Vector (CGAL::ORIGIN, pt2);
93
94 }
95 else
96 return n;
97 }
98
99 template <typename Traits>
regularize_normals_from_prior(const typename Traits::Vector_3 & np,const typename Traits::Vector_3 & n,const typename Traits::Vector_3 & symmetry_direction,typename Traits::FT cos_symmetry)100 typename Traits::Vector_3 regularize_normals_from_prior
101 (const typename Traits::Vector_3& np,
102 const typename Traits::Vector_3& n,
103 const typename Traits::Vector_3& symmetry_direction,
104 typename Traits::FT cos_symmetry)
105 {
106 typedef typename Traits::FT FT;
107 typedef typename Traits::Point_3 Point;
108 typedef typename Traits::Vector_3 Vector;
109 typedef typename Traits::Line_3 Line;
110 typedef typename Traits::Plane_3 Plane;
111
112 Plane plane_orthogonality (CGAL::ORIGIN, np);
113 Point pt_symmetry = CGAL::ORIGIN + cos_symmetry* symmetry_direction;
114
115 Plane plane_symmetry (pt_symmetry, symmetry_direction);
116
117 Line line;
118 CGAL::Object ob_1 = CGAL::intersection (plane_orthogonality, plane_symmetry);
119 if (!assign(line, ob_1))
120 return regularize_normal<Traits> (n, symmetry_direction, cos_symmetry);
121
122 Point projected_origin = line.projection (CGAL::ORIGIN);
123 FT R = CGAL::squared_distance (Point (CGAL::ORIGIN), projected_origin);
124
125 if (R <= 1) // 2 (or 1) possible points intersecting the unit sphere and line
126 {
127 FT delta = std::sqrt ((FT)1. - R);
128 Vector line_vector(line);
129 line_vector = line_vector / std::sqrt (line_vector * line_vector);
130 Point pt1 = projected_origin + delta * line_vector;
131 Point pt2 = projected_origin - delta * line_vector;
132
133 Point pt_n = CGAL::ORIGIN + n;
134 if (CGAL::squared_distance (pt_n, pt1) <= CGAL::squared_distance (pt_n, pt2))
135 return Vector (CGAL::ORIGIN, pt1);
136 else
137 return Vector (CGAL::ORIGIN, pt2);
138 }
139 else //no point intersecting the unit sphere and line
140 return regularize_normal<Traits> (n,symmetry_direction, cos_symmetry);
141
142 }
143
144 template <typename Traits,
145 typename PointRange,
146 typename PointMap,
147 typename IndexMap>
compute_centroids_and_areas(const PointRange & points,PointMap point_map,std::size_t nb_planes,IndexMap index_map,std::vector<typename Traits::Point_3> & centroids,std::vector<typename Traits::FT> & areas)148 void compute_centroids_and_areas (const PointRange& points,
149 PointMap point_map,
150 std::size_t nb_planes,
151 IndexMap index_map,
152 std::vector<typename Traits::Point_3>& centroids,
153 std::vector<typename Traits::FT>& areas)
154 {
155 typedef typename Traits::FT FT;
156 typedef typename Traits::Point_3 Point;
157
158 std::vector<std::vector<Point> > listp (nb_planes);
159
160 for (std::size_t i = 0; i < points.size(); ++ i)
161 {
162 int idx = get (index_map, i);
163 if (idx != -1)
164 listp[std::size_t(idx)].push_back (get(point_map, *(points.begin() + i)));
165 }
166
167 centroids.reserve (nb_planes);
168 areas.reserve (nb_planes);
169 for (std::size_t i = 0; i < nb_planes; ++ i)
170 {
171 centroids.push_back (CGAL::centroid (listp[i].begin (), listp[i].end ()));
172 areas.push_back ((FT)(listp[i].size() / (FT)100.));
173 }
174 }
175
176
177 template <typename Traits,
178 typename PlaneRange,
179 typename PlaneMap>
compute_parallel_clusters(PlaneRange & planes,PlaneMap plane_map,std::vector<Plane_cluster<Traits>> & clusters,std::vector<typename Traits::FT> & areas,typename Traits::FT tolerance_cosangle,const typename Traits::Vector_3 & symmetry_direction)180 void compute_parallel_clusters (PlaneRange& planes,
181 PlaneMap plane_map,
182 std::vector<Plane_cluster<Traits> >& clusters,
183 std::vector<typename Traits::FT>& areas,
184 typename Traits::FT tolerance_cosangle,
185 const typename Traits::Vector_3& symmetry_direction)
186 {
187
188 typedef typename Traits::FT FT;
189 typedef typename Traits::Vector_3 Vector;
190
191 // find pairs of epsilon-parallel primitives and store them in parallel_planes
192 std::vector<std::vector<std::size_t> > parallel_planes (planes.size ());
193 for (std::size_t i = 0; i < std::size_t(planes.size ()); ++ i)
194 {
195 typename PlaneRange::iterator it = planes.begin() + i;
196 Vector v1 = get(plane_map, *it).orthogonal_vector();
197
198 for (std::size_t j = 0; j < std::size_t(planes.size()); ++ j)
199 {
200 if (i == j)
201 continue;
202
203 typename PlaneRange::iterator it2 = planes.begin() + j;
204 Vector v2 = get(plane_map, *it2).orthogonal_vector();
205
206 if (std::fabs (v1 * v2) > 1. - tolerance_cosangle)
207 parallel_planes[i].push_back (j);
208 }
209 }
210
211
212 std::vector<bool> is_available (planes.size (), true);
213
214 for (std::size_t i = 0; i < std::size_t(planes.size()); ++ i)
215 {
216
217 if(is_available[i])
218 {
219 const typename Traits::Plane_3& plane = get(plane_map, *(planes.begin() + i));
220
221 is_available[i] = false;
222
223 clusters.push_back (Plane_cluster<Traits>());
224 Plane_cluster<Traits>& clu = clusters.back ();
225
226 //initialization containers
227 clu.planes.push_back (i);
228
229 std::vector<std::size_t> index_container_former_ring_parallel;
230 index_container_former_ring_parallel.push_back(i);
231
232 std::list<std::size_t> index_container_current_ring_parallel;
233
234 //propagation over the pairs of epsilon-parallel primitives
235 bool propagation=true;
236 clu.normal = plane.orthogonal_vector ();
237 clu.area = areas[i];
238
239 do
240 {
241 propagation = false;
242
243 for (std::size_t k = 0; k < index_container_former_ring_parallel.size(); ++ k)
244 {
245
246 std::size_t plane_index = index_container_former_ring_parallel[k];
247
248 for (std::size_t l = 0; l < parallel_planes[plane_index].size(); ++ l)
249 {
250 std::size_t it = parallel_planes[plane_index][l];
251
252 Vector normal_it = get(plane_map, *(planes.begin() + it)).orthogonal_vector ();
253
254 if(is_available[it]
255 && std::fabs (normal_it*clu.normal) > 1. - tolerance_cosangle )
256 {
257 propagation = true;
258 index_container_current_ring_parallel.push_back(it);
259 is_available[it]=false;
260
261 if(clu.normal * normal_it <0)
262 normal_it = -normal_it;
263
264 clu.normal = (FT)clu.area * clu.normal
265 + (FT)areas[it] * normal_it;
266 FT norm = (FT)1. / std::sqrt (clu.normal.squared_length());
267 clu.normal = norm * clu.normal;
268 clu.area += areas[it];
269 }
270 }
271 }
272
273 //update containers
274 index_container_former_ring_parallel.clear();
275 for (std::list<std::size_t>::iterator it = index_container_current_ring_parallel.begin();
276 it != index_container_current_ring_parallel.end(); ++it)
277 {
278 index_container_former_ring_parallel.push_back(*it);
279 clu.planes.push_back(*it);
280 }
281 index_container_current_ring_parallel.clear();
282
283 }
284 while(propagation);
285
286 if (symmetry_direction != CGAL::NULL_VECTOR)
287 {
288 clu.cosangle_symmetry = symmetry_direction * clu.normal;
289 if (clu.cosangle_symmetry < 0.)
290 {
291 clu.normal = -clu.normal;
292 clu.cosangle_symmetry = -clu.cosangle_symmetry;
293 }
294 }
295 }
296 }
297
298 is_available.clear();
299 }
300
301 template <typename Traits>
cluster_symmetric_cosangles(std::vector<Plane_cluster<Traits>> & clusters,typename Traits::FT tolerance_cosangle,typename Traits::FT tolerance_cosangle_ortho)302 void cluster_symmetric_cosangles (std::vector<Plane_cluster<Traits> >& clusters,
303 typename Traits::FT tolerance_cosangle,
304 typename Traits::FT tolerance_cosangle_ortho)
305 {
306 typedef typename Traits::FT FT;
307
308 std::vector < FT > cosangle_centroids;
309 std::vector < std::size_t> list_cluster_index;
310 for( std::size_t i = 0; i < clusters.size(); ++ i)
311 list_cluster_index.push_back(static_cast<std::size_t>(-1));
312
313 std::size_t mean_index = 0;
314 for (std::size_t i = 0; i < clusters.size(); ++ i)
315 {
316 if(list_cluster_index[i] == static_cast<std::size_t>(-1))
317 {
318 list_cluster_index[i] = mean_index;
319 FT mean = clusters[i].area * clusters[i].cosangle_symmetry;
320 FT mean_area = clusters[i].area;
321
322 for (std::size_t j = i+1; j < clusters.size(); ++ j)
323 {
324 if (list_cluster_index[j] == static_cast<std::size_t>(-1)
325 && std::fabs (clusters[j].cosangle_symmetry -
326 mean / mean_area) < tolerance_cosangle_ortho)
327 {
328 list_cluster_index[j] = mean_index;
329 mean_area += clusters[j].area;
330 mean += clusters[j].area * clusters[j].cosangle_symmetry;
331 }
332 }
333 ++ mean_index;
334 mean /= mean_area;
335 cosangle_centroids.push_back (mean);
336 }
337 }
338
339 for (std::size_t i = 0; i < cosangle_centroids.size(); ++ i)
340 {
341 if (cosangle_centroids[i] < tolerance_cosangle_ortho)
342 cosangle_centroids[i] = 0;
343 else if (cosangle_centroids[i] > 1. - tolerance_cosangle)
344 cosangle_centroids[i] = 1;
345 }
346 for (std::size_t i = 0; i < clusters.size(); ++ i)
347 clusters[i].cosangle_symmetry = cosangle_centroids[list_cluster_index[i]];
348 }
349
350
351 template <typename Traits>
subgraph_mutually_orthogonal_clusters(std::vector<Plane_cluster<Traits>> & clusters,const typename Traits::Vector_3 & symmetry_direction)352 void subgraph_mutually_orthogonal_clusters (std::vector<Plane_cluster<Traits> >& clusters,
353 const typename Traits::Vector_3& symmetry_direction)
354 {
355 typedef typename Traits::FT FT;
356 typedef typename Traits::Vector_3 Vector;
357
358 std::vector < std::vector < std::size_t> > subgraph_clusters;
359 std::vector < std::size_t> subgraph_clusters_max_area_index;
360
361 for (std::size_t i = 0; i < clusters.size(); ++ i)
362 clusters[i].is_free = true;
363
364 for (std::size_t i = 0; i < clusters.size(); ++ i)
365 {
366 if(clusters[i].is_free)
367 {
368 clusters[i].is_free = false;
369 FT max_area = clusters[i].area;
370 std::size_t index_max_area = i;
371
372 //initialization containers
373 std::vector < std::size_t > index_container;
374 index_container.push_back(i);
375 std::vector < std::size_t > index_container_former_ring;
376 index_container_former_ring.push_back(i);
377 std::list < std::size_t > index_container_current_ring;
378
379 //propagation
380 bool propagation=true;
381 do
382 {
383 propagation=false;
384
385 //neighbors
386 for (std::size_t k=0;k<index_container_former_ring.size();k++)
387 {
388
389 std::size_t cluster_index=index_container_former_ring[k];
390
391 for (std::size_t j = 0; j < clusters[cluster_index].orthogonal_clusters.size(); ++ j)
392 {
393 std::size_t cluster_index_2 = clusters[cluster_index].orthogonal_clusters[j];
394 if(clusters[cluster_index_2].is_free)
395 {
396 propagation = true;
397 index_container_current_ring.push_back(cluster_index_2);
398 clusters[cluster_index_2].is_free = false;
399
400 if(max_area < clusters[cluster_index_2].area)
401 {
402 max_area = clusters[cluster_index_2].area;
403 index_max_area = cluster_index_2;
404 }
405 }
406 }
407 }
408
409 //update containers
410 index_container_former_ring.clear();
411 for(std::list < std::size_t>::iterator it = index_container_current_ring.begin();
412 it != index_container_current_ring.end(); ++it)
413 {
414 index_container_former_ring.push_back(*it);
415 index_container.push_back(*it);
416 }
417 index_container_current_ring.clear();
418
419 }
420 while(propagation);
421 subgraph_clusters.push_back(index_container);
422 subgraph_clusters_max_area_index.push_back(index_max_area);
423 }
424 }
425
426 //create subgraphs of mutually orthogonal clusters in which the
427 //largest cluster is excluded and store in
428 //subgraph_clusters_prop
429 std::vector < std::vector < std::size_t> > subgraph_clusters_prop;
430 for (std::size_t i=0;i<subgraph_clusters.size(); i++)
431 {
432 std::size_t index=subgraph_clusters_max_area_index[i];
433 std::vector < std::size_t> subgraph_clusters_prop_temp;
434 for (std::size_t j=0;j<subgraph_clusters[i].size(); j++)
435 if(subgraph_clusters[i][j]!=index)
436 subgraph_clusters_prop_temp.push_back(subgraph_clusters[i][j]);
437
438 subgraph_clusters_prop.push_back(subgraph_clusters_prop_temp);
439 }
440
441 //regularization of cluster normals : in each subgraph, we start
442 //from the largest area cluster and we propagate over the subgraph
443 //by regularizing the normals of the clusters according to
444 //orthogonality and cos angle to symmetry direction
445
446 for (std::size_t i = 0; i < clusters.size(); ++ i)
447 clusters[i].is_free = true;
448
449 for (std::size_t i = 0; i < subgraph_clusters_prop.size(); ++ i)
450 {
451
452 std::size_t index_current=subgraph_clusters_max_area_index[i];
453
454 Vector vec_current=regularize_normal<Traits>
455 (clusters[index_current].normal,
456 symmetry_direction,
457 clusters[index_current].cosangle_symmetry);
458 clusters[index_current].normal = vec_current;
459 clusters[index_current].is_free = false;
460
461 //initialization containers
462 std::vector < std::size_t> index_container;
463 index_container.push_back(index_current);
464 std::vector < std::size_t> index_container_former_ring;
465 index_container_former_ring.push_back(index_current);
466 std::list < std::size_t> index_container_current_ring;
467
468 //propagation
469 bool propagation=true;
470 do
471 {
472 propagation=false;
473
474 //neighbors
475 for (std::size_t k=0;k<index_container_former_ring.size();k++)
476 {
477
478 std::size_t cluster_index=index_container_former_ring[k];
479
480 for (std::size_t j = 0; j < clusters[cluster_index].orthogonal_clusters.size(); ++ j)
481 {
482 std::size_t cluster_index_2 = clusters[cluster_index].orthogonal_clusters[j];
483 if(clusters[cluster_index_2].is_free)
484 {
485 propagation = true;
486 index_container_current_ring.push_back(cluster_index_2);
487 clusters[cluster_index_2].is_free = false;
488
489 Vector new_vect=regularize_normals_from_prior<Traits>
490 (clusters[cluster_index].normal,
491 clusters[cluster_index_2].normal,
492 symmetry_direction,
493 clusters[cluster_index_2].cosangle_symmetry);
494 clusters[cluster_index_2].normal = new_vect;
495 }
496 }
497 }
498
499 //update containers
500 index_container_former_ring.clear();
501 for(std::list < std::size_t>::iterator it = index_container_current_ring.begin();
502 it != index_container_current_ring.end(); ++it)
503 {
504 index_container_former_ring.push_back(*it);
505 index_container.push_back(*it);
506 }
507 index_container_current_ring.clear();
508 }while(propagation);
509 }
510 }
511
512
513 } // namespace PlaneRegularization
514 } // namespace internal
515 /// \endcond
516
517
518 // ----------------------------------------------------------------------------
519 // Public section
520 // ----------------------------------------------------------------------------
521
522 /// \ingroup PkgShapeDetectionRef
523
524 /*!
525 Given a set of detected planes with their corresponding inlier sets,
526 this function enables to regularize the planes:
527
528 - Planes near parallel can be made exactly parallel;
529
530 - Planes near orthogonal can be made exactly orthogonal;
531
532 - Planes parallel and near coplanar can be made exactly coplanar;
533
534 - Planes near symmetrical with a user-defined axis can be made
535 exactly symmetrical.
536
537 Planes are directly modified. Points are left unaltered, as well as
538 their relationships to planes (no transfer of point from a primitive
539 plane to another).
540
541 The implementation follows \cgalCite{cgal:vla-lod-15}.
542
543 \tparam PointRange must be a model of `ConstRange` with points.
544
545 \tparam PointPMap must be a model of `ReadablePropertyMap` with value type `Kernel::Point_3`.
546 It can be omitted if the value type of the iterator of `PointRange` is convertible to `Point_3<Kernel>`.
547
548 \tparam PlaneRange must be a model of `Range` with planes.
549
550 \tparam PlaneMap must be a model of `WritablePropertyMap` with value type `Kernel::Plane_3`.
551 It can be omitted if the value type of the iterator of `PlaneRange` is convertible to `Plane_3<Kernel>`.
552
553 \tparam IndexMap must be a model of `ReadablePropertyMap` with value type `int`.
554
555 \tparam Kernel must be a geometric traits class.
556 It can be omitted and deduced automatically from the value type of `PointMap`.
557
558 \param points `ConstRange` of points
559 \param point_map property map: value_type of `typename PointRange::const_iterator` -> `Point_3`
560 \param planes range of planes
561 \param plane_map property map: value_type of `typename PlaneRange::iterator` -> `Plane_3`
562 \param index_map property map: index of point `std::size_t` -> index of plane `int` (-1 if the point is not assigned to a plane)
563
564 \param regularize_parallelism select whether parallelism is
565 regularized or not
566
567 \param regularize_orthogonality select whether orthogonality is
568 regularized or not
569
570 \param regularize_coplanarity select whether coplanarity is
571 regularized or not
572
573 \param regularize_axis_symmetry select whether axis symmetry is
574 regularized or not
575
576 \param tolerance_angle tolerance of deviation between normal
577 vectors of planes (in degrees) used for parallelism, orthogonality,
578 and axis symmetry. %Default value is 25 degrees.
579
580 \param tolerance_coplanarity maximal distance between two parallel
581 planes such that they are considered coplanar. %Default value is
582 0.01.
583
584 \param symmetry_direction chosen axis for symmetry
585 regularization. %Default value is the Z axis.
586 */
587
588 // This variant requires all parameters
589 template <typename PointRange,
590 typename PointMap,
591 typename PlaneRange,
592 typename PlaneMap,
593 typename IndexMap,
594 typename Kernel>
595 void regularize_planes(const PointRange& points,
596 PointMap point_map,
597 PlaneRange& planes,
598 PlaneMap plane_map,
599 IndexMap index_map,
600 const Kernel&,
601 bool regularize_parallelism,
602 bool regularize_orthogonality,
603 bool regularize_coplanarity,
604 bool regularize_axis_symmetry,
605 double tolerance_angle = 25.0,
606 double tolerance_coplanarity = 0.01,
607 typename Kernel::Vector_3 symmetry_direction
608 = typename Kernel::Vector_3 (0.0, 0.0, 1.0))
609 {
610 typedef typename Kernel::FT FT;
611 typedef typename Kernel::Point_3 Point;
612 typedef typename Kernel::Vector_3 Vector;
613 typedef typename Kernel::Plane_3 Plane;
614
615 typedef typename internal::PlaneRegularization::Plane_cluster<Kernel>
616 Plane_cluster;
617
618 /*
619 * Compute centroids and areas
620 */
621 std::vector<Point> centroids;
622 std::vector<FT> areas;
623 internal::PlaneRegularization::compute_centroids_and_areas<Kernel>
624 (points, point_map, planes.size(), index_map, centroids, areas);
625
626 tolerance_angle = tolerance_angle * (FT)CGAL_PI / (FT)(180);
627 FT tolerance_cosangle = (FT)(1. - std::cos (tolerance_angle));
628 FT tolerance_cosangle_ortho = (FT)(std::cos ((FT)0.5 * (FT)CGAL_PI - (FT)tolerance_angle));
629
630 // clustering the parallel primitives and store them in clusters
631 // & compute the normal, size and cos angle to the symmetry
632 // direction of each cluster
633 std::vector<Plane_cluster> clusters;
634 internal::PlaneRegularization::compute_parallel_clusters<Kernel>
635 (planes, plane_map, clusters, areas,
636 (regularize_parallelism ? tolerance_cosangle : (FT)0.0),
637 (regularize_axis_symmetry ? symmetry_direction : CGAL::NULL_VECTOR));
638
639 if (regularize_orthogonality)
640 {
641 //discovery orthogonal relationship between clusters
642 for (std::size_t i = 0; i < clusters.size(); ++ i)
643 {
644 for (std::size_t j = i + 1; j < clusters.size(); ++ j)
645 {
646 if (std::fabs (clusters[i].normal * clusters[j].normal) < tolerance_cosangle_ortho)
647 {
648 clusters[i].orthogonal_clusters.push_back (j);
649 clusters[j].orthogonal_clusters.push_back (i);
650 }
651 }
652 }
653 }
654
655 if (regularize_axis_symmetry)
656 {
657 //clustering the symmetry cos angle and store their centroids in
658 //cosangle_centroids and the centroid index of each cluster in
659 //list_cluster_index
660 internal::PlaneRegularization::cluster_symmetric_cosangles<Kernel>
661 (clusters, tolerance_cosangle, tolerance_cosangle_ortho);
662 }
663
664 //find subgraphs of mutually orthogonal clusters (store index of
665 //clusters in subgraph_clusters), and select the cluster of
666 //largest area
667 if (regularize_orthogonality || regularize_axis_symmetry)
668 internal::PlaneRegularization::subgraph_mutually_orthogonal_clusters<Kernel>
669 (clusters, (regularize_axis_symmetry ? symmetry_direction : CGAL::NULL_VECTOR));
670
671 //recompute optimal plane for each primitive after normal regularization
672 for (std::size_t i=0; i < clusters.size(); ++ i)
673 {
674 Vector vec_reg = clusters[i].normal;
675 for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j)
676 {
677 std::size_t index_prim = clusters[i].planes[j];
678 const Plane& plane = get(plane_map, *(planes.begin() + index_prim));
679
680 Point pt_reg = plane.projection (centroids[index_prim]);
681 if(plane.orthogonal_vector () * vec_reg < 0)
682 vec_reg=-vec_reg;
683 Plane plane_reg(pt_reg,vec_reg);
684
685 if( std::fabs(plane.orthogonal_vector () * vec_reg) > 1. - tolerance_cosangle)
686 put(plane_map, *(planes.begin() + index_prim), plane_reg);
687 }
688 }
689
690 if (regularize_coplanarity)
691 {
692 //detecting co-planarity and store in list_coplanar_prim
693 for (std::size_t i = 0; i < clusters.size(); ++ i)
694 {
695 Vector vec_reg = clusters[i].normal;
696
697 for (std::size_t ip = 0; ip < clusters[i].planes.size(); ++ ip)
698 clusters[i].coplanar_group.push_back (static_cast<std::size_t>(-1));
699
700 std::size_t cop_index=0;
701
702 for (std::size_t j = 0; j < clusters[i].planes.size(); ++ j)
703 {
704 std::size_t index_prim = clusters[i].planes[j];
705
706 if (clusters[i].coplanar_group[j] == static_cast<std::size_t>(-1))
707 {
708 const Plane& plane = get(plane_map, *(planes.begin() + index_prim));
709
710 clusters[i].coplanar_group[j] = cop_index;
711
712 Point pt_reg = plane.projection(centroids[index_prim]);
713 Plane plan_reg(pt_reg,vec_reg);
714
715 for (std::size_t k = j + 1; k < clusters[i].planes.size(); ++ k)
716 {
717 if (clusters[i].coplanar_group[k] == static_cast<std::size_t>(-1))
718 {
719 std::size_t index_prim_next = clusters[i].planes[k];
720 const Plane& plane_next = get(plane_map, *(planes.begin() + index_prim_next));
721 Point pt_reg_next = plane_next.projection(centroids[index_prim_next]);
722 Point pt_proj=plan_reg.projection(pt_reg_next);
723 FT distance = std::sqrt (CGAL::squared_distance(pt_reg_next,pt_proj));
724
725 if (distance < tolerance_coplanarity)
726 clusters[i].coplanar_group[k] = cop_index;
727 }
728 }
729 cop_index++;
730 }
731 }
732 //regularize primitive position by computing barycenter of coplanar planes
733 std::vector<Point> pt_bary (cop_index, Point ((FT)0., (FT)0., (FT)0.));
734 std::vector<FT> area (cop_index, 0.);
735
736 for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j)
737 {
738 std::size_t index_prim = clusters[i].planes[j];
739 std::size_t group = clusters[i].coplanar_group[j];
740
741 Point pt_reg = get(plane_map, *(planes.begin()+index_prim)).projection(centroids[index_prim]);
742
743 pt_bary[group] = CGAL::barycenter (pt_bary[group], area[group], pt_reg, areas[index_prim]);
744 area[group] += areas[index_prim];
745 }
746
747
748 for (std::size_t j = 0; j < clusters[i].planes.size (); ++ j)
749 {
750 std::size_t index_prim = clusters[i].planes[j];
751 std::size_t group = clusters[i].coplanar_group[j];
752 Plane plane_reg (pt_bary[group], vec_reg);
753
754 if (get(plane_map, *(planes.begin() + index_prim)).orthogonal_vector ()
755 * plane_reg.orthogonal_vector() < 0)
756 put(plane_map, *(planes.begin() + index_prim), plane_reg.opposite());
757 else
758 put(plane_map, *(planes.begin() + index_prim), plane_reg);
759 }
760 }
761 }
762 }
763
764
765 /// \cond SKIP_IN_MANUAL
766
767 // Workaround for bug reported here:
768 // https://developercommunity.visualstudio.com/content/problem/340310/unaccepted-typename-that-other-compilers-require.html
769 #if _MSC_VER == 1915
770 #define CGAL_TYPENAME_FOR_MSC
771 #else
772 #define CGAL_TYPENAME_FOR_MSC typename
773 #endif
774
775 // This variant deduces the kernel from the point property map.
776 template <typename PointRange,
777 typename PointMap,
778 typename PlaneRange,
779 typename PlaneMap,
780 typename IndexMap>
781 void regularize_planes (const PointRange& points,
782 PointMap point_map,
783 PlaneRange& planes,
784 PlaneMap plane_map,
785 IndexMap index_map,
786 bool regularize_parallelism,
787 bool regularize_orthogonality,
788 bool regularize_coplanarity,
789 bool regularize_axis_symmetry,
790 double tolerance_angle = 25.0,
791 double tolerance_coplanarity = 0.01,
792 typename Kernel_traits
793 <typename boost::property_traits
794 <PointMap>::value_type>::Kernel::Vector_3 symmetry_direction
795 = CGAL_TYPENAME_FOR_MSC Kernel_traits
796 <typename boost::property_traits
797 <PointMap>::value_type>::Kernel::Vector_3
798 (
799 CGAL_TYPENAME_FOR_MSC Kernel_traits
800 <typename boost::property_traits
801 <PointMap>::value_type>::Kernel::FT(0.),
802 CGAL_TYPENAME_FOR_MSC Kernel_traits
803 <typename boost::property_traits
804 <PointMap>::value_type>::Kernel::FT(0.),
805 CGAL_TYPENAME_FOR_MSC Kernel_traits
806 <typename boost::property_traits
807 <PointMap>::value_type>::Kernel::FT(1.)))
808 {
809 typedef typename boost::property_traits<PointMap>::value_type Point;
810 typedef typename Kernel_traits<Point>::Kernel Kernel;
811
812 regularize_planes (points, point_map, planes, plane_map, index_map, Kernel(),
813 regularize_parallelism, regularize_orthogonality,
814 regularize_coplanarity, regularize_axis_symmetry,
815 tolerance_angle, tolerance_coplanarity, symmetry_direction);
816 }
817
818 #ifdef CGAL_TYPENAME_FOR_MSC
819 #undef CGAL_TYPENAME_FOR_MSC
820 #endif
821
822
823 /// \endcond
824
825 } // namespace CGAL
826
827 #endif // CGAL_REGULARIZE_PLANES_H
828