WIP does not even compile

This commit is contained in:
Andreas Fabri 2019-11-08 16:21:20 +01:00
parent e18413c9b5
commit b6bc521c22
5 changed files with 99 additions and 60 deletions

View File

@ -24,7 +24,8 @@
#include <boost/random/uniform_int.hpp>
#include <boost/random/variate_generator.hpp>
#include <tbb/parallel_sort.h>
#include <tbb/task_group.h>
#include <algorithm>
#include <iterator>
#include <functional>
@ -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<int> 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

View File

@ -25,6 +25,8 @@
#include <CGAL/Box_intersection_d/Box_traits_d.h>
#include <CGAL/Box_intersection_d/box_limits.h>
#include <tbb/task_group.h>
#include <vector>
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<NT>::inf();
const NT sup = Box_intersection_d::box_limits<NT>::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();
}

View File

@ -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." )

View File

@ -2,6 +2,8 @@
#include <CGAL/Surface_mesh.h>
#include <CGAL/Polygon_mesh_processing/self_intersections.h>
#include <CGAL/Timer.h>
#include <CGAL/Real_timer.h>
#include <fstream>
@ -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<std::pair<face_descriptor, face_descriptor> > 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;

View File

@ -40,12 +40,17 @@
#include <CGAL/Polygon_mesh_processing/internal/named_function_params.h>
#include <CGAL/Polygon_mesh_processing/internal/named_params_helper.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#ifdef DOXYGEN_RUNNING
#define CGAL_PMP_NP_TEMPLATE_PARAMETERS NamedParameters
#define CGAL_PMP_NP_CLASS NamedParameters
#endif
namespace CGAL {
namespace internal {
template <class TM,//TriangleMesh
class Kernel,
@ -100,6 +105,38 @@ struct Intersect_facets
do_intersect_3_functor(kernel.do_intersect_3_object())
{ }
void operator()(const tbb::blocked_range<std::size_t> &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 <class TriangleMesh
, class OutputIterator
#ifdef DOXYGEN_RUNNING
, class NamedParameters
, class NamedParameters>
#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<P,T,R>& np)
, const Named_function_parameters<P,T,R>& 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<TM>::face_descriptor face_descriptor;
typedef typename boost::graph_traits<TM>::vertex_descriptor vertex_descriptor;
typedef typename CGAL::Box_intersection_d::Box_with_info_d<double, 3, face_descriptor> 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<TM, NamedParameters>::type GeomTraits;
Emptyset_iterator dev0;
CGAL::internal::Intersect_facets<TM,GeomTraits,Box,Emptyset_iterator,VertexPointMap>
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<std::size_t>(0, num_vertices(tmesh)), intersect_facets_parallel);
// compute self-intersections filtered out by boxes
CGAL::internal::Intersect_facets<TM,GeomTraits,Box,OutputIterator,VertexPointMap>
intersect_facets(tmesh, out, vpmap,
parameters::choose_parameter(parameters::get_parameter(np, internal_np::geom_traits), GeomTraits()));