CGAL 6.1 - 3D Isosurfacing

Loading...
Searching...
No Matches
User Manual

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

Figure 64.1 Generating a surface from a 3D gray level image using Marching Cubes (3D input image from 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, a 3D image, 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 are 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) : a simple and efficient method that generates a triangle mesh, with almost no guarantees.
- Topologically Correct Marching Cubes (TMC) : an extension to Marching Cubes that provides additional guarantees for the output.
- Dual Contouring (DC) : a method that generates a polygon mesh, with a focus on sharp features recovery.

## Marching Cubes (MC)

Marching Cubes (MC) [6] 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 using 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 64.2 illustrates the configurations in 2D. In 3D, there are no less than 33 configurations (not shown) [1].

Figure 64.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 or cap-shaped triangles.

Furthermore, Marching Cubes does not preserve the sharp features present in the isovalue of the input value field (see Figure 64.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 [4]. 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 64.3 Marching Cubes vs Topologically Correct Marching Cubes. Some vertex values can represent complex surfaces within a single cell: on the left, the values at the vertices of the cell have been interpolated within the cell and correspond to a tunnel-like surface. This surface is not correctly recovered when applying the cases of the Marching Cubes algorithm to the cell: 2 independent disk-like sheets are created (middle). However, by inserting new vertices, TMC correctly captures the topology within a single cell (right).

## Dual Contouring (DC)

Dual Contouring (DC) [5] 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 64.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 [3] : 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 using a named parameter.

## 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 | Quads** | no | no | no | yes (not guaranteed)

(* assuming the isosurface does not exit the specified bounding box of the input 3D domain)
(** which can be triangulated on-the-fly)

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

Figure 64.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 of this package:

- Marching Cubes: `CGAL::Isosurfacing::marching_cubes()`, using the named parameter: `use_topologically_correct_marching_cubes` set to `false`;
- Topologically Correct Marching Cubes: `CGAL::Isosurfacing::marching_cubes()`, using the named parameter: `use_topologically_correct_marching_cubes` set to `true`;
- Dual Contouring: `CGAL::Isosurfacing::dual_contouring()`.

All these 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,
...);
CGAL::Sequential_tag

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 (i.e., 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).

All isosurfacing algorithms can run either sequentially in one thread or in parallel using 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:

- Partition: this must be a class that describes the partition of the 3D volume into cells. The most basic example of such a class is `CGAL::Isosurfacing::Cartesian_grid_3`, but users can pass their own partition, provided it meets the requirements described by the concept `IsosurfacingPartition_3`.
- ValueField: this must be a class that provides the 3D value field at the vertices of the partition. A few classes are provided, such as `CGAL::Isosurfacing::Value_function_3` and `CGAL::Isosurfacing::Interpolated_discrete_values_3`. Users can pass their own value class, provided it meets the requirements described by the concept `IsosurfacingValueField_3`.
- GradientField: (`CGAL::Isosurfacing::Dual_contouring_domain_3` only) this must be a class that provides the gradient of the value field at the vertices of the partition. A few classes are provided by default, such as `CGAL::Isosurfacing::Finite_difference_gradient_3` and `CGAL::Isosurfacing::Interpolated_discrete_gradients_3`. Users can pass their own gradient class, provided it meets the requirements described by the concept `IsosurfacingGradientField_3`.
- EdgeIntersectionOracle: this must be a class that provides a function to compute the intersection between an edge and the isosurface. The default is linear interpolation for `CGAL::Isosurfacing::Marching_cubes_domain_3`, and a dichotomy for `CGAL::Isosurfacing::Dual_contouring_domain_3`. This parameter should be adjusted depending on how the value field is defined: there is for example no point incorporating a dichotomy in Dual Contouring if the value field is defined through linear interpolation. Users can pass their own edge intersection oracle, provided it meets the requirements described by the concept `IsosurfacingEdgeIntersectionOracle_3`.

# 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. The resulting triangle soup is converted to a triangle mesh, and is remeshed using the isotropic remeshing algorithm.

File Isosurfacing_3/marching_cubes.cpp

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh.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/Bbox_3.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <CGAL/Polygon_mesh_processing/remesh.h>
#include <cmath>
#include <iostream>
#include <vector>
using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
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>;
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
using Mesh = CGAL::Surface_mesh<Point>;
auto value_fn = [](const Point& p) -> FT
{
const FT& x = p.x(), y = p.y(), z = p.z();
// "Klein Bottle" - https://www-sop.inria.fr/galaad/surface/
return -1-4*y*z*z*x*x-2*y+6*z*z*x*x*y*y-16*x*z+16*x*z*y*y+3*x*x+7*y*y+11*z*z-11*z*z*z*z+z*z*z*z*z*z-3*x*x*x*x-7*y*y*y*y+x*x*x*x*x*x+y*y*y*y*y*y-14*z*z*x*x-18*z*z*y*y+3*z*z*z*z*x*x+3*z*z*z*z*y*y-10*x*x*y*y-4*y*y*y*z*z+3*z*z*x*x*x*x+3*z*z*y*y*y*y+16*x*x*x*z+3*x*x*x*x*y*y+3*x*x*y*y*y*y+4*x*x*y-12*z*z*y-2*x*x*x*x*y-4*x*x*y*y*y-2*z*z*z*z*y+16*x*z*z*z+12*y*y*y-2*y*y*y*y*y-32*x*z*y;
};
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.1;
const FT box_c = (argc > 2) ? std::abs(std::stod(argv[2])) : 5.;
const std::size_t grid_n = (argc > 3) ? std::stoi(argv[3]) : 50;
// create bounding box and grid
const CGAL::Bbox_3 bbox { -box_c, -box_c, -box_c, box_c, box_c, box_c };
Grid grid { bbox, CGAL::make_array<std::size_t>(grid_n, grid_n, grid_n) };
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
Values values { value_fn, grid };
// Below is equivalent to:
// Domain domain { grid, values };
Domain domain = CGAL::Isosurfacing::create_marching_cubes_domain_3(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 << "Soup #vertices: " << points.size() << std::endl;
std::cout << "Soup #triangles: " << triangles.size() << std::endl;
if(!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)) {
std::cerr << "Warning: the soup is not a 2-manifold surface, non-manifoldness?..." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Convert the soup to a triangle mesh..." << std::endl;
Mesh mesh;
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, mesh);
CGAL::IO::write_polygon_mesh("marching_cubes.off", mesh, CGAL::parameters::stream_precision(17));
// Let's remesh it to something nicer looking
std::cout << "Remeshing..." << std::endl;
const FT target_edge_length = box_c / 50;
CGAL::Polygon_mesh_processing::isotropic_remeshing(faces(mesh), target_edge_length, mesh,
CGAL::parameters::number_of_iterations(5)
.number_of_relaxation_steps(5));
CGAL::IO::write_polygon_mesh("marching_cubes-remeshed.off", mesh, CGAL::parameters::stream_precision(17));
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
CGAL::Bbox_3
CGAL::Exact_predicates_inexact_constructions_kernel
CGAL::Isosurfacing::Cartesian_grid_3
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
CGAL::Isosurfacing::Marching_cubes_domain_3
A domain that can be used with the Marching Cubes algorithm.
Definition: Marching_cubes_domain_3.h:50
CGAL::Isosurfacing::Value_function_3
The class Value_function_3 represents a field of scalars computed using a user-provided unary functio...
Definition: Value_function_3.h:40
Kernel::Point_3
Kernel::Vector_3
CGAL::Isosurfacing::create_marching_cubes_domain_3
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
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh
void polygon_soup_to_polygon_mesh(const PointRange &points, const PolygonRange &polygons, PolygonMesh &out, const NamedParameters_PS &np_ps=parameters::default_values(), const NamedParameters_PM &np_pm=parameters::default_values())
CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh
bool is_polygon_soup_a_polygon_mesh(const PolygonRange &polygons)
CGAL::Polygon_mesh_processing::isotropic_remeshing
void isotropic_remeshing(const FaceRange &faces, SizingFunction &sizing, PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
CGAL::IO::write_polygon_mesh
bool write_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
CGAL::Polygon_mesh_processing::bbox
CGAL::Bbox_3 bbox(const PolygonMesh &pmesh, const NamedParameters &np=parameters::default_values())
Kernel

Figure 64.6 Results of the Marching Cubes algorithm, and the final result after remeshing.

## 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/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/Gradient_function_3.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/IO/polygon_mesh_io.h>
#include <CGAL/Polygon_mesh_processing/polygon_soup_to_polygon_mesh.h>
#include <cmath>
#include <iostream>
#include <vector>
using Kernel = CGAL::Simple_cartesian<double>;
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>;
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Grid>;
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
using Mesh = CGAL::Surface_mesh<Point>;
// "Devil" - 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;
const FT box_c = (argc > 2) ? std::abs(std::stod(argv[2])) : 1.;
const std::size_t grid_n = (argc > 3) ? std::stoi(argv[3]) : 50;
// create bounding box and grid
const CGAL::Bbox_3 bbox { -box_c, -box_c, -box_c, box_c, box_c, box_c };
Grid grid { bbox, CGAL::make_array<std::size_t>(grid_n, grid_n, grid_n) };
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 & gradients
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 << "Soup #vertices: " << points.size() << std::endl;
std::cout << "Soup #triangles: " << triangles.size() << std::endl;
if(!CGAL::Polygon_mesh_processing::is_polygon_soup_a_polygon_mesh(triangles)) {
std::cerr << "Warning: the soup is not a 2-manifold surface, non-manifoldness?..." << std::endl;
return EXIT_FAILURE;
}
std::cout << "Create the mesh..." << std::endl;
Mesh mesh;
CGAL::Polygon_mesh_processing::polygon_soup_to_polygon_mesh(points, triangles, mesh);
CGAL::IO::write_polygon_mesh("dual_contouring.off", mesh, CGAL::parameters::stream_precision(17));
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
CGAL::Isosurfacing::Dual_contouring_domain_3
A domain that can be used as input in the Dual Contouring algorithm.
Definition: Dual_contouring_domain_3.h:52
CGAL::Isosurfacing::Gradient_function_3
The class Gradient_function_3 represents a field of vectors computed using a user-provided unary func...
Definition: Gradient_function_3.h:39
CGAL::Isosurfacing::create_dual_contouring_domain_3
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::Simple_cartesian

Figure 64.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).

## 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 Kernel = CGAL::Simple_cartesian<double>;
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 Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel, CGAL::Isosurfacing::Cache_vertex_locations>;
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Grid>;
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)
{
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
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)
{
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
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;
}
CGAL::IO::write_polygon_soup
bool write_polygon_soup(const std::string &fname, const PointRange &points, const PolygonRange &polygons, const NamedParameters &np=parameters::default_values())

## 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 Kernel = CGAL::Simple_cartesian<double>;
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>;
using Values = CGAL::Isosurfacing::Interpolated_discrete_values_3<Grid>;
using Gradients = CGAL::Isosurfacing::Interpolated_discrete_gradients_3<Grid>;
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;
}
CGAL::Isosurfacing::Interpolated_discrete_gradients_3
Class template for a gradient field that is computed using discrete values and interpolation.
Definition: Interpolated_discrete_gradients_3.h:40
CGAL::Isosurfacing::Interpolated_discrete_values_3
Class template for a field of values that are calculated using discrete values and interpolation.
Definition: Interpolated_discrete_values_3.h:41
CGAL::sqrt
NT sqrt(const NT &x)
CGAL::NULL_VECTOR
const CGAL::Null_vector NULL_VECTOR
CGAL::ORIGIN
const CGAL::Origin ORIGIN
CGAL::squared_distance
Kernel::FT squared_distance(Type1< Kernel > obj1, Type2< Kernel > obj2)
CGAL::Isosurfacing
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 Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
using Values = CGAL::Isosurfacing::Interpolated_discrete_values_3<Grid>;
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
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
CGAL::Image_3 image;
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;
}
CGAL::Image_3
CGAL::Image_3::read
bool read(const char *file)
CGAL::Isosurfacing::Finite_difference_gradient_3
Class template for a gradient that is calculated using finite differences.
Definition: Finite_difference_gradient_3.h:39

Figure 64.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 Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
using Point = typename Kernel::Point_3;
using Vector = typename Kernel::Vector_3;
using Mesh = CGAL::Surface_mesh<Point>;
using Grid = CGAL::Isosurfacing::Cartesian_grid_3<Kernel>;
using Values = CGAL::Isosurfacing::Value_function_3<Grid>;
using Gradients = CGAL::Isosurfacing::Finite_difference_gradient_3<Kernel>;
using Point_range = std::vector<Point>;
using Polygon_range = std::vector<std::vector<std::size_t> >;
struct Offset_oracle
{
using Primitive = CGAL::AABB_face_graph_triangle_primitive<Mesh>;
using Traits = CGAL::AABB_traits_3<Kernel, Primitive>;
using Tree = CGAL::AABB_tree<Traits>;
private:
const bool is_closed;
const Tree tree;
CGAL::Side_of_triangle_mesh<Mesh, Kernel> sotm;
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)
{
using Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Grid, Values>;
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)
{
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Grid, Values, Gradients>;
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
CGAL::Bbox_3 bbox = CGAL::Polygon_mesh_processing::bbox(mesh);
const FT diag_length = sqrt(CGAL::square(bbox.xmax() - bbox.xmin()) +
CGAL::square(bbox.ymax() - bbox.ymin()) +
CGAL::square(bbox.zmax() - bbox.zmin()));
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;
}
CGAL::AABB_face_graph_triangle_primitive
CGAL::AABB_traits_3
CGAL::AABB_tree
CGAL::Bbox_3::ymin
double ymin() const
CGAL::Bbox_3::xmax
double xmax() const
CGAL::Bbox_3::zmin
double zmin() const
CGAL::Bbox_3::zmax
double zmax() const
CGAL::Bbox_3::ymax
double ymax() const
CGAL::Bbox_3::xmin
double xmin() const
CGAL::Side_of_triangle_mesh
CGAL::square
NT square(const NT &x)
CGAL::is_empty
bool is_empty(const FaceGraph &g)
CGAL::is_closed
bool is_closed(const FaceGraph &g)
CGAL::IO::read_polygon_mesh
bool read_polygon_mesh(const std::string &fname, Graph &g, const NamedParameters &np=parameters::default_values())
CGAL::ON_BOUNDED_SIDE
ON_BOUNDED_SIDE
CGAL

## Contouring using Octrees

The following example shows the use of an octree as discretization instead of a Cartesian grid for both Dual Contouring and Marching Cubes. Note that in this configuration, all methods lose guarantees, and cracks can easily appear in the surface for Marching Cubes if the octree is not adapted to the surface being extracted.

File Isosurfacing_3/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/Bbox_3.h>
#include <CGAL/IO/polygon_soup_io.h>
#include <cmath>
#include <iostream>
#include <vector>
using Kernel = CGAL::Simple_cartesian<double>;
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> >;
using Octree = CGAL::Octree<Kernel, std::vector<Point> >;
using Values = CGAL::Isosurfacing::Value_function_3<Octree>;
using Gradients = CGAL::Isosurfacing::Gradient_function_3<Octree>;
using MC_Domain = CGAL::Isosurfacing::Marching_cubes_domain_3<Octree, Values>;
using Domain = CGAL::Isosurfacing::Dual_contouring_domain_3<Octree, Values, Gradients>;
namespace IS = CGAL::Isosurfacing;
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());
};
auto blobby_function = [](const Point& p) -> FT
{
return std::exp(-1.5 * ((p.x() - 0.2) * (p.x() - 0.2) + (p.y() - 0.2) * (p.y() - 0.2) + (p.z() - 0.2) * (p.z() - 0.2))) +
std::exp(-1.5 * ((p.x() + 0.2) * (p.x() + 0.2) + (p.y() + 0.2) * (p.y() + 0.2) + (p.z() + 0.2) * (p.z() + 0.2))) +
std::exp(-1.5 * ((p.x() - 0.4) * (p.x() - 0.4) + (p.y() + 0.4) * (p.y() + 0.4) + (p.z() - 0.4) * (p.z() - 0.4))) +
std::exp(-6 * ((p.x() - 0.1) * (p.x() - 0.1) + (p.y() - 0.1) * (p.y() - 0.1))) +
std::exp(-6 * ((p.y() + 0.1) * (p.y() + 0.1) + (p.z() + 0.1) * (p.z() + 0.1))) +
std::exp(-6 * ((p.x() + 0.1) * (p.x() + 0.1) + (p.z() - 0.1) * (p.z() - 0.1))) -
0.3;
};
auto blobby_gradient = [](const Point& p) -> Vector
{
const FT g1 = -3 * std::exp(-1.5 * ((p.x() - 0.2) * (p.x() - 0.2) + (p.y() - 0.2) * (p.y() - 0.2) + (p.z() - 0.2) * (p.z() - 0.2)));
const FT g2 = -3 * std::exp(-1.5 * ((p.x() + 0.2) * (p.x() + 0.2) + (p.y() + 0.2) * (p.y() + 0.2) + (p.z() + 0.2) * (p.z() + 0.2)));
const FT g3 = -3 * std::exp(-1.5 * ((p.x() - 0.4) * (p.x() - 0.4) + (p.y() + 0.4) * (p.y() + 0.4) + (p.z() - 0.4) * (p.z() - 0.4)));
const FT g4 = -12 * std::exp(-6 * ((p.x() - 0.1) * (p.x() - 0.1) + (p.y() - 0.1) * (p.y() - 0.1)));
const FT g5 = -12 * std::exp(-6 * ((p.y() + 0.1) * (p.y() + 0.1) + (p.z() + 0.1) * (p.z() + 0.1)));
const FT g6 = -12 * std::exp(-6 * ((p.x() + 0.1) * (p.x() + 0.1) + (p.z() - 0.1) * (p.z() - 0.1)));
return Vector(g1 * (p.x() - 0.2) + g2 * (p.x() + 0.2) + g3 * (p.x() - 0.4) + g4 * (p.x() - 0.1) + g6 * (p.x() + 0.1),
g1 * (p.y() - 0.2) + g2 * (p.y() + 0.2) + g3 * (p.y() + 0.4) + g4 * (p.y() - 0.1) + g5 * (p.y() + 0.1),
g1 * (p.z() - 0.2) + g2 * (p.z() + 0.2) + g3 * (p.z() - 0.4) + g5 * (p.z() + 0.1) + g6 * (p.z() - 0.1));
};
// This is a naive refinement that is adapted to the isosurface:
// This refines:
// - at the minimum till minimum depth
// - at the maximum till maximum depth
// - we split if the the isovalue goes through the voxel, i.e. if not all vertices of the cell
// are on the same side of the isosurface defined by a function
// It's not a great refinement technique because the surface can enter and leave a cell
// without involving the cell's vertex. In practice, that means a hole if at nearby adjacent
// cells the voxels did get refined and registered the surface.
struct Refine_around_isovalue
{
std::size_t min_depth_;
std::size_t max_depth_;
std::function<FT(const Point&)> function_;
FT isovalue_;
Refine_around_isovalue(std::size_t min_depth,
std::size_t max_depth,
std::function<FT(const Point&)> function,
FT isovalue)
: min_depth_(min_depth),
max_depth_(max_depth),
function_(function),
isovalue_(isovalue)
{}
bool operator()(const Octree::Node_index& ni, const Octree& octree) const
{
// Ensure minimum depth refinement
if (octree.depth(ni) < min_depth_)
return true;
// Stop refinement at maximum depth
if (octree.depth(ni) >= max_depth_)
return false;
// Get the bounding box of the node
auto bbox = octree.bbox(ni);
// Evaluate the function at the corners of the bounding box
std::array<FT, 8> corner_values;
int index = 0;
for (FT x : {bbox.xmin(), bbox.xmax()})
for (FT y : {bbox.ymin(), bbox.ymax()})
for (FT z : {bbox.zmin(), bbox.zmax()})
corner_values[index++] = function_(Point(x, y, z));
// Check if the function values cross the isovalue
bool has_positive = false, has_negative = false;
for (const auto& value : corner_values)
{
if (value > isovalue_)
has_positive = true;
if (value < isovalue_)
has_negative = true;
if (has_positive && has_negative)
return true; // Refine if the isosurface intersects the voxel
}
return false; // No refinement needed
}
};
// This is a refinement that is NOT adapted to the isosurface
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;
}
};
template <typename Splitter>
void run_DC_octree(const CGAL::Bbox_3 bbox,
const Splitter& split_predicate,
const std::function<FT(const Point&)> function,
const std::function<Vector(const Point&)> gradient,
const FT isovalue,
const std::string& name)
{
std::vector<Point> bbox_points { { bbox.xmin(), bbox.ymin(), bbox.zmin() },
{ bbox.xmax(), bbox.ymax(), bbox.zmax() } };
Octree octree(bbox_points);
octree.refine(split_predicate);
auto leaf_range = octree.traverse(CGAL::Orthtrees::Leaves_traversal<Octree>(octree));
std::size_t leaf_counter = std::distance(leaf_range.begin(), leaf_range.end());
std::cout << "octree has " << leaf_counter << " leaves" << std::endl;
// fill up values and gradients
Values values { function, octree };
Gradients gradients { gradient, octree };
Domain domain { octree, values, gradients };
// output containers
Point_range points;
Polygon_range triangles;
std::cout << "Running Dual Contouring with isovalue = " << isovalue << std::endl;
// run Dual Contouring
IS::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 (DC): " << points.size() << std::endl;
std::cout << "Output #triangles (DC): " << triangles.size() << std::endl;
std::ofstream oo("octree_DC_" + name + ".polylines.txt");
oo.precision(17);
octree.dump_to_polylines(oo);
CGAL::IO::write_polygon_soup("DC_" + name + ".off", points, triangles);
}
template <typename Splitter>
void run_MC_octree(const CGAL::Bbox_3 bbox,
const Splitter& split_predicate,
const std::function<FT(const Point&)> function,
const FT isovalue,
const std::string& name)
{
std::vector<Point> bbox_points { { bbox.xmin(), bbox.ymin(), bbox.zmin() },
{ bbox.xmax(), bbox.ymax(), bbox.zmax() } };
Octree octree(bbox_points);
octree.refine(split_predicate);
Values values { function, octree };
Point_range points;
Polygon_range triangles;
MC_Domain mcdomain { octree, values };
std::cout << "Running MC" << std::endl;
CGAL::Isosurfacing::marching_cubes<CGAL::Parallel_if_available_tag>(mcdomain, isovalue, points, triangles);
std::cout << "Output #vertices (MC): " << points.size() << std::endl;
std::cout << "Output #triangles (MC): " << triangles.size() << std::endl;
std::ofstream oo("octree_MC_" + name + ".polylines.txt");
oo.precision(17);
octree.dump_to_polylines(oo);
CGAL::IO::write_polygon_soup("MC_" + name + ".off", points, triangles);
}
// Whether you are using MC, TMC, or DC, there is no guarantee for an octree:
// it should behave well if your nodes are split with a uniform size around the surface,
// but it is sure to produce cracks if you have varying depths around the isosurface.
int main(int argc, char** argv)
{
const FT isovalue = (argc > 1) ? std::stod(argv[1]) : 0.3;
const CGAL::Bbox_3 bbox { -1., -1., -1., 1., 1., 1. };
Refine_one_eighth one_eight_splitter(3, 5);
run_DC_octree(bbox, one_eight_splitter, sphere_function, sphere_gradient, isovalue, "one_eight");
// This is
Refine_around_isovalue isovalue_splitter(3, 5, sphere_function, isovalue);
run_DC_octree(bbox, isovalue_splitter, sphere_function, sphere_gradient, isovalue, "sphere_adapted");
Refine_around_isovalue isvalue_splitter_2(3, 5, blobby_function, isovalue);
run_DC_octree(bbox, isvalue_splitter_2, blobby_function, blobby_gradient, isovalue, "blobby_adapted");
// to illustrate cracks
run_MC_octree(bbox, isovalue_splitter, sphere_function, isovalue, "adapted");
run_MC_octree(bbox, one_eight_splitter, sphere_function, isovalue, "one_eight");
std::cout << "Done" << std::endl;
return EXIT_SUCCESS;
}
CGAL::Orthtree
CGAL::Octree
Orthtree< Orthtree_traits_point< GeomTraits, PointRange, PointMap, cubic_nodes, 3 > > Octree
CGAL::Orthtrees::Leaves_traversal

Figure 64.9 Running Marching Cubes with an octree: left, the octree (blue wireframe) has been applied arbitrary refinement in one eighth of the domain and the resulting mesh exhibits cracks (boundaries) at the regions where octree levels are discontinuous. Right, the octree was constructed to be dense around the isosurface, and the surface is manifold.

# 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.