High-order algorithms for computing closest points on implicitly-defined surfaces

A powerful technique for representing curves in two dimensions (N = 2) and surfaces in three dimension (N = 3) is to define them implicitly as a level set of a continuous, N-dimensional scalar function. This technique of embedding the surface geometry in a higher-dimensional function (which is often called the "level set function") leads to a wide array of mathematical and computational advantages, as exemplified by the level set method for moving interface problems, and the Voronoi implicit interface method for computing multiphase interface evolution.

Numerical methods making use of this idea often require accurate approximations of minimum distances to and closest points on implicitly defined surfaces. I have developed high-order accurate algorithms for this purpose and these are described in the paper:

Provided on this page is C++ code implementing the algorithms developed in the paper. The code is freely available (subject to the license) and can be used to implement, for example, high-order accurate reinitialisation/redistancing algorithms in level set methods. The code mainly applies to the case that the level set function is defined on a rectangular Cartesian grid. However, as discussed in the paper, it is possible to extend the algorithms to the case of unstructured grids – one could use the code as a starting point. For example, the implementations of the k-d tree and Newton's method could be used as-is without modification in such an adaptation.

C++ code

Most of the code is in C++ header files (.hpp). Except for two demonstration programs, there is only one .cpp file (poly.cpp) that needs to be compiled. This file contains the precomputed pseudoinverses of the interpolation matrices and the various stencils for performing least-squares polynomial interpolation on a rectangular Cartesian grid.

The package contains:

Instructions and code dependencies

Much of the code is templated on the dimension N. To assist with this functionality, the code makes use of blitz++, an open-source implementation of N-dimensional arrays and fixed-length vectors in C++ with convenient expression template techniques.

To compile and test the code:

  1. Download and install blitz++.
  2. Download (from below) the high-order closest point code.
  3. Compile:
    1. If you are using blitz-0.9, uncomment the line "#include <blitz/tinyvec-et.h>" in hocp_util.hpp. If you are using blitz-0.10, you do not need to do this.
    2. Modify the Makefile to correctly point to the location of the blitz directory in the include command.
    3. make, or compile according to your favourite methods.
  4. Run:
    1. The program "testSimple" demonstrates a very simple problem of reinitialising a level set function in 2D and 3D. Type in a dimension (2 or 3) and the number of grid points in each axis to have the program reinitialise a spherical interface level set function and report the maximum-norm error of the computed signed distance function.
    2. The program "testConvergence" implements a series of convergence tests on different interface geometries and can be used to replicate the results in the above mentioned paper. Use it by specifying the test problem via command-line parameters – for example, ./testConvergence 2 sphere or ./testConvergence 3 ellipsoid.


The source code is being made available according to a license similar to the BSD 3-Clause license. By downloading or using this code you are agreeing to the license in the file License.txt that is also contained in the packages. Please read it before downloading or using the code.

Download as a zip file (66 KB).

Comments or suggestions

Feel free to contact me if you have any comments or suggestions: rsaye {at} lbl {dot} gov