From e2e155fdb409ea31cbcc3ec28bccdc9bb429c07a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Loriot?= Date: Wed, 11 Jul 2012 16:21:54 +0000 Subject: [PATCH] add include files for the demo --- .gitattributes | 9 + .../include/CGAL/corefinement_operations.h | 485 ++++ .../Combinatorial_map_for_corefinement.h | 110 + .../corefinement/Polyhedron_constness_types.h | 57 + .../Polyhedron_subset_extraction.h | 200 ++ .../intersection_coplanar_triangles_3.h | 327 +++ .../intersection_triangle_segment_3.h | 166 ++ ...intersection_triangle_segment_3_coplanar.h | 395 +++ .../CGAL/intersection_of_Polyhedra_3.h | 1958 +++++++++++++++ ...ection_of_Polyhedra_3_refinement_visitor.h | 2149 +++++++++++++++++ 10 files changed, 5856 insertions(+) create mode 100644 Polyhedron/include/CGAL/corefinement_operations.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/Polyhedron_constness_types.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/Polyhedron_subset_extraction.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/intersection_coplanar_triangles_3.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3.h create mode 100644 Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3_coplanar.h create mode 100644 Polyhedron/include/CGAL/intersection_of_Polyhedra_3.h create mode 100644 Polyhedron/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h diff --git a/.gitattributes b/.gitattributes index 8ef556f9783..f44c136b2df 100644 --- a/.gitattributes +++ b/.gitattributes @@ -3279,6 +3279,15 @@ Polyhedron/examples/Polyhedron/corner.off -text svneol=unset#application/octet-s Polyhedron/examples/Polyhedron/corner_with_hole.off -text svneol=unset#application/octet-stream Polyhedron/examples/Polyhedron/corner_with_sharp_edge.off -text svneol=unset#application/octet-stream Polyhedron/examples/Polyhedron/cross.off -text svneol=unset#application/octet-stream +Polyhedron/include/CGAL/corefinement_operations.h -text +Polyhedron/include/CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h -text +Polyhedron/include/CGAL/internal/corefinement/Polyhedron_constness_types.h -text +Polyhedron/include/CGAL/internal/corefinement/Polyhedron_subset_extraction.h -text +Polyhedron/include/CGAL/internal/corefinement/intersection_coplanar_triangles_3.h -text +Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3.h -text +Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3_coplanar.h -text +Polyhedron/include/CGAL/intersection_of_Polyhedra_3.h -text +Polyhedron/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h -text Polynomial/doc_tex/Polynomial/Polynomial.png -text Polynomial/doc_tex/Polynomial_ref/PolynomialTraits_d_ConstructCoefficientConstIteratorRange.tex -text Polynomial/doc_tex/Polynomial_ref/PolynomialTraits_d_ConstructInnermostCoefficientConstIteratorRange.tex -text diff --git a/Polyhedron/include/CGAL/corefinement_operations.h b/Polyhedron/include/CGAL/corefinement_operations.h new file mode 100644 index 00000000000..96e68d19b90 --- /dev/null +++ b/Polyhedron/include/CGAL/corefinement_operations.h @@ -0,0 +1,485 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_COREFINEMENT_OPERATIONS_H +#define CGAL_COREFINEMENT_OPERATIONS_H + +#include +#include + +namespace CGAL{ +/** \cond */ +namespace internal{ + +template +class Import_volume_as_polyhedron : public CGAL::Modifier_base { + typedef typename HDS::Halfedge_handle Halfedge_handle; + typedef typename HDS::Vertex_handle Vertex_handle; + typedef typename HDS::Face_handle Face_handle; + typedef typename HDS::Vertex Vertex; + typedef typename HDS::Halfedge Halfedge; + typedef typename HDS::Face Face; + + //data members + Face_handle current_face; + std::vector< typename Vertex::Point > points; + unsigned nb_edges; + std::vector > faces; + + typename HDS::Halfedge::Base* + unlock_halfedge(Halfedge_handle h){ + return static_cast(&(*h)); + } + + typename HDS::Face::Base* + unlock_face(Face_handle f){ + return static_cast(&(*f)); + } + +public: + + //to import each piece individually + template + Import_volume_as_polyhedron(const Combinatorial_map_3& map,typename Combinatorial_map_3::Dart_const_handle dart):nb_edges(0) + { + typedef Combinatorial_map_3 CMap; + typedef std::map Vertex_map; + Vertex_map vertex_map; + unsigned int index=0; + //recover all the vertices in the current volumeiterator over all the point + //map the vertex to the index of the point in the vector + for (typename CMap::template One_dart_per_incident_cell_const_range<0,3>::const_iterator + it=map.template one_dart_per_incident_cell<0,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<0,3>(dart).end(); + it!=itend; ++it) + { + points.push_back(it->template attribute<0>()->point()); + vertex_map.insert(std::make_pair(&it->template attribute<0>()->point(),index++)); + } + + //count the number of edges + nb_edges+=map.template one_dart_per_incident_cell<1,3>(dart).size(); + + //recover one dart per face + for (typename CMap::template One_dart_per_incident_cell_const_range<2,3>::const_iterator + it=map.template one_dart_per_incident_cell<2,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<2,3>(dart).end(); + it!=itend; ++it) + { + //warning: the convention used into a polyhedron is that the normal + // of a triangle indicates the outside of the object; thus + // we need to reverse the orientation of the faces of the + // combinatorial map. + unsigned int i=vertex_map[&it->template attribute<0>()->point()]; + unsigned int j=vertex_map[&it->beta(0)->template attribute<0>()->point()]; + unsigned int k=vertex_map[&it->beta(1)->template attribute<0>()->point()]; + faces.push_back(CGAL::cpp0x::make_tuple(i,j,k)); + } + } + + //for intersection and symetric difference + template + Import_volume_as_polyhedron(const Combinatorial_map_3& map, + Iterator dart_begin, + Iterator dart_end):nb_edges(0) + { + typedef Combinatorial_map_3 CMap; + typedef std::map Vertex_map; + Vertex_map vertex_map; + unsigned int index=0; + + for (Iterator it=dart_begin;it!=dart_end;++it) + { + typename Combinatorial_map_3::Dart_const_handle dart=*it; + //recover all the vertices in the current volumeiterator over all the point + //map the vertex to the index of the point in the vector + for (typename CMap::template One_dart_per_incident_cell_const_range<0,3>::const_iterator + it=map.template one_dart_per_incident_cell<0,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<0,3>(dart).end(); + it!=itend; ++it) + { + if ( vertex_map.insert(std::make_pair(&it->template attribute<0>()->point(),index)).second ) + { + points.push_back(it->template attribute<0>()->point()); + ++index; + } + } + + //count the number of edges + nb_edges+=map.template one_dart_per_incident_cell<1,3>(dart).size(); + + //recover one dart per face + for (typename CMap::template One_dart_per_incident_cell_const_range<2,3>::const_iterator + it=map.template one_dart_per_incident_cell<2,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<2,3>(dart).end(); + it!=itend; ++it) + { + //warning: the convention used into a polyhedron is that the normal + // of a triangle indicates the outside of the object; thus + // we need to reverse the orientation of the faces of the + // combinatorial map. + unsigned int i=vertex_map[&it->template attribute<0>()->point()]; + unsigned int j=vertex_map[&it->beta(0)->template attribute<0>()->point()]; + unsigned int k=vertex_map[&it->beta(1)->template attribute<0>()->point()]; + faces.push_back(CGAL::cpp0x::make_tuple(i,j,k)); + } + } + } + + //for union : use the inverse of the complementary + template + Import_volume_as_polyhedron(const Combinatorial_map_3& map, + Iterator dart_begin, + Iterator dart_end,bool):nb_edges(0) + { + typedef Combinatorial_map_3 CMap; + typedef std::map Vertex_map; + Vertex_map vertex_map; + unsigned int index=0; + + for (Iterator it=dart_begin;it!=dart_end;++it) + { + typename Combinatorial_map_3::Dart_const_handle dart=*it; + //recover all the vertices in the current volumeiterator over all the point + //map the vertex to the index of the point in the vector + for (typename CMap::template One_dart_per_incident_cell_const_range<0,3>::const_iterator + it=map.template one_dart_per_incident_cell<0,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<0,3>(dart).end(); + it!=itend; ++it) + { + if (vertex_map.insert(std::make_pair(&it->template attribute<0>()->point(),index)).second ) + { + points.push_back(it->template attribute<0>()->point()); + ++index; + } + } + + //count the number of edges + nb_edges+=map.template one_dart_per_incident_cell<1,3>(dart).size(); + + //recover one dart per face + for (typename CMap::template One_dart_per_incident_cell_const_range<2,3>::const_iterator + it=map.template one_dart_per_incident_cell<2,3>(dart).begin(), + itend=map.template one_dart_per_incident_cell<2,3>(dart).end(); + it!=itend; ++it) + { + //warning: the convention used into a polyhedron is that the normal + // of a triangle indicates the outside of the object; thus + // we need to reverse the orientation of the faces of the + // combinatorial map. Since to get the complementary we + // also need to reverse the orientation, we finally do + // not change it. + unsigned int i=vertex_map[&it->template attribute<0>()->point()]; + unsigned int j=vertex_map[&it->beta(1)->template attribute<0>()->point()]; + unsigned int k=vertex_map[&it->beta(0)->template attribute<0>()->point()]; + faces.push_back(CGAL::cpp0x::make_tuple(i,j,k)); + } + } + } + + + void operator()( HDS& hds) + { + CGAL::Polyhedron_incremental_builder_3 B(hds, true); + B.begin_surface( points.size(), faces.size(),2*nb_edges); + + //insert vertices + for (typename std::vector::iterator it=points.begin();it!=points.end();++it) + B.add_vertex(*it); + + //create faces + for (std::vector >::iterator it=faces.begin();it!=faces.end();++it) + { + B.begin_facet(); + B.add_vertex_to_facet(CGAL::cpp0x::get<0>(*it)); + B.add_vertex_to_facet(CGAL::cpp0x::get<1>(*it)); + B.add_vertex_to_facet(CGAL::cpp0x::get<2>(*it)); + B.end_facet(); + } + B.end_surface(); + } +}; + +} +/** \endcond */ + +/*! \class Polyhedron_corefinement corefinement_operations.h CGAL/corefinement_operations.h + * Function object to compute the decomposition of the space induced by two polyhedra. + * @tparam Polyhedron must be an instantiation of CGAL::Polyhedron_3. + * @tparam Kernel must be a CGAL Kernel compatible with the underlying kernel of Polyhedron. + */ +template +class Polyhedron_corefinement +{ + typedef internal::Import_volume_as_polyhedron Volume_import_modifier; + +public: + /** Enumeration of the different feature tags, listing all kind of space decomposition the functor can compute given two input polyhedra P and Q.*/ + enum Boolean_operation_tag + { + Join_tag=1, /*!< the union of P and Q. */ + Intersection_tag=2, /*!< the intersection of P and Q. */ + P_minus_Q_tag=4, /*!< P minus Q. */ + Q_minus_P_tag=8, /*!< Q minus P. */ + Parts_of_P_tag=32, /*!< decomposition of the volume bounded by P induced by Q. */ + Parts_of_Q_tag=64, /*!< decomposition of the volume bounded by Q induced by P. */ + Decomposition_tag=16 /*!< Both decompositions of P and Q. */ + }; + + /** + * This computes different polyhedra according to the value of features. + * Each input polyhedron is expected to bound a volume: the volume is bounded by the surface polyhedron with facet + * normals pointing outside the volume. Each facet is supposed to be counterclockwise oriented, that is its vertex + * sequence (induced by halfedges) is seen counterclockwise from the side of the facet pointed by its normal. If one or + * both polyhedra are not closed, the algorithm will end up correctly if each intersection polyline separates each surface + * in two components. In that case for each triangle of each surface path boundary, the interior of the volume is considered + * to be on the side opposite of that pointed by it normals (but the orientation must be consistent on each patch). + * the surface the volume bounded by an open surface is considered to be an infinite volume above + * or below the surface (the side not pointed by normal vectors). + * @tparam Polyline_output_iterator must be an output iterator of std::vector. + * @tparam Polyhedron_ptr_and_type_output_iterator an output iterator of std::pair. + * @param P is the first input triangulated polyhedron. Note that a reference is taken and P will be updated to contain the 1D intersection between the two surfaces P and Q. + * @param Q is the first input triangulated polyhedron. Note that a reference is taken and Q will be updated to contain the 1D intersection between the two surfaces P and Q. + * @param polyline_output is an output iterator that collects intersection polylines between P and Q. + * @param poly_output is an output iterator that collects output polyhedra. Note that each polyhedron is allocated within this function (thus must be explicitly deleted when no longer used). The integer is a feature tag corresponding to the volume represented by the polyhedron of the pair. + * @param features is an integer indicating what polyhedra the function must compute. If several kind of polyhedra are expected, feature tags must be combined by |. For example if features = Polyhedron_corefinement::Join_tag | Polyhedron_corefinement::Intersection_tag, then poly_output will collect two polyhedra, the union and the intersection of P and Q. + */ + template + void operator()( Polyhedron& P, Polyhedron& Q, + Polyline_output_iterator polyline_output, + Polyhedron_ptr_and_type_output_iterator poly_output, + int features) const + { + typedef CGAL::Node_visitor_refine_polyhedra Split_visitor; + Split_visitor visitor; + CGAL::Intersection_of_Polyhedra_3 polyline_intersections(visitor); + + polyline_intersections(P, Q, polyline_output); + + typedef typename Split_visitor::Combinatorial_map_3 Combinatorial_map_3; + typedef typename Split_visitor::Volume_info Volume_info; + typedef typename Combinatorial_map_3::Dart_const_handle Dart_const_handle; + + const typename Split_visitor::Combinatorial_map_3& final_map=visitor.combinatorial_map(); + + typename Combinatorial_map_3::template One_dart_per_cell_const_range<3> cell_range=final_map.template one_dart_per_cell<3>(); + + std::list intersection; + std::list union_; + std::list P_minus_Q; + std::list Q_minus_P; + + for (typename Combinatorial_map_3::template One_dart_per_cell_const_range<3>::const_iterator + it = cell_range.begin(), it_end= cell_range.end(); + it!= it_end; + ++it ) + { + + const Volume_info& info=it->template attribute<3>()->info(); + std::size_t inside_size=info.inside.size(); + std::size_t outside_size=info.outside.size(); + + if ( inside_size + outside_size != 2){ + std::cerr << "Error: One volume cannot be represented using a polyhedron. Aborted.\n"; + break; + } + + switch (outside_size) + { + case 2: + if (features & Join_tag) union_.push_back(it); + break; + case 0: + if (features & Intersection_tag) intersection.push_back(it); + break; + default: + if ( *info.inside.begin() == &P ) + { if (features & P_minus_Q_tag) P_minus_Q.push_back(it); } + else + { if (features & Q_minus_P_tag) Q_minus_P.push_back(it); } + } + + if ( features&Decomposition_tag ) + { + Volume_import_modifier modifier(final_map,it); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++ = std::make_pair( new_poly,static_cast(Decomposition_tag) ); + } + + if ( features&Parts_of_P_tag && info.inside.find(&P)!=info.inside.end() ) + { + Volume_import_modifier modifier(final_map,it); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++ = std::make_pair( new_poly,static_cast(Parts_of_P_tag) ); + } + + if ( features&Parts_of_Q_tag && info.inside.find(&Q)!=info.inside.end() ) + { + Volume_import_modifier modifier(final_map,it); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++ = std::make_pair( new_poly,static_cast(Parts_of_Q_tag) ); + } + } + + if (!intersection.empty()) + { + Volume_import_modifier modifier(final_map,intersection.begin(),intersection.end()); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++=std::make_pair( new_poly,static_cast(Intersection_tag) ); + } + + if (!P_minus_Q.empty()) + { + Volume_import_modifier modifier(final_map,P_minus_Q.begin(),P_minus_Q.end()); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++=std::make_pair( new_poly,static_cast(P_minus_Q_tag) ); + } + + if (!Q_minus_P.empty()) + { + Volume_import_modifier modifier(final_map,Q_minus_P.begin(),Q_minus_P.end()); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++=std::make_pair( new_poly,static_cast(Q_minus_P_tag) ); + } + + if (!union_.empty()) + { + Volume_import_modifier modifier(final_map,union_.begin(),union_.end(),true); + Polyhedron* new_poly=new Polyhedron(); + new_poly->delegate(modifier); + *poly_output++=std::make_pair( new_poly,static_cast(Join_tag) ); + } + } + /** + * This updates P according to the value of features. + * Each input polyhedron is expected to bound a volume: the volume is bounded by the surface polyhedron with facet + * normals pointing outside the volume. Each facet is supposed to be counterclockwise oriented, that is its vertex + * sequence (induced by halfedges) is seen counterclockwise from the side of the facet pointed by its normal. If one or + * both polyhedra are not closed, the algorithm will end up correctly if each intersection polyline separates each surface + * in two components. In that case for each triangle of each surface path boundary, the interior of the volume is considered + * to be on the side opposite of that pointed by it normals (but the orientation must be consistent on each patch). + * the surface the volume bounded by an open surface is considered to be an infinite volume above + * or below the surface (the side not pointed by normal vectors). + * @tparam Polyline_output_iterator must be an output iterator of std::vector. + * @param P is the first input triangulated polyhedron. Note that a reference is taken and P will be updated to contain the 1D intersection between the two surfaces P and Q. + * @param Q is the first input triangulated polyhedron. Note that a reference is taken and Q will be updated to contain the 1D intersection between the two surfaces P and Q. + * @param polyline_output is an output iterator that collects intersection polylines between P and Q. + * @param features is either Join_tag, Intersection_tag, P_minus_Q_tag or Q_minus_P_tag + */ + template + void operator()( Polyhedron& P, Polyhedron& Q, + Polyline_output_iterator polyline_output, + Boolean_operation_tag features) const + { + typedef CGAL::Node_visitor_refine_polyhedra Split_visitor; + Split_visitor visitor; + CGAL::Intersection_of_Polyhedra_3 polyline_intersections(visitor); + + polyline_intersections(P, Q, polyline_output); + + typedef typename Split_visitor::Combinatorial_map_3 Combinatorial_map_3; + typedef typename Split_visitor::Volume_info Volume_info; + typedef typename Combinatorial_map_3::Dart_const_handle Dart_const_handle; + + const typename Split_visitor::Combinatorial_map_3& final_map=visitor.combinatorial_map(); + + typename Combinatorial_map_3::template One_dart_per_cell_const_range<3> cell_range=final_map.template one_dart_per_cell<3>(); + + std::list darts; + + for (typename Combinatorial_map_3::template One_dart_per_cell_const_range<3>::const_iterator + it = cell_range.begin(), it_end= cell_range.end(); + it!= it_end; + ++it ) + { + + const Volume_info& info=it->template attribute<3>()->info(); + std::size_t inside_size=info.inside.size(); + std::size_t outside_size=info.outside.size(); + + if ( inside_size + outside_size != 2){ + std::cerr << "Error: One volume cannot be represented using a polyhedron. Aborted.\n"; + break; + } + + switch (outside_size) + { + case 2: + if ( features==Join_tag ) darts.push_back(it); + break; + case 0: + if ( features==Intersection_tag ) darts.push_back(it); + break; + default: + if ( *info.inside.begin() == &P ) + { if (features == P_minus_Q_tag) darts.push_back(it); } + else + { if (features == Q_minus_P_tag) darts.push_back(it); } + } + } + + P.clear(); + + if (!darts.empty()) + { + Volume_import_modifier modifier= + features==Join_tag? + Volume_import_modifier(final_map,darts.begin(),darts.end(),true) + : Volume_import_modifier(final_map,darts.begin(),darts.end()); + P.delegate(modifier); + } + } + + /** \cond */ + static std::string get_type_str(const std::string& Pname,const std::string& Qname,int i) + { + switch (i) + { + case Join_tag: + return Pname+std::string("_union_")+Qname; + case P_minus_Q_tag: + return Pname+std::string("_minus_")+Qname; + case Q_minus_P_tag: + return Qname+std::string("_minus_")+Pname; + case Intersection_tag: + return Pname+std::string("_inter_")+Qname; + case Decomposition_tag: + return std::string("Decomposition"); + } + return std::string("Unknow"); + } + /** \endcond */ +}; + + + + +} + +#endif // CGAL_COREFINEMENT_OPERATIONS_H diff --git a/Polyhedron/include/CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h b/Polyhedron/include/CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h new file mode 100644 index 00000000000..d28dff4768d --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/Combinatorial_map_for_corefinement.h @@ -0,0 +1,110 @@ +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// 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. +// +// $UR$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_COMBINATORIAL_MAP_FOR_COREFINEMENT_H +#define CGAL_INTERNAL_COMBINATORIAL_MAP_FOR_COREFINEMENT_H + +#include +#include +#include +#include + +namespace CGAL{ + namespace internal_IOP{ + +template +struct Volume_info{ + std::set outside; + std::set inside; + bool is_empty; + Volume_info():is_empty(false){} +}; + +struct Volume_on_merge +{ + template + void operator() (Attribute& a1,const Attribute& a2) const + { + CGAL_assertion(!a1.info().is_empty && !a2.info().is_empty); + std::copy(a2.info().outside.begin(),a2.info().outside.end(),std::inserter(a1.info().outside,a1.info().outside.begin())); + std::copy(a2.info().inside.begin(),a2.info().inside.end(),std::inserter(a1.info().inside,a1.info().inside.begin())); + } +}; + +#ifndef NDEBUG +struct Point_on_merge +{ + template + void operator() (Attribute& a1,const Attribute& a2) const + { + CGAL_assertion(a1.point()==a2.point() ); + } +}; +#endif + + +template < class Refs, class T, class Point_, + class Functor_on_merge_=CGAL::Null_functor, + class Functor_on_split_=CGAL::Null_functor > +class My_cell_attribute_with_point : + public CGAL::Cell_attribute_without_info +{ + Point_ mpoint; +public: + typedef Point_ Point; + typedef Functor_on_merge_ Functor_on_merge; + typedef Functor_on_split_ Functor_on_split; + + My_cell_attribute_with_point(){} + My_cell_attribute_with_point(const Point& apoint) : mpoint(apoint) {} + Point& point() { return mpoint; } + const Point& point() const { return mpoint; } + +}; + +template +class Item_with_points_and_volume_info +{ +public: + public: + static const unsigned int dimension = 3; + static const unsigned int NB_MARKS = 32; + template + struct Dart_wrapper + { + typedef CGAL::Dart<3, Refs > Dart; + typedef Traits_ Traits; + typedef typename Traits::FT FT; + typedef typename Traits::Point_3 Point; + typedef typename Traits::Vector_3 Vector; + #ifndef NDEBUG + typedef My_cell_attribute_with_point< Refs,CGAL::Tag_true,Point,Point_on_merge> Vertex_attribute; + #else + typedef My_cell_attribute_with_point< Refs,CGAL::Tag_true,Point> Vertex_attribute; + #endif + typedef CGAL::Cell_attribute< Refs,Volume_info,CGAL::Tag_true,Volume_on_merge > Volume_attribute; + typedef CGAL::cpp0x::tuple< Vertex_attribute, + void, + void, + Volume_attribute> Attributes; + }; +}; + +} } //namespace CGAL::internal_IOP + +#endif //CGAL_INTERNAL_COMBINATORIAL_MAP_FOR_COREFINEMENT_H diff --git a/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_constness_types.h b/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_constness_types.h new file mode 100644 index 00000000000..0dfdf0a95f9 --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_constness_types.h @@ -0,0 +1,57 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_POLYHEDRON_CONSTNESS_TYPES_H +#define CGAL_INTERNAL_POLYHEDRON_CONSTNESS_TYPES_H + +namespace CGAL { +namespace internal_IOP{ + +template +struct Polyhedron_types; + +template +struct Polyhedron_types{ + typedef Polyhedron& Polyhedron_ref; + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_iterator Halfedge_iterator; + typedef typename Polyhedron::Facet_iterator Facet_iterator; + typedef typename Polyhedron::Facet_handle Facet_handle; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + typedef typename Polyhedron::Vertex Vertex; + typedef typename Polyhedron::Halfedge Halfedge; + typedef typename Polyhedron::Facet Facet; +}; + +template +struct Polyhedron_types{ + typedef const Polyhedron& Polyhedron_ref; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_const_iterator Halfedge_iterator; + typedef typename Polyhedron::Facet_const_iterator Facet_iterator; + typedef typename Polyhedron::Facet_const_handle Facet_handle; + typedef typename Polyhedron::Vertex_const_handle Vertex_handle; + typedef const typename Polyhedron::Vertex Vertex; + typedef const typename Polyhedron::Halfedge Halfedge; + typedef const typename Polyhedron::Facet Facet; +}; + +} } //namespace CGAL::internal_IOP + +#endif //CGAL_INTERNAL_POLYHEDRON_CONSTNESS_TYPES_H diff --git a/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_subset_extraction.h b/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_subset_extraction.h new file mode 100644 index 00000000000..54867a5853f --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/Polyhedron_subset_extraction.h @@ -0,0 +1,200 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_POLYHEDRON_SUBSET_EXTRACTION_H +#define CGAL_INTERNAL_POLYHEDRON_SUBSET_EXTRACTION_H + +#include +#include +#include +#include + +namespace CGAL { + namespace internal{ + +template +struct Compare_handle_ptr{ + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle; + + bool operator()(Facet_const_handle f1,Facet_const_handle f2) const { + return &(*f1) < &(*f2); + } + + bool operator()(Vertex_const_handle v1,Vertex_const_handle v2) const { + return &(*v1) < &(*v2); + } +}; + +struct Dummy_true{ + template + bool operator()(T) const {return true;} +}; + +template +class Build_polyhedron_subset : public ::CGAL::Modifier_base { + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef typename Polyhedron::Vertex_const_handle Vertex_const_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + + typedef typename HDS::Vertex::Point Point; + std::list points; + std::list< std::vector > facets; + + template + typename Polyhedron::Halfedge_const_handle get_facet_halfedge(Facet_iterator facet_it) const + { + return (*facet_it)->halfedge(); + } + + typename Polyhedron::Halfedge_const_handle get_facet_halfedge(typename Polyhedron::Facet_const_handle facet) const + { + return facet->halfedge(); + } + +public: + template + Build_polyhedron_subset(const Polyhedron&,Facets_const_iterator begin,Facets_const_iterator end) + { + typedef std::map > Vertices; + Vertices vertices; + unsigned int index=0; + //get vertices and get face description relatively to the restricted set of vertices + for (Facets_const_iterator it=begin;it!=end;++it) + { + Halfedge_const_handle start=get_facet_halfedge(it); + Halfedge_const_handle curr=start; + facets.push_back(std::vector()); + std::vector& indices = facets.back(); + do{ + bool is_new_vertex; + typename Vertices::iterator it_vertex; + ::CGAL::cpp0x::tie(it_vertex,is_new_vertex)=vertices.insert(std::make_pair(curr->vertex(),index)); + if (is_new_vertex) { + ++index; + points.push_back(curr->vertex()); + } + indices.push_back(it_vertex->second); + curr=curr->next(); + }while(curr!=start); + } + } + + void operator()( HDS& hds) { + ::CGAL::Polyhedron_incremental_builder_3 B( hds, true); + B.begin_surface( points.size(), facets.size() ); + for (typename std::list::iterator it=points.begin();it!=points.end();++it) + B.add_vertex((*it)->point()); + for (typename std::list< std::vector >::iterator + it=facets.begin();it!=facets.end();++it) + { + B.begin_facet(); + for (std::vector::iterator it_i=it->begin();it_i!=it->end();++it_i) + B.add_vertex_to_facet(*it_i); + B.end_facet(); + } + B.end_surface(); + } +}; + + + +template +void extract_connected_components( + const Polyhedron& P, + const Adjacency_criterium& adjacent, + CGAL::Union_find& uf, + Face_to_UF_handle_map& map_f2h, + Result& result + ) +{ + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + typedef ::CGAL::Union_find UF; + typedef typename UF::handle UF_handle; + typedef typename UF::iterator UF_iterator; + + CGAL_precondition(P.is_pure_triangle()); + +//init union-find: each facet is in its own set + for (typename Polyhedron::Facet_const_iterator it=P.facets_begin();it!=P.facets_end();++it){ + map_f2h.insert(std::make_pair(it,uf.make_set(it))); + } +//merge 2 facets if they share a common edge + for (typename Polyhedron::Facet_const_iterator it=P.facets_begin();it!=P.facets_end();++it){ + Facet_const_handle facet=it; + + UF_handle current=map_f2h.find(it)->second; + Halfedge_const_handle neighbors[3]; + neighbors[0]=facet->halfedge()->opposite(); + neighbors[1]=facet->halfedge()->next()->opposite(); + neighbors[2]=facet->halfedge()->next()->next()->opposite(); + + for (int i=0;i!=3;++i){ + if ( neighbors[i]->is_border_edge() ) continue; + UF_handle neigh=map_f2h.find(neighbors[i]->facet())->second; + if ( adjacent(neighbors[i]) && !uf.same_set(current,neigh) ){ + uf.unify_sets(current,neigh); + } + } + } + +//recover merged sets + for (UF_iterator it=uf.begin();it!=uf.end();++it){ + UF_handle master=uf.find(it); + result[*master].push_back(*it); + } +} + +template +void extract_connected_components(const Polyhedron& P,const Adjacency_criterium& adjacent,Output_iterator out) +{ + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef ::CGAL::Union_find UF; + typedef typename UF::handle UF_handle; + typedef typename UF::iterator UF_iterator; + typedef std::map,Compare_handle_ptr > Result; + typedef std::map > Facet_to_handle_map; + + UF uf; + Facet_to_handle_map map_f2h; + Result result; + + extract_connected_components(P,adjacent,uf,map_f2h,result); + + for (typename Result::iterator it=result.begin();it!=result.end();++it) + { + typedef std::list Facets; + const Facets& facets=it->second; + Polyhedron new_poly; + Build_polyhedron_subset modifier(new_poly,facets.begin(),facets.end()); + new_poly.delegate(modifier); + *out++=new_poly; + } +} + +template +void extract_connected_components(const Polyhedron& P,Output_iterator out) +{ + extract_connected_components(P,Dummy_true(),out); +} + +} } //namespace CGAL::internal + +#endif //CGAL_INTERNAL_POLYHEDRON_SUBSET_EXTRACTION_H diff --git a/Polyhedron/include/CGAL/internal/corefinement/intersection_coplanar_triangles_3.h b/Polyhedron/include/CGAL/internal/corefinement/intersection_coplanar_triangles_3.h new file mode 100644 index 00000000000..64341765ca0 --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/intersection_coplanar_triangles_3.h @@ -0,0 +1,327 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_INTERSECTION_COPLANAR_TRIANGLES_3_H +#define CGAL_INTERNAL_INTERSECTION_COPLANAR_TRIANGLES_3_H + +#include //for Intersection_type +#include +#include +#include + +//TODO rename this file when doing proper integration +namespace CGAL{ +namespace internal_IOP{ + + +//intersection point of two coplanar triangles that keeps track of +//the location of that point onto the triangles. +template +struct Intersection_point_with_info +{ + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + typedef IK Input_kernel; + typedef CGAL::Cartesian_converter Converter; + + Intersection_type type_1,type_2; //intersection type for 1st and 2nd facets + Halfedge_handle info_1,info_2; //halfedge providing primitive indicated by type_1 and type_2 + typename Exact_kernel::Point_3 point; //the geometric embedding of the intersection + + //constructor from a vertex of first triangle initialized inside the second triangle + Intersection_point_with_info(Halfedge_handle info1,Halfedge_handle info2): + type_1(VERTEX),type_2(FACET),info_1(info1),info_2(info2) + { + Converter converter; + point=converter(info_1->vertex()->point()); + } + + //constructor for intersection of edges. prev and curr are two points on an edge of the first facet (preserving the + //orientation of the facet). This edge is intersected by info2 from the second facet. + // + //The rational is the following: we first check whether curr and prev are on the same edge. I so we create + //an intersection point between two edges. Otherwise, the point is a vertex of the second facet included into + //the first facet. + // + //(V,F) : point initialy constructed + //(V,E) : (V,F) updated by get_orientation_and_update_info_2 (i.e lies on one edge) + //(V,V) : (V,E) updated by get_orientation_and_update_info_2 (i.e lies on two edges) + //(E,E) : created in the following function when prev and curr lie on the same edge + //(E,V) : (E,E) updated by get_orientation_and_update_info_2 (always done as lies on two edges) + //(E,F) : impossible + //(F,V) : detected when curr and prev and not on the same edge + //(F,E) : impossible + //(F,F) : impossible + // + Intersection_point_with_info(Intersection_point_with_info prev,Intersection_point_with_info curr, + Halfedge_handle info1,Halfedge_handle info2): + type_2(EDGE),info_2(info2) + { + #ifdef CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + std::cout << "prev: "; prev.print_debug(); + std::cout << "curr: "; curr.print_debug(); std::cout << std::endl; + #endif //CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + Converter converter; + if (prev.type_1==VERTEX && prev.info_1->next() == curr.info_1){ + CGAL_assertion(curr.type_1!=FACET); + type_1=EDGE; + info_1=curr.info_1; + } + else{ + if(curr.type_1==VERTEX && prev.info_1 == curr.info_1){ + CGAL_assertion(prev.type_1!=FACET); + type_1=EDGE; + info_1=curr.info_1; + } + else{ + if (curr.type_1==EDGE && prev.type_1==EDGE && curr.info_1==prev.info_1){ + type_1=EDGE; + info_1=curr.info_1; + } + else{ + //curr and prev are not on the same edge of the first facet. + //The intersection point to be computed is a VERTEX of the second facet + type_1=FACET; + info_1=info1; + type_2=VERTEX; + + //this is used to select the correct endpoint of the edge of the second facet + typename Exact_kernel::Collinear_3 is_collinear = Exact_kernel().collinear_3_object(); + if ( !is_collinear(prev.point,curr.point,converter(info_2->vertex()->point()) ) ){ + info_2=info_2->next()->next(); + CGAL_assertion( is_collinear(prev.point,curr.point,converter(info_2->vertex()->point()) ) ); + } + point = converter( info_2->vertex()->point() ); + return; + } + } + } + + //handle degenerate case when two edges overlap + //at least one of the two vertex has already been found as a vertex of a facet. Here we set it for the second point + if(prev.type_2!=FACET && curr.type_2!=FACET && (prev.type_1==VERTEX || prev.type_2==VERTEX) && (curr.type_1==VERTEX || curr.type_2==VERTEX)){ + typename Exact_kernel::Collinear_3 is_collinear = Exact_kernel().collinear_3_object(); + if ( is_collinear(prev.point,curr.point,converter(info_2->opposite()->vertex()->point()) ) ){ + info_2=info_2->next()->next(); + type_2=VERTEX; + point = converter( info_2->vertex()->point() ); + return; + } + if ( is_collinear(prev.point,curr.point,converter(info_2->vertex()->point()) ) ){ + type_2=VERTEX; + point = converter( info_2->vertex()->point() ); + return; + } + } + + //handle regular intersection of two edges + typename Exact_kernel::Construct_line_3 line_3=Exact_kernel().construct_line_3_object(); + typename Exact_kernel::Line_3 l1= + line_3(converter(info_2->vertex()->point()),converter(info_2->opposite()->vertex()->point())); + typename Exact_kernel::Line_3 l2=line_3(prev.point,curr.point); + CGAL::Object res=Exact_kernel().intersect_3_object()(l1,l2); + const typename Exact_kernel::Point_3* ptptr=CGAL::object_cast(&res); + CGAL_assertion(ptptr!=NULL); + point=*ptptr; + } + + void print_debug() const{ + std::cout << " ("; + if (type_1==VERTEX) std::cout << "V"; + if (type_1==EDGE) std::cout << "E"; + if (type_1==FACET) std::cout << "F"; + if (type_1==EMPTY) std::cout << "?"; + std::cout << "-" << &(*info_1); + std::cout << ";"; + if (type_2==VERTEX) std::cout << "V"; + if (type_2==EDGE) std::cout << "E"; + if (type_2==FACET) std::cout << "F"; + if (type_2==EMPTY) std::cout << "?"; + std::cout << ")" << "[" << CGAL::to_double(point.x()) << "," << CGAL::to_double(point.y()) << "," << CGAL::to_double(point.z()) << "]"; + } + + int debug_unique_type_int() const{ + int res=0; + switch (type_1){ + case VERTEX: res+=1; break; + case EDGE: res+=3; break; + case FACET: res+=7; break; + default: break; + } + switch (type_2){ + case VERTEX: res+=1; break; + case EDGE: res+=3; break; + case FACET: res+=7; break; + default: break; + } + return res; + } + + bool is_valid(Intersection_type type,Halfedge_handle info){ + bool valid=true; + Converter converter; + switch (type){ + case VERTEX: valid&= converter(info->vertex()->point())==point; break; + case EDGE:{ + typename Exact_kernel::Segment_3 seg= + Exact_kernel().construct_segment_3_object()( converter(info->vertex()->point()), + converter(info->opposite()->vertex()->point()) ); + valid&= Exact_kernel().has_on_3_object()(seg,point); + } + break; + case FACET:{ + typename Exact_kernel::Coplanar_orientation_3 orient=Exact_kernel().coplanar_orientation_3_object(); + typename Exact_kernel::Point_3 p=converter(info->vertex()->point()); + typename Exact_kernel::Point_3 q=converter(info->next()->vertex()->point()); + typename Exact_kernel::Point_3 r=converter(info->opposite()->vertex()->point()); + valid &= orient(p,q,r,point)==POSITIVE; + valid &= orient(q,r,p,point)==POSITIVE; + valid &= orient(r,p,q,point)==POSITIVE; + } + break; + default: valid=false; + } + return valid; + } + + bool is_valid(){ + return is_valid(type_1,info_1) && is_valid(type_2,info_2); + } + +}; + +template +CGAL::Orientation get_orientation_and_update_info_2(Halfedge_handle h,Inter_pt& p) +{ + typename Inter_pt::Exact_kernel::Coplanar_orientation_3 orient= + typename Inter_pt::Exact_kernel().coplanar_orientation_3_object(); + typename Inter_pt::Converter converter; + + CGAL::Orientation res = orient(converter(h->opposite()->vertex()->point()), + converter(h->vertex()->point()), + converter(h->next()->vertex()->point()), + p.point); + + if ( (p.type_1==VERTEX || p.type_1==EDGE) && res==COLLINEAR){ + if (p.type_2==FACET){ //detect a case (VERTEX,EDGE) + p.type_2=EDGE; + p.info_2=h; + } + else{ + //detect a case (VERTEX,VERTEX) or (EDGE,VERTEX) + CGAL_assertion(p.type_2==EDGE); + p.type_2=VERTEX; + if (p.info_2->next()!=h){ + CGAL_assertion(h->next()==p.info_2); + p.info_2=h; + } + } + } + + return res; +} + +template +void intersection_coplanar_facets_cutoff(Facet_handle f,std::list& inter_pts,Facet_handle other) +{ + #ifdef CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + std::cout << "cutoff: " << f->opposite()->vertex()->point() << " " << f->vertex()->point() << std::endl; + #endif //CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + if ( inter_pts.empty() ) return; + typedef typename std::list::iterator Iterator; + + std::map orientations; + for (Iterator it=inter_pts.begin();it!=inter_pts.end();++it){ + #ifdef CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + it->print_debug(); + #endif //CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + orientations[ &(*it) ]=get_orientation_and_update_info_2(f,*it); + } + #ifdef CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + std::cout << std::endl; + #endif //CGAL_DEBUG_COPLANAR_TRIANGLE_INTERSECTION + + int pt_added=0; + + Inter_pt* prev = &(*boost::prior(inter_pts.end())); + bool inter_pts_size_g_2 = inter_pts.size() > 2; + Iterator stop = inter_pts_size_g_2 ? inter_pts.end() : boost::prior(inter_pts.end()); + for (Iterator it=inter_pts.begin();it!=stop;++it) + { + Inter_pt* curr=&(*it); + if (!inter_pts_size_g_2) std::swap(prev,curr); + Orientation or_prev=orientations[prev],or_curr=orientations[curr]; + if ( (or_prev==POSITIVE && or_curr==NEGATIVE) || (or_prev==NEGATIVE && or_curr==POSITIVE) ) + { + Iterator it_curr = inter_pts_size_g_2 ? it:boost::next(it); + prev=&(* inter_pts.insert( it_curr,Inter_pt(*prev,*curr,other,f) ) ); + orientations[prev]=COLLINEAR; + ++pt_added; + } + prev=&(*it); + } + + CGAL_kernel_assertion(pt_added<3); + Iterator it=inter_pts.begin(); + std::size_t nb_interpt=inter_pts.size(); + //this boolean allows to reverse order of intersection points in case there were 3 remaining intersection points + //and the point in the middle was removed. In that case the order must be reversed to preserve the orientations + //of the last edge: + // A---X---B --> AB to be consistent with the other cases this should be BA! + // X---B---A --> BA + // B---A---X --> BA + // + bool should_revert_list=false; + + while(it!=inter_pts.end()) + { + if (orientations[&(*it)]==NEGATIVE){ + inter_pts.erase(it++); + if (--nb_interpt == 2 && it!=inter_pts.end() && boost::next(it)==inter_pts.end()) should_revert_list=true; + } + else + ++it; + } + if (should_revert_list && nb_interpt==2) inter_pts.reverse(); +} + +template +void +intersection_coplanar_facets( + Halfedge_handle f1, + Halfedge_handle f2, + std::list >& output ) +{ + typedef Intersection_point_with_info Inter_pt; + output.push_back( Inter_pt(f1,f2) ); + output.push_back( Inter_pt(f1->next(),f2) ); + output.push_back( Inter_pt(f1->next()->next(),f2) ); + + //intersect f2 with the three half planes which intersection defines f1 + intersection_coplanar_facets_cutoff(f2,output,f1); + intersection_coplanar_facets_cutoff(f2->next(),output,f1); + intersection_coplanar_facets_cutoff(f2->next()->next(),output,f1); +} + + + +} } //namespace CGAL::internal_IOP + + +#endif //CGAL_INTERNAL_INTERSECTION_COPLANAR_TRIANGLES_3_H + diff --git a/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3.h b/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3.h new file mode 100644 index 00000000000..f9ff0fccab5 --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3.h @@ -0,0 +1,166 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_H +#define CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_H + +//TODO rename this file when doing proper integration +#include +#include +namespace CGAL{ +namespace internal_IOP{ + +enum Intersection_type {FACET,EDGE,VERTEX,EMPTY,COPLNR}; + + +//define shortcut for intersection result types +template +struct Intersection_types{ + typedef typename Polyhedron_types::Halfedge_handle Intersection_info; + typedef cpp0x::tuple Intersection_result; +}; + + +template +typename Intersection_types::Intersection_result +find_intersection(const typename Kernel::Point_3& p, const typename Kernel::Point_3& q, + const typename Kernel::Point_3& a, const typename Kernel::Point_3& b, const typename Kernel::Point_3& c, + typename Polyhedron_types::Halfedge_handle /*hh*/, + typename Polyhedron_types::Facet_handle fh, + bool is_vertex_coplanar=false,bool is_vertex_opposite_coplanar=false) +{ + typedef typename Intersection_types::Intersection_info Intersection_info; + typedef typename Intersection_types::Intersection_result Intersection_result; + + Orientation ab=orientation(p,q,a,b); + Orientation bc=orientation(p,q,b,c); + Orientation ca=orientation(p,q,c,a); + + if ( ab==POSITIVE || bc==POSITIVE || ca==POSITIVE ) + return Intersection_result(EMPTY,Intersection_info(),false,false); + + int nb_coplanar=(ab==COPLANAR?1:0) + (bc==COPLANAR?1:0) + (ca==COPLANAR?1:0); + + if ( nb_coplanar==0 ) + return cpp0x::make_tuple(FACET,Intersection_info(fh->halfedge()),is_vertex_coplanar,is_vertex_opposite_coplanar); + + if (nb_coplanar==1){ + if (ab==COPLANAR) + return cpp0x::make_tuple(EDGE,Intersection_info(fh->halfedge()->next()),is_vertex_coplanar,is_vertex_opposite_coplanar); + if (bc==COPLANAR) + return cpp0x::make_tuple(EDGE,Intersection_info(fh->halfedge()->next()->next()),is_vertex_coplanar,is_vertex_opposite_coplanar); + CGAL_assertion(ca==COPLANAR); + return cpp0x::make_tuple(EDGE,Intersection_info(fh->halfedge()),is_vertex_coplanar,is_vertex_opposite_coplanar); + } + + CGAL_assertion(nb_coplanar==2); + + if (ab!=COPLANAR) + return cpp0x::make_tuple(VERTEX,Intersection_info(fh->halfedge()->next()->next()),is_vertex_coplanar,is_vertex_opposite_coplanar); + if (bc!=COPLANAR) + return cpp0x::make_tuple(VERTEX,Intersection_info(fh->halfedge()),is_vertex_coplanar,is_vertex_opposite_coplanar); + CGAL_assertion(ca!=COPLANAR); + return cpp0x::make_tuple(VERTEX,Intersection_info(fh->halfedge()->next()),is_vertex_coplanar,is_vertex_opposite_coplanar); +} + + +template +typename Intersection_types::Intersection_result +do_intersect(typename Polyhedron_types::Halfedge_handle hh, + typename Polyhedron_types::Facet_handle fh) +{ + typedef typename Intersection_types::Intersection_info Intersection_info; + typedef typename Intersection_types::Intersection_result Intersection_result; + + const typename Kernel::Point_3 & a = fh->halfedge()->vertex()->point(); + const typename Kernel::Point_3 & b = fh->halfedge()->next()->vertex()->point(); + const typename Kernel::Point_3 & c = fh->halfedge()->next()->next()->vertex()->point(); + const typename Kernel::Point_3 & p = hh->vertex()->point(); + const typename Kernel::Point_3 & q = hh->opposite()->vertex()->point(); + + + const Orientation abcp = orientation(a,b,c,p); + const Orientation abcq = orientation(a,b,c,q); + + + switch ( abcp ) { + case POSITIVE: + switch ( abcq ) { + case POSITIVE: + // the segment lies in the positive open halfspaces defined by the + // triangle's supporting plane + return Intersection_result(EMPTY,Intersection_info(),false,false); + case NEGATIVE: + // p sees the triangle in counterclockwise order + return find_intersection(p,q,a,b,c,hh,fh); + case COPLANAR: + // q belongs to the triangle's supporting plane + // p sees the triangle in counterclockwise order + return find_intersection(p,q,a,b,c,hh,fh,false,true); + + default: // should not happen. + CGAL_assertion(false); + return Intersection_result(EMPTY,Intersection_info(),false,false); + } + case NEGATIVE: + switch ( abcq ) { + case POSITIVE: + // q sees the triangle in counterclockwise order + return find_intersection(q,p,a,b,c,hh,fh); + case NEGATIVE: + // the segment lies in the negative open halfspaces defined by the + // triangle's supporting plane + return Intersection_result(EMPTY,Intersection_info(),false,false); + case COPLANAR: + // q belongs to the triangle's supporting plane + // p sees the triangle in clockwise order + return find_intersection(q,p,a,b,c,hh,fh,false,true); + default: // should not happen. + CGAL_assertion(false); + return Intersection_result(EMPTY,Intersection_info(),false,false); + } + case COPLANAR: // p belongs to the triangle's supporting plane + switch ( abcq ) { + case POSITIVE: + // q sees the triangle in counterclockwise order + return find_intersection(q,p,a,b,c,hh,fh,true,false); + case NEGATIVE: + // q sees the triangle in clockwise order + return find_intersection(p,q,a,b,c,hh,fh,true,false); + case COPLANAR: + // the segment is coplanar with the triangle's supporting plane + // we test whether the segment intersects the triangle in the common + // supporting plane + if ( ::CGAL::internal::do_intersect_coplanar(a,b,c,p,q,Kernel()) ) + return Intersection_result(COPLNR,Intersection_info(),true,true); + return Intersection_result(EMPTY,Intersection_info(),true,true); + + default: // should not happen. + CGAL_assertion(false); + return Intersection_result(EMPTY,Intersection_info(),false,false); + } + default: // should not happen. + CGAL_assertion(false); + return Intersection_result(EMPTY,Intersection_info(),false,false); + } +} + +}} //namespace CGAL::internal_IOP + +#endif //CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_H diff --git a/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3_coplanar.h b/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3_coplanar.h new file mode 100644 index 00000000000..21e7dbbab88 --- /dev/null +++ b/Polyhedron/include/CGAL/internal/corefinement/intersection_triangle_segment_3_coplanar.h @@ -0,0 +1,395 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_COPLANAR_H +#define CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_COPLANAR_H + + +namespace CGAL{ + namespace internal_IOP{ + +//enum Intersection_type {FACET,EDGE,VERTEX,EMPTY,COPLNR}; + +template +struct Intersection_point_coplanar{ + typedef typename Polyhedron_types::Halfedge_handle Halfedge_handle; + + Intersection_type type; + Halfedge_handle info; + Halfedge_handle segment_vertex; + + Intersection_point_coplanar(){}; + + Intersection_point_coplanar(Intersection_type type_, + Halfedge_handle info_, + Halfedge_handle segment_vertex_ + ) :type(type_),info(info_),segment_vertex(segment_vertex_){} +}; + +//test q vs segment bc +template +void point_vs_segment(Inter_pt_coplanar& pt, + Halfedge_handle bc, Halfedge_handle ab, + const Orientation& pqc, const Orientation& pqb) +{ + if (pqb==COLLINEAR){ + pt.type=VERTEX; + pt.info=ab; + } + else{ + if (pqc==COLLINEAR){ + pt.type=VERTEX; + pt.info=bc; + } + else{ + pt.type=EDGE; + pt.info=bc; + } + } +} + + +/* +// c +// / \ +// ___/___\____ with q-------p +// / \ +// a_______b +//supporting line of segment pq intersects ac and bc +*/ + +template +std::pair +decision_tree(const Point_3* a,const Point_3* b,const Point_3* c, + const Point_3* p,const Point_3* q, + const Orientation& pqa,const Orientation& pqb,const Orientation& pqc, + Halfedge_handle pq, + Halfedge_handle ca,Halfedge_handle ab,Halfedge_handle bc) +{ + CGAL_precondition(pqa!=NEGATIVE); + CGAL_precondition(pqb!=NEGATIVE); + CGAL_precondition(pqc!=POSITIVE); + + Inter_pt_coplanar pt1(EMPTY,NULL,NULL),pt2(EMPTY,NULL,NULL); + + //the segment supporting line intersects ac and bc + const Orientation bcq = coplanar_orientation(*b,*c,*q); + switch(bcq){ + case NEGATIVE: //segment does not intersect the triangle + return std::make_pair(pt1,pt2); + case COLLINEAR: //q \in [bc] + point_vs_segment(pt1,bc,ab,pqc,pqb); + pt1.segment_vertex=pq; + return std::make_pair(pt1,pt2); + default: + break; + } + + //===> q is to the left of bc + + const Orientation cap = coplanar_orientation(*c,*a,*p); + switch(cap){ + case NEGATIVE: //segment does not intersect the triangle + return false; + case COLLINEAR: //p \in [ca] + point_vs_segment(pt1,ca,bc,pqa,pqc); + pt1.segment_vertex=pq->opposite(); + return std::make_pair(pt1,pt2); + default: + break; + } + //===> p is to the right of ac + + + const Orientation cqa = coplanar_orientation(*c,*q,*a); + const Orientation cbp = coplanar_orientation(*c,*b,*p); + + if (pqc==COLLINEAR && cqa==POSITIVE && cbp==POSITIVE){ + //special case when c is inside the segment pq + pt1.type=VERTEX; + pt1.info=bc; + return std::make_pair(pt1,pt2); + } + + //where is located q? + switch (cqa){ + case POSITIVE: //q is outside the triangle + point_vs_segment(pt1,ca,bc,pqa,pqc); + break; + case NEGATIVE: //q is inside the triangle + pt1.type=FACET; + pt1.info=pq; + break; + case COLLINEAR: //q \in [ca] + point_vs_segment(pt1,ca,bc,pqa,pqc); + pt1.segment_vertex=pq; + break; + } + + //where is located p? + switch (cbp){ + case POSITIVE: //p is outside the triangle + point_vs_segment(pt2,bc,ab,pqc,pqb); + break; + case NEGATIVE: //p is inside the triangle + pt2.type=FACET; + pt2.info=pq->opposite(); + break; + case COLLINEAR: //p \in [bc] + point_vs_segment(pt2,bc,ab,pqc,pqb); + pt2.segment_vertex=pq->opposite(); + break; + } + return std::make_pair(pt1,pt2); +} + + +/* +// c +// / \ +// / \ +// / \ +// ____a_______b______ with q-------p +//supporting lines of segments pq and ab are the same. +*/ + +template +std::pair +collinear_decision_tree(const Point_3* a,const Point_3* b,const Point_3* c, + const Point_3* p,const Point_3* q, + const Orientation& pqa,const Orientation& pqb,const Orientation& pqc, + Halfedge_handle pq, + Halfedge_handle ca,Halfedge_handle ab,Halfedge_handle bc) +{ + CGAL_precondition(pqa==COLLINEAR); + CGAL_precondition(pqb==COLLINEAR); + CGAL_precondition(pqc==NEGATIVE); + + Inter_pt_coplanar pt1(EMPTY,NULL,NULL),pt2(EMPTY,NULL,NULL); + + //the segment supporting line intersects ac and bc + const Orientation bcq = coplanar_orientation(*b,*c,*q); + switch(bcq){ + case NEGATIVE: //segment does not intersect the triangle + return std::make_pair(pt1,pt2); + case COLLINEAR: // q = b + pt1.type=VERTEX; + pt1.info=ab; + pt1.segment_vertex=pq; + return std::make_pair(pt1,pt2); + default: + break; + } + + //===> q is to the left of b + + const Orientation cap = coplanar_orientation(*c,*a,*p); + switch(cap){ + case NEGATIVE: //segment does not intersect the triangle + return false; + case COLLINEAR: //p = a + pt1.type=VERTEX; + pt1.info=ca; + pt1.segment_vertex=pq->opposite(); + return std::make_pair(pt1,pt2); + default: + break; + } + + //===> p is to the right of a + + + const Orientation cqa = coplanar_orientation(*c,*q,*a); + const Orientation cbp = coplanar_orientation(*c,*b,*p); + + //where is located q? + switch (cqa){ + case POSITIVE: //q is to the left of a + pt1.type=VERTEX; + pt1.info=ca; + break; + case NEGATIVE: //q is to the right of a + pt1.type=EDGE; + pt1.info=ab; + pt1.segment_vertex=pq; + break; + case COLLINEAR: //q = a + pt1.type=VERTEX; + pt1.info=ca; + pt1.segment_vertex=pq; + break; + } + + //where is located p? + switch (cbp){ + case POSITIVE: //p is to the right of b + { + pt2.type=VERTEX; + pt2.info=ab; + } + break; + case NEGATIVE: //p is to the left of b + { + pt2.type=EDGE; + pt2.info=ab; + pt2.segment_vertex=pq->opposite(); + } + break; + case COLLINEAR: //p = b + pt2.type=VERTEX; + pt2.info=ab; + pt2.segment_vertex=pq->opposite(); + break; + } + + return std::make_pair(pt1,pt2); +} + +//std::pair,Intersection_point_coplanar > +template +std::pair,Intersection_point_coplanar > +do_intersect_coplanar(typename Polyhedron_types::Halfedge_handle pq, + typename Polyhedron_types::Face_handle fh) +{ + typedef Intersection_point_coplanar Inter_pt_coplanar; + typedef std::pair Return_type; + + + typedef typename Polyhedron_types::Halfedge_handle Halfedge_handle; + + const typename Kernel::Point_3 & A = fh->halfedge()->vertex()->point(); + const typename Kernel::Point_3 & B = fh->halfedge()->next()->vertex()->point(); + const typename Kernel::Point_3 & C = fh->halfedge()->next()->next()->vertex()->point(); + const typename Kernel::Point_3 & P = pq->vertex()->point(); + const typename Kernel::Point_3 & Q = pq->opposite()->vertex()->point(); + + const typename Kernel::Point_3 * a = &A; + const typename Kernel::Point_3 * b = &B; + const typename Kernel::Point_3 * c = &C; + + const typename Kernel::Point_3* p = &P; + const typename Kernel::Point_3* q = &Q; + + Halfedge_handle ca=fh->halfedge(); + Halfedge_handle ab=fh->halfedge()->next(); + Halfedge_handle bc=fh->halfedge()->next()->next(); + + + // Determine the orientation of the triangle in the common plane + + if (coplanar_orientation(A,B,C) != POSITIVE){ + // The triangle is not counterclockwise oriented swap two vertices. + b = &C; + c = &B; + std::swap(bc,ab); + } + + // Test whether the segment's supporting line intersects the + // triangle in the common plane + + Orientation pqa = coplanar_orientation(*p,*q,*a); + + //ensure pqa >= 0 + if (pqa == NEGATIVE){ + std::swap(p,q); + pqa=POSITIVE; + pq=pq->opposite(); + } + + const Orientation pqb = coplanar_orientation(*p,*q,*b); + const Orientation pqc = coplanar_orientation(*p,*q,*c); + + + + //Handle first case pq collinear with a triangle edge + if (pqa == COLLINEAR){ + if (pqb == COLLINEAR){ + //ab and pq are on the same line + if (pqc == NEGATIVE) + return collinear_decision_tree + (a,b,c,p,q,pqa,pqb,pqc,pq,ca,ab,bc); + else + return collinear_decision_tree + (a,b,c,q,p,pqa,pqb,NEGATIVE,pq->opposite(),ca,ab,bc); + } + if (pqc == COLLINEAR){ + //ac and pq are on the same line + if (pqb == NEGATIVE) + return collinear_decision_tree + (c,a,b,p,q,pqc,pqa,pqb,pq,bc,ca,ab); + else + return collinear_decision_tree + (c,a,b,q,p,pqc,pqa,NEGATIVE,pq->opposite(),bc,ca,ab); + } + } + else + if(pqb ==COLLINEAR && pqc == COLLINEAR){ + //bc and pq are on the same line + return collinear_decision_tree + (b,c,a,q,p,pqb,pqc,NEGATIVE,pq->opposite(),ab,bc,ca); + } + + + CGAL_assertion(pqa!=NEGATIVE); + + switch ( pqa ) { + case POSITIVE: + switch ( pqb ) { + case POSITIVE: + if (pqc == POSITIVE) + return std::make_pair(Inter_pt_coplanar(EMPTY,NULL,NULL),Inter_pt_coplanar(EMPTY,NULL,NULL)); + return decision_tree(a,b,c,p,q,pqa,pqb,pqc,pq,ca,ab,bc); + + case COLLINEAR: + case NEGATIVE: + if (pqc == POSITIVE) // b is isolated on the negative side + return decision_tree(c,a,b,p,q,pqc,pqa,pqb,pq,bc,ca,ab); + return decision_tree + (b,c,a,q,p,POSITIVE,opposite(pqc),NEGATIVE,pq->opposite(),ab,bc,ca); + default:// should not happen. + CGAL_assertion(false); + return std::make_pair(Inter_pt_coplanar(EMPTY,NULL,NULL),Inter_pt_coplanar(EMPTY,NULL,NULL)); + } + case COLLINEAR: + switch ( pqb ) { + case POSITIVE: + if (pqc == POSITIVE) // a is isolated on the negative side + return decision_tree(b,c,a,p,q,pqb,pqc,pqa,pq,ab,bc,ca); + return decision_tree(a,b,c,p,q,pqa,pqb,pqc,pq,ca,ab,bc); + + case NEGATIVE: + if (pqc == NEGATIVE) // a is isolated on the positive side + return decision_tree + (b,c,a,q,p,POSITIVE,POSITIVE,pqa,pq->opposite(),ab,bc,ca); + return decision_tree(c,a,b,p,q,pqc,pqa,pqb,pq,bc,ca,ab); + + default:// should not happen. + CGAL_assertion(false); + return std::make_pair(Inter_pt_coplanar(EMPTY,NULL,NULL),Inter_pt_coplanar(EMPTY,NULL,NULL)); + } + default:// should not happen. + CGAL_assertion(false); + return std::make_pair(Inter_pt_coplanar(EMPTY,NULL,NULL),Inter_pt_coplanar(EMPTY,NULL,NULL)); + + } +} + +}}//namespace CGAL::internal_IOP + +#endif //CGAL_INTERNAL_INTERSECTION_TRIANGLE_SEGMENT_3_COPLANAR_H diff --git a/Polyhedron/include/CGAL/intersection_of_Polyhedra_3.h b/Polyhedron/include/CGAL/intersection_of_Polyhedra_3.h new file mode 100644 index 00000000000..85182f11c56 --- /dev/null +++ b/Polyhedron/include/CGAL/intersection_of_Polyhedra_3.h @@ -0,0 +1,1958 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// Licensees holding a valid commercial license may use this file in +// accordance with the commercial license agreement provided with the software. +// +// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE +// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +// +// $URL$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERSECTION_OF_POLYHEDRA_3_H +#define CGAL_INTERSECTION_OF_POLYHEDRA_3_H + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#ifdef CGAL_COREFINEMENT_DEBUG +#warning look at CGAL/Mesh_3/Robust_intersection_traits.h and the statically filtered decision tree +#endif + +namespace CGAL{ + +//This functor computes the pairwise intersection of polyhedral surfaces. +//Intersection are given as a set of polylines +//The algorithm works as follow: +//From each polyhedral surface we can get it as a set of segments or as a set of triangles. +//We first use Box_intersection_d to filter intersection between all polyhedral +//surface segments and polyhedral triangles. +//From this filtered set, for each pair (segment,triangle), we look at the +//intersection type. If not empty, we can have three different cases +// 1)the segment intersect the interior of the triangle: +// We compute the intersection point and for each triangle incident +// to the segment, we write the fact that the point belong to the intersection +// of these two triangles. +// 2)the segment intersect the triangle on an edge +// We do the same thing as described above but +// for all triangle incident to the edge intersected +// 3)the segment intersect the triangle at a vertex +// for each edge incident to the vertex, we do +// the same operations as in 2) +// +//In case the segment intersect the triangle at one of the segment endpoint, +//we repeat the same procedure for each segment incident to this +//endpoint. +// +//Note that given a pair (segment,triangle)=(S,T), if S belongs +//to the plane of T, we have nothing to do in the following cases: +// -- no triangle T' contains S such that T and T' are coplanar +// -- at least one triangle contains S +// Indeed, the intersection points of S and T will be found using segments +// of T or segments adjacent to S. + + +namespace internal_IOP { + //an enum do decide which kind of intersection points are needed + struct No_predicates_on_constructions{}; + struct Predicates_on_constructions{}; +} // namespace internal_IOP + +template +struct Empty_node_visitor{ + typedef typename Polyhedron::Halfedge_const_handle Halfedge_handle; + void new_node_added(int,internal_IOP::Intersection_type,Halfedge_handle,Halfedge_handle,bool,bool){} + template + void annotate_graph(Iterator,Iterator){} + void update_terminal_nodes(std::vector&){} + void set_number_of_intersection_points_from_coplanar_facets(int){}; + void add_filtered_intersection(Halfedge_handle,Halfedge_handle,const Polyhedron&,const Polyhedron&){} + void start_new_polyline(int,int){} + void add_node_to_polyline(int){} + void new_input_polyhedron(const Polyhedron&){} + template + void finalize(T&){} + typedef internal_IOP::No_predicates_on_constructions Node_storage_type; + typedef Tag_true Is_polyhedron_const; + static const bool do_need_vertex_graph = false; +}; + +namespace internal_IOP{ + +template +struct Compare_handles{ + typedef typename Polyhedron_types::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron_types::Halfedge Halfedge; + typedef typename Polyhedron_types::Vertex Vertex; + typedef typename Polyhedron_types::Facet Facet; + typedef typename Polyhedron_types::Facet_handle Facet_handle; + typedef std::pair Vertex_handle_pair; + + static inline Vertex_handle_pair + make_sorted_pair_of_vertices(Halfedge_handle h) { + const Vertex* v1=&(* h->vertex() ); + const Vertex* v2=&(* h->opposite()->vertex() ); + if ( v1 < v2 ) + return Vertex_handle_pair(v1,v2); + return Vertex_handle_pair(v2,v1); + } + + bool operator()(Halfedge_handle h1,Halfedge_handle h2) const { + Vertex_handle_pair p1=make_sorted_pair_of_vertices(h1); + Vertex_handle_pair p2=make_sorted_pair_of_vertices(h2); + return p1 < p2; + } + + bool operator()(Facet_handle f1,Facet_handle f2) const { + return &(*f1) < &(*f2); + } + + bool operator()(const std::pair& p1, const std::pair& p2) const{ + Halfedge* h1= (std::min) ( &(*(p1.first)), &(*(p1.first->opposite())) ); + Halfedge* h2= (std::min) ( &(*(p2.first)), &(*(p2.first->opposite())) ); + return h1 +struct Compare_handle_pairs{ + typedef typename Polyhedron_types::Facet Facet; + typedef typename Polyhedron_types::Facet_handle Facet_handle; + typedef std::pair Facet_pair; + typedef std::pair Facet_pair_and_int; + + bool operator()(const Facet_pair& p1, const Facet_pair& p2) const{ + Facet* f1=&(*p1.first); + Facet* f2=&(*p2.first); + if (f1==f2){ + f1=&(*p1.second); + f2=&(*p2.second); + } + return f1f2) return false; + return p1.second +struct Order_along_a_halfedge{ + typedef typename Polyhedron_types::Halfedge_handle Halfedge_handle; + const Node_vector& nodes; + Halfedge_handle hedge; + + Order_along_a_halfedge(Halfedge_handle hedge_,const Node_vector& nodes_):nodes(nodes_),hedge(hedge_){} + bool operator()(int i,int j) const { + //returns true, iff q lies strictly between p and r. + try{ + typename Node_vector::Protector p; + return CGAL::collinear_are_strictly_ordered_along_line(nodes.to_interval(hedge->vertex()->point()), + nodes.interval_node(j), + nodes.interval_node(i)); + } + catch(CGAL::Uncertain_conversion_exception&){ + return CGAL::collinear_are_strictly_ordered_along_line(nodes.to_exact(hedge->vertex()->point()), + nodes.exact_node(j), + nodes.exact_node(i)); + } + } +}; + + +template +class Split_halfedge_at_point : public CGAL::Modifier_base { + typedef typename HDS::Halfedge_handle Halfedge_handle; + typedef typename HDS::Vertex_handle Vertex_handle; + typedef typename HDS::Vertex Vertex; + Halfedge_handle hedge; + Vertex vertex; + + typename HDS::Halfedge::Base* + unlock_halfedge(Halfedge_handle h){ + return static_cast(&(*h)); + } + +public: + + template + Split_halfedge_at_point( Halfedge_handle h,const Point_3& point):hedge(h),vertex(point){} + + void operator()( HDS& hds) { + + Vertex_handle v=hds.vertices_push_back(vertex); + Halfedge_handle opposite=hedge->opposite(); + + Halfedge_handle new_hedge=hds.edges_push_back(*hedge); + Halfedge_handle new_opposite=new_hedge->opposite(); + + //update next relations + unlock_halfedge(new_hedge)->set_next(hedge); + unlock_halfedge(new_hedge->prev())->set_next(new_hedge); + unlock_halfedge(hedge)->set_prev(new_hedge); + + unlock_halfedge(opposite)->set_next(new_opposite); + unlock_halfedge(new_opposite)->set_prev(opposite); + unlock_halfedge(new_opposite->next())->set_prev(new_opposite); + + unlock_halfedge(opposite)->set_vertex(v); + unlock_halfedge(new_hedge)->set_vertex(v); + + v->set_halfedge(new_hedge); + new_opposite->vertex()->set_halfedge(new_opposite); + } +}; + + + +} //namespace internal_IOP + + + +//WARNING THIS IS DONE ONLY FOR POLYHEDRON +template +class Node_visitor_for_polyline_split{ +//typedefs + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge Halfedge; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + //Info stores information about a particular intersection on an + //edge or a vertex of a polyhedron. The two first elements in + //the template describe the intersected simplex of the considered + //polyhedron; the two last elements describe the element of the + //second polyhedron (can be either a vertex, an edge of a facet) + //involved in the intersection + typedef CGAL::cpp0x::tuple Info; + typedef std::map Hedge_to_polyhedron_map; + typedef std::vector Infos; + typedef std::map Node_to_infos_map; +//data members + Node_to_infos_map node_infos; + Hedge_to_polyhedron_map hedge_to_polyhedron; + Halfedge_predicate is_on_polyline; + Set_vertex_corner set_as_corner; +//functions + void handle_principal_edge(int node_id, + internal_IOP::Intersection_type type, + Halfedge_handle principal_edge, + Halfedge_handle additional_edge, + bool is_vertex_coplanar, + bool is_vertex_opposite_coplanar) + { + bool coplanar_v=false; + if (is_vertex_coplanar) coplanar_v=true; + else if (is_vertex_opposite_coplanar){ + principal_edge=principal_edge->opposite(); + coplanar_v=true; + } + + if (coplanar_v) + handle_on_vertex(node_id,principal_edge,type,additional_edge); + else{ + if ( is_on_polyline(principal_edge) ){ + typename Node_to_infos_map::iterator it_res= + node_infos.insert(std::make_pair(node_id,Infos())).first; + it_res->second.push_back( Info(internal_IOP::EDGE,principal_edge,type,additional_edge) ); + } + } + } + + void handle_on_vertex(int node_id,Halfedge_handle edge,internal_IOP::Intersection_type type,Halfedge_handle additional_edge){ + Halfedge_handle current=edge; + do{ + if (is_on_polyline(current)){ + typename Node_to_infos_map::iterator it_res= + node_infos.insert(std::make_pair(node_id,Infos())).first; + it_res->second.push_back( Info(internal_IOP::VERTEX,current,type,additional_edge) ); + break; + } + current=current->next()->opposite(); + } + while(current!=edge); + } + + Halfedge* make_unique_key(Halfedge_handle h){ + if (&(*h) < &(*h->opposite())) + return &(*h); + else + return &(*h->opposite()); + } + + // new_hedge hedge + // -----------> -----------> + // v + // <----------- <----------- + // new_opposite opposite + // + void split_edge_and_retriangulate(Halfedge_handle hedge,const typename Kernel::Point_3& point,Polyhedron& P){ + internal_IOP::Split_halfedge_at_point delegated(hedge,point); + P.delegate( delegated ); + CGAL_assertion(P.is_valid()); + //triangulate the two adjacent facets + if (!hedge->is_border()) + P.split_facet(hedge->prev(),hedge->next()); + if (!hedge->opposite()->is_border()) + P.split_facet(hedge->opposite(),hedge->opposite()->next()->next()); + CGAL_assertion(P.is_valid()); + } + + //sort node ids so that we can split the hedge + //consecutively + template + void sort_vertices_along_hedge(std::vector& node_ids,Halfedge_handle hedge,const Node_vector& nodes) + { + std::sort(node_ids.begin(), + node_ids.end(), + internal_IOP::Order_along_a_halfedge(hedge,nodes) + ); + } + +public: + static const bool do_need_vertex_graph = false; + typedef internal_IOP::Predicates_on_constructions Node_storage_type; + typedef Tag_false Is_polyhedron_const; + + Node_visitor_for_polyline_split(){} + Node_visitor_for_polyline_split(const Halfedge_predicate& getting, + const Set_vertex_corner& setting) + :is_on_polyline(getting),set_as_corner(setting){} + + void new_node_added(int node_id, + internal_IOP::Intersection_type type, + Halfedge_handle principal_edge, + Halfedge_handle additional_edge, + bool is_vertex_coplanar, + bool is_vertex_opposite_coplanar) + { + switch(type) + { + case internal_IOP::FACET: //Facet intersected by an edge + handle_principal_edge(node_id,type,principal_edge,additional_edge,is_vertex_coplanar,is_vertex_opposite_coplanar); + break; + case internal_IOP::EDGE: //Edge intersected by an edge + handle_principal_edge(node_id,type,principal_edge,additional_edge,is_vertex_coplanar,is_vertex_opposite_coplanar); + if ( is_on_polyline(additional_edge) ){ + typename Node_to_infos_map::iterator it_res= + node_infos.insert(std::make_pair(node_id,Infos())).first; + it_res->second.push_back( Info(type,additional_edge, + ( is_vertex_coplanar||is_vertex_opposite_coplanar ) ? internal_IOP::VERTEX:internal_IOP::EDGE, + is_vertex_opposite_coplanar?principal_edge->opposite():principal_edge) ); + } + break; + case internal_IOP::VERTEX://Vertex intersected by an edge + handle_principal_edge(node_id,type,principal_edge,additional_edge,is_vertex_coplanar,is_vertex_opposite_coplanar); + handle_on_vertex( node_id,additional_edge, + ( is_vertex_coplanar||is_vertex_opposite_coplanar ) ? internal_IOP::VERTEX:internal_IOP::EDGE, + is_vertex_opposite_coplanar?principal_edge->opposite():principal_edge); + break; + default: + break; + } + } + + template + void annotate_graph(Iterator begin,Iterator end){ + for(Iterator it=begin;it!=end;++it){ + typename Node_to_infos_map::iterator it_res=node_infos.find(it->first); + if (it_res!=node_infos.end()) + it->second.make_terminal(); + } + } + + void update_terminal_nodes(std::vector& terminal_bools){ + for (typename Node_to_infos_map::iterator it=node_infos.begin();it!=node_infos.end();++it){ + terminal_bools[it->first]=true; + } + } + + void new_input_polyhedron(const Polyhedron&){} + void start_new_polyline(int,int){} + void add_node_to_polyline(int){} + void set_number_of_intersection_points_from_coplanar_facets(int){} + + void add_filtered_intersection(Halfedge_handle eh,Halfedge_handle fh,Polyhedron& Pe,Polyhedron& Pf){ + hedge_to_polyhedron.insert(std::make_pair(make_unique_key(eh),&Pe)); + hedge_to_polyhedron.insert(std::make_pair(make_unique_key(fh),&Pf)); + hedge_to_polyhedron.insert(std::make_pair(make_unique_key(fh->next()),&Pf)); + hedge_to_polyhedron.insert(std::make_pair(make_unique_key(fh->next()->next()),&Pf)); + } + + //split_halfedges + template + void finalize(const Node_vector& nodes){ + typedef std::map, + std::vector,internal_IOP::Compare_handles > Halfedges_to_split; + + Halfedges_to_split halfedges_to_split; + + for (typename Node_to_infos_map::iterator it=node_infos.begin();it!=node_infos.end();++it){ + int node_id=it->first; + const Infos& infos=it->second; + std::map > hedges_to_split; + + //collect information about halfedge to split + typename Infos::const_iterator it_info=infos.begin(); + for (;it_info!=infos.end();++it_info) + { + typename Hedge_to_polyhedron_map::iterator it_poly= + hedge_to_polyhedron.find(make_unique_key(CGAL::cpp0x::get<1>(*it_info))); + CGAL_assertion(it_poly!=hedge_to_polyhedron.end()); + //associate information to an intersection point: + //we give which simplex of the other polyhedron intersect the simplex considered + set_as_corner.add_info_to_node(node_id,it_poly->second,*it_info); + switch(CGAL::cpp0x::get<0>(*it_info)) + { + case internal_IOP::EDGE: + { + halfedges_to_split.insert( + std::make_pair( std::make_pair(CGAL::cpp0x::get<1>(*it_info),&(*(it_poly->second))),std::vector() ) + ).first->second.push_back(node_id); + break; + } + case internal_IOP::VERTEX: + set_as_corner(CGAL::cpp0x::get<1>(*it_info)->vertex(),node_id,it_poly->second); + break; + default: + CGAL_assertion(false); + //should never be here + } + } + } + + + //do the split + for(typename Halfedges_to_split::iterator it=halfedges_to_split.begin();it!=halfedges_to_split.end();++it){ + Halfedge_handle hedge=it->first.first; + Polyhedron* P=it->first.second; + std::vector& node_ids=it->second; + + sort_vertices_along_hedge(node_ids,hedge,nodes); + for (std::vector::iterator it_id=node_ids.begin();it_id!=node_ids.end();++it_id){ + split_edge_and_retriangulate(hedge,nodes[*it_id],*P); + set_as_corner(hedge->opposite()->vertex(),*it_id,(Polyhedron*)(0)); + } + } + } +}; + + +namespace internal_IOP{ + + template + typename Exact_kernel::Point_3 + compute_triangle_segment_intersection_point( + typename Polyhedron::Vertex_const_handle vh1,typename Polyhedron::Vertex_const_handle vh2, + typename Polyhedron::Vertex_const_handle vf1,typename Polyhedron::Vertex_const_handle vf2,typename Polyhedron::Vertex_const_handle vf3, + const Exact_kernel& ek) + { + CGAL::Cartesian_converter to_exact; + typename Exact_kernel::Triangle_3 t(to_exact( vf1->point() ), + to_exact( vf2->point() ), + to_exact( vf3->point() ) + ); + + typename Exact_kernel::Segment_3 s (to_exact( vh1->point() ), + to_exact( vh2->point() ) + ); + + typename Exact_kernel::Intersect_3 exact_intersect=ek.intersect_3_object(); + CGAL::Object inter=exact_intersect(t,s); + CGAL_assertion(CGAL::do_intersect(t,s)); + const typename Exact_kernel::Point_3* e_pt=CGAL::object_cast(&inter); + CGAL_assertion(e_pt!=NULL); + return *e_pt; + } + + template + typename Exact_kernel::Point_3 + compute_triangle_segment_intersection_point( + typename Polyhedron::Halfedge_const_handle edge, + typename Polyhedron::Facet_const_handle facet, + const Exact_kernel& ek) + { + return compute_triangle_segment_intersection_point( + edge->vertex(),edge->opposite()->vertex(), + facet->halfedge()->vertex(),facet->halfedge()->next()->vertex(),facet->halfedge()->opposite()->vertex(), + ek); + + } + + + //A class containing a vector of the intersection points. + //The third template parameter indicates whether an + //exact representation is required + template ::value> + class Triangle_segment_intersection_point; + + + //Store only the double version of the intersection points. + template + class Triangle_segment_intersection_point + { + //typedefs + typedef std::vector Node_vector; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_handle; + typedef typename Polyhedron::Facet_const_handle Facet_handle; + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + typedef CGAL::Cartesian_converter Exact_to_double; + //members + Node_vector nodes; + Exact_kernel ek; + Exact_to_double exact_to_double; + public: + typedef CGAL::Interval_nt::Protector Protector; + + const typename Kernel::Point_3& + operator[](int i) const { + return nodes[i]; + } + + const typename Kernel::Point_3& exact_node(int i) const {return nodes[i];} + const typename Kernel::Point_3& interval_node(int i) const {return nodes[i];} + const typename Kernel::Point_3& to_exact(const typename Kernel::Point_3& p) const {return p;} + const typename Kernel::Point_3& to_interval(const typename Kernel::Point_3& p) const {return p;} + + size_t size() const {return nodes.size();} + + //add a new node in the final graph. + //it is the intersection of the triangle with the segment + void add_new_node(Halfedge_handle edge,Facet_handle facet) + { + nodes.push_back ( exact_to_double( + compute_triangle_segment_intersection_point(edge,facet,ek) + )); + } + + void add_new_node(const typename Exact_kernel::Point_3& p) + { + nodes.push_back( exact_to_double(p) ); + } + + void add_new_node(const typename Kernel::Point_3& p) + { + nodes.push_back(p); + } + }; + + + //second specializations: store an exact copy of the points so that we can answer exactly predicates + //FYI, it used to have two specializations (one in the case the polyhedron + //can be edited and on if it cannot) building exact representation on demand. + //In the former case, we were using facet and halfedge while in the latter + //triple of vertex_handle and pair of vertex_handle + template + class Triangle_segment_intersection_point + { + //typedefs + public: + typedef CGAL::Simple_cartesian > Ikernel; + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + private: + typedef CGAL::Cartesian_converter Interval_to_double; + typedef CGAL::Cartesian_converter Double_to_interval; + typedef CGAL::Cartesian_converter Exact_to_interval; + typedef CGAL::Cartesian_converter Double_to_exact; + + typedef typename Polyhedron::Vertex_const_handle Vertex_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_handle; + typedef typename Polyhedron::Facet_const_handle Facet_handle; + + typedef std::vector Interval_nodes; + typedef std::vector Exact_nodes; + + + //members + Interval_nodes inodes; + Exact_nodes enodes; + + Interval_to_double interval_to_double; + Exact_to_interval exact_to_interval; + Double_to_interval double_to_interval; + Double_to_exact double_to_exact; + Exact_kernel ek; + + public: + typedef CGAL::Interval_nt::Protector Protector; + + typename Kernel::Point_3 + operator[](int i) const { + return interval_to_double(inodes[i]); + } + + const typename Ikernel::Point_3& + interval_node(int i) const { + return inodes[i]; + } + + typename Ikernel::Point_3 + to_interval(const typename Kernel::Point_3& p) const { + return double_to_interval(p); + } + + const Exact_kernel::Point_3 + exact_node(int i) const { + return enodes[i]; + } + + typename Exact_kernel::Point_3 + to_exact(const typename Kernel::Point_3& p) const { + return double_to_exact(p); + } + + + size_t size() const {return enodes.size();} + + void add_new_node(Halfedge_handle edge,Facet_handle facet) + { + enodes.push_back(compute_triangle_segment_intersection_point(edge,facet,ek) ); + inodes.push_back( exact_to_interval(enodes.back()) ); + } + + void add_new_node(const Exact_kernel::Point_3& p){ + enodes.push_back(p); + inodes.push_back( exact_to_interval(p) ); + } + + //the point is an input + void add_new_node(const typename Kernel::Point_3& p){ + enodes.push_back(to_exact(p)); + inodes.push_back( double_to_interval(p) ); + } + }; + + //Third specialization: The kernel already has exact constructions. + template + class Triangle_segment_intersection_point + { + //typedefs + typedef std::vector Node_vector; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_handle; + typedef typename Polyhedron::Facet_const_handle Facet_handle; + //members + Node_vector nodes; + Kernel k; + public: + typedef Kernel Ikernel; + typedef Kernel Exact_kernel; + typedef void* Protector; + const typename Kernel::Point_3& + operator[](int i) const { + return nodes[i]; + } + + size_t size() const {return nodes.size();} + const typename Kernel::Point_3& exact_node(int i) const {return nodes[i];} + const typename Kernel::Point_3& interval_node(int i) const {return nodes[i];} + + //add a new node in the final graph. + //it is the intersection of the triangle with the segment + void add_new_node(Halfedge_handle edge,Facet_handle facet) + { + nodes.push_back ( + compute_triangle_segment_intersection_point(edge,facet,k) + ); + } + + void add_new_node(const typename Kernel::Point_3& p) + { + nodes.push_back(p); + } + + const typename Kernel::Point_3& to_interval(const typename Kernel::Point_3& p) const { return p; } + const typename Kernel::Point_3& to_exact(const typename Kernel::Point_3& p) const { return p; } + + }; + +} + +//TODO an important requirement is that the Polyhedron should be based on a list-based +//HDS. We use a lot of maps that use the address of Facet,Halfedge and a reallocation would +//be dramatic. + +template< class Polyhedron, + class Kernel=typename Polyhedron::Traits::Kernel, + class Node_visitor=Empty_node_visitor, + class Node_storage_type=typename Node_visitor::Node_storage_type, + class Use_const_polyhedron=typename Node_visitor::Is_polyhedron_const + > +class Intersection_of_Polyhedra_3{ + +//typedefs + typedef typename Kernel::Triangle_3 Triangle; + typedef typename Kernel::Segment_3 Segment; + typedef internal_IOP:: + Polyhedron_types Polyhedron_types; + + typedef typename Polyhedron_types::Polyhedron_ref Polyhedron_ref; + typedef typename Polyhedron_types::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron_types::Halfedge_iterator Halfedge_iterator; + typedef typename Polyhedron_types::Facet_iterator Facet_iterator; + typedef typename Polyhedron_types::Facet_handle Facet_handle; + typedef typename Polyhedron_types::Vertex_handle Vertex_handle; + typedef typename Polyhedron_types::Vertex Vertex; + typedef typename Polyhedron_types::Facet Facet; + typedef CGAL::Box_intersection_d::Box_with_handle_d< + double, 3, Halfedge_handle> Box; + typedef std::pair Facet_pair; + typedef std::pair Facet_pair_and_int; + + typedef internal_IOP:: + Compare_handles Compare_handles; + + typedef internal_IOP:: + Compare_handle_pairs Compare_handle_pairs; + + + typedef std::map,Compare_handle_pairs> Facets_to_nodes_map;//Indeed the boundary of the intersection of two coplanar triangles may contain several segments. + typedef std::set Coplanar_facets_set;//any insertion should be done with make_sorted_pair_of_facets + typedef typename Kernel::Point_3 Node; + typedef internal_IOP::Triangle_segment_intersection_point + Node_vector; + + typedef typename internal_IOP:: + Intersection_types + ::Intersection_result Intersection_result; + + typedef std::set Facet_set; + + typedef std::map + Edge_to_intersected_facets; + #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES + typedef std::set Coplanar_duplicated_intersection_set; + #endif +//helper functions + static inline Facet_pair + make_sorted_pair_of_facets(Facet_handle fh1,Facet_handle fh2) { + const Facet* f1=&(*fh1); + const Facet* f2=&(*fh2); + if (f1 < f2) + return Facet_pair(fh1,fh2); + return Facet_pair(fh2,fh1); + } + + static inline Facet_pair_and_int + make_sorted_pair_of_facets_with_int(Facet_handle fh1,Facet_handle fh2,int i) { + return std::make_pair(make_sorted_pair_of_facets(fh1,fh2),i); + } + + static inline std::pair make_sorted_void_pair(void* v1,void* v2){ + if (v1opposite()) ) return h; + return h->opposite(); + } + +//member variables + Edge_to_intersected_facets edge_to_sfacet; //Associate a segment to a filtered set of facets that may be intersected + Facets_to_nodes_map f_to_node; //Associate a pair of triangle to their intersection points + Coplanar_facets_set coplanar_facets;//Contains all pairs of triangular facets intersecting that are coplanar + Node_vector nodes; //Contains intersection points of polyhedra + Node_visitor* visitor; + bool is_default_visitor; //indicates whether the visitor need to be deleted + #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES + //this does not occur only when one extremity of an edge is inside a face. + // The problem occur every time an edge or a part of an edge with two incident triangles + // is on the intersection polyline, I choose to direcly filter the output by removing duplicated edges + Coplanar_duplicated_intersection_set coplanar_duplicated_intersection;//Set containing edges that are duplicated because of edges (partially) included in a triangle + #endif + +//functions that should come from a traits class + bool has_at_least_two_incident_faces(Halfedge_handle edge) + { + return !edge->is_border_edge(); + } + + template + void get_incident_facets(Halfedge_handle edge,Output_iterator out){ + if (!edge->is_border()) *out++=edge->facet(); + if (!edge->opposite()->is_border()) *out++=edge->opposite()->facet(); + } + + template + void get_incident_edges_to_vertex(Halfedge_handle edge,Output_iterator out){ + Halfedge_handle current=edge; + do{ + *out++=current; + current=current->next()->opposite(); + } + while(current!=edge); + } + +//internal functions + + class Map_edge_facet_bbox_intersection { + Edge_to_intersected_facets& edge_to_sfacet; + Polyhedron_ref polyhedron_triangle; + Polyhedron_ref polyhedron_edge; + Node_visitor& visitor; + public: + Map_edge_facet_bbox_intersection(Edge_to_intersected_facets& map_, + Polyhedron_ref P, + Polyhedron_ref Q, + Node_visitor& visitor_) + :edge_to_sfacet(map_),polyhedron_triangle(P),polyhedron_edge(Q),visitor(visitor_){} + + void operator()( const Box* fb, const Box* eb) const { + Halfedge_handle fh = fb->handle(); + Halfedge_handle eh = eb->handle(); + + typename Edge_to_intersected_facets::iterator res= + edge_to_sfacet.insert(std::make_pair(eh,Facet_set())).first; + res->second.insert(fh->facet()); + visitor.add_filtered_intersection(eh,fh,polyhedron_edge,polyhedron_triangle); + } + }; + + class Map_edge_facet_bbox_intersection_extract_coplanar { + Edge_to_intersected_facets& edge_to_sfacet; + Coplanar_facets_set& coplanar_facets; + Polyhedron_ref polyhedron_triangle; + Polyhedron_ref polyhedron_edge; + Node_visitor& visitor; + public: + Map_edge_facet_bbox_intersection_extract_coplanar( + Edge_to_intersected_facets& map_, + Coplanar_facets_set& coplanar_facets_, + Polyhedron_ref P, + Polyhedron_ref Q, + Node_visitor& visitor_) + :edge_to_sfacet(map_),coplanar_facets(coplanar_facets_),polyhedron_triangle(P),polyhedron_edge(Q),visitor(visitor_) + {} + + void operator()( const Box* fb, const Box* eb) const { + Halfedge_handle fh = fb->handle(); + Halfedge_handle eh = eb->handle(); + + + //check if the segment intersects the plane of the facet or if it is included in the plane + const typename Kernel::Point_3 & a = fh->vertex()->point(); + const typename Kernel::Point_3 & b = fh->next()->vertex()->point(); + const typename Kernel::Point_3 & c = fh->next()->next()->vertex()->point(); + const Orientation abcp = orientation(a,b,c,eh->vertex()->point()); + const Orientation abcq = orientation(a,b,c,eh->opposite()->vertex()->point()); + if (abcp==abcq){ + if (abcp!=COPLANAR){ +// std::cout << "rejected " << &(*fh->facet()) << "{" << &(*eh->facet()) << " " <<&(*eh->opposite()->facet()) << " "<< eh->vertex()->point() << " " << eh->opposite()->vertex()->point() << "}" <is_border() && orientation(a,b,c,eh->next()->vertex()->point())==COPLANAR){ + coplanar_facets.insert(make_sorted_pair_of_facets(eh->facet(),fh->facet())); + } + if (!eh->opposite()->is_border() && orientation(a,b,c,eh->opposite()->next()->vertex()->point())==COPLANAR){ + coplanar_facets.insert(make_sorted_pair_of_facets(eh->opposite()->facet(),fh->facet())); + } + visitor.add_filtered_intersection(eh,fh,polyhedron_edge,polyhedron_triangle); + //in case only the edge is coplanar, the intersection points will be detected using an incident facet + //(see remark at the beginning of the file) + return; + } + + typename Edge_to_intersected_facets::iterator res= + edge_to_sfacet.insert(std::make_pair(eh,Facet_set())).first; + res->second.insert(fh->facet()); + visitor.add_filtered_intersection(eh,fh,polyhedron_edge,polyhedron_triangle); + } + }; + + // This function tests the intersection of the faces of P with the edges of Q + void filter_intersections( Polyhedron_ref P, Polyhedron_ref Q) { + std::vector facet_boxes, edge_boxes; + facet_boxes.reserve( P.size_of_facets()); + for ( Facet_iterator i = P.facets_begin(); i != P.facets_end(); ++i){ + facet_boxes.push_back( + Box( i->halfedge()->vertex()->point().bbox() + + i->halfedge()->next()->vertex()->point().bbox() + + i->halfedge()->next()->next()->vertex()->point().bbox(), + i->halfedge())); + } + std::vector facet_box_ptr; + facet_box_ptr.reserve( P.size_of_facets()); + for ( typename std::vector::iterator j = facet_boxes.begin(); j != facet_boxes.end(); ++j){ + facet_box_ptr.push_back( &*j); + } + + for ( Halfedge_iterator i = Q.halfedges_begin(); i != Q.halfedges_end(); ++i){ + if(&*i < &*(i->opposite())){ + edge_boxes.push_back( + Box( i->vertex()->point().bbox() + + i->opposite()->vertex()->point().bbox(), + i)); + } + } + + std::vector edge_box_ptr; + edge_box_ptr.reserve( Q.size_of_halfedges()/2); + for ( typename std::vector::iterator j = edge_boxes.begin(); j != edge_boxes.end(); ++j){ + edge_box_ptr.push_back( &*j); + } + + CGAL::box_intersection_d( facet_box_ptr.begin(), facet_box_ptr.end(), + edge_box_ptr.begin(), edge_box_ptr.end(), + #ifdef DO_NOT_HANDLE_COPLANAR_FACETS + Map_edge_facet_bbox_intersection(edge_to_sfacet,P,Q,*visitor), + #else + Map_edge_facet_bbox_intersection_extract_coplanar(edge_to_sfacet,coplanar_facets,P,Q,*visitor), + #endif + std::ptrdiff_t(2000) + ); + } + + + void add_intersection_point_to_facet_and_all_edge_incident_facets(Facet_handle facet, + Halfedge_handle edge, + int node_id) + { + std::vector incident_facets; + get_incident_facets(edge,std::back_inserter(incident_facets)); + for (typename std::vector::iterator it=incident_facets.begin(); + it!=incident_facets.end();++it) + { + CGAL_assertion(cgal_do_intersect_debug(facet,*it)); + + Facet_pair facet_pair = make_sorted_pair_of_facets(facet,*it); + if ( !coplanar_facets.empty() && coplanar_facets.find(facet_pair)!=coplanar_facets.end() ) continue; + typename Facets_to_nodes_map::iterator it_list= + f_to_node.insert( std::make_pair( Facet_pair_and_int(facet_pair,0),std::set()) ).first; + it_list->second.insert(node_id); + } + } + + void cip_handle_case_edge(int node_id, + Facet_set* fset, + Halfedge_handle edge, + Halfedge_handle edge_intersected) + { + //associate the intersection point to all facets incident to the intersected edge using edge + std::vector incident_facets; + get_incident_facets(edge_intersected,std::back_inserter(incident_facets)); + for (typename std::vector::iterator it=incident_facets.begin(); + it!=incident_facets.end();++it) + { + add_intersection_point_to_facet_and_all_edge_incident_facets(*it,edge,node_id); + if (fset!=NULL) fset->erase(*it); + } + incident_facets.clear(); + + //associate the intersection point to all facets incident to edge using the intersected edge + //at least one pair of facets is already handle above + + typename Edge_to_intersected_facets::iterator it_fset=edge_to_sfacet.find(edge_intersected); + if (it_fset==edge_to_sfacet.end()) return; + Facet_set& fset_bis=it_fset->second; + get_incident_facets(edge,std::back_inserter(incident_facets)); + for (typename std::vector::iterator it=incident_facets.begin(); + it!=incident_facets.end();++it) + { +// add_intersection_point_to_facet_and_all_edge_incident_facets(*it,edge_intersected,node_id); //this call is not needed, already done in the first loop + fset_bis.erase(*it); + } + } + + void cip_handle_case_vertex(int node_id, + Facet_set* fset, + Halfedge_handle edge, + Halfedge_handle vertex_intersected) + { + std::vector incident_halfedges; + get_incident_edges_to_vertex(vertex_intersected,std::back_inserter(incident_halfedges)); + for (typename std::vector::iterator + it=incident_halfedges.begin();it!=incident_halfedges.end();++it) + { + cip_handle_case_edge(node_id,fset,edge,*it); + } + } + + //add a new node in the final graph. + //it is the intersection of the triangle with the segment + void add_new_node(Halfedge_handle edge, + Facet_handle facet, + const Intersection_result& inter_res, + Node_vector& nodes) + { + bool is_vertex_coplanar = CGAL::cpp0x::get<2>(inter_res); + if (is_vertex_coplanar) + nodes.add_new_node(edge->vertex()->point()); + else{ + bool is_opposite_vertex_coplanar = CGAL::cpp0x::get<3>(inter_res); + if (is_opposite_vertex_coplanar) + nodes.add_new_node(edge->opposite()->vertex()->point()); + else + nodes.add_new_node(edge,facet); + } + } + + //either the exact or input point can be used to create a node + //with this function + template + void add_new_node(const Point& pt) + { + nodes.add_new_node(pt); + } + + + #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES + void check_coplanar_edge(Halfedge_handle hedge,Facet_handle facet) + { + const typename Kernel::Point_3& p0=facet->halfedge()->vertex()->point(); + const typename Kernel::Point_3& p1=facet->halfedge()->next()->vertex()->point(); + const typename Kernel::Point_3& p2=facet->halfedge()->opposite()->vertex()->point(); + CGAL_precondition( orientation( p0,p1,p2,hedge->vertex()->point() ) == COPLANAR ); + + if ( has_at_least_two_incident_faces(hedge) && orientation( p0,p1,p2,hedge->opposite()->vertex()->point() ) == COPLANAR ) + { + //In case two facets are incident along such this edge, the intersection + //will be reported twice. We keep track of this so that at the end, we can remove one intersecting edge out of the two + //choose the smaller of the two faces (only one need to de deleted) + Facet_handle smaller=make_sorted_pair_of_facets(hedge->face(),hedge->opposite()->face()).first; + coplanar_duplicated_intersection.insert(make_sorted_pair_of_facets(smaller,facet)); + } + } + + bool are_incident_facets_coplanar(Halfedge_handle hedge){ + const typename Kernel::Point_3& p0=hedge->vertex()->point(); + const typename Kernel::Point_3& p1=hedge->next()->vertex()->point(); + const typename Kernel::Point_3& p2=hedge->opposite()->vertex()->point(); + const typename Kernel::Point_3& p3=hedge->opposite()->next()->vertex()->point(); + return orientation( p0,p1,p2,p3 ) == COPLANAR; + } + + void check_coplanar_edge(Halfedge_handle hedge,Halfedge_handle additional_edge,internal_IOP::Intersection_type type) + { + switch(type){ + case internal_IOP::FACET: + check_coplanar_edge(hedge,additional_edge->face()); + break; + + case internal_IOP::EDGE: + if ( !additional_edge->is_border() ){ + check_coplanar_edge(hedge,additional_edge->face()); + } + if (!additional_edge->opposite()->is_border()) + check_coplanar_edge(hedge,additional_edge->opposite()->face()); + break; + case internal_IOP::VERTEX: + { + //consider all incident faces + Halfedge_handle current=additional_edge; + do{ + if( !current->is_border() ) + check_coplanar_edge(hedge,current->face()); + current=current->next()->opposite(); + } + while(current!=additional_edge); + } + break; + + default: + CGAL_assertion(type==internal_IOP::COPLNR); + break; + } + } + + void check_coplanar_edge_old(Halfedge_handle hedge,Halfedge_handle additional_edge,internal_IOP::Intersection_type type) + { + switch(type){ + case internal_IOP::FACET: + check_coplanar_edge(hedge,additional_edge->face()); + break; + + case internal_IOP::EDGE: + { + if ( !additional_edge->is_border() ){ + if (!additional_edge->opposite()->is_border()){ + if ( are_incident_facets_coplanar(additional_edge) ) + { + Facet_handle facet=additional_edge->face(); + const typename Kernel::Point_3& p0=facet->halfedge()->vertex()->point(); + const typename Kernel::Point_3& p1=facet->halfedge()->next()->vertex()->point(); + const typename Kernel::Point_3& p2=facet->halfedge()->opposite()->vertex()->point(); + CGAL_precondition( orientation( p0,p1,p2,hedge->vertex()->point() ) == COPLANAR ); + + if ( has_at_least_two_incident_faces(hedge) && orientation( p0,p1,p2,hedge->opposite()->vertex()->point() ) == COPLANAR ) + { + //In case two facets are incident along a common edge of two coplanar triangles. + //We need to remove three out of the four reported pair + Facet_handle smaller=make_sorted_pair_of_facets(hedge->face(),hedge->opposite()->face()).first; + coplanar_duplicated_intersection.insert(make_sorted_pair_of_facets(hedge->face(),facet)); + coplanar_duplicated_intersection.insert(make_sorted_pair_of_facets(hedge->opposite()->face(),facet)); + coplanar_duplicated_intersection.insert(make_sorted_pair_of_facets(hedge->opposite()->face(),additional_edge->opposite()->face())); + } + } + else + { + check_coplanar_edge(hedge,additional_edge->face()); + check_coplanar_edge(hedge,additional_edge->opposite()->face()); + } + } + else + check_coplanar_edge(hedge,additional_edge->face()); + } + else{ + CGAL_assertion(!additional_edge->opposite()->is_border()); + check_coplanar_edge(hedge,additional_edge->opposite()->face()); + } + } + break; + case internal_IOP::VERTEX: + + break; + + default: + CGAL_assertion(type==internal_IOP::COPLNR); + break; + } + } + + template + void check_coplanar_edges(Hedge_iterator begin,Hedge_iterator end,Halfedge_handle additional_edge,internal_IOP::Intersection_type type) + { + for (Hedge_iterator it=begin;it!=end;++it) + check_coplanar_edge(*it,additional_edge,type); + } + #endif + + void print_type_debug(internal_IOP::Intersection_type type,bool cpl,bool opp_cpl) + { + switch(type){ + case internal_IOP::COPLNR: + std::cout << "COPLNR " << cpl << " " << opp_cpl << std::endl; + break; + case internal_IOP::EMPTY: + std::cout << "EMPTY " << cpl << " " << opp_cpl << std::endl; + break; + + case internal_IOP::FACET: + std::cout << "FACET " << cpl << " " << opp_cpl << std::endl; + break; + + case internal_IOP::EDGE: + std::cout << "EDGE " << cpl << " " << opp_cpl << std::endl; + break; + case internal_IOP::VERTEX: + std::cout << "VERTEX " << cpl << " " << opp_cpl << std::endl; + break; + } + } + + + void handle_coplanar_case_VERTEX_FACET(Halfedge_handle vertex,Halfedge_handle facet,int node_id){ + visitor->new_node_added(node_id,internal_IOP::FACET,vertex,facet,true,false); + std::vector all_edges; + get_incident_edges_to_vertex(vertex,std::back_inserter(all_edges)); + typename std::vector::iterator it_edge=all_edges.begin(); + for (;it_edge!=all_edges.end();++it_edge){ + add_intersection_point_to_facet_and_all_edge_incident_facets(facet->facet(),*it_edge,node_id); + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + if (it_ets!=edge_to_sfacet.end()) it_ets->second.erase(facet->facet()); + } + } + + void handle_coplanar_case_VERTEX_EDGE(Halfedge_handle vertex,Halfedge_handle edge,int node_id){ + visitor->new_node_added(node_id,internal_IOP::VERTEX,edge,vertex,false,false); + std::vector all_edges; + get_incident_edges_to_vertex(vertex,std::back_inserter(all_edges)); + typename std::vector::iterator it_edge=all_edges.begin(); + for (;it_edge!=all_edges.end();++it_edge){ + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + Facet_set* fset = (it_ets!=edge_to_sfacet.end())?&(it_ets->second):NULL; + cip_handle_case_edge(node_id,fset,*it_edge,edge); + } + } + + void handle_coplanar_case_VERTEX_VERTEX(Halfedge_handle vertex1,Halfedge_handle vertex2,int node_id){ + visitor->new_node_added(node_id,internal_IOP::VERTEX,vertex2,vertex1,true,false); + std::vector all_edges; + get_incident_edges_to_vertex(vertex1,std::back_inserter(all_edges)); + typename std::vector::iterator it_edge=all_edges.begin(); + for (;it_edge!=all_edges.end();++it_edge){ + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + Facet_set* fset = (it_ets!=edge_to_sfacet.end())?&(it_ets->second):NULL; + cip_handle_case_vertex(node_id,fset,*it_edge,vertex2); + } + } + + + template + int get_or_create_node(Cpl_inter_pt& ipt,int& current_node,Coplanar_node_map& coplanar_node_map){ + void *v1, *v2; + switch(ipt.type_1){ + case internal_IOP::VERTEX: v1=(void*)( &(*(ipt.info_1->vertex())) ); break; + case internal_IOP::EDGE : v1=(void*)( &(*smaller_handle(ipt.info_1)) ); break; + case internal_IOP::FACET : v1=(void*)( &(*(ipt.info_1->facet())) ); break; + default: CGAL_assertion(!"Should not get there!"); + } + + switch(ipt.type_2){ + case internal_IOP::VERTEX: v2=(void*)( &(*(ipt.info_2->vertex())) ); break; + case internal_IOP::EDGE : v2=(void*)( &(*smaller_handle(ipt.info_2)) ); break; + case internal_IOP::FACET : v2=(void*)( &(*(ipt.info_2->facet())) ); break; + default: CGAL_assertion(!"Should not get there!"); + } + + std::pair key=make_sorted_void_pair(v1,v2); + + std::pair res=coplanar_node_map.insert(std::make_pair(key,current_node+1)); + if (res.second){ //insert a new node + + if (ipt.type_1==internal_IOP::VERTEX) + add_new_node(ipt.info_1->vertex()->point()); + else{ + if(ipt.type_2==internal_IOP::VERTEX) + add_new_node(ipt.info_2->vertex()->point()); + else + add_new_node(ipt.point); + } + return ++current_node; + } + return res.first->second; + } + + void compute_intersection_of_coplanar_facets(int& current_node){ + typedef std::map,int> Coplanar_node_map; + Coplanar_node_map coplanar_node_map; + + for (typename Coplanar_facets_set::iterator it=coplanar_facets.begin();it!=coplanar_facets.end();++it){ + Facet_handle f1=it->first; + Facet_handle f2=it->second; + typedef internal_IOP::Intersection_point_with_info Cpl_inter_pt; + std::list inter_pts; + internal_IOP::intersection_coplanar_facets(f1->halfedge(),f2->halfedge(),inter_pts); +// std::cout << "found " << inter_pts.size() << " inter pts: "; + std::size_t nb_pts=inter_pts.size(); + std::vector cpln_nodes; cpln_nodes.reserve(nb_pts); + for (typename std::list::iterator iti=inter_pts.begin();iti!=inter_pts.end();++iti){ + #ifdef CGAL_COREFINEMENT_DEBUG + //iti->print_debug(); + #endif + CGAL_assertion(iti->is_valid()); + int node_id=get_or_create_node(*iti,current_node,coplanar_node_map); + cpln_nodes.push_back(node_id); + + switch(iti->type_1){ + case internal_IOP::VERTEX: + { + switch(iti->type_2){ + case internal_IOP::VERTEX: handle_coplanar_case_VERTEX_VERTEX(iti->info_1,iti->info_2,node_id); break; + case internal_IOP::EDGE : handle_coplanar_case_VERTEX_EDGE(iti->info_1,iti->info_2,node_id); break; + case internal_IOP::FACET : handle_coplanar_case_VERTEX_FACET(iti->info_1,iti->info_2,node_id); break; + default: CGAL_assertion(!"Should not get there!"); + } + } + break; + case internal_IOP::EDGE:{ + switch(iti->type_2){ + case internal_IOP::VERTEX:handle_coplanar_case_VERTEX_EDGE(iti->info_2,iti->info_1,node_id);break; + case internal_IOP::EDGE: + { + visitor->new_node_added(node_id,internal_IOP::EDGE,iti->info_1,iti->info_2,false,false); + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(iti->info_1); + Facet_set* fset = (it_ets!=edge_to_sfacet.end())?&(it_ets->second):NULL; + cip_handle_case_edge(node_id,fset,iti->info_1,iti->info_2); + } + break; + default: CGAL_assertion(!"Should not get there!"); + } + } + break; + + case internal_IOP::FACET: + { + CGAL_assertion(iti->type_2==internal_IOP::VERTEX); + handle_coplanar_case_VERTEX_FACET(iti->info_2,iti->info_1,node_id); + } + break; + + default: CGAL_assertion(!"Should not get there!"); + } + } + switch (nb_pts){ + case 0: break; + case 1: + { + typename Facets_to_nodes_map::iterator it_list= + f_to_node.insert( std::make_pair( Facet_pair_and_int(*it,1),std::set()) ).first; + it_list->second.insert(cpln_nodes[0]); + } + break; + default: + { + int i=0; + std::size_t stop=nb_pts + (nb_pts<3?-1:0); + for (std::size_t k=0;k()) ).first; + it_list->second.insert( cpln_nodes[k] ); + it_list->second.insert( cpln_nodes[(k+1)%nb_pts] ); + } + } + } +// std::cout << std::endl; + } + } + + void compute_intersection_points(int& current_node){ + for(typename Edge_to_intersected_facets::iterator it=edge_to_sfacet.begin();it!=edge_to_sfacet.end();++it){ + Halfedge_handle edge=it->first; + Facet_set& fset=it->second; + while (!fset.empty()){ + Facet_handle facet=*fset.begin(); + + Intersection_result res=internal_IOP::do_intersect(edge,facet); + internal_IOP::Intersection_type type=CGAL::cpp0x::get<0>(res); + + //handle degenerate case: one extremity of edge below to facet + std::vector all_edges; + if ( CGAL::cpp0x::get<2>(res) ) + get_incident_edges_to_vertex(edge,std::back_inserter(all_edges)); + else{ + if ( CGAL::cpp0x::get<3>(res) ) + get_incident_edges_to_vertex(edge->opposite(),std::back_inserter(all_edges)); + else + all_edges.push_back(edge); + } + + CGAL_precondition(*all_edges.begin()==edge || *all_edges.begin()==edge->opposite()); +// print_type_debug(type,CGAL::cpp0x::get<2>(res),CGAL::cpp0x::get<3>(res)); + + #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES + check_coplanar_edges(boost::next(all_edges.begin()),all_edges.end(),CGAL::cpp0x::get<1>(res),type); + #endif + + typename std::vector::iterator it_edge=all_edges.begin(); + switch(type){ + case internal_IOP::COPLNR: + #ifndef DO_NOT_HANDLE_COPLANAR_FACETS + assert(!"COPLNR : this point should never be reached!"); + #else + //nothing need to be done, cf. comments at the beginning of the file + #endif + break; + case internal_IOP::EMPTY: + fset.erase(fset.begin()); + CGAL_assertion(!cgal_do_intersect_debug(edge,facet)); + break; + + case internal_IOP::FACET: + { + CGAL_assertion(cgal_do_intersect_debug(edge,facet)); + CGAL_assertion(facet==CGAL::cpp0x::get<1>(res)->face()); + + int node_id=++current_node; + add_new_node(edge,facet,res,nodes); + visitor->new_node_added(node_id,internal_IOP::FACET,edge,facet->halfedge(),CGAL::cpp0x::get<2>(res),CGAL::cpp0x::get<3>(res)); + for (;it_edge!=all_edges.end();++it_edge){ + add_intersection_point_to_facet_and_all_edge_incident_facets(facet,*it_edge,node_id); + //erase facet from the list to test intersection with it_edge + if ( it_edge==all_edges.begin() ) + fset.erase(fset.begin()); + else + { + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + if(it_ets!=edge_to_sfacet.end()) it_ets->second.erase(facet); + } + } + } + break; + + case internal_IOP::EDGE: + { + CGAL_assertion(cgal_do_intersect_debug(edge,facet)); + int node_id=++current_node; + add_new_node(edge,facet,res,nodes); + Halfedge_handle edge_intersected=CGAL::cpp0x::get<1>(res); + visitor->new_node_added(node_id,internal_IOP::EDGE,edge,edge_intersected,CGAL::cpp0x::get<2>(res),CGAL::cpp0x::get<3>(res)); + for (;it_edge!=all_edges.end();++it_edge){ + if ( it_edge!=all_edges.begin() ){ + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + Facet_set* fset_bis = (it_ets!=edge_to_sfacet.end())?&(it_ets->second):NULL; + cip_handle_case_edge(node_id,fset_bis,*it_edge,edge_intersected); + } + else + cip_handle_case_edge(node_id,&fset,*it_edge,edge_intersected); + } + } + break; + + case internal_IOP::VERTEX: + { + CGAL_assertion(cgal_do_intersect_debug(edge,facet)); + int node_id=++current_node; + Halfedge_handle vertex_intersected=CGAL::cpp0x::get<1>(res); + add_new_node(vertex_intersected->vertex()->point()); //we use the original vertex to create the node + //before it was internal_IOP::FACET but do not remember why, probably a bug... + visitor->new_node_added(node_id,internal_IOP::VERTEX,edge,vertex_intersected,CGAL::cpp0x::get<2>(res),CGAL::cpp0x::get<3>(res)); + for (;it_edge!=all_edges.end();++it_edge){ + if ( it_edge!=all_edges.begin() ){ + typename Edge_to_intersected_facets::iterator it_ets=edge_to_sfacet.find(*it_edge); + Facet_set* fset_bis = (it_ets!=edge_to_sfacet.end())?&(it_ets->second):NULL; + cip_handle_case_vertex(node_id,fset_bis,*it_edge,vertex_intersected); + } + else + cip_handle_case_vertex(node_id,&fset,*it_edge,vertex_intersected); + } + } + break; + + } + } + } + CGAL_assertion(nodes.size()==unsigned(current_node+1)); + } + + #ifdef USE_DETECTION_MULTIPLE_DEFINED_EDGES + void remove_duplicated_intersecting_edges() + { + for (typename Coplanar_duplicated_intersection_set::iterator + it=coplanar_duplicated_intersection.begin(); + it!=coplanar_duplicated_intersection.end(); + ++it) + { + typename Facets_to_nodes_map::iterator res=f_to_node.find(*it); + //CGAL_assertion(res!=f_to_node.end()); + //we cannot use a precondition here: in some cases the coplanarity test is not sufficient. + //if on surface1 we have to coplanar triangle incident along edge e. Then in surface2, take a triangle + //with one edge inside one triangle of surface1 such that one vertex in on e. Then the collinearity test + //return true for both faces but only one implies a duplicated edge + if(res!=f_to_node.end()) + f_to_node.erase(res); + } + } + #else + void remove_duplicated_intersecting_edges() + { + std::set< std::pair > already_seen; + std::list to_erase; + for (typename Facets_to_nodes_map::iterator + it=f_to_node.begin(); + it!=f_to_node.end(); + ++it + ) + { + if (it->second.size()==1) continue; + CGAL_precondition(it->second.size()==2); + //it->second is a set so int are already sorted + bool is_new=already_seen.insert( std::make_pair( + *(it->second.begin()), + *boost::next(it->second.begin()) ) + ).second; + + if (!is_new) + to_erase.push_back(it); + } + + for ( typename std::list::iterator + it=to_erase.begin();it!=to_erase.end();++it + ) + f_to_node.erase(*it); + } + #endif + + + struct Graph_node{ + std::set neighbors; + unsigned size; + + Graph_node():size(0){} + + void insert(int i){ + ++size; + CGAL_assertion(neighbors.find(i)==neighbors.end()); + neighbors.insert(i); + } + + void erase(int i){ + CGAL_assertion(neighbors.find(i)!=neighbors.end()); + neighbors.erase(i); + } + void make_terminal() {size=45;} + bool is_terminal()const {return size!=2;} + bool empty() const {return neighbors.empty();} + int top() const {return *neighbors.begin();} + void pop() { + CGAL_assertion(!neighbors.empty()); + neighbors.erase(neighbors.begin()); + } + }; + + + template + void construct_polylines(Node_vector& nodes,Output_iterator out){ + typedef std::map Graph; + Graph graph; + + //counts the number of time each node has been seen + std::size_t nb_nodes=nodes.size(); + std::vector node_mult(nb_nodes,0); + bool isolated_point_seen=false; + for (typename Facets_to_nodes_map::iterator it=f_to_node.begin();it!=f_to_node.end();++it){ + const std::set& segment=it->second; + CGAL_assertion(segment.size()==2 || segment.size()==1); + if (segment.size()==2){ + int i=*segment.begin(); + int j=*boost::next(segment.begin()); + typename Graph::iterator ins_res=graph.insert(std::make_pair(i,Graph_node())).first; + ins_res->second.insert(j); + ins_res=graph.insert(std::make_pair(j,Graph_node())).first; + ins_res->second.insert(i); + ++(node_mult[i]); + ++(node_mult[j]); + } + else{ + CGAL_assertion(segment.size()==1); + isolated_point_seen=true; + } + } + + //add isolated points + if (isolated_point_seen){ + for (unsigned index=0;index(1,nodes[index]); + visitor->start_new_polyline(index,index); + } + } + + //visitor call + visitor->annotate_graph(graph.begin(),graph.end()); + + bool only_cycle=false; + while (!graph.empty()){ + typename Graph::iterator it=graph.begin(); + for (;!only_cycle && it!=graph.end();++it){ + if (it->second.is_terminal()) break; + } + + std::vector polyline; + + if(!only_cycle && it!=graph.end()){ + //this is a polyline + int i=it->first; + int j=it->second.top(); + visitor->start_new_polyline(i,j); + CGAL_assertion(i!=j); + it->second.pop(); + if (it->second.empty()) + graph.erase(it); + polyline.push_back(nodes[i]); + visitor->add_node_to_polyline(i); + while(true){ + it=graph.find(j); + CGAL_assertion(it!=graph.end()); + it->second.erase(i); + i=j; + polyline.push_back(nodes[i]); + visitor->add_node_to_polyline(i); + if (it->second.empty()){ + graph.erase(it); + break; + } + if (it->second.is_terminal()) break; + j=it->second.top(); + it->second.pop(); + if (it->second.empty()) + graph.erase(it); + } + *out++=polyline; + } + else{ + //it remains only cycles + only_cycle=true; + it=graph.begin(); + int i=it->first; + int j=it->second.top(); + visitor->start_new_polyline(i,j); + graph.erase(it); + polyline.push_back(nodes[i]); + visitor->add_node_to_polyline(i); + int first=i; + do{ + it=graph.find(j); + it->second.erase(i); + i=j; + polyline.push_back(nodes[i]); + visitor->add_node_to_polyline(i); + j=it->second.top(); + graph.erase(it); + }while(j!=first); + polyline.push_back(nodes[j]);// we duplicate first point for cycles + visitor->add_node_to_polyline(j); + *out++=polyline; + } + } + } + + int get_other_int(const std::set& s,int i) const { + if (*s.begin()!=i) return *s.begin(); + return *boost::next(s.begin()); + } + + template + bool is_grabbed(const Dispatch_out_it&) const{ + return CGAL::Is_in_tuple::value; + } + + + template + struct Is_dispatch_based_ouput_iterator{ + typedef boost::false_type type; + }; + + template class Dispatch_based_output_it,class V,class O> + struct Is_dispatch_based_ouput_iterator >{ + typedef typename boost::is_base_of< Dispatch_output_iterator, + Dispatch_based_output_it >::type type; + }; + + template + inline void construct_polylines_with_info(Node_vector& nodes,Output_iterator out){ + return construct_polylines_with_info(nodes,out,typename Is_dispatch_based_ouput_iterator::type()); + } + + template + void construct_polylines_with_info(Node_vector& nodes,Output_iterator out,boost::false_type){ + construct_polylines_with_info(nodes, + dispatch_or_drop_output >(out),boost::true_type()); + } + + template class Dispatch_based_output_it,class V,class O> + void construct_polylines_with_info(Node_vector& nodes, + Dispatch_based_output_it out,boost::true_type) + { + typedef typename Facets_to_nodes_map::value_type Edge; + typedef std::list Polyline_info; + + + std::size_t nb_nodes=nodes.size(); + std::vector node_mult(nb_nodes,0); + std::vector terminal_bools(nb_nodes,false); + std::vector< std::set > connections(nb_nodes); + // --counts the number of time each node has been seen + // --associate to each node its incident edges. An edge = a pair of Facet_handle+2 indices of intersection points + for (typename Facets_to_nodes_map::iterator it=f_to_node.begin();it!=f_to_node.end();++it){ + const std::set& segment=it->second; + CGAL_assertion(segment.size()==2 || segment.size()==1); + if (segment.size()==2){ + int i=*segment.begin(); + int j=*boost::next(segment.begin()); + connections[i].insert(&(*it)); + connections[j].insert(&(*it)); + ++(node_mult[i]); + ++(node_mult[j]); + } + } + + //detect terminal nodes and isolated nodes + for (unsigned k=0;k(1,nodes[k]); + if ( is_grabbed >(out)) + *out++=std::vector(); + } + } + + //visitor call + visitor->update_terminal_nodes(terminal_bools); + + //We start from a node N and recursively walk one edge to find other + // node. If this is a cycle we stop when we are back at the node N; + //otherwise we stop at a terminal node and restart a walk from N + //following another edge (if N is not terminal) until finding a terminal + //node. With this we can associate to each edge the pair of facet_handle + //intersecting that define this edge. + unsigned current_node=0; + while (current_node!=nb_nodes){ + if (connections[current_node].empty()){ + ++current_node; + continue; + } + + Edge* edge=*connections[current_node].begin(); + connections[current_node].erase(connections[current_node].begin()); + + Polyline_info polyline_info; + std::list polyline_embedding; + + if ( is_grabbed >(out)) + polyline_info.push_back(edge->first.first); + polyline_embedding.push_back(nodes[current_node]); + + unsigned i=current_node; + unsigned start_node=current_node; + bool reverse=false; + while (true){ + i=get_other_int(edge->second,i); + connections[i].erase(edge); + + if (reverse) polyline_embedding.push_front(nodes[i]); + else polyline_embedding.push_back(nodes[i]); + + if (i==start_node) break; + if (terminal_bools[i]){ + if (reverse || terminal_bools[current_node]) break; + reverse=true; + i=current_node; + if ( connections[i].empty() ) break; + } + + edge=*connections[i].begin(); + connections[i].erase(connections[i].begin()); + + if ( is_grabbed >(out)){ + if (reverse) polyline_info.push_front(edge->first.first); + else polyline_info.push_back(edge->first.first); + } + + } + + *out++=std::vector(polyline_embedding.begin(),polyline_embedding.end()); + if ( is_grabbed >(out)){ + CGAL_assertion(polyline_embedding.size()==polyline_info.size()+1); + *out++=std::vector(polyline_info.begin(),polyline_info.end()); + } + } + } + +//debug functions + + bool cgal_do_intersect_debug(Halfedge_handle eh,Facet_handle fh){ + Triangle t( fh->halfedge()->vertex()->point(), + fh->halfedge()->next()->vertex()->point(), + fh->halfedge()->next()->next()->vertex()->point()); + + Segment s( eh->vertex()->point(), + eh->opposite()->vertex()->point()); + + return CGAL::do_intersect( s, t); + } + + bool cgal_do_intersect_debug(Facet_handle fh1,Facet_handle fh2){ + Triangle t1( fh1->halfedge()->vertex()->point(), + fh1->halfedge()->next()->vertex()->point(), + fh1->halfedge()->next()->next()->vertex()->point()); + Triangle t2( fh2->halfedge()->vertex()->point(), + fh2->halfedge()->next()->vertex()->point(), + fh2->halfedge()->next()->next()->vertex()->point()); + + + return CGAL::do_intersect( t1, t2); + } + + void print_f_to_node_debug(){ + std::cout << "print_f_to_node_debug " << &f_to_node << std::endl; + for (typename Facets_to_nodes_map::iterator it=f_to_node.begin();it!=f_to_node.end();++it){ + std::cout << &(*(it->first.first.first)) << " " << &(*(it->first.first.second)) << " " << it->first.second << " -> {"; + std::copy(it->second.begin(),it->second.end(),std::ostream_iterator(std::cout,",")); + std::cout << "}" <& graph){ + for (typename std::map::const_iterator it=graph.begin();it!=graph.end();++it){ + std::cout << it->first << " -> {"; + std::copy(it->second.neighbors.begin(),it->second.neighbors.end(),std::ostream_iterator(std::cout,",")); + std::cout << "}" <second; + std::cout << &fset << " fset size " << fset.size() << std::endl; + } + } + + + template + OutputIterator main_run(OutputIterator out,bool build_polylines=true){ + // std::cout << edge_to_sfacet.size() << std::endl; + int current_node=-1; + + //print_edge_to_sfacet_debug(); + #ifndef DO_NOT_HANDLE_COPLANAR_FACETS + //first handle coplanar triangles + compute_intersection_of_coplanar_facets(current_node); + visitor->set_number_of_intersection_points_from_coplanar_facets(current_node+1); + #endif + //print_edge_to_sfacet_debug(); + //compute intersection points of segments and triangles. + //build node of the graph + //build connectivity info + compute_intersection_points(current_node); + + if (!build_polylines){ + visitor->finalize(nodes); + return out; + } + //remove duplicated intersecting edges: + // In case two facets are incident along such an edge coplanar in a facet of another polyhedron (and one extremity inside the facet), the intersection + // will be reported twice. We kept track (check_coplanar_edge(s)) of this so that, we can remove one intersecting edge out of the two + //print_f_to_node_debug(); + remove_duplicated_intersecting_edges(); + + //std::cout << "f_to_node "<< f_to_node.size() << std::endl; + //print_f_to_node_debug(); + //collect connectivity infos and create polylines + if ( Node_visitor::do_need_vertex_graph ) + construct_polylines(nodes,out); //using the graph approach (at some point we know all connections between intersection points) + else + construct_polylines_with_info(nodes,out); //direct construction by propagation + + visitor->finalize(nodes); + + return out; + } + +public: + Intersection_of_Polyhedra_3():visitor(new Node_visitor()),is_default_visitor(true){} + Intersection_of_Polyhedra_3(Node_visitor& v):visitor(&v),is_default_visitor(false){} + ~Intersection_of_Polyhedra_3(){if (is_default_visitor) delete visitor;} + //pairwise intersection between all elements in the range + template + OutputIterator + operator()(InputIterator begin, InputIterator end, OutputIterator out) { + for(InputIterator it1=begin;it1!=end;++it1){ + CGAL_precondition( it1->is_pure_triangle() ); + Polyhedron_ref P=*it1; + visitor->new_input_polyhedron(P); + for(InputIterator it2=boost::next(it1);it2!=end;++it2){ + CGAL_precondition( it2->is_pure_triangle() ); + Polyhedron_ref Q=*it2; + filter_intersections(P, Q); + filter_intersections(Q, P); + } + } + + return main_run(out); + } + + //pairwise intersection between all elements in the range + //(pointers version) + template + OutputIterator + operator()(InputIterator begin, InputIterator end, OutputIterator out, int) { + for(InputIterator it1=begin;it1!=end;++it1){ + CGAL_precondition( (*it1)->is_pure_triangle() ); + Polyhedron_ref P=**it1; + visitor->new_input_polyhedron(P); + for(InputIterator it2=boost::next(it1);it2!=end;++it2){ + CGAL_precondition( (*it2)->is_pure_triangle() ); + Polyhedron_ref Q=**it2; + filter_intersections(P, Q); + filter_intersections(Q, P); + } + } + + return main_run(out); + } + + //intersection between P and each element in the range + template + OutputIterator + operator()( Polyhedron_ref P, InputIterator begin, InputIterator end, OutputIterator out) { + CGAL_precondition( P.is_pure_triangle() ); + visitor->new_input_polyhedron(P); + for(InputIterator it=begin;it!=end;++it){ + CGAL_precondition( it->is_pure_triangle() ); + Polyhedron_ref Q=*it; + visitor->new_input_polyhedron(Q); + filter_intersections(P, Q); + filter_intersections(Q, P); + } + return main_run(out); + } + + //intersection between P and each element in the range + //(pointers version) + template + OutputIterator + operator()(Polyhedron_ref P, InputIterator begin, InputIterator end, OutputIterator out, int) { + CGAL_precondition( P.is_pure_triangle() ); + visitor->new_input_polyhedron(P); + for(InputIterator it=begin;it!=end;++it){ + CGAL_precondition( (*it)->is_pure_triangle() ); + Polyhedron_ref Q=**it; + visitor->new_input_polyhedron(Q); + filter_intersections(P, Q); + filter_intersections(Q, P); + } + return main_run(out); + } + + //intersection between P and Q + template + OutputIterator + operator()(Polyhedron_ref P, Polyhedron_ref Q, OutputIterator out) { + visitor->new_input_polyhedron(P); + visitor->new_input_polyhedron(Q); + filter_intersections(P, Q); + filter_intersections(Q, P); + return main_run(out); + } + + //intersection between P and Q, only visitor called not polyline is constructed + void operator()(Polyhedron_ref P, Polyhedron_ref Q) { + visitor->new_input_polyhedron(P); + visitor->new_input_polyhedron(Q); + filter_intersections(P, Q); + filter_intersections(Q, P); + main_run(Emptyset_iterator(),false); + } +}; + +template +OutputIterator +intersection_Polyhedron_3_Polyhedron_3(const Polyhedron& P, const Polyhedron& Q, OutputIterator out) +{ + return Intersection_of_Polyhedra_3()(P,Q,out); +} + +}// namespace CGAL + +#endif //CGAL_INTERSECTION_OF_POLYHEDRA_3_H diff --git a/Polyhedron/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h b/Polyhedron/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h new file mode 100644 index 00000000000..f95ebcbf667 --- /dev/null +++ b/Polyhedron/include/CGAL/intersection_of_Polyhedra_3_refinement_visitor.h @@ -0,0 +1,2149 @@ +// Copyright (c) 2011 GeometryFactory (France). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org); you may redistribute it under +// the terms of the Q Public License version 1.0. +// See the file LICENSE.QPL distributed with CGAL. +// +// 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. +// +// $UR$ +// $Id$ +// +// +// Author(s) : Sebastien Loriot + +#ifndef CGAL_INTERSECTION_OF_POLYHEDRA_3_REFINEMENT_VISITOR_H +#define CGAL_INTERSECTION_OF_POLYHEDRA_3_REFINEMENT_VISITOR_H + +#include +#include +#include +#include +#include + +#include + +#include + +#include + +#include +#include + +// TODOCUMENT +// --We suppose that the two input polyhedra are triangulated orientable surfaces. +// --Any polyhedron defines two bounding volumes: one inside and one outside. +// The convention used is the following: the normal of a triangle always indicates +// the outside of the object. +// --The input polyhedra should not touch at only one point locally. If so, the current +// implementation just ignore it (TAG SL001) +// --Polyhedron type should be list-based or should guarantee no reallocation. We use maps +// on pointer of halfedges,facets and vertices +// --Polyhedral_mesh_domain requires the domain tp be closed: do not provided as input +// an open polyhedral surface and a polyhedron with a connected component free from intersection +//OPTIMIZATIONS +// --cdt: try using intervals? in that case, only points inside the face should be considered +// and points on edge should be handled by hand (simply start using the point opposite to the edge) +// --filtered_order_around_edge: can be done using the original supporting planes +// --in intersection_of_Polyhedra_3: upon call to Triangle_segment_intersection_point::add_new_node, interval and exact nodes are +// inserted into a vector. Since we do not know the final size of vector this lead to reallocation of data. +// --in Triangle_segment_intersection_point, try using EPEC instead of Interval_nt+SC +// --use a sorted pair of indices in edge_to_hedge+simplify the code TAG_SLXX1 +// --in sew_2_marked_darts arrange how darts are passed to avoid comparing to a Point_3 +//TODO: +// --validity of the embedding: points inserted in polyhedron are approximation of the real +// intersection points. It may happen that because of the approximation, the embedding gets +// wrong. To avoid this, for each new triangle created, we should make an orientation test +// with the approximated point to check if this is correct. If not, points must be moved +// within their double interval so that all triangles incident to each of these points are correctly +// oriented. This is probably an expensive test that can be activated only with a template parameter +// of something similar. +namespace CGAL +{ + + namespace internal_IOP + { + template + struct Compare_unik_address{ + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + typedef typename Polyhedron::Halfedge Halfedge; + + bool operator()(Halfedge_handle h1,Halfedge_handle h2) const { + Halfedge* ph1=&(*h1) < &(*h1->opposite()) ? &(*h1) : &(*h1->opposite()); + Halfedge* ph2=&(*h2) < &(*h2->opposite()) ? &(*h2) : &(*h2->opposite()); + return ph1 < ph2; + } + + bool operator()(Halfedge_const_handle h1,Halfedge_const_handle h2) const { + const Halfedge* ph1=&(*h1) < &(*h1->opposite()) ? &(*h1) : &(*h1->opposite()); + const Halfedge* ph2=&(*h2) < &(*h2->opposite()) ? &(*h2) : &(*h2->opposite()); + return ph1 < ph2; + } + }; + + template + struct Compare_address{ + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + typedef typename Polyhedron::Halfedge Halfedge; + + bool operator()(Halfedge_handle h1,Halfedge_handle h2) const { + return &(*h1) < &(*h2); + } + + bool operator()(Halfedge_const_handle h1,Halfedge_const_handle h2) const { + return &(*h1) < &(*h2); + } + }; + + template + class Non_intersection_halfedge{ + typedef std::map< typename Polyhedron::Halfedge_const_handle, + std::pair, + Compare_unik_address + > Intersection_hedges_set; + Intersection_hedges_set intersection_hedges_; + public: + Non_intersection_halfedge(const Intersection_hedges_set& the_set) : intersection_hedges_(the_set){} + + + bool operator()(typename Polyhedron::Halfedge_const_handle h) const + { + return intersection_hedges_.find(h)==intersection_hedges_.end(); + } + }; + + + template + class Triangulate_a_face : public CGAL::Modifier_base { + typedef typename HDS::Halfedge_handle Halfedge_handle; + typedef typename HDS::Vertex_handle Vertex_handle; + typedef typename HDS::Face_handle Face_handle; + typedef typename HDS::Vertex Vertex; + typedef typename HDS::Halfedge Halfedge; + typedef typename HDS::Face Face; + + //data members + Face_handle current_face; + std::map nodes_; + std::map& node_to_polyhedron_vertex_; + std::map,Halfedge_handle>& edge_to_hedge_; + std::vector > edges_to_create_; + std::vector > faces_to_create_; + + typename HDS::Halfedge::Base* + unlock_halfedge(Halfedge_handle h){ + return static_cast(&(*h)); + } + + typename HDS::Face::Base* + unlock_face(Face_handle f){ + return static_cast(&(*f)); + } + + public: + + template + Triangulate_a_face( Face_handle face, + const Node_vector& nodes, + const std::vector& node_ids, + std::map& node_to_polyhedron_vertex, + std::map,Halfedge_handle>& edge_to_hedge, + const Triangulation& triangulation) + :current_face(face),node_to_polyhedron_vertex_(node_to_polyhedron_vertex),edge_to_hedge_(edge_to_hedge) + { + //grab vertices to be inserted to copy them from the vector + for (std::vector::const_iterator it=node_ids.begin();it!=node_ids.end();++it) + { + nodes_.insert(std::make_pair(*it,nodes[*it])); + } + //grab edges that are not on the convex hull (these have already been created) + for (typename Triangulation::Finite_edges_iterator + it=triangulation.finite_edges_begin(); + it!=triangulation.finite_edges_end(); + ++it) + { + typename Triangulation::Vertex_handle v0=it->first->vertex((it->second+1)%3); + typename Triangulation::Vertex_handle v1=it->first->vertex((it->second+2)%3); + //warning in degenerate cases you can insert outsite expected convex hull edges: need exact here. + //an alternative is to test if one the incident faces are infinite (cf assertion below) + if ( edge_to_hedge_.find(std::make_pair(v0->info(),v1->info()))==edge_to_hedge_.end() && + edge_to_hedge_.find(std::make_pair(v1->info(),v0->info()))==edge_to_hedge_.end() ) + { + edges_to_create_.push_back( std::make_pair(v0->info(),v1->info()) ); + } + else + CGAL_assertion( triangulation.is_infinite(it->first->vertex(it->second)) || triangulation.is_infinite( triangulation.mirror_vertex(it->first,it->second)) ); + } + //grab triangles. + for (typename Triangulation::Finite_faces_iterator + it=triangulation.finite_faces_begin(); + it!=triangulation.finite_faces_end(); + ++it) + { + typename Triangulation::Vertex_handle v0=it->vertex(0); + typename Triangulation::Vertex_handle v1=it->vertex(1); + typename Triangulation::Vertex_handle v2=it->vertex(2); + //warning in degenerate case we can have non wanted triangles: need exact here + faces_to_create_.push_back( CGAL::cpp0x::make_tuple( v0->info(),v1->info(),v2->info() ) ); + } + } + + + + void operator()( HDS& hds) { +// std::cerr << "node_to_polyhedron_vertex_"<< std::endl; +// for (typename std::map::iterator it=node_to_polyhedron_vertex_.begin();it!=node_to_polyhedron_vertex_.end();++it) +// std::cerr << it->first << " " << &(*(it->second)) << std::endl; + + //insert the intersection point interior to the face inside the polyhedron and + //save their Polyhedron::vertex_handle + for (typename std::map::iterator it=nodes_.begin();it!=nodes_.end();++it) + { + Vertex_handle v=hds.vertices_push_back(Vertex(it->second)); + CGAL_assertion( node_to_polyhedron_vertex_.find( it->first ) == node_to_polyhedron_vertex_.end()); + node_to_polyhedron_vertex_.insert( std::make_pair(it->first,v) ); +// std::cerr << "vertices " << it->first << " " << &(*v) << std::endl; + } + + //insert the new halfedge and set their incident vertex + for (typename std::vector >::iterator + it=edges_to_create_.begin();it!=edges_to_create_.end();++it) + { + Halfedge_handle he=hds.edges_push_back(Halfedge(),Halfedge()); + + //associate edge to halfedge going from i to j with j as incident vertex + CGAL_assertion(node_to_polyhedron_vertex_.find(it->second)!= node_to_polyhedron_vertex_.end()); + Vertex_handle v=node_to_polyhedron_vertex_.find(it->second)->second; + unlock_halfedge(he)->set_vertex( v ); + v->set_halfedge(he); +// std::cerr << " --in edge " << &(*v) << std::endl; + edge_to_hedge_.insert( std::make_pair(*it,he) ); + v=node_to_polyhedron_vertex_.find(it->first)->second; +// std::cerr << " --in edge " << &(*v) << std::endl; + unlock_halfedge( he->opposite() )->set_vertex( v ); + v->set_halfedge(he->opposite()); + edge_to_hedge_.insert( std::make_pair(std::make_pair(it->second,it->first),he->opposite()) ); +// std::cerr << "edges " << it->first << " " << it->second << std::endl; + } + + std::vector >::iterator it=faces_to_create_.begin(); + + //create the new faces and update adjacencies + while (true) + { + int i=cpp0x::get<0>(*it),j=cpp0x::get<1>(*it),k=cpp0x::get<2>(*it); +// std::cerr << "faces " << i << " " << j << " " << k<< std::endl; + Halfedge_handle current = edge_to_hedge_.find(std::make_pair(i,j))->second; + Halfedge_handle next = edge_to_hedge_.find(std::make_pair(j,k))->second; + Halfedge_handle previous = edge_to_hedge_.find(std::make_pair(k,i))->second; + + + CGAL_assertion (edge_to_hedge_.find(std::make_pair(i,j))!=edge_to_hedge_.end()); + CGAL_assertion (edge_to_hedge_.find(std::make_pair(j,k))!=edge_to_hedge_.end()); + CGAL_assertion (edge_to_hedge_.find(std::make_pair(k,i))!=edge_to_hedge_.end()); + + CGAL_assertion(current->vertex()==node_to_polyhedron_vertex_.find(j)->second); + CGAL_assertion(next->vertex()==node_to_polyhedron_vertex_.find(k)->second); + CGAL_assertion(previous->vertex()==node_to_polyhedron_vertex_.find(i)->second); + + unlock_halfedge(current)->set_next(next); + unlock_halfedge(next)->set_next(previous); + unlock_halfedge(previous)->set_next(current); + + unlock_halfedge(current)->set_prev(previous); + unlock_halfedge(next)->set_prev(current); + unlock_halfedge(previous)->set_prev(next); + + //update face halfedge + unlock_face(current_face)->set_halfedge(current); + + //update face of halfedges + unlock_halfedge(current) ->set_face(current_face); + unlock_halfedge(next) ->set_face(current_face); + unlock_halfedge(previous) ->set_face(current_face); + + if ( ++it!=faces_to_create_.end() ) + current_face=hds.faces_push_back(Face()); + else + break; + } + } + }; + + } //namespace internal_IOP + + +//Considering the plane with normal vector [O_prime,O] and containing O. +//We define the counterclockwise order around O when looking from the side of the plane +//into which the vector [O_prime,O] is pointing. +//We consider the portion of the plane defined by rotating a ray starting at O +//from the planar projection of P1 to the planar projection of P2 in counterclockwise order. +//The predicates indicates whether the planar projection of point Q lies in this portion of the plane. +//Preconditions: +// O_prime,O,P1 are not collinear +// O_prime,O,P2 are not collinear +// O_prime,O,Q are not collinear +// O_prime,O,P1,Q are not coplanar or coplanar_orientation(O,O_prime,P1,Q)==NEGATIVE +// O_prime,O,P2,Q are not coplanar or coplanar_orientation(O,O_prime,P2,Q)==NEGATIVE +template +bool is_in_interior_of_object( + const typename Kernel::Point_3& O_prime,const typename Kernel::Point_3& O, + const typename Kernel::Point_3& P1,const typename Kernel::Point_3& P2, + const typename Kernel::Point_3& Q) +{ + //guarantee to have non-flat triangles + CGAL_precondition( !collinear(O_prime,O,P1) ); + CGAL_precondition( !collinear(O_prime,O,P2) ); + CGAL_precondition( !collinear(O_prime,O,Q) ); + + //no two triangles are coplanar and on the same side of their common edge + CGAL_precondition( !coplanar(O_prime,O,P1,Q) || coplanar_orientation(O,O_prime,P1,Q)==NEGATIVE ); + CGAL_precondition( !coplanar(O_prime,O,P2,Q) || coplanar_orientation(O,O_prime,P2,Q)==NEGATIVE ); + + Sign s0 = sign( determinant(O-O_prime,P1-O,P2-O) ); + + if ( s0==ZERO ){ + //O, O_prime, P1 and P2 are coplanar + Orientation o=orientation(O_prime,O,P1,Q); + CGAL_precondition(o!=COPLANAR); + return o==POSITIVE; + } + + //O, O_prime, P1 and P2 are not coplanar + Sign s1 = sign( determinant(O-O_prime,P1-O,Q -O) ); + Sign s2 = sign( determinant(O-O_prime,Q -O,P2-O) ); + + if (s0 == POSITIVE) // the angle P1,O,P2 is smaller that Pi. + return ( s1 == POSITIVE ) && ( s2 ==POSITIVE ); //true if the angles P1,O,Q and Q,O,P2 are smaller than Pi + else + return ( s1 != NEGATIVE ) || ( s2 != NEGATIVE ); //true if the angle P1,O,Q or the angle Q,O,P2 is smaller than or equal to Pi +} + +//import into the combinatorial map facets in the given range. +//they are supposed to be in the same connected component. +//two volume are created (each facets gives two opposite orientation 2-cell in the map) +template +typename Map::Dart_handle import_from_polyhedron_subset( Map& amap, + Face_iterator faces_begin, + Face_iterator faces_end, + const Non_special_edge_predicate& is_non_special_edge, + Halfedge_to_dart_map_& selected_hedge_to_dart, + int mark_index + ) +{ + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + typedef std::map < Halfedge_const_handle, typename Map::Dart_handle,internal_IOP::Compare_address > Halfedge_to_dart_map; + + Halfedge_to_dart_map hedge_to_dart; + typename Map::Dart_handle first_dart = NULL; + // First traversal to build the darts and link them. + for (Face_iterator it_face = faces_begin; it_face != faces_end; ++it_face) + { + Halfedge_const_handle start=(*it_face)->halfedge(); + + CGAL_precondition(start->next()!=start); + + Halfedge_const_handle current=start; + typename Map::Dart_handle prev = NULL; + typename Map::Dart_handle first_dart_of_face = NULL; + do + { + typename Map::Dart_handle d = amap.create_dart(); + amap.template link_beta<3>(d,amap.create_dart()); //for opposite volume + hedge_to_dart[current] = d; + + if (prev != NULL){ + amap.template link_beta<1>(prev, d); + amap.template link_beta<1>(d->beta(3),prev->beta(3));//for opposite volume + } + else + { + first_dart_of_face = d; + if (first_dart==NULL) first_dart=d; + } + + if ( is_non_special_edge (current) ){ + if ( !current->is_border_edge() ){ + CGAL_assertion(current != current->opposite()); + typename Halfedge_to_dart_map::iterator it = hedge_to_dart.find(current->opposite()); + if (it != hedge_to_dart.end()){ //link the opposites halfedges only when both corresponding darts have been created + amap.template link_beta<2>(d, it->second); + amap.template link_beta<2>(d->beta(3), it->second->beta(3));//for opposite volume + } + } + } + else{ + typename Halfedge_to_dart_map_::iterator it_hedge_map=selected_hedge_to_dart.find(current); + //all marked hedges are not the selected one for its polyline + if ( it_hedge_map!=selected_hedge_to_dart.end() ) it_hedge_map->second=d; + //darts d and d->beta(3) are special edges + amap.mark(d,mark_index); + amap.mark(d->beta(3),mark_index); + } + prev = d; + current=current->next(); + } + while (current != start); + amap.template link_beta<1>(prev, first_dart_of_face); + amap.template link_beta<1>(first_dart_of_face->beta(3),prev->beta(3));//for opposite volume + } + + // Second traversal to update the geometry. + // We run one again through the facets of the HDS. + for (Face_iterator it_face = faces_begin; it_face != faces_end; ++it_face) + { + Halfedge_const_handle start=(*it_face)->halfedge(); + Halfedge_const_handle current=start; + do + { + typename Map::Dart_handle d = hedge_to_dart[current]; // Get the dart associated to the Halfedge + if (d->template attribute<0>() == NULL) + { + amap.template set_attribute<0>(d, + amap.template create_attribute<0>(current->opposite()->vertex()->point())); + } + current=current->next(); + } + while (current != start); + } + + return first_dart; +} + + //turn around the target vertex of dart to find a marked dart +template +boost::optional +next_marked_dart_around_target_vertex( + const Combinatorial_map_3& final_map, + typename Combinatorial_map_3::Dart_handle dart, + int mark_index) +{ + CGAL_precondition(final_map.is_marked(dart,mark_index)); + typename Combinatorial_map_3::Dart_handle next=dart->beta(1); + while ( ! final_map.is_marked(next,mark_index) ){ + if (next->is_free(2) )//we reach a boundary + return boost::optional(); + next=next->beta(2)->beta(1); + } + if (next == dart) //no new dart have been found + return boost::optional(); + CGAL_precondition(&dart->beta(1)->template attribute<0>()->point() == &next->template attribute<0>()->point()); + return boost::optional (next); +} + +//turn around the target vertex of dart to find a marked dart +//with expected_target as target vertex +template +typename Combinatorial_map_3::Dart_handle +get_next_marked_dart_around_target_vertex( + const Combinatorial_map_3& final_map, + typename Combinatorial_map_3::Dart_handle dart, + int mark_index) +{ + CGAL_precondition(final_map.is_marked(dart,mark_index)); + typename Combinatorial_map_3::Dart_handle next=dart->beta(1); + while ( !final_map.is_marked(next,mark_index) ){ + CGAL_assertion( !next->is_free(2) ); + next=next->beta(2)->beta(1); + CGAL_assertion(next != dart); + } + CGAL_precondition(&dart->beta(1)->template attribute<0>()->point() == &next->template attribute<0>()->point()); + return next; +} + +//turn around the source vertex of dart to find a marked dart +//with expected_source as source vertex +template +typename Combinatorial_map_3::Dart_handle +get_next_marked_dart_around_source_vertex( + const Combinatorial_map_3& final_map, + typename Combinatorial_map_3::Dart_handle dart, + int mark_index) +{ + CGAL_precondition(final_map.is_marked(dart,mark_index)); + typename Combinatorial_map_3::Dart_handle next=dart->beta(0); + while ( ! final_map.is_marked(next,mark_index) ){ + CGAL_assertion( !next->is_free(2) ); + next=next->beta(2)->beta(0); + CGAL_assertion(next != dart); + } + CGAL_precondition(&dart->template attribute<0>()->point() == &next->beta(1)->template attribute<0>()->point()); + return next; +} + +//given two marked darts, this function links these two darts with beta<2> +//but in addition it follows the marked darts connected to the same vertex +//(there should be only one) to connect them all together +//( this function is a kind of zipper ;) ) +template +void sew_2_marked_darts( Combinatorial_map_3& final_map, + typename Combinatorial_map_3::Dart_handle dart_1 , + typename Combinatorial_map_3::Dart_handle dart_2 , + int mark_index, + const Node_vector& nodes, + const std::pair& indices, + const std::pair& polyline_info) +{ + typedef boost::optional O_Dart_handle; + + CGAL_precondition( dart_1->is_free(2) ); + CGAL_precondition( dart_2->is_free(2) ); + CGAL_precondition( final_map.is_marked(dart_1,mark_index) ); + CGAL_precondition( final_map.is_marked(dart_2,mark_index) ); + CGAL_precondition( dart_1->template attribute<0>()->point() == dart_2->beta(1)->template attribute<0>()->point() ); + CGAL_precondition( dart_1->beta(1)->template attribute<0>()->point() == dart_2->template attribute<0>()->point() ); + + int src_index = ( ( indices.first < indices.second) == polyline_info.first ) + ? indices.second:indices.first; + + if ( dart_1->template attribute<0>()->point() != nodes[ src_index ] ) std::swap(dart_1,dart_2); + + int nb_segs=polyline_info.second-1,k=1; + + do{ + CGAL_precondition( final_map.template is_sewable<2>(dart_1,dart_2) ); + final_map.template sew<2>(dart_1,dart_2); + + if (k==nb_segs) break; + + dart_1=get_next_marked_dart_around_target_vertex(final_map,dart_1,mark_index); + dart_2=get_next_marked_dart_around_source_vertex(final_map,dart_2,mark_index); + } + while(++k); +} + +//not_top and not_down are two darts from volumes that get merged with an existing +//other one because of a set of identical coplanar triangles. +//top and down is the dart of the volumes "replacing" that of not_top and not down respectively, +//The function is considering all triangles that are bounded by a cycle of marked edges. +//The volume not_top and not_down are part of are those that will disappear at the +//end of the main algorithm. +//( this function is a kind of facet gluer ;) ) +template +void sew_3_marked_darts( Combinatorial_map_3& final_map, + typename Combinatorial_map_3::Dart_handle not_top , + typename Combinatorial_map_3::Dart_handle not_down , + typename Combinatorial_map_3::Dart_handle top , + typename Combinatorial_map_3::Dart_handle down , + int mark_index, + std::set& darts_to_remove) +{ + typedef boost::optional O_Dart_handle; + + if ( not_top->template attribute<3>()->info().is_empty ){ + CGAL_assertion(not_down->template attribute<3>()->info().is_empty); + return; + } + + CGAL_assertion(!not_down->template attribute<3>()->info().is_empty); + + //merge attribute of the two volumes: + internal_IOP::Volume_on_merge merge_attributes; + merge_attributes(*top->template attribute<3>(),*not_top->template attribute<3>()); + merge_attributes(*down->template attribute<3>(),*not_down->template attribute<3>()); + + //set volume attributes as empty to avoid double sew_3 of the same topological disk of triangles + not_top->template attribute<3>()->info().is_empty=true; + not_down->template attribute<3>()->info().is_empty=true; + + CGAL_precondition( final_map.is_marked(not_top,mark_index) && final_map.is_marked(top,mark_index) ); + CGAL_precondition( final_map.is_marked(not_down,mark_index) && final_map.is_marked(down,mark_index) ); + CGAL_precondition( not_top->template attribute<0>()->point() == not_down->beta(1)->template attribute<0>()->point() ); + CGAL_precondition( not_top->beta(1)->template attribute<0>()->point() == not_down->template attribute<0>()->point() ); + CGAL_precondition( not_top->template attribute<0>()->point() == top->template attribute<0>()->point() ); + CGAL_precondition( not_down->template attribute<0>()->point() == down->template attribute<0>()->point() ); + + CGAL_assertion( top->beta(3)==down ); + + //set to be removed the darts of the two no longer used volumes + typename Combinatorial_map_3::Dart_handle start=not_top; + do + { + CGAL_assertion(!not_top->is_free(3)); + darts_to_remove.insert(not_top); darts_to_remove.insert(not_top->beta(1)); darts_to_remove.insert(not_top->beta(1)->beta(1)); + darts_to_remove.insert(not_top->beta(3)); darts_to_remove.insert(not_top->beta(3)->beta(1)); darts_to_remove.insert(not_top->beta(3)->beta(1)->beta(1)); + O_Dart_handle current_1=next_marked_dart_around_target_vertex(final_map,not_top,mark_index); + CGAL_precondition(current_1); + not_top=*current_1; + } + while(not_top!=start); +} + +template +struct Halfedge_marker{ + template + static void mark(Halfedge_handle){} +}; + +template<> +struct Halfedge_marker{ + template + static void mark(Halfedge_handle h){h->set_mark();} +}; + +template +class Node_visitor_refine_polyhedra{ +//typedefs + typedef typename Polyhedron::Halfedge_handle Halfedge_handle; + typedef typename Polyhedron::Halfedge_const_handle Halfedge_const_handle; + typedef typename Polyhedron::Face_handle Face_handle; + typedef typename Polyhedron::Halfedge Halfedge; + typedef typename Polyhedron::Vertex_handle Vertex_handle; + typedef internal_IOP::Compare_handles Cmp_handle; //This ensures uniqueness of edges when comparing halfedges + typedef internal_IOP::Compare_unik_address Cmp_unik_ad; //This ensures uniqueness of edges when comparing halfedges + + //constrained triangulation used for triangulation interior of faces + #ifdef DO_NO_USE_EXACT_CDT + typedef CGAL::Triangulation_vertex_base_with_info_2 Vbi; + typedef CGAL::Constrained_triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS_2; + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; //DO WE NEED DELAUNAY???? + #else + typedef CGAL::Exact_predicates_exact_constructions_kernel Exact_kernel; + typedef CGAL::Triangulation_vertex_base_with_info_2 Vbi; + typedef CGAL::Constrained_triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS_2; + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; //DO WE NEED DELAUNAY???? + #endif + + typedef std::map Hedge_to_polyhedron_map; + + typedef std::vector Node_ids; + typedef std::set Node_id_set; //avoid having duplicated node on edge of coplanar triangles + typedef std::map< Face_handle,Node_ids,Cmp_handle > In_face_map; + typedef std::map< Halfedge_handle,Node_id_set,Cmp_unik_ad > In_halfedge_map; + //to keep the correspondance between node_id and vertex_handle in each polyhedron + typedef std::map Node_to_polyhedron_vertex_map; + typedef std::map Poly_to_map_node; + //to maintain an polyhedron halfedge on each polyline + pair + //with first = "is the key (pair) was reversed?" and second is the number of edges +1 in the polyline + typedef std::map< std::pair, std::pair< std::map,std::pair > > An_edge_per_polyline_map; + //to handle coplanar halfedge of polyhedra that are full in the intersection + typedef std::map< int,Halfedge_handle > Node_to_target_of_hedge_map; + typedef std::map< Polyhedron*,Node_to_target_of_hedge_map> Poly_to_vertices_on_intersection_map; + + //Combinatorial map typedefs + typedef internal_IOP::Item_with_points_and_volume_info Items; + typedef CGAL::Combinatorial_map<3,Items> Combinatorial_map_3_; + typedef typename Combinatorial_map_3_::Dart_handle Dart_handle; + +//data members + Hedge_to_polyhedron_map hedge_to_polyhedron; + In_face_map in_face; + In_halfedge_map in_hedge; + std::map< int,std::set > graph_of_constraints; + std::map< int,std::set > coplanar_constraints; + An_edge_per_polyline_map an_edge_per_polyline; + typename An_edge_per_polyline_map::iterator last_polyline; + Poly_to_vertices_on_intersection_map poly_to_vertices_on_inter; + Poly_to_map_node polyhedron_to_map_node_to_polyhedron_vertex; + std::set non_manifold_nodes; //contain nodes that are original vertices of input polyhedron and that neighborhood is not a topological disk + std::map nodes_that_are_original_vertices;//to keep the correspondance between original polyhedron vertices that are also nodes + + Combinatorial_map_3_* final_map_ptr; + Combinatorial_map_3_& final_map() {return *final_map_ptr;} + bool final_map_comes_from_outside; + // new_hedge hedge + // -----------> -----------> + // v + // <----------- <----------- + // new_opposite opposite + // + Vertex_handle split_edge( Halfedge_handle hedge, + const typename Kernel::Point_3& point, + Polyhedron& P) + { + internal_IOP::Split_halfedge_at_point delegated(hedge,point); + P.delegate( delegated ); + CGAL_assertion(P.is_valid()); + + Vertex_handle v=boost::prior(P.vertices_end()); + CGAL_assertion(v->point()==point); + return v; + } + + //sort node ids so that we can split the hedge + //consecutively + template + void sort_vertices_along_hedge(std::vector& node_ids,Halfedge_handle hedge,const Node_vector& nodes) + { + std::sort(node_ids.begin(), + node_ids.end(), + internal_IOP::Order_along_a_halfedge(hedge,nodes) + ); + } + + //insert intersection as constrained edges in a CDT triangulation + template + void insert_constrained_edges_coplanar_case(int node_id, + CDT& triangulation, + std::map& id_to_CDT_vh) + { + if (node_id < number_coplanar_vertices){ + //XSL_TAG_CPL_VERT + //Insert constrained edges from coplanar facets that have been retriangulated. This ensure that triangulations are compatible + std::map< int,std::set >::iterator it_neighbors=coplanar_constraints.find(node_id); + if (it_neighbors!=coplanar_constraints.end()) + { + typename CDT::Vertex_handle vh=id_to_CDT_vh.find(node_id)->second; + for (std::set::iterator it_n=it_neighbors->second.begin();it_n!=it_neighbors->second.end();++it_n){ + typename std::map::iterator it_vh=id_to_CDT_vh.find(*it_n); + // this condition ensures to consider only graph edges that are in the same triangle (not in a neighbor one when involving node on a triangle edge) + // here we can't make the difference between a point on the interior or the boundary, so points_on_triangle is not used. + if ( it_vh!=id_to_CDT_vh.end() ){ + triangulation.insert_constraint(vh,id_to_CDT_vh.find(*it_n)->second); + } + } + } + } + } + //insert intersection as constrained edges in a CDT triangulation + template + void insert_constrained_edges(Node_ids& node_ids, //index of vertices we are interested in + CDT& triangulation, + std::map& id_to_CDT_vh, + Constrained_edges_map& constrained_edges, //list of pair of int to indicate edges that are constrained + bool points_on_triangle=false) + { + for (Node_ids::iterator it_node_id=node_ids.begin();it_node_id!=node_ids.end();++it_node_id){ + std::map< int,std::set >::iterator it_neighbors=graph_of_constraints.find(*it_node_id); + if (it_neighbors!=graph_of_constraints.end()) + { + typename CDT::Vertex_handle vh=id_to_CDT_vh.find(*it_node_id)->second; + for (std::set::iterator it_n=it_neighbors->second.begin();it_n!=it_neighbors->second.end();++it_n){ + typename std::map::iterator it_vh=id_to_CDT_vh.find(*it_n); + // this condition ensures to consider only graph edges that are in the same triangle (not in a neighbor one when involving node on a triangle edge) + if ( !points_on_triangle || it_vh!=id_to_CDT_vh.end() ){ + CGAL_assertion(it_vh!=id_to_CDT_vh.end()); + triangulation.insert_constraint(vh,id_to_CDT_vh.find(*it_n)->second); + constrained_edges.push_back(std::make_pair(*it_node_id,*it_n)); + } + } + } + #ifndef NDEBUG + else + { + std::cout << "X0: Found an isolated point" << std::endl; + } + #endif + + insert_constrained_edges_coplanar_case(*it_node_id,triangulation,id_to_CDT_vh); + } + } + + std::pair make_sorted_pair(int i,int j) const {return i indices,typename Polyhedron::Halfedge_handle hedge) + { + std::pair sorted_pair=make_sorted_pair(indices.first,indices.second); + typename An_edge_per_polyline_map::iterator it=an_edge_per_polyline.find(sorted_pair); + if (it!=an_edge_per_polyline.end()){ + it->second.first.insert(std::make_pair( P,sorted_pair.first==indices.first?hedge:hedge->opposite() )); + } + } + + int node_index_of_incident_vertex(Halfedge_const_handle h, + const std::map,Cmp_unik_ad >& border_halfedges) + { + //WARNING this may be expensive + typedef std::map,Cmp_unik_ad > Border_halfedges_map; + + Halfedge_const_handle start=h; + Halfedge_const_handle curr=start; + do{ + typename Border_halfedges_map::const_iterator it_border=border_halfedges.find( curr ); + if (it_border!=border_halfedges.end()) + return it_border->first==curr?it_border->second.second:it_border->second.first; + curr=curr->next()->opposite(); + }while(curr!=start); + + return -1; + } + + template + bool filtered_order_around_edge(int O_prime_index, + int O_index, + int P1_index, + int P2_index, + int Q_index, + Vertex_handle P1, + Vertex_handle P2, + Vertex_handle Q, + const Node_vector& nodes) + { + try{ + return is_in_interior_of_object( + nodes.interval_node(O_prime_index), + nodes.interval_node(O_index), + P1_index == -1 ? nodes.to_interval(P1->point()): nodes.interval_node(P1_index), + P2_index == -1 ? nodes.to_interval(P2->point()): nodes.interval_node(P2_index), + Q_index == -1 ? nodes.to_interval(Q->point()) : nodes.interval_node(Q_index ) + ); + } + catch(Uncertain_conversion_exception&){ + return is_in_interior_of_object( + nodes.exact_node(O_prime_index), + nodes.exact_node(O_index), + P1_index == -1 ? nodes.to_exact(P1->point()): nodes.exact_node(P1_index), + P2_index == -1 ? nodes.to_exact(P2->point()): nodes.exact_node(P2_index), + Q_index == -1 ? nodes.to_exact(Q->point()) : nodes.exact_node(Q_index ) + ); + } + } + +//keep track of the fact that a polyhedron original vertex is a node +void all_incident_faces_got_a_node_as_vertex(Halfedge_handle incident_to_vertex_edge,int node_id) +{ + nodes_that_are_original_vertices.insert(std::make_pair(incident_to_vertex_edge->vertex(),node_id)); +} + +//if an original polyhedron vertex is also a node, do no use a fake id +void set_triangle_boundary_indices( + Vertex_handle* triangle_boundary, + int* triangle_boundary_indices) +{ + triangle_boundary_indices[0]=-1; + triangle_boundary_indices[1]=-2; + triangle_boundary_indices[2]=-3; + + for (int k=0;k<3;++k){ + typename std::map::iterator it=nodes_that_are_original_vertices.find(triangle_boundary[k]); + if (it!=nodes_that_are_original_vertices.end()) + triangle_boundary_indices[k]=it->second; + } +} + +//======================================================================// +//functions internally used to glue piece of the final combinatorial map// +//======================================================================// + + + //-----first polyhedron +template +inline Dart_handle get_associated_dart(Halfedge_handle hedge,Halfedge_to_dart_map& selected_hedge_to_dart){ + typename Halfedge_to_dart_map::iterator it_saved_dart=selected_hedge_to_dart.find(hedge); + CGAL_assertion(it_saved_dart!=selected_hedge_to_dart.end()); + return it_saved_dart->second; +} + +//first_hedge defines four volumes, second_hedge only two +//first_poly is not needed as inside/outside volume is update during the merge +//of the sew. Only second_poly is needed +template +void sew_2_three_volumes_case( Halfedge_handle first_hedge, Halfedge_handle second_hedge, + const std::pair& indices, + const Node_vector& nodes, + Border_halfedges_map& border_halfedges, + Halfedge_to_dart_map& selected_hedge_to_dart, + Polyhedron* /*first_poly*/, Polyhedron* second_poly, + int mark_index, + std::set& darts_to_remove, + const std::pair& polyline_info) +{ + bool took_opposite=second_hedge->is_border(); + if (took_opposite) second_hedge=second_hedge->opposite(); + + Vertex_handle P1=first_hedge->opposite()->next()->vertex(); + Vertex_handle P2=first_hedge->next()->vertex(); + // when looking from the side of indices.second, the interior of the first polyhedron is described + // by turning counterclockwise from P1 to P2 + + Vertex_handle Q = second_hedge->next()->vertex(); + + //check if the third point of each triangular face is an original point (stay -1) + //or a intersection point (in that case we need the index of the corresponding node to + //have the exact value of the point) + int index_p1=node_index_of_incident_vertex(first_hedge->opposite()->next(),border_halfedges); + int index_p2=node_index_of_incident_vertex(first_hedge->next(),border_halfedges); + int index_q =node_index_of_incident_vertex(second_hedge->next(),border_halfedges); + + //Recover the dart that will be the start point of the different sewing + // dof_X_outside = dart of face of , meaning the triangle containing the + // point X and part of the volume outside of the corresponding polyhedron + //-----first polyhedron + Dart_handle dof_P1_outside = get_associated_dart(first_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_P2_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + if (index_p1!=-1 && index_p1==index_q){ + Dart_handle top=dof_P1_outside->beta(3), not_top=took_opposite?dof_Q_outside->beta(3):dof_Q_outside; + Dart_handle down=dof_P1_outside, not_down=took_opposite?dof_Q_outside:dof_Q_outside->beta(3); + + if ( top->template attribute<3>()->info().is_empty ) std::swap(not_top,top); + if ( down->template attribute<3>()->info().is_empty ) std::swap(not_down,down); + CGAL_assertion( !top->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down->template attribute<3>()->info().is_empty ); + + sew_2_marked_darts( final_map(),top , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1P2 or QP2 + sew_2_marked_darts( final_map(),dof_P2_outside , down ,mark_index, nodes, indices, polyline_info); //P2Q or P2P1 + sew_3_marked_darts( final_map(),not_top,not_down,top,down,mark_index,darts_to_remove); + + return; + } + + if (index_p2!=-1 && index_p2==index_q){ + Dart_handle top=dof_P2_outside->beta(3), not_top=took_opposite?dof_Q_outside:dof_Q_outside->beta(3); + Dart_handle down=dof_P2_outside, not_down=took_opposite?dof_Q_outside->beta(3):dof_Q_outside; + + if ( top->template attribute<3>()->info().is_empty ) std::swap(not_top,top); + if ( down->template attribute<3>()->info().is_empty ) std::swap(not_down,down); + CGAL_assertion( !top->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down->template attribute<3>()->info().is_empty ); + + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , top ,mark_index, nodes, indices, polyline_info); //P1Q or P1P2 + sew_2_marked_darts( final_map(),down , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //QP1 or P2P1 + sew_3_marked_darts( final_map(),not_top,not_down,top,down,mark_index,darts_to_remove); + + return; + } + + bool Q_is_between_P1P2 = filtered_order_around_edge(indices.first,indices.second,index_p1,index_p2,index_q,P1,P2,Q,nodes); + + + if (Q_is_between_P1P2) + { + // poly_first - poly_second = took_opposite?P1Q:QP2 + // poly_second - poly_first = {0} + // poly_first \cap poly_second = took_opposite?QP2:P1Q + // opposite( poly_first U poly_second ) = P2P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , took_opposite?dof_Q_outside:dof_Q_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1Q + sew_2_marked_darts( final_map(),took_opposite?dof_Q_outside->beta(3):dof_Q_outside , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //QP2 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //P2P1 + dof_P1_outside->template attribute<3>()->info().outside.insert(second_poly); //update P2P1 outside poly + } + else + { + // poly_first - poly_second = P1P2 + // poly_second - poly_first = took_opposite?QP1:P2Q + // poly_first \cap poly_second = {0} + // opposite( poly_first U poly_second ) = took_opposite?P2Q:QP1 + sew_2_marked_darts( final_map(),dof_P2_outside , took_opposite?dof_Q_outside:dof_Q_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P2Q + sew_2_marked_darts( final_map(),took_opposite?dof_Q_outside->beta(3):dof_Q_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //QP1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1P2 + dof_P1_outside->beta(3)->template attribute<3>()->info().outside.insert(second_poly); //update P1P2 outside poly + } +} + +//first_hedge defines two volumes, second_hedge only two +template +void sew_2_two_volumes_case( Halfedge_handle first_hedge, Halfedge_handle second_hedge, + Border_halfedges_map& border_halfedges, + Halfedge_to_dart_map& selected_hedge_to_dart, + int mark_index, + std::set& darts_to_remove, + const Node_vector& nodes, + const std::pair& indices, + const std::pair& polyline_info) +{ + bool first_took_opposite=first_hedge->is_border(); + if (first_took_opposite) first_hedge=first_hedge->opposite(); + bool second_took_opposite=second_hedge->is_border(); + if (second_took_opposite) second_hedge=second_hedge->opposite(); + + //-----first polyhedron + Dart_handle dof_P_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + + + + int index_p =node_index_of_incident_vertex(first_hedge->next(),border_halfedges); + int index_q =node_index_of_incident_vertex(second_hedge->next(),border_halfedges); + + if (index_p!=-1 && index_q!=-1 && index_p==index_q){ + Dart_handle top=dof_P_outside, not_top=dof_Q_outside->beta(3); + Dart_handle down=dof_P_outside->beta(3), not_down=dof_Q_outside; + + if (first_took_opposite==second_took_opposite) + { + top=dof_P_outside->beta(3); not_top=dof_Q_outside->beta(3); + down=dof_P_outside; not_down=dof_Q_outside; + } + + if ( top->template attribute<3>()->info().is_empty ) std::swap(not_top,top); + if ( down->template attribute<3>()->info().is_empty ) std::swap(not_down,down); + CGAL_assertion( !top->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down->template attribute<3>()->info().is_empty ); + + sew_3_marked_darts( final_map(),not_top,not_down,top,down,mark_index,darts_to_remove); + + return; + } + + + + //since the edge is shared, the inside of each polyhedron must be on opposite orientation halfedges + if (first_took_opposite==second_took_opposite) + { + //sew out with in + sew_2_marked_darts( final_map(),dof_P_outside->beta(3) , dof_Q_outside ,mark_index, nodes, indices, polyline_info); //PQ + sew_2_marked_darts( final_map(),dof_Q_outside->beta(3) , dof_P_outside ,mark_index, nodes, indices, polyline_info); //QP + } + else + { + //sew in with in + sew_2_marked_darts( final_map(),dof_P_outside , dof_Q_outside ,mark_index, nodes, indices, polyline_info); //PQ + sew_2_marked_darts( final_map(),dof_Q_outside->beta(3), dof_P_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //QP + } +} + +//4 volume case with 2 identical volume +//Q2 is supposed to be identical to P2 +template +void sew_2_four_volumes_case_1( Halfedge_handle first_hedge, Halfedge_handle second_hedge, + const std::pair& indices, + const Node_vector& nodes, + int index_p1, int index_p2, int index_q1, + Halfedge_to_dart_map& selected_hedge_to_dart, + int mark_index, + std::set& darts_to_remove, + const std::pair& polyline_info, + bool swap_in_out_Q=false) +{ + Vertex_handle P1=first_hedge->opposite()->next()->vertex(); + Vertex_handle P2=first_hedge->next()->vertex(); + // when looking from the side of indices.second, the interior of the first polyhedron is described + // by turning counterclockwise from P1 to P2 + Vertex_handle Q1=second_hedge->opposite()->next()->vertex(); + // Vertex_handle Q2=second_hedge->next()->vertex(); + bool Q1_is_between_P1P2 = filtered_order_around_edge(indices.first,indices.second,index_p1,index_p2,index_q1,P1,P2,Q1,nodes); + + + //Recover the dart that will be the start point of the different sewing + // dof_X_outside = dart of face of , meaning the triangle containing the + // point X and part of the volume outside of the corresponding polyhedron + //-----first polyhedron + Dart_handle dof_P1_outside = get_associated_dart(first_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_P2_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q1_outside = get_associated_dart(second_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_Q2_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + if( swap_in_out_Q ){ + dof_Q1_outside=dof_Q1_outside->beta(3); + dof_Q2_outside=dof_Q2_outside->beta(3); + } + + if (Q1_is_between_P1P2){ + Dart_handle top=dof_Q2_outside->beta(3), not_top=dof_P2_outside->beta(3); + Dart_handle down=dof_Q2_outside, not_down=dof_P2_outside; + if ( top->template attribute<3>()->info().is_empty ) std::swap(not_top,top); + if ( down->template attribute<3>()->info().is_empty ) std::swap(not_down,down); + CGAL_assertion( !top->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down->template attribute<3>()->info().is_empty ); + + // poly_first - poly_second = P1Q1 + // poly_second - poly_first = {0} + // poly_first \cap poly_second = Q1P2 or Q1Q2 + // opposite( poly_first U poly_second ) = Q2P1 or P2P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //P1Q1 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , top ,mark_index, nodes, indices, polyline_info); //Q1P2 or Q1Q2 + sew_2_marked_darts( final_map(),down , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q2P1 or P2P1 + sew_3_marked_darts( final_map(),not_top,not_down,top,down,mark_index,darts_to_remove); + + } + else{ + Dart_handle top=dof_Q2_outside->beta(3), not_top=dof_P2_outside->beta(3); + Dart_handle down=dof_Q2_outside, not_down=dof_P2_outside; + if ( top->template attribute<3>()->info().is_empty ) std::swap(not_top,top); + if ( down->template attribute<3>()->info().is_empty ) std::swap(not_down,down); + CGAL_assertion( !top->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down->template attribute<3>()->info().is_empty ); + + // poly_first - poly_second = {0} + // poly_second - poly_first = Q1P1 + // poly_first \cap poly_second = P1P2 or P1Q2 + // opposite( poly_first U poly_second ) = Q2Q1 or P2Q1 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q1P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , top ,mark_index, nodes, indices, polyline_info); //P1P2 or P1Q2 + sew_2_marked_darts( final_map(),down , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //Q2Q1 or P2Q1 + sew_3_marked_darts( final_map(),not_top,not_down,top,down,mark_index,darts_to_remove); + } +} + +template +bool coplanar_triangles_case_handled(Halfedge_handle first_hedge,Halfedge_handle second_hedge, + const std::pair& indices, + const Node_vector& nodes, + int index_p1, int index_p2, int index_q1,int index_q2, + Halfedge_to_dart_map& selected_hedge_to_dart, + int mark_index, + std::set& darts_to_remove, + const std::pair& polyline_info) +{ + if( index_p1!=-1 ){ + if (index_p1==index_q1){ + if(index_p2!=-1){ + CGAL_assertion(index_p2!=index_q1); + if(index_p2==index_q2){ + //-----first polyhedron + Dart_handle dof_P1_outside = get_associated_dart(first_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_P2_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q1_outside = get_associated_dart(second_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_Q2_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + Dart_handle top_1=dof_P1_outside->beta(3), not_top_1=dof_Q1_outside->beta(3); + Dart_handle top_2=dof_P2_outside->beta(3), not_top_2=dof_Q2_outside->beta(3); + Dart_handle down_1=dof_P1_outside, not_down_1=dof_Q1_outside; + Dart_handle down_2=dof_P2_outside, not_down_2=dof_Q2_outside; + if ( top_1->template attribute<3>()->info().is_empty ) std::swap(top_1,not_top_1); + if ( top_2->template attribute<3>()->info().is_empty ) std::swap(top_2,not_top_2); + if ( down_1->template attribute<3>()->info().is_empty ) std::swap(down_1,not_down_1); + if ( down_2->template attribute<3>()->info().is_empty ) std::swap(down_2,not_down_2); + CGAL_assertion( !top_1->template attribute<3>()->info().is_empty ); + CGAL_assertion( !top_2->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down_1->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down_2->template attribute<3>()->info().is_empty ); + + // poly_first - poly_second = {0} + // poly_second - poly_first = {0} + // poly_first \cap poly_second = P1P2 or Q1Q2 or P1Q1 or Q1P2 + // opposite( poly_first U poly_second ) = P2P1 or Q2Q1 or P2Q1 or Q2P1 + sew_2_marked_darts( final_map(),top_1 , top_2 ,mark_index, nodes, indices, polyline_info); //P1P2 or Q1Q2 or P1Q1 or Q1P2 + sew_2_marked_darts( final_map(),down_2 , down_1 ,mark_index, nodes, indices, polyline_info); //P2P1 or Q2Q1 or P2Q1 or Q2P1 + sew_3_marked_darts( final_map(),not_top_1, not_down_1 ,top_1,down_1,mark_index,darts_to_remove); + sew_3_marked_darts( final_map(),not_top_2, not_down_2 ,top_2,down_2,mark_index,darts_to_remove); + return true; + } + } + sew_2_four_volumes_case_1(first_hedge->opposite(),second_hedge->opposite(),std::make_pair(indices.second,indices.first),nodes,index_p2,index_p1,index_q2,selected_hedge_to_dart,mark_index,darts_to_remove,polyline_info); + return true; + } + if (index_p1==index_q2){ + if(index_p2!=-1){ + CGAL_assertion(index_p2!=index_q2); + if(index_p2==index_q1){ + //-----first polyhedron + Dart_handle dof_P1_outside = get_associated_dart(first_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_P2_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q1_outside = get_associated_dart(second_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_Q2_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + Dart_handle top_1=dof_P1_outside->beta(3), not_top_1=dof_Q2_outside; + Dart_handle top_2=dof_P2_outside->beta(3), not_top_2=dof_Q1_outside; + Dart_handle down_1=dof_P1_outside, not_down_1=dof_Q2_outside->beta(3); + Dart_handle down_2=dof_P2_outside, not_down_2=dof_Q1_outside->beta(3); + if ( top_1->template attribute<3>()->info().is_empty ) std::swap(top_1,not_top_1); + if ( top_2->template attribute<3>()->info().is_empty ) std::swap(top_2,not_top_2); + if ( down_1->template attribute<3>()->info().is_empty ) std::swap(down_1,not_down_1); + if ( down_2->template attribute<3>()->info().is_empty ) std::swap(down_2,not_down_2); + CGAL_assertion( !top_1->template attribute<3>()->info().is_empty ); + CGAL_assertion( !top_2->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down_1->template attribute<3>()->info().is_empty ); + CGAL_assertion( !down_2->template attribute<3>()->info().is_empty ); + + // poly_first - poly_second = P1P2 or P1Q1 or Q2P2 or Q2Q1 + // poly_second - poly_first = Q1Q2 or Q1P1 or P2P1 or P2Q2 + // poly_first \cap poly_second = {0} + // opposite( poly_first U poly_second ) = all space + sew_2_marked_darts( final_map(),top_1 , top_2 ,mark_index, nodes, indices, polyline_info); //P1P2 or Q1Q2 or P1Q1 or Q1P2 + sew_2_marked_darts( final_map(),down_2 , down_1 ,mark_index, nodes, indices, polyline_info); //P2P1 or Q2Q1 or P2Q1 or Q2P1 + sew_3_marked_darts( final_map(),not_top_1, not_down_1 ,top_1,down_1,mark_index,darts_to_remove); + sew_3_marked_darts( final_map(),not_top_2, not_down_2 ,top_2,down_2,mark_index,darts_to_remove); + return true; + } + } + sew_2_four_volumes_case_1(first_hedge->opposite(),second_hedge,std::make_pair(indices.second,indices.first),nodes,index_p2,index_p1,index_q1,selected_hedge_to_dart,mark_index,darts_to_remove,polyline_info,true); + return true; + } + } + + if(index_p2!=-1){ + if (index_p2==index_q1){ + sew_2_four_volumes_case_1(first_hedge,second_hedge->opposite(),indices,nodes,index_p1,index_p2,index_q2,selected_hedge_to_dart,mark_index,darts_to_remove,polyline_info,true); + return true; + } + if(index_p2==index_q2){ + sew_2_four_volumes_case_1(first_hedge,second_hedge,indices,nodes,index_p1,index_p2,index_q1,selected_hedge_to_dart,mark_index,darts_to_remove,polyline_info); + return true; + } + } + + + return false; +} + +//===// +//end// +//===// + bool do_not_build_cmap; //set to true in the case only the corefinement must be done + int number_coplanar_vertices; //number of intersection points between coplanar facets, see fixes XSL_TAG_CPL_VERT +public: + Node_visitor_refine_polyhedra (Combinatorial_map_3_* ptr=NULL,bool do_not_build_cmap_=false):do_not_build_cmap(do_not_build_cmap_) + { + if (ptr!=NULL){ + final_map_comes_from_outside=true; + final_map_ptr=ptr; + } + else{ + final_map_comes_from_outside=false; + final_map_ptr=new Combinatorial_map_3_(); + } + } + + ~Node_visitor_refine_polyhedra(){ + if(!final_map_comes_from_outside) delete final_map_ptr; + } + + + typedef internal_IOP::Predicates_on_constructions Node_storage_type; + typedef Tag_false Is_polyhedron_const; + static const bool do_need_vertex_graph = true; //because we need to know which edges are constrained + + typedef Combinatorial_map_3_ Combinatorial_map_3; + typedef internal_IOP::Volume_info Volume_info; + + const Combinatorial_map_3& combinatorial_map() const {return *final_map_ptr;} + + void set_number_of_intersection_points_from_coplanar_facets(int n){ + number_coplanar_vertices=n; + } + + void check_node_on_non_manifold_vertex(int node_id,Halfedge_handle hedge){ + //we turn around the hedge and check no halfedge is a border halfedge + Halfedge_handle curr=hedge; + do{ + if ( curr->is_border_edge() ){ + non_manifold_nodes.insert(node_id); + return; + } + curr=curr->next()->opposite(); + } + while(curr!=hedge); + } + + void check_node_on_non_manifold_edge(int node_id,Halfedge_handle hedge){ + if ( hedge->is_border_edge() ) non_manifold_nodes.insert(node_id); + } + + void new_node_added(int node_id, + internal_IOP::Intersection_type type, + Halfedge_handle principal_edge, + Halfedge_handle additional_edge, + bool is_vertex_coplanar, + bool is_vertex_opposite_coplanar) + { + switch(type) + { + case internal_IOP::FACET: //Facet intersected by an edge + { + typename In_face_map::iterator it_fmap=in_face.insert(std::make_pair(additional_edge->face(), Node_ids())).first; + it_fmap->second.push_back(node_id); + } + break; + case internal_IOP::EDGE: //Edge intersected by an edge + { + typename In_halfedge_map::iterator it_hedge_map=in_hedge.insert(std::make_pair(additional_edge,Node_id_set())).first; + it_hedge_map->second.insert(node_id); + check_node_on_non_manifold_edge(node_id,additional_edge); + } + break; + case internal_IOP::VERTEX: + { + //grab original vertex that is on commom intersection + typename Hedge_to_polyhedron_map::iterator it=hedge_to_polyhedron.find(additional_edge->facet()->halfedge()); + CGAL_assertion(it!=hedge_to_polyhedron.end()); + poly_to_vertices_on_inter[it->second].insert(std::make_pair(node_id,additional_edge)); + polyhedron_to_map_node_to_polyhedron_vertex[it->second].insert(std::make_pair(node_id,additional_edge->vertex())); + all_incident_faces_got_a_node_as_vertex(additional_edge,node_id); + check_node_on_non_manifold_vertex(node_id,additional_edge); + } + break; + default: + return; + } + + CGAL_assertion(!is_vertex_coplanar || !is_vertex_opposite_coplanar); //coplanar edge are not forwarded + + + if ( is_vertex_coplanar ) + { + //grab original vertex that is on commom intersection + typename Hedge_to_polyhedron_map::iterator it=hedge_to_polyhedron.find(principal_edge->facet()->halfedge()); + poly_to_vertices_on_inter[it->second].insert(std::make_pair(node_id,principal_edge)); + polyhedron_to_map_node_to_polyhedron_vertex[it->second].insert(std::make_pair(node_id,principal_edge->vertex())); + all_incident_faces_got_a_node_as_vertex(principal_edge,node_id); + check_node_on_non_manifold_vertex(node_id,principal_edge); + } + else{ + if ( is_vertex_opposite_coplanar ){ + //grab original vertex that is on commom intersection + typename Hedge_to_polyhedron_map::iterator it=hedge_to_polyhedron.find(principal_edge->facet()->halfedge()); + poly_to_vertices_on_inter[it->second].insert(std::make_pair(node_id,principal_edge->opposite())); + polyhedron_to_map_node_to_polyhedron_vertex[it->second].insert(std::make_pair(node_id,principal_edge->opposite()->vertex())); + all_incident_faces_got_a_node_as_vertex(principal_edge->opposite(),node_id); + check_node_on_non_manifold_vertex(node_id,principal_edge->opposite()); + } + else{ + //handle intersection on principal edge + typename In_halfedge_map::iterator it_hedge_map=in_hedge.insert(std::make_pair(principal_edge,Node_id_set())).first; + it_hedge_map->second.insert(node_id); + check_node_on_non_manifold_edge(node_id,principal_edge); + } + } + } + + template + void annotate_graph(Iterator begin,Iterator end) + { +// std::cout << "Annotation graph..." << std::endl; + for (Iterator it=begin;it!=end;++it) + { + int node_id=it->first; + if (non_manifold_nodes.find(node_id)!=non_manifold_nodes.end()) it->second.make_terminal(); + const std::set& neighbors = it->second.neighbors; + graph_of_constraints.insert(std::make_pair(node_id,neighbors)); + } + } + + void update_terminal_nodes(std::vector&) + { + CGAL_assertion(!"Must not call this function"); + } + + void add_filtered_intersection(Halfedge_handle eh,Halfedge_handle fh,Polyhedron& Pe,Polyhedron& Pf){ + //use the representant halfedge of the facet as key + //--set polyhedron for the two facets incident to the edge + hedge_to_polyhedron.insert(std::make_pair(eh->facet()->halfedge(),&Pe)); + if ( !eh->opposite()->is_border() ) + hedge_to_polyhedron.insert(std::make_pair(eh->opposite()->facet()->halfedge(),&Pe)); + //--set polyhedron for the facet intersected by the edge + hedge_to_polyhedron.insert(std::make_pair(fh->facet()->halfedge(),&Pf)); + } + + + struct Polyhedron_face_boundary{ + std::vector node_ids_array[3]; // the node_ids on each halfedges + std::map hedges_ids; + Halfedge_handle halfedges[3]; //the three halfedges of the original face + Vertex_handle vertices[3]; //the three vertices of the original face + //node_ids_array[0] corresponds to the original edge vertices[0],vertices[1] = halfedges[0] + //node_ids_array[1] corresponds to the original edge vertices[1],vertices[2] = halfedges[1] + //node_ids_array[2] corresponds to the original edge vertices[2],vertices[0] = halfedges[2] + Polyhedron_face_boundary(Halfedge_handle first) + { + CGAL_assertion(first->next()->next()->next()==first); //the face is a triangle + hedges_ids.insert(std::make_pair(first,0)); + hedges_ids.insert(std::make_pair(first->next(),1)); + hedges_ids.insert(std::make_pair(first->next()->next(),2)); + halfedges[0]=first; + halfedges[1]=first->next(); + halfedges[2]=first->next()->next(); + + vertices[0]=halfedges[0]->opposite()->vertex(); + vertices[1]=halfedges[1]->opposite()->vertex(); + vertices[2]=halfedges[2]->opposite()->vertex(); + } + + //used when object was create with hedge but opposite was used to split the original face + void update_original_halfedge(Halfedge_handle original,Halfedge_handle new_hedge) + { + typename std::map::iterator it_id=hedges_ids.find(original); + CGAL_assertion(it_id!=hedges_ids.end()); + int index=it_id->second; + CGAL_assertion(halfedges[index]==original); + hedges_ids.erase(it_id); + hedges_ids.insert(std::make_pair(new_hedge,index)); + halfedges[index]=new_hedge; + } + + template + void copy_node_ids(Halfedge_handle hedge,Iterator begin,Iterator end) + { + typename std::map::iterator it_id=hedges_ids.find(hedge); + CGAL_assertion(it_id!=hedges_ids.end()); + std::copy(begin,end,std::back_inserter(node_ids_array[it_id->second])); + } + }; + + + void start_new_polyline(int i, int j) + { + if ( i==j ) //case of a single point + { + //TAG SL001 + //nothing is done + return; + } + std::pair res= + an_edge_per_polyline.insert( + std::make_pair( make_sorted_pair(i,j), + std::make_pair( std::map(),std::make_pair(false,0)) ) + ); + CGAL_assertion(res.second); + last_polyline=res.first; + if ( i !=last_polyline->first.first ) + last_polyline->second.second.first=true; + } + + void add_node_to_polyline(int){ + ++(last_polyline->second.second.second); + } + + void new_input_polyhedron(Polyhedron& P) + { + #ifndef NDEBUG + std::pair res = + #endif + polyhedron_to_map_node_to_polyhedron_vertex.insert(std::make_pair( &P,Node_to_polyhedron_vertex_map() )); + #ifndef NDEBUG + CGAL_assertion(res.second == true); + #endif + } + + //1) split_halfedges and retriangulate faces with no intersection point interior to the facet + //2) retriangulate using a constrained Delaunay triangulation each triangle in each Polyhedron that contains at least + // one intersection point inside the facet + //3) mark polyhedron edges that are on the intersection + //4) create one output polyhedron per connected component of polyhedron, connected by an edge which is not an intersection edge + //5) import each piece into a common combinatorial map + //6) glue all the pieces together + template + void finalize(const Node_vector& nodes){ + //mark halfedge that are on the intersection + //SL: I needed to use a map because to get the orientation around the edge, + // I need to know in the case the third vertex is a node its index (for exact construction) + typedef std::map,Cmp_unik_ad > Border_halfedges_map; + Border_halfedges_map border_halfedges; + + //Additionnal stuct to mark halfedge of the original polyhedra that are on the intersection + typedef Halfedge_marker Marker; + + //store for each triangle facet which boundary is intersected by the other surface, + //original vertices (and halfedges in the refined mesh pointing on these vertices) + typedef std::map Faces_boundary; + Faces_boundary faces_boundary; + + //0) For each polyhedron, collect original vertices that belongs to the intersection. + // From the graph of constaints, extract intersection edges that are incident to such vertices. In case + // there exists another original vertex adjacent to the first one found, this halfedge must be + // marked on the boundary (and possibly update an_edge_per_polyline). + // This is done first to avoid halfedges stored to be modified in the steps following. + for (typename Poly_to_vertices_on_intersection_map::iterator + it=poly_to_vertices_on_inter.begin(); + it!=poly_to_vertices_on_inter.end(); + ++it) + { + Polyhedron* poly=it->first; + std::set > already_done; + Node_to_target_of_hedge_map& nodes_to_hedge=it->second; + for(typename Node_to_target_of_hedge_map::iterator + it_node_2_hedge=nodes_to_hedge.begin(); + it_node_2_hedge!=nodes_to_hedge.end(); + ++it_node_2_hedge) + { + int node_id_of_first=it_node_2_hedge->first; + std::map< int,std::set >::iterator it_neighbors=graph_of_constraints.find(node_id_of_first); + if ( it_neighbors!=graph_of_constraints.end() ) + { + std::set& neighbors=it_neighbors->second; + for (std::set::iterator it_id=neighbors.begin();it_id!=neighbors.end();++it_id){ + if ( already_done.find(std::make_pair(*it_id,node_id_of_first))!=already_done.end() ) continue;//already done for the opposite + typename Node_to_target_of_hedge_map::iterator it_node_2_hedge_two=nodes_to_hedge.find(*it_id); + if ( it_node_2_hedge_two!=nodes_to_hedge.end() ) //a full edge is on intersection + { + //get the corresponding halfedge with vertex corresponding to node_id_of_first + Halfedge_handle hedge=it_node_2_hedge->second; + #ifndef NDEBUG + Halfedge_handle start=hedge; + #endif + while ( hedge->opposite()->vertex()!=it_node_2_hedge_two->second->vertex() ){ + hedge=hedge->next()->opposite(); + #ifndef NDEBUG + CGAL_assertion(hedge!=start); + #endif + } + std::pair edge_pair(*it_id,node_id_of_first); + border_halfedges.insert( std::make_pair(hedge,edge_pair) ); + Marker::mark(hedge); + update_edge_per_polyline(poly,edge_pair,hedge); + //save the fact that we already handle this edge + already_done.insert(std::make_pair(node_id_of_first,*it_id)); + } + } + } + #ifdef CGAL_COREFINEMENT_DEBUG + else + { + std::cout << "X1: Found an isolated point" << std::endl; + } + #endif + } + } + + //1) First split halfedges cut by the intersection polyline(s) + for (typename In_halfedge_map::iterator it=in_hedge.begin();it!=in_hedge.end();++it) + { + Halfedge_handle hedge=it->first; //the halfedge to be split (and its opposite too) + Node_ids node_ids; //indices of the intersection points to be inserted + //we used a set to avoid having duplicated nodes reported on an edge of two coplanar triangles + std::copy(it->second.begin(),it->second.end(),std::back_inserter(node_ids)); + typename Hedge_to_polyhedron_map::iterator it_poly=hedge_to_polyhedron.find( hedge->facet()->halfedge() ); + CGAL_assertion(it_poly!=hedge_to_polyhedron.end()); + Polyhedron* P=it_poly->second; //the polyhedron in which vertices should be added + + sort_vertices_along_hedge(node_ids,hedge,nodes); + + //save original face and nodes for face of hedge (1) + if ( !hedge->is_border() ){ + typename Faces_boundary::iterator it_face=faces_boundary.find(hedge->face()); + if (it_face==faces_boundary.end()) + it_face=faces_boundary.insert(std::make_pair(hedge->face(),Polyhedron_face_boundary(hedge))).first; + it_face->second.copy_node_ids(hedge,node_ids.begin(),node_ids.end()); + } + + //save original face and nodes for face of hedge->opposite (2) + typename Faces_boundary::iterator opposite_original_info=faces_boundary.end(); + if ( !hedge->opposite()->is_border() ){ + opposite_original_info=faces_boundary.find(hedge->opposite()->face()); + if (opposite_original_info==faces_boundary.end()) + opposite_original_info=faces_boundary.insert(std::make_pair(hedge->opposite()->face(),Polyhedron_face_boundary(hedge->opposite()))).first; + opposite_original_info->second.copy_node_ids(hedge->opposite(),node_ids.rbegin(),node_ids.rend()); + } + + typename Poly_to_map_node::iterator it_map=polyhedron_to_map_node_to_polyhedron_vertex.find(P); + CGAL_assertion(it_map!=polyhedron_to_map_node_to_polyhedron_vertex.end()); + //a map to identify the vertex in the polyhedron corresponding to an intersection point + Node_to_polyhedron_vertex_map& node_to_polyhedron_vertex=it_map->second; + + #ifndef NDEBUG + Vertex_handle original_vertex=hedge->opposite()->vertex(); + #endif + + //We need an edge incident to the source vertex of hedge. This is the first opposite edge created. + bool first=true; Halfedge_handle hedge_incident_to_src; + //do split the edges + for (std::vector::const_iterator it_id=node_ids.begin();it_id!=node_ids.end();++it_id){ + Vertex_handle v=split_edge(hedge,nodes[*it_id],*P); + node_to_polyhedron_vertex.insert(std::make_pair(*it_id,v)); + if (first){ + first=false; + hedge_incident_to_src=hedge->opposite()->next(); + } + } + + #ifndef NDEBUG + CGAL_assertion(hedge_incident_to_src->vertex()==original_vertex); + CGAL_assertion(hedge_incident_to_src->face()==hedge->opposite()->face()); + #endif + + //save original face and nodes for face of hedge->opposite (2) + if ( !hedge->opposite()->is_border() ){ + CGAL_assertion(opposite_original_info!=faces_boundary.end()); + opposite_original_info->second.update_original_halfedge(hedge->opposite(),hedge_incident_to_src); + } + + //insert the two incident faces in in_face map so that they will be triangulated. + if (!hedge->is_border()) in_face.insert(std::make_pair(hedge->face(),Node_ids())); + if (!hedge->opposite()->is_border()) in_face.insert(std::make_pair(hedge->opposite()->face(),Node_ids())); + } + + //2)triangulation of the triangle faces containing intersection point in their interior + // and also those with intersection points only on the boundary. + for (typename In_face_map::iterator it=in_face.begin();it!=in_face.end();++it) + { + Face_handle f = it->first; //the face to be retriangulated + Node_ids& node_ids = it->second; //the index of the intersection point that are interior to the face + typename Faces_boundary::iterator it_fb=faces_boundary.find(f); + + + typename Hedge_to_polyhedron_map::iterator it_polyhedron = hedge_to_polyhedron.find (f->halfedge()); //we can do this because the halfedge is still the same (at least its address)+no Face::set_halfedge called + CGAL_assertion(it_polyhedron != hedge_to_polyhedron.end()); + Polyhedron* P=it_polyhedron->second; + typename Poly_to_map_node::iterator it_map=polyhedron_to_map_node_to_polyhedron_vertex.find(P); + CGAL_assertion(it_map!=polyhedron_to_map_node_to_polyhedron_vertex.end()); + //a map to identify the vertex in the polyhedron corresponding to an intersection point + Node_to_polyhedron_vertex_map& node_to_polyhedron_vertex=it_map->second; + + std::map id_to_CDT_vh; + + //associate an edge of the triangulation to a halfedge in a given polyhedron + std::map,Halfedge_handle> edge_to_hedge; + + Vertex_handle triangle_boundary[3]; + int triangle_boundary_indices[3]; //the node_id of the triangle original vertex or a fake id + if (it_fb!=faces_boundary.end()){ //the boundary of the triangle face was refined + triangle_boundary[0]=it_fb->second.vertices[0]; + triangle_boundary[1]=it_fb->second.vertices[1]; + triangle_boundary[2]=it_fb->second.vertices[2]; + set_triangle_boundary_indices(triangle_boundary,triangle_boundary_indices); + } + else{ + triangle_boundary[0]=f->halfedge()->vertex(); //-1 + triangle_boundary[1]=f->halfedge()->next()->vertex(); //-2 + triangle_boundary[2]=f->halfedge()->next()->next()->vertex(); //-3 + CGAL_assertion(f->halfedge()->next()->next()->next()==f->halfedge());//check this is a triangle + set_triangle_boundary_indices(triangle_boundary,triangle_boundary_indices); + edge_to_hedge.insert (std::make_pair( std::make_pair( triangle_boundary_indices[2],triangle_boundary_indices[0] ) , f->halfedge() ) ); + edge_to_hedge.insert (std::make_pair( std::make_pair( triangle_boundary_indices[0],triangle_boundary_indices[1] ) , f->halfedge()->next() ) ); + edge_to_hedge.insert (std::make_pair( std::make_pair( triangle_boundary_indices[1],triangle_boundary_indices[2] ) , f->halfedge()->next()->next() ) ); + } + + + #ifdef DO_NO_USE_EXACT_CDT + typename Kernel::Plane_3 plane(triangle_boundary[0]->point(),triangle_boundary[1]->point(),triangle_boundary[2]->point()); + #else + CGAL::Cartesian_converter convert; + typename Exact_kernel::Plane_3 plane(convert(triangle_boundary[0]->point()),convert(triangle_boundary[1]->point()),convert(triangle_boundary[2]->point())); + #endif + CDT triangulation; + //insert point inside face + for (std::vector::iterator it_node_id=node_ids.begin();it_node_id!=node_ids.end();++it_node_id){ + #ifdef DO_NO_USE_EXACT_CDT + typename CDT::Vertex_handle vh=triangulation.insert(plane.to_2d(nodes[*it_node_id])); + #else + typename CDT::Vertex_handle vh=triangulation.insert(plane.to_2d(nodes.exact_node(*it_node_id))); + #endif + vh->info()=*it_node_id; + id_to_CDT_vh.insert(std::make_pair(*it_node_id,vh)); + } + + + typename CDT::Vertex_handle triangle_vertices[3]; + #ifdef DO_NO_USE_EXACT_CDT + triangle_vertices[0]=triangulation.insert(plane.to_2d(triangle_boundary[0]->point())); + triangle_vertices[1]=triangulation.insert(plane.to_2d(triangle_boundary[1]->point())); + triangle_vertices[2]=triangulation.insert(plane.to_2d(triangle_boundary[2]->point())); + #else + //we can do this because these are input points. + triangle_vertices[0]=triangulation.insert(plane.to_2d(convert(triangle_boundary[0]->point()))); + triangle_vertices[1]=triangulation.insert(plane.to_2d(convert(triangle_boundary[1]->point()))); + triangle_vertices[2]=triangulation.insert(plane.to_2d(convert(triangle_boundary[2]->point()))); + #endif + + triangle_vertices[0]->info()=triangle_boundary_indices[0]; + triangle_vertices[1]->info()=triangle_boundary_indices[1]; + triangle_vertices[2]->info()=triangle_boundary_indices[2]; + //insert face_extremities: we use operator[] because indice -1,-2,-3 are used in each loop and are specific to the current face + node_to_polyhedron_vertex[-1]=triangle_boundary[0]; + node_to_polyhedron_vertex[-2]=triangle_boundary[1]; + node_to_polyhedron_vertex[-3]=triangle_boundary[2]; + + //if one of the triangle original vertex is also a node + for (int ik=0;ik<3;++ik){ + if ( triangle_boundary_indices[ik]>=0 ) + id_to_CDT_vh.insert(std::make_pair(triangle_boundary_indices[ik],triangle_vertices[ik])); + } + //insert points on edges + #ifdef DO_NO_USE_EXACT_CDT + //and constrains these edges + #endif + if (it_fb!=faces_boundary.end()) //is there at least one intersection point on the boundary of the face? + { + //in the following loop, for each original edge of the triangle, we insert the constrained edges + // and we recover the halfedge_handle corresponding to these constrained (they are already in the polyhedron) + for (int i=0;i<3;++i){ +// std::cerr << "Boundary edges" << std::endl; +// std::cerr << " " << -1-i <second.node_ids_array[i]; + typename CDT::Vertex_handle previous=triangle_vertices[i]; + int previous_index=triangle_boundary_indices[i]; //index of original Polyhedron vertex + Halfedge_handle hedge = it_fb->second.halfedges[ (i+2) % 3]->next(); + CGAL_assertion( hedge->opposite()->vertex()==it_fb->second.vertices[i] ); + if (!bounding_ids.empty()){ //is there al least one intersection point on this edge? + for (Node_ids::iterator it_id=bounding_ids.begin();it_id!=bounding_ids.end();++it_id){ +// std::cerr << " "<< *it_id << std::endl; + #ifdef DO_NO_USE_EXACT_CDT + typename CDT::Vertex_handle vh=triangulation.insert(plane.to_2d(nodes[*it_id])); + #else + typename CDT::Vertex_handle vh=triangulation.insert(plane.to_2d(nodes.exact_node(*it_id))); + #endif + vh->info()=*it_id; + id_to_CDT_vh.insert(std::make_pair(*it_id,vh)); + #ifdef DO_NO_USE_EXACT_CDT + triangulation.insert_constraint(previous,vh); + #endif + edge_to_hedge.insert (std::make_pair( std::make_pair(previous_index,*it_id),hedge) ); + previous=vh; + hedge=hedge->next(); + previous_index=*it_id; + } + } + else{ + CGAL_assertion( it_fb->second.halfedges[i]->vertex() == it_fb->second.vertices[ (i+1) % 3 ] ); + CGAL_assertion( it_fb->second.halfedges[i]->opposite()->vertex() == it_fb->second.vertices[ i ] ); + } + CGAL_assertion(hedge==it_fb->second.halfedges[i]); + edge_to_hedge.insert (std::make_pair( std::make_pair(previous_index,triangle_boundary_indices[(i+1) % 3]) , it_fb->second.halfedges[i] ) ); +// std::cerr << " " << -1 - ( (i+1) % 3 ) < > constrained_edges; + + //insert constraints that are interior to the triangle (in the case no edges are collinear in the meshes) + insert_constrained_edges(node_ids,triangulation,id_to_CDT_vh,constrained_edges); + + //insert constraints between points that are on the boundary (not a contrained on the triangle boundary) + if (it_fb!=faces_boundary.end()) //is there at least one intersection point on the boundary of the face? + { + for (int i=0;i<3;++i){ + Node_ids& bounding_ids=it_fb->second.node_ids_array[i]; + insert_constrained_edges(bounding_ids,triangulation,id_to_CDT_vh,constrained_edges,true); + } + } + + //insert coplanar edges for endpoints of triangles + for (int i=0;i<3;++i){ + int nindex=triangle_vertices[i]->info(); + if ( nindex >=0 ) + insert_constrained_edges_coplanar_case(nindex,triangulation,id_to_CDT_vh); + } + + //XSL_TAG_CPL_VERT + //collect edges incident to a point that is the intersection of two coplanar faces. + //This ensure that triangulations are compatible. + if (it_fb!=faces_boundary.end()) //is there at least one intersection point on the boundary of the face? + { + for (typename CDT::Finite_vertices_iterator vit=triangulation.finite_vertices_begin(), + vit_end=triangulation.finite_vertices_end();vit_end!=vit;++vit) + { + //skip original vertices (that are not nodes) and non-coplanar facet issued vertices + //(this is working because intersection points between coplanar facets are the first inserted) + if ( vit->info() < 0 || vit->info() >= number_coplanar_vertices) continue; + std::map< int,std::set >::iterator res=coplanar_constraints.insert(std::make_pair(vit->info(),std::set())).first; + //turn around the vertex and get incident edge + typename CDT::Edge_circulator start=triangulation.incident_edges(vit); + typename CDT::Edge_circulator curr=start; + do{ + if (triangulation.is_infinite(*curr) ) continue; + typename CDT::Edge mirror_edge=triangulation.mirror_edge(*curr); + if ( triangulation.is_infinite( curr->first->vertex(curr->second) ) || + triangulation.is_infinite( mirror_edge.first->vertex(mirror_edge.second) ) ) + continue; //skip edges that are on the boundary of the triangle (these are already constrained) + //insert edges in the set of constraints + int nindex = + curr->first->vertex( (curr->second+1)%3 )==static_cast(vit)? + (curr->second+2)%3:(curr->second+1)%3; + typename CDT::Vertex_handle vn=curr->first->vertex(nindex); + if ( vit->info() > vn->info() ) continue; //take only one out of the two edges + skip negative vn->info() + CGAL_assertion(vn->info()>=0); + res->second.insert( vn->info() ).second; + }while(start!=++curr); + } + +// this is a working alternative that should be slower +// for (typename CDT::Finite_edges_iterator eit=triangulation.finite_edges_begin(), +// eit_end=triangulation.finite_edges_end();eit_end!=eit;++eit) +// { +// typename CDT::Edge mirror_edge=triangulation.mirror_edge(*eit); +// if ( triangulation.is_infinite( eit->first->vertex(eit->second) ) || +// triangulation.is_infinite( mirror_edge.first->vertex(mirror_edge.second) ) ) +// continue; //skip edges that are on the boundary of the triangle (these are already constrained) +// typename CDT::Vertex_handle v1=eit->first->vertex( (eit->second+1)%3 ), +// v2=eit->first->vertex( (eit->second+2)%3 ); +// if (v1->info()<0 || v2->info()<0) continue; +// if ( v1->info() > v2->info() ) std::swap(v1,v2); +// coplanar_constraints.insert(std::make_pair(v1->info(),std::set())).first->second.insert(v2->info()); +// } + } + + + //create a modifier to insert nodes and copy the triangulation of the face + //inside the polyhedron + internal_IOP::Triangulate_a_face modifier(f,nodes,node_ids,node_to_polyhedron_vertex,edge_to_hedge,triangulation); + + CGAL_assertion(P->is_valid()); + P->delegate(modifier); + CGAL_assertion(P->is_valid()); + + //3) mark halfedges that are common to two polyhedral surfaces + //recover halfedges inserted that are on the intersection + for (std::list >::iterator it_cst=constrained_edges.begin();it_cst!=constrained_edges.end();++it_cst) + { + typename std::map,Halfedge_handle>::iterator it_poly_hedge=edge_to_hedge.find(*it_cst); + //we cannot have an assertion here in the case an edge or part of an edge is a constraints. + //Indeed, the graph_of_constraints report an edge 0,1 and 1,0 for example while only one of the two + //is defined as one of them defines an adjacent face + //CGAL_assertion(it_poly_hedge!=edge_to_hedge.end()); + if( it_poly_hedge!=edge_to_hedge.end() ){ + border_halfedges.insert( std::make_pair(Halfedge_const_handle(it_poly_hedge->second),*it_cst) ); + Marker::mark(it_poly_hedge->second); + update_edge_per_polyline(P,it_poly_hedge->first,it_poly_hedge->second); + } + else{ + //WARNING: in few case this is needed if the marked edge is on the border + //to optimize it might be better to only use sorted pair. TAG_SLXX1 + std::pair opposite_pair(it_cst->second,it_cst->first); + it_poly_hedge=edge_to_hedge.find(opposite_pair); + CGAL_assertion( it_poly_hedge!=edge_to_hedge.end() ); + + border_halfedges.insert( std::make_pair(Halfedge_const_handle(it_poly_hedge->second),opposite_pair) ); + Marker::mark(it_poly_hedge->second); + update_edge_per_polyline(P,it_poly_hedge->first,it_poly_hedge->second); + } + } + } + + if (do_not_build_cmap) return; + + //4) create one output polyhedron per connected component of polyhedron, + // connected by an edge which is not an intersection edge + //5) import into a Combinatorial map + #ifdef CGAL_COREFINEMENT_DEBUG + std::cout << "Nb marked edges " << border_halfedges.size() << std::endl; +// for (typename Border_halfedges_map::iterator it=border_halfedges.begin();it!=border_halfedges.end();++it) +// std::cout << it->first->opposite()->vertex()->point() << " " << it->first->vertex()->point() << " is constrained " << std::endl; + std::cout << "Nb polylines " << an_edge_per_polyline.size() << std::endl; + #endif + + internal_IOP::Non_intersection_halfedge criterium(border_halfedges); + + int mark_index=final_map().get_new_mark(); //mark used to tag dart that are on an intersection + + //define a map that will contain the correspondance between selected halfedges of the boundary and + //their corresponding Dart_handle in the cmap. + typedef std::map > Halfedge_to_dart_map; + Halfedge_to_dart_map selected_hedge_to_dart; + for (typename An_edge_per_polyline_map::iterator it=an_edge_per_polyline.begin();it!=an_edge_per_polyline.end();++it) + { + CGAL_assertion(it->second.first.size()==2); + //orientation of faces around the edge (to be sure we can do it) + Halfedge_handle first_hedge=it->second.first.begin()->second; + Halfedge_handle second_hedge=boost::next(it->second.first.begin())->second; + + if (!first_hedge->is_border()) selected_hedge_to_dart.insert(std::make_pair(first_hedge,Dart_handle(NULL))); + if (!first_hedge->opposite()->is_border()) selected_hedge_to_dart.insert(std::make_pair(first_hedge->opposite(),Dart_handle(NULL))); + if (!second_hedge->is_border()) selected_hedge_to_dart.insert(std::make_pair(second_hedge,Dart_handle(NULL))); + if (!second_hedge->opposite()->is_border()) selected_hedge_to_dart.insert(std::make_pair(second_hedge->opposite(),Dart_handle(NULL))); + } + + #ifdef CGAL_COREFINEMENT_DEBUG + int polynb=0; + #endif + for (typename Poly_to_map_node::iterator + it=polyhedron_to_map_node_to_polyhedron_vertex.begin(); + it!=polyhedron_to_map_node_to_polyhedron_vertex.end(); + ++it + ) + { + typedef typename Polyhedron::Facet_const_handle Facet_const_handle; + typedef ::CGAL::Union_find UF; + typedef typename UF::handle UF_handle; + typedef typename UF::iterator UF_iterator; + typedef std::map,internal::Compare_handle_ptr > Result; + typedef std::map > Facet_to_handle_map; + + UF uf; + Facet_to_handle_map map_f2h; + Result result; + Polyhedron* current_poly=it->first; + + #ifdef CGAL_COREFINEMENT_DEBUG + std::cout << "writing poly debug"<< std::endl; + std::stringstream ss; + ss << "output_debug-" << ++polynb << ".off"; + std::ofstream output_debug(ss.str().c_str()); + output_debug << *current_poly; + #endif + + extract_connected_components(*(current_poly),criterium,uf,map_f2h,result); + + + //add each connected component in the map with 2 volumes per component. + for (typename Result::iterator it_res=result.begin();it_res!=result.end();++it_res) + { + //create in the final Cmap a 2D component containing faces of a connected component + //(twice: one with same orientation and one with the opposite orientation to model the other volume) + Dart_handle d = import_from_polyhedron_subset(final_map(),it_res->second.begin(),it_res->second.end(),criterium,selected_hedge_to_dart,mark_index); + //set an attribute to one volume represented by this component to indicates + //a part outside of the polyhedron current_poly + typename Combinatorial_map_3_::template Attribute_range<3>::type::iterator attrib=final_map().template create_attribute<3>(); + attrib->info().outside.insert(current_poly); + final_map().template set_attribute<3>(d,attrib); + //set the attribute for the opposite volume: represent a part inside current_poly + attrib=final_map().template create_attribute<3>(); + attrib->info().inside.insert(current_poly); + final_map().template set_attribute<3>(d->beta(3),attrib); + + #ifdef CGAL_COREFINEMENT_DEBUG + final_map().display_characteristics(std::cout); + std::cout << std::endl; + #endif + } + } + #ifndef NDEBUG + for(typename Halfedge_to_dart_map::iterator it=selected_hedge_to_dart.begin();it!=selected_hedge_to_dart.end();++it) + CGAL_assertion(it->second!=Dart_handle(NULL)); + #endif + + CGAL_assertion(final_map().is_valid()); + + std::set darts_to_remove; + + //6) Glue pieces together + // using one edge per intersection polyline, we merge the different volumes + for (typename An_edge_per_polyline_map::iterator it=an_edge_per_polyline.begin();it!=an_edge_per_polyline.end();++it) + { + CGAL_assertion(it->second.first.size()==2); + //orientation of faces around the edge (to be sure we can do it) + std::pair indices = it->first; + const std::pair& polyline_info=it->second.second; + + //get the two halfedges incident to the edge [indices.first,indices.second] + Halfedge_handle first_hedge=it->second.first.begin()->second; + Halfedge_handle second_hedge=boost::next(it->second.first.begin())->second; + + CGAL_assertion(nodes[indices.second]==first_hedge->vertex()->point()); + CGAL_assertion(nodes[indices.first]==first_hedge->opposite()->vertex()->point()); + CGAL_assertion(nodes[indices.second]==second_hedge->vertex()->point()); + CGAL_assertion(nodes[indices.first]==second_hedge->opposite()->vertex()->point()); + + Polyhedron* first_poly = it->second.first.begin()->first; + Polyhedron* second_poly = boost::next(it->second.first.begin())->first; + + //different handling depending on the number of incident triangles to the edge. + //After sewing there are two,three or four volumes if there are two,three or four incident triangles respectively + if ( first_hedge->is_border() || first_hedge->opposite()->is_border() ){ + if (second_hedge->is_border() || second_hedge->opposite()->is_border()) + sew_2_two_volumes_case(first_hedge,second_hedge,border_halfedges,selected_hedge_to_dart,mark_index,darts_to_remove,nodes,indices,polyline_info); + else + sew_2_three_volumes_case(second_hedge, first_hedge,indices,nodes,border_halfedges,selected_hedge_to_dart,second_poly,first_poly,mark_index,darts_to_remove,polyline_info); + } + else + if (second_hedge->is_border() || second_hedge->opposite()->is_border()) + sew_2_three_volumes_case(first_hedge, second_hedge,indices,nodes,border_halfedges,selected_hedge_to_dart,first_poly,second_poly,mark_index,darts_to_remove,polyline_info); + else + { + //Sort the four triangle facets around their common edge + // we suppose that the exterior of the polyhedron is indicated by + // counterclockwise oriented facets. + Vertex_handle P1=first_hedge->opposite()->next()->vertex(); + Vertex_handle P2=first_hedge->next()->vertex(); + // when looking from the side of indices.second, the interior of the first polyhedron is described + // by turning counterclockwise from P1 to P2 + Vertex_handle Q1=second_hedge->opposite()->next()->vertex(); + Vertex_handle Q2=second_hedge->next()->vertex(); + // when looking from the side of indices.second, the interior of the second polyhedron is described + // by turning from Q1 to Q2 + + //check if the third point of each triangular face is an original point (stay -1) + //or a intersection point (in that case we need the index of the corresponding node to + //have the exact value of the point) + int index_p1=node_index_of_incident_vertex(first_hedge->opposite()->next(),border_halfedges); + int index_p2=node_index_of_incident_vertex(first_hedge->next(),border_halfedges); + int index_q1=node_index_of_incident_vertex(second_hedge->opposite()->next(),border_halfedges); + int index_q2=node_index_of_incident_vertex(second_hedge->next(),border_halfedges); + + #ifdef CGAL_COREFINEMENT_DEBUG + std::cout << index_p1 << " " << index_p2 << " " << index_q1 << " " <point() << " | " << P2->point() << " | " << Q1->point() << " | " <point() << std::endl; + #endif + + if ( coplanar_triangles_case_handled(first_hedge,second_hedge,indices,nodes,index_p1,index_p2,index_q1,index_q2,selected_hedge_to_dart,mark_index,darts_to_remove,polyline_info) ) + continue; + + CGAL_assertion(P1->point() !=Q1->point() && P1->point()!=Q2->point() && P2->point() !=Q1->point() && P2->point()!=Q2->point()); + + bool Q1_is_between_P1P2 = filtered_order_around_edge(indices.first,indices.second,index_p1,index_p2,index_q1,P1,P2,Q1,nodes); + bool Q2_is_between_P1P2 = filtered_order_around_edge(indices.first,indices.second,index_p1,index_p2,index_q2,P1,P2,Q2,nodes); + + + //Recover the dart that will be the start point of the different sewing + // dof_X_outside = dart of face of , meaning the triangle containing the + // point X and part of the volume outside of the corresponding polyhedron + //-----first polyhedron + Dart_handle dof_P1_outside = get_associated_dart(first_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_P2_outside = get_associated_dart(first_hedge,selected_hedge_to_dart); + //-----second polyhedron + Dart_handle dof_Q1_outside = get_associated_dart(second_hedge->opposite(),selected_hedge_to_dart); + Dart_handle dof_Q2_outside = get_associated_dart(second_hedge,selected_hedge_to_dart); + + if ( Q1_is_between_P1P2 ){ + if( Q2_is_between_P1P2 ) + { + bool P1_is_between_Q1Q2 = filtered_order_around_edge(indices.first,indices.second,index_q1,index_q2,index_p1,Q1,Q2,P1,nodes); + if (!P1_is_between_Q1Q2){ + // poly_first - poly_second = P1Q1 U Q2P2 + // poly_second - poly_first = {0} + // poly_first \cap poly_second = Q1Q2 + // opposite( poly_first U poly_second ) = P2P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //P1Q1 + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q2P2 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q1Q2 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //P2P1 + //update inside outside info (because darts from the same volume have been merged) + dof_Q1_outside->beta(3)->template attribute<3>()->info().inside.insert(first_poly); //update Q1Q2 inside poly + dof_P2_outside->template attribute<3>()->info().outside.insert(second_poly);//update P2P1 outside poly + } + else{ + // poly_first - poly_second = Q2Q1 + // poly_second - poly_first = P2P1 + // poly_first \cap poly_second = P1Q2 U Q1P2 + // opposite( poly_first U poly_second ) = {O} + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //Q2Q1 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //P2P1 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q1P2 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1Q2 + //update inside outside info (because darts from the same volume have been merged) + dof_Q2_outside->template attribute<3>()->info().inside.insert(first_poly); //update Q2Q1 inside poly + dof_P2_outside->template attribute<3>()->info().inside.insert(second_poly);//update P2P1 inside poly + } + } + else + { + // poly_first - poly_second = P1Q1 + // poly_second - poly_first = P2Q2 + // poly_first \cap poly_second = Q1P2 + // opposite( poly_first U poly_second ) = Q2P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //P1Q1 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P2Q2 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q1P2 + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q2P1 + } + } + else + { + if( Q2_is_between_P1P2 ) + { + // poly_first - poly_second = Q2P2 + // poly_second - poly_first = Q1P1 + // poly_first \cap poly_second = P1Q2 + // opposite( poly_first U poly_second ) = P2Q1 + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q2P2 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q1P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1Q2 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //P2Q1 + } + else + { + bool P1_is_between_Q1Q2 = filtered_order_around_edge(indices.first,indices.second,index_q1,index_q2,index_p1,Q1,Q2,P1,nodes); + if (!P1_is_between_Q1Q2){ + // poly_first - poly_second = P1P2 + // poly_second - poly_first = Q1Q2 + // poly_first \cap poly_second = {0} + // opposite( poly_first U poly_second ) = P2Q1 U Q2P1 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1P2 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //Q1Q2 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //P2Q1 + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q2P1 + //update inside outside info (because darts from the same volume have been merged) + dof_Q1_outside->beta(3)->template attribute<3>()->info().outside.insert(first_poly); //update Q1Q2 outside poly + dof_P1_outside->beta(3)->template attribute<3>()->info().outside.insert(second_poly);//update P2P1 outside poly + } + else{ + // poly_first - poly_second = {0} + // poly_second - poly_first = Q1P1 U P2Q2 + // poly_first \cap poly_second = P1P2 + // opposite( poly_first U poly_second ) = Q2Q1 + sew_2_marked_darts( final_map(),dof_Q1_outside->beta(3) , dof_P1_outside ,mark_index, nodes, indices, polyline_info); //Q1P1 + sew_2_marked_darts( final_map(),dof_P2_outside , dof_Q2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P2Q2 + sew_2_marked_darts( final_map(),dof_P1_outside->beta(3) , dof_P2_outside->beta(3) ,mark_index, nodes, indices, polyline_info); //P1P2 + sew_2_marked_darts( final_map(),dof_Q2_outside , dof_Q1_outside ,mark_index, nodes, indices, polyline_info); //Q2Q1 + //update inside outside info (because darts from the same volume have been merged) + dof_P1_outside->beta(3)->template attribute<3>()->info().inside.insert(second_poly); //update P1P2 inside poly + dof_Q2_outside->template attribute<3>()->info().outside.insert(first_poly);//update Q2Q1 outside poly + } + } + } + } + } + + #ifdef CGAL_COREFINEMENT_DEBUG + std::cout << "number of darts to remove: " << darts_to_remove.size() <::iterator itdart=darts_to_remove.begin(),end=darts_to_remove.end();itdart!=end;++itdart){ + final_map().erase_dart(*itdart); + } + + //remove empty volumes + typedef typename Combinatorial_map_3_::template Attribute_range<3>::type Volume_attribute_range; + Volume_attribute_range& ratrib=final_map().template attributes<3>(); + typename Volume_attribute_range::iterator curr=ratrib.begin(),end=ratrib.end(); + do{ + if (curr->info().is_empty) + final_map().template erase_attribute<3>(curr++); + else + ++curr; + } + while(curr!=end); + + CGAL_assertion(final_map().is_valid()); + + //update the info of each volume knowing about only one polyhedron: + //this happens when one polyhedron has a connected component + //that do not intersect the other polyhedron + typedef CGAL::Polyhedral_mesh_domain_3 Mesh_domain; + CGAL_precondition(polyhedron_to_map_node_to_polyhedron_vertex.size()==2); + Polyhedron* Poly_A = polyhedron_to_map_node_to_polyhedron_vertex.begin()->first; + Polyhedron* Poly_B = boost::next(polyhedron_to_map_node_to_polyhedron_vertex.begin())->first; + Mesh_domain* domain_A_ptr=NULL; + Mesh_domain* domain_B_ptr=NULL; + + #ifdef CGAL_COREFINEMENT_DEBUG + final_map().display_characteristics(std::cout); std::cout << "\n"; + #endif + + typename Combinatorial_map_3::template One_dart_per_cell_range<3> cell_range=final_map().template one_dart_per_cell<3>(); + for (typename Combinatorial_map_3::template One_dart_per_cell_range<3>::iterator + it = cell_range.begin(), it_end=cell_range.end(); + it_end!= it; + ++it ) + { + internal_IOP::Volume_info& info=it->template attribute<3>()->info(); + unsigned inside_size=info.inside.size(); + unsigned outside_size=info.outside.size(); + if ( inside_size + outside_size == 1) + { + bool is_inside = (inside_size==1); + Polyhedron* current_poly= is_inside? (*info.inside.begin()):(*info.outside.begin()); + Polyhedron* test_poly; + Mesh_domain* domain_ptr; + if ( current_poly==Poly_A) + { + test_poly=Poly_B; + if (domain_B_ptr == NULL) domain_B_ptr=new Mesh_domain(*Poly_B); + domain_ptr=domain_B_ptr; + } + else + { + test_poly=Poly_A; + if (domain_A_ptr == NULL) domain_A_ptr=new Mesh_domain(*Poly_A); + domain_ptr=domain_A_ptr; + } + typename Mesh_domain::Is_in_domain is_in_domain(*domain_ptr); + typename Kernel::Point_3 query=it->template attribute<0>()->point(); + if ( is_in_domain(query) ) + info.inside.insert(test_poly); + else + info.outside.insert(test_poly); + } + + #ifdef CGAL_COREFINEMENT_DEBUG + std::cout << "This volume has inside: "; + for (typename std::set::iterator itpoly=info.inside.begin();itpoly!=info.inside.end();++itpoly) + std::cout << " " << *itpoly; + std::cout << " and outside: "; + for (typename std::set::iterator itpoly=info.outside.begin();itpoly!=info.outside.end();++itpoly) + std::cout << " " << *itpoly; + std::cout << std::endl; + #endif + } + if (domain_A_ptr!=NULL) delete domain_A_ptr; + if (domain_B_ptr!=NULL) delete domain_B_ptr; + + } +}; + +}//namespace CGAL + + + +#endif //CGAL_INTERSECTION_OF_POLYHEDRA_3_REFINEMENT_VISITOR_H