Eikonal Solver for the Distance Function to Triangular Meshes

A complete distance field contours (white) updated in 8 sweeps (red) of the Fast Sweeping Algorithm, with grid resolution 38 x 38 x 34.

I would start with a warning, that this post is going to be slightly more mathematical than the previous ones, so for the “copy-pasters” out there: feel free to skip to the code snippets which you can adjust according to your needs. If you’re not in a hurry and curious about how this method works, do read through my take on the underlying theory.

After implementing a fast enough approach to find the voxel outline of a polygonal mesh, we have everything we need to proceed to the next step of completing the distance field computation. To elaborate on the current state of our scalar grid with domain G = [G_{x, \min}, G_{x, \max}] \times [G_{y, \min}, G_{y, \max}] \times [G_{z, \min}, G_{z, \max}] with a subset of grid points \textbf{x}_{i,j,k}^\Gamma corresponding to voxel cells intersecting the mesh (boundary) \Gamma = \overline{X} .

Mathematically, the distance function: d^{+} : \mathbb{R}^3 \rightarrow \mathbb{R}^{+} : \textbf{x} \mapsto \inf \big\{ \| \textbf{p} - \textbf{x} \| \big| \ \textbf{p} \in \Gamma , \ \Gamma \subset \mathbb{R}^3 \big\} is a solution to a hyperbolic partial differential equation with a Dirichlet boundary condition

\| \nabla d^{+}(\textbf{x}) \| = 1 \ , \ \ d^{+}(\textbf{x})\big|_{\Gamma} \equiv 0

which is a specific case for a general n-dimensional Eikonal equation: \| \nabla u(\textbf{x}) \| = f(\textbf{x}) with general boundary condition u(\textbf{x})\big|_{\Gamma} = \phi(\textbf{x}) , \Gamma \subset \mathbb{R}^n . The equation basically says that the distance function d^{+} to a set \Gamma \subset \mathbb{R}^3 is linear with maximum slope 1 everywhere (except at the boundary \Gamma and the object’s skeleton where its first partial derivatives are discontinuous in at least one direction), which is exactly what distance should behave like.

The boundary condition says that the distance to an object at any of the object’s points is zero, and it increases linearly with slope 1 in the direction “away from” the source. The gradient operator \nabla computes a vector at a given point which points in the direction of the steepest growth of a function and its magnitude corresponds to the rate of change in this direction.

Eikonal Solvers

Since the Eikonal equation is a non-linear hyperbolic PDE of second order, its exact solutions are hard to find. In fact, there are only a few exact solutions for the distance function equation (e.g.: the 2 -sphere with radius r centered at the origin: d^{+}(x,y,z) = x^2 + y^2 + z^2 - r^2 ). This means that the solution d^{+} needs to be obtained numerically for most shapes the source domain (boundary) \Gamma can acquire.

Hyperbolic equations describe “conservation laws“, that is, some quantity is conserved in the system’s evolution. For example the solutions to the wave equation \partial_{tt} u(\textbf{x}, t) - c^2 \Delta u(\textbf{x}, t) = \phi(\textbf{x}) conserve mechanical energy as waves with given parameters (as initial and boundary conditions) propagate through a medium. Hyperbolic systems of first-order PDE’s can describe various flow problems. A key property is that their solution carry information along a set of curves called characteristics. First-order quasilinear and hyperbolic systems can be reduced to a system of ordinary differential equations and solved for a solution along characteristics which follows the conservation law.

In particular, the characteristics of the Eikonal distance equation \| \nabla d^{+}(\textbf{x}) \| = 1 are straight lines along which their derivative is always 1 : \partial_\textbf{x} d^{+} (\textbf{x}(s)) \cdot \textbf{x}'(s)= 1 . By these means the characteristics of the distance function are perpendicular to the outline \Gamma .

Characteristics cannot intersect (and for linear hyperbolic systems they follow this property) because otherwise the system of ODE’s would not have a unique solution at these intersection points. The reader then notices that for the distance function d^{+} these lines surely intersect at some points, particularly when they’re emitted inward into the object’s interior where they form the object’s skeleton \sigma , or in the outward direction, where they form shockwaves \rho as they meet from non-convex parts of \Gamma . The formation of shockwaves or rarefaction waves is typical of non-linear hyperbolic systems, so even though d^{+} is bound to not have continuous first derivatives in these regions the conservation law holds. The flow of information in shock and rarefaction wave regions is determined by the flow along incoming or outflowing characteristics (see the image below).

A contour plot of the distance function to a star-shaped outline \Gamma emanating characteristics (thin black) which form shockwave regions \rho and peak values of the distance field at the object’s skeleton \sigma

The property of outward propagation of information is used in algorithms computing numerical solutions to the Eikonal problem. In this section, we mention two possibilities, the latter of which will be used.

The front of “neighbours” in the Fast Marching algorithm

One possibility is to use the Fast Marching Method (FMM) which utilizes the same principle as Dijkstra’s algorithm, namely that the information flows outward from a given seeding area by updating the neighbor with the smallest value. The worst case complexity of the FMM is \mathcal{O}(M \log{M}) where M is the grid size. The logarithm orginates from a heap sort step applied for each iteration for acceleration . The total heap size corresponds to the possible size of the heap data structure in which we store the neighboring values of the updated front.

In a 2005 paper by Hongkai Zhao the 90’s FMM scheme by Sethian was dwarfed by the Fast Sweeping Algorithm with linear complexity \mathcal{O}(M) with respect to the grid size. The speed of the method lies in the fact that the Fast Sweeping is applied to Eikonal problems with characteristics that do not change very much, which is the case for the distance function whose characteristics are rays emitted from a source (which might collide in shock/rarefaction fronts).

Zhao’s Fast Sweeping Method

The input data for the Fast Sweeping Method is a domain G with a pixel/voxel outline with exact values as boundary condition (the source of the distance field) and a “large enough” number(such that no possible distance value would be larger) everywhere else (e.g.: DBL_MAX ). In n dimensions the method needs to perform 2^n “sweeps“, that is: iterations with alternating index directions across =domain G , which is a min-max (axis aligned bounding box) in our case. In particular, for n=2 we sweep 4 times from each corner of a rectangular domain:

i_x = 1, ..., N_x \ , \ \ i_y = 1, ..., N_y
i_x = N_x, ..., 1 \ , \ \ i_y = 1, ..., N_y
i_x = N_x, ..., 1 \ , \ \ i_y = N_y, ..., 1
i_x = 1, ..., N_x \ , \ \ i_y = N_y, ..., 1

A C++ implementation of the method, given a grayscale image giving rise to the input distance grid distGrid with marked outline pixels as “frozen”, and DBL_MAX everywhere else would look something like this:

void FastSweep2D(double* distGrid, bool* frozenCells, ImageParams* imgPar)
{
	const int height = imgPar->height, width = imgPar->width;
	const int row = imgPar->row;

	const int NSweeps = 4;
	// sweep directions { start, end, step }
	const int dirX[NSweeps][3] = { {0, width - 1, 1} , {width - 1, 0, -1}, {width - 1, 0, -1}, {0, width - 1, 1} };
	const int dirY[NSweeps][3] = { {0, height - 1, 1}, {0, height - 1, 1}, {height - 1, 0, -1}, {height - 1, 0, -1} };
	double aa[2], eps = 1e-6;
	double d_new, a, b;
	int s, ix, iy, gridPos;
	const double h = 1.0, f = 1.0;

	for (s = 0; s < NSweeps; s++) {

		for (iy = dirY[s][0]; dirY[s][2] * iy <= dirY[s][1]; iy += dirY[s][2]) {
			for (ix = dirX[s][0]; dirX[s][2] * ix <= dirX[s][1]; ix += dirX[s][2]) {

				gridPos = iy * row + ix;

				if (!frozenCells[gridPos]) {

					// === neighboring cells (Upwind Godunov) ===
					if (iy == 0 || iy == (height - 1)) {
						if (iy == 0) {
							aa[1] = distGrid[gridPos] < distGrid[(iy + 1) * row + ix] ? distGrid[gridPos] : distGrid[(iy + 1) * row + ix];
						}
						if (iy == (height - 1)) {
							aa[1] = distGrid[(iy - 1) * row + ix] < distGrid[gridPos] ? distGrid[(iy - 1) * row + ix] : distGrid[gridPos];
						}
					}
					else {
						aa[1] = distGrid[(iy - 1) * row + ix] < distGrid[(iy + 1) * row + ix] ? distGrid[(iy - 1) * row + ix] : distGrid[(iy + 1) * row + ix];
					}

					if (ix == 0 || ix == (width - 1)) {
						if (ix == 0) {
							aa[0] = distGrid[gridPos] < distGrid[iy * row + (ix + 1)] ? distGrid[gridPos] : distGrid[iy * row + (ix + 1)];
						}
						if (ix == (width - 1)) {
							aa[0] = distGrid[iy * row + (ix - 1)] < distGrid[gridPos] ? distGrid[iy * row + (ix - 1)] : distGrid[gridPos];
						}
					}
					else {
						aa[0] = distGrid[iy * row + (ix - 1)] < distGrid[iy * row + (ix + 1)] ? distGrid[iy * row + (ix - 1)] : distGrid[iy * row + (ix + 1)];
					}

					a = aa[0]; b = aa[1];
					d_new = (fabs(a - b) < f * h ? (a + b + sqrt(2.0 * f * f * h * h - (a - b) * (a - b))) * 0.5 : std::fminf(a, b) + f * h);

					distGrid[gridPos] = distGrid[gridPos] < d_new ? distGrid[gridPos] : d_new;
				}
			}
		}
	}
}

Drawing the result onto an image with pixel values from [0, 255] we simply normalize the output distance field between min and max values [d_{\min}, d_{\max}] and get:

Note that the zero-value outline was inverted to white color to be distinguishable from dark gray colors close to the outline \Gamma .

The specific details of each step in each sweep rely on mathematics provided by Zhao’s paper. Of course, only cells that are not frozen are processed. And only cells in \Gamma possess such property. Moreover, f = 1 is the value of the right-hand side of the Eikonal equation (i.e.: the velocity function) and h = 1 is the (index) step size of the numerical method corresponding to the distance s as used in my previous post for cell size.

The Godunov upwind difference scheme corresponds to a finite volume scheme (for pixels/voxels) for conservation laws (given a hyperbolic PDE) such that it is sensitive to the direction of propagation in a flow problem (along characteristics). In Zhao’s 2D formulation for a grid point \textbf{x}_{i,j} with value d_{i,j}^h :

[(d_{i,j}^h - d_{x, \min}^h)^{+}]^2 + [(d_{i,j}^h - d_{y, \min}^h)^{+}]^2 = f_{i,j}^2 h^2
i = 2, ..., N_x - 1 \ , \ \ j = 2 ,..., N_y - 1

Where d_{x,\min}^h = \min \{ d_{i-1, j}^h, d_{i+1,j}^h \} \ , \ \ d_{y,\min}^h = \min \{ d_{i, j - 1}^h, d_{i,j + 1}^h \} and operator ( \cdot )^{+} corresponds to:

The individual (Gauss-Seidel) iterations look for the unique solution to a quadratic equation

[(x-a)^{+}]^2 + [(x-b)^{+}]^2 = f_{i,j}^2 h^2

where a = d_{x,\min}^h and b = d_{y,\min}^h , that is:

Values of a and b at the boundary of G use a one-sided difference, for example:

For reference, see line 39 of the code snippet. We use a full array of coefficients a = a_1 \ , \ \ b = a_2 to be able to extend the algorithm to n dimensions with a quadratic equation:

[(x-a_1)^{+}]^2 + [(x-a_2)^{+}]^2 + ... + [(x-a_n)^{+}]^2= f_{i,j}^2 h^2

The condition that holds for n = 2 has to be generalized for higher dimensions. In particular, we need to order the neighbor values: a_1 \leq a_2 \leq ... \leq a_n (by sorting array aa in our implementation). Then we recursively add values a_k \ , \ k = 2,...,n until we find a solution \overline{x} of quadratic equation

(x-a_1)^2 + (x-a_2)^2 + ... + (x-a_k)^2= f_{i,j}^2 h^2

such that a_k < \overline{x} \leq a_{k+1} ( a_n < \overline{x} < \infty ). The Gauss-Seidel iterations need to be carried out from each of the 2^n corners of the grid. See this applied to a 3D version of the method:

void FastSweep3D(double* grid, bool* frozenCells, GridParams* gridPar) {
   /* ... */
   /* initialization & iterating from 8 corners */
   
   uint gridPos = ((iz * Ny + iy) * Nx + ix);
   
   if (!frozenCells[gridPos]) {
    // === neighboring cells (Upwind Godunov) ...
    // ... same as for 2D image, but with all domain walls & corners ===
    // ...
        // simple bubble sort
        if (aa[0] > aa[1]) { tmp = aa[0]; aa[0] = aa[1]; aa[1] = tmp; }
        if (aa[1] > aa[2]) { tmp = aa[1]; aa[1] = aa[2]; aa[2] = tmp; }
        if (aa[0] > aa[1]) { tmp = aa[0]; aa[0] = aa[1]; aa[1] = tmp; }

        double d_curr = aa[0] + h * f; // just a linear equation with the first neighbor value
        double d_new;
        if (d_curr <= (aa[1] + eps)) {
            d_new = d_curr; // accept the solution
        }
        else {
            // quadratic equation with coefficients involving 2 neighbor values aa
            double a = 2.0; 
            double b = -2.0 * (aa[0] + aa[1]);
            double c = aa[0] * aa[0] + aa[1] * aa[1] - h * h * f * f;
            double D = sqrt(b * b - 4.0 * a * c);
            // choose the minimal root
            d_curr = ((-b + D) > (-b - D) ? (-b + D) : (-b - D)) / (2.0 * a);

            if (d_curr <= (aa[2] + eps))
                d_new = d_curr; // accept the solution
            else {
                // quadratic equation with coefficients involving all 3 neighbor values aa
                a = 3.0;
                b = -2.0 * (aa[0] + aa[1] + aa[2]);
                c = aa[0] * aa[0] + aa[1] * aa[1] + aa[2] * aa[2] - h * h * f * f;
                D = sqrt(b * b - 4.0 * a * c);
                // choose the minimal root
                d_new = ((-b + D) > (-b - D) ? (-b + D) : (-b - D)) / (2.0 * a);
            }
        }
        // update if d_new is smaller
        grid[gridPos] = grid[gridPos] < d_new ? grid[gridPos] : d_new
   }
       
   /* for-loops scope ends */
   /* ... */
}

The grid corners with step s as triples (e.g: (-1,1,1) for corner 3 ) can be seen in the first image below:

The procedure of updating field voxels according to radiating characteristics of the Eikonal problem relies on a simple geometric fact that the solution gets updated by the velocity * distance step, i.e.: f_{i,j}h weighed by the degree of coincidence of the field gradient (characteristic velocity) with the sweeping direction at each voxel.

After a closer look at the 2D algorithm’s logic applied to each voxel, we notice that if a_1 and a_2 are coordinates of the center of a circle, and r = f_{i,j}h its radius, then whenever the identity line y = x intersects this circle, we have |a_1 - a_2| < f_{i,j}h = r and we look for the farthest solution of the quadratic for intersection (I.). If the circle gets too far from the identity line (II.), we simply update the minimal coordinate by the radius r = f_{i,j}h .

By sorting neighbor cell values a_1, a_2, ..., a_n in higher dimensions, we secure the positivity of the difference a_{k} - a_{k-1} and hence we first consider an intersection of a k -sphere with the identity line before moving on to higher (k+1) -sphere. We can think of this as testing whether the possible gradient vector at grid point \textbf{x}_{i,j,k} differs too much from the characteristic velocity and in that case we move on to a higher-dimensional intersection.

The updating procedure can move against the direction of characteristics when it updates voxels with LARGE_VAL value (those that have not yet been processed). A closer look at the animation for the star image shows that during the first sweep the grid values in the interior \text{Int}{\Gamma} of contour \Gamma are updated as if there was no “other side” of the contour. The following sweep then corrects the discrepancies in the negative x direction.

The “closeness” of the k -sphere (k = 1, 2, ..., n ) centered at (a_1, a_2,..., a_k)^\top to the identity line \textbf{y} = \textbf{x} , corresponds to the balance of the neighboring cell values in the direction of sweeping. In other words, if we sweep through i = 1, ..., N_x \ , \ \ j = 1 ,..., N_y (from the top left corner), and the values d^{h}_{i-1, j} (previous row) and d^{h}_{i, j - 1} (previous column) are chosen as minimal (e.g.: less than LARGE_VAL) and balance out: d^{h}_{i-1, j} \approx d^{h}_{i, j - 1} , then current cell value d^{h}_{i, j} is computed as the maximal intersection of the circle with the identity line. If, for example, d^{h}_{i-1, j} > d^{h}_{i, j - 1} then the flow is biased in the y-direction and the minimal value d^{h}_{i, j - 1} from the previous column incremented by f_{i,j}h = r is considered.

When the same cell is re-visited in the next sweep it gets updated when it is necessary. To make sure that distance values are not biased in the direction of the last sweep, we verify whether the computed new value candidates are lower than the current cell value. If that is the case, we should, of course, consider the minimum distance value, otherwise we need to keep the old distance value. The latter should happen in the interior of \Gamma where a gradient discontinuity forms at the object’s skeleton. It would make no sense to overwrite correctly computed distance values when sweeping from the opposite direction.

  • 8 sweeps animation
  • Complete distance field
  • absError = | d_FS – d_Exact |

Final Remarks

The complete procedure for computing a distance field (without computing a sign dependent on position in \text{Int}{\Gamma} or \text{Ext}{\Gamma} – interior/exterior of the contour) can be summarized by the following pseudocode:

The computation times (in seconds) on my machine for different grid resolutions and individual steps of the overall algorithm are listed in the table below:

Clearly, the demands get higher for larger grids, but the overall time is orders of magnitude better than the brute force approach. I still computed the exact distance field d^{\text{exact}} to the bunny mesh for the 60^3 grid resolution (taking over 20 mins) in order to compute the global error of the numerical method. The type of error I used for evaluating the method’s precision is referred to as the L_2 -error (based on the L_2 Lebesque integral norm on function space L_2(G, \mathbb{R}) ). In discrete terms:

\epsilon_{N_x, N_y, N_z} = \sqrt{\frac{1}{N_x N_y N_z} \sum_{i,j,k = 1}^{N_x, N_y, N_z} (d_{i,j,k}^{\text{FS}} - d_{i,j,k}^{\text{exact}})^2 }

where d^{\text{FS}} is the distance field computed by the Fast Sweeping method. The pointwise absolute error |d_{i,j,k}^{\text{FS}} - d_{i,j,k}^{\text{exact}} | usually accumulates towards the grid corners, as can be seen in the image for the bunny mesh.

To verify the computation one should compare L_2 -errors from different grid resolutions, e.g.: N_A^3 and N_B^3 such that N_A < N_B using the experimental order of convergence:

\text{EOC} = \log_{N_B / N_A} \big( \frac{\epsilon_{N_A}}{\epsilon_{N_B}} \big) \approx 2

In the next post I will finalize the signed distance computation by introducing an easy-to-implement method using voxel flood fill for sign computation.

Leave a comment