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 aaace81b7ba..80e3a9d10ff 100644 --- a/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h +++ b/Poisson_surface_reconstruction_3/include/CGAL/Poisson_reconstruction_function.h @@ -205,6 +205,67 @@ 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) { } + + // Copy operator to satisfy vector, shouldn't be used + Cached_bary_coord(const Cached_bary_coord&) + { + CGAL_error(); + } + + 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]; } + }; + + // 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) { } + + // 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 + { + 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. @@ -213,14 +274,13 @@ 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; // 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; @@ -284,7 +344,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))) @@ -304,7 +364,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))) { @@ -323,7 +383,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, @@ -527,21 +587,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 @@ -552,19 +614,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 @@ -580,11 +644,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 @@ -627,7 +688,13 @@ public: void initialize_matrix_entry(Cell_handle ch) const { - boost::array & entry = (*m_Bary)[ch->info()]; + Cached_bary_coord& bary = (*m_bary)[ch->info()]; + + if (bary.is_initialized()) + return; + + // If the cache was uninitialized, this thread is in charge of + // initializing it const Point& pa = ch->vertex(0)->point(); const Point& pb = ch->vertex(1)->point(); const Point& pc = ch->vertex(2)->point(); @@ -638,9 +705,13 @@ public: Vector vc = pc - pd; 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]); + 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_initialized(); } /// \endcond @@ -847,10 +918,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(); 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 b8cb454be73..efd62081391 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 @@ -29,6 +29,15 @@ if ( CGAL_FOUND ) # Executables that require Eigen 3.1 create_single_source_cgal_program( "poisson_reconstruction_test.cpp" ) target_link_libraries(poisson_reconstruction_test PUBLIC CGAL::Eigen_support) + + find_package(TBB) + include(CGAL_TBB_support) + if (TBB_FOUND) + create_single_source_cgal_program( "poisson_and_parallel_mesh_3.cpp" ) + target_link_libraries(poisson_and_parallel_mesh_3 PUBLIC CGAL::Eigen_support CGAL::TBB_support) + 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.") 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..c053152f5f5 --- /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; +}