<< . .

( 23)

. . >>

|u0 (x) ’ C1 |2 dx + |u0 (x) ’ C2 |2 dx,
E1 (“) = (12.4)
image. The parameter µ de¬nes the scale for the method in the sense that
inside “ outside “
if µ = ∞, then the length of the boundaries should be minimized. So we
where “ is any curve and C1 , C2 are the averages of u0 inside “ and take u to be the average of u0 over the whole image; this is the coarsest
outside “. In this simple case it is obvious that “0 , the boundary of the scale. If µ = 0, then there is no penalty for length; each grid point (pixel) is
object, is the minimizer of the ¬tting term. See Figure 12.1. the average of u0 , or just the value of u0 ; and the segmentation u is equal
In the active contour model proposed in [35] and [34] the ¬tting term to the original image. As µ increases the segmentation coarsens. This is the
plus some regularizing terms will be minimized. The regularizing terms will idea behind the segmentation of [99]; it is a split and merge method, not a
involve the length of the boundary “ and the area of „¦, the region inside “. partial di¬erential equations-based approach, as in [36, 37].
This is in the spirit of the Mumford-Shah functional [117]. Thus, using the Returning to the model in equation (12.5), it is easy to see that with
variational level set formulation [175], the energy can be written, with φ respect to the constants C1 and C2 it is easy to express these two in terms
12.2. Active Contours Without Edges 123 124 12. Snakes, Active Contours, and Segmentation

of φ as forcing function on all the level contours of φ, forcing some to go through 0
u0 (x)H(φ(x)) dx
The scale parameter should be small if we have to detect many objects,
C1 (φ) = , (12.8)
H(φ(x)) dx including small objects. If we have to detect only large objects (for example
u0 (x)(1 ’ H(φ(x)) dx a cluster of lights), the scale parameter µ should be larger. An extremely
C2 (φ) = . (12.9)
trivial but slightly instructive analytic example is the following. Let µ = 0
(1 ’ H(φ(x)) dx
and »1 = »2 = », so the ¬nest-scale segmentation occurs; i.e., every point
This expresses the fact that the best constant value for the segment u is should be a boundary point and u(x) = u0 (x). That this does occur follows
just the average of u0 over the subregion. from the evolution
In order to compute the Euler-Lagrange equations we use the variational
‚φ C1 + C2
level set approach and arrive at = δ (φ) »(C2 ’ C1 ) u0 ’ , (12.14)
‚t 2
‚φ φ
’ ν ’ »1 (u0 ’ C1 )2 + »2 (u0 ’ C2 )2
= | φ| µ · so steady state can occur only if the average of u0 over „¦ equals its aver-
| φ|
age over „¦c (which is an unstable equilibrium) or if φ ≡ 0, which is the
desired equilibrium. If we take µ > 0 and take the limit as µ ’ 0, we
φ(x, 0) = φ0 (x). (12.11) believe intuitively that the u = u0 , φ ≡ 0 solution is the stable limit, since
the curvature term will tend to move the boundaries when these averages
However, it was found in [35, 34], that the nonmorphological approach
happen to be equal. Numerical experiments indicate that this is true, and
was more e¬ective; i.e., | φ| is replaced by δ (φ) in the term multiplying
hence an in¬nite (actually as many as there are grid points) number of new
the brackets in equation (12.10). Here
zero-level contours develop in a stable fashion.
‚ ‚1 2 z
δ (z) = H (z) = 1+ (12.12)
‚z ‚z 2 π
for > 0 and small, which gives a globally positive approximation to the 12.3 Results
delta function. This is necessary, as we shall discuss below. Thus the model
de¬ned in [35, 34] is We illustrate in Figures 12.2, 12.3, 12.4 and 12.5 the main advantages of
this active contour model without edges [34]: detection of cognitive contours
‚φ φ
’ ν ’ »1 (u0 ’ C1 )2 + »2 (u0 ’ C2 )2
= δ (φ) µ (12.13) (which are not de¬ned by gradients) in Figure 12.2, detection of contours
| φ|
in a noisy image in Figure 12.3, detection of interior contours automatically
with C1 and C2 de¬ned in equations (12.8) and (12.9). Generally, the pa- and extension to three dimensions in Figures 12.4 and 12.5. Also, note that
rameters are taken to be ν = 0, »1 = »2 = 1, and µ > 0 is the scale the initial curve does not need to enclose the objects, as in the classical
parameter. Although only two regions „¦ and „¦c can be constructed, they snakes and as in active contour models based on the gradient-edge detector.
can, and generally will, be disconnected into numerous components in the
¬ne-scale case, with each component having one of two constant values
for u.
12.4 Extensions
One important remark concerning this model as opposed to other level
set evolutions is its global nature. All level sets of φ have the potential to
be important. This means that other isocontours corresponding to nonzero As in the original Mumford-Shah functional [117] and the implementations
values of φ might evolve so as to push through the φ = 0 barrier and create of [99] and [34], one may propose the use of other channels, e.g., replacing u0
by the curvature of its level sets · ( u0 /| u0 |), or by their orientations
new segmented regions. Thus reinitialization to the distance function is
u0 = tan’1 ((u0 )x /(u0 )y ), to do texture segmentation.
not a good idea here (as pointed out by Fedkiw). One can even begin with
φ > 0 or φ < 0 throughout the region and watch new zeros develop. Of Chan et al. [32] extended the method to vector-valued images as fol-
course, this also explains why we need δ (z) > 0 in equation (12.12). Again, lows. Let u0,i be the ith channel of an image on the usual square region
the goal here is to detect interior contours. The technical reason why this with N channels and “ the evolving curve. See [32] for examples of these
works is that the image u0 acts in a nontrivial and nonlinear way as a channels, which include color images. The extension to the vector case is
12.4. Extensions 125 126 12. Snakes, Active Contours, and Segmentation

straightforward. The level set evolution becomes in arti¬cial time, as usual, is
‚φ φ
’ |u+ ’ u0 |2 + |u’ ’ u0 |2
‚φ φ 1 ·
= δ (φ) µ (12.22)
»+ (u0,i +
C i )2
· ’ ’
= δ (φ) µ (12.15) | φ|
| φ|
‚ N i=1
’ ν| u+ |2 + ν| u’ |2 .
»’ (u0,i ’ Ci )2 .

+ i
N The new numerical challenge is to obtain the numerical solution of the set
of elliptic boundary value problems in equations (12.18) to (12.21) for u’
and u+ in multiply connected regions. This is done by ¬rst extending u’ on
Here the Ci are the averages of u0,i on φ > 0 and φ < 0, respectively.
{x | φ(x) > 0} while retaining boundary conditions, and similarly for u+ .
Another extension, again using only one level set function, involves re-
There are several methods suggested; see [37] for a brief description. (One
moving the piecewise constant assumption and allowing piecewise-smooth
is the ghost ¬‚uid method of Fedkiw et al. [63].) See [37] for some interesting
solutions to the variational problem, smooth inside each zero isocontour
of φ, with jumps across the edges, as described in Chan and Vese [37]. This
The last extension is to get several, or indeed many, di¬erent regions
is in the spirit of the original Mumford-Shah functional, although multi-
corresponding to di¬erent level set functions. The idea is as follows. Based
ple junctions are not yet allowed; see below for that. The minimization
on the four color theorem, we can “color” all regions in a partition using
procedure is
only four colors such that any two adjacent regions have di¬erent colors.
F (u+ , u’ φ),
inf (12.16) Therefore, using two level set functions we can identify the four colors by
u+ ,u’ ,φ
the four possibilities φi > 0, φi < 0, i = 1, 2. This automatically gives a
segmentation of the image. However, as we shall see below, this modi¬es
the minimization problem a bit.
As above, the link between the four regions can be made by introducing
F (u+ , u’ , φ) = |u+ ’ u0 |2 H(φ)dx (12.17)
four functions u++ , u+’ , u’+ , and u++ in an obvious fashion:
± ++
|u’ ’ u0 |2 (1 ’ H(φ))dx
+ u (x) if φ1 (x), φ2 (x) > 0,

 +’
u (x) if φ (x) > 0 > φ (x),
1 2
u(x) = (12.23)
| u+ |2 H(φ)dx
+ν ’+
u (x) if φ1 (x) < 0 < φ2 (x),

 ’’
u (x) if φ1 (x), φ2 (x) < 0.
| u’ |2 (1 ’ H(φ))dx

This gives us
| H(φ)|dx.

u =u++ H(φ1 )H(φ2 ) + u+’ H(φ1 )(1 ’ H(φ2 )) (12.24)
+ u’+ (1 ’ H(φ1 ))H(φ2 ) + u’’ (1 ’ H(φ1 )(1 ’ H(φ2 )).
The Euler-Lagrange equations for u+ and u’ are
The energy in level set formulation based on the Mumford-Shah functional
u+ ’ u0 = ν∆u+ on {x | φ(x) > 0} (12.18)
= 0 on {x | φ(x) = 0} (12.19)
|u++ ’ u0 |2 H(φ1 )H(φ2 ) dx
‚n F (u, φ) = (12.25)
u’ ’ u0 = ν∆u’ on {x | φ(x) < 0} (12.20)
| u++ |2 H(φ1 )H(φ2 ) dx
‚u’ +ν
= 0 on {x | φ(x) = 0}. (12.21)
|u+’ ’ u0 |2 H(φ1 )(1 ’ H(φ2 )) dx
These two sets of elliptic boundary value problems will have a smoothing
and denoising e¬ect on the image, but only inside homogeneous regions, not
| u+’ |2 H(φ1 )(1 ’ H(φ2 )) dx

across edges. The Euler-Lagrange equations for φ, using gradient descent
12.4. Extensions 127 128 12. Snakes, Active Contours, and Segmentation

|u’+ ’ u0 |2 (1 ’ H(φ1 ))H(φ2 ) dx

| u’+ |2 (1 ’ H(φ1 ))H(φ2 ) dx

|u’’ ’ u0 |2 (1 ’ H(φ1 ))(1 ’ H(φ2 )) dx

| u’’ |2 (1 ’ H(φ1 ))(1 ’ H(φ2 )) dx

| H(φ1 )| dx + µ | H(φ2 )| dx.

As the authors themselves note in [37], the last expression | H(φ1 )| dx+
| H(φ2 )| dx is not the length of the free boundary. However, it is certainly
between one and two times that quantity. Some segments are counted once,
some twice. However, this releases the Mumford-Shah functional from the
well-known restriction that only 120—¦ -angle functions are possible, within
the class of multiple junctions. If one wishes to minimize the precise term
proportional to length in a multiphase problem, one can use the technique
involving constraints and more level set functions that was introduced by
Zhao et al. [175]. For the segmentation active contours problem described
here, this technique involving 2 (or [log2 n] if one does not use the four color
result and n is the number of separate regions desired) level set functions
seems to work quite well for the piecewise constant case (see [36]).
The Euler-Lagrange equations for the four u functions are as in the two-
phase case which means that they decouple. The time-dependent coupled
gradient descent equations for φ1 and φ2 are easily solved with very simple
changes over the two-phase, one-φ case.
In Figures 12.6 and 12.7 the vector-valued active contour model from
[32] is used, where objects are recovered from combined channels with
missing information in each channel. In Figures 12.9, 12.10, and 12.11,
the piecewise-constant four-phase segmentation model from [36] is used, as
a particular case of the piecewise-smooth four-phase model from [37]. A
similar result from [36] is shown in Figure 12.12, where triple junctions are
also detected in a color image, with the piecewise-constant segmentation
model using three level set functions.

Figure 12.2. Europe night-lights [34].
12.4. Extensions 129 130 12. Snakes, Active Contours, and Segmentation

Figure 12.4. Evolution of an active surface using the 3D version of the active
contour without edges from [33] on volumetric MRI brain data. We show here
only a 61 — 61 — 61 cube from the 3D calculations performed on a larger domain
containing the brain.

Figure 12.3. Detection of the contours of a plane in a noisy environment [34].
12.4. Extensions 131 132 12. Snakes, Active Contours, and Segmentation

Figure 12.5. Cross-sections of the previous 3D calculations showing the evolv-
ing contour and the ¬nal segmentation on a slice of the volumetric image. We
illustrate here how interior boundaries are automatically detected.

Figure 12.6. Numerical results using the multichannel version of the active con-
tour model without edges (from [32]) to detect the full contour of an airplane from
two channels. Note that channel 1 has an occlusion, while channel 2 is noisy.
12.4. Extensions 133 134 12. Snakes, Active Contours, and Segmentation

Figure 12.7. Color image, its gray-level version, and the three RGB channels. (See
also color ¬gure, Plate 3.)

Figure 12.8. Recovered objects without well-de¬ned boundaries, using the multi-
channel version of the active contour model without edges from [32]. The three
objects could not be recovered using only one channel or the intensity image.
(See also color ¬gure, Plate 4.)
12.4. Extensions 135 136 12. Snakes, Active Contours, and Segmentation

Figure 12.10. Evolution of the four-phase segmentation model using two level
set functions. Left: the evolving curves. Right: corresponding piecewise-constant
Figure 12.9. Original and segmented images (top row); ¬nal segments (second
segmentations. Initially, we seed the image with small circles to speed up the
and third rows) [36].
numerical calculation [36].
12.4. Extensions 137 138 12. Snakes, Active Contours, and Segmentation

Figure 12.11. Segmentation of an outdoor picture using two level set functions
and four phases. In the bottom row we show the four segments obtained [36].
Figure 12.12. Color picture with junctions. Three level set functions representing
up to eight regions. Six segments are detected. We show the ¬nal zero level sets
[36]. (See also color ¬gure, Plate 5.)
140 13. Reconstruction of Surfaces from Unorganized Data Points

and Voronoi diagrams. Although this approach is more versatile in that
it can deal with more general data sets, the constructed surface is only
piecewise linear, and it is di¬cult to handle nonuniform and noisy data.
Furthermore, the tracking of large deformations and topological changes
13 can be di¬cult using explicit surfaces.
Recently, implicit surfaces, or volumetric representations, have attracted
signi¬cant attention. The traditional approach (see Bloomenthal et al. [16])
Reconstruction of Surfaces from uses a combination of smooth basis function primitives such as blobs to
Unorganized Data Points ¬nd a scalar function such that all data points are close to an isocontour
of that scalar function. This isocontour represents the constructed implicit
surface. However, computation costs are very high for large data sets, since
the construction is global, which results in solving a large linear system;
i.e., the basis functions are coupled together, and a single data point change
can result in globally di¬erent coe¬cients. This makes human interaction,
incremental updates, and deformation di¬cult. However, recently, Carr et
al. [26] used polyharmonic radial basis functions (RBF) to model large
data sets by a single RBF. The key new idea here is the use of the Fast
Multipole Method (FMM) of Greengard and Rokhlin [76] to greatly reduce
the storage and computational costs of the method. Another crucial idea is
the use of o¬-surface points on both sides of the point cloud. However, the
ability to interpolate curves and surface patches, the ability to do dynamic
deformation, the performance on coarse data sets, and the speed of the
13.1 Introduction method all seem to be inferior to the method we describe below. On the
other hand, the method proposed by [26] does give an analytic, grid-free
Surface reconstruction from an unorganized data set is very challenging.
expression and exact control of the ¬ltering error.
The problem is ill-posed, i.e., there is no unique solution. Furthermore, the
Zhao et al. [177, 176] proposed a new weighted minimal surface model
ordering or connectivity of the data set and the topology of the real surface
based on variational formulations. Only the unsigned distance function to
can be rather complicated. A desirable reconstruction procedure should be
the data set is used, and the reconstructed surface is smoother than piece-
able to deal with complicated topology and geometry as well as noise and
wise linear. The formulation is a regularization that is adaptive to the local
nonuniformity of the data to construct a surface that is a good approxi-
sampling density that can keep sharp features if a local sampling condi-
mation of the data set and has some smoothness (regularity). Moreover,
tion is satis¬ed. The method handles noisy as well as nonuniform data and
the reconstructed surface should have a representation and data structure
works well in three spatial dimensions.
that is not only good for static rendering but also good for deformation,
animation, and other dynamic operations on surfaces.
For parametric surfaces such as NURBS (see Peigl and Tiller [128] or
13.2 The Basic Model
Rogers [137]), the reconstructed surface is smooth, and the data set can
be nonuniform. However, this requires one to parameterize the data set
in a nice way such that the reconstructed surface is a graph in the pa- Let S denote a general data set, which can include data points, curves, and
rameter space. The parameterization and patching can be di¬cult for pieces of surfaces. De¬ne d(x) = dist(x, S) to be the distance function to S.
surface reconstruction from an arbitrary data set. Also, noise in the data is The following surface energy is de¬ned for the variational formulation:
di¬cult to deal with. Another popular approach is to reconstruct a trian- 1
gulated surface using Delaunay triangulations and Voronoi diagrams. The p
1 ¤ p ¤ ∞,
E(“) = d (x) ds , (13.1)
reconstructed surface is typically a subset of the faces of the Delaunay tri- “
angulations. A lot of work has been done along these lines (see, for example,
where “ is an arbitrary surface and ds is the surface area. The energy
Amenta and Bern [6], Boissonat and Cazals [17], and Edelsbrunner [58]),
functional is independent of parameterization and is invariant under ro-
and e¬cient algorithms are available to compute Delaunay triangulations
13.2. The Basic Model 141 142 13. Reconstruction of Surfaces from Unorganized Data Points

13.3 The Convection Model
tation and translation. When p = ∞, E(“) is the value of the distance
of the point x on “ furthest from S. We take the local minimizer of our
energy functional, which mimics a weighted minimal surface or an elastic The evolution equation (13.2) involves the mean curvature of the surface,
membrane attached to the data set, to be the reconstructed surface. and it is a nonlinear parabolic equation. A time-implicit scheme is not
The gradient ¬‚ow of the energy functional in equation (13.1) is currently available. A stable time-explicit scheme requires a restrictive-time
step size, ∆t = O( x2 ). Thus it is desirable to have an e¬cient algorithm
p ’1
d“ 1 to ¬nd a good approximation before we start the gradient ¬‚ow for the
dp (x) ds dp’1 (x) (13.2)
=’ d(x)· N + d(x)κ N ,
dt p minimal surface. We propose the following physically motivated convection

model for this purpose.
and the minimizer or steady-state solution of the gradient ¬‚ow satis¬es the
If a velocity ¬eld is created by a potential ¬eld F , then v = ’ F . In
Euler-Lagrange equation
our convection model the potential ¬eld is the distance function d(x) to
1 the data set S. This leads to the convection equation
dp’1 (x) d(x) · N + d(x)κ = 0. (13.3)
p ‚φ
= d(x) · φ. (13.6)
We see a balance between the attraction d(x) · N and the surface tension
For a general data set S, a particle will be attracted to its closest point
d(x)κ. Moreover, the nonlinear regularization due to surface tension has a
in S unless the particle is located an equal distance from two or more data
desirable scaling d(x). Thus the reconstructed surface is more ¬‚exible in the
points. The set of equal distance points has measure zero. Similarly, points
region where the sampling density is high and is more rigid in the region
on our surface, except those equal distance points, are attracted by their
where the sampling density is low. We start with an initial surface that
closest points in the data set. The ambiguity at those equal distance points
encloses all the data and follow the gradient ¬‚ow in equation (13.2). When
is resolved by adding a small surface tension force, which automatically
p = 1, the surface energy de¬ned in equation (13.1) has the dimension
exists as numerical viscosity in our ¬nite di¬erence schemes. Those equal
of volume and the gradient ¬‚ow in equation (13.2) is scale-invariant. In
distance points on the curve or surface are dragged by their neighbors, and
practice we ¬nd that p = 1 is a good choice.
the whole curve or surface is attracted to the data set until it reaches a
We use the same motion law for all level sets of the level set function,
local equilibrium, which is a polygon or polyhedron whose vertices belong
which results in a morphological partial di¬erential equation. The level set
to the data set as the viscosity tends to zero.
formulation becomes
The convection equation can be solved using a time step ∆t = O( x),
‚φ φ φ leading to signi¬cant computational savings over typical parabolic ∆t =
d· · | φ|.
= +d (13.4)
| φ| | φ|
‚t O( x2 ) time-step restrictions. The convection model by itself very often

<< . .

( 23)

. . >>

Copyright Design by: Sunlight webdesign