1 // Copyright (c) 2001,2011 Max-Planck-Institute Saarbruecken (Germany).
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/Convex_hull_3/include/CGAL/convex_hull_3.h $
7 // $Id: convex_hull_3.h 61d0fb5 2021-01-05T12:06:20+01:00 Sébastien Loriot
8 // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
9 //
10 //
11 // Author(s) : Susan Hert <hert@mpi-sb.mpg.de>
12 // : Amol Prakash <prakash@mpi-sb.mpg.de>
13 // : Andreas Fabri
14
15 #ifndef CGAL_CONVEX_HULL_3_H
16 #define CGAL_CONVEX_HULL_3_H
17
18 #include <CGAL/license/Convex_hull_3.h>
19
20 #include <CGAL/disable_warnings.h>
21
22 #include <CGAL/basic.h>
23 #include <CGAL/algorithm.h>
24 #include <CGAL/convex_hull_2.h>
25 #include <CGAL/Projection_traits_xy_3.h>
26 #include <CGAL/Projection_traits_xz_3.h>
27 #include <CGAL/Projection_traits_yz_3.h>
28 #include <CGAL/Convex_hull_traits_3.h>
29 #include <CGAL/Convex_hull_2/ch_assertions.h>
30 #include <CGAL/Triangulation_data_structure_2.h>
31 #include <CGAL/Triangulation_vertex_base_with_info_2.h>
32 #include <CGAL/Cartesian_converter.h>
33 #include <CGAL/Simple_cartesian.h>
34
35 #include <CGAL/internal/Exact_type_selector.h>
36 #include <CGAL/boost/graph/copy_face_graph.h>
37 #include <CGAL/boost/graph/Named_function_parameters.h>
38 #include <CGAL/boost/graph/graph_traits_Triangulation_data_structure_2.h>
39 #include <CGAL/boost/graph/properties_Triangulation_data_structure_2.h>
40 #include <CGAL/Polyhedron_3_fwd.h>
41 #include <CGAL/boost/graph/Euler_operations.h>
42 #include <CGAL/boost/iterator/transform_iterator.hpp>
43 #include <CGAL/boost/graph/named_params_helper.h>
44 #include <CGAL/is_iterator.h>
45
46 #include <boost/next_prior.hpp>
47 #include <boost/type_traits/is_floating_point.hpp>
48 #include <boost/type_traits/is_same.hpp>
49 #include <boost/mpl/has_xxx.hpp>
50 #include <boost/graph/graph_traits.hpp>
51
52 #ifndef CGAL_CH_NO_POSTCONDITIONS
53 #include <CGAL/convexity_check_3.h>
54 #endif // CGAL_CH_NO_POSTCONDITIONS
55
56 #include <algorithm>
57 #include <iostream>
58 #include <list>
59 #include <memory>
60 #include <vector>
61 #include <type_traits>
62 #include <utility>
63
64 // first some internal stuff to avoid using a true Face_graph model for extreme_points_3
65 namespace CGAL {
66
67 // Forward declaration
68 template<class VertexPointMap,class Base_traits> class Extreme_points_traits_adapter_3;
69
70 namespace Convex_hull_3 {
71 namespace internal {
72
73 // wrapper used as a MutableFaceGraph to extract extreme points
74 template <class OutputIterator>
75 struct Output_iterator_wrapper
76 {
77 OutputIterator out;
Output_iterator_wrapperOutput_iterator_wrapper78 Output_iterator_wrapper(OutputIterator out)
79 : out(out)
80 {}
81 };
82
83 template <class Point_3, class OutputIterator>
add_isolated_points(const Point_3 & point,Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> & w)84 void add_isolated_points(const Point_3& point,
85 Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator>& w)
86 {
87 *w.out++ = point;
88 }
89
90 template <class Point_3, class PolygonMesh>
add_isolated_points(const Point_3 & point,PolygonMesh & P)91 void add_isolated_points(const Point_3& point, PolygonMesh& P)
92 {
93 put(get(CGAL::vertex_point, P), add_vertex(P), point);
94 }
95
96 template <class Point_3, class OutputIterator>
copy_ch2_to_face_graph(const std::list<Point_3> & CH_2,Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> & w)97 void copy_ch2_to_face_graph(const std::list<Point_3>& CH_2,
98 Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator>& w)
99 {
100 for(const Point_3& p : CH_2)
101 *w.out++ = p;
102 }
103
104 } // namespace internal
105 } // namespace Convex_hull_3
106
107 template <class TDS, class OutputIterator>
copy_face_graph(const TDS & tds,Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> & wrapper)108 void copy_face_graph(const TDS& tds, Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator>& wrapper)
109 {
110 typedef typename boost::graph_traits<TDS>::vertex_descriptor vertex_descriptor;
111 typename boost::property_map<TDS, boost::vertex_point_t >::const_type vpm = get(boost::vertex_point, tds);
112 for(vertex_descriptor vh : vertices(tds))
113 {
114 *wrapper.out++ = get(vpm, vh);
115 }
116 }
117
118 template <class OutputIterator>
clear(Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> &)119 void clear(Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator>&)
120 {}
121
122 template <class Point, class OutputIterator>
make_tetrahedron(const Point & p0,const Point & p1,const Point & p2,const Point & p3,Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> & w)123 void make_tetrahedron(const Point& p0, const Point& p1, const Point& p2, const Point& p3,
124 Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator>& w)
125 {
126 *w.out++ = p0;
127 *w.out++ = p1;
128 *w.out++ = p2;
129 *w.out++ = p3;
130 }
131
132 } // CGAL
133
134 namespace boost {
135
136 // needed so that the regular make_tetrahedron of CGAL does not complain when tried to be instantiated
137 template <class OutputIterator>
138 struct graph_traits<CGAL::Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> >
139 {
140 typedef void* halfedge_descriptor;
141 };
142
143 } // namespace boost
144
145 namespace CGAL {
146 namespace Convex_hull_3 {
147 namespace internal {
148
149 //struct to select the default traits class for computing convex hull
150 template< class Point_3,
151 class PolygonMesh = Default,
152 class Is_floating_point=typename boost::is_floating_point<typename Kernel_traits<Point_3>::Kernel::FT>::type,
153 class Has_filtered_predicates_tag=typename Kernel_traits<Point_3>::Kernel::Has_filtered_predicates_tag >
154 struct Default_traits_for_Chull_3{
155 typedef typename Kernel_traits<Point_3>::Kernel type;
156 };
157
158 //FT is a floating point type and Kernel is a filtered kernel
159 template <class Point_3, class PolygonMesh>
160 struct Default_traits_for_Chull_3<Point_3, PolygonMesh, boost::true_type,Tag_true>{
161 typedef Convex_hull_traits_3< typename Kernel_traits<Point_3>::Kernel, PolygonMesh, Tag_true > type;
162 };
163
164 template <class Traits>
165 struct Default_polyhedron_for_Chull_3{
166 typedef CGAL::Polyhedron_3<Traits> type;
167 };
168
169 template <class K, class P, class Tag>
170 struct Default_polyhedron_for_Chull_3<Convex_hull_traits_3<K, P, Tag> >{
171 typedef typename Convex_hull_traits_3<K, P, Tag>::Polygon_mesh type;
172 };
173
174 template <class T>
175 struct Is_cartesian_kernel
176 {
177 typedef boost::false_type type;
178 };
179
180 template <class Kernel, class PolygonMesh>
181 struct Is_cartesian_kernel< Convex_hull_traits_3<Kernel, PolygonMesh, Tag_true> >
182 {
183 // Rational here is that Tag_true can only be passed by us since it is not documented
184 // so we can assume that Kernel is a CGAL Kernel
185 typedef typename boost::is_same<typename Kernel::Kernel_tag, Cartesian_tag>::type type;
186 };
187
188 // Predicate internally used as a wrapper around has_on_positive_side
189 // We provide a partial specialization restricted to the case of CGAL Cartesian Kernels with inexact constructions below
190 //template <class Traits,class Tag_use_advanced_filtering=typename Use_advanced_filtering<Traits>::type >
191 template <class Traits, class Is_CK = typename Is_cartesian_kernel<Traits>::type >
192 class Is_on_positive_side_of_plane_3{
193 typedef typename Traits::Point_3 Point_3;
194 typename Traits::Plane_3 plane;
195 typename Traits::Has_on_positive_side_3 has_on_positive_side;
196 public:
197 typedef Protect_FPU_rounding<false> Protector;
198
199 Is_on_positive_side_of_plane_3(const Traits& traits,const Point_3& p,const Point_3& q,const Point_3& r)
200 :plane(traits.construct_plane_3_object()(p,q,r)),has_on_positive_side(traits.has_on_positive_side_3_object()) {}
201
202 bool operator() (const Point_3& s) const
203 {
204 return has_on_positive_side(plane,s);
205 }
206 };
207
208 template <class Base_traits, class VPM, class Is_CK>
209 class Is_on_positive_side_of_plane_3< Extreme_points_traits_adapter_3<VPM,Base_traits>, Is_CK>
210 : public Is_on_positive_side_of_plane_3< Base_traits >
211 {
212 typedef Extreme_points_traits_adapter_3<VPM, Base_traits> Traits;
213 typedef Is_on_positive_side_of_plane_3< Base_traits > Base;
214 typedef typename Traits::Point_3 Point_3;
215 const Traits& m_traits;
216 public:
217 typedef typename Base::Protector Protector;
218
219 Is_on_positive_side_of_plane_3(const Traits& traits,
220 const Point_3& p,const Point_3& q,const Point_3& r)
221 : Base(static_cast<const Base_traits&>(traits), traits.get_point(p), traits.get_point(q), traits.get_point(r))
222 , m_traits(traits)
223 {}
224
225 bool operator() (const Point_3& s) const
226 {
227 return static_cast<const Base*>(this)->operator()(m_traits.get_point(s));
228 }
229 };
230
231 //This predicate uses copy of the code from the statically filtered version of
232 //Orientation_3. The rational is that the plane is a member of the functor
233 //so optimization are done to avoid doing several time operations on the plane.
234 //The main operator() first tries the static version of the predicate, then uses
235 //interval arithmetic (the protector must be created before using this predicate)
236 //and in case of failure, exact arithmetic is used.
237 template <class Kernel, class P>
238 class Is_on_positive_side_of_plane_3<Convex_hull_traits_3<Kernel, P, Tag_true>, boost::true_type >{
239 typedef Simple_cartesian<CGAL::internal::Exact_field_selector<double>::Type> Exact_K;
240 typedef Simple_cartesian<Interval_nt_advanced > Approx_K;
241 typedef Convex_hull_traits_3<Kernel, P, Tag_true> Traits;
242 typedef typename Traits::Point_3 Point_3;
243
244 Cartesian_converter<Kernel,Approx_K> to_AK;
245 Cartesian_converter<Kernel,Exact_K> to_EK;
246
247 template <typename K>
248 struct Vector_plus_point {
249 typename K::Vector_3 vector;
250 typename K::Point_3 point;
251 };
252
253 const Point_3& p,q,r;
254 mutable Vector_plus_point<Approx_K> ak_plane;
255 mutable Vector_plus_point<Exact_K>* ek_plane_ptr;
256
257 double m10,m20,m21,Maxx,Maxy,Maxz;
258
259 static const int STATIC_FILTER_FAILURE = 555;
260
261 //this function is a made from the statically filtered version of Orientation_3
262 int static_filtered(double psx,double psy, double psz) const{
263
264 // Then semi-static filter.
265 double apsx = CGAL::abs(psx);
266 double apsy = CGAL::abs(psy);
267 double apsz = CGAL::abs(psz);
268
269 double maxx = (Maxx < apsx)? apsx : Maxx;
270 double maxy = (Maxy < apsy)? apsy : Maxy;
271 double maxz = (Maxz < apsz)? apsz : Maxz;
272
273 double det = psx*m10 - m20*psy + m21*psz;
274
275 // Sort maxx < maxy < maxz.
276 if (maxx > maxz)
277 std::swap(maxx, maxz);
278 if (maxy > maxz)
279 std::swap(maxy, maxz);
280 else if (maxy < maxx)
281 std::swap(maxx, maxy);
282
283 // Protect against underflow in the computation of eps.
284 if (maxx < 1e-97) /* cbrt(min_double/eps) */ {
285 if (maxx == 0)
286 return 0;
287 }
288 // Protect against overflow in the computation of det.
289 else if (maxz < 1e102) /* cbrt(max_double [hadamard]/4) */ {
290 double eps = 5.1107127829973299e-15 * maxx * maxy * maxz;
291 if (det > eps) return 1;
292 if (det < -eps) return -1;
293 }
294 return STATIC_FILTER_FAILURE;
295 }
296
297 public:
298 typedef typename Interval_nt_advanced::Protector Protector;
299
300 Is_on_positive_side_of_plane_3(const Traits&,const Point_3& p_,const Point_3& q_,const Point_3& r_)
301 : p(p_),q(q_),r(r_)
302 , ak_plane()
303 , ek_plane_ptr(nullptr)
304 {
305 ak_plane.vector =
306 typename Approx_K::Vector_3(Interval_nt_advanced(0., std::numeric_limits<double>::infinity()),
307 0., 0.);
308 double pqx = q.x() - p.x();
309 double pqy = q.y() - p.y();
310 double pqz = q.z() - p.z();
311 double prx = r.x() - p.x();
312 double pry = r.y() - p.y();
313 double prz = r.z() - p.z();
314
315
316 m10 = pqy*prz - pry*pqz;
317 m20 = pqx*prz - prx*pqz;
318 m21 = pqx*pry - prx*pqy;
319
320 double aprx = CGAL::abs(prx);
321 double apry = CGAL::abs(pry);
322 double aprz = CGAL::abs(prz);
323
324 Maxx = CGAL::abs(pqx);
325 if (Maxx < aprx) Maxx = aprx;
326 Maxy = CGAL::abs(pqy);
327 if (Maxy < apry) Maxy = apry;
328 Maxz = CGAL::abs(pqz);
329 if (Maxz < aprz) Maxz = aprz;
330 }
331
332 ~Is_on_positive_side_of_plane_3(){
333 if (ek_plane_ptr!=nullptr) delete ek_plane_ptr;
334 }
335
336 bool operator() (const Point_3& s) const
337 {
338 double psx = s.x() - p.x();
339 double psy = s.y() - p.y();
340 double psz = s.z() - p.z();
341
342 int static_res = static_filtered(psx,psy,psz);
343 if (static_res != STATIC_FILTER_FAILURE)
344 return static_res == 1;
345
346 try{
347 // infinity() is the sentinel for uninitialized `ak_plane`
348 if (ak_plane.vector.x().sup() == std::numeric_limits<double>::infinity())
349 {
350 const typename Approx_K::Point_3 ap = to_AK(p);
351 ak_plane.vector = cross_product(to_AK(q)-ap, to_AK(r)-ap);
352 ak_plane.point = ap;
353 }
354 Uncertain<Sign> res =
355 sign(scalar_product(to_AK(s) - ak_plane.point,
356 ak_plane.vector));
357 if(is_certain(res)) {
358 return (get_certain(res) == POSITIVE);
359 }
360 }
361 catch (Uncertain_conversion_exception&){}
362 if (ek_plane_ptr==nullptr) {
363 const typename Exact_K::Point_3 ep = to_EK(p);
364 ek_plane_ptr = new Vector_plus_point<Exact_K>;
365 ek_plane_ptr->vector = cross_product(to_EK(q)-ep, to_EK(r)-ep);
366 ek_plane_ptr->point = ep;
367 }
368 return sign(scalar_product(to_EK(s) - ek_plane_ptr->point,
369 ek_plane_ptr->vector)) == POSITIVE;
370 }
371 };
372
373
374 BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Traits_has_typedef_Traits_xy_3,Traits_xy_3,false)
375 BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Traits_has_typedef_Traits_yz_3,Traits_xy_3,false)
376 BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(Traits_has_typedef_Traits_xz_3,Traits_xy_3,false)
377
378 template <class T,bool has_projection_traits=
379 Traits_has_typedef_Traits_xy_3<T>::value &&
380 Traits_has_typedef_Traits_yz_3<T>::value &&
381 Traits_has_typedef_Traits_xz_3<T>::value
382 >
383 struct Projection_traits{
384 Projection_traits(const T&){}
385 typedef typename Kernel_traits<typename T::Point_3>::Kernel K;
386 typedef CGAL::Projection_traits_xy_3<K> Traits_xy_3;
387 typedef CGAL::Projection_traits_yz_3<K> Traits_yz_3;
388 typedef CGAL::Projection_traits_xz_3<K> Traits_xz_3;
389
390 Traits_xy_3 construct_traits_xy_3_object()const
391 {return Traits_xy_3();}
392 Traits_yz_3 construct_traits_yz_3_object()const
393 {return Traits_yz_3();}
394 Traits_xz_3 construct_traits_xz_3_object()const
395 {return Traits_xz_3();}
396 };
397
398 template <class T>
399 struct Projection_traits<T,true>{
400 const T& traits;
401 Projection_traits(const T& t):traits(t){}
402 typedef typename T::Traits_xy_3 Traits_xy_3;
403 typedef typename T::Traits_yz_3 Traits_yz_3;
404 typedef typename T::Traits_xz_3 Traits_xz_3;
405 Traits_xy_3 construct_traits_xy_3_object()const
406 {return traits.construct_traits_xy_3_object();}
407 Traits_yz_3 construct_traits_yz_3_object()const
408 {return traits.construct_traits_yz_3_object();}
409 Traits_xz_3 construct_traits_xz_3_object()const
410 {return traits.construct_traits_xz_3_object();}
411 };
412
413 template <class Point_3, class PolygonMesh>
414 void copy_ch2_to_face_graph(const std::list<Point_3>& CH_2, PolygonMesh& P)
415 {
416 typename boost::property_map<PolygonMesh, CGAL::vertex_point_t>::type vpm
417 = get(CGAL::vertex_point, P);
418 typedef boost::graph_traits<PolygonMesh> Graph_traits;
419 typedef typename Graph_traits::vertex_descriptor vertex_descriptor;
420 typedef typename Graph_traits::halfedge_descriptor halfedge_descriptor;
421 typedef typename Graph_traits::face_descriptor face_descriptor;
422 std::vector<vertex_descriptor> vertices;
423 vertices.reserve(CH_2.size());
424 for(const Point_3& p : CH_2){
425 vertices.push_back(add_vertex(P));
426 put(vpm, vertices.back(),p);
427 }
428 face_descriptor f = Euler::add_face(vertices, P);
429
430 // Then triangulate that face
431 const halfedge_descriptor he = halfedge(f, P);
432 halfedge_descriptor other_he = next(next(he, P), P);
433 for(std::size_t i = 3, end = vertices.size(); i < end; ++i) {
434 const halfedge_descriptor next_he = next(other_he, P);
435 Euler::split_face(other_he, he, P);
436 other_he = next_he;
437 }
438 }
439
440
441 template <class InputIterator, class Point_3, class PolygonMesh, class Traits>
442 void coplanar_3_hull(InputIterator first, InputIterator beyond,
443 const Point_3& p1, const Point_3& p2, const Point_3& p3,
444 PolygonMesh& P, const Traits& traits )
445 {
446 typedef typename Convex_hull_3::internal::Projection_traits<Traits> PTraits;
447 typedef typename PTraits::Traits_xy_3 Traits_xy_3;
448 typedef typename PTraits::Traits_yz_3 Traits_yz_3;
449 typedef typename PTraits::Traits_xz_3 Traits_xz_3;
450
451 PTraits ptraits(traits);
452
453 std::list<Point_3> CH_2;
454
455 Traits_xy_3 traits_xy = ptraits.construct_traits_xy_3_object();
456 typename Traits_xy_3::Left_turn_2 left_turn_in_xy = traits_xy.left_turn_2_object();
457 if ( left_turn_in_xy(p1,p2,p3) || left_turn_in_xy(p2,p1,p3) )
458 convex_hull_points_2( first, beyond,
459 std::back_inserter(CH_2),
460 traits_xy );
461 else{
462 Traits_yz_3 traits_yz = ptraits.construct_traits_yz_3_object();
463 typename Traits_yz_3::Left_turn_2 left_turn_in_yz = traits_yz.left_turn_2_object();
464 if ( left_turn_in_yz(p1,p2,p3) || left_turn_in_yz(p2,p1,p3) )
465 convex_hull_points_2( first, beyond,
466 std::back_inserter(CH_2),
467 traits_yz );
468 else{
469 Traits_xz_3 traits_xz = ptraits.construct_traits_xz_3_object();
470 CGAL_assertion_code( typename Traits_xz_3::Left_turn_2 left_turn_in_xz = traits_xz.left_turn_2_object(); )
471 CGAL_assertion( left_turn_in_xz(p1,p2,p3) || left_turn_in_xz(p2,p1,p3) );
472 convex_hull_points_2( first, beyond,
473 std::back_inserter(CH_2),
474 traits_xz );
475 }
476 }
477 copy_ch2_to_face_graph(CH_2, P);
478 }
479
480
481 //
482 // visible is the set of facets visible from point and reachable from
483 // start_facet.
484 //
485 template <class TDS_2, class Traits>
486 void
487 find_visible_set(TDS_2& tds,
488 const typename Traits::Point_3& point,
489 typename TDS_2::Face_handle start,
490 std::list<typename TDS_2::Face_handle>& visible,
491 std::map<typename TDS_2::Vertex_handle, typename TDS_2::Edge>& outside,
492 const Traits& traits)
493 {
494 typedef typename TDS_2::Face_handle Face_handle;
495 typedef typename TDS_2::Vertex_handle Vertex_handle;
496
497 std::vector<Vertex_handle> vertices;
498 vertices.reserve(10);
499 int VISITED=1, BORDER=2;
500 visible.clear();
501 typename std::list<Face_handle>::iterator vis_it;
502 visible.push_back(start);
503 start->info() = VISITED;
504 vertices.push_back(start->vertex(0));
505 vertices.push_back(start->vertex(1));
506 vertices.push_back(start->vertex(2));
507 start->vertex(0)->info() = start->vertex(1)->info() = start->vertex(2)->info() = VISITED;
508
509 for (vis_it = visible.begin(); vis_it != visible.end(); vis_it++)
510 {
511 // check all the neighbors of the current face to see if they have
512 // already been visited or not and if not whether they are visible
513 // or not.
514
515 for(int i=0; i < 3; i++) {
516 // the facet on the other side of the current halfedge
517 Face_handle f = (*vis_it)->neighbor(i);
518 // if haven't already seen this facet
519 if (f->info() == 0) {
520 f->info() = VISITED;
521 Is_on_positive_side_of_plane_3<Traits> is_on_positive_side(
522 traits,f->vertex(0)->point(),f->vertex(2)->point(),f->vertex(1)->point());
523 int ind = f->index(*vis_it);
524 if ( !is_on_positive_side(point) ){ // is visible
525 visible.push_back(f);
526 Vertex_handle vh = f->vertex(ind);
527 if(vh->info() == 0){ vertices.push_back(vh); vh->info() = VISITED;}
528 } else {
529 f->info() = BORDER;
530 f->vertex(TDS_2::cw(ind))->info() = BORDER;
531 f->vertex(TDS_2::ccw(ind))->info() = BORDER;
532 outside.insert(std::make_pair(f->vertex(TDS_2::cw(ind)),
533 typename TDS_2::Edge(f,ind)));
534 }
535 } else if(f->info() == BORDER) {
536 int ind = f->index(*vis_it);
537 f->vertex(TDS_2::cw(ind))->info() = BORDER;
538 f->vertex(TDS_2::ccw(ind))->info() = BORDER;
539 outside.insert(std::make_pair(f->vertex(TDS_2::cw(ind)),
540 typename TDS_2::Edge(f,ind)));
541 }
542 }
543 }
544
545 for(typename std::vector<Vertex_handle>::iterator vit = vertices.begin();
546 vit != vertices.end();
547 ++vit){
548 if((*vit)->info() != BORDER){
549 tds.delete_vertex(*vit);
550 } else {
551 (*vit)->info() = 0;
552 }
553 }
554
555 }
556
557 // using a third template parameter for the point instead of getting it from
558 // the traits class as it should be is required by M$VC6
559 template <class Face_handle, class Traits, class Point>
560 typename std::list<Point>::iterator
561 farthest_outside_point(Face_handle f, std::list<Point>& outside_set,
562 const Traits& traits)
563 {
564
565 typedef typename std::list<Point>::iterator Outside_set_iterator;
566 CGAL_ch_assertion(!outside_set.empty());
567
568 typename Traits::Plane_3 plane =
569 traits.construct_plane_3_object()(f->vertex(0)->point(),
570 f->vertex(1)->point(),
571 f->vertex(2)->point());
572
573 typename Traits::Less_signed_distance_to_plane_3 less_dist_to_plane =
574 traits.less_signed_distance_to_plane_3_object();
575 Outside_set_iterator farthest_it =
576 std::max_element(outside_set.begin(),
577 outside_set.end(),
578 [&less_dist_to_plane,&plane](const Point& p1, const Point& p2)
579 { return less_dist_to_plane(plane, p1, p2); });
580 return farthest_it;
581 }
582
583 template <class Face_handle, class Traits, class Point>
584 void
585 partition_outside_sets(const std::list<Face_handle>& new_facets,
586 std::list<Point>& vis_outside_set,
587 std::list<Face_handle>& pending_facets,
588 const Traits& traits)
589 {
590 typename std::list<Face_handle>::const_iterator f_list_it;
591 typename std::list<Point>::iterator point_it, to_splice;
592
593 // walk through all the new facets and check each unassigned outside point
594 // to see if it belongs to the outside set of this new facet.
595 for (f_list_it = new_facets.begin(); (f_list_it != new_facets.end()) && (! vis_outside_set.empty());
596 ++f_list_it)
597 {
598 Face_handle f = *f_list_it;
599 Is_on_positive_side_of_plane_3<Traits> is_on_positive_side(
600 traits,f->vertex(0)->point(),f->vertex(1)->point(),f->vertex(2)->point());
601 std::list<Point>& point_list = f->points;
602
603 for (point_it = vis_outside_set.begin();point_it != vis_outside_set.end();){
604 if( is_on_positive_side(*point_it) ) {
605 to_splice = point_it;
606 ++point_it;
607 point_list.splice(point_list.end(), vis_outside_set, to_splice);
608 } else {
609 ++point_it;
610 }
611 }
612 if(! point_list.empty()){
613 pending_facets.push_back(f);
614 f->it = boost::prior(pending_facets.end());
615 } else {
616 f->it = pending_facets.end();
617 }
618 }
619
620
621 for (; f_list_it != new_facets.end();++f_list_it)
622 (*f_list_it)->it = pending_facets.end();
623 }
624
625
626
627 template <class TDS_2, class Traits>
628 void
629 ch_quickhull_3_scan(TDS_2& tds,
630 std::list<typename TDS_2::Face_handle>& pending_facets,
631 const Traits& traits)
632 {
633 typedef typename TDS_2::Edge Edge;
634 typedef typename TDS_2::Face_handle Face_handle;
635 typedef typename TDS_2::Vertex_handle Vertex_handle;
636 typedef typename Traits::Point_3 Point_3;
637 typedef std::list<Point_3> Outside_set;
638 typedef typename std::list<Point_3>::iterator Outside_set_iterator;
639 typedef std::map<typename TDS_2::Vertex_handle, typename TDS_2::Edge> Border_edges;
640
641 std::list<Face_handle> visible_set;
642 typename std::list<Face_handle>::iterator vis_set_it;
643 Outside_set vis_outside_set;
644 Border_edges border;
645
646 while (!pending_facets.empty())
647 {
648 vis_outside_set.clear();
649
650 Face_handle f_handle = pending_facets.front();
651
652 Outside_set_iterator farthest_pt_it = farthest_outside_point(f_handle, f_handle->points, traits);
653 Point_3 farthest_pt = *farthest_pt_it;
654 f_handle->points.erase(farthest_pt_it);
655 find_visible_set(tds, farthest_pt, f_handle, visible_set, border, traits);
656
657 // for each visible facet
658 for (vis_set_it = visible_set.begin(); vis_set_it != visible_set.end();
659 vis_set_it++)
660 {
661
662 // add its outside set to the global outside set list
663 std::list<Point_3>& point_list = (*vis_set_it)->points;
664 if(! point_list.empty()){
665 vis_outside_set.splice(vis_outside_set.end(), point_list, point_list.begin(), point_list.end());
666 }
667
668 if((*vis_set_it)->it != pending_facets.end()){
669 pending_facets.erase((*vis_set_it)->it);
670 }
671 (*vis_set_it)->info() = 0;
672 }
673
674 std::vector<Edge> edges;
675 edges.reserve(border.size());
676 typename Border_edges::iterator it = border.begin();
677 Edge e = it->second;
678 e.first->info() = 0;
679 edges.push_back(e);
680 border.erase(it);
681 while(! border.empty()){
682 it = border.find(e.first->vertex(TDS_2::ccw(e.second)));
683 CGAL_ch_assertion(it != border.end());
684 e = it->second;
685 e.first->info() = 0;
686 edges.push_back(e);
687 border.erase(it);
688 }
689
690 // If we want to reuse the faces we must only pass |edges| many, and call delete_face for the others.
691 // Also create facets if necessary
692 std::ptrdiff_t diff = visible_set.size() - edges.size();
693 if(diff < 0){
694 for(int i = 0; i<-diff;i++){
695 visible_set.push_back(tds.create_face());
696 }
697 } else {
698 for(int i = 0; i<diff;i++){
699 tds.delete_face(visible_set.back());
700 visible_set.pop_back();
701 }
702 }
703 Vertex_handle vh = tds.star_hole(edges.begin(), edges.end(), visible_set.begin(), visible_set.end());
704 vh->point() = farthest_pt;
705 vh->info() = 0;
706
707 // now partition the set of outside set points among the new facets.
708
709 partition_outside_sets(visible_set, vis_outside_set,
710 pending_facets, traits);
711
712 }
713 }
714
715 template <class TDS_2, class Traits>
716 void non_coplanar_quickhull_3(std::list<typename Traits::Point_3>& points,
717 TDS_2& tds, const Traits& traits)
718 {
719 typedef typename Traits::Point_3 Point_3;
720
721 typedef typename TDS_2::Face_handle Face_handle;
722 typedef typename TDS_2::Face_iterator Face_iterator;
723 typedef typename std::list<Point_3>::iterator P3_iterator;
724
725 std::list<Face_handle> pending_facets;
726
727 typename Is_on_positive_side_of_plane_3<Traits>::Protector p;
728
729 // for each facet, look at each unassigned point and decide if it belongs
730 // to the outside set of this facet.
731 for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){
732 Is_on_positive_side_of_plane_3<Traits> is_on_positive_side(
733 traits,fit->vertex(0)->point(),fit->vertex(1)->point(),fit->vertex(2)->point() );
734 for (P3_iterator point_it = points.begin() ; point_it != points.end(); )
735 {
736 if( is_on_positive_side(*point_it) ) {
737 P3_iterator to_splice = point_it;
738 ++point_it;
739 fit->points.splice(fit->points.end(), points, to_splice);
740 } else {
741 ++point_it;
742 }
743 }
744 }
745 // add all the facets with non-empty outside sets to the set of facets for
746 // further consideration
747 for(Face_iterator fit = tds.faces_begin(); fit != tds.faces_end(); ++fit){
748 if (! fit->points.empty()){
749 pending_facets.push_back(fit);
750 fit->it = boost::prior(pending_facets.end());
751 } else {
752 fit->it = pending_facets.end();
753 }
754 }
755
756
757 ch_quickhull_3_scan(tds, pending_facets, traits);
758
759 //std::cout << "|V(tds)| = " << tds.number_of_vertices() << std::endl;
760 // CGAL_ch_expensive_postcondition(all_points_inside(points.begin(),
761 // points.end(),P,traits));
762 // CGAL_ch_postcondition(is_strongly_convex_3(P, traits));
763 }
764
765 template <class InputIterator, class PolygonMesh, class Traits>
766 void
767 ch_quickhull_face_graph(std::list<typename Traits::Point_3>& points,
768 InputIterator point1_it, InputIterator point2_it, InputIterator point3_it,
769 PolygonMesh& P,
770 const Traits& traits)
771 {
772 typedef typename Traits::Point_3 Point_3;
773 typedef typename Traits::Plane_3 Plane_3;
774 typedef typename std::list<Point_3>::iterator P3_iterator;
775
776 typedef Triangulation_data_structure_2<
777 Triangulation_vertex_base_with_info_2<int, GT3_for_CH3<Traits> >,
778 Convex_hull_face_base_2<int, Traits> > Tds;
779
780 typedef typename Tds::Vertex_handle Vertex_handle;
781 typedef typename Tds::Face_handle Face_handle;
782
783 // found three points that are not collinear, so construct the plane defined
784 // by these points and then find a point that has maximum distance from this
785 // plane.
786 typename Traits::Construct_plane_3 construct_plane =
787 traits.construct_plane_3_object();
788 Plane_3 plane = construct_plane(*point3_it, *point2_it, *point1_it);
789 typedef typename Traits::Less_signed_distance_to_plane_3 Dist_compare;
790 Dist_compare compare_dist = traits.less_signed_distance_to_plane_3_object();
791
792 typename Traits::Coplanar_3 coplanar = traits.coplanar_3_object();
793 // find both min and max here since using signed distance. If all points
794 // are on the negative side of the plane, the max element will be on the
795 // plane.
796 std::pair<P3_iterator, P3_iterator> min_max;
797 min_max = CGAL::min_max_element(points.begin(), points.end(),
798 [&compare_dist, &plane](const Point_3& p1, const Point_3& p2)
799 { return compare_dist(plane, p1, p2); },
800 [&compare_dist, &plane](const Point_3& p1, const Point_3& p2)
801 { return compare_dist(plane, p1, p2); });
802 P3_iterator max_it;
803 if (coplanar(*point1_it, *point2_it, *point3_it, *min_max.second))
804 {
805 max_it = min_max.first;
806 // want the orientation of the points defining the plane to be positive
807 // so have to reorder these points if all points were on negative side
808 // of plane
809 std::swap(*point1_it, *point3_it);
810 }
811 else
812 max_it = min_max.second;
813
814 // if the maximum distance point is on the plane then all are coplanar
815 if (coplanar(*point1_it, *point2_it, *point3_it, *max_it)) {
816 coplanar_3_hull(points.begin(), points.end(), *point1_it, *point2_it, *point3_it, P, traits);
817 } else {
818 Tds tds;
819 Vertex_handle v0 = tds.create_vertex(); v0->set_point(*point1_it);
820 Vertex_handle v1 = tds.create_vertex(); v1->set_point(*point2_it);
821 Vertex_handle v2 = tds.create_vertex(); v2->set_point(*point3_it);
822 Vertex_handle v3 = tds.create_vertex(); v3->set_point(*max_it);
823
824 v0->info() = v1->info() = v2->info() = v3->info() = 0;
825 Face_handle f0 = tds.create_face(v0,v1,v2);
826 Face_handle f1 = tds.create_face(v3,v1,v0);
827 Face_handle f2 = tds.create_face(v3,v2,v1);
828 Face_handle f3 = tds.create_face(v3,v0,v2);
829 tds.set_dimension(2);
830 v0->set_face(f0);
831 v1->set_face(f0);
832 v2->set_face(f0);
833 v3->set_face(f1);
834 f0->set_neighbors(f2, f3, f1);
835 f1->set_neighbors(f0, f3, f2);
836 f2->set_neighbors(f0, f1, f3);
837 f3->set_neighbors(f0, f2, f1);
838
839 points.erase(point1_it);
840 points.erase(point2_it);
841 points.erase(point3_it);
842 points.erase(max_it);
843 if (!points.empty()){
844 non_coplanar_quickhull_3(points, tds, traits);
845 copy_face_graph(tds,P);
846 }
847 else{
848 CGAL_assertion( traits.has_on_positive_side_3_object()(
849 construct_plane(v2->point(),v1->point(),v0->point()),
850 v3->point()) );
851 make_tetrahedron(v0->point(),v1->point(),v3->point(),v2->point(),P);
852 }
853 }
854
855 }
856
857 } // namespace internal
858 } // namespace Convex_hull_3
859
860 template <class InputIterator, class Traits>
861 void
862 convex_hull_3(InputIterator first, InputIterator beyond,
863 Object& ch_object,
864 const Traits& traits)
865 {
866 typedef typename Traits::Point_3 Point_3;
867 typedef std::list<Point_3> Point_3_list;
868 typedef typename Point_3_list::iterator P3_iterator;
869 typedef std::pair<P3_iterator,P3_iterator> P3_iterator_pair;
870
871 if (first == beyond) // No point
872 return;
873
874 // If the first and last point are equal the collinearity test some lines below will always be true.
875 Point_3_list points(first, beyond);
876 std::size_t size = points.size();
877 while((size > 1) && (points.front() == points.back())){
878 points.pop_back();
879 --size;
880 }
881
882 typename Traits::Collinear_3 collinear = traits.collinear_3_object();
883
884 if ( size == 1 ) // 1 point
885 {
886 ch_object = make_object(*points.begin());
887 return;
888 }
889 else if ( size == 2 ) // 2 points
890 {
891 typedef typename Traits::Segment_3 Segment_3;
892 typename Traits::Construct_segment_3 construct_segment =
893 traits.construct_segment_3_object();
894 Segment_3 seg = construct_segment(*points.begin(), *(++points.begin()));
895 ch_object = make_object(seg);
896 return;
897 }
898 else if ( ( size == 3 ) && (! collinear(*(points.begin()),
899 *(++points.begin()),
900 *(--points.end()) ) ) ) // 3 points
901 {
902 typedef typename Traits::Triangle_3 Triangle_3;
903 typename Traits::Construct_triangle_3 construct_triangle =
904 traits.construct_triangle_3_object();
905 Triangle_3 tri = construct_triangle(*(points.begin()),
906 *(++points.begin()),
907 *(--points.end()));
908 ch_object = make_object(tri);
909 return;
910 }
911
912 // at least 4 points
913
914 P3_iterator point1_it = points.begin();
915 P3_iterator point2_it = points.begin();
916 point2_it++;
917 P3_iterator point3_it = points.end();
918 point3_it--;
919
920 // find three that are not collinear
921 while (point2_it != points.end() &&
922 collinear(*point1_it,*point2_it,*point3_it))
923 point2_it++;
924
925
926 // all are collinear, so the answer is a segment
927 if (point2_it == points.end())
928 {
929 typedef typename Traits::Less_distance_to_point_3 Less_dist;
930
931 Less_dist less_dist = traits.less_distance_to_point_3_object();
932 P3_iterator_pair endpoints =
933 min_max_element(points.begin(), points.end(),
934 [&points, &less_dist](const Point_3& p1, const Point_3& p2)
935 { return less_dist(*points.begin(), p1, p2); },
936 [&points, &less_dist](const Point_3& p1, const Point_3& p2)
937 { return less_dist(*points.begin(), p1, p2); });
938
939 typename Traits::Construct_segment_3 construct_segment =
940 traits.construct_segment_3_object();
941 typedef typename Traits::Segment_3 Segment_3;
942
943 Segment_3 seg = construct_segment(*endpoints.first, *endpoints.second);
944 ch_object = make_object(seg);
945 return;
946 }
947
948 // result will be a polyhedron
949 typedef typename Convex_hull_3::internal::Default_polyhedron_for_Chull_3<Traits>::type Polyhedron;
950 Polyhedron P;
951
952 P3_iterator minx, maxx, miny, it;
953 minx = maxx = miny = it = points.begin();
954 ++it;
955 for(; it != points.end(); ++it){
956 if(it->x() < minx->x()) minx = it;
957 if(it->x() > maxx->x()) maxx = it;
958 if(it->y() < miny->y()) miny = it;
959 }
960 if(! collinear(*minx, *maxx, *miny) ){
961 Convex_hull_3::internal::ch_quickhull_face_graph(points, minx, maxx, miny, P, traits);
962 } else {
963 Convex_hull_3::internal::ch_quickhull_face_graph(points, point1_it, point2_it, point3_it, P, traits);
964 }
965 CGAL_assertion(num_vertices(P)>=3);
966 typename boost::graph_traits<Polyhedron>::vertex_iterator b,e;
967 boost::tie(b,e) = vertices(P);
968 if (num_vertices(P) == 3){
969 typename boost::property_map<Polyhedron, vertex_point_t>::type vpmap = get(CGAL::vertex_point, P);
970 typedef typename Traits::Triangle_3 Triangle_3;
971 typename Traits::Construct_triangle_3 construct_triangle =
972 traits.construct_triangle_3_object();
973 typedef typename Traits::Point_3 Point_3;
974 Point_3 p = get(vpmap, *b); ++b;
975 Point_3 q = get(vpmap, *b); ++b;
976 Point_3 r = get(vpmap, *b);
977 Triangle_3 tri = construct_triangle(p,q,r);
978 ch_object = make_object(tri);
979 }
980 else
981 ch_object = make_object(P);
982 }
983
984
985 template <class InputIterator>
986 void convex_hull_3(InputIterator first, InputIterator beyond,
987 Object& ch_object)
988 {
989 typedef typename std::iterator_traits<InputIterator>::value_type Point_3;
990 typedef typename Convex_hull_3::internal::Default_traits_for_Chull_3<Point_3>::type Traits;
991 convex_hull_3(first, beyond, ch_object, Traits());
992 }
993
994 template <class InputIterator, class PolygonMesh, class Traits>
995 void convex_hull_3(InputIterator first, InputIterator beyond,
996 PolygonMesh& polyhedron,
997 const Traits& traits)
998 {
999 typedef typename Traits::Point_3 Point_3;
1000 typedef std::list<Point_3> Point_3_list;
1001 typedef typename Point_3_list::iterator P3_iterator;
1002 typedef std::pair<P3_iterator,P3_iterator> P3_iterator_pair;
1003
1004 if(first == beyond)
1005 return;
1006
1007 Point_3_list points(first, beyond);
1008 clear(polyhedron);
1009
1010 typename Traits::Collinear_3 collinear = traits.collinear_3_object();
1011 typename Traits::Equal_3 equal = traits.equal_3_object();
1012 P3_iterator point1_it = points.begin();
1013 P3_iterator point2_it = points.begin();
1014 point2_it++;
1015
1016 while (point2_it != points.end() && equal(*point1_it,*point2_it))
1017 ++point2_it;
1018
1019 // if there is only one point or all points are equal
1020 if(point2_it == points.end())
1021 {
1022 Convex_hull_3::internal::add_isolated_points(*point1_it, polyhedron);
1023 return;
1024 }
1025
1026 P3_iterator point3_it = point2_it;
1027 ++point3_it;
1028
1029 while (point3_it != points.end() && collinear(*point1_it,*point2_it,*point3_it))
1030 ++point3_it;
1031
1032 if (point3_it == points.end())
1033 {
1034 // only 2 points or all points are collinear
1035 typedef typename Traits::Less_distance_to_point_3 Less_dist;
1036 Less_dist less_dist = traits.less_distance_to_point_3_object();
1037 P3_iterator_pair endpoints =
1038 min_max_element(points.begin(), points.end(),
1039 [&points, &less_dist](const Point_3& p1, const Point_3& p2)
1040 { return less_dist(*points.begin(), p1, p2); },
1041 [&points, &less_dist](const Point_3& p1, const Point_3& p2)
1042 { return less_dist(*points.begin(), p1, p2); });
1043 Convex_hull_3::internal::add_isolated_points(*endpoints.first, polyhedron);
1044 Convex_hull_3::internal::add_isolated_points(*endpoints.second, polyhedron);
1045 return;
1046 }
1047
1048 Convex_hull_3::internal::ch_quickhull_face_graph(points, point1_it, point2_it, point3_it,
1049 polyhedron, traits);
1050 }
1051
1052 template <class InputIterator, class PolygonMesh>
1053 void convex_hull_3(InputIterator first, InputIterator beyond,
1054 PolygonMesh& polyhedron,
1055 // workaround to avoid ambiguity with next overload.
1056 typename std::enable_if<CGAL::is_iterator<InputIterator>::value>::type* = 0)
1057 {
1058 typedef typename std::iterator_traits<InputIterator>::value_type Point_3;
1059 typedef typename Convex_hull_3::internal::Default_traits_for_Chull_3<Point_3, PolygonMesh>::type Traits;
1060 convex_hull_3(first, beyond, polyhedron, Traits());
1061 }
1062
1063 template <class VertexListGraph, class PolygonMesh, class NamedParameters>
1064 void convex_hull_3(const VertexListGraph& g,
1065 PolygonMesh& pm,
1066 const NamedParameters& np)
1067 {
1068 using CGAL::parameters::choose_parameter;
1069 using CGAL::parameters::get_parameter;
1070
1071 typedef typename GetVertexPointMap<VertexListGraph, NamedParameters>::const_type Vpmap;
1072 typedef CGAL::Property_map_to_unary_function<Vpmap> Vpmap_fct;
1073 Vpmap vpm = choose_parameter(get_parameter(np, internal_np::vertex_point),
1074 get_const_property_map(boost::vertex_point, g));
1075
1076 Vpmap_fct v2p(vpm);
1077 convex_hull_3(boost::make_transform_iterator(vertices(g).begin(), v2p),
1078 boost::make_transform_iterator(vertices(g).end(), v2p), pm);
1079 }
1080
1081 template <class VertexListGraph, class PolygonMesh>
1082 void convex_hull_3(const VertexListGraph& g,
1083 PolygonMesh& pm)
1084 {
1085 convex_hull_3(g,pm,CGAL::parameters::all_default());
1086 }
1087
1088 template <class InputRange, class OutputIterator, class Traits>
1089 OutputIterator
1090 extreme_points_3(const InputRange& range,
1091 OutputIterator out,
1092 const Traits& traits)
1093 {
1094 Convex_hull_3::internal::Output_iterator_wrapper<OutputIterator> wrapper(out);
1095 convex_hull_3(range.begin(), range.end(), wrapper, traits);
1096 return out;
1097 }
1098
1099 template <class InputRange, class OutputIterator>
1100 OutputIterator
1101 extreme_points_3(const InputRange& range, OutputIterator out)
1102 {
1103 typedef typename InputRange::const_iterator Iterator_type;
1104 typedef typename std::iterator_traits<Iterator_type>::value_type Point_3;
1105 typedef typename Convex_hull_3::internal::Default_traits_for_Chull_3<Point_3>::type Traits;
1106
1107 return extreme_points_3(range, out, Traits());
1108 }
1109
1110 } // namespace CGAL
1111
1112 #include <CGAL/enable_warnings.h>
1113
1114 #endif // CGAL_CONVEX_HULL_3_H
1115