13.6 A Fast and Precise Level Set Algorithm

The idea leading to fast level set algorithms stems from observing that only the values of the level set function near its zero level set are essential, and thus only the values at the grid points in a narrow band around the zero level set have to be calculated. As the zero level set moves, the signed distance function in the narrow band must be maintained (cf. Section 13.3.2).

Both extending the speed function and narrow banding require constructing the distance function from the zero level set in the order of increasing distance. But calculating the exact distance function from a curve or surface consisting of a large number of small line segments or triangles is computationally expensive and can only be justified for initialization. An approximation to the distance function can be computed by a special fast marching method (cf. Section 13.3.1).

By intertwining both, extending the speed function and narrow banding, the following algorithm keeps expensive calculations to a minimum. Although the level set method is a seemingly computationally expensive method, since it requires solving a partial differential equation for describing surface evolutions, the computation time consumed for the surface evolution by narrow banding is negligible compared to that required for the physical models, i.e., radiosity or diffusion.

Combining narrow banding and extending the speed function into one algorithm provides several benefits. First the speed function is retained as the signed distance function throughout the simulation, which assures good accuracy until the end of the simulation. Second narrow banding reduces the number of active points that have to be updated from $O(n^2)$ to $O(n)$ in two dimensions and from $O(n^3)$ to $O(n^2)$ in three dimensions. By retaining the signed distance function the width of the narrow band can be kept down to two or three points on each side without decreasing accuracy. Third time consuming calculations are reduced to a minimum by intertwining the computations necessary for narrow banding and extending the speed function. Finally it is noted that the width of the narrow band can be adjusted if desired.

The algorithm works as follows (cf. Section 13.3.1). First the initial points near the zero level set, where the speed function is known, and the neighboring trial points are determined. In the main loop it is checked if there is still a trial point to be considered in the narrow band. All trial points are stored in a heap ordered by their distance to the zero level set. If there is a point to be considered, both its distance is approximated and its extension speed calculated, and its neighbors are updated accordingly. Finally after the main loop, bookkeeping information for the narrow band points is updated using distance information just computed.

In more detail the steps of this algorithm are as follows, where several auxiliary functions are used. Concerning the data structures, the information about the level set grid, the distance function, the extended speed function, and the tags for the fast marching method are stored in arrays. The trial points are stored in a heap [113].

1. First find the grid points whose speed function values are initially known. These values are computed in the physical simulation step and translated to the grid. Then compute the distance for the initial points, and tag them as known, and initialize the corresponding grid points of the speed function. Find the trial points which are the neighbors of the initially known points, and compute their tentative distance values. All other points are far points.
2. While there are trial points, do the following:
   1. Remove the first point from the heap and call it $a$. It has the smallest distance from the zero level set of all points in the heap.
   2. If narrow banding is used and the maximum width of the narrow band is smaller than the distance of $a$, return from the loop.
   3. Mark $a$ as known.
   4. For all neighbors $b$ of $a$, do the following:
      1. If $b$ is a far point, recompute its distance and speed function values and mark it as a trial point.
      2. If $b$ is a trial point, recompute its tentative distance and speed function values, unless it was computed in the previous step.
3. If narrow banding shall be used, set the distance values of all points which are not marked as known to the width of the narrow band. Set the sign of the distance function.
4. Finally return two values, namely the grid containing the new signed distance function and the grid containing the extended speed function.

Then a finite difference scheme from Section 13.5 is used for updating the level set grid. The main level set function must perform the bookkeeping for the narrow band considering the old narrow band from the previous iteration and the new narrow band. The points outside the old narrow band but inside the current narrow band are initialized to the signed distance function just computed.