songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Start rework and testing

This commit is contained in:
Julian Stahl 2023-08-09 00:55:59 +02:00
parent 6cacbee79f
commit 499e0a8a2c
6 changed files with 139 additions and 79 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp
View File

@ -9,6 +9,8 @@
#include <vector>
#include <tbb/tick_count.h>
using Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
using Vector = typename Kernel::Vector_3;
@ -20,7 +22,7 @@ using Polygon_range = std::vector<std::vector<std::size_t> >;
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);

9
Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp
View File

@ -17,7 +17,7 @@ using Polygon_range = std::vector<std::vector<std::size_t> >;
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);

19
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/Grid_topology_3.h
View File

@ -20,6 +20,7 @@
#ifdef CGAL_LINKED_WITH_TBB
#include <tbb/parallel_for.h>
#include <tbb/blocked_range3d.h>
#endif // CGAL_LINKED_WITH_TBB
#include <array>
@ -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<std::size_t>& r) {
for(std::size_t i = r.begin(); i != r.end(); ++i)
for(std::size_t j=0; j<sj - 1; ++j)
for(std::size_t k=0; k<sk - 1; ++k)
auto iterator = [&f, sj, sk](const tbb::blocked_range3d<std::size_t>& 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<std::size_t>(0, size_i - 1), iterator);
tbb::blocked_range3d<std::size_t> range(0, size_i - 1, 0, size_j - 1, 0, size_k - 1);
tbb::parallel_for(range, iterator);
}
#endif // CGAL_LINKED_WITH_TBB

42
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/marching_cubes_functors.h
View File

@ -48,7 +48,7 @@
#include <CGAL/assertions.h>
#ifdef CGAL_LINKED_WITH_TBB
#include <tbb/concurrent_vector.h>
#include <tbb/enumerable_thread_specific.h>
#else
#include <vector>
#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<std::array<Point_3, 3> >;
using Triangles = tbb::enumerable_thread_specific<std::vector<std::array<Point_3, 3>>>;
#else
using Triangles = std::vector<std::array<Point_3, 3> >;
#endif
@ -245,7 +261,7 @@ public:
{ }
// gets the created triangle list
const Triangles& triangles() const
Triangles& triangles()
{
return m_triangles;
}

137
Isosurfacing_3/include/CGAL/Isosurfacing_3/internal/topologically_correct_marching_cubes_functors.h
View File

@ -46,11 +46,16 @@
#include <CGAL/Isosurfacing_3/internal/marching_cubes_functors.h>
#include <CGAL/Isosurfacing_3/internal/tables.h>
#include <tbb/concurrent_vector.h>
#include <tbb/concurrent_hash_map.h>
#include <array>
#include <atomic>
#include <cmath>
#include <map>
#include <mutex>
#include <atomic>
#include <optional>
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<Edge_index, Point_index>;
private:
const Domain& m_domain;
FT m_isovalue;
Point_range& m_points;
Polygon_range& m_polygons;
std::atomic<Point_index> m_point_counter;
tbb::concurrent_vector<Point_3> m_points;
std::mutex mutex;
Edge_point_map m_edges;
tbb::concurrent_vector<std::array<Point_index, 3>> 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<Point_3, 8> 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<std::mutex> 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<std::mutex> 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<FT, 8>& values,
const std::array<Point_3, 8>& 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;
}
};

2
Isosurfacing_3/include/CGAL/Isosurfacing_3/marching_cubes_3.h
View File

@ -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<Domain, PointRange, TriangleRange> functor(domain, isovalue, points, triangles);
internal::TMC_functor<Domain, PointRange, TriangleRange> functor(domain, isovalue);
domain.template iterate_cells<ConcurrencyTag>(functor);
}
else