diff --git a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt b/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt index a689d15bcc0..c5405002e97 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/Isosurfacing_3.txt @@ -7,48 +7,61 @@ namespace CGAL { \cgalAutoToc \author Julian Stahl and Daniel Zint +gif with different isosurfaces zooming in or out \section secmyintroduction Introduction This package provides functions to compute a surface mesh representing an isosurface. -The data structure from which an isosurface can be extracted is a 3-dimensional scalar function. -An isosurface is defined as the surface on which the value of this function is equal to a given constant. -This constant value is called the isovalue. -The representation that is used to store the isosurface is a triangular or quadrilateral indexed face set. Or triangle soup? +The data structure from which an isosurface is extracted is a 3-dimensional scalar function. +An isosurface is defined as the surface on which the value of this function is equal to a given constant, i.e. the isovalue. +The isosurface is stored as a triangular or quadrilateral indexed face set. \section secmyalgorithms Algorithms -There are multiple algorithms to extract an isosurface from a uniform grid or octree. -This package contains marching cubes, topologically correct marching cubes, dual contouring and octree marching?. +There are multiple algorithms to extract isosurfaces. +This package contains Marching Cubes, topologically correct Marching Cubes, Dual Contouring and Octree Marching. (References?) -\subsection subsecmc Marching cubes -The marching cubes algorithm iterates over all cells of the input domain and processes each cell individually. -Firstly, the values of all eight corners of a cell are compared to the isovalue. -Each corner gets a sign (+/-) to indicate if it is outside or inside the isosurface. -On every edge of the cube that connects corners with different signs a new vertex is created via linear interpolation. +\subsection subsecmc Marching Cubes +The Marching Cubes algorithm 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. -
- -
+Marching Cubes is an old algorithm that is fast but approximate. It can be run on any domain that has cubical cells. +The algorithm creates a manifold triangle mesh if the domain is a regular grid. +This guarantees that every edge is incident to at most two faces and every vertex is surrounded by one triangle fan. +As long as the isosurface does not intersect the domains boundaries the resulting mesh is also watertight. +If the input domain contains multiple isosurfaces with the same value that are not connected Marching Cubes will also generate multiple meshes. +Compared to other approaches the algorithm often generates more triangles and many thin or acute triangles. +Sharp edges are not well preserved in a mesh generated with Marching Cubes. -\subsection subsectmc Topologically correct marching cubes +\image html mc_cases.png -\subsection subsecdc Dual contouring +\subsection subsectmc Topologically correct Marching Cubes +This algorithm is an extension to Marching Cubes that provides additional guarantees. +It generates a mesh that is homeomorph to the trilinear interpolant of the input function inside each cube. ... + +TODO examples + +\subsection subsecdc Dual Contouring +The Dual Contouring algorithm first iterates over all cells. If a cell intersects the isosurface a vertex inside this cell is created. +Then it checks for each edge if it intersects the isosurface. If that is the case the vertices of all incident cells are connects to form a face. + +This algorithm needs a functor that returns the normals of vertices as an additional parameter. +The placement of a vertex position inside a cell is also configurable with an optional parameter. ... + +Dual Contouring works on any domain but does not guarantee a manifold or watertight mesh. +It creates less faces than Marching Cubes. +The algorithms can represent sharp edges better than Marching Cubes. \subsection subseccomparison Comparison table with pros/cons or applicability and guarantees -\subsection subsecparameters Parameters -gif with different isosurfaces zooming in or out - -\section secmydomains Domains -necessary? - \section secmyinterface Interface -The algorithms can be called with their respective functions. The parameters are always the same: +Each algorithm is represented by a single functions. The function signature is the same for each algorithm (except for some additional parameters for some of them): \code{.cpp} template @@ -56,22 +69,43 @@ void make_triangle_mesh_using_marching_cubes(const Domain_& domain, const typena PointRange& points, PolygonRange& polygons); \endcode -The `domain` is an object that provides functions to access the input data and its topology. -There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm. -Examples are `Cartesian_grid_domain`, which takes a `Cartesian_grid_3` as parameter, -and `Implicit_domain`, which represents the grid as an implicit function. +The input is provided in the form of a `domain` (see \ref secmydomains). The `iso_value` parameter describes the grid value the isosurface should represent. -The output is in the form of a polygon soup that is written to the two collections `points` and `polygons`. -The vertices are stored as `Point_3` in `points`. Each face in `polygons` is a list of indices pointing into the `points` collection. +The output is in the form of an indexed face set that is written to the two collections `points` and `polygons`. +The vertex positions are stored as `Point_3` in `points`. Each face in `polygons` is a list of indices pointing into the `points` collection. +Depending on the algorithm, the indexed face set may either store a polygon soup or a topological mesh. Algorithms can run sequentially on one CPU core or in parallel. -The Concurrency_tag can be Sequential_tag or Parallel_tag and is used to specify how the algorithm is executed. +The Concurrency_tag is used to specify how the algorithm is executed and is either Sequential_tag or Parallel_tag. To enable parallelism, CGAL needs to be linked with the Intel TBB library. If the parallel version is not availible the sequential version will always be used as a fallback. +\section secmydomains Domains + +A domain is an object that provides functions to access the input data, its geometry, topology, and optionally its gradient. +There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm. +Users can also define new domains by implementing the `Isosurfacing_domain` concept. + +\subsection mysubsecimplicitdomain Implicit domain +The `Implicit_domain` represents the input function in an implicit form without storing any values. +It takes a functor that computes the value of the function from the position of a vertex as parameter. +Additionally, the bounding box and spacing between grid points have to be specified. + +\subsection mysubseccartesiandomain Cartesian grid domain +The `Cartesian_grid_domain` only takes a `Cartesian_grid_3` as parameter. +It represents a uniform grid of values that are either computed by the user or read from an `Image_3`. +The constructor of `Cartesian_grid_3` needs the number of grid points in each dimension and the bounding box. +The values are read and written with `value(x, y, z)` where x, y, and z are the coordinates of a grid point. +Alternatively, all required data can be copied from an `Image_3`. + +\subsection mysubsecoctreedomain Octree domain +The `Octree_domain` wraps an octree to be used by isosurfacing algorithms. +The octree has to be already refined. ... + + \section secmyexamples Examples \subsection myExampleImplicit_domain Implicit sphere @@ -81,17 +115,33 @@ The domain is an `Implicit_domain` that describes a sphere by the distance to it \cgalExample{Isosurfacing_3/marching_cubes_implicit_sphere.cpp} -\subsection myExampleCartesian_grid_domain Cartesian_grid_3 sphere +\subsection myExampleAll_cube Cartesian_grid_3 cube comparison -The following example shows the use of the marching cubes algorithm to extract an isosurface. -The domain is an `Cartesian_grid_domain` that describes a sphere by storing the distance to its origin in a `Cartesian_grid_3`. +The following example compares all provided algorithms to extract an isosurface. +The domain is an `Cartesian_grid_domain` that describes a cube by storing the manhattan distance to its origin in a `Cartesian_grid_3`. -\cgalExample{Isosurfacing_3/dual_contouring_cartesian_grid.cpp} +\cgalExample{Isosurfacing_3/all_cartesian_cube.cpp} + +\subsection myExampleImage Image_3 + +The following example shows how to load data from an `Image_3`. \cgalExample{Isosurfacing_3/marching_cubes_inrimage.cpp} +\subsection myExampleMeshOffset Offset mesh + +The following example shows how to compute an offset mesh. +The original mesh is passed to an `AABB_tree` to allow fast distance queries. +With the use of `Side_of_triangle_mesh` the sign of the distance function is flipped inside the mesh. + \cgalExample{Isosurfacing_3/marching_cubes_mesh_offset.cpp} +\subsection myExampleOctree Octree Dual Contouring + +The following example shows how to extract an isosurface from an octree using Dual Contouring. +The domain is an `Octree_domain` that describes a sphere by the distance to its origin stored in an octree. +The octree is highly refined in one octant and only coarse in the others. + \cgalExample{Isosurfacing_3/dual_contouring_octree.cpp} */ diff --git a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt b/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt index 3a09957b302..a4d7ac7bc2d 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/PackageDescription.txt @@ -22,7 +22,7 @@ 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 output is an indexed face set that stores an isosurface in the form of a surface mesh. -Available algorithms include marching cubes, dual contouring, and others. +Available algorithms include Marching Cubes, Dual Contouring, and others. \cgalClassifedRefPages diff --git a/Isosurfacing_3/doc/Isosurfacing_3/examples.txt b/Isosurfacing_3/doc/Isosurfacing_3/examples.txt index a7e560287f2..b5b64b2a4e0 100644 --- a/Isosurfacing_3/doc/Isosurfacing_3/examples.txt +++ b/Isosurfacing_3/doc/Isosurfacing_3/examples.txt @@ -6,4 +6,5 @@ \example Isosurfacing_3/dual_contouring_cartesian_grid.cpp \example Isosurfacing_3/dual_contouring_mesh_offset.cpp \example Isosurfacing_3/dual_contouring_octree.cpp +\example Isosurfacing_3/all_cartesian_cube.cpp */ diff --git a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt index 831ac0b0528..2c54d61b84c 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt +++ b/Isosurfacing_3/examples/Isosurfacing_3/CMakeLists.txt @@ -13,6 +13,7 @@ create_single_source_cgal_program( "marching_cubes_inrimage.cpp" ) create_single_source_cgal_program( "dual_contouring_cartesian_grid.cpp" ) create_single_source_cgal_program( "dual_contouring_mesh_offset.cpp" ) create_single_source_cgal_program( "dual_contouring_octree.cpp" ) +create_single_source_cgal_program( "all_cartesian_cube.cpp" ) find_package(TBB) include(CGAL_TBB_support) @@ -24,4 +25,5 @@ if(TARGET CGAL::TBB_support) target_link_libraries(dual_contouring_cartesian_grid PUBLIC CGAL::TBB_support) target_link_libraries(dual_contouring_mesh_offset PUBLIC CGAL::TBB_support) target_link_libraries(dual_contouring_octree PUBLIC CGAL::TBB_support) + target_link_libraries(all_cartesian_cube PUBLIC CGAL::TBB_support) endif() diff --git a/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp new file mode 100644 index 00000000000..41e8acc2e2a --- /dev/null +++ b/Isosurfacing_3/examples/Isosurfacing_3/all_cartesian_cube.cpp @@ -0,0 +1,63 @@ +#include +#include +#include +#include +#include +#include +#include + +typedef CGAL::Simple_cartesian Kernel; +typedef typename Kernel::FT FT; +typedef typename Kernel::Point_3 Point; +typedef typename Kernel::Vector_3 Vector; + +typedef CGAL::Cartesian_grid_3 Grid; + +typedef std::vector Point_range; +typedef std::vector> Polygon_range; + +int main() { + // create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3 + Grid grid(100, 100, 100, {-1, -1, -1, 1, 1, 1}); + + // calculate the value 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++) { + + const FT pos_x = x * grid.get_spacing()[0] + grid.get_bbox().xmin(); + const FT pos_y = y * grid.get_spacing()[1] + grid.get_bbox().ymin(); + const FT pos_z = z * grid.get_spacing()[2] + grid.get_bbox().zmin(); + + // manhattan distance to the origin + grid.value(x, y, z) = std::max({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)}); + + // the normal depends on the side of the cube + if (grid.value(x, y, z) == std::abs(pos_x)) { + grid.gradient(x, y, z) = Vector(std::copysign(1.0, pos_x), 0, 0); + } else if (grid.value(x, y, z) == std::abs(pos_y)) { + grid.gradient(x, y, z) = Vector(0, std::copysign(1.0, pos_y), 0); + } else if (grid.value(x, y, z) == std::abs(pos_z)) { + grid.gradient(x, y, z) = Vector(0, 0, std::copysign(1.0, pos_z)); + } + } + } + } + + // create a domain from the grid + CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); + + // prepare collections for the results + Point_range points_mc, points_tmc, points_dc; + Polygon_range polygons_mc, polygons_tmc, polygons_dc; + + // execute marching cubes, topologically correct marching cubes and dual contouring with an isovalue of 0.8 + CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 0.8, points_mc, polygons_mc); + CGAL::Isosurfacing::make_triangle_mesh_using_tmc(domain, 0.8, points_tmc, polygons_tmc); + CGAL::Isosurfacing::make_quad_mesh_using_dual_contouring(domain, 0.8, points_dc, polygons_dc); + + // save the results in the OFF format + CGAL::IO::write_OFF("result_mc.off", points_mc, polygons_mc); + CGAL::IO::write_OFF("result_tmc.off", points_tmc, polygons_tmc); + CGAL::IO::write_OFF("result_dc.off", points_dc, polygons_dc); +} 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 c8e608e4833..83cd79015f7 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_implicit_sphere.cpp @@ -19,7 +19,7 @@ int main() { // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02 CGAL::Isosurfacing::Implicit_domain(sphere_function))> domain( - {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient(sphere_function)); // TODO + {-1, -1, -1, 1, 1, 1}, Vector(0.02f, 0.02f, 0.02f), sphere_function, CGAL::Isosurfacing::Default_gradient(sphere_function)); // TODO: this is ugly // prepare collections for the result Point_range points; diff --git a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp index 5ca98e4393a..acf29468518 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_inrimage.cpp @@ -16,21 +16,27 @@ typedef std::vector> Polygon_range; int main() { const std::string fname = "../data/skull_2.9.inr"; - // Load image + + // load the image CGAL::Image_3 image; if (!image.read(fname)) { std::cerr << "Error: Cannot read file " << fname << std::endl; return EXIT_FAILURE; } + // convert it to a cartesian grid const Grid grid(image); + // create a domain from the grid CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); + // prepare collections for the result Point_range points; Polygon_range polygons; + // execute marching cubes with an isovalue of 0.8 CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, 2.9f, 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/marching_cubes_mesh_offset.cpp b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp index ab7a8f05ec3..d9ea07d24d7 100644 --- a/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp +++ b/Isosurfacing_3/examples/Isosurfacing_3/marching_cubes_mesh_offset.cpp @@ -28,6 +28,7 @@ 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 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()); @@ -38,13 +39,14 @@ int main() { const int n_voxels = 20; const FT offset_value = 0.01; + // load the original mesh Mesh mesh_input; if (!CGAL::IO::read_OFF(input_name, mesh_input)) { std::cerr << "Could not read mesh" << std::endl; exit(-1); } - // compute bounding box + // compute the new bounding box 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); aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox(); @@ -55,10 +57,9 @@ int main() { CGAL::Side_of_triangle_mesh::type> sotm(mesh_input); + // create the grid Grid grid(n_voxels, n_voxels, n_voxels, aabb_grid); - CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); - for (std::size_t z = 0; z < grid.zdim(); z++) { for (std::size_t y = 0; y < grid.ydim(); y++) { for (std::size_t x = 0; x < grid.xdim(); x++) { @@ -68,8 +69,10 @@ 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 grid.value(x, y, z) = distance_to_mesh(tree, p); + // flip the sign, so the distance is negative inside the mesh const bool is_inside = (sotm(p) == CGAL::ON_BOUNDED_SIDE); if (is_inside) { grid.value(x, y, z) *= -1; @@ -78,10 +81,16 @@ int main() { } } + // create a domain from the grid + CGAL::Isosurfacing::Cartesian_grid_domain domain(grid); + + // prepare collections for the result Point_range points; Polygon_range polygons; + // execute marching cubes with an isovalue equal to the offset CGAL::Isosurfacing::make_triangle_mesh_using_marching_cubes(domain, offset_value, points, polygons); + // save the result in the OFF format CGAL::IO::write_OFF("result.off", points, polygons); }