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 b814c056603..f875dabfbbf 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -5,6 +5,8 @@ #include #include +#include + using Kernel = CGAL::Simple_cartesian; using FT = typename Kernel::FT; using Point = typename Kernel::Point_3; @@ -32,10 +34,16 @@ int main(int, char**) Point_range points; Triangle_range triangles; + const tbb::tick_count start = tbb::tick_count::now(); + // execute marching cubes with a given isovalue const FT isovalue = 0.8; CGAL::Isosurfacing::marching_cubes(domain, isovalue, points, triangles); + 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("output.off", points, triangles); 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..030e7c8d383 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,45 @@ 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::array; + + struct Hash_compare + { + static size_t hash(const Edge_index& key) + { + std::size_t res = 17; + res = res * 31 + std::hash()(key[0]); + res = res * 31 + std::hash()(key[1]); + res = res * 31 + std::hash()(key[2]); + res = res * 31 + std::hash()(key[3]); + return res; + } + + static bool equal(const Edge_index& key1, const Edge_index& key2) + { + return key1[0] == key2[0] && key1[1] == key2[1] && key1[2] == key2[2] && key1[3] == key2[3]; + } + }; + + 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 +125,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 +147,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 +159,90 @@ 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(); + const Point_index p0 = add_point(vertices[eg0], compute_edge_index(cell, eg0)); + const Point_index p1 = add_point(vertices[eg1], compute_edge_index(cell, eg1)); + const Point_index p2 = add_point(vertices[eg2], compute_edge_index(cell, eg2)); - triangle.push_back(p0_idx + 2); - triangle.push_back(p0_idx + 1); - triangle.push_back(p0_idx + 0); + add_triangle(p2, p1, p0); + } + } + + // gets the created triangle list + template + void to_triangle_soup(PointRange& points, TriangleRange& triangles) const + { + points.insert(points.begin(), m_points.begin(), m_points.end()); + for (const std::array& tri : m_triangles) { + triangles.push_back({ tri[0], tri[1], tri[2] }); } } 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); + + return { 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 + 1); + m_points[i] = p; + + return i; + } + + Point_index add_point_unchecked(const Point_3& p) + { + const Point_index i = m_point_counter++; + + m_points.grow_to_at_least(i + 1); + m_points[i] = p; + + return i; + } + + 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 +253,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}; @@ -191,25 +265,18 @@ private: // A hexahedron has twelve edges, save the intersection of the isosurface with the edge // save global edge and global vertex index of isosurface - std::vector vertices(12); + std::vector vertices(12); // save local coordinate along the edge of intersection point std::vector ecoord{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; // 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; @@ -223,22 +290,9 @@ private: const FT py = (1 - l) * y_coord(corners[v0]) + l * y_coord(corners[v1]); const FT pz = (1 - l) * z_coord(corners[v0]) + l * z_coord(corners[v1]); - // set vertex in map - // set vertex index - // auto s_index = m_vertices.find(vertices[eg].g_edg); - // if(s_index == m_vertices.end()) - // { - const int g_idx = static_cast(m_points.size()); - vertices[eg] = g_idx; - // m_vertices[vertices[eg].g_edg] = g_idx; - m_points.push_back(point(px, py, pz)); - //} else { - // vertices[eg] = s_index->second; - //} + // add vertex and insert to map + vertices[eg] = add_point(point(px, py, pz), compute_edge_index(cell, eg)); } - /*else { - e_set[eg] = false; - }*/ // next edge flag <<= 1; @@ -409,8 +463,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 +535,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; } } } @@ -815,7 +869,7 @@ private: // compute vertices of inner hexagon, save new vertices in list and compute and keep // global vertices index to build triangle connectivity later on. - std::size_t tg_idx[6]; + Point_index tg_idx[6]; for(int i=0; i<6; ++i) { const FT u = hvt[i][0]; @@ -834,8 +888,7 @@ private: w * ((1 - v) * (z_coord(corners[4]) + u * (z_coord(corners[5]) - z_coord(corners[4]))) + v * (z_coord(corners[6]) + u * (z_coord(corners[7]) - z_coord(corners[6])))); - tg_idx[i] = m_points.size(); - m_points.push_back(point(px, py, pz)); + tg_idx[i] = add_point_unchecked(point(px, py, pz)); } // triangulate contours with inner hexagon @@ -1113,10 +1166,10 @@ private: wcoord * ((1 - vcoord) * (z_coord(corners[4]) + ucoord * (z_coord(corners[5]) - z_coord(corners[4]))) + vcoord * (z_coord(corners[6]) + ucoord * (z_coord(corners[7]) - z_coord(corners[6])))); - const std::size_t g_index = m_points.size(); + bool pt_used = false; + Point_index g_index = 0; // loop over the contours - bool pt_used = false; for(int i=0; i functor(domain, isovalue, points, triangles); + internal::TMC_functor functor(domain, isovalue); domain.template iterate_cells(functor); + functor.to_triangle_soup(points, triangles); } else {