The remove is now fully parallel.

The remove operations are done in parallel as long as the dimension stays = 3, then the few remaining points are removed sequentially.
This commit is contained in:
Clement Jamin 2013-03-27 13:47:40 +01:00
parent 57d84ecb24
commit de1f43e30b
3 changed files with 373 additions and 59 deletions

View File

@ -1,6 +1,8 @@
#ifndef TYPEDEFS_H
#define TYPEDEFS_H
#define CGAL_CONCURRENT_TRIANGULATION_3_PROFILING
#include <vector> //dynamic array
#include <list> //linked list
@ -84,17 +86,21 @@ private:
* and point location is then O(n^(1/3)) time
*/
#ifdef CONCURRENT_TRIANGULATION_3
typedef CGAL::Spatial_grid_lock_data_structure_3<
CGAL::Tag_priority_blocking_with_atomics> Lock_ds;
typedef CGAL::Triangulation_data_structure_3<
Vertex_base<Kernel>,
CGAL::Triangulation_ds_cell_base_3<>,
CGAL::Compact_container_strategy_base,
CGAL::Compact_container_strategy_base,
CGAL::Parallel_tag > Tds;
typedef CGAL::Delaunay_triangulation_3<Kernel, Tds> DT3;
typedef CGAL::Delaunay_triangulation_3<
Kernel, Tds, CGAL::Default, Lock_ds> DT3;
#else
typedef CGAL::Triangulation_data_structure_3< Vertex_base<Kernel> > Tds;
typedef CGAL::Delaunay_triangulation_3<
Kernel, Tds, CGAL::Fast_location> DT3;
Kernel, Tds/*, CGAL::Fast_location*/> DT3;
#endif
typedef DT3::Object Object_3;

View File

@ -50,6 +50,7 @@
#ifdef CGAL_LINKED_WITH_TBB
# include <tbb/parallel_for.h>
# include <tbb/enumerable_thread_specific.h>
# include <tbb/concurrent_vector.h>
#endif
#ifdef CGAL_DELAUNAY_3_OLD_REMOVE
@ -281,9 +282,10 @@ public:
{
size_t num_points = points.size();
int i = 0;
// Sequential until dim = 3 (or more)
// Insert "num_points_seq" points sequentially
// (or more if dim < 3 after that)
Vertex_handle hint;
size_t num_points_seq = (std::min)(num_points, (size_t)500); // CJTODO: 100, 1000... ?
size_t num_points_seq = (std::min)(num_points, (size_t)500);
while (dimension() < 3 || i < num_points_seq)
{
hint = insert(points[i], hint);
@ -650,7 +652,10 @@ public:
}
// REMOVE
void remove(Vertex_handle v, bool *p_could_lock_zone = 0);
void remove(Vertex_handle v);
// Concurrency-safe
// See Triangulation_3::remove for more information
bool remove(Vertex_handle v, bool *p_could_lock_zone);
// return new cells (internal)
template <class OutputItCells>
@ -670,8 +675,9 @@ public:
if (is_parallel())
{
std::vector<Vertex_handle> vertices(first, beyond);
tbb::concurrent_vector<Vertex_handle> vertices_to_remove_sequentially;
tbb::task_scheduler_init init(10); // CJTODO TEMP
tbb::task_scheduler_init init(2); // CJTODO TEMP
// CJTODO: lambda functions OK?
tbb::parallel_for(
tbb::blocked_range<size_t>( 0, vertices.size()),
@ -679,21 +685,29 @@ public:
{
for( size_t i_vertex = r.begin() ; i_vertex != r.end() ; ++i_vertex)
{
bool could_lock_zone;
Vertex_handle v = vertices[i_vertex];
bool could_lock_zone, needs_to_be_done_sequentially;
do
{
remove(vertices[i_vertex], &could_lock_zone);
needs_to_be_done_sequentially =
remove(v, &could_lock_zone);
this->unlock_all_elements();
/*if (!could_lock_zone)
{
if (r.end() - i_vertex > 2)
std::swap(vertices[i_vertex], vertices[i_vertex + 1 + (std::rand()%(r.end() - i_vertex - 2))]);
else if (i_vertex != r.end()-1)
std::swap(vertices[i_vertex], vertices[i_vertex + 1]);
}*/
} while (!could_lock_zone);
if (needs_to_be_done_sequentially)
vertices_to_remove_sequentially.push_back(v);
}
});
// Do the rest sequentially
for ( tbb::concurrent_vector<Vertex_handle>::const_iterator
it = vertices_to_remove_sequentially.begin(),
it_end = vertices_to_remove_sequentially.end()
; it != it_end
; ++it)
{
remove(*it);
}
}
// Sequential
else
@ -1122,13 +1136,26 @@ public:
template < class Gt, class Tds, class Lds >
void
Delaunay_triangulation_3<Gt,Tds,Default,Lds>::
remove(Vertex_handle v)
{
Self tmp;
Vertex_remover<Self> remover (tmp);
Tr_Base::remove(v, remover);
CGAL_triangulation_expensive_postcondition(is_valid());
}
template < class Gt, class Tds, class Lds >
bool
Delaunay_triangulation_3<Gt,Tds,Default,Lds>::
remove(Vertex_handle v, bool *p_could_lock_zone)
{
Self tmp;
Vertex_remover<Self> remover (tmp);
Tr_Base::remove(v, remover, p_could_lock_zone);
bool ret = Tr_Base::remove(v, remover, p_could_lock_zone);
CGAL_triangulation_expensive_postcondition(is_valid());
return ret;
}
template < class Gt, class Tds, class Lds >

View File

@ -630,7 +630,7 @@ public:
// Create the 3D triangulation of p0, p1, p3 and p4
// Precondition: p0, p1, p3 and p4 MUST BE positively oriented
Triangulation_3(const Point &p0, const Point &p1, // CJTODO: add "const"
Triangulation_3(const Point &p0, const Point &p1,
const Point &p3, const Point &p4,
const GT & gt = GT(), Lock_data_structure *p_lock_ds = 0)
: Base(p_lock_ds), _gt(gt)
@ -1360,12 +1360,21 @@ protected:
bool test_dim_down_using_incident_cells_3(
Vertex_handle v, std::vector<Cell_handle> &incident_cells,
std::vector<Vertex_handle> &adj_vertices) const;
std::vector<Vertex_handle> &adj_vertices,
bool *p_could_lock_zone = 0) const;
// REMOVAL
template < class VertexRemover >
void remove(Vertex_handle v, VertexRemover &remover,
bool *p_could_lock_zone = 0);
void remove(Vertex_handle v, VertexRemover &remover);
template < class VertexRemover >
// Concurrency-safe version
// Pre-condition: dimension = 3
// The return value is only meaningful if *p_could_lock_zone = true:
// * returns true if the vertex was removed
// * returns false if the vertex wasn't removed since it would decrease
// the dimension => needs to be done sequentially
bool remove(Vertex_handle v, VertexRemover &remover,
bool *p_could_lock_zone);
template < class VertexRemover, class OutputItCells >
void remove_and_give_new_cells(Vertex_handle v, VertexRemover &remover,
@ -1513,6 +1522,10 @@ private:
template < class VertexRemover >
VertexRemover& remove_2D(Vertex_handle v, VertexRemover &remover);
template < class VertexRemover >
VertexRemover& remove_3D(Vertex_handle v, VertexRemover &remover);
// Version of remove_3D if the incident cells and the adjacent vertices
// are already known
template < class VertexRemover >
VertexRemover& remove_3D(Vertex_handle v, VertexRemover &remover,
const std::vector<Cell_handle> &inc_cells,
std::vector<Vertex_handle> &adj_vertices);
@ -1847,7 +1860,70 @@ public:
}
return true;
}
template <class OutputIterator>
bool
try_lock_and_get_adjacent_vertices_and_cells_3(
Vertex_handle v, OutputIterator vertices,
std::vector<Cell_handle> &cells) const
{
// We need to lock v individually first, to be sure v->cell() is valid
if (!try_lock_vertex(v))
return false;
Cell_handle d = v->cell();
if (!try_lock_cell(d)) // LOCK
{
return false;
}
cells.push_back(d);
d->tds_data().mark_in_conflict();
int head=0;
int tail=1;
do {
Cell_handle c = cells[head];
for (int i=0; i<4; ++i) {
if (c->vertex(i) == v)
continue;
Cell_handle next = c->neighbor(i);
if (!try_lock_cell(next)) // LOCK
{
BOOST_FOREACH(Cell_handle& ch,
std::make_pair(cells.begin(), cells.end()))
{
ch->tds_data().clear();
}
cells.clear();
return false;
}
if (! next->tds_data().is_clear())
continue;
cells.push_back(next);
++tail;
next->tds_data().mark_in_conflict();
}
++head;
} while(head != tail);
std::set<Vertex_handle> tmp_vertices;
BOOST_FOREACH(Cell_handle& ch, std::make_pair(cells.begin(), cells.end()))
{
ch->tds_data().clear();
for (int i = 0; i < 4; ++i)
{
Vertex_handle w = ch->vertex(i);
if (w != v && tmp_vertices.insert(w).second)
{
*vertices = w;
}
}
}
return true;
}
template <class OutputIterator>
OutputIterator
@ -1889,7 +1965,6 @@ public:
return _tds.adjacent_vertices(v, vertices);
}
// correct name
template <class OutputIterator>
OutputIterator
adjacent_vertices_and_cells_3(Vertex_handle v, OutputIterator vertices,
@ -3963,15 +4038,26 @@ bool
Triangulation_3<GT,Tds,Lds>::
test_dim_down_using_incident_cells_3(
Vertex_handle v, std::vector<Cell_handle> &incident_cells,
std::vector<Vertex_handle> &adj_vertices) const
std::vector<Vertex_handle> &adj_vertices,
bool *p_could_lock_zone) const
{
CGAL_triangulation_precondition(dimension() == 3);
CGAL_triangulation_precondition(! is_infinite(v) );
// collect all vertices on the boundary
// Collect all vertices on the boundary
// and all incident cells
adjacent_vertices_and_cells_3(
v, std::back_inserter(adj_vertices), incident_cells);
if (p_could_lock_zone)
{
*p_could_lock_zone = try_lock_and_get_adjacent_vertices_and_cells_3(
v, std::back_inserter(adj_vertices), incident_cells);
if (*p_could_lock_zone == false)
return false;
}
else
{
adjacent_vertices_and_cells_3(
v, std::back_inserter(adj_vertices), incident_cells);
}
typedef Filter_iterator<
std::vector<Vertex_handle>::const_iterator,
@ -3981,12 +4067,10 @@ test_dim_down_using_incident_cells_3(
adj_vertices.end(), Infinite_tester(this), adj_vertices.begin());
Finite_vertex_iterator vit_end(
adj_vertices.end(), Infinite_tester(this));
//std::vector<Vertex_handle>::const_iterator vit = adj_vertices.begin();
const Point &p1 = (*vit++)->point();
const Point &p2 = (*vit++)->point();
const Point &p3 = (*vit++)->point();
//for (; vit != adj_vertices.end(); ++vit )
for ( ; vit != vit_end ; ++vit )
{
if (!coplanar(p1, p2, p3, (*vit)->point()))
@ -4466,6 +4550,188 @@ remove_2D(Vertex_handle v, VertexRemover &remover)
return remover;
}
template <class Gt, class Tds, class Lds>
template < class VertexRemover >
VertexRemover&
Triangulation_3<Gt,Tds,Lds>::
remove_3D(Vertex_handle v, VertexRemover &remover)
{
std::vector<Cell_handle> hole;
hole.reserve(64);
// Construct the set of vertex triples on the boundary
// with the facet just behind
typedef std::map<Vertex_triple,Facet> Vertex_triple_Facet_map;
Vertex_triple_Facet_map outer_map;
Vertex_triple_Facet_map inner_map;
make_hole_3D(v, outer_map, hole);
CGAL_assertion(remover.hidden_points_begin() ==
remover.hidden_points_end() );
// Output the hidden points.
for (typename std::vector<Cell_handle>::iterator
hi = hole.begin(), hend = hole.end(); hi != hend; ++hi)
remover.add_hidden_points(*hi);
bool inf = false;
// collect all vertices on the boundary
std::vector<Vertex_handle> vertices;
vertices.reserve(64);
adjacent_vertices(v, std::back_inserter(vertices));
// create a Delaunay triangulation of the points on the boundary
// and make a map from the vertices in remover.tmp towards the vertices
// in *this
unsigned int i = 0;
Unique_hash_map<Vertex_handle,Vertex_handle> vmap;
Cell_handle ch = Cell_handle();
#ifdef CGAL_TRIANGULATION_3_USE_THE_4_POINTS_CONSTRUCTOR
size_t num_vertices = vertices.size();
if (num_vertices >= 5)
{
//std::random_shuffle(vertices.begin(), vertices.end());
for (int j = 0 ; j < 4 ; ++j)
{
if (is_infinite(vertices[j]))
{
std::swap(vertices[j], vertices[4]);
break;
}
}
Orientation o = orientation(
vertices[0]->point(),
vertices[1]->point(),
vertices[2]->point(),
vertices[3]->point());
if (o == NEGATIVE)
std::swap(vertices[0], vertices[1]);
if (o != ZERO)
{
Vertex_handle vh1, vh2, vh3, vh4;
remover.tmp.init_tds(
vertices[0]->point(), vertices[1]->point(),
vertices[2]->point(), vertices[3]->point(),
vh1, vh2, vh3, vh4);
ch = vh1->cell();
vmap[vh1] = vertices[0];
vmap[vh2] = vertices[1];
vmap[vh3] = vertices[2];
vmap[vh4] = vertices[3];
i = 4;
}
}
#endif
for(; i < vertices.size(); i++){
if(! is_infinite(vertices[i])){
Vertex_handle vh = remover.tmp.insert(vertices[i]->point(), ch);
ch = vh->cell();
vmap[vh] = vertices[i];
}else {
inf = true;
}
}
if(remover.tmp.dimension()==2){
Vertex_handle fake_inf = remover.tmp.insert(v->point());
vmap[fake_inf] = infinite_vertex();
} else {
vmap[remover.tmp.infinite_vertex()] = infinite_vertex();
}
CGAL_triangulation_assertion(remover.tmp.dimension() == 3);
// Construct the set of vertex triples of remover.tmp
// We reorient the vertex triple so that it matches those from outer_map
// Also note that we use the vertices of *this, not of remover.tmp
if(inf){
for(All_cells_iterator it = remover.tmp.all_cells_begin(),
end = remover.tmp.all_cells_end(); it != end; ++it){
for(i=0; i < 4; i++){
Facet f = std::pair<Cell_handle,int>(it,i);
Vertex_triple vt_aux = make_vertex_triple(f);
Vertex_triple vt(vmap[vt_aux.first],vmap[vt_aux.third],vmap[vt_aux.second]);
make_canonical(vt);
inner_map[vt]= f;
}
}
} else {
for(Finite_cells_iterator it = remover.tmp.finite_cells_begin(),
end = remover.tmp.finite_cells_end(); it != end; ++it){
for(i=0; i < 4; i++){
Facet f = std::pair<Cell_handle,int>(it,i);
Vertex_triple vt_aux = make_vertex_triple(f);
Vertex_triple vt(vmap[vt_aux.first],vmap[vt_aux.third],vmap[vt_aux.second]);
make_canonical(vt);
inner_map[vt]= f;
}
}
}
// Grow inside the hole, by extending the surface
while(! outer_map.empty()){
typename Vertex_triple_Facet_map::iterator oit = outer_map.begin();
while(is_infinite(oit->first.first) ||
is_infinite(oit->first.second) ||
is_infinite(oit->first.third)){
++oit;
// otherwise the lookup in the inner_map fails
// because the infinite vertices are different
}
typename Vertex_triple_Facet_map::value_type o_vt_f_pair = *oit;
Cell_handle o_ch = o_vt_f_pair.second.first;
unsigned int o_i = o_vt_f_pair.second.second;
typename Vertex_triple_Facet_map::iterator iit =
inner_map.find(o_vt_f_pair.first);
CGAL_triangulation_assertion(iit != inner_map.end());
typename Vertex_triple_Facet_map::value_type i_vt_f_pair = *iit;
Cell_handle i_ch = i_vt_f_pair.second.first;
unsigned int i_i = i_vt_f_pair.second.second;
// create a new cell and glue it to the outer surface
Cell_handle new_ch = tds().create_cell();
new_ch->set_vertices(vmap[i_ch->vertex(0)], vmap[i_ch->vertex(1)],
vmap[i_ch->vertex(2)], vmap[i_ch->vertex(3)]);
o_ch->set_neighbor(o_i,new_ch);
new_ch->set_neighbor(i_i, o_ch);
// for the other faces check, if they can also be glued
for(i = 0; i < 4; i++){
if(i != i_i){
Facet f = std::pair<Cell_handle,int>(new_ch,i);
Vertex_triple vt = make_vertex_triple(f);
make_canonical(vt);
std::swap(vt.second,vt.third);
typename Vertex_triple_Facet_map::iterator oit2 = outer_map.find(vt);
if(oit2 == outer_map.end()){
std::swap(vt.second,vt.third);
outer_map[vt]= f;
} else {
// glue the faces
typename Vertex_triple_Facet_map::value_type o_vt_f_pair2 = *oit2;
Cell_handle o_ch2 = o_vt_f_pair2.second.first;
int o_i2 = o_vt_f_pair2.second.second;
o_ch2->set_neighbor(o_i2,new_ch);
new_ch->set_neighbor(i, o_ch2);
outer_map.erase(oit2);
}
}
}
outer_map.erase(oit);
}
tds().delete_vertex(v);
tds().delete_cells(hole.begin(), hole.end());
return remover;
}
template <class Gt, class Tds, class Lds>
template < class VertexRemover >
VertexRemover&
@ -4647,57 +4913,70 @@ template <class Gt, class Tds, class Lds>
template < class VertexRemover >
void
Triangulation_3<Gt, Tds, Lds>::
remove(Vertex_handle v, VertexRemover &remover,
bool *p_could_lock_zone) {
remove(Vertex_handle v, VertexRemover &remover) {
CGAL_triangulation_precondition( v != Vertex_handle());
CGAL_triangulation_precondition( !is_infinite(v));
CGAL_triangulation_expensive_precondition( tds().is_vertex(v) );
if (test_dim_down (v)) {
remove_dim_down (v, remover);
}
else {
switch (dimension()) {
case 1: remove_1D (v, remover); break;
case 2: remove_2D (v, remover); break;
case 3: remove_3D (v, remover); break;
default:
CGAL_triangulation_assertion (false);
}
}
}
template <class Gt, class Tds, class Lds>
template < class VertexRemover >
bool
Triangulation_3<Gt, Tds, Lds>::
remove(Vertex_handle v, VertexRemover &remover, bool *p_could_lock_zone)
{
// N.B.: dimension doesn't need to be atomic since the parallel removal
// will never decrease the dimension (the last few removes are done
// sequentially)
CGAL_triangulation_precondition( v != Vertex_handle());
CGAL_triangulation_precondition( !is_infinite(v));
CGAL_triangulation_precondition( dimension() == 3);
CGAL_triangulation_expensive_precondition( tds().is_vertex(v) );
#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING
static Profile_branch_counter_3 bcounter(
"early withdrawals / late withdrawals / successes [Delaunay_tri_3::remove]");
#endif
if (p_could_lock_zone && !try_lock_vertex(v))
bool needs_to_be_done_sequentially = false;
// Locking vertex v is a good start
if (!try_lock_vertex(v))
{
*p_could_lock_zone = false;
#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING
bcounter.increment_branch_2(); // THIS is an early withdrawal!
#endif
return;
}
switch (dimension())
else
{
case 1:
if (test_dim_down (v))
remove_dim_down (v, remover);
else
remove_1D (v, remover);
break;
case 2:
if (test_dim_down (v))
remove_dim_down (v, remover);
else
remove_2D (v, remover);
break;
case 3:
std::vector<Cell_handle> incident_cells;
incident_cells.reserve(64);
std::vector<Vertex_handle> adj_vertices;
adj_vertices.reserve(64);
bool dim_down = test_dim_down_using_incident_cells_3(
v, incident_cells, adj_vertices, p_could_lock_zone);
if (*p_could_lock_zone)
{
std::vector<Cell_handle> incident_cells;
incident_cells.reserve(64);
std::vector<Vertex_handle> adj_vertices;
adj_vertices.reserve(64);
if (test_dim_down_using_incident_cells_3(v, incident_cells, adj_vertices))
remove_dim_down (v, remover);
if (dim_down)
needs_to_be_done_sequentially = true;
else
remove_3D (v, remover, incident_cells, adj_vertices);
break;
remove_3D (v, remover, incident_cells, adj_vertices);
}
default:
if (test_dim_down (v))
remove_dim_down (v, remover);
else
CGAL_triangulation_assertion (false);
}
#ifdef CGAL_CONCURRENT_TRIANGULATION_3_PROFILING
@ -4709,6 +4988,8 @@ remove(Vertex_handle v, VertexRemover &remover,
bcounter.increment_branch_1(); // THIS is a late withdrawal!
}
#endif
return needs_to_be_done_sequentially;
}
// The remove here uses the remover, but