LINEBURG

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

spontaneously.

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-

| φ|

‚t

age over „¦c (which is an unstable equilibrium) or if φ ≡ 0, which is the

(12.10)

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

tan’1

δ (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

| φ|

‚t

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

‚φ φ

N

’ |u+ ’ u0 |2 + |u’ ’ u0 |2

‚φ φ 1 ·

= δ (φ) µ (12.22)

»+ (u0,i +

C i )2

· ’ ’

= δ (φ) µ (12.15) | φ|

‚t

i

| φ|

‚ N i=1

’ ν| u+ |2 + ν| u’ |2 .

N

1

»’ (u0,i ’ Ci )2 .

’

+ i

N The new numerical challenge is to obtain the numerical solution of the set

i=1

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

results.

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

where

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)

is

‚u+

= 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)

‚n

|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

p

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

1

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)

‚t

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

Copyright Design by: Sunlight webdesign