User Manual

Author

Julian Stahl and Daniel Zint

Figure 60.1 Different isovalues of a bunny.

Introduction

This package provides methods to compute a surface mesh representing an isosurface of a 3-dimensional scalar field. An isosurface is defined as the surface on which the value of this field is equal to a given constant, i.e. the isovalue. Depending on the isosurfacing method and the input data structure, the result is either a triangular, quadrilateral, or higher order polygonal indexed face set.

Algorithms

There are multiple algorithms to extract isosurfaces. This package contains Marching Cubes, topologically correct Marching Cubes, and Dual Contouring.

Marching Cubes (MC)

MC processes all cells of the input domain individually. Each cell corner gets a sign (+/-) to indicate if it is above or below the isovalue. A new vertex is created on every cell edge where the sign changes, i.e. the isosurface is intersected. The vertex position is computed via linear interpolation of the scalar values of the incident corners. Depending on the configuration of signs at the corners the resulting vertices are connected to form triangles within the cell.

MC can process any input data structure that consists of hexahedral cells. In case of a conforming grid, MC produces a triangle mesh that is manifold in most scenarios. If the mesh is manifold and the isosurface does not intersect the domain boundaries, the mesh is also watertight. Compared to other approaches the algorithm often generates more and many thin triangles with acute angles. MC does not preserve sharp features of the input data.

Figure 60.2 Some example cases of Marching Cubes.

Topologically correct Marching Cubes (TMC)

This algorithm is an extension to MC and provides additional guarantees for the output. It generates a mesh that is homeomorphic to the trilinear interpolant of the input function inside each cube. This means that the mesh can accurately represent small complex features. For example a tunnel of the isosurface within a single cell can be topologically resolved. Furthermore, the mesh is guaranteed to be manifold and watertight, as long as the isosurface does not intersect the domain boundaries.

This algorithm is based on the following paper:

Roberto Grosso: Construction of Topologically Correct and Manifold Isosurfaces. Computer Graphics Forum 35(5):187-196 - August 2016

Dual Contouring (DC)

DC creates one vertex in every cell that is intersected by the isosurface. Next, a face 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 mesh.

The classical DC method requires the gradient of the scalar field. The provided domain therefore has to implement the concept IsosurfacingDomainWithGradient. All default domain implementations do this but assume the gradient to be zero if it is not provided as an additional parameter. Thus, for using the classical DC, the gradient has to be defined by the user. Alternatively, there are also some default gradient functions implemented, such as Finite_difference_gradient and Explicit_cartesian_grid_gradient.

Different versions of DC compute the vertex positions differently. Therefore, the vertex positioning is configurable with an optional parameter passed to the function. Some of them do not require the gradient and therefore even work with the zero gradient.

Dual Contouring works on any domain but does not guarantee a manifold or watertight mesh. It creates less faces than Marching Cubes. Another advantage of DC over MC is the ability to represent sharp edges.

Figure 60.3 Isosurface of the IWP function generated by Dual Contouring.

Comparison

Algorithm	Domains	Faces	Manifold	Watertight*	Topologically correct
MC	Hexahedral	Triangles	no	no	no
TMC	Hexahedral	Triangles	yes	yes	yes
DC	All	Polygons	no	no	no

(* assuming the isosurface does not leave the given bounding box of the domain)

Figure 60.4 Comparison between a cube generated by Dual Contouring (left) and Marching Cubes (right).

Interface

Each algorithm is represented by a single functions. The function signature is the same for all algorithms:

template <typename Concurrency_tag = Sequential_tag, class Domain_, class PointRange, class PolygonRange>
void marching_cubes(const Domain_& domain, const typename Domain_::FT iso_value,
                    PointRange& points, PolygonRange& polygons);

The input is provided in the form of a domain (see Domains).

The iso_value parameter describes the grid value the isosurface should represent.

The output is in the form of an indexed face set that is written to the two collections points and polygons. The vertex positions are stored as Point_3 in points. Each face in polygons is a list of indices pointing into the points collection. Depending on the algorithm, the indexed face set may either store a polygon soup or a topological mesh.

Algorithms can run sequentially on one CPU core or in parallel. The Concurrency_tag is used to specify how the algorithm is executed and is either Sequential_tag or Parallel_tag. To enable parallelism, CGAL needs to be linked with the Intel TBB library. If the parallel version is not availible the sequential version will always be used as a fallback.

Domains

A domain is an object that provides functions to access the input data, its geometry, topology, and optionally its gradient. There are some predefined domain classes that wrap the input data and provide a generalized interface for the algorithm. Users can also define new domains by implementing the Isosurfacing_domain or IsosurfacingDomainWithGradient concept.

Implicit cartesian grid domain

The Implicit_cartesian_grid_domain represents the input function in an implicit form without storing any values. It takes a functor or lambda that computes the value of the function from the position of a vertex as parameter. Additionally, the bounding box and spacing between grid points have to be specified.

Explicit cartesian grid domain

The Explicit_cartesian_grid_domain only takes a Cartesian_grid_3 as parameter. It represents a uniform grid of values that are either computed by the user or read from an Image_3. The constructor of Cartesian_grid_3 needs the number of grid points in each dimension and the bounding box. The values are read and written with value(x, y, z) where x, y, and z are the coordinates of a grid point. Alternatively, all required data can be copied from an Image_3.

Performance

Figure 60.5 Scaling of Marching Cubes with more threads.

Examples

Implicit sphere

The following example shows the usage of the Marching Cubes algorithm to extract an isosurface. The domain is an Implicit_domain that describes a sphere by the distance to its origin as an implicit function.

File Isosurfacing_3/marching_cubes_implicit_sphere.cpp

#include <CGAL/Bbox_3.h>
#include <CGAL/Implicit_cartesian_grid_domain.h>
#include <CGAL/Marching_cubes_3.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/boost/graph/IO/OFF.h>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::Vector_3 Vector;
typedef typename Kernel::Point_3 Point;

typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;

int main() {
    const CGAL::Bbox_3 bbox{-1, -1, -1, 1, 1, 1};
    const Vector spacing(0.04f, 0.04f, 0.04f);

    auto sphere_function = [&](const Point& p) { return std::sqrt(p.x() * p.x() + p.y() * p.y() + p.z() * p.z()); };

    // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
    auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(bbox, spacing, sphere_function);

    // prepare collections for the result
    Point_range points;
    Polygon_range polygons;

    // execute marching cubes with an isovalue of 0.8
    CGAL::Isosurfacing::marching_cubes(domain, 0.8f, points, polygons);

    // save the result in the OFF format
    CGAL::IO::write_OFF("result.off", points, polygons);
}

Cartesian_grid_3 cube comparison

The following example compares all provided algorithms to extract an isosurface. The domain is an Cartesian_grid_domain that describes a cube by storing the manhattan distance to its origin in a Cartesian_grid_3.

File Isosurfacing_3/all_cartesian_cube.cpp

#include <CGAL/Cartesian_grid_3.h>
#include <CGAL/Dual_contouring_3.h>
#include <CGAL/Explicit_cartesian_grid_domain.h>
#include <CGAL/Marching_cubes_3.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/boost/graph/IO/OFF.h>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;

typedef CGAL::Cartesian_grid_3<Kernel> Grid;

typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;

FT sign(FT value) {
    return (value > 0) - (value < 0);
}

int main() {
    // create a cartesian grid with 100^3 grid points and the bounding box [-1, 1]^3
    const CGAL::Bbox_3 bbox(-1, -1, -1, 1, 1, 1);
    std::shared_ptr<Grid> grid = std::make_shared<Grid>(7, 7, 7, bbox);

    // calculate the value at all grid points
    for (std::size_t x = 0; x < grid->xdim(); x++) {
        for (std::size_t y = 0; y < grid->ydim(); y++) {
            for (std::size_t z = 0; z < grid->zdim(); z++) {

                const FT pos_x = x * grid->get_spacing()[0] + bbox.xmin();
                const FT pos_y = y * grid->get_spacing()[1] + bbox.ymin();
                const FT pos_z = z * grid->get_spacing()[2] + bbox.zmin();

                // manhattan distance to the origin
                grid->value(x, y, z) = std::max({std::abs(pos_x), std::abs(pos_y), std::abs(pos_z)});
            }
        }
    }

    auto cube_gradient = [](const Point& p) {
        // the normal depends on the side of the cube
        const FT max_value = std::max({std::abs(p.x()), std::abs(p.y()), std::abs(p.z())});

        Vector g(0, 0, 0);
        if (max_value == std::abs(p.x())) {
            g += Vector(sign(p.x()), 0, 0);
        }
        if (max_value == std::abs(p.y())) {
            g += Vector(0, sign(p.y()), 0);
        }
        if (max_value == std::abs(p.z())) {
            g += Vector(0, 0, sign(p.z()));
        }
        const FT length_sq = g.squared_length();
        if (length_sq > 0.00001) {
            g /= CGAL::approximate_sqrt(length_sq);
        }
        return g;
    };

    // create a domain from the grid
    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid, cube_gradient);

    // prepare collections for the results
    Point_range points_mc, points_dc;
    Polygon_range polygons_mc, polygons_dc;

    // execute topologically correct marching cubes and dual contouring with an isovalue of 0.8
    CGAL::Isosurfacing::marching_cubes(domain, 0.88, points_mc, polygons_mc);
    CGAL::Isosurfacing::dual_contouring(domain, 0.88, points_dc, polygons_dc);

    // save the results in the OFF format
    CGAL::IO::write_OFF("result_mc.off", points_mc, polygons_mc);
    CGAL::IO::write_OFF("result_dc.off", points_dc, polygons_dc);
}

Image_3

The following example shows how to load data from an Image_3.

File Isosurfacing_3/marching_cubes_inrimage.cpp

#include <CGAL/Cartesian_grid_3.h>
#include <CGAL/Explicit_cartesian_grid_domain.h>
#include <CGAL/Marching_cubes_3.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/boost/graph/IO/OFF.h>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::Point_3 Point;

typedef CGAL::Cartesian_grid_3<Kernel> Grid;

typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;

int main() {

    const std::string fname = CGAL::data_file_path("images/skull_2.9.inr");

    // load the image
    CGAL::Image_3 image;
    if (!image.read(fname)) {
        std::cerr << "Error: Cannot read file " << fname << std::endl;
        return EXIT_FAILURE;
    }

    // convert it to a cartesian grid
    std::shared_ptr<Grid> grid = std::make_shared<Grid>(image);

    // create a domain from the grid
    auto domain = CGAL::Isosurfacing::create_explicit_cartesian_grid_domain<Kernel>(grid);

    // prepare collections for the result
    Point_range points;
    Polygon_range polygons;

    // execute marching cubes with an isovalue of 2.9
    CGAL::Isosurfacing::marching_cubes(domain, 2.9, points, polygons);

    // save the result in the OFF format
    CGAL::IO::write_OFF("result.off", points, polygons);
}

Offset mesh

The following example shows how to compute an offset mesh. The original mesh is passed to an AABB_tree to allow fast distance queries. With the use of Side_of_triangle_mesh the sign of the distance function is flipped inside the mesh.

File Isosurfacing_3/dual_contouring_mesh_offset.cpp

#include <CGAL/AABB_face_graph_triangle_primitive.h>
#include <CGAL/AABB_traits.h>
#include <CGAL/AABB_tree.h>
#include <CGAL/Bbox_3.h>
#include <CGAL/Dual_contouring_3.h>
#include <CGAL/Implicit_cartesian_grid_domain.h>
#include <CGAL/Side_of_triangle_mesh.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Surface_mesh.h>
#include <CGAL/boost/graph/IO/OFF.h>

#include <iostream>

typedef CGAL::Simple_cartesian<double> Kernel;
typedef typename Kernel::FT FT;
typedef typename Kernel::Point_3 Point;
typedef typename Kernel::Vector_3 Vector;

typedef CGAL::Surface_mesh<Point> Mesh;

typedef CGAL::AABB_face_graph_triangle_primitive<Mesh> Primitive;
typedef CGAL::AABB_traits<Kernel, Primitive> Traits;
typedef CGAL::AABB_tree<Traits> Tree;

typedef std::vector<Point> Point_range;
typedef std::vector<std::vector<std::size_t>> Polygon_range;

int main() {
    const std::string input_name = CGAL::data_file_path("meshes/cross.off");
    const Vector grid_spacing(0.1, 0.1, 0.1);
    const FT offset_value = 0.2;

    Mesh mesh_input;
    if (!CGAL::IO::read_OFF(input_name, mesh_input)) {
        std::cerr << "Could not read mesh" << std::endl;
        exit(-1);
    }

    // compute bounding box
    CGAL::Bbox_3 aabb_grid = CGAL::Polygon_mesh_processing::bbox(mesh_input);
    Vector aabb_increase_vec = Vector(offset_value + 0.01, offset_value + 0.01, offset_value + 0.01);
    aabb_grid += (Point(aabb_grid.xmax(), aabb_grid.ymax(), aabb_grid.zmax()) + aabb_increase_vec).bbox();
    aabb_grid += (Point(aabb_grid.xmin(), aabb_grid.ymin(), aabb_grid.zmin()) - aabb_increase_vec).bbox();

    // construct AABB tree
    Tree tree(mesh_input.faces_begin(), mesh_input.faces_end(), mesh_input);

    CGAL::Side_of_triangle_mesh<Mesh, CGAL::GetGeomTraits<Mesh>::type> sotm(mesh_input);


    auto mesh_distance = [&tree](const Point& p) {
        const Point& x = tree.closest_point(p);
        return std::sqrt((p - x).squared_length());
    };

    auto mesh_normal = [&tree](const Point& p) {
        const Point& x = tree.closest_point(p);
        const Vector n = p - x;
        return n / std::sqrt(n.squared_length());
    };

    // create a domain with bounding box [-1, 1]^3 and grid spacing 0.02
    auto domain = CGAL::Isosurfacing::create_implicit_cartesian_grid_domain<Kernel>(aabb_grid, grid_spacing,
                                                                                    mesh_distance, mesh_normal);

    Point_range points;
    Polygon_range polygons;

    CGAL::Isosurfacing::dual_contouring(domain, offset_value, points, polygons);

    CGAL::IO::write_OFF("result.off", points, polygons);
}