# Technical Explanation of Level Set Methods

Given a moving closed hypersurface G(t), we wish to produce an Eulerian formulation for the motion of the hypersurface propagating along its normal direction with speed F, where F can be a function of various arguments, including the curvature, normal direction, etc. The main idea is to embed this propagating interface as the zero level set of a higher dimensional function phi. Let phi(x,t=0), where x is in n-dimensional space, be defined by

phi(x,t=0) = d(signed)

where d(signed) is the distance from x to G(t=0), and the plus (minus) sign is chosen if the point x is outside (inside) the initial hypersurface G(t=0). Thus, we have an initial function phi(x,t=0) with the property that

G(t=0) = ( x | phi( x, t= 0) = 0 )

Our goal is to now produce an equation for the evolving function phi(x,t) which contains the embedded motion of G(t) as the level set phi = 0. Let x(t) be the path of a point on the propagating front. That is, x (t=0) is a point on the initial front G(t=0), and dx/dt = F (x(t)) with the vector dx/dx normal to the front at x(t). Since the evolving function phi is always zero on the propagating hypersurface, we must have

phi( x ( t) , t ) = 0

By the chain rule,
d(phi)/dt + grad ( x(t,t) ) * dx/dt= 0

where grad is the gradient operator, and the * denotes the dot product. Since F already gives the speed in the outward normal direction, then dx/dt * n = F, where n = grad phi /|grad phi|. Thus, we then have the evolution equation for phi(x,t), namely

d(phi)/dt + F | grad phi | = 0

with
phi (x,t=0 ) given

We refer to this as a Hamilton-Jacobi ``type'' equation because, for certain forms of the speed function F, we obtain the standard Hamilton-Jacobi equation. There are four major advantages to this Eulerian Hamilton-Jacobi formulation.

• First,the evolving function phi(x,t) always remains a function as long as F is smooth. However, the level surface phi = 0, and hence the propagating hypersurface G(t), may change topology, break, merge, and form sharp corners as the function phi evolves.
• The second advantage of this Eulerian formulation concerns numerical approximation. Because phi(x,t) remains a function as it evolves, we may use a discrete grid in the domain of x and substitute finite difference approximations for the spatial and temporal derivatives. For example, using a uniform mesh of spacing h on an orthogonal grid, with grid nodes (i,j), and employing the standard notation that phi(i,j,n) is the approximation to the solution phi(i h , j h , n k), where k is the time step, we might write

(phi(i,j,n+1)- phi(i,j,n))/k + ( F ) ( grad(i,j) phi(i,j,n)) = 0

where we have used forward differences in time, and grad(ij) phi(i,j,n) represents some appropriate finite difference operator for the spatial derivative. As discussed above, the correct entropy-satisfying approximation to the difference operator comes from exploiting the technology of hyperbolic conservation laws. Following (Sethian, J.A.. Level Set Methods, Cambridge University Press, 1996), given a speed function F(K), we update the front by the following scheme. First, separate F(K) into a constant advection term F0 and the remainder F1 (K), that is,

F(K) = F0 + F1 (K)

The advection component F0 of the speed function is then approximated using upwind schemes, while the remainder is approximated using central differences. In one space dimension with positive F0, we have

phi(i,n+1) = phi(i,n) - (k) (F0) sqrt [max (DM phi(i),0)^2 + min (DP phi(i),0)^2 ]
- | F1(K) grad phi(i,n) |
where DM and DP are the backwards and forwards difference operators. Extension to higher dimensions are straightforward.
• The third major advantage of the above formulation is that intrinsic geometric properties of the front may be easily determined from the level function phi. Both the normal derivative and the curvature may be easily calculated by finite difference approximations.
• Finally, the fourth major advantage of the above level set approach is that there are no significant differences in following fronts in three space dimensions. By simply extending the array structures and gradients operators, propagating surfaces are easily handled.