Problem: Reconstructing meshes with sharp features from signed distance data

A geometric problem that I've been interested in for some time is how to get an accurate 3D surface mesh (triangulation) back from a volumetric signed distance field. Since I'm not likely to tackle this in the near future, in the spirit of "open science" I thought I'd just post some thoughts, and see how the world responds.

Level sets in a nutshell

A common implicit surface representation is a signed distance field stored on a volumetric regular grid. At each gridpoint the sign indicates whether the gridpoint is inside or outside of the surface, and the value at the gridpoint is the distance to the nearest point on the surface. The interpolated zero isocontour represents the position of the surface implicitly. Although technically I believe level sets refer to any implicit function, signed distance is by the most commonly used, at least for graphics applications.

Mesh to level set conversion

Because one of the de facto standard surface representations in graphics is the triangle mesh, there are a fairly large number of different algorithms for converting such meshes into grid-based signed distance data efficiently and accurately. The simplest (though terribly inefficient) is to brute force it: compute the closest point on the triangle mesh for every grid point and use some form of raycasting/parity count to determine the signs. Doing the latter in a geometrically robust way can be a bit tricky, but is certainly doable. (Robert Bridson has a small code called mesh_query on his website for doing just that.) Since level sets are useful for various geometry processing operations, for fast proximity queries, and for representing deformable surfaces such as liquids, converting from standard triangle meshes to the level set representation is a useful operation. For a survey of signed distance fields and much better ways to compute them from meshes, see the article 3D Distance Fields: A Survey of Techniques and Applications.

Level set to mesh conversion - What are the issues?

Going from a level set representation back to a triangle mesh seems (to me) to be a much less well-studied problem. The simplest approach here is to assume a bi/tri-linear interpolant of the data within each grid cell, and then use marching cubes or some other standard surface extraction method to get a triangle mesh out. Unfortunately, that performs rather poorly. Sharp features present in the data are rounded away, joints between cells are visually apparent, and overall the reconstruction is of poor quality. This is highlighted for example in the teaser image (Figure 1) of the paper Feature Sensitive Surface Extraction from Volume Data. That paper also goes on to illustrate that regardless of the distance field resolution, the extracted surface never converges to the correct sharp features (Figure 6).

A first thought might be to use a higher order interpolant, say tri-cubic interpolation. Our good friend Wikipedia illustrates that this yields a smoother approximation of our distance field between cells, which is good, but it doesn't necessarily resolve the problem of capturing sharp features. It can also introduce overshoots that cause the reconstructed surface to diverge further from the data.

Another thought is to use methods designed with reconstructing sharp features in mind, such as Dual Contouring or Feature Sensitive Surface Extraction from Volume Data mentioned above. Unfortunately these methods require normal vectors (in addition to distance data) in order to find and extract sharp features. We could construct some approximate normals by taking finite difference gradients of our signed distance field, however in regions with sharp features where the grid resolution is comparatively low, those gradients are likely going to be of very poor quality.

There's been a fair amount of work on reconstructing sharp features from point set surfaces using moving least squares, but this representation assumes you have point samples directly on the surface, rather than distance values lying away from the surface. Perhaps these could be adapted to work with signed distances.

As a general observation, my impression is that people think of signed distance fields as somehow inherently smoothed or rounded, whereas this isn't necessarily the case. The interpolants (or reconstructions) most commonly used with level sets have this effect, but the underlying data often contains more information that I hope will allow us to do a much better job with a bit more effort.

What can we do?

For the purposes of this problem, I am going to assume we have essentially exact signed distance data to work with, not corrupted with noise. This tells us a few things that we might potentially exploit. One is that our grid-based distance values are a sampling from a continuous signed distance function with unit length gradients everywhere, other than at local extrema. ie. if phi is our (true) distance field, then ||grad phi|| = 1. (eg. in 1D, the distance function always has slope +/-1) Perhaps we can use this to determine improved normal estimates, and then use Dual Contouring or one of its descendants...

The second helpful fact is that each grid sample can be thought of as representing a sphere centred at the grid point with radius equal to its distance value. It will be tangent to the unknown surface at (a minimum of) one point, ie. the closest point. In 2D we can visualize the problem in this way, and it becomes apparent that information regarding the presence of sharp features is not always entirely destroyed in the conversion to a signed distance field.

Examples

Input 2D polygon mesh: Marching triangles reconstruction (at the grid resolution):Visualized signed distance field (blue = inside, red = outside):
The far right image illustrates that visually at least, sharp features can be identified quite readily, assuming there is sufficient grid resolution around them. In particular, in areas of the grid where there is no sign change (eg. top left) it is still quite obvious where the sharp feature should be.

Observations:
Since the signed distance disks are always tangent to the closest surface point, sharp features often occur where many of the disks intersect at the same point. In 3D, a further challenge will be distinguishing sharp points from sharp edges.
Any point at which that two spheres of opposite sign are tangent to one another represents a surface point that is known exactly.
Somewhat non-local information is likely important to achieving a good reconstruction. If we only consider points within one cell length of the surface, sharp features are harder to visually identify with certainty:

Code

Here's the starter code with which I generated the above 2D images. dist_field_visualize.zip

Caveats

Many operations that we might typically apply to a signed distance field will destroy the signed distance property to some extent, making a surface reconstruction method based on this property less effective. Hence this problem is primarily of theoretical value for the moment, though I suspect applications could be found with a bit of effort (eg. improving liquid surface tracking).
Furthermore, my examples assumed that our distance data came from a well-resolved polygon in the first place. Presumably similar behaviour would occur for curved surfaces with sharp features..

Related research

An interesting loosely related problem is trying to assign signs to unsigned distance data: Signing the Unsigned: Robust Surface Reconstruction from Raw Pointsets.

So that's where it stands for the moment. I'd be interested to know if anyone has an obvious solution, if there's any important references I missed or points I should add, or if you have thoughts about this problem in general. Feel free to drop me a line about it at christopherbatty@yahoo.com.

-Christopher Batty, March 5, 2011