From e3a5f7c2bf5b4e9ceb2e5a89aeb926036d01f201 Mon Sep 17 00:00:00 2001 From: Ning Xu Date: Mon, 18 Aug 2014 19:49:02 -0400 Subject: [PATCH] Add parallel triangular expansion algorithm. Compile ok, still contain bugs, waiting for debug --- ...rallel_triangular_expansion_visibility_2.h | 589 ++++++++++++++++++ 1 file changed, 589 insertions(+) create mode 100644 Visibility_2/include/CGAL/Parallel_triangular_expansion_visibility_2.h diff --git a/Visibility_2/include/CGAL/Parallel_triangular_expansion_visibility_2.h b/Visibility_2/include/CGAL/Parallel_triangular_expansion_visibility_2.h new file mode 100644 index 00000000000..e779ed7668f --- /dev/null +++ b/Visibility_2/include/CGAL/Parallel_triangular_expansion_visibility_2.h @@ -0,0 +1,589 @@ +// Copyright (c) 2013 Technical University Braunschweig (Germany). +// All rights reserved. +// +// This file is part of CGAL (www.cgal.org). +// You can redistribute it and/or modify it under the terms of the GNU +// General Public License as published by the Free Software Foundation, +// either version 3 of the License, or (at your option) any later version. +// +// 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): Michael Hemmer +// Ning Xu +// + +#ifndef CGAL_PARALLEL_TRIANGULAR_EXPANSION_VISIBILITY_2__H +#define CGAL_PARALLEL_TRIANGULAR_EXPANSION_VISIBILITY_2__H + +#include +#include +#include +#include + +#ifdef CGAL_LINKED_WITH_TBB + #include "tbb/parallel_do.h" + #include "tbb/concurrent_vector.h" +#endif + +namespace CGAL { + +template < class Arrangement_2_, + class RegularizationTag = CGAL::Tag_true, + class ConcurrencyTag = CGAL::Parallel_tag > +class Parallel_triangular_expansion_visibility_2 +{ + typedef typename Arrangement_2_::Geometry_traits_2 Geometry_traits_2; + typedef typename Geometry_traits_2::Kernel K; +public: + // Currently only consider with same type for both + typedef Arrangement_2_ Arrangement_2; + typedef typename Arrangement_2::Traits_2 Traits_2; + typedef typename Arrangement_2::Halfedge_const_handle Halfedge_const_handle; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator + Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Face_const_handle Face_const_handle; + typedef typename Arrangement_2::Face_handle Face_handle; + typedef typename Arrangement_2::Vertex_const_handle Vertex_const_handle; + typedef typename Arrangement_2::Vertex_handle Vertex_handle; + + typedef typename K::Point_2 Point_2; + typedef typename Geometry_traits_2::Ray_2 Ray_2; + typedef typename Geometry_traits_2::Segment_2 Segment_2; + typedef typename Geometry_traits_2::Line_2 Line_2; + typedef typename Geometry_traits_2::Vector_2 Vector_2; + typedef typename Geometry_traits_2::Direction_2 Direction_2; + typedef typename Geometry_traits_2::FT Number_type; + typedef typename Geometry_traits_2::Object_2 Object_2; + + // TODO + typedef RegularizationTag Regularization_tag; + + typedef CGAL::Tag_true Supports_general_polygon_tag; + typedef CGAL::Tag_true Supports_simple_polygon_tag; + +private: + typedef CGAL::Triangulation_vertex_base_2 Vb; + typedef CGAL::Constrained_triangulation_face_base_2 Fb; + typedef CGAL::Triangulation_data_structure_2 TDS; + typedef CGAL::No_intersection_tag Itag; + typedef CGAL::Constrained_Delaunay_triangulation_2 CDT; + + typedef typename CDT::Face_handle CDT_Face_handle; + typedef std::pair CDT_Edge; + typedef typename CDT::Vertex_handle CDT_Vertex_handle; + +private: + class Edge + { + private: + CDT_Face_handle fh; + int idx; + CDT_Vertex_handle lvh; + CDT_Vertex_handle rvh; + public: + Edge () + : fh( NULL ), idx( -1 ), lvh( NULL ), rvh( NULL ) + {} + Edge ( CDT_Face_handle f, int i, CDT_Vertex_handle l, CDT_Vertex_handle r ) + : fh( f ), idx( i ), lvh( l ), rvh( r ) + {} + CDT_Face_handle face () const + { return fh; } + int index () const + { return idx; } + CDT_Vertex_handle left_source () const + { return lvh; } + CDT_Vertex_handle right_source () const + { return rvh; } + CDT_Vertex_handle left_vertex () const + { return fh->vertex( fh->cw( idx ) ); } + CDT_Vertex_handle right_vertex () const + { return fh->vertex( fh->ccw( idx ) ); } + + CDT_Edge edge () const + { return CDT_Edge( fh, idx ); } + + void trace ( ostream& os, int level ) const + { + for ( int i = 0; i < level; i++ ) os << " "; + os << "lvh = " << lvh->point() << ", rvh = " << rvh->point() << " , left = " << left_vertex()->point() << " , right = " << right_vertex()->point() << endl; + } + }; + + class Edge_bundle + { + private: + std::queue data; + std::vector constrained; + public: + Edge_bundle () + {} + + int capacity () const + { return 128; } + int size () const + { return data.size(); } + bool empty () const + { return data.empty(); } + bool full () const + { return ( data.size() == capacity() ); } + int constrained_size () const + { return constrained.size(); } + const Edge& get_constrained( int pos ) const + { return constrained[pos]; } + + const Edge& front () const + { return data.front(); } + void push( const Edge& e ) + { data.push( e ); } + void pop () + { data.pop(); } + void clear () + { + while ( !data.empty() ) data.pop(); + constrained.clear(); + } + void push_constrained ( const Edge& e ) + { + constrained.push_back( e ); + } + + void trace ( ostream& os, int level ) const + { + std::queue tmp = data; + int i = 0; + while ( !tmp.empty() ) { + for ( int j = 0; j < level; j++ ) os << " "; + os << "edge[" << i << "] = "; + tmp.front().trace( os, level+1 ); + i++; + tmp.pop(); + } + } + }; + +private: + typedef std::vector EB_vector; +#ifdef CGAL_LINKED_WITH_TBB + typedef tbb::concurrent_vector Result_container; +#else + typedef std::vector Result_container; +#endif + +private: + +#ifdef CGAL_LINKED_WITH_TBB + class Parallel_expand + { + private: + Parallel_triangular_expansion_visibility_2 * algo; + Result_container * p_result; + public: + Parallel_expand( Parallel_triangular_expansion_visibility_2 * a, + Result_container * pr ) + : algo( a ), p_result( pr ) + {} + void operator () ( Edge_bundle& bundle, tbb::parallel_do_feeder& feeder ) const + { + algo->compute_visibility_bundle( bundle ); + for ( int i = 0; i < bundle.constrained_size(); i++ ) + p_result->push_back( bundle.get_constrained( i ) ); + if ( bundle.empty() ) + return; + Edge_bundle res; + while ( !bundle.empty() ) { + res.push( bundle.front() ); + bundle.pop(); + if ( res.full() ) { + feeder.add( res ); + res.clear(); + } + } + if ( !res.empty() ) + feeder.add( res ); + } + }; +#endif + +private: + Edge left_edge ( CDT_Face_handle f, int idx ) + { return Edge( f, f->cw(idx), f->vertex(f->ccw(idx)), f->vertex(idx) ); } + Edge right_edge ( CDT_Face_handle f, int idx ) + { return Edge( f, f->ccw(idx), f->vertex(idx), f->vertex(f->cw(idx)) ); } + + CDT_Face_handle find_face_having_halfedge( CDT_Face_handle fh, int index ) + { + CDT_Face_handle curr, circ, nfh; + int curr_idx, nindex; + + curr = circ = fh; + curr_idx = index; + do { + if ( curr->vertex( curr->cw( curr_idx ) )->point() == query_edge->source()->point() ) + return curr; + nfh = curr->neighbor( curr->cw(curr_idx) ); + nindex = nfh->index( curr ); + curr = nfh; + curr_idx = nfh->cw( nindex ); + } while ( curr != circ && !p_cdt->is_infinite( curr ) ); + + if ( p_cdt->is_infinite( curr ) ) { + curr = circ = fh; + curr_idx = index; + do { + if ( curr->vertex( curr->cw( curr_idx ) )->point() == query_edge->source()->point() ) + return curr; + nfh = curr->neighbor( curr->ccw(curr_idx) ); + nindex = nfh->index( curr ); + curr = nfh; + curr_idx = nfh->ccw( nindex ); + } while ( curr != circ && !p_cdt->is_infinite( curr ) ); + } + // Something must be wrong here! +cout << "must be wrong" << endl; + return NULL; + } + int find_index_of_halfedge( CDT_Face_handle fh ) + { + for ( int i = 0; i < 3; i++ ) { + if ( fh->vertex( i )->point() == query_pt ) + return i; + } + // Something must be wrong here! + return -1; + } + void insert_edges_of_vertex_query( Edge_bundle& edges, CDT_Face_handle fh, int index ) + { + CDT_Face_handle curr = fh; + int curr_idx = index; + do { + edges.push( left_edge( curr, curr->ccw( curr_idx ) ) ); + CDT_Face_handle nfh = curr->neighbor( curr->cw(curr_idx) ); + int nindex = nfh->index( curr ); + curr = nfh; + curr_idx = nfh->cw( nindex ); + } while ( curr->vertex(curr->cw(curr_idx))->point() != query_edge->next()->target()->point() ); + } + + int expand_threshold () + { return 1024; } + + void compute_visibility_bundle( Edge_bundle& bundle ) + { + int count = 0; + int threshold = expand_threshold(); + while ( !bundle.empty() && count < threshold ) { + const Edge& e = bundle.front(); + + if ( p_cdt->is_constrained( e.edge() ) ) { + // A constrained edge + bundle.push_constrained( e ); + } else { + // An unconstrained edge, should be expanded + // deal with the left edge + + CDT_Face_handle nfh = e.face()->neighbor( e.index() ); + int nindex = nfh->index( e.face() ); + int rindex = p_cdt->ccw( nindex ); + int lindex = p_cdt->cw( nindex ); + + CDT_Vertex_handle nvh = nfh->vertex( nindex ); + CDT_Vertex_handle rvh = nfh->vertex( lindex ); + CDT_Vertex_handle lvh = nfh->vertex( rindex ); + + CGAL::Orientation ro = CGAL::orientation( query_pt, + e.right_source()->point(), + nvh->point() ); + CGAL::Orientation lo = CGAL::orientation( query_pt, + e.left_source()->point(), + nvh->point() ); + + if ( ro == CGAL::LEFT_TURN ) { + // The right edge is seen + if ( lo == CGAL::LEFT_TURN ) { + // No split is needed + bundle.push( Edge( nfh, rindex, e.left_source(), e.right_source() ) ); + count++; + } else { + // split at the new vertex + bundle.push( Edge( nfh, rindex, nvh, e.right_source() ) ); + count++; + } + } + + if ( lo == CGAL::RIGHT_TURN ) { + // The left edge is seen + if ( ro == CGAL::RIGHT_TURN ) { + // No split is needed + bundle.push( Edge( nfh, lindex, e.left_source(), e.right_source() ) ); + count++; + } else { + // split at the new vertex + bundle.push( Edge( nfh, lindex, e.left_source(), nvh ) ); + count++; + } + } + } + + bundle.pop(); + } + } + + void compute_visibility_parallel ( Edge_bundle& bundle, + Result_container& result, + CGAL::Sequential_tag ) + { + } + + void compute_visibility_parallel ( Edge_bundle& bundle, + Result_container& result, + CGAL::Parallel_tag ) + { +#ifdef CGAL_LINKED_WITH_TBB + std::vector bundles; + bundles.push_back( bundle ); + Parallel_expand exp( this, &result ); + tbb::parallel_do( bundles.begin(), bundles.end(), exp ); +#else + compute_visibility_parallel( bundle, CGAL::Sequential_tag() ); +#endif + } + + + Point_2 ray_seg_intersection ( const Point_2& dp, + const Point_2& s, + const Point_2& t ) + { + Ray_2 ray( query_pt, dp ); + Segment_2 seg( s, t ); + assert( CGAL::do_intersect( ray, seg ) ); + CGAL::Object obj = CGAL::intersection( ray, seg ); + const Point_2 * ipoint = CGAL::object_cast( &obj ); + assert( ipoint != NULL ); + return *ipoint; + } + + template < class VARR > + void compute_visibility_impl ( VARR& arr_out) + { + arr_out.clear(); + // Locate the query point + typename CDT::Locate_type location; + int index; + CDT_Face_handle fh = p_cdt->locate( query_pt, location, index ); + + Edge_bundle edges; + + switch( location ) { + case CDT::FACE: + // The query point locates in a face + // It only occurs in FACE_QUERY + assert( query_type == FACE_QUERY ); + edges.push( Edge( fh, 0, fh->vertex(2), fh->vertex(1) ) ); + edges.push( Edge( fh, 1, fh->vertex(0), fh->vertex(2) ) ); + edges.push( Edge( fh, 2, fh->vertex(1), fh->vertex(0) ) ); + break; + case CDT::VERTEX: + // The query point locates at a vertex + // It only occurs in VERTEX_QUERY + // Find the face containing the halfedge in query + assert( query_type == VERTEX_QUERY ); + fh = find_face_having_halfedge( fh, index ); + assert( fh != NULL ); + index = find_index_of_halfedge( fh ); + assert( index >= 0 ); + + insert_edges_of_vertex_query( edges, fh, index ); + break; + case CDT::EDGE: + // The query point locates on an edge + // It may occur in FACE_QUERY and EDGE_QUERY + assert( query_type != VERTEX_QUERY ); + if ( !p_cdt->is_constrained( std::make_pair(fh, index) ) ) { + // the edge is not contrained + CDT_Face_handle nfh = fh->neighbor( index ); + int nindex = nfh->index( fh ); + assert( !p_cdt->is_infinite( fh ) ); + assert( !p_cdt->is_infinite( nfh ) ); + edges.push( left_edge( fh, index ) ); + edges.push( right_edge( fh, index ) ); + edges.push( left_edge( nfh, nindex ) ); + edges.push( right_edge( nfh, nindex ) ); + } else { + // the edge is contrained + if ( fh->vertex(fh->cw(index))->point() == query_edge->target()->point() ) { + // the locating edge is the halfedge in query + assert( !p_cdt->is_infinite( fh ) ); + edges.push( left_edge( fh, index ) ); + edges.push( right_edge( fh, index ) ); + } else { + // the locating edge is not the halfedge in query + CDT_Face_handle nfh = fh->neighbor( index ); + int nindex = nfh->index( fh ); + assert( !p_cdt->is_infinite( nfh ) ); + edges.push( left_edge( nfh, nindex ) ); + edges.push( right_edge( nfh, nindex ) ); + } + } + break; + } + +#ifndef NDEBUG + cout << "... Initialize edges ..." << endl; + edges.trace( cout, 0 ); + cout << endl; +#endif + + Result_container results; + + compute_visibility_parallel( edges, results, ConcurrencyTag() ); + +#ifndef NDEBUG + for ( int i = 0; i < results.size(); i++ ) { + results[i].trace( cout, 0 ); + } +#endif + + // Generate results + std::vector segs; + for ( int i = 0; i < results.size(); i++ ) { + // push segments on the edge + Point_2 lp = ray_seg_intersection( results[i].left_source()->point(), + results[i].left_vertex()->point(), + results[i].right_vertex()->point() ); + Point_2 rp = ray_seg_intersection( results[i].right_source()->point(), + results[i].left_vertex()->point(), + results[i].right_vertex()->point() ); + if ( lp != rp ) + segs.push_back( Segment_2( lp, rp ) ); + + //push segments on the left side + if ( results[i].left_source() != results[i].left_vertex() ) { + segs.push_back( Segment_2( results[i].left_source()->point(), lp ) ); + } + + //push segments on the right side + if ( results[i].right_source() != results[i].right_vertex() ) { + segs.push_back( Segment_2( results[i].right_source()->point(), rp ) ); + } + } + +#ifndef NDEBUG + cout << "Output segments ... " << endl; + for ( int i = 0; i < segs.size(); i++ ) { + cout << segs[i].source() << " --> " << segs[i].target() << endl; + } +#endif + + //CGAL::insert_non_intersecting_curves( arr_out, segs.begin(), segs.end() ); + CGAL::insert( arr_out, segs.begin(), segs.end() ); + } + + void init_cdt () + { + std::vector< std::pair > constraints; + typename Arrangement_2::Edge_const_iterator eit; + for ( eit = p_arr->edges_begin(); eit != p_arr->edges_end(); eit++ ) { + constraints.push_back( std::make_pair( eit->source()->point(), + eit->target()->point() ) ); + } + CDT * p = new CDT( constraints.begin(), constraints.end() ); + assert( p != NULL ); + p_cdt = boost::shared_ptr( p ); + } + +public: + Parallel_triangular_expansion_visibility_2 () + : p_arr( NULL ) + {} + Parallel_triangular_expansion_visibility_2 ( const Arrangement_2& arr ) + : p_arr( &arr ) + { init_cdt(); } + void attach ( const Arrangement_2& arr ) + { + if ( p_arr != &arr ) { + p_arr = arr; + init_cdt(); + } + } + void detach () + { + p_arr = NULL; + p_cdt = boost::shared_ptr(); + } + bool is_attached () const + { return ( p_arr != NULL ); } + const Arrangement_2& arr () const + { return &p_arr; } + const std::string name () const + { return std::string( "PT_visibility_2" ); } + + + template < class VARR > + typename VARR::Face_handle + compute_visibility ( const Point_2& q, + const Face_const_handle& fh, + VARR& arr_out ) + { + query_pt = q; + query_type = FACE_QUERY; + + + compute_visibility_impl( arr_out ); + + assert( arr_out.number_of_faces() == 2 ); + + if ( arr_out.faces_begin()->is_unbounded() ) + return ++arr_out.faces_begin(); + else + return arr_out.faces_begin(); + } + + template < class VARR > + typename VARR::Face_handle + compute_visibility ( const Point_2& q, + const Halfedge_const_handle& eh, + VARR& arr_out ) + { + if ( q == eh->source()->point() ) + return compute_visibility( q, eh->prev(), arr_out ); + + query_pt = q; + query_edge = eh; + if ( query_pt == eh->target()->point() ) { + query_type = VERTEX_QUERY; + } else { + query_type = EDGE_QUERY; + } + compute_visibility_impl( arr_out ); + + assert( arr_out.number_of_faces() == 2 ); + + if ( arr_out.faces_begin()->is_unbounded() ) + return ++arr_out.faces_begin(); + else + return arr_out.faces_begin(); + } + + +private: + const Arrangement_2* p_arr; + boost::shared_ptr p_cdt; + + Point_2 query_pt; + enum { VERTEX_QUERY, EDGE_QUERY, FACE_QUERY } query_type; + Halfedge_const_handle query_edge; +}; + +} // namespace CGAL + +#endif //CGAL_TRIANGULAR_EXPANSION_VISIBILITY_2__H