CGAL 6.0 - 3D Isosurfacing
Loading...
Searching...
No Matches
User Manual

Author
Mael Rouxel-Labbé, Julian Stahl, Daniel Zint, and Pierre Alliez

Figure 61.1 Generating a surface from a 3D gray level image using Marching Cubes (3D input image: qim.dk)


Introduction

Given a field of a scalar values, an isosurface is defined as the locus of points where the field has a given constant value; in other words, it is a level set. This constant value is referred to as the "isovalue", and, for well-behaved fields, the level set forms a surface. In the following, we shall refer to the the field of scalar values as the value field. "Isosurfacing", also known as "isosurface extraction" or "contouring", is the process of constructing the isosurface corresponding to a given value field and isovalue. Isosurfacing is often needed for volume visualization and for the simulation of physical phenomena.

This CGAL package provides methods to extract isosurfaces from 3D value fields. These contouring techniques rely on the partition of the 3D space and of the field to construct an approximate representation of the isosurface. The 3D value field can be described through various representations: an implicit function, an interpolated set of discrete sampling values, etc. (see Examples). The isovalue is user-defined. The output is a polygon soup, made either of triangles or quads depending on the method, and may consist of a single connected component, or multiple, disjoint components. Note that due to the inherent approximate nature of these discrete methods that only sample 3D value fields, parts of the "true" isosurface may be missing from the output, and the output may contain artifacts that re not present in the true isosurface.

Isosurfacing Methods

The scientific literature abounds with algorithms for extracting isosurfaces, each coming with different properties for the output and requirements for the input [2]. This package offers the following methods

Marching Cubes (MC)

Marching Cubes (MC) [5] uses a volumetric grid, i.e., a 3D iso-cuboid partitioned into hexahedral cells. All cells of the grid are processed individually using values of the value field sampled at the grid corners. Each cell corner is assigned a sign (+/-) to indicate whether its value field value is above or below the user-defined isovalue. A vertex is created along each grid edge where a sign change occurs, i.e., where the edge intersects the isosurface. More specifically, the vertex location is computed via linear interpolation of the value field values evaluated at the cell corners forming the edge. These vertices are connected to form triangles within the cell, depending on the configuration of signs at the cell corners. Figure 61.2 illustrates the configurations in 2D. In 3D, there are no less than 33 configurations (not shown) [1].

Figure 61.2 Examples of some configurations for 2D Marching Cubes.


The implementation within CGAL is generic in the sense that it can process any grid-like data structure that consists of hexahedral cells. When the hexahedral grid is a conforming grid (meaning that the intersection of two hexahedral cells is a face, an edge, or a vertex), the Marching Cubes algorithm generates as output a surface triangle mesh that is almost always combinatorially 2-manifold, yet can still exhibit topological issues such as holes or non-conforming edges.

If the mesh is 2-manifold and the isosurface does not intersect the domain boundary, then the output mesh is watertight. As the Marching Cubes algorithm uses linear interpolation of the sampled value field along the grid edges, it can miss details or components that are not captured by the sampling of the value field.

Compared to other meshing approaches such as Delaunay refinement, Marching Cubes is substantially faster, but often tends to generate more triangle facets for an equivalent desired sizing field. In addition, the quality of the triangle facets is in general poor, with many needle/cap shaped triangles.

Furthermore, Marching Cubes does not preserve the sharp features present in the isovalue of the input value field (see Figure 61.4).

Topologically Correct Marching Cubes (TMC)

Topologically Correct Marching Cubes is an extension to the Marching Cubes algorithm which provides additional guarantees for the output [3]. More specifically, it generates as output a mesh that is homeomorphic to the trilinear interpolant of the input value field inside each cube. This means that the output mesh can accurately represent small complex features. For example, a tunnel of the isosurface within a single cell is topologically resolved. To achieve this, the algorithm can insert additional vertices within cells. Furthermore, the mesh is guaranteed to be 2-manifold and watertight, as long as the isosurface does not intersect the domain boundaries.

Figure 61.3 Marching Cubes vs Topologically Correct Marching Cubes.


Dual Contouring (DC)

Dual Contouring (DC) [4] is a method that does not generate vertices on the grid edges, but within cells instead. A facet is created for each edge that intersects the isosurface by connecting the vertices of the incident cells. For a uniform hexahedral grid, this results in a quadrilateral surface mesh. Dual Contouring can deal with any domain but guarantees neither a 2-manifold nor a watertight mesh. On the other hand, it generates fewer faces and higher quality faces than Marching Cubes, in general. Finally, its main advantage over Marching Cubes is its ability to recover sharp creases and corners.

Figure 61.4 Comparison between a mesh of a CSG shape generated by Marching Cubes (left) and Dual Contouring (right).


In addition to the 3D value field, Dual Contouring requires knowledge about the gradient of the value field.

The CGAL implementation uses a vertex positioning strategy based on Quadric Error Metrics: for a cell, the vertex location is computed by minimizing the error to the sum of the quadrics defined at each edge-isosurface intersection. Using this approach, the vertex may be located outside the cell, which is a desirable property to improve the odds of recovering sharp features, but it might also create self-intersections. Users can choose to constrain the vertex location inside the cell.

By default, Dual Contouring generates quads, but using edge-isosurface intersections, one can "star" these quads to generate four triangles. Triangulating the quads is the default behavior in CGAL, but this can be changed via named parameters.

Comparisons

The following table summarizes the differences between the algorithms in terms of constraints over the input 3D domain, the facets of the output surface mesh, and the properties of the output surface mesh.

Algorithm Facets 2-Manifold Watertight* Topologically correct Recovery of Sharp Features
MC Triangles no no no no
TMC Triangles yes yes yes no
DC Polygons no no no yes (not guaranteed)

(* assuming the isosurface does not exit the specified bounding box of the input 3D domain)

Note that the output mesh has boundaries when the isosurface intersects the domain boundaries, regardless of the method (see Figure 61.5).

Figure 61.5 Outputs of Marching Cubes (left) and Dual Contouring (right) for an implicit sphere of radius 1.1 and a domain of size 2x2x2, both centered at the origin. Output meshes can have boundaries when the isosurface intersects the domain boundary.


Interface

The following functions are the main entry points to the isosurfacing algorithms:

These three free functions share the same signature:

template <typename ConcurrencyTag = CGAL::Sequential_tag,
typename Domain,
typename PointRange,
typename PolygonRange>
void ...(const Domain& domain,
const typename Domain::FT isovalue,
PointRange& points,
PolygonRange& polygons);

The input (space partition, value field, gradient field) is provided in the form of a domain, see Domains for a complete description.

The isovalue scalar parameter is the value that defines the isosurface being approximated.

The output discrete surface is provided in the form of a polygon soup, which is stored into two containers: points and polygons. Depending on the algorithm, the polygon soup may store either unorganized polygons with no relationship to one another (ie, no connectivity is shared between them) or polygons sharing points (the same point in adjacent polygons will be the same point in the point range).

The isosurfacing algorithms can run either sequentially in one thread or in parallel via multithreading. The template parameter ConcurrencyTag is used to specify how the algorithm is executed. To enable parallelism, CGAL must be linked with the Intel TBB library (see the CMakeLists.txt file in the examples folder).

Domains

A domain is an object that provides functions to access the partition of the 3D volume, the value field, and, optionally, the gradient field at a given point query. These requirements are described through two concepts: IsosurfacingDomain_3 and IsosurfacingDomainWithGradient_3.

Two domains, CGAL::Isosurfacing::Marching_cubes_domain_3 and CGAL::Isosurfacing::Dual_contouring_domain_3, are provided as the respective default class models that fulfill the requirements of the concepts. Both these domain models have template parameters enabling the user to customize the domain:

Performance

Due to their cell-based nature, the isosurfacing algorithms are well-suited for parallel execution.

Examples

The first two examples are very basic examples for Marching Cubes and Dual Contouring. Afterwards, the focus is shifted from the method to the type of input data, and examples run both methods on different types of input data.

Marching Cubes

The following example illustrates a basic run of the Marching Cubes algorithm, and in particular the free function to create a domain from a Cartesian grid, and the named parameter that enables the user to switch from Marching Cubes to Topologically Correct Marching Cubes.


File Isosurfacing_3/marching_cubes.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/Value_function_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
// create bounding box and grid
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
Grid grid { bbox, CGAL::make_array<std::size_t>(30, 30, 30) };
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
// fill up values
auto sphere_value_fn = [](const Point& p) -> FT
{
return sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
};
Values values { sphere_value_fn, grid };
// Below is equivalent to:
// Domain domain { grid, values };
Point_range points;
Polygon_range triangles;
// run marching cubes isosurfacing
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles,
CGAL::parameters::use_topologically_correct_marching_cubes(true));
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("marching_cubes.off", points, triangles);
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
The class Cartesian_grid_3 represents a 3D Cartesian grid, that is the partition of an iso-cuboid int...
Definition: Cartesian_grid_3.h:148
A domain that can be used with the Marching Cubes algorithm.
Definition: Marching_cubes_domain_3.h:50
The class Value_function_3 represents a field of scalars computed using a user-provided unary functio...
Definition: Value_function_3.h:40
bool write_polygon_soup(const std::string &fname, const PointRange &points, const PolygonRange &polygons, const NamedParameters &np=parameters::default_values())
Marching_cubes_domain_3< Partition, ValueField, EdgeIntersectionOracle > create_marching_cubes_domain_3(const Partition &partition, const ValueField &values, const EdgeIntersectionOracle &intersection_oracle=EdgeIntersectionOracle())
creates a new instance of a domain that can be used with the Marching Cubes algorithm.
Definition: Marching_cubes_domain_3.h:92
NT sqrt(const NT &x)

Dual Contouring

The following example illustrates a basic run of the Dual Contouring algorithm, and in particular the free function to create a domain from a Cartesian grid, and the named parameters that enable (or disable) triangulation of the output, and to constrain the vertex location within the cell.


File Isosurfacing_3/dual_contouring.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/Value_function_3.h>
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
// https://www-sop.inria.fr/galaad/surface/
auto devil_value = [](const Point& point)
{
const FT x = point.x(), y = point.y(), z = point.z();
return x*x*x*x + 2*x*x*z*z - 0.36*x*x - y*y*y*y + 0.25*y*y + z*z*z*z;
};
auto devil_gradient = [](const Point& point)
{
const FT x = point.x(), y = point.y(), z = point.z();
const FT gx = 4*x*x*x + 4*x*z*z - 0.72*x;
const FT gy = -4*y*y*y + 0.5*y;
const FT gz = 4*x*x*z + 4*z*z*z;
Vector g(gx, gy, gz);
return g / std::sqrt(gx*gx + gy*gy + gz*gz);
};
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.;
// create bounding box and grid
const CGAL::Bbox_3 bbox { -1, -1, -1, 1, 1, 1 };
Grid grid { bbox, CGAL::make_array<std::size_t>(50, 50, 50) };
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
Values values { devil_value, grid };
Gradients gradients { devil_gradient, grid };
// Below is equivalent to:
// Domain domain { grid, values, gradients };
Domain domain = CGAL::Isosurfacing::create_dual_contouring_domain_3(grid, values, gradients);
Point_range points;
Polygon_range triangles;
// run dual contouring isosurfacing
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
CGAL::Isosurfacing::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles,
CGAL::parameters::do_not_triangulate_faces(true)
.constrain_to_cell(false));
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring.off", points, triangles);
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
A domain that can be used as input in the Dual Contouring algorithm.
Definition: Dual_contouring_domain_3.h:52
The class Gradient_function_3 represents a field of vectors computed using a user-provided unary func...
Definition: Gradient_function_3.h:39
Dual_contouring_domain_3< Partition, ValueField, GradientField, EdgeIntersectionOracle > create_dual_contouring_domain_3(const Partition &partition, const ValueField &values, const GradientField &gradients, const EdgeIntersectionOracle &intersection_oracle=EdgeIntersectionOracle())
creates a new instance of a domain that can be used with the Dual Contouring algorithm.
Definition: Dual_contouring_domain_3.h:99
CGAL::Bbox_3 bbox(const PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())

Figure 61.7 Results of the Dual Contouring algorithm: untriangulated (left column) or triangulated (right column), unconstrained vertex location (top row) or constrained vertex location (bottom row).


Dual Contouring using Octree

The following example shows the use of an octree for dual contouring or marching cubes.


File Isosurfacing_3/dual_contouring_octree.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/Isosurfacing_3/Value_function_3.h>
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
#include <CGAL/Isosurfacing_3/Octree_partition.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <CGAL/Real_timer.h>
#include <cmath>
#include <iostream>
#include <vector>
using FT = typename Kernel::FT;
using Vector = typename Kernel::Vector_3;
using Point = typename Kernel::Point_3;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
// Refine one of the octant
struct Refine_one_eighth
{
std::size_t min_depth_;
std::size_t max_depth_;
std::size_t octree_dim_;
Refine_one_eighth(std::size_t min_depth,
std::size_t max_depth)
: min_depth_(min_depth),
max_depth_(max_depth)
{
octree_dim_ = std::size_t(1) << max_depth_;
}
Octree::Global_coordinates uniform_coordinates(const Octree::Node_index& node_index, const Octree& octree) const
{
auto coords = octree.global_coordinates(node_index);
const std::size_t depth_factor = std::size_t(1) << (max_depth_ - octree.depth(node_index));
for(int i=0; i < 3; ++i)
coords[i] *= uint32_t(depth_factor);
return coords;
}
bool operator()(const Octree::Node_index& ni, const Octree& octree) const
{
if(octree.depth(ni) < min_depth_)
return true;
if(octree.depth(ni) == max_depth_)
return false;
auto leaf_coords = uniform_coordinates(ni, octree);
if(leaf_coords[0] >= octree_dim_ / 2)
return false;
if(leaf_coords[1] >= octree_dim_ / 2)
return false;
if(leaf_coords[2] >= octree_dim_ / 2)
return false;
return true;
}
};
auto sphere_function = [](const Point& p) -> FT
{
return std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
};
auto sphere_gradient = [](const Point& p) -> Vector
{
const Vector g = p - CGAL::ORIGIN;
return g / std::sqrt(g.squared_length());
};
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
const CGAL::Bbox_3 bbox{-1., -1., -1., 1., 1., 1.};
std::vector<Kernel::Point_3> bbox_points { {bbox.xmin(), bbox.ymin(), bbox.zmin()},
{ bbox.xmax(), bbox.ymax(), bbox.zmax() } };
CGAL::Real_timer timer;
timer.start();
Octree octree(bbox_points);
Refine_one_eighth split_predicate(3, 5);
octree.refine(split_predicate);
// fill up values and gradients
Values values { sphere_function, octree };
Gradients gradients { sphere_gradient, octree };
Domain domain { octree, values, gradients };
// output containers
Point_range points;
Polygon_range triangles;
// run Dual Contouring
CGAL::Isosurfacing::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles, CGAL::parameters::do_not_triangulate_faces(true));
// run Marching Cubes
// ToDo: Does not yet work with topologically correct marching cubes
// MC_Domain mcdomain{ octree, values };
// CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_if_available_tag>(mcdomain, isovalue, points, triangles);
timer.stop();
std::ofstream oo("octree2.polylines.txt");
oo.precision(17);
octree.dump_to_polylines(oo);
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
std::cout << "Output #vertices (DC): " << points.size() << std::endl;
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring_octree.off", points, triangles);
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
double ymin() const
double xmax() const
double zmin() const
double zmax() const
double ymax() const
double xmin() const
Orthtree< Orthtree_traits_point< GeomTraits, PointRange, PointMap, cubic_nodes, 3 > > Octree
const CGAL::Origin ORIGIN

Implicit Data

The following example shows the usage of Marching Cubes and Dual Contouring algorithms to extract an isosurface. The domain is implicit and describes the unit sphere by the distance to its center (set to the origin) as an implicit 3D value field.


File Isosurfacing_3/contouring_implicit_data.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/Value_function_3.h>
#include <CGAL/Isosurfacing_3/Gradient_function_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/Real_timer.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <array>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
// using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel, CGAL::Isosurfacing::Do_not_cache_vertex_locations>;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
// ---
const FT alpha = 5.01;
auto iwp_value = [](const Point& point)
{
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
return cos(x)*cos(y) + cos(y)*cos(z) + cos(z)*cos(x) - cos(x)*cos(y)*cos(z); // isovalue = 0
};
auto iwp_gradient = [](const Point& point)
{
const FT x = alpha * (point.x() + FT(1.0)) * CGAL_PI;
const FT y = alpha * (point.y() + FT(1.0)) * CGAL_PI;
const FT z = alpha * (point.z() + FT(1.0)) * CGAL_PI;
const FT gx = CGAL_PI * alpha * sin(x) * (cos(y) * (cos(z) - FT(1.0)) - cos(z));
const FT gy = CGAL_PI * alpha * sin(y) * (cos(x) * (cos(z) - FT(1.0)) - cos(z));
const FT gz = CGAL_PI * alpha * sin(z) * (cos(x) * (cos(y) - FT(1.0)) - cos(y));
return Vector(gx, gy, gz);
};
void run_marching_cubes(const Grid& grid,
const FT isovalue)
{
std::cout << "\n ---- " << std::endl;
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
CGAL::Real_timer timer;
timer.start();
// fill up values
Values values { iwp_value, grid };
Domain domain { grid, values };
// output containers
Point_range points;
Polygon_range triangles;
// run Marching Cubes
CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles);
timer.stop();
std::cout << "Output #vertices (MC): " << points.size() << std::endl;
std::cout << "Output #triangles (MC): " << triangles.size() << std::endl;
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
CGAL::IO::write_polygon_soup("marching_cubes_implicit.off", points, triangles);
}
void run_dual_contouring(const Grid& grid,
const FT isovalue)
{
std::cout << "\n ---- " << std::endl;
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
CGAL::Real_timer timer;
timer.start();
// fill up values and gradients
Values values { iwp_value, grid };
Gradients gradients { iwp_gradient, grid };
Domain domain { grid, values, gradients };
// output containers
Point_range points;
Polygon_range triangles;
// run Dual Contouring
CGAL::Isosurfacing::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles,
CGAL::parameters::do_not_triangulate_faces(true));
timer.stop();
std::cout << "Output #vertices (DC): " << points.size() << std::endl;
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
std::cout << "Elapsed time: " << timer.time() << " seconds" << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring_implicit.off", points, triangles);
}
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.;
const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
const FT step = 0.02;
const Vector spacing { step, step, step };
Grid grid { bbox, spacing };
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
run_marching_cubes(grid, isovalue);
run_dual_contouring(grid, isovalue);
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}

Discrete Data

In the following example, the input data is sampled at the vertices of a grid, and interpolated.


File Isosurfacing_3/contouring_discrete_data.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/Interpolated_discrete_gradients_3.h>
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
namespace IS = CGAL::Isosurfacing;
void run_marching_cubes(const Grid& grid,
const FT isovalue)
{
using Domain = IS::Marching_cubes_domain_3<Grid, Values, IS::Linear_interpolation_edge_intersection>;
std::cout << "\n ---- " << std::endl;
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
// fill up values
Values values { grid };
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)
{
const Point& p = grid.point(i,j,k);
const FT d = sqrt(CGAL::squared_distance(p, Point(CGAL::ORIGIN)));
values(i,j,k) = d;
}
}
}
Domain domain { grid, values };
Point_range points;
Polygon_range triangles;
// run marching cubes isosurfacing
IS::marching_cubes<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles);
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("marching_cubes_discrete.off", points, triangles);
}
void run_dual_contouring(const Grid& grid,
const FT isovalue)
{
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Linear_interpolation_edge_intersection>;
std::cout << "\n ---- " << std::endl;
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
// fill up values and gradients
Values values { grid };
Gradients gradients { grid };
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)
{
const Point& p = grid.point(i,j,k);
const FT d = sqrt(CGAL::squared_distance(p, Point(CGAL::ORIGIN)));
values(i,j,k) = d;
if(d != 0)
gradients(i,j,k) = Vector(CGAL::ORIGIN, p) / d;
else
gradients(i,j,k) = CGAL::NULL_VECTOR;
}
}
}
// If the gradient was not known analytically, we could use:
// - Finite_difference_gradient_3 to use finite difference to compute it from values;
// - the function below to use finite difference _and_ store the values at the grid vertices
// gradients.compute_discrete_gradients(values);
Domain domain { grid, values, gradients };
Point_range points;
Polygon_range triangles;
// run dual contouring isosurfacing
IS::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles);
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring_discrete.off", points, triangles);
}
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.8;
// create bounding box and grid
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
Grid grid { bbox, CGAL::make_array<std::size_t>(30, 30, 30) };
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
run_marching_cubes(grid, isovalue);
run_dual_contouring(grid, isovalue);
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
Class template for a gradient field that is computed using discrete values and interpolation.
Definition: Interpolated_discrete_gradients_3.h:40
Class template for a field of values that are calculated using discrete values and interpolation.
Definition: Interpolated_discrete_values_3.h:41
const CGAL::Null_vector NULL_VECTOR
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)
Definition: partition_traits.h:2

3D Image

The following example shows how to load data from an Image_3, and generates an isosurface from this voxel data.


File Isosurfacing_3/contouring_inrimage.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
#include <CGAL/Isosurfacing_3/Interpolated_discrete_values_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/Image_3.h>
#include <CGAL/Isosurfacing_3/IO/Image_3.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <iostream>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
namespace IS = CGAL::Isosurfacing;
void run_marching_cubes(const Grid& grid,
const FT isovalue,
const Values& values)
{
using Domain = IS::Marching_cubes_domain_3<Grid, Values, IS::Linear_interpolation_edge_intersection>;
std::cout << "\n ---- " << std::endl;
std::cout << "Running Marching Cubes with isovalue = " << isovalue << std::endl;
// fill up values
// create a domain from the grid
Domain domain { grid, values };
// prepare collections for the output indexed soup
Point_range points;
Polygon_range triangles;
// execute marching cubes
IS::marching_cubes<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles);
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
// save output indexed mesh to a file, in the OFF format
CGAL::IO::write_polygon_soup("marching_cubes_inrimage.off", points, triangles);
}
void run_dual_contouring(const Grid& grid,
const FT isovalue,
const Values& values)
{
using Domain = IS::Dual_contouring_domain_3<Grid, Values, Gradients, IS::Linear_interpolation_edge_intersection>;
std::cout << "\n ---- " << std::endl;
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
// fill up values and gradients
const FT step = CGAL::approximate_sqrt(grid.spacing().squared_length()) * 0.01; // finite difference step
Gradients gradients { values, step };
Domain domain { grid, values, gradients };
Point_range points;
Polygon_range triangles;
// run dual contouring isosurfacing
IS::dual_contouring<CGAL::Parallel_if_available_tag>(domain, isovalue, points, triangles);
std::cout << "Output #vertices: " << points.size() << std::endl;
std::cout << "Output #triangles: " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring_inrimage.off", points, triangles);
}
int main(int argc, char* argv[])
{
const std::string fname = (argc > 1) ? argv[1] : CGAL::data_file_path("images/skull_2.9.inr");
const FT isovalue = (argc > 2) ? std::stod(argv[2]) : - 2.9;
// load volumetric image from a file
if(!image.read(fname))
{
std::cerr << "Error: Cannot read image file " << fname << std::endl;
return EXIT_FAILURE;
}
// convert image to a Cartesian grid
Grid grid;
Values values { grid }; // 'values' keeps a reference to the grid
if(!IS::IO::convert_image_to_grid(image, grid, values))
{
std::cerr << "Error: Cannot convert image to Cartesian grid" << std::endl;
return EXIT_FAILURE;
}
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
run_marching_cubes(grid, isovalue, values);
run_dual_contouring(grid, isovalue, values);
return EXIT_SUCCESS;
}
bool read(const char *file)
Class template for a gradient that is calculated using finite differences.
Definition: Finite_difference_gradient_3.h:39
std::string data_file_path(const std::string &filename)

Figure 61.8 Results of the Topologically Correct Marching Cubes algorithm for different isovalues (1, 2, and 2.9) on the skull model.


Offset Mesh

The following example illustrates how to generate a mesh approximating a signed offset to an input closed surface mesh. The input mesh is stored into an AABB_tree data structure to provide fast distance queries. Via the Side_of_triangle_mesh functor, the sign of the distance field is made negative inside the mesh.


File Isosurfacing_3/contouring_mesh_offset.cpp

#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/Isosurfacing_3/Cartesian_grid_3.h>
#include <CGAL/Isosurfacing_3/dual_contouring_3.h>
#include <CGAL/Isosurfacing_3/Dual_contouring_domain_3.h>
#include <CGAL/Isosurfacing_3/Value_function_3.h>
#include <CGAL/Isosurfacing_3/Finite_difference_gradient_3.h>
#include <CGAL/Isosurfacing_3/marching_cubes_3.h>
#include <CGAL/Isosurfacing_3/Marching_cubes_domain_3.h>
#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_traits_3.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/Polygon_mesh_processing/bbox.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <CGAL/boost/graph/IO/polygon_mesh_io.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <iostream>
#include <string>
#include <vector>
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point>;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
struct Offset_oracle
{
private:
const bool is_closed;
const Tree tree;
public:
Offset_oracle(const Mesh& mesh)
: is_closed(CGAL::is_closed(mesh)), tree(mesh.faces_begin(), mesh.faces_end(), mesh), sotm(mesh)
{ }
FT distance(const Point& p) const
{
const Point cp = tree.closest_point(p);
FT d = sqrt((p - cp).squared_length());
if(is_closed && sotm(p) == (CGAL::ON_BOUNDED_SIDE))
d *= -1.0;
return d;
}
};
void run_marching_cubes(const Grid& grid,
const FT offset_value,
const Offset_oracle& offset_oracle)
{
std::cout << "\n ---- " << std::endl;
std::cout << "Running Marching Cubes with offset value = " << offset_value << std::endl;
// fill up values
auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); };
Values values { mesh_distance, grid };
Domain domain { grid, values };
Point_range points;
Polygon_range triangles;
// run marching cubes
std::cout << "Running Marching Cubes with isovalue = " << offset_value << std::endl;
CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_if_available_tag>(domain, offset_value, points, triangles);
std::cout << "Output #vertices (MC): " << points.size() << std::endl;
std::cout << "Output #triangles (MC): " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("marching_cubes_offsets.off", points, triangles);
}
void run_dual_contouring(const Grid& grid,
const FT offset_value,
const Offset_oracle& offset_oracle)
{
std::cout << "\n ---- " << std::endl;
std::cout << "Running Dual Contouring with offset value = " << offset_value << std::endl;
// fill up values and gradients
auto mesh_distance = [&offset_oracle](const Point& p) { return offset_oracle.distance(p); };
const FT step = CGAL::approximate_sqrt(grid.spacing().squared_length()) * 0.01;
Values values { mesh_distance, grid };
Gradients gradients { values, step };
Domain domain { grid, values, gradients };
// output containers
Point_range points;
Polygon_range triangles;
// run dual contouring
std::cout << "Running Dual Contouring with isovalue = " << offset_value << std::endl;
CGAL::Isosurfacing::dual_contouring<CGAL::Parallel_if_available_tag>(
domain, offset_value, points, triangles,
CGAL::parameters::do_not_triangulate_faces(true));
std::cout << "Output #vertices (DC): " << points.size() << std::endl;
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
CGAL::IO::write_polygon_soup("dual_contouring_mesh_offset.off", points, triangles);
}
int main(int argc, char** argv)
{
const std::string filename = (argc > 1) ? argv[1] : CGAL::data_file_path("meshes/cross.off");
const FT offset_value = (argc > 2) ? std::stod(argv[2]) : 0.2;
if(offset_value < 0)
{
std::cerr << "Offset value must be positive" << std::endl;
return EXIT_FAILURE;
}
Mesh mesh;
if(!CGAL::IO::read_polygon_mesh(filename, mesh) || is_empty(mesh))
{
std::cerr << "Could not read input mesh" << std::endl;
return EXIT_FAILURE;
}
if(CGAL::is_closed(mesh))
std::cout << "Input mesh is closed - using signed distance offset" << std::endl;
else
std::cout << "Input mesh is not closed - using unsigned distance offset" << std::endl;
// construct loose bounding box from input mesh
const FT diag_length = sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
const FT loose_offset = offset_value + 0.1 * diag_length;
Vector aabb_increase_vec = Vector(loose_offset, loose_offset, loose_offset);
bbox += (Point(bbox.xmax(), bbox.ymax(), bbox.zmax()) + aabb_increase_vec).bbox();
bbox += (Point(bbox.xmin(), bbox.ymin(), bbox.zmin()) - aabb_increase_vec).bbox();
const int nv = 15;
Grid grid { bbox, CGAL::make_array<std::size_t>(nv, nv, nv) };
std::cout << "Span: " << grid.span() << std::endl;
std::cout << "Cell dimensions: " << grid.spacing()[0] << " " << grid.spacing()[1] << " " << grid.spacing()[2] << std::endl;
std::cout << "Cell #: " << grid.xdim() << ", " << grid.ydim() << ", " << grid.zdim() << std::endl;
Offset_oracle offset_oracle(mesh);
run_marching_cubes(grid, offset_value, offset_oracle);
run_dual_contouring(grid, offset_value, offset_oracle);
return EXIT_SUCCESS;
}
NT square(const NT &x)
bool is_empty(const FaceGraph &g)
bool is_closed(const FaceGraph &g)
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
ON_BOUNDED_SIDE

Design and Implementation History

The development of this package started during the 2022 Google Summer of Code, with the contribution of Julian Stahl, mentored by Daniel Zint and Pierre Alliez, providing a first implementation of Marching Cubes, Topologically Correct Marching Cubes, and Dual Contouring. Marching Cubes tables were provided by Roberto Grosso (FAU Erlangen-Nürnberg). Mael Rouxel-Labbé worked on improving the initial Dual Contouring implementation, and on the first complete version of the package.