diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp index 05ad64b0d8b..e6cab966bfb 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -9,6 +9,8 @@ #include +#include + using Kernel = CGAL::Simple_cartesian; using FT = typename Kernel::FT; using Vector = typename Kernel::Vector_3; @@ -20,7 +22,7 @@ using Polygon_range = std::vector >; int main(int, char**) { const CGAL::Bbox_3 bbox { -1.0, -1.0, -1.0, 1.0, 1.0, 1.0 }; - const FT spacing = 0.04; + const FT spacing = 0.002; const Vector vec_spacing(spacing, spacing, spacing); // Euclidean distance function to the origin @@ -36,9 +38,14 @@ int main(int, char**) Point_range points; Polygon_range polygons; + const tbb::tick_count start = tbb::tick_count::now(); + // execute marching cubes with an isovalue of 0.8 CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, polygons); + const tbb::tick_count end = tbb::tick_count::now(); + std::cout << (end - start).seconds() << std::endl; + // save ouput indexed mesh to a file, in the OFF format CGAL::IO::write_OFF("result.off", points, polygons); diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp index 5acb6f59f7c..7015db3a914 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp @@ -17,7 +17,7 @@ using Polygon_range = std::vector >; int main(int, char**) { - const std::string fname = CGAL::data_file_path("images/skull_2.9.inr"); + const std::string fname = "../examples/Isosurfacing_3/FullHead.inr";//CGAL::data_file_path("images/skull_2.9.inr"); // load volumetric image from a file CGAL::Image_3 image; @@ -30,6 +30,11 @@ int main(int, char**) // convert image to a Cartesian grid Grid grid{image}; + for (std::size_t i = 0; i < grid.xdim(); i++) + for (std::size_t j = 0; j < grid.ydim(); j++) + for (std::size_t k = 0; k < grid.zdim(); k++) + grid.value(i, j, k) = 2 * 1120 - grid.value(i, j, k); + // create a domain from the grid auto domain = CGAL::Isosurfacing::create_explicit_Cartesian_grid_domain(grid); @@ -38,7 +43,7 @@ int main(int, char**) Polygon_range polygons; // execute marching cubes - CGAL::Isosurfacing::marching_cubes(domain, 2.9 /*isovalue*/, points, polygons); + CGAL::Isosurfacing::marching_cubes(domain, 1120 /*isovalue*/, points, polygons); // save output indexed mesh to a file, in the OFF format CGAL::IO::write_OFF("result.off", points, polygons); diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h index 92770bddb6a..ed1887886e6 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h @@ -20,6 +20,7 @@ #ifdef CGAL_LINKED_WITH_TBB #include +#include #endif // CGAL_LINKED_WITH_TBB #include @@ -213,14 +214,22 @@ public: const std::size_t sk = size_k; // for now only parallelize outer loop - auto iterator = [&f, sj, sk](const tbb::blocked_range& r) { - for(std::size_t i = r.begin(); i != r.end(); ++i) - for(std::size_t j=0; j& r) { + const std::size_t i_begin = r.pages().begin(); + const std::size_t i_end = r.pages().end(); + const std::size_t j_begin = r.rows().begin(); + const std::size_t j_end = r.rows().end(); + const std::size_t k_begin = r.cols().begin(); + const std::size_t k_end = r.cols().end(); + + for(std::size_t i = i_begin; i != i_end; ++i) + for(std::size_t j = j_begin; j != j_end; ++j) + for(std::size_t k = k_begin; k != k_end; ++k) f({i, j, k}); }; - tbb::parallel_for(tbb::blocked_range(0, size_i - 1), iterator); + tbb::blocked_range3d range(0, size_i - 1, 0, size_j - 1, 0, size_k - 1); + tbb::parallel_for(range, iterator); } #endif // CGAL_LINKED_WITH_TBB diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h index 386c68d1e39..ea5e15da42f 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h @@ -48,7 +48,7 @@ #include #ifdef CGAL_LINKED_WITH_TBB -#include +#include #else #include #endif @@ -187,7 +187,12 @@ void mc_construct_triangles(const int i_case, const int eg2 = Cube_table::triangle_cases[t_index + 2]; // insert new triangle in list - triangles.push_back({vertices[eg0], vertices[eg1], vertices[eg2]}); + #ifdef CGAL_LINKED_WITH_TBB + auto& tris = triangles.local(); + #else + auto& tris = triangles; + #endif + tris.push_back({vertices[eg0], vertices[eg1], vertices[eg2]}); } } @@ -198,17 +203,28 @@ void triangles_to_polygon_soup(const TriangleRange& triangles, PointRange& points, PolygonRange& polygons) { - for(auto& triangle : triangles) - { - const std::size_t id = points.size(); + #ifdef CGAL_LINKED_WITH_TBB + for(const auto& triangle_list : triangles) + { + #else + const auto& triangle_list = triangles; + #endif - points.push_back(triangle[0]); - points.push_back(triangle[1]); - points.push_back(triangle[2]); + for(const auto& triangle : triangle_list) + { + const std::size_t id = points.size(); - // simply use increasing indices - polygons.push_back({id + 2, id + 1, id + 0}); - } + points.push_back(triangle[0]); + points.push_back(triangle[1]); + points.push_back(triangle[2]); + + // simply use increasing indices + polygons.push_back({id + 2, id + 1, id + 0}); + } + + #ifdef CGAL_LINKED_WITH_TBB + } + #endif } // Marching Cubes implemented as a functor that runs on every cell of the grid @@ -225,7 +241,7 @@ public: using Cell_descriptor = typename Domain::Cell_descriptor; #ifdef CGAL_LINKED_WITH_TBB - using Triangles = tbb::concurrent_vector >; + using Triangles = tbb::enumerable_thread_specific>>; #else using Triangles = std::vector >; #endif @@ -245,7 +261,7 @@ public: { } // gets the created triangle list - const Triangles& triangles() const + Triangles& triangles() { return m_triangles; } diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h index 492ef74eb16..68ea824741d 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h @@ -46,11 +46,16 @@ #include #include +#include +#include + #include #include #include #include #include +#include +#include namespace CGAL { namespace Isosurfacing { @@ -73,26 +78,27 @@ private: using Edge_descriptor = typename Domain::Edge_descriptor; using Cell_descriptor = typename Domain::Cell_descriptor; - using uint = unsigned int; + using Point_index = std::size_t; + using Edge_index = std::size_t; + + using Edge_point_map = tbb::concurrent_hash_map; private: const Domain& m_domain; FT m_isovalue; - Point_range& m_points; - Polygon_range& m_polygons; + std::atomic m_point_counter; + tbb::concurrent_vector m_points; - std::mutex mutex; + Edge_point_map m_edges; + + tbb::concurrent_vector> m_triangles; public: TMC_functor(const Domain& domain, - const FT isovalue, - Point_range& points, - Polygon_range& polygons) + const FT isovalue) : m_domain(domain), - m_isovalue(isovalue), - m_points(points), - m_polygons(polygons) + m_isovalue(isovalue) { } void operator()(const Cell_descriptor& cell) @@ -101,18 +107,19 @@ public: std::array corners; const int i_case = get_cell_corners(m_domain, cell, m_isovalue, corners, values); - const int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1 - if(Cube_table::intersected_edges[i_case] == 0 || - Cube_table::intersected_edges[i_case] == all_bits_set) - { - return; - } - // this is the only difference to mc - int tcm = int(Cube_table::t_ambig[i_case]); + const int tcm = Cube_table::t_ambig[i_case]; if(tcm == 105) { - p_slice(cell, m_isovalue, values, corners, i_case); + if (p_slice(cell, m_isovalue, values, corners, i_case)) + return; + else + std::cerr << "WARNING: the result might not be topologically correct" << std::endl; + } + + constexpr int all_bits_set = (1 << (8 + 1)) - 1; // last 8 bits are 1 + if(i_case == 0 || i_case == all_bits_set) + { return; } @@ -122,7 +129,6 @@ public: // @todo improve triangle generation // construct triangles - std::lock_guard lock(mutex); for(int t=0; t<16; t += 3) { const int t_index = i_case * 16 + t; @@ -135,38 +141,63 @@ public: const int eg1 = Cube_table::triangle_cases[t_index + 1]; const int eg2 = Cube_table::triangle_cases[t_index + 2]; - const std::size_t p0_idx = m_points.size(); - - m_points.push_back(vertices[eg0]); - m_points.push_back(vertices[eg1]); - m_points.push_back(vertices[eg2]); // insert new triangle into list - m_polygons.emplace_back(); - auto& triangle = m_polygons.back(); - - triangle.push_back(p0_idx + 2); - triangle.push_back(p0_idx + 1); - triangle.push_back(p0_idx + 0); } } private: - void add_triangle(const std::size_t p0, - const std::size_t p1, - const std::size_t p2) + Edge_index compute_edge_index(const Cell_descriptor& cell, int edge) { - std::lock_guard lock(mutex); + // edge is in 0 - 11 - m_polygons.emplace_back(); - auto& triangle = m_polygons.back(); + // there are 12 edges, assign to each vertex three edges, the global edge numbering + // consists of 3*global_vertex_id + edge_offset. + const unsigned long long gei_pattern_ = 670526590282893600ull; - triangle.push_back(p0); - triangle.push_back(p1); - triangle.push_back(p2); + // the edge global index is given by the vertex global index + the edge offset + const std::size_t shift = 5 * edge; + const std::size_t ix = cell[0] + ((gei_pattern_ >> shift) & 1); // global_edge_id[edge][0]; + const std::size_t iy = cell[1] + ((gei_pattern_ >> (shift + 1)) & 1); // global_edge_id[edge][1]; + const std::size_t iz = cell[2] + ((gei_pattern_ >> (shift + 2)) & 1); // global_edge_id[edge][2]; + const std::size_t off_val = ((gei_pattern_ >> (shift + 3)) & 3); + + int g_edg = int(m_cell_shift_factor * m_ugrid.global_index(ix, iy, iz) + off_val); } - void p_slice(const Cell_descriptor& cell, + bool find_point(const Edge_index e, Point_index& i) + { + Edge_point_map::const_accessor acc; + if (m_edges.find(acc, e)) + { + i = acc->second; + return true; + } + return false; + } + + Point_index add_point(const Point_3& p, const Edge_index e) + { + Edge_point_map::accessor acc; + if (!m_edges.insert(acc, e)) + return acc->second; + + const Point_index i = m_point_counter++; + acc->second = i; + acc.release(); + + m_points.grow_to_at_least(i); + m_points[i] = p; + } + + void add_triangle(const Point_index p0, + const Point_index p1, + const Point_index p2) + { + m_triangles.push_back({p0, p1, p2}); + } + + bool p_slice(const Cell_descriptor& cell, const FT i0, const std::array& values, const std::array& corners, @@ -177,9 +208,7 @@ private: typename Geom_traits::Compute_z_3 z_coord = m_domain.geom_traits().compute_z_3_object(); typename Geom_traits::Construct_point_3 point = m_domain.geom_traits().construct_point_3_object(); - // there are 12 edges, assign to each vertex three edges, the global edge numbering - // consists of 3*global_vertex_id + edge_offset. - const unsigned long long gei_pattern_ = 670526590282893600ull; + using uint = unsigned int; // code edge end vertices for each of the 12 edges const unsigned char l_edges_[12] = {16, 49, 50, 32, 84, 117, 118, 100, 64, 81, 115, 98}; @@ -198,18 +227,11 @@ private: // collect vertices unsigned short flag{1}; - for(int eg=0; eg<12; ++eg) + for(int eg = 0; eg < 12; ++eg) { if(flag & Cube_table::intersected_edges[i_case]) { - // the edge global index is given by the vertex global index + the edge offset - // uint shift = 5 * eg; - // const int ix = i_index + (int)((gei_pattern_ >> shift) & 1); // global_edge_id[eg][0]; - // const int iy = j_index + (int)((gei_pattern_ >> (shift + 1)) & 1); // global_edge_id[eg][1]; - // const int iz = k_index + (int)((gei_pattern_ >> (shift + 2)) & 1); // global_edge_id[eg][2]; - // const int off_val = (int)((gei_pattern_ >> (shift + 3)) & 3); - - // int g_edg = int(m_cell_shift_factor * m_ugrid.global_index(ix, iy, iz) + off_val); + // generate vertex here, do not care at this point if vertex already exists uint v0, v1; @@ -409,8 +431,8 @@ private: } else { - std::cerr << "ERROR: can't correctly triangulate cell's face\n"; - return; + // std::cerr << "ERROR: can't correctly triangulate cell's face\n"; + return false; } } } @@ -481,8 +503,8 @@ private: } else { - std::cerr << "ERROR: can't correctly triangulate cell's face\n"; - return; + // std::cerr << "ERROR: can't correctly triangulate cell's face\n"; + return false; } } } @@ -1138,6 +1160,7 @@ private: m_points.emplace_back(px, py, pz); } } + return true; } }; diff --git a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h index 2d0a93356aa..b86c9c39fc0 100644 --- a/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h +++ b/Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h @@ -74,7 +74,7 @@ void marching_cubes(const Domain& domain, { // run topologically correct marching cubes // and directly write the result to points and triangles - internal::TMC_functor functor(domain, isovalue, points, triangles); + internal::TMC_functor functor(domain, isovalue); domain.template iterate_cells(functor); } else