Recent Advances in the Level Set Method - Biomedical Image Analysis - page 215

the velocity field in the first place. This approach was introduced by Adalsteins-
son and Sethian [3] as an alternative to the use of reinitialization. As noted in
the introduction, velocity extensions also serve the purpose of constructing a
velocity field for the entire domain of φ even when the speed, F , is defined only
on the interface itself.
For velocity extensions, the objective is to force the velocity field, F ,tobe
such that the signed distance function is preserved, i.e.

∇ φ ·∇ φ ≡ 1 .

(4.31)

Differentiating Eq. 4.31 with respect to t , and using Eqs. 4.5 and 4.31, gives

∇ φ ·∇ F = 0 .

(4.32)

Viewed geometrically, Eq. 4.32 makes sense because it requires the speed func-
tion normal to the interface to be constant along that normal. This effectively
keeps the level sets of φ evenly spaced.
To solve Eq. 4.32, assume the function F is given on the zero level set of φ .
The goal is to construct an extension velocity F ext , such that
F ext φ = 0 = F φ = 0

∇ F ext

and ·∇ φ = 0 . (4.33)

The solution of Eq. 4.33 is done in a manner very similar to the fast marching
method. The discretization of Eq. 4.33 is given by

min( D + x φ i , j , 0) D + x F ext

i , j + max( D − x φ i , j , 0) D − x F ext i , j

+ min( D + y φ i , j , 0) D + y F ext

i , j + max( D − y φ i , j , 0) D − y F ext

i , j = 0 .

(4.34)

This is a linear equation in terms of the unknown F ext

i , j and is easily solved. Note
that Eq. 4.34 must be solved at the grid points x i , j in the order of increasing mag-
nitude of φ i , j similar to the fast marching method. This is easily accomplished
using the same heap-sort strategy described in Section 4.2.3.
The initialization of F ext

i , j on the grid points near the interface φ − 1 (0) is done
using the bicubic interpolation method discussed in Section 4.2.3. Given a grid
point x i , j , the point, y , on the interface φ − 1 (0) nearest to x i , j is computed using
the bicubic interpolant. The value of F ext

i , j must be the same as F ( y ), because the
vector x i , j − y is orthogonal to the interface, and hence parallel to ∇ φ ,so F ext
must be constant along that vector. This populates the grid points adjacent to
the interface, and the velocity extension algorithm can then proceed.
The algorithm for velocity extensions is therefore given by