From 07b6cb75c798210afc87b9b2d7c88ba9bb6af58c Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 29 Sep 2020 14:02:34 +0200 Subject: [PATCH 1/5] thread safe cached barycentric coordinates --- .../CGAL/Poisson_reconstruction_function.h | 97 ++++++++++++++----- 1 file changed, 71 insertions(+), 26 deletions(-) diff --git a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h index 10db057f40b..6b71b63d517 100644 --- a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h @@ -216,6 +216,32 @@ private: typedef typename Triangulation::All_cells_iterator All_cells_iterator; typedef typename Triangulation::Locate_type Locate_type; + enum Cache_state { UNINITIALIZED, BUSY, INITIALIZED }; + // Thread-safe cache for barycentric coordinates of a cell + class Cached_bary_coord + { + private: + std::atomic m_state; + std::array m_bary; + public: + Cached_bary_coord() : m_state (UNINITIALIZED) { } + Cached_bary_coord(const Cached_bary_coord& other) + : m_state (other.m_state.load()) // not atomic + { } + + Cached_bary_coord& operator= (const Cached_bary_coord& other) + { + m_state.store (other.m_state.load()); // not atomic + } + + Cache_state exchange_state(const Cache_state& s) { return m_state.exchange(s); } + Cache_state get_state() { return m_state; } + void set_state(const Cache_state& s) { m_state = s; } + + const double& operator[] (const std::size_t& idx) const { return m_bary[idx]; } + double& operator[] (const std::size_t& idx) { return m_bary[idx]; } + }; + // Data members. // Warning: the Surface Mesh Generation package makes copies of implicit functions, // thus this class must be lightweight and stateless. @@ -224,8 +250,7 @@ private: // operator() is pre-computed on vertices of *m_tr by solving // the Poisson equation Laplacian(f) = divergent(normals field). boost::shared_ptr m_tr; - - mutable boost::shared_ptr > > m_Bary; + mutable std::shared_ptr > m_bary; mutable std::vector Dual; mutable std::vector Normal; @@ -295,7 +320,7 @@ public: PointPMap point_pmap, ///< property map: `value_type of InputIterator` -> `Point` (the position of an input point). NormalPMap normal_pmap ///< property map: `value_type of InputIterator` -> `Vector` (the *oriented* normal of an input point). ) - : m_tr(new Triangulation), m_Bary(new std::vector > ) + : m_tr(new Triangulation), m_bary(new std::vector) , average_spacing(CGAL::compute_average_spacing (CGAL::make_range(first, beyond), 6, CGAL::parameters::point_map(point_pmap))) @@ -315,7 +340,7 @@ public: PointPMap point_pmap, ///< property map: `value_type of InputIterator` -> `Point` (the position of an input point). NormalPMap normal_pmap, ///< property map: `value_type of InputIterator` -> `Vector` (the *oriented* normal of an input point). Visitor visitor) - : m_tr(new Triangulation), m_Bary(new std::vector > ) + : m_tr(new Triangulation), m_bary(new std::vector) , average_spacing(CGAL::compute_average_spacing(CGAL::make_range(first, beyond), 6, CGAL::parameters::point_map(point_pmap))) { @@ -334,7 +359,7 @@ public: boost::is_convertible::value_type, Point> >::type* = 0 ) - : m_tr(new Triangulation), m_Bary(new std::vector > ) + : m_tr(new Triangulation), m_bary(new std::vector) , average_spacing(CGAL::compute_average_spacing(CGAL::make_range(first, beyond), 6)) { forward_constructor(first, beyond, @@ -590,11 +615,8 @@ public: void initialize_barycenters() const { - m_Bary->resize(m_tr->number_of_cells()); - - for(std::size_t i=0; i< m_Bary->size();i++){ - (*m_Bary)[i][0]=-1; - } + m_bary->clear(); + m_bary->resize(m_tr->number_of_cells()); } void initialize_cell_normals() const @@ -637,20 +659,44 @@ public: void initialize_matrix_entry(Cell_handle ch) const { - boost::array & entry = (*m_Bary)[ch->info()]; - const Point& pa = ch->vertex(0)->point(); - const Point& pb = ch->vertex(1)->point(); - const Point& pc = ch->vertex(2)->point(); - const Point& pd = ch->vertex(3)->point(); + Cached_bary_coord& bary = (*m_bary)[ch->info()]; + Cache_state s = bary.exchange_state(BUSY); - Vector va = pa - pd; - Vector vb = pb - pd; - Vector vc = pc - pd; + // If the cache was already initialized, then let's just restore + // the value + if (s == INITIALIZED) + bary.set_state (INITIALIZED); - internal::invert(va.x(), va.y(), va.z(), - vb.x(), vb.y(), vb.z(), - vc.x(), vc.y(), vc.z(), - entry[0],entry[1],entry[2],entry[3],entry[4],entry[5],entry[6],entry[7],entry[8]); + // If the cache was uninitialized, this thread is in charge of + // initializing it + else if (s == UNINITIALIZED) + { + const Point& pa = ch->vertex(0)->point(); + const Point& pb = ch->vertex(1)->point(); + const Point& pc = ch->vertex(2)->point(); + const Point& pd = ch->vertex(3)->point(); + + Vector va = pa - pd; + Vector vb = pb - pd; + Vector vc = pc - pd; + + internal::invert(va.x(), va.y(), va.z(), + vb.x(), vb.y(), vb.z(), + vc.x(), vc.y(), vc.z(), + bary[0],bary[1],bary[2], + bary[3],bary[4],bary[5], + bary[6],bary[7],bary[8]); + + bary.set_state (INITIALIZED); + } + + // Else, the cache is busy, let's just wait for it to be + // initialized + else + { + CGAL_assertion (s == BUSY); + while (bary.get_state() != INITIALIZED) { } + } } /// \endcond @@ -857,10 +903,9 @@ private: // vb.x(), vb.y(), vb.z(), // vc.x(), vc.y(), vc.z(), // i00, i01, i02, i10, i11, i12, i20, i21, i22); - const boost::array & i = (*m_Bary)[cell->info()]; - if(i[0]==-1){ - initialize_matrix_entry(cell); - } + initialize_matrix_entry(cell); + const Cached_bary_coord& i = (*m_bary)[cell->info()]; + // UsedBary[cell->info()] = true; a = i[0] * vp.x() + i[3] * vp.y() + i[6] * vp.z(); b = i[1] * vp.x() + i[4] * vp.y() + i[7] * vp.z(); From c005b37e2dae01852cc81a76f79e047959b6756a Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 29 Sep 2020 14:03:22 +0200 Subject: [PATCH 2/5] Thread safe cell hint --- .../CGAL/Poisson_reconstruction_function.h | 70 +++++++++++++------ 1 file changed, 49 insertions(+), 21 deletions(-) diff --git a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h index 6b71b63d517..3e7e3d38e9f 100644 --- a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h @@ -242,6 +242,30 @@ private: double& operator[] (const std::size_t& idx) { return m_bary[idx]; } }; + // Wrapper for thread safety of maintained cell hint for fast + // locate, with conversions atomic/Cell_handle + class Cell_hint + { + std::atomic m_cell; + public: + + Cell_hint() : m_cell(nullptr) { } + Cell_hint(const Cell_hint& other) + : m_cell (other.m_cell.load()) // not atomic + { } + + Cell_hint& operator= (const Cell_hint& other) + { + m_cell.store (other.m_cell.load()); // not atomic + } + + Cell_handle get() const + { + return Triangulation_data_structure::Cell_range::s_iterator_to(*m_cell); + } + void set (Cell_handle ch) { m_cell = ch.operator->(); } + }; + // Data members. // Warning: the Surface Mesh Generation package makes copies of implicit functions, // thus this class must be lightweight and stateless. @@ -256,7 +280,7 @@ private: // contouring and meshing Point m_sink; // Point with the minimum value of operator() - mutable Cell_handle m_hint; // last cell found = hint for next search + mutable Cell_hint m_hint; // last cell found = hint for next search FT average_spacing; @@ -562,21 +586,23 @@ public: boost::tuple special_func(const Point& p) const { - m_hint = m_tr->locate(p ,m_hint ); // no hint when we use hierarchy + Cell_handle hint = m_hint.get(); + hint = m_tr->locate(p, hint); // no hint when we use hierarchy + m_hint.set(hint); - if(m_tr->is_infinite(m_hint)) { - int i = m_hint->index(m_tr->infinite_vertex()); - return boost::make_tuple(m_hint->vertex((i+1)&3)->f(), - m_hint, true); + if(m_tr->is_infinite(hint)) { + int i = hint->index(m_tr->infinite_vertex()); + return boost::make_tuple(hint->vertex((i+1)&3)->f(), + hint, true); } FT a,b,c,d; - barycentric_coordinates(p,m_hint,a,b,c,d); - return boost::make_tuple(a * m_hint->vertex(0)->f() + - b * m_hint->vertex(1)->f() + - c * m_hint->vertex(2)->f() + - d * m_hint->vertex(3)->f(), - m_hint, false); + barycentric_coordinates(p,hint,a,b,c,d); + return boost::make_tuple(a * hint->vertex(0)->f() + + b * hint->vertex(1)->f() + + c * hint->vertex(2)->f() + + d * hint->vertex(3)->f(), + hint, false); } /// \endcond @@ -587,19 +613,21 @@ public: */ FT operator()(const Point& p) const { - m_hint = m_tr->locate(p ,m_hint); + Cell_handle hint = m_hint.get(); + hint = m_tr->locate(p, hint); + m_hint.set(hint); - if(m_tr->is_infinite(m_hint)) { - int i = m_hint->index(m_tr->infinite_vertex()); - return m_hint->vertex((i+1)&3)->f(); + if(m_tr->is_infinite(hint)) { + int i = hint->index(m_tr->infinite_vertex()); + return hint->vertex((i+1)&3)->f(); } FT a,b,c,d; - barycentric_coordinates(p,m_hint,a,b,c,d); - return a * m_hint->vertex(0)->f() + - b * m_hint->vertex(1)->f() + - c * m_hint->vertex(2)->f() + - d * m_hint->vertex(3)->f(); + barycentric_coordinates(p,hint,a,b,c,d); + return a * hint->vertex(0)->f() + + b * hint->vertex(1)->f() + + c * hint->vertex(2)->f() + + d * hint->vertex(3)->f(); } /// \cond SKIP_IN_MANUAL From 2e56893cf0d1830bd88a41ccd1746d06076c11a6 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 29 Sep 2020 14:33:20 +0200 Subject: [PATCH 3/5] Add test for poisson and parallel Mesh_3 --- .../CMakeLists.txt | 8 +- .../poisson_and_parallel_mesh_3.cpp | 78 +++++++++++++++++++ 2 files changed, 85 insertions(+), 1 deletion(-) create mode 100644 Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp diff --git a/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/CMakeLists.txt b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/CMakeLists.txt index e1d9326840b..f3e847dce16 100644 --- a/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/CMakeLists.txt +++ b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/CMakeLists.txt @@ -30,6 +30,13 @@ if ( CGAL_FOUND ) # Executables that require Eigen 3.1 create_single_source_cgal_program( "poisson_reconstruction_test.cpp" ) + find_package(TBB) + if (TBB_FOUND) + create_single_source_cgal_program( "poisson_and_parallel_mesh_3.cpp" ) + CGAL_target_use_TBB(poisson_and_parallel_mesh_3) + else() + message(STATUS "NOTICE: test with parallel Mesh_3 needs TBB and will not be compiled.") + endif() else() message(STATUS "NOTICE: Some of the executables in this directory need Eigen 3.1 (or greater) and will not be compiled.") @@ -41,4 +48,3 @@ else() message(STATUS "NOTICE: This program requires the CGAL library, and will not be compiled.") endif() - diff --git a/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp new file mode 100644 index 00000000000..9f9619201a1 --- /dev/null +++ b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp @@ -0,0 +1,78 @@ +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel; +typedef Kernel::FT FT; +typedef Kernel::Point_3 Point_3; +typedef Kernel::Vector_3 Vector_3; + +typedef std::pair Pwn; +typedef CGAL::First_of_pair_property_map Point_map; +typedef CGAL::Second_of_pair_property_map Vector_map; + +typedef CGAL::Poisson_reconstruction_function Poisson; + +typedef CGAL::Labeled_mesh_domain_3 Implicit_domain; +typedef CGAL::Mesh_triangulation_3::type Mesh_tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +int main(int, char**) +{ + std::vector points; + + std::ifstream stream("data/oni.pwn"); + if (!stream || + !CGAL::read_xyz_points + (stream, std::back_inserter(points), + CGAL::parameters:: + point_map(Point_map()). + normal_map(Vector_map()))) + { + std::cerr << "Error: cannot read file" << std::endl; + return EXIT_FAILURE; + } + + Poisson poisson (points.begin(), points.end(), Point_map(), Vector_map()); + if (!poisson.compute_implicit_function()) + { + std::cerr << "Error: cannot compute implicit function" << std::endl; + return EXIT_FAILURE; + } + + CGAL::Bbox_3 bbox + = CGAL::bbox_3 (boost::make_transform_iterator + (points.begin(), + CGAL::Property_map_to_unary_function()), + boost::make_transform_iterator + (points.end(), + CGAL::Property_map_to_unary_function())); + + Implicit_domain domain + = Implicit_domain::create_implicit_mesh_domain + (poisson, bbox); + + Mesh_criteria criteria (CGAL::parameters::facet_angle = 30, + CGAL::parameters::facet_size = 4, + CGAL::parameters::facet_distance = 0.1); + C3t3 c3t3 = CGAL::make_mesh_3 (domain, criteria, + CGAL::parameters::manifold()); + + return EXIT_SUCCESS; +} From 3e0fde968320788784bdadddcc1c9f0725bcc760 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Tue, 29 Sep 2020 15:56:13 +0200 Subject: [PATCH 4/5] Improve thread-safety structures from review --- .../CGAL/Poisson_reconstruction_function.h | 87 +++++++++---------- 1 file changed, 42 insertions(+), 45 deletions(-) diff --git a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h index 3e7e3d38e9f..ade184edc37 100644 --- a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h @@ -225,18 +225,33 @@ private: std::array m_bary; public: Cached_bary_coord() : m_state (UNINITIALIZED) { } - Cached_bary_coord(const Cached_bary_coord& other) - : m_state (other.m_state.load()) // not atomic - { } - Cached_bary_coord& operator= (const Cached_bary_coord& other) + // Copy operator to satisfy vector, shouldn't be used + Cached_bary_coord(const Cached_bary_coord&) { - m_state.store (other.m_state.load()); // not atomic + CGAL_error(); } - Cache_state exchange_state(const Cache_state& s) { return m_state.exchange(s); } - Cache_state get_state() { return m_state; } - void set_state(const Cache_state& s) { m_state = s; } + bool is_initialized() + { + Cache_state s = m_state; + if (s == UNINITIALIZED) + { + // If the following line successfully replaces UNINITIALIZED + // by BUSY, then the current thread in charge of initialization + if (m_state.compare_exchange_weak(s, BUSY)) + return false; + } + // Otherwise, either the thread is BUSY by another thread, or + // it's already INITIALIZED. Either way, we way until it's INITIALIZED + else + while (m_state != INITIALIZED) { } + + // At this point, it's always INITIALIZED + return true; + } + + void set_initialized() { m_state = INITIALIZED; } const double& operator[] (const std::size_t& idx) const { return m_bary[idx]; } double& operator[] (const std::size_t& idx) { return m_bary[idx]; } @@ -250,14 +265,10 @@ private: public: Cell_hint() : m_cell(nullptr) { } - Cell_hint(const Cell_hint& other) - : m_cell (other.m_cell.load()) // not atomic - { } - Cell_hint& operator= (const Cell_hint& other) - { - m_cell.store (other.m_cell.load()); // not atomic - } + // Poisson_reconstruction_function should be copyable, although we + // don't need to copy that + Cell_hint(const Cell_hint&) : m_cell(nullptr) { } Cell_handle get() const { @@ -688,43 +699,29 @@ public: void initialize_matrix_entry(Cell_handle ch) const { Cached_bary_coord& bary = (*m_bary)[ch->info()]; - Cache_state s = bary.exchange_state(BUSY); - // If the cache was already initialized, then let's just restore - // the value - if (s == INITIALIZED) - bary.set_state (INITIALIZED); + if (bary.is_initialized()) + return; // If the cache was uninitialized, this thread is in charge of // initializing it - else if (s == UNINITIALIZED) - { - const Point& pa = ch->vertex(0)->point(); - const Point& pb = ch->vertex(1)->point(); - const Point& pc = ch->vertex(2)->point(); - const Point& pd = ch->vertex(3)->point(); + const Point& pa = ch->vertex(0)->point(); + const Point& pb = ch->vertex(1)->point(); + const Point& pc = ch->vertex(2)->point(); + const Point& pd = ch->vertex(3)->point(); - Vector va = pa - pd; - Vector vb = pb - pd; - Vector vc = pc - pd; + Vector va = pa - pd; + Vector vb = pb - pd; + Vector vc = pc - pd; - internal::invert(va.x(), va.y(), va.z(), - vb.x(), vb.y(), vb.z(), - vc.x(), vc.y(), vc.z(), - bary[0],bary[1],bary[2], - bary[3],bary[4],bary[5], - bary[6],bary[7],bary[8]); + internal::invert(va.x(), va.y(), va.z(), + vb.x(), vb.y(), vb.z(), + vc.x(), vc.y(), vc.z(), + bary[0],bary[1],bary[2], + bary[3],bary[4],bary[5], + bary[6],bary[7],bary[8]); - bary.set_state (INITIALIZED); - } - - // Else, the cache is busy, let's just wait for it to be - // initialized - else - { - CGAL_assertion (s == BUSY); - while (bary.get_state() != INITIALIZED) { } - } + bary.set_initialized(); } /// \endcond From 37d2002b2c3e2c6a3ccc5ff29ec3698cef9cdb11 Mon Sep 17 00:00:00 2001 From: Simon Giraudot Date: Wed, 30 Sep 2020 15:22:07 +0200 Subject: [PATCH 5/5] Fix tabs --- .../poisson_and_parallel_mesh_3.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp index 9f9619201a1..c053152f5f5 100644 --- a/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp +++ b/Poisson_surface_reconstruction_3/test/Poisson_surface_reconstruction_3/poisson_and_parallel_mesh_3.cpp @@ -49,14 +49,14 @@ int main(int, char**) return EXIT_FAILURE; } - Poisson poisson (points.begin(), points.end(), Point_map(), Vector_map()); - if (!poisson.compute_implicit_function()) + Poisson poisson (points.begin(), points.end(), Point_map(), Vector_map()); + if (!poisson.compute_implicit_function()) { std::cerr << "Error: cannot compute implicit function" << std::endl; return EXIT_FAILURE; } - CGAL::Bbox_3 bbox + CGAL::Bbox_3 bbox = CGAL::bbox_3 (boost::make_transform_iterator (points.begin(), CGAL::Property_map_to_unary_function()), @@ -64,14 +64,14 @@ int main(int, char**) (points.end(), CGAL::Property_map_to_unary_function())); - Implicit_domain domain + Implicit_domain domain = Implicit_domain::create_implicit_mesh_domain (poisson, bbox); - Mesh_criteria criteria (CGAL::parameters::facet_angle = 30, + Mesh_criteria criteria (CGAL::parameters::facet_angle = 30, CGAL::parameters::facet_size = 4, CGAL::parameters::facet_distance = 0.1); - C3t3 c3t3 = CGAL::make_mesh_3 (domain, criteria, + C3t3 c3t3 = CGAL::make_mesh_3 (domain, criteria, CGAL::parameters::manifold()); return EXIT_SUCCESS;