The aim of this document is to provide a documentation of the methods implemented by this library. We limit the discussion to
Note: In case math rendering looks odd, you might want to read this page in raw mode.
Volumentric data is commonly found in many scientific, engineering, and medical applications. Such volumentric data can be efficiently encoded using signed distance fields (SDFs). A SDF is a vector valued scalar function,
Its useful to note that the tuple
Still, this property gives raise to efficient ray maching schemes for visualizing SDF volumes. Another useful property of SDFs: the gradient
For many primitive shapes in
Isosurface extraction is the task of finding a suitable tesselation of the boundary (or any constant offset value) of a SDF. The methods considered in this library rely a regular SDF sampling volume from which the resulting mesh is generated. The two schemes that dominate the isosurface extraction field differ in the way they generate the tesselated topology. The following table lists the differences:
Intersecting grid element | |||
Method | Edge | Face | Voxel |
---|---|---|---|
Primal | Vertex | Edge | Face(s) |
Dual | Face | Edge | Vertex |
The table above lists, for each type of grid element that intersects the boundary of the SDF surface, which topological element is generated using either primary or dual approaches. For example, first row, second column means that in primal methods a vertex is created for each sampling edge that crosses the SDF boundary. See [1] for a more elaborate discussion.
For detailed discussions we need to define the set of coordinate systems involved. These coordinate systems match the implementation of the library.
The following image shows a sampling grid of (3,3,3)
(blue) points. The sampling coordinates are with respect to the data coordinate system (x,y,z). The sampling points are indexed by a grid coordinate system (i,j,k) that has its origin at the mininium sampling point. Each voxel in the grid is indexed by the integer grid coordinates of its minimum (i,j,k) corner. Note that the next voxel in the i-direction of voxel (i,j,k) is (i+1,j,k).
We index edges by its source voxel index plus a label e
Having three (forward) edges per voxel index allows us to easily enumerate all edges without duplicates and without missing any edges. Note, at the outer positive border faces we get a set of invalid edges (for example (2,0,2,0)
is invalid, while (2,0,2,1)
is valid).
This library implements a generic dual isosurface extraction method based. This blueprint allows the user to set different behavioral aspects to implement varios approaches (or hybrids thereof) proposed in literature.
Given a SDF and grid defining the sampling locations, the basic dual isosurface algorithm works as follows
- Active edges: For each edge in the sampling grid, determine if it intersects the boundary of the SDF. We call those edges with intersections active edges.
- Edge intersection: For each active edge find the intersection point with the boundary of the surface along the edge.
- Vertex placement: For each grid (active) voxel with at least one active edge, determine a single vertex location.
- Face generation: For each active edge create a quadliteral connecting the vertices of the four active voxels sharing this active edge.
See [3],[1] for more information.
The library implements this recipe in vectorized form. That is, all steps of the algorithms are capable to work with multiple elements at once. This allows for a better resource usage and generally speeds up algorithmic runtime dramatically. It is also the reason that you will hardly find for-loops sprinkled all over the code.
The recipe above gives raise to different behavioral aspects that are implemented as exchangable modules:
- edge strategies: determines the specific location of the surface intersection along an active edge.
- vertex strategies: determines the way the vertex position is calculated for all active voxels.
Edge strategies implement different methods to determine the edge/surface crossing. The following strategies are implemented
This method determines the intersection by finding the root of a linear equation guided by the SDF values at the two edge endpoints. This is most commonly found method in literature. It makes the following two assumptions
- Surface Smoothness: the SDF is assumed to be smooth
- Small edges: the edge lengths are supposed to be small compared to size of the shape of the SDF
Together, this two assumptions lead to linearity of the surface close to edges. Hence, modelling the surface boundary using a linear equation leads to accurate estimations. Also, in case you do not have access to analytic SDFs (e.g discretized volume of SDF values) it is the best you can do.
If the assumptions of the linear strategy are wrong, this leads to misplaced intersections that affect the quality of the final mesh. If you happen to have access to an analytic SDF, you might do better: We drop the assumption of surface linearity and instead find the root of the SDF along the edge iteratively. One algorithm with quadric convergence is Newton's method (it requires access to the gradient of the SDF). For our usecase (vector valued scalar function) we need a variant of it, the so called directional Newton method.
The bisection method is useful when a) linearity is not given, b) you have access to an analytic SDF and c) the gradient does not convey information along the edge direction (e.g. for some points in the SDF of a box).
Here is a diagram comparing the edge strategies on a spherical cross-section.
One notices, that for the linear estimator only the smaller edge seems to yield an accurate fit. For the larger edge, the main assumptions of the linear method break and hence the estimated root is off. Newtons method as well as the bisection method do not expose this issue at the cost of additional computational steps.
Shown below, is a similar plot for the cross section of a box.
The linear method fails for both edges because its main assumptions are violated. For Newton's method, the intersection for only one of the edges is computed correctly. No intersection is found for the other edge, since the gradient is orthogonal to the edge direction (no information along the edge dir). Only the bisection method is capable for producing an accurate result for both cases.
These strategies determine the final vertex locations active voxels.
The simplest strategy that places the vertex in the center the voxel. The resulting meshes are Mincraft-like box worlds.
Naive SurfaceNets [1] determine the the voxel's vertex location as the centroid of all 12 active edge intersection points. This strategy works well in practice but has the following downsides: a) it cannot reproduce sharp features like corners (like primal Marching Cubes [4]) and b) its prone to shrink the surface in regions of sharp features.
Dual Contouring [5] attempts to restore sharp features by determining the vertex position as the point that has the maximum support from all active edges of the voxel. Consider the plane formed by an active edge intersection point
Additionally, many edge configurations convey too little information to uniquely determine the location (as an example consider a xy-plane intersecting the four k-directed voxel edges). These configurations lead to an underdetermined system of linear equations. From a geometrical standpoint only the intersection of at least three planes (a corner) will uniquely determine the position of the vertex. In terms of the null-space, we can see that in the case of a single (independent) plane we end up with two dimensional null-space (all points in the plane). For two independent equations, we end up with one dimensional null-space (all points along the edge). For more than three-planes we resort to a least squares solution that minimizes
To deal null-spaces of dimension > 0, we add additional equations (multi objective linear least squares) that directly encode a preferred vertex location. These take the following
The following plot compares the three vertex placement strategies using the SDF of a union of two (offsetted) axis aligned boxes given by
boxes = sdftoolbox.sdfs.Union([
sdftoolbox.sdfs.Box().transform(trans=(0.5, 0.5, 0.5)),
sdftoolbox.sdfs.Box(),
]).transform(trans=(-0.25, -0.25, -0.25))
We use a low resolution grid of resolution 10x10x10
.
![](/cheind/sdftoolbox/raw/main/doc/vertex_strategies_aligned_box.gif)
Since the boxes are aligned with the world coordinate system, the midpoint strategy generates visually pleasing reconstruction. However, due to the placement of the vertices in the voxel centers, the resulting surface model has too little volume. Similarily, the naive SurfaceNets variant deforms the shape of the model due to a contraction of the average. This causes a loss of sharp features and a shape volume that is too small. Only the Dual Contouring strategy is capable of reconstructing sharp features and placing the vertices at locations that give rise to a volumetric matching reconstruction.
Shown below are the reconstructions of the same boxes, but this time rotated around the axis (1,1,1)
by 45°.
![](/cheind/sdftoolbox/raw/main/doc/vertex_strategies_rot_box.gif)
This time the midpoint placement strategy fails to capture the shape of the object (like aliasing effects in rendering). The naive SurfaceNets method has the same issues as in the aligned case (non sharp features, volume contraction). The Dual Contouring method manages to capture the shape correctly but induces a few non-manifold vertices (hard to see from the plot).
All plots are generated by doc_plots.py.
- [1] mikolalysenko. https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/
- [2] Inigio Quilez's. https://iquilezles.org/articles/distfunctions/
- [3] Gibson, Sarah FF. "Constrained elastic surface nets: Generating smooth surfaces from binary segmented data." International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, Berlin, Heidelberg, 1998.
- [4]Lorensen, William E., and Harvey E. Cline. "Marching cubes: A high resolution 3D surface construction algorithm." ACM siggraph computer graphics 21.4 (1987): 163-169.
- [5] Ju, Tao, et al. "Dual contouring of hermite data." Proceedings of the 29th annual conference on Computer graphics and interactive techniques. 2002.