diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/self_intersections_example.cpp b/Polygon_mesh_processing/examples/Polygon_mesh_processing/self_intersections_example.cpp index 440d7738bd5..179eeec96ca 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/self_intersections_example.cpp +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/self_intersections_example.cpp @@ -1,3 +1,5 @@ +struct Parallel_tag_bis {}; + #include #include diff --git a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h index 909a00381f3..c40df55bae6 100644 --- a/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h +++ b/Polygon_mesh_processing/include/CGAL/Polygon_mesh_processing/self_intersections.h @@ -41,9 +41,11 @@ #include #include +#ifdef CGAL_LINKED_WITH_TBB #include #include #include +#endif #ifdef DOXYGEN_RUNNING #define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters @@ -81,6 +83,7 @@ struct Intersect_facets typedef typename Kernel::Triangle_3 Triangle; typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; typedef typename boost::property_map::const_type Ppmap; // members @@ -108,7 +111,7 @@ struct Intersect_facets { } void operator()(const Box* b, const Box* c) const - { + { halfedge_descriptor h = halfedge(b->info(), m_tmesh); halfedge_descriptor g = halfedge(c->info(),m_tmesh); @@ -206,7 +209,10 @@ struct Intersect_facets } // end operator () }; // end struct Intersect_facets + +#ifdef CGAL_LINKED_WITH_TBB +// The functor for testing only triangles that do not share an edge or vertex in parallel template @@ -269,6 +275,146 @@ struct TriangleTriangle { }; +// The functor for doing all geometric tests in parallel +template +struct AllPairs { + const TM& m_tmesh; + const VertexPointMap m_vpmap; + mutable OutputIterator m_iterator; + typename Kernel::Construct_segment_3 segment_functor; + typename Kernel::Construct_triangle_3 triangle_functor; + typename Kernel::Do_intersect_3 do_intersect_3_functor; + typedef typename Kernel::Segment_3 Segment; + typedef typename Kernel::Triangle_3 Triangle; + + typedef typename boost::graph_traits::halfedge_descriptor halfedge_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + const std::vector >& seq_faces; + + + AllPairs(const std::vector >& seq_faces, + OutputIterator it, + const TM& tmesh, VertexPointMap vpmap, const Kernel& kernel) + : seq_faces(seq_faces), + m_tmesh(tmesh), + m_vpmap(vpmap), + m_iterator(it), + triangle_functor(kernel.construct_triangle_3_object()), + do_intersect_3_functor(kernel.do_intersect_3_object()) + {} + + void operator()(const tbb::blocked_range &r) const + { + for (std::size_t ri = r.begin(); ri != r.end(); ++ri) { + const std::pair& ff = seq_faces[ri]; + this->operator()(ff); + } + } + + void operator()(const std::pair& ff) const + { + halfedge_descriptor h = halfedge(ff.first, m_tmesh), g = halfedge(ff.second, m_tmesh); + + vertex_descriptor hv[3], gv[3]; + hv[0] = target(h, m_tmesh); + hv[1] = target(next(h, m_tmesh), m_tmesh); + hv[2] = source(h, m_tmesh); + + gv[0] = target(g, m_tmesh); + gv[1] = target(next(g, m_tmesh), m_tmesh); + gv[2] = source(g, m_tmesh); + + halfedge_descriptor opp_h; + + // check for shared egde + for(unsigned int i=0; i<3; ++i){ + opp_h = opposite(h, m_tmesh); + if(face(opp_h, m_tmesh) == ff.second){ + // there is an intersection if the four points are coplanar and + // the triangles overlap + get(m_vpmap, hv[i]); + get(m_vpmap, hv[(i + 1) % 3]); + get(m_vpmap, hv[(i + 2) % 3]); + get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh)); + + if(CGAL::coplanar(get(m_vpmap, hv[i]), + get(m_vpmap, hv[(i+1)%3]), + get(m_vpmap, hv[(i+2)%3]), + get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh))) && + CGAL::coplanar_orientation(get(m_vpmap, hv[(i+2)%3]), + get(m_vpmap, hv[i]), + get(m_vpmap, hv[(i+1)%3]), + get(m_vpmap, target(next(opp_h, m_tmesh), m_tmesh))) + == CGAL::POSITIVE){ + *m_iterator++ = ff; + return; + } else { // there is a shared edge but no intersection + return; + } + } + h = next(h, m_tmesh); + } + + // check for shared vertex --> maybe intersection, maybe not + + halfedge_descriptor v; + + int i(0),j(0); + bool shared = false; + for(; i < 3 && (! shared); ++i){ + for(j = 0; j < 3 && (! shared); ++j){ + if(hv[i]==gv[j]){ + shared = true; + break; + } + } + if (shared) { + break; + } + } + if(shared){ + // found shared vertex: + assert(hv[i] == gv[j]); + // geometric check if the opposite segments intersect the triangles + Triangle t1 = triangle_functor( get(m_vpmap, hv[0]), + get(m_vpmap, hv[1]), + get(m_vpmap, hv[2])); + Triangle t2 = triangle_functor( get(m_vpmap, gv[0]), + get(m_vpmap, gv[1]), + get(m_vpmap, gv[2])); + + Segment s1 = segment_functor( get(m_vpmap, hv[(i+1)%3]), + get(m_vpmap, hv[(i+2)%3]) ); + Segment s2 = segment_functor( get(m_vpmap, gv[(j+1)%3]), + get(m_vpmap, gv[(j+2)%3])); + + if(do_intersect_3_functor(t1,s2)){ + *m_iterator ++ = ff; + } else if(do_intersect_3_functor(t2,s1)){ + *m_iterator ++ = ff; + } + return; + } + + // check for geometric intersection + Triangle t1 = triangle_functor( get(m_vpmap, hv[0]), + get(m_vpmap, hv[1]), + get(m_vpmap, hv[2])); + Triangle t2 = triangle_functor( get(m_vpmap, gv[0]), + get(m_vpmap, gv[1]), + get(m_vpmap, gv[2])); + if(do_intersect_3_functor(t1, t2)){ + *m_iterator ++ = ff; + } + } +}; + + +// The functor for filtering pairs of faces that share an edge or vertex template @@ -339,6 +485,27 @@ struct Incident_faces_filter } // end operator () }; +// The functor that filters nothing +template +struct All_faces_filter +{ + mutable OutputIterator m_iterator; + + + All_faces_filter(OutputIterator it) + : m_iterator(it) + {} + + template + void operator()(const Box* b, const Box* c) const + { + *m_iterator ++ = std::make_pair(b->info(), c->info()); + + } // end operator () +}; + + +// The functor that computes intersections for faces incident to a vertex template box_ptr; - box_ptr.reserve(num_faces(tmesh)); + box_ptr.reserve(boxes.size()); for(Box& b : boxes) box_ptr.push_back(&b); @@ -630,6 +798,7 @@ self_intersections( const FaceRange& face_range, return intersect_facets.m_iterator; } +#ifdef CGAL_LINKED_WITH_TBB template +OutputIterator +self_intersections( const FaceRange& face_range, + const TriangleMesh& tmesh, + OutputIterator out, + const NamedParameters& np, + Parallel_tag) +{ + CGAL_precondition(CGAL::is_triangle_mesh(tmesh)); + + typedef TriangleMesh TM; + typedef typename boost::graph_traits::face_descriptor face_descriptor; + typedef typename boost::graph_traits::vertex_descriptor vertex_descriptor; + typedef typename CGAL::Box_intersection_d::Box_with_info_d Box; + + // make one box per facet + std::vector boxes; + boxes.reserve( + std::distance( boost::begin(face_range), boost::end(face_range) ) + ); + + typedef typename GetVertexPointMap::const_type VertexPointMap; + VertexPointMap vpmap = parameters::choose_parameter(parameters::get_parameter(np, internal_np::vertex_point), + get_const_property_map(boost::vertex_point, tmesh)); + + for(face_descriptor f : face_range) + { + typename boost::property_traits::reference + p = get(vpmap, target(halfedge(f,tmesh),tmesh)), + q = get(vpmap, target(next(halfedge(f, tmesh), tmesh), tmesh)), + r = get(vpmap, target(next(next(halfedge(f, tmesh), tmesh), tmesh), tmesh)); + + if ( collinear(p, q, r) ) + *out++= std::make_pair(f,f); + else + boxes.push_back(Box(p.bbox() + q.bbox() + r.bbox(), f)); + } + // generate box pointers + std::vector box_ptr; + box_ptr.reserve(boxes.size()); + + for(Box& b : boxes) + box_ptr.push_back(&b); + + // (A) Sequentially write all pairs of faces with intersecting bbox into a std::vector + std::ptrdiff_t cutoff = 2000; + typedef std::vector >SeqV; + typedef std::back_insert_iterator SeqVI; + SeqV face_pairs; + internal::All_faces_filter all_faces_filter(std::back_inserter(face_pairs)); + CGAL::box_self_intersection_d(box_ptr.begin(), box_ptr.end(),all_faces_filter,cutoff); + + + // (B) Parallel: perform the geometric tests + typedef typename GetGeomTraits::type GeomTraits; + + typedef tbb::concurrent_vector > CV; + CV cv_faces; + typedef std::back_insert_iterator CVI; + + + CGAL::internal::AllPairs + all_pairs(face_pairs, + std::back_inserter(cv_faces), // result + tmesh, + vpmap, + parameters::choose_parameter(parameters::get_parameter(np, internal_np::geom_traits), GeomTraits())); + + Real_timer rt; + rt.start(); + + tbb::parallel_for(tbb::blocked_range(0, face_pairs.size()), all_pairs); + + // (C) Sequentially: Copy from the concurent container to the output iterator + for(CV::iterator it = cv_faces.begin(); it != cv_faces.end(); ++it){ + *out ++= *it; + } + return out; +} +#endif // CGAL_LINKED_WITH_TBB + + /// \cond SKIP_IN_MANUAL template