the Level–Set representation is, in fact, a volume representation.
Figure 7.2: Figure illustrates how the exterior part of the source volume (oval) shrinks while the interior expands to approach the target volume.
The work of Breen and Whitaker rests on the observation that if the shape in the source volume is deformed by a speed function that is the signed distance function of the target volume, the shape in the source volume will deform to the target volume. When using this method, the behaviour of the deforming surface is more clearly defined: The part of the source solid that is outside the target solid will shrink, and the part that is exterior expands. Genus change is not a problem since the expanding and collapsing surfaces may merge and break. However, the initial surfaces must overlap.
The Level–Set Method can be used for many things beside morphing. For in-stance, if curvature is used to specify the speed function, smoothing is possible as well as filleting and blending [180]. The add blob deformation can be imple-mented through a speed function with the shape of a bump, say a 3D Gaußian.
Outside the volume graphics community, 2D Level–Set methods have been used to implement things like morphological manipulations [144], and shearing [163]
that can easily be extended to 3D.
7.3 The Level–Set Method
The Level–Set method [156] is a technique for tracking the evolution of a deform-ing interface or surface. Assume that we are dealdeform-ing with a surface B(t)⊂R3 wheretis the time parameterization. Bis assumed to change according to some speed function that pushesB in the normal direction. The speed function may depend on the geometry ofBor be completely independent ofB. A good exam-ple of the latter is a speed function that is always constant causingBto grow by constantly moving in the normal direction. A good example of a speed function that does depend on B is one that depends on the curvature ofB and pushes
the surface in the direction of the curvature centre. Such a speed function will smooth the surface and can be very useful.
The Level–Set method tracks the motion ofBin the normal direction, and this is expressed by a relationship with an embedding function Φ :R3×R+ →R.
For all points onB the value of Φ must be zero. This leads to the equation
Φ(B(t), t) = 0 (7.1)
where B(t) denotes a given point onB at time t. (7.1) simply says that B(t) is an isosurface (here called a level–set) of Φ(·, t). Because this holds for any point in time, both B and Φ may evolve but the Level–Set equation continues to hold implying that also
dΦ(B(t), t)/dt= 0 (7.2)
To see how the change of Φ andB are coupled, we take the derivative of (7.1) using the chain rule
i. Because all motion is in the normal direction, we can write the change ofB in terms of a speed functionF times the normal ||∇∇ΦΦ||
dB(t)
dt =F ∇Φ
||∇Φ|| (7.5)
Plugging this equation back into (7.4), we obtain the Level–Set equation
∂Φ
∂t +F||∇Φ||= 0 (7.6)
The Level–Set Method works on a discrete grid representation of Φ, that is3 Φn[i, j, k] = Φ(i∆x, j∆y, k∆z, n∆t)
This is a 4D discrete function, but, in general, only one time step is stored. In other words, Φ is really represented as a 3D voxel grid. Moreover, the initial value, Φ0is typically a distance field. In other words, the voxel gridsGthat are used throughout this thesis are precisely the same type of representation as the discretized embedding function Φ that the Level–Set Method works on.
3For simplicity we will assume in the following that unit time step is used and that the grid spacing is unit.
7.3 The Level–Set Method 139
The Level–Set Method is, essentially, a solver for an initial value problem: Given a Φ0what is the value at time stepn. With a little care, it is possible to ensure that Φn remains a distance field. That is important since the distance field property is one of the things we wish to be preserved by manipulations.
7.3.1 Discrete Implementation
The basic Level–Set Method as introduced by Osher and Sethian is to approxi-mate the time derivative with the forward difference operator
∂Φ
∂t ≈D+tΦ = Φn+1[i, j, k]−Φn[i, j, k] (7.7) This leads to an algorithm for computing the solution to (7.6) from an initial value Φ0
Φn+1[i, j, k] = Φn[i, j, k]−F||∇Φn|| (7.8) where the gradient||∇Φn||must be computed in the upwind direction. IfF ≥0
||∇Φn||2=
At first this upwinding seems to be a bit odd; why not simply approximate the gradient with central differences? The answer is that F indicates which way information propagates, and the gradient should be approximated using only voxels that lie in the direction whence information comes. If this principle is not obeyed, the numerical solution can easily become unstable in the presence of discontinuities. A more mathematical explanation is that the upwinding scheme is necessary because Φ might have discontinuities in which case the differential equation doesn’t have a normal solution. However, an integral form of the equation can have a weak solution and the upwinding is a part of this weak solution. This is explained in [156] but the discussion of weak solutions
is sketchy. To fully appreciate the issues, insight into the field ofconservation laws [96] is required.
What time–step is appropriate? A condition known as the CFL (Courant Friederichs Lewy) condition asserts that given a first order scheme like the one discussed above, the speed function must obey
maxF ≤∆x/∆t (7.11)
If we consider only grids with unit spacing and unit time step, this reduces to the simple condition that the speed function should not exceed 1. The CFL condition is mentioned by Sethian [156] and explained more deeply by LeVeque [96].
If F depends on the curvature of the evolving surface, discontinuities do not occur because the rˆole of curvature is to keep things smooth. Sethian suggests using central differences both for the gradient and for the computation of the first and second order partial derivatives involved in computing the curvature.
In 3D there is more than one type of curvature, but in the context of the Level–
Set Method mean curvature [80, 32] is almost always used.
7.3.2 Narrow Banding
A simple implementation of the Level–Set Method would update all voxels ac-cording to (7.8) but that is highly inefficient. In this context, we are only in-terested in the zero–set of Φ which corresponds to the boundary of an evolving solid. To address this problem,narrow band methods are used [156].
When using the narrow band approach, the values of Φ are updated only in a narrow band around the zero–set of Φ. Since the zero–set moves, we have to reinitialize the narrow band at times. This is done by tagging certain voxels as being “land mines”. When the evolving zero–set crosses a land mine, the narrow band is reinitialized. This is typically done by using the fast marching method to recompute distances to the zero–set.
Whitaker introduced a variation of this approach called the sparse–field tech-nique [178]. The idea is to keep track of a so–calledactive set which is the set of all voxels that are immediately adjacent to the zero–set – in the sense that one of their 6–neighbours are on the other side of the zero–set. From the active set the distances are propagated to neighbouring voxels using the city block distance. In practical terms, for each voxel that is 6–neighbour to a voxel in the active set, we set the distance to 1 plus the value of the active set neighbour.
7.3 The Level–Set Method 141
Whitaker keeps track of the voxels that are within 1/2vudistance to the surface, and these voxels are updated using the Level–Set Method. The rest are updated layer–wise. If a voxel,a, that is not in the active set has a neighbour,b, in the active set, the distance value ofabecomesa=b+1. Further layers are computed similarly. This method is faster than the Fast Marching Method but less precise.
7.3.3 Extension Velocities
It is important to note one problem with (7.6). The speedF is defined as the speed of the evolution of the surface B, and it is not always trivial to decide what valueF should take away from the surface.
Adalsteinsson and Sethian demonstrate in [1] (and chapter 11 of [156]) that if the initial Φ is a distance field and the gradient of the speed function is everywhere orthogonal to∇Φ i.e.∇F·∇Φ = 0 thendtd||∇Φ||= 0. In other words Φ remains a distance field. Certain other conditions must hold: For instance, the gradient is not well defined if the surface has a sharp edge. Never the less, this gives us a good indication that the speed function should be constant along the direction of∇φ. Because, the gradient of the speed function∇F is only orthogonal to∇Φ ifF is constant in the direction of ∇Φ. Consequently, a good way of extending the speed function from the surface to R3 is by finding the closest point onB and evaluating the result.
How to do this in practice is a new problem which Adalsteinsson solves by propagating the speed function using a variation of the fast marching method that propagates not only distances but also the speed function. However, a simpler approach and the one taken here is to simply find the foot point4 of a given voxel and evaluate the speed function at the foot point.
7.3.4 The CIR Scheme
The CIR (Courant Isaacson Rees) scheme has recently been used to solve the Level–Set equation by John Strain [162]. Say we are following the characteristic curves(t) defined by
In other words, Φ is constant alongs. At any given point, we can approximate a step alongs by the speed function times the gradient, and that leads to the CIR scheme which is, essentially, to track the characteristic curve from a voxel position one time–step back and then assign the value at that point.
The algorithm as implemented by Strain consists of three steps carried out for all grid points. Let the grid point be x. First we evaluate the speed function F(x). A step back along the characteristic is approximated by
s=x−F(x)∇Φ (7.14)
where (as usual) unit time step is assumed. The value of Φ at sis computed.
Strain uses the so–called ENO scheme [96] to find the value ats(which is not in general a grid point), but trilinear interpolation is also a good choice. Finally, the interpolated value Φ(s) is assigned to the grid pointx.
When all grid points have been updated, the entire grid is redistanced by re-placing the values at all grid points with the computed distance to the zero–set.
Strain claims that several features of his method makes it possible to use larger time steps than with other methods. The reasoning rests on showing heuristi-cally that the CFL condition is satisfied. Also, experiments indicate that the method converges to the exact solution when the time step is refined.