Isosurface Countouring

This notebook shows a very simple example of isosurface contouring using GEL. An isosurface (aka level set) is a set of points, $S$ in $\mathbb{R}^3$ that have the same value, $\tau$, for a given function: $f: \mathbb{R}^3 \rightarrow \mathbb R$. In other words, $S = \{ \mathbf{x} \in \mathbb{R}^3 : f(\mathbf{x}) = \tau \}$. Of course, it is easy to find examples where $S$ is not a surface, e.g. $f(\mathbf x) = \tau$. However, if we assume that the gradient of $f$ does not vanish, then $S$ is, in fact, a surface.

By isosurface contouring we understand finding a polygonal approximation of $S$. There is a number of ways of doing that, the most famous being the _Marching Cubes_ method due to Lorensen and Cline. More recently, the _Dual Contouring_ method was proposed by Ju et al. who based their work (in part) on the idea of _Surface Nets_ due to Frisken. Common to these methods is that space is divided into a grid of (often cuboid) cells and $f$ is usually only known at the grid points. The grid points are often called voxels although, confusingly, the word voxel is also frequently used to denote a rectangular cell in 3D. It is perhaps best to think of a voxel as the generalization of pixel to 3D and pixels can also be seen both as little squares and point samples of a continuous image.

For some illustrations and comparisons of dual contouring and marching cubes, please see my old page on isosurface polygonization. Read on for a presentation of the method implemented in GEL.

The GEL library implements a variation of the Dual Contouring method. Dual Contouring is arguably simpler than Marching Cubes and it often produces better triangle meshes in the sense that there are fewer triangles with close to zero area. That being said, it is not as clear how to place vertices in Dual Contouring. The scheme used in GEL is the following: each vertex is initially placed at the barycenter of the eight grid points (i.e. voxels) that are at the corners of the cube to which it belongs. The value of $f$ at the vertex is simply the average of those eight voxels. We now consider the edges to the eight corners. For each such edge where the values at the end points are on opposite sides of $\tau$, we find the intersection (i.e. the point on the edge where linear interpolation would produce $\tau$ exactly), and the final vertex position is the average of these.

For vertices which lie on the boundary of the domain, the initial position is the barycenter of the grid points which belong the voxel grid, and we only consider intersections of the edges to those grid points. This ensures that the output vertices lie within the voxel grid.

The API for isosurface contouring is as follows:

def volumetric_isocontouring(data, bbox_min = None, bbox_max = None, tau=0.0, make_triangles=True, high_is_inside=True)

The first argument, data contains the volume which we pass as a numpy 3D array. The next two arguments define the bounding box through the coordinates of its smallest and biggest corner. If not provided, the bounding box is set to the extents of the volumetric grid. tau is the isovalue, make_triangles (defaults to True) indicates whether the raw quads should be output or the mesh triangulated. Finally, high_is_inside indicates whether values greater than the isovalue are considered to be inside the contoured shape. This makes a difference if the contour intersects the boundary of the volume since it tells us how to close the surface.

Creating a volume

A 3D grid of voxels is often called a volume. Below, we use a 3D numpy array as our representation of the volume and fill it with the values of the trigonometric approximation of the gyroid

Contouring the isosurface

Most settings are at default, but high_is_inside is set to false meaning that values smaller than the isovalue are inside and values greater are outside. The isovalue is left at 0.

The method used for isocontouring is the dual contouring method which is perhaps slightly less well-known than marching cubes but tends to produce nicer triangles. The price for the nicer triangles is that vertex placement is a bit less straight forward.

A Plane

We definitely want the polygonizer to be able to reconstruct a planar surface exactly. We also want the polygonizer to behave reasonably on the boundary of the domain. Both are tested below where we compute an isocontour from an implicit representation of a plane. For this test, we set make_triangles=False which only omits the triangulation step since the polygonizer produces quads initially.

A Distance Field

Finally, we create a distance field of the Stanford Bunny and use the isocontouring to produce a new mesh. Initially, the mesh is loaded, and we obtain the bounding box.

Next, we sample the distance field at the points of a regular voxel grid. We adjust the bounding box corner corresponding to the highest indexed (along all axes) grid point. We will pass the corners of the bounding box to the volumetric_isocontouring function. This makes it possible to map the vertices such that the isocontour has the same scale and position as the input model.

Finally, we extract the isocontour. We map the distance field onto the resulting mesh to show any deviation rom the distance field. Note that the deviations are small since the mesh was created as a polygonized isocontour of the same distance field, but they are not zero - especially not in high-curvature regions.

Discussion

There is a lot more to say about isosurface contouring. Methods have been developed which provide various guarantees and which contour analytic functions or irregular grids or hierarchical grids. However, this brief tutorial aims only to demonstrate the capabilities of GEL and the Pythin bindings in PyGEL.