From b6bc521c22dfd24e6393740eb9d896d2a1d964d1 Mon Sep 17 00:00:00 2001 From: Andreas Fabri Date: Fri, 8 Nov 2019 16:21:20 +0100 Subject: [PATCH] WIP does not even compile --- .../CGAL/Box_intersection_d/segment_tree.h | 65 ++++------------- .../include/CGAL/box_intersection_d.h | 9 ++- .../Polygon_mesh_processing/CMakeLists.txt | 1 + .../self_intersections_example.cpp | 12 +++- .../self_intersections.h | 72 +++++++++++++++++-- 5 files changed, 99 insertions(+), 60 deletions(-) diff --git a/Box_intersection_d/include/CGAL/Box_intersection_d/segment_tree.h b/Box_intersection_d/include/CGAL/Box_intersection_d/segment_tree.h index 74a113c46b8..d7514750c9b 100644 --- a/Box_intersection_d/include/CGAL/Box_intersection_d/segment_tree.h +++ b/Box_intersection_d/include/CGAL/Box_intersection_d/segment_tree.h @@ -24,7 +24,8 @@ #include #include - +#include +#include #include #include #include @@ -93,8 +94,8 @@ void one_way_scan( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, bool in_order = true ) { typedef typename Traits::Compare Compare; - std::sort( p_begin, p_end, Compare( 0 ) ); - std::sort( i_begin, i_end, Compare( 0 ) ); + tbb::parallel_sort( p_begin, p_end, Compare( 0 ) ); + tbb::parallel_sort( i_begin, i_end, Compare( 0 ) ); // for each box viewed as interval i for( RandomAccessIter2 i = i_begin; i != i_end; ++i ) { @@ -133,8 +134,8 @@ void modified_two_way_scan( { typedef typename Traits::Compare Compare; - std::sort( p_begin, p_end, Compare( 0 ) ); - std::sort( i_begin, i_end, Compare( 0 ) ); + tbb::parallel_sort( p_begin, p_end, Compare( 0 ) ); + tbb::parallel_sort( i_begin, i_end, Compare( 0 ) ); // for each box viewed as interval while( i_begin != i_end && p_begin != p_end ) { @@ -323,7 +324,8 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, RandomAccessIter2 i_begin, RandomAccessIter2 i_end, T lo, T hi, Callback callback, Predicate_traits traits, - std::ptrdiff_t cutoff, int dim, bool in_order ) + std::ptrdiff_t cutoff, int dim, bool in_order, + tbb::task_group& tg) { typedef typename Predicate_traits::Spanning Spanning; typedef typename Predicate_traits::Lo_less Lo_less; @@ -332,37 +334,10 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, const T inf = box_limits< T >::inf(); const T sup = box_limits< T >::sup(); -#if CGAL_BOX_INTERSECTION_DEBUG - CGAL_STATIC_THREAD_LOCAL_VARIABLE(int, level, -1); - Counter bla( level ); - CGAL_BOX_INTERSECTION_DUMP("range: [" << lo << "," << hi << ") dim " - << dim << std::endl ) - CGAL_BOX_INTERSECTION_DUMP("intervals: " ) - //dump_box_numbers( i_begin, i_end, traits ); - dump_intervals( i_begin, i_end, traits, dim ); - CGAL_BOX_INTERSECTION_DUMP("points: " ) - //dump_box_numbers( p_begin, p_end, traits ); - dump_points( p_begin, p_end, traits, dim ); -#endif - -#if CGAL_SEGMENT_TREE_CHECK_INVARIANTS - { - // first: each point is inside segment [lo,hi) - for( RandomAccessIter1 it = p_begin; it != p_end; ++it ) { - CGAL_assertion( Lo_less( hi, dim )(*it) ); - CGAL_assertion( Lo_less( lo, dim )(*it) == false ); - } - // second: each interval intersects segment [lo,hi) - for( RandomAccessIter2 it = i_begin; it != i_end; ++it ) - CGAL_assertion( Hi_greater(lo,dim)(*it) && Lo_less(hi,dim)(*it)); - } -#endif - if( p_begin == p_end || i_begin == i_end || lo >= hi ) return; if( dim == 0 ) { - CGAL_BOX_INTERSECTION_DUMP( "dim = 0. scanning ... " << std::endl ) one_way_scan( p_begin, p_end, i_begin, i_end, callback, traits, dim, in_order ); return; @@ -371,7 +346,6 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, if( std::distance( p_begin, p_end ) < cutoff || std::distance( i_begin, i_end ) < cutoff ) { - CGAL_BOX_INTERSECTION_DUMP( "scanning ... " << std::endl ) modified_two_way_scan( p_begin, p_end, i_begin, i_end, callback, traits, dim, in_order ); return; @@ -381,23 +355,18 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, std::partition( i_begin, i_end, Spanning( lo, hi, dim ) ); if( i_begin != i_span_end ) { - CGAL_BOX_INTERSECTION_DUMP( "checking spanning intervals ... " - << std::endl ) + // make two calls for roots of segment tree at next level. segment_tree( p_begin, p_end, i_begin, i_span_end, inf, sup, - callback, traits, cutoff, dim - 1, in_order ); + callback, traits, cutoff, dim - 1, in_order, tg ); segment_tree( i_begin, i_span_end, p_begin, p_end, inf, sup, - callback, traits, cutoff, dim - 1, !in_order ); + callback, traits, cutoff, dim - 1, !in_order, tg ); } T mi; RandomAccessIter1 p_mid = split_points( p_begin, p_end, traits, dim, mi ); if( p_mid == p_begin || p_mid == p_end ) { - CGAL_BOX_INTERSECTION_DUMP( "unable to split points! ") - //dump_points( p_begin, p_end, traits, dim ); - CGAL_BOX_INTERSECTION_DUMP( "performing modified two_way_san ... " - << std::endl ) modified_two_way_scan( p_begin, p_end, i_span_end, i_end, callback, traits, dim, in_order ); return; @@ -407,21 +376,17 @@ void segment_tree( RandomAccessIter1 p_begin, RandomAccessIter1 p_end, // separate left intervals. // left intervals have a low point strictly less than mi i_mid = std::partition( i_span_end, i_end, Lo_less( mi, dim ) ); - CGAL_BOX_INTERSECTION_DUMP("->left" << std::endl ) + segment_tree( p_begin, p_mid, i_span_end, i_mid, lo, mi, - callback, traits, cutoff, dim, in_order ); + callback, traits, cutoff, dim, in_order, tg ); // separate right intervals. // right intervals have a high point strictly higher than mi i_mid = std::partition( i_span_end, i_end, Hi_greater( mi, dim ) ); - CGAL_BOX_INTERSECTION_DUMP("->right"<< std::endl ) + segment_tree( p_mid, p_end, i_span_end, i_mid, mi, hi, - callback, traits, cutoff, dim, in_order ); + callback, traits, cutoff, dim, in_order, tg ); } -#if CGAL_BOX_INTERSECTION_DEBUG - #undef CGAL_BOX_INTERSECTION_DUMP -#endif -#undef CGAL_BOX_INTERSECTION_DEBUG } // end namespace Box_intersection_d diff --git a/Box_intersection_d/include/CGAL/box_intersection_d.h b/Box_intersection_d/include/CGAL/box_intersection_d.h index d7952400617..d2e3415ed8e 100644 --- a/Box_intersection_d/include/CGAL/box_intersection_d.h +++ b/Box_intersection_d/include/CGAL/box_intersection_d.h @@ -25,6 +25,8 @@ #include #include +#include + #include namespace CGAL { @@ -40,6 +42,7 @@ void box_intersection_custom_predicates_d( std::ptrdiff_t cutoff = 10, Box_intersection_d::Setting setting = Box_intersection_d::BIPARTITE) { + tbb::task_group tg; typedef BoxPredicateTraits Traits; typedef typename Traits::NT NT; CGAL_assertion( Traits::dimension() > 0 ); @@ -47,10 +50,12 @@ void box_intersection_custom_predicates_d( const NT inf = Box_intersection_d::box_limits::inf(); const NT sup = Box_intersection_d::box_limits::sup(); Box_intersection_d::segment_tree(begin1, end1, begin2, end2, - inf, sup, callback, traits, cutoff, dim, true); + inf, sup, callback, traits, cutoff, dim, true, tg); if(setting == Box_intersection_d::BIPARTITE) Box_intersection_d::segment_tree(begin2, end2, begin1, end1, - inf, sup, callback, traits, cutoff, dim, false); + inf, sup, callback, traits, cutoff, dim, false, tg); + + tg.wait(); } diff --git a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt index 88701b9e196..eeebbd3e6cc 100644 --- a/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt +++ b/Polygon_mesh_processing/examples/Polygon_mesh_processing/CMakeLists.txt @@ -108,6 +108,7 @@ endif(OpenMesh_FOUND) find_package( TBB ) if( TBB_FOUND ) + CGAL_target_use_TBB(self_intersections_example) CGAL_target_use_TBB(hausdorff_distance_remeshing_example) else() message( STATUS "NOTICE: Intel TBB was not found. Sequential code will be used." ) 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 93a2de84858..403d22786a0 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 @@ -2,6 +2,8 @@ #include #include +#include +#include #include @@ -22,17 +24,23 @@ int main(int argc, char* argv[]) std::cerr << "Not a valid input file." << std::endl; return 1; } - + /* bool intersecting = PMP::does_self_intersect(mesh, PMP::parameters::vertex_point_map(get(CGAL::vertex_point, mesh))); std::cout << (intersecting ? "There are self-intersections." : "There is no self-intersection.") << std::endl; + */ + CGAL::Real_timer rtimer; + CGAL::Timer timer; + rtimer.start(); + timer.start(); + std::vector > intersected_tris; PMP::self_intersections(mesh, std::back_inserter(intersected_tris)); - + std::cout << rtimer.time() << " " << timer.time() << " sec." << std::endl; std::cout << intersected_tris.size() << " pairs of triangles intersect." << std::endl; return 0; 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 decbb718f93..e18be4a35d4 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 @@ -40,12 +40,17 @@ #include #include +#include +#include + #ifdef DOXYGEN_RUNNING #define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters #define CGAL_PMP_NP_CLASS NamedParameters #endif namespace CGAL { + + namespace internal { template &r) const + { + for (std::size_t i = r.begin(); i != r.end(); ++i) { + vertex_descriptor vd(i); + if(degree(vd, m_tmesh)<= 3){ + return; + } + halfedge_descriptor hd = halfedge(vd, m_tmesh), done(hd); + + do { //for each hd with vd as target + // no need to look at the neighbor faces + // later do the shared edge test with them + Triangle th = triangle_functor( get(m_vpmap, vd), + get(m_vpmap, source(hd,m_tmesh)), + get(m_vpmap, target(next(hd,m_tmesh), m_tmesh)) ); + halfedge_descriptor start = prev(opposite(prev(opposite(hd,m_tmesh),m_tmesh),m_tmesh),m_tmesh); + halfedge_descriptor stop = opposite(next(hd,m_tmesh),m_tmesh); + while(start != stop){ + if(start < h){ + Segment ss = segment_functor(get(m_vpmap, source(start,m_tmesh)), + get(m_vpmap, target(next(start,m_tmesh), m_tmesh))); + if(do_intersect_3_functor(th, ss)){ + *m_iterator_wrapper++ = std::make_pair(face(h,m_tmesh), face(start, m_tmesh)); + } + } + start = prev(opposite(hd,m_tmesh),m_tmesh); + } + ++hd; + }while(hd != done); + } + } + void operator()(const Box* b, const Box* c) const { halfedge_descriptor h = halfedge(b->info(), m_tmesh); @@ -163,6 +200,7 @@ struct Intersect_facets } } if(shared){ +#if 0 // found shared vertex: assert(hv[i] == gv[j]); // geometric check if the opposite segments intersect the triangles @@ -183,6 +221,7 @@ struct Intersect_facets } else if(do_intersect_3_functor(t2,s1)){ *m_iterator_wrapper++ = std::make_pair(b->info(), c->info()); } +#endif return; } @@ -211,6 +250,8 @@ struct Throw_at_output { }// namespace internal + + namespace Polygon_mesh_processing { #ifndef DOXYGEN_RUNNING @@ -226,6 +267,7 @@ self_intersections( const FaceRange& face_range, const NamedParameters& np); #endif + /** * \ingroup PMP_intersection_grp * detects and records self-intersections of a triangulated surface mesh. @@ -253,22 +295,27 @@ self_intersections( const FaceRange& face_range, * * @return `out` */ + + template #else //avoid ambiguity with self_intersections(faces, tmesh, out) - , class P, class T, class R + , class P, class T, class R> #endif -> + + OutputIterator self_intersections(const TriangleMesh& tmesh , OutputIterator out #ifdef DOXYGEN_RUNNING - , const NamedParameters& np) + , const NamedParameters& np #else - , const Named_function_parameters& np) + , const Named_function_parameters& np #endif + ) + { return self_intersections(faces(tmesh), tmesh, out, np); } @@ -281,6 +328,8 @@ self_intersections(const TriangleMesh& tmesh, OutputIterator out) return self_intersections(tmesh, out, CGAL::Polygon_mesh_processing::parameters::all_default()); } + + /// \endcond /*! @@ -326,6 +375,7 @@ self_intersections( const FaceRange& face_range, 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 @@ -357,8 +407,18 @@ self_intersections( const FaceRange& face_range, for(Box& b : boxes) box_ptr.push_back(&b); - // compute self-intersections filtered out by boxes + typedef typename GetGeomTraits::type GeomTraits; + + Emptyset_iterator dev0; + CGAL::internal::Intersect_facets + intersect_facets_parallel(tmesh, dev0, vpmap, + parameters::choose_parameter(parameters::get_parameter(np, internal_np::geom_traits), GeomTraits())); + + std::cout << "do it"<< std::endl; + tbb::parallel_for(tbb::blocked_range(0, num_vertices(tmesh)), intersect_facets_parallel); + + // compute self-intersections filtered out by boxes CGAL::internal::Intersect_facets intersect_facets(tmesh, out, vpmap, parameters::choose_parameter(parameters::get_parameter(np, internal_np::geom_traits), GeomTraits()));