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.
# Boilerplate code
from pygel3d import hmesh, jupyter_display as jd
import numpy as np
from math import cos, sin
jd.set_export_mode(True)
array = np.array