diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt b/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt index c6fdbeba7a0..4e3e7a7af1d 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt @@ -17,34 +17,38 @@ Different isovalues of a bunny. \section secmyintroduction Introduction -This package provides methods to compute a surface mesh representing an isosurface of a 3-dimensional scalar field. -An isosurface is defined as the surface on which the value of this field is equal to a given constant, i.e. the isovalue. -Depending on the isosurfacing method and the input data structure, the result is either a triangular, quadrilateral, or higher order polygonal indexed face set. - +This package provides algorithms to compute a surface mesh approximating an isosurface of a 3-dimensional scalar field defined over an input 3D domain. +An isosurface is defined as the surface on which the value of this scalar field equates a given constant, i.e. a user-defined isovalue. +Depending on the isosurfacing method and the input data structure, the ouput is either a triangular, quadrilateral, +or higher order polygonal surface mesh represented as an indexed face set. The output can be empty when the isovalue is absent in the input domain. \section secmyalgorithms Algorithms -There are multiple algorithms to extract isosurfaces. -This package contains Marching Cubes, topologically correct Marching Cubes, and Dual Contouring. +There is a variety of algorithms to extract isosurfaces. +This package offers Marching Cubes, topologically correct Marching Cubes and Dual Contouring. \subsection subsecmc Marching Cubes (MC) -MC processes all cells of the input domain individually. +MC runs over a 3D grid, i.e. a 3D domain partitioned into hexahedral cells. +It processes all cells of the input domain individually. Each cell corner gets a sign (+/-) to indicate if it is above or below the isovalue. -A new vertex is created on every cell edge where the sign changes, i.e. the isosurface is intersected. -The vertex position is computed via linear interpolation of the scalar values of the incident corners. -Depending on the configuration of signs at the corners the resulting vertices are connected to form triangles within the cell. +A vertex is created on every cell edge where the sign changes, i.e. where the edge intersects the isosurface. +The vertex position is computed via linear interpolation of the scalar values evaluated at the cell corners forming the edge. +Depending on the configuration of signs at the cell corners, the resulting vertices are connected to form triangles within the cell. -MC can process any input data structure that consists of hexahedral cells. -In case of a conforming grid, MC produces a triangle mesh that is manifold in most scenarios. -If the mesh is manifold and the isosurface does not intersect the domain boundaries, the mesh is also watertight. -Compared to other approaches the algorithm often generates more and many thin triangles with acute angles. -MC does not preserve sharp features of the input data. +The proposed implementation is generic in that it can process any grid-based data structure that consists of hexahedral cells. +In case of a conforming grid, MC generates as output a surface triangle mesh that is 2-manifold in most scenarios. [TODO: reformulate, as +it is either strict or not - precise in which cases it is not, and say whether you talk about just combinatorially 2-manifold, or truly 2-manifold] +If the mesh is manifold and the isosurface does not intersect the domain boundaries, then the output mesh is also watertight. + +Compared to other approaches such a Delaunay refinement (TODO: add link to component), the MC algorithm often generates more triangles, +and triangles with small angles. +MC does not preserve the sharp features present in the isovalue of the input scalar field. \cgalFigureAnchor{isosurfacing_mc_cases}
\cgalFigureCaptionBegin{isosurfacing_mc_cases} -Some example cases of Marching Cubes. +Ouputs of Marching Cubes over a variety of input domains. \cgalFigureCaptionEnd \subsection subsectmc Topologically correct Marching Cubes (TMC) diff --git a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt b/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt index 20d8c107719..80bfd5a25a7 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt @@ -8,7 +8,7 @@ \cgalPkgPicture{isosurfacing3_ico.png} \cgalPkgSummaryBegin \cgalPkgAuthor{Julian Stahl and Daniel Zint} -\cgalPkgDesc{This package implements isosurfacing algorithms to generate a mesh from a scalar field. The algorithms provide different guarantees for the output and should be chosen depending on the use-case. The representations of input data is flexible and includes explicit and implicit formats. } +\cgalPkgDesc{This package implements isosurfacing algorithms to generate a surface mesh from a 3D scalar field. The algorithms provide different guarantees for the output and should be chosen depending on each use case. The possible format for the input data is flexible and includes explicit as well as implicit representations. } \cgalPkgManuals{Chapter_Isosurfacing3,PkgIsosurfacing3Ref} \cgalPkgSummaryEnd \cgalPkgShortInfoBegin @@ -20,9 +20,9 @@ \cgalPkgDescriptionEnd This package provides algorithms to extract isosurfaces from different inputs. -The input is represented as a domain and can be an implicit function, a cartesion grid, or an octree. +The input is represented as a 3D domain and can be an implicit function, a cartesion grid or an octree. The output is an indexed face set that stores an isosurface in the form of a surface mesh. -Available algorithms include Marching Cubes, topologically correct Marching Cubes, and Dual Contouring. +Available algorithms include Marching Cubes, topologically correct Marching Cubes and Dual Contouring. \cgalClassifedRefPages diff --git a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp index 9066ec6366a..6456a5af55f 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp @@ -15,13 +15,14 @@ typedef CGAL::Cartesian_grid_3 Grid; typedef std::vector Point_range; typedef std::vector> Polygon_range; +// return 1.0 if value has positive sign, and -1.0 otherwise FT sign(FT value) { - return (value > 0) - (value < 0); + return (value > 0.0) - (value < 0.0); } int main() { - // create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3 - const CGAL::Bbox_3 bbox(-1, -1, -1, 1, 1, 1); + // create a cartesian grid with 7^3 grid points and the bounding box [-1, 1]^3 + const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0); std::shared_ptr grid = std::make_shared(7, 7, 7, bbox); // calculate the value at all grid points @@ -39,19 +40,20 @@ int main() { } } + // compute function gradient auto cube_gradient = [](const Point& p) { // the normal depends on the side of the cube const FT max_value = std::max({std::abs(p.x()), std::abs(p.y()), std::abs(p.z())}); - Vector g(0, 0, 0); + Vector g(0.0, 0.0, 0.0); if (max_value == std::abs(p.x())) { - g += Vector(sign(p.x()), 0, 0); + g += Vector(sign(p.x()), 0.0, 0.0); } if (max_value == std::abs(p.y())) { - g += Vector(0, sign(p.y()), 0); + g += Vector(0.0, sign(p.y()), 0.0); } if (max_value == std::abs(p.z())) { - g += Vector(0, 0, sign(p.z())); + g += Vector(0.0, 0.0, sign(p.z())); } const FT length_sq = g.squared_length(); if (length_sq > 0.00001) { @@ -60,18 +62,19 @@ int main() { return g; }; - // create a domain from the grid + // create domain from given grid and gradient auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid, cube_gradient); - // prepare collections for the results + // containers for output indexed surface meshes Point_range points_mc, points_dc; Polygon_range polygons_mc, polygons_dc; - // execute topologically correct marching cubes and dual contouring with an isovalue of 0.8 - CGAL::Isosurfacing::marching_cubes(domain, 0.88, points_mc, polygons_mc); - CGAL::Isosurfacing::dual_contouring(domain, 0.88, points_dc, polygons_dc); + // run topologically correct marching cubes and dual contouring with given isovalue + const FT isovalue = 0.88; + CGAL::Isosurfacing::marching_cubes(domain, isovalue, points_mc, polygons_mc, true); + CGAL::Isosurfacing::dual_contouring(domain, isovalue, points_dc, polygons_dc); - // save the results in the OFF format + // save output indexed meshes to files, in the OFF format CGAL::IO::write_OFF("result_mc.off", points_mc, polygons_mc); CGAL::IO::write_OFF("result_dc.off", points_dc, polygons_dc); } diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp index 87d648191d8..31513fbae5d 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_implicit_iwp.cpp @@ -37,19 +37,20 @@ int main() { return Vector(gx, gy, gz); }; - const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1}; - const Vector spacing(0.05f, 0.05f, 0.05f); + const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; + const FT spacing = 0.5; + const Vector vec_spacing(spacing, spacing, spacing); - // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02 + // create a domain with given bounding box and grid spacing auto domain = - CGAL::Isosurfacing::create_implicit_cartesian_grid_domain(bbox, spacing, iwp_value, iwp_gradient); + CGAL::Isosurfacing::create_implicit_cartesian_grid_domain(bbox, vec_spacing, iwp_value, iwp_gradient); // prepare collections for the result Point_range points; Polygon_range polygons; - // execute marching cubes with an isovalue of 0 - CGAL::Isosurfacing::dual_contouring(domain, 0.0f, points, polygons); + // execute marching cubes with an isovalue of 0.0 + CGAL::Isosurfacing::dual_contouring(domain, 0.0, points, polygons); // save the result in the OFF format CGAL::IO::write_OFF("result.off", points, polygons); diff --git a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp index b2bfccf3653..cdd580c1ee6 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/dual_contouring_mesh_offset.cpp @@ -23,7 +23,7 @@ typedef CGAL::AABB_traits Traits; typedef CGAL::AABB_tree Tree; typedef std::vector Point_range; -typedef std::vector> Polygon_range; +typedef std::vector > Polygon_range; int main() { const std::string input_name = CGAL::data_file_path("meshes/cross.off"); @@ -32,7 +32,7 @@ int main() { Mesh mesh_input; if (!CGAL::IO::read_OFF(input_name, mesh_input)) { - std::cerr << "Could not read mesh" << std::endl; + std::cerr << "Could not read input mesh" << std::endl; exit(-1); } @@ -47,7 +47,7 @@ int main() { CGAL::Side_of_triangle_mesh::type> sotm(mesh_input); - + // functors for addressing distance and normal queries auto mesh_distance = [&tree](const Point& p) { const Point& x = tree.closest_point(p); return std::sqrt((p - x).squared_length()); @@ -55,18 +55,20 @@ int main() { auto mesh_normal = [&tree](const Point& p) { const Point& x = tree.closest_point(p); - const Vector n = p - x; - return n / std::sqrt(n.squared_length()); + const Vector n = p - x; // TODO: address case where norm = zero + return n / std::sqrt(n.squared_length()); // normalize output vector }; - // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02 + // create a domain with given bounding box and grid spacing auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain(aabb_grid, grid_spacing, mesh_distance, mesh_normal); - + // containers for output indexed surface mesh Point_range points; Polygon_range polygons; + // dual contouring CGAL::Isosurfacing::dual_contouring(domain, offset_value, 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/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp index 1811eea7218..57c660f05b9 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_cartesian_grid_sphere.cpp @@ -16,10 +16,10 @@ typedef std::vector> Polygon_range; int main() { // create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3 - const CGAL::Bbox_3 bbox(-1, -1, -1, 1, 1, 1); + const CGAL::Bbox_3 bbox(-1.0, -1.0, -1.0, 1.0, 1.0, 1.0); std::shared_ptr grid = std::make_shared(50, 50, 50, bbox); - // calculate the value at all grid points + // calculate function values at all grid points for (std::size_t x = 0; x < grid->xdim(); x++) { for (std::size_t y = 0; y < grid->ydim(); y++) { for (std::size_t z = 0; z < grid->zdim(); z++) { @@ -28,7 +28,7 @@ int main() { const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin(); const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin(); - // distance to the origin + // Euclidean distance to the origin grid->value(x, y, z) = std::sqrt(pos_x * pos_x + pos_y * pos_y + pos_z * pos_z); } } @@ -44,6 +44,6 @@ int main() { // execute marching cubes with an isovalue of 0.8 CGAL::Isosurfacing::marching_cubes(domain, 0.8f, points, polygons); - // save the result in the OFF format + // save the result 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_implicit_sphere.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp index 7c6215e2f4d..284e0c4ce0e 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -1,4 +1,3 @@ - #include #include #include @@ -6,28 +5,31 @@ #include typedef CGAL::Simple_cartesian Kernel; +typedef typename Kernel::FT FT; typedef typename Kernel::Vector_3 Vector; typedef typename Kernel::Point_3 Point; typedef std::vector Point_range; -typedef std::vector> Polygon_range; +typedef std::vector > Polygon_range; int main() { - const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1}; - const Vector spacing(0.04f, 0.04f, 0.04f); + const CGAL::Bbox_3 bbox{-1.0, -1.0, -1.0, 1.0, 1.0, 1.0}; + const FT spacing = 0.04; + const Vector vec_spacing(spacing, spacing, spacing); + // Euclidean distance function to the origin auto sphere_function = [&](const Point& p) { return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); }; - // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02 - auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain(bbox, spacing, sphere_function); + // create a domain with given bounding box and grid spacing + auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain(bbox, vec_spacing, sphere_function); - // prepare collections for the result + // prepare collections for the output indexed mesh Point_range points; Polygon_range polygons; // execute marching cubes with an isovalue of 0.8 - CGAL::Isosurfacing::marching_cubes(domain, 0.8f, points, polygons); + CGAL::Isosurfacing::marching_cubes(domain, 0.8, points, polygons); - // save the result in the OFF format + // 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 a17d2390bc0..bccfee1ad80 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp @@ -6,7 +6,6 @@ typedef CGAL::Simple_cartesian Kernel; typedef typename Kernel::Point_3 Point; - typedef CGAL::Cartesian_grid_3 Grid; typedef std::vector Point_range; @@ -16,26 +15,26 @@ int main() { const std::string fname = CGAL::data_file_path("images/skull_2.9.inr"); - // load the image + // load volumetric image from a file CGAL::Image_3 image; if (!image.read(fname)) { - std::cerr << "Error: Cannot read file " << fname << std::endl; + std::cerr << "Error: Cannot read image file " << fname << std::endl; return EXIT_FAILURE; } - // convert it to a cartesian grid + // convert image to a cartesian grid std::shared_ptr grid = std::make_shared(image); // create a domain from the grid auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid); - // prepare collections for the result + // prepare collections for the output indexed mesh Point_range points; Polygon_range polygons; // execute marching cubes with an isovalue of 2.9 CGAL::Isosurfacing::marching_cubes(domain, 2.9, points, polygons); - // save the result in the OFF format + // save output 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_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp index 0f13251d035..44906dd347a 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp @@ -28,7 +28,8 @@ typedef std::vector Point_range; typedef std::vector> Polygon_range; -// computes the distance of a point p from the mesh with the use of a AABB_tree +// computes the Euclidean distance from query point p to the mesh +// via the AABB tree data structure inline Kernel::FT distance_to_mesh(const Tree& tree, const Point& p) { const Point& x = tree.closest_point(p); return std::sqrt((p - x).squared_length()); @@ -39,25 +40,25 @@ int main() { const int n_voxels = 20; const FT offset_value = 0.2; - // load the original mesh + // load input mesh Mesh mesh_input; if (!CGAL::IO::read_OFF(input_name, mesh_input)) { - std::cerr << "Could not read mesh" << std::endl; + std::cerr << "Could not read input mesh" << std::endl; exit(-1); } - // compute the new bounding box + // compute loose bounding box of the mesh CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input); - Vector aabb_increase_vec = Vector(offset_value + 0.01, offset_value + 0.01, offset_value + 0.01); + const FT loose_offset = offset_value + 0 .01; + Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset); aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox(); aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox(); - // construct AABB tree + // construct AABB tree and functor to address inside/outside point queries Tree tree(mesh_input.faces_begin(), mesh_input.faces_end(), mesh_input); - CGAL::Side_of_triangle_mesh::type> sotm(mesh_input); - // create the grid + // create grid std::shared_ptr grid = std::make_shared(n_voxels, n_voxels, n_voxels, aabb_grid); for (std::size_t z = 0; z < grid->zdim(); z++) { @@ -69,28 +70,28 @@ int main() { const FT pos_z = z * grid->get_spacing()[2] + grid->get_bbox().zmin(); const Point p(pos_x, pos_y, pos_z); - // compute the distance + // compute unsigned distance to input mesh grid->value(x, y, z) = distance_to_mesh(tree, p); - // flip the sign, so the distance is negative inside the mesh + // sign distance so that it is negative inside the mesh const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE); if (is_inside) { - grid->value(x, y, z) *= -1; + grid->value(x, y, z) *= -1.0; } } } } - // create a domain from the grid + // create domain from the grid auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain(grid); - // prepare collections for the result + // prepare collections for output indexed surface mesh Point_range points; Polygon_range polygons; - // execute marching cubes with an isovalue equal to the offset + // execute marching cubes with an isovalue equating offset CGAL::Isosurfacing::marching_cubes(domain, offset_value, points, polygons); - // save the result in the OFF format + // save output indexed surface mesh to a file, in the OFF format CGAL::IO::write_OFF("result.off", points, polygons); } diff --git a/Isosurfacing_3/package_info/Isosurfacing_3/copyright b/Isosurfacing_3/package_info/Isosurfacing_3/copyright index 434d5723b62..fd4b0243d94 100644 --- a/Isosurfacing_3/package_info/Isosurfacing_3/copyright +++ b/Isosurfacing_3/package_info/Isosurfacing_3/copyright @@ -1 +1 @@ -INRIA Sophia-Antipolis (France) +Inria Sophia-Antipolis (France) diff --git a/Isosurfacing_3/package_info/Isosurfacing_3/description.txt b/Isosurfacing_3/package_info/Isosurfacing_3/description.txt index bfc5002e1e7..a5b1968e688 100644 --- a/Isosurfacing_3/package_info/Isosurfacing_3/description.txt +++ b/Isosurfacing_3/package_info/Isosurfacing_3/description.txt @@ -1 +1 @@ -3D Isosurfacing Algorithms \ No newline at end of file +The 3D Isosurfacing package provides several isosurfacing algorithms (marching cubes, dual contouring) to generate surface triangle and quadrangle surface meshes from an input 3D domain and a user-defined isovalue. \ No newline at end of file diff --git a/Isosurfacing_3/package_info/Isosurfacing_3/long_description.txt b/Isosurfacing_3/package_info/Isosurfacing_3/long_description.txt index e947047a135..6650e546139 100644 --- a/Isosurfacing_3/package_info/Isosurfacing_3/long_description.txt +++ b/Isosurfacing_3/package_info/Isosurfacing_3/long_description.txt @@ -1 +1 @@ -Various implementation of the marching cube algorithms to generate triangle surfaces +This component takes a 3D domain as input and a user-specified isovalue, and generates a surface mesh approximating the specified isovalue. The meshing algorithm can be chosen among two isosurfacing algorithms: marching cubes and dual contouring. Two variants of the marching cubes algorithm are offererd: with or without topological guarantees. The domain can be created from an explicit grid, an implicit grid or a volumetric image.