(PDF) Isophotes Selection and Reaction-Diffusion Model for Object Boundaries Estimation

Isophotes Selection and Reaction-Diffusion Modelfor Object Boundaries Estimation

CHARLES KERVRANN, MARK HOEBEKE AND ALAIN TRUBUILINRA—Biometrie, Domaine de Vilvert, 78352 Jouy-en-Josas, France

[emailprotected]

Received February 16, 2001; Revised May 13, 2002; Accepted May 13, 2002

Abstract. This paper investigates generic region-based segmentation schemes using area-minimization constraintand background modeling, and develops a computationally efficient framework based on level lines selection coupledwith biased anisotropic diffusion. A common approach to image segmentation is to construct a cost functionwhose minima yield the segmented image. This is generally achieved by competition of two terms in the costfunction, one that punishes deviations from the original image and another that acts as a regularization term. Wepropose a variational framework for characterizing global minimizers of a particular segmentation energy that cangenerates irregular object boundaries in image segmentation. Our motivation comes from the observation that energyfunctionals are traditionally complex, for which it is usually difficult to precise global minimizers correspondingto “best” segmentations. In this paper, we prove that the set of curves that minimizes the basic energy model underconcern is a subset of level lines or isophotes, i.e. the boundaries of image level sets. The connections of ourapproach with region-growing techniques, snakes and geodesic active contours are also discussed. Moreover, it isabsolutely necessary to regularize isophotes delimiting object boundaries and to determine piecewise smooth orconstant approximations of the image data inside the objects boundaries for vizualization and pattern recognitionpurposes. Thus, we have constructed a reaction-diffusion process based on the Perona-Malik anisotropic diffusionequation. In particular, a reaction term has been added to force the solution to remain close to the data inside objectboundaries and to be constant in non-informative regions, that is the background region. In the overall approach,diffusion requires the design of the background and foreground regions obtained by segmentation, and segmentationof the adaptively smoothed image is performed after each iteration of the diffusion process. From an applicationpoint of view, the sound initialization-free algorithm is shown to perform well in a variety of imaging contextswith variable texture, noise and lighting conditions, including optical imaging, medical imaging and meteorologicalimaging. Depending on the context, it yields either a reliable segmentation or a good pre-segmentation that can beused as initialization for more sophisticated, application-dependent segmentation models.

Keywords: grouping and segmentation, energy minimization, level sets, level lines, isophotes, connectedcomponents, anisotropic diffusion

1. Introduction

One of the primary goals of early vision is to segmentthe domain of an image into regions ideally correspond-ing to distinct physical objects in the scene. While it hasbeen clear that image segmentation is a critical problem

in a wide range of computer vision applications, it hasproven difficult to express segmentation criteria thatcapture non-local properties of an image and to developefficient algorithms for computing segmentations. Thisis especially true if the images were obtained underdifferent conditions and have different content. This

64 Kervrann, Hoebeke and Trubuil

work is concerned with the segmentation of salientstructures from natural images by integrating two ap-proaches based on the minimization of a segmentationfunctional and the technique of anisotropic diffusion,to form an overall approach that is robust to noise andpoor initialization.

In general, the goal of image segmentation is to findgroups which are both hom*ogeneous in the same groupand well separated (Deriche, 1987; Pavlidis and Liow,1990; Chu and Agarwal, 1993; Geman and Geman,1984; Blake and Zisserman, 1987; Mumford and Shah,1989; Leclerc, 1989). Many of segmentation tech-niques rely on the design and minimization of an en-ergy functional which captures the interaction betweenmodels and image data. We note that one of the mostwidely studied mathematical models in image process-ing is the active contour model (or “snakes”) to ex-tract the boundaries of hom*ogeneous regions withinan image (Kass et al., 1987). In practical imaging, theinteractive initialization of the energy-based active con-tour near the desired boundary significantly reduces thedifficulty of segmentation. On the other hand, becausethey rely on the image gradients, these models can failin presence of strong noise and may be sensitive to thestarting position. In addition, most algorithms optimiz-ing the energy function associated to active contours,find only local minima, and thus have no measure of thesignificance of the extracted boundary for the image asa whole. Several improvements have been proposed tosupport several levels of user initialization/interaction:a “balloon force” which controls the interior area ofthe curve has been introduced in Cohen (1991) to passover weak edges and extract salient edges; topologi-cally adaptive active contour methods involve solvingthe energy-based active contours minimization by thecomputation of geodesics (Kichenesamy et al., 1996;Caselles et al., 1997). In this approach, a curve is em-bedded as a zero level set of higher dimensional surface(Osher and Sethian, 1988; Sethian, 1996). The entiresurface is evolved to minimize an intrinsic Riemannianmetric depending on image gradient. These models aremore adaptable to detect an arbitrary number of objectsand less sensitive to initialization but use unreliable lo-cal information to make a hard premature decision.

In other respects, several works were conductedin the field of regions-based approaches which aremore suitable for unsupervised image segmentation.Segmentation consists in finding partitions of the im-age pixels into zones the most hom*ogeneous possi-ble corresponding to coherent image properties such

as brightness, color and texture. hom*ogeneity may bemeasured by a given global objective function and harddecisions are made only when information from thewhole image is examined at the same time. Thus, pastapproaches have centered on formulating the problemas the minimization of an energy functional involvingthe image intensity and edge functions. Some energymodels are based on a discrete model of the image, suchas Markov random fields (Geman and Geman, 1984;Blake and Zisserman, 1987; Leclerc, 1989; Wang,1998) or markovian deformable templates (Grenanderand Miller, 1994) whereas variational models are basedon a continuous model of the image (Mumford andShah, 1989; Geiger and Yuille, 1991; Mumford, 1994;Morel and Solimini, 1994; Schnorr, 1998; Jermyn andIshikawa, 1999). Finally, Blake and Zisserman (1987)and Mumford and Shah (1989) have written about mostaspects of this approach to segmentation and have pro-posed various complex functionals whose minima cor-respond to segmented images. In a recent review, Moreland Solimini (1994) have, indeed, shown that most ap-proaches aim at optimizing a cost functional which isthe combination of three terms: one which ensures thatthe smoothed image approximates the observed one,another which states that the gradient of the smoothedimage should be small, except on a discontinuity set,and a last one which ensures that the discontinuity sethas a small length. Such underlying a priori constraintsencourage the emergence of few regions whose bound-aries are regular (Mumford and Shah, 1989; Leclerc,1989). The integration of boundary and region informa-tion sources is also currently achieved through a single-objective function (Zhu and Yuille, 1996; Paragios andDeriche, 2000). These functionals include both localand global information, adding robustness to noise andweak boundaries.

In other respects, while these different approachesoffer powerful theoretical frameworks and minimizersexist (Morel and Solimini, 1994; Schnorr, 1998), it isoften computationally difficult to optimize the associ-ated functionals. Depending on the formulation of thecost function with either continuous or discrete vari-ables, different minimization techniques can be ap-plied. For cost functions of continuous variables, itis in many cases possible to compute a set of par-tial differential equations, the Euler Lagrange equa-tions, that are solved by the minima. For cost func-tions of discrete variables one must generally rely ondirect minimization techniques. Often graduate-non-convexity (Blake and Zisserman, 1987) or Monte Carlo

Isophotes Selection and Reaction-Diffusion Model 65

Markov Chain algorithms such as the well knownsimulated annealing (Geman and Geman, 1984) arethe methods of choice. Typically, they estimate thecurves that maximally separate unknown statistics in-side and outside the curves and avoid bad local min-ima of cost functionals (Blake and Zisserman, 1987;Leclerc, 1989; Grenander and Miller, 1994). Unlikeprevious segmentation techniques, other methods usespecific a priori knowledge to ease the segmentationtask: they may assume that the number of objects(Chesnaud et al., 2000; Figueiredo et al., 2000), lay-ers (Darrell and Pentland, 1995), classes (Paragiosand Deriche, 2000; Yezzi et al., 1999), or the statis-tics inside region boundaries (Grenander and Miller,1994; Samson et al., 1999; Chan and Vese, 1999;Amadieu et al., 1999) are known or estimated us-ing Expectation-Maximization procedures or ad-hocmethods (Darrell and Pentland, 1995; Paragios andDeriche, 2000; Yezzi et al., 1999). The region bound-aries propagation can be then implemented using thelevel set theory (Osher and Sethian, 1988; Sethian,1996; Kichenesamy et al., 1996) but these supervisedsegmentation methods may be sensitive to the initialcurve conditions (Amadieu et al., 1999; Samson et al.,1999; Chan and Vese, 1999; Paragios and Deriche,2000).

In this paper, we present a general approach forimage segmentation that incorporates edge-preservingsmoothing and energy-based regions extraction. Notmuch research has been done which analyzes thecombination of non-linear diffusion methods withsegmentation techniques. However, extracting signifi-cant connected regions from a piecewise smooth imageestimated using non-linear diffusion methods (Peronaand Malik, 1990; Black and Rangarajan, 1996) still re-quires considerable effort. In order to create a true seg-mentation, the regularized image must be processed bya ad-hoc method leading to a partitioning of the entireimage into a finite number of regions. This topic hasbeen already addressed in Weickert (1998b) by combin-ing partial differential equations (PDE)-based regular-ization methods with a watershed algorithm (Vincentand Soille, 1991): the algorithm estimates a piecewisealmost constant image, and performs the partitioning ofthe entire image afterwards. In our work, we shall seethat regions can be extracted at each discretized timeof the diffusion equation as connected components oflevel sets of the PDE-based regularized image. From amethodology point of view, the main contributions ofthe paper are the following:

• We aim at characterizing the minimizers of a ba-sic energy functional combining area-minimizationconstraint with Gaussian priors for luminance dis-tributions over object and background regions, asa subset image level sets. The energy model underconcern, introduced in a discrete setting by Beaulieuand Goldberg (1989) and reviewed by Morel andSolimini (1994), allows to partition the image intoregions, though in a more restrictive manner thanwidely studied approaches (Mumford and Shah,1989; Blake and Zisserman, 1987; Leclerc, 1989;Geiger and Yuille, 1991; Zhu and Yuille, 1996)or more recent approaches (Chan and Vese, 1999;Yezzi et al., 1999; Samson et al., 1999; Paragios andDeriche, 2000) since it can generate irregular objectboundaries. The studied energy model is composedof two terms, one that punishes deviations fromthe image and another that encourages the emer-gence of a small number of regions with irregularboundaries. Nevertheless, unlike related classifica-tion works (Chan and Vese, 1999; Yezzi et al., 1999;Samson et al., 1999; Paragios and Deriche, 2000), thestatistics or the number of classes are not estimatedbefore performing the segmentation: the works pro-posed by Chan and Vese (1999) and Yezzi et al.(1999) with an implementation based on the levelset formulation concerns respectively the bi-modaland tri-modal cases. More generally, a characteriza-tion of minimizers of energy functionals involving aconstant-piecewise model to describe images are notalways available; a fairly complete analysis is avail-able only for a simplified version of the Mumfordand Shah model that approximates a given imagewith piecewise constant functions (Koepfler et al.,1994; Morel and Solimini, 1994). The present inves-tigation is based on a variational approach. We provehereafter that the set of curves that minimizes the en-ergy model under concern is a subset of level linesor isophotes defined from level sets of the image.In other words, minimizing the energy is equivalentto select a subset of connected components corre-sponding to image regions that optimally partitionthe image domain.

• An additional difference with previous related workslies in the algorithmic implementation. Instead ofa commonly used level set formulation (Chan andVese, 1999; Yezzi et al., 1999; Paragios and Deriche,2000), we derive an efficient and initialization-freealgorithm, based on incremental level lines selectionover uniformly quantized version of original data,

66 Kervrann, Hoebeke and Trubuil

which yields a sub-optimal solution to the globalminimizer.

• In some respects, it is absolutely necessary to reg-ularize object boundaries and to determine piece-wise smooth or constant approximations of the imagedata inside the boundaries of objects for vizualiza-tion and pattern recognition purposes. Hence, wehave coupled the level line selection procedure witha robust anisotropic diffusion scheme with sepa-rate object/background processing and robust criticalscale estimation to provide a regularized segmenta-tion. The spatially varying “edge-stopping” param-eter is herein influenced by the segmentation of theforeground (composed of objects) and backgroundregions, at each discretized time of the diffusionprocess.

The rest of the paper is organized as follows: inSection 2, we describe the energy-based model forimage segmentation and we outline its general fea-tures and minimizers. In this section, we describea numerical implementation of the initialization-freesegmentation algorithm. In Section 3, we present exper-imental results obtained by applying the isophotes se-lection algorithm to synthetic as well as real images. InSection 4, we propose a robust statistical measure of thegradient variation and use this in an adaptive reaction-diffusion framework to regularize level lines delimitingobject boundaries. The qualitative and quantitative as-sessment of the accuracy of the combined approachis provided in Section 5. Conclusions are presented inSection 6.

2. Isophotes Selection for Image Segmentation

Our model is the minimizer of a segmentation energy.Most of the time, the energy is designed as a combina-tion of several terms, each of them corresponding to aprecise property which much be satisfied by the optimalsolution. These energy models usually have two parts: adata model and a prior model also called “regularizer”since it was initially conceived to make the problem ofminimizing the data model well-posed (in the sense ofHadamard).

2.1. Description of the Model

Our theoretical setting is the following. Let S be anopen subset of R

2 and f a grey-scale image treatedas a function defined on S. In practical imaging S

is a collection of pixels within a discretized rectan-gle, and possible values of f are given by integers[0, 256[ ∩ N, so both the domain and the range of fare discrete sets. Below we will work in the contin-uous setup, where S is a subset of a Euclidian spaceand f : S → R represents the observed data func-tion. The continuous setup allows us to refer to ana-lytic tools, while leaving always a possibility to “dis-cretize” the problem. We use the notation f to denoteboth the function and the corresponding grey-scaleimage. We use the terminology “site” or “pixel” todenote a point of the image, even in the continu-ous case. Each point x ∈ S is assigned a grey valuef (x). Without further notice we assume that f is up-per semi-continuous. This ensures the measurabilityof f . According to Matheron (1975), we interpretthe image f as a nested family of sets defined byLγ ( f ) = {x ∈ S : f (x) ≥ γ }, γ ∈ R. Each level setLγ ( f ) is assumed to be of finite perimeter (Caselleset al., 1999).

We define the solution to the segmentation problemas the global minimum of a regularized criterion overall regions. Let �i ⊂ S, i = 1, . . . , P , be a non-emptyimage domain or object and ∂�i its boundary. A par-titioning of S consists in finding a set {�i }P

i=1 and abackground � defined as the complementary subset ofthe union of objects (see Fig. 1):

S\� =P⋃

i=1

�i , �i

⋂i �= j

� j = ∅ and

(1)�i

⋂� = ∅.

We seek a strong segmentation, that is, a parti-tion of the rectangle S into a finite set of patches,each of which corresponding to a part of the im-age where f is approximated by a constant value.Moreover, we just wish to control the number of re-gions in the image. The regularization theory leadsto associate with the unknown domains �i and� the following regularized objective function, in-spired from Beaulieu and Goldberg (1989) and Istas(1997):

Eλ( f, �1, . . . , �P , �)

=P∑

i=1

∫�i

(f (x) − f �i

)2dx

+∫

�

( f (x) − f �)2 dx + λ

P∑i=1

|�i | (2)

Isophotes Selection and Reaction-Diffusion Model 67

Figure 1. A minimal energy-based partition is given by the set of (P = 2) objects {�2, �3} and a background � = �1 ∪ �4.

where f �i and f � denote respectively the unknownaverage values of f over �i and � = S\⋃P

i=1 �i , |�i |is the two-dimensional Lebesgue measure which corre-sponds to the area of �i since S is defined as a subset ofthe Euclidian space, and λ ≥ 0 is the regularization pa-rameter. We point out that the two first terms in Eq. (2)are coming from the Mumford-Shah energy (Mumfordand Shah, 1989) and the last term has been already usedin active contours (Zhu and Yuille, 1996), balloonsmodels (Cohen, 1991; Zhu and Yuille, 1996; Rougonand Preteux, 1998) and also by Nitzberg and Mumford(1990). We shall see in Section 2.5 that objects canbe extracted as connected components of image levelssets.

Instead of fixing a priori the cardinality of the seg-mentation, which is a highly arbitrary choice, it seemsmore natural to control the emergence of regions byan object area-based penalty weighted by a param-eter λ. The regularization term E p(�1, . . . , �P ) =∑P

i=1 |�i | is herein introduced to encourage the emer-gence of a large background � and gives no con-trol on the local smoothness of boundaries. The con-stant λ can be interpreted as a scale parameter of

the functional that only tunes the number of regions(Beaulieu and Goldberg, 1989; Morel and Solimini,1994). If λ is large, the background is encouragedto be a large part of the image. As λ decreases, alot of regions are allowed and the segmentation isfine. If λ = 0, each point is potentially a region and� = ∅ ; the global minimum value coincides withzero and this segmentation is called the “trivial seg-mentation” (Koepfler et al., 1994; Morel and Solimini,1994).

In this modeling, implicitly, a Gaussian distributionfor the noise is assumed. This means that one observes acorrupted function f = ftrue+ε, where ε is a zero-meanGaussian white noise. The true image ftrue is supposedpiecewise constant:

ftrue(x) =P∑

i=1

f �i 1x∈�i + f �1x∈�, (3)

where, for a set E ⊂ S, the set indicator function1x∈E = 1 if x ∈ E and 1x∈E = 0 otherwise. The stan-dard deviation is assumed to be constant over the entireimage (Zhu and Yuille, 1996).

68 Kervrann, Hoebeke and Trubuil

Our aim is now to define objects in f . Therefore, wedefine the following collection CP of P ≥ 0 admissibleobjects

CP ={

{�1, . . . , �P} ⊂ S are connected;

S \ � =P⋃

i=1

�i ; �i

⋂1≤i �= j≤P

� j = ∅}

(4)

where the subsets {�1, . . . , �P} are the objects of theimage and � is the background. When P = 0, thereis no object in the image. An optimal segmentationof image f is by definition a global minimum of theenergy (when exists)

(�1, . . . , �

P , �

)

= argmin0≤P≤T

argmin{�1,...,�P }∈CP

Eλ( f, �1, . . . , �P , �). (5)

where the collection CT denotes the bank of all ad-missible objects, i.e. CP ⊆ CT , ∀P ≤ T , and T is themaximum number of registered objects in the image.In Fig. 1, T = 4 and a minimal energy-based partitionis given by the set of (P = 2) objects {�2, �3} and abackground � = �1 ∪ �4.

A joint minimization with respect to all unknowndomains �i and parameters f �i is an intricate prob-lem (Blake and Zisserman, 1987; Mumford and Shah,1989; Morel and Solimini, 1994; Zhu and Yuille, 1996;Cohen, 1996), even if T is low, because of the largenumber of possibilities of placing objects inside S andP is unknown. In addition, the set of unknown vari-ables (sets and functions) are, by definition, clearly notindependent. In Section 2.4, we prove that the objectboundaries that minimize (2) are level lines of functionf , which makes the problem tractable.

2.2. Related Works

The Mumford-Shah functional (Mumford and Shah,1989) is known as the most synthetic criterion for thesegmentation process (Morel and Solimini, 1994), de-pending on two variables, the unknown image functionf and the set of unknown boundaries denoted D. Thisfunctional is the combination of three terms: the firstterm ensures that the estimated function f approxi-mates the observed function f ; the second term statesthat the gradient ∇ f is small, except on a discontinu-ity set D; the last term leads to a discontinuity set Dhaving a small length. The simplest energy functional

associated with this generic description of Mumfordand Shah is the restriction to piecewise-constant func-tions (Koepfler et al., 1994; Morel and Solimini, 1994;Zhu and Yuille, 1996):

EM S( f , D) =∫

S\D( f (x) − f (x))2 dx + ν�(D), (6)

where D is a union of boundaries in S with length �(D)and f is piecewise constant on S\D. The constant ν

represents the scale parameter of the functional andmeasures the amount of boundaries: if ν is low a lotof boundaries are allowed; as ν increases, the segmen-tation gets coarser. This functional represents the sim-plest compromise between accuracy of the regions andparcimony of the boundaries. If we fix the boundariesD, then the corresponding minimal f is completelydefined by the fact that its value on each region �i ofS\D is completely defined by the fact that its value oneach region �i of S\D is equal to the average value off �i of f over �i . This corresponds to the “cartoon”model (Zhu and Yuille, 1996)

Ecartoon( f, �1, . . . , �P ) =P∑

i=1

∫�i

(f (x) − f �i

)2dx

+ ν�

(P⋃

i=1

∂�i

). (7)

which also includes the definition of a connected “back-ground” region since, by definition, the union of �i =S. Unlike this model, the “background” region definedin our segmentation energy (2) is intentionally a non-connected region. In addition, from a numerical point ofview, it is not easy to compute a minimizer of (7) sincethe set of boundaries is unknown. We mention the levelset method (Osher and Sethian, 1988; Sethian, 1996)methods and the merging regions method (Koepfleret al., 1994; Morel and Solimini, 1994) to solve theso-called “minimal partition problem”.

We note also that our energy functional (2) is also re-lated to balloons models (Cohen, 1991; Zhu and Yuille,1996; Rougon and Preteux, 1998). Snakes and balloonshave been typically used with creative external forcesto segment various anatomical structures. An activecontour is a closed contour of a region ∂�, definedby ∂�(s) where s can be the arc length of the con-tour. The two-dimensional balloon model extends thesnake energy to include a force λ �N (s) which pushesthe contour out (or in) along its normal �N (s). Chang-ing the sign of λ makes this normal force inflating or

Isophotes Selection and Reaction-Diffusion Model 69

deflating. Contrasting the balloon and snake models,we note that incorporation of the normal force in theballoon model allows the initial position of the contourto be further from the intended final position, while stillenabling convergence. As well, in the balloon model,the initial position can lie either inside or outside the in-tended contour, while the snake model requires the ini-tial position to surround the intended contour (Cohen,1991; Rougon and Preteux, 1998). The additional in-flation or deflation force can be derived from an energyfunctional. It corresponds to the following energy term|�| = λ

∫�

dx which is a reduced form of our priormodel E p(�1, . . . , �P ) = λ

∑Pi=1 |�i | where P = 1.

Thus, the balloon model maximizes its area whilemaximizing the intensity gradient along the contour.Because of the isoperimetric inequality, this results al-most everywhere into a smoothing behavior.

We want also to point out that the approach pre-sented in this paper shares common aspects with theregion competition approach of Zhu and Yuille (1996).The Zhu-Yuille’s algorithm is derived by minimizing ageneralized Bayes/MDL criterion using the variationalprinciple. In this part of the section, we propose to deter-mine the motion for a point xb at the common boundary∂�i (parametrized by s ∈ [0, 1]) of a region �i andthe background � by computing the gradient descenton (2). The time dependent position of the boundary∂�i can be expressed parametrically by xb(s, t). Themotion of the boundary ∂�i is governed by the Euler-Lagrange differential equation. For any point xb(s, t)on the boundary we obtain

dxb(s, t)

dt= −δEλ( f, �i , �)

δxb(s)

= λ �N (xb(s, t)) + [(f (x) − f �i

− ( f (x) − f �)2] �N (xb(s, t)) (8)

where �N (xb(s, t)) is the unit normal to ∂�i at pointxb(s, t). There are two forces acting on the contour, bothpointing along the normal: the first term is analogous toa pressure term; the second term is the statistics force.The better the point xb(s, t) satisfies the hom*ogene-ity requirement, and the weaker is the statistics force.This equation can be seen as a degenerate case of theregion competition algorithm (Zhu and Yuille, 1996).The Euler-Lagrange equations solving for each regioncan be complex and the region competition algorithmfinds a local minima. Notice the Eq. (7) does not holdat T-junctions between domains boundaries ∂�i where

the normal vector is probably ill-defined. The level setframework has been recently investigated for treatingthis issue (Merriman et al., 1994). In particular, us-ing the level-set formulation (Osher and Sethian, 1988;Chan and Vese, 1999; Samson et al., 1999; Paragiosand Deriche, 2000; Yezzi et al., 1999), suitable numer-ical schemes have been derived for solving propagatingequations. However, in both cases, seed regions mustbe provided by the user or randomly put across theimage, and average values f �i and f � are updated ateach step of the iterative algorithm. In this paper, weshall see that the steady solution of a set of P equa-tions of the form (8) can be characterized (P regionssuperimposed on a background �).

2.3. Upper Bound of the Objects Number

As pointed out by Zhu and Yuille (1996), regularizationparameters calibrating energy terms in segmentationmodels can be interpreted in the sense of statistics.1 Itappears, most of the time, that variations in the valuesof these parameters λ have significant effects on thequalitative properties of the minimizer (Younes, 2000).In this section, we show that the maximum number ofobjects is explicitly influenced by λ.

Lemma 1. If there exists an optimal segmentationdefined by (2) and |�i | ≥ |�min|, i = 1, . . . , P, thenthe optimal number P of objects is upper bounded by

Pmax = (λ|�min|)−1∫

S( f (x) − f S)2 dx if

E p(�1, . . . , �P ) =P∑

i=1

|�i |.

Proof of Lemma 1. Suppose that the segmentationenergy is maximum when there is no object in theimage. We have

P∑i=1

|�i | ≤ Eλ( f, �

1, . . . , �P , �

) ≤ Eλ( f, S)

=∫

S( f (x) − f S)2 dx .

If |�i | ≥ |�min|, we have P|�min| ≤ ∑P

i=1 |�i | ≤

λ−1∫

S( f (x) − f S)2 dx . Finally, we obtain

P ≤ (λ|�min|)−1∫

S( f (x) − f S)2 dx .

70 Kervrann, Hoebeke and Trubuil

2.4. Minimizer Description and Isophotes

In this paper, our aim is not to investigate conditionsfor having global minima of the energy under con-cern and discuss their existence, which is quite dif-ficult (Morel and Solimini, 1994). In what follows,we just assume existence of minimizers of the energyEλ( f, �1, . . . , �P , �) among functions of sets of fi-nite perimeter (or of bounded variation (noted BV)).Our estimator is defined by (when exists)

(�1, . . . , �P , ˆ�)

= arg min0≤P≤T

arg min{�1,...,�P }∈CP

Eλ( f, �1, . . . , �P , �) (9)

Hence, we propose the following lemma:

Lemma 2. If there exist minimizers and that infinites-imal variations in the neighborhood of the minimal so-lution do not introduce topological changes, then theset of curves that globally minimizes the energy is asubset of level lines of f :

f|∂�i≡ µi , i = 1, . . . , P.

i.e. the boundary ∂�i of each �i is a boundary of aconnected component of a level set of f .

We defer the Proof of Lemma 2 to the appendix andstudy a situation where the perturbation is connected.Finally, Eq. (35) (see Appendix) states a necessary con-dition which is essential to prove that a subset of levellines is a minimizer of the energy. In this approach,the family of possible partitions, based on level lineswhich do not cross each other’s, includes local andglobal minimizers. Indeed, placing a level curve insideS corresponds to select a local or a global minimumof the energy. A global minimizer (not unique) can becomputed by supervising all combinations of bound-aries of connected components. The rest of the paperis devoted to the computation a subset of level curvesthat minimizes the segmentation energy. In particular,the computing of energies for all possible partitionsmay produce an expensive computational cost. Accord-ingly, we propose a fast stepwise greedy algorithm (seeSection 2.5) to select a suboptimal configuration of ob-jects corresponding to a local minima of the energy. Webasically consider further that a connected componentis an object �i : each bounded connected componenthas a topological border that is composed of edgelscalled lines lines; a level line separates the plane into

two disjoint connected parts, it bounded interior and itsunbounded exterior; �i is comprised in the interior ofone of the level lines.

2.5. A Stepwise Algorithm for Image Segmentation

This section describes our algorithmic procedure forobject boundaries estimation using the result describedabove. We discuss issues that have arisen in convertingthe theory to practice for our applications. We presenta formal description of the method starting with a de-scription of the input parameters of the algorithm. Ourrecommendations for the concrete choice of these pa-rameters as well as the default choices used in our sim-ulations and applications are collected in this section.The algorithm we propose is automatic and does re-quire neither the number of regions nor any initial av-erage values for regions and background.

2.5.1. Level Sets and Object Boundaries. The key in-gredient of the procedure is the construction of objectswhose boundaries are isophotes in the image (Caselleset al., 1999). In practical imaging, both the domain Sand the range of f are discrete sets. We recall thata discrete gray-scale image can be interpreted as anested family of level sets (Matheron, 1975), computedby simple thresholding. In the today’s technology, wecan traditionally associate with an image 256 level sets{Lγ ( f )}, γ ∈ {0, 1, 2, . . . , 255}. Let γ be a fixed levelof the image, 0 ≤ γ ≤ 255, and let uγ be the bi-nary image at level γ of the discretized original imagef , defined by uγ (x) = 1 if f (x) ≥ γ and uγ (x) = 0otherwise. A crude way to build pixels sets correspond-ing to objects would be to proceed to a connected com-ponents labeling (Monasse and Guichard, 2000) of bi-nary images {uγ }, 0 ≤ γ ≤ 255, and to associate eachlabel with an object �i . Its boundary ∂�i would be theborder of the connected component within the imagelevel sets (Caselles et al., 1999).

Instead of computing all the 256 level sets, we restrictthis computation to a small number of L(<256) levelsets and uniformly quantize the image histogram. Now,we consider the scenario where a point x belongs to onesingle connected component at once within the imagelevel sets. We take into account this fact and define thebilevel sets of f with levels v and w, 0 ≤ v ≤ w,as the set of pixels x ∈ S such as v ≤ f (x) ≤ w.For l ∈ N varying from 1 to L , let bl be the binaryimage with bl(x) = 1 if f (x) ∈ [tl−1, tl[ and bl(x) = 0otherwise, where tl is a threshold (tl ∈ [ fmin, fmax]). We

Isophotes Selection and Reaction-Diffusion Model 71

call those images L-bilevel sets of f (Alvarez et al.,1999). Each bilevel image represents a quantizationlevel of the original image. In general, each bilevel set ismade up of n(tl) connected components, where n(tl) isa function of the threshold tl . Notice that the connectedcomponents {�tl ,1, �tl ,2, . . . , �tl ,n(tl )}, 1 ≤ l ≤ L aredisjoint and their union is the image domain S:

l=L⋃l=1

[�tl ,1

⋃�tl ,2

⋃· · ·

⋃�tl ,n(tl )

]= S. (10)

The connected components {�i } of level sets canbe characterized by their surrounding curves {∂�i },that is the level lines (Caselles et al., 1999; Alvarezet al., 1999). If we map these level lines for a given setof L levels, we get a segmentation of the image alsocalled topographic map (Caselles et al., 1999; Froment,2000). Recently, Monasse and Guichard proposed a fastdiscrete algorithm to compute a topographic map usinga contrast invariant tree representation of connectedcomponents (Monasse and Guichard, 2000). As madeclear in Caselles et al. (1999), the topographic mapis the basic structure of the image. More generally,one can consider a segmentation achieved using onlysome connected components of level sets, which is thephilosophy of our approach. The most perceptible levellines can be determined by an isoperimetric criterion(Froment, 2000) or the detection of T-junctions of levellines (Caselles et al., 1999). Both criteria are strongindicators of region boundaries.

Instead, we define herein perceptually significantlevel lines as the level sets boundaries of an uniformlyquantized image by using L quantizers. As a conse-quence, the detection of meaningful level lines willdepend on the quantization parameter L . Unlike pre-vious criteria (Caselles et al., 1999; Froment, 2000),this quantization operation is not invariant to contrastchanges. Nevertheless, we shall see that, in practice,L = {4, . . . , 8} seems sufficient to detect physicallymeaningful objects with large areas in the image.

2.5.2. The Procedure. The proposed algorithm is nota region growing algorithm as described in Koepfleret al. (1994), Morel and Solimini (1994), Beaulieu andGoldberg (1989), and Pavlidis and Liow (1990) sinceall objects are built once and for all. It differs fromthe wathershed approach since regions that emergefrom the watershed segmentation are not necessar-ily connected components within the image level sets(Vincent and Soille, 1991). In addition, the watershed

approach may include an over-segmentation and com-putational expense. Nevertheless, over-segmentationcan be avoided by using non-linear diffusion as a pre-processing step (Weickert, 1998b) and parallel real-time implementations of the watershed algorithm exit.In our approach, we prevent over-segmentation bychoosing a small number L of quantization levels andpost-processing the connected components to removeany components whose surface area |�i | is less thansome threshold |�min| (a parameter of the method).This parameter is commonly used in image segmenta-tion (Salembier and Serra, 1995; Acton and Mukherjee,2000) to eliminate regions corresponding to noise andartifacts in the original image. To implement our levelset image segmentation based on energy minimization,a four stages method is used. Let L , λ, |�min| be theinput parameters set by the user.

1. Bilevel Set Construction. The first step completesa crude mapping of each image pixel on a givenbilevel set. At present, we uniformly quantize thefunction f in L = 4, 8, 16 or 32 equal-sized andnon-overlapping intervals [tl−1, tl[, l = {1, . . . , L}.Given this set of intervals, let bl be the bilevel setimage with bl(x) = 1 if f (x) ∈ [tl−1, tl[ and bl(x) =0 otherwise.

2. Object Extraction. A crude way to build pixels setscorresponding to objects is to proceed to a connectedcomponents labeling of bilevels image bl and to as-sociate each label with an object �i . Though thisprocess may work in the noise-free case, in gen-eral we would also need some smoothing effectof the connected components labeling. So we con-sider a size-oriented morphological operator actingon sets that consists in keeping all connected com-ponents of the output of area larger than a limit|�min|. This area operator does not introduce newfeatures or edges and boundaries of connected com-ponents are preserved (Salembier and Serra, 1995;Acton and Mukherjee, 2000). The list of remainedconnected components then forms the bank CT ofadmissible T objects {�1, . . . , �T } (P ≤ T ) suchas |�i | ≥ |�min|. The connected components of arealower than |�min| are a part of the background �.

3. Configuration Determination. The connected com-ponents are then combined during the third step toform object configurations. For instance, these con-figurations can be built by enumeration of all pos-sible object combinations, i.e. 2T configurations.Each configuration is made of a subset of objectstaken in the bank {�1, . . . , �T } The background

72 Kervrann, Hoebeke and Trubuil

� corresponds to the complementary set of objectsselected for each configuration. Each possible con-figuration can then be indexed by a binary numberbi which is the binary expansion of i(0 ≤ i ≤ 2T ).The binary value of each bit in bi determines thepresence or absence of a given object in the config-uration (see Fig. 1).

4. Energy Computation and Object ConfigurationSelection. Energy calculations take the image in-tensities of the original (not quantized) image toestablish piecewise-constant approximation errors.Energies of the form {∫

�i( f (x)− f �i )

2 dx} are com-puted once and stored in RAM memory. The energyterm

∫�

( f (x) − f �)2 dx is efficiently updated foreach configuration since � is the complementarysubset of the union of objects {�i }P

i=1. The con-figuration that globally minimizes the energy func-tional corresponds to the optimal segmentation. Thetime necessary to perform image segmentation es-sentially depends on the size of the object bank CT ,i.e. the number T of registered connected compo-nents. Nevertheless, all configurations are indepen-dent and could be potentially evaluated on suitableparallel architectures.

2.5.3. Computational Issues. Now we discuss howsome parameters of the procedure can be selected andindicate one possible choice used in our experimentalresults.

2.5.3.1. Connectivity. We define object boundaries asconnected components borders within the image levelsets. Two pixels are said to be connected when theyare neighbors according to some neighborhood rela-tionship. On the discrete domain S (with rectangulartessellation) the neighborhoods of a pixel x are typi-cally defined via 4-connectivity or 8-connectivity.

2.5.3.2. Number of Bilevel Sets. The value of L ismainly determined by the number of meaningful ob-jects that one wishes to extract and the computationaleffort one is able to spend. Decreases L allows to re-duce the number of connected components and pro-cess a small bank of T objects. We propose herein amethod for mapping a set of pixels to a small set oflevels such that each connected component forms arelatively large and meaningful region. We select theoptimal configuration of objects by supervising a smallset of levels (Ishikawa and Geiger, 1998). In practice,our approach successfully segmented various images

into only 4 or 8 levels. Presenting typically level lineswith levels multiple of a fixed amount say fmax− fmin

L willpreserve all edges of an uniformly quantized imagewith L quantizers.

2.5.3.3. Minimal Size of Objects. This area opera-tor affects the image by remaining connected compo-nents within the image level sets that do not satisfy theminimum area criterion (Salembier and Serra, 1995;Acton and Mukherjee, 2000). Boundaries of connectcomponents are not distorted by this operator sinceit does not incorporate any shape of structuring ele-ment on the processed image. We post-process the con-nected components to eliminate patches correspond-ing to noise in the original date. Our default choice is|�min| ∈ [0.0001 − 0.015] × |S|.

2.5.3.4. Hyperparameter λ. The choice of this pa-rameter determines mostly the properties of the seg-mentation result. Increasing this parameter reduces thefinal number of objects to be extracted. If f is a functionfrom S to [0, 255], a default choice for the hyperparam-eter is λ ∈ [0.001−1.]×2552. Of course larger valuesof λ lead to extraction of only one object. In practice,λ = 0.01×2552 provides a reasonable compromise formost cases. However, we keep a possibility to tuningthis parameter in some specific situations depending onwhat is important in each particular case.

2.5.3.5. Energy Minimization. For a fixed bank CT ={�1, . . . , �T } of T objects, one way to choose theoptimal set of objects {�1, . . . , �P}, P ≤ T definedby (9), is to search for all possible combinationsof P objects and compute the corresponding energyEλ( f, �1, . . . , �P , �). Then, by comparing the ener-gies, we can see which collection of objects is the best.Enumerating all possible sets of objects in the objectbank and comparing their energies is computationallytoo expensive if T is large (typically, it is infeasible ifT > 32). Instead of a such a brute force search, onlyused in experiments when T ≤ 20, we propose thefollowing stepwise greedy algorithm for minimizingEλ( f, �1, . . . , �P , �).

We start from P = 0 and introduce one object � j

at a time. Energies of all objects are assumed to bealready stored in a RAM memory. At the first step, wecompute the T energies with one single object � j atonce against the complementary subset � = S\∪T

j �=i=1

�i . Let �1 the estimated object that best lowers Eλ.This object is stored on a RAM memory as an object of

Isophotes Selection and Reaction-Diffusion Model 73

the optimal configuration. It is removed from the initialbank of objects CT . At any steps of the algorithm, a newestimated object is chosen to maximally decrease theenergy Eλ.

Suppose that at the P-th step, P and ˆ� are notknown but we have estimated P objects {�1, . . . , �P}and a current background � = S\{�1, . . . , �P}. LetEλ( f, �1, . . . , �P , �) the current computed energy.Then at the (P + 1)-th step, for each object � j ∈CT \{�1, . . . , �P}, we calculate the following energydifference:

�Eλ( f, � j , �, �◦

) = Eλ( f, �1, . . . , �P , �)

− Eλ( f, �1, . . . , �P , � j , �◦

)

(11)

where � = �◦ ∪ � j . Intuitively, we choose the object

which has the maximal difference, i.e.,

�P+1 = arg max� j ∈CT \{�1,...,�P }

�Eλ( f, � j , �, �◦

(12)

The algorithm stops at P-th step when the addingof any object does not decrease Eλ. This means thatthe optimal number of objects is P = P and the re-mained objects of the bank are a part of the estimatedbackground, i.e. ˆ� = S\{�1, . . . , �P}. In summary,we propose the algorithm of object selection describedin Fig. 2.

This algorithm selects a suboptimal configurationof objects corresponding to a local minima of the en-

Figure 2. Object selection algorithm.

ergy functional. Using this algorithm, T ×(T +1)2 object

configurations are examined at the most, whereas thesupervision of all the configurations correspond to 2T

iterations. The CPU-time taken by our implementationusing the parameter settings is reported in Section 3.

3. Experimental Results in Image Segmentation

This section presents experiments on synthetic imagesas well as real-world images. We are interested in theuse of the technique in the context of optical, medicaland meteorological imagery. Our system successfullysegmented various images into a few regions. For thebulk of the experiments, the algorithm parameters wereset in these experiments as follows: L = 4, 8, 16 or 32,and regions which areas |�i | < |�min| are discarded.The cleaning threshold |�min| ∈ [0.0005 − 0.015]×|S|indicates the minimum area of connected componentsto be retained in the segmentation. Recall that obtain-ing the most significant objects is the goal of this work.For this reason, L was set fairly low in the experi-ments to obtain large regions and to improve robust-ness to noise and artifacts in the image. For our method,λ ∈ [0.001 − 1.] × 2552 varies across the images de-pending on the image content. It is set empirically andvalues that gave visually better results were chosen. Inour experiments, only a few values (3 or 4 values) witha crude precision were tried for this parameter. Thechoice of λ is crucial and a special attention must begranted by the user. Nevertheless, a small perturbationof the chosen parameter does not drastically alter the

74 Kervrann, Hoebeke and Trubuil

segmentation results from experiments. Most segmen-tations took approximately about 5–30 seconds on a296 MHz workstation. Simulations were conducted onsynthetic as well as real-world images to evaluate theperformance of the algorithm. We start by examiningthe influence of the penalization parameter λ on results.

3.1. Influence of the Penalization Parameter λ

Figure 3(a) shows an artificially computed 256 × 256image representing the superposition of two bidi-mensional Gaussian functions located respectively atx0 = (64, 128) and x1 = (160, 128) with variance ofσ0 = 792 and σ1 = 1024.

Figure 3(b) shows the result of the uniform quanti-zation operation applied on Fig. 3(a) (L = 32). Thelevels lines associated with the quantized image are

Figure 3. Segmentation of a synthetic image (L = 32). (a) Original image; (b) uniformly quantized image (L = 32); (c) level lines superimposedon the quantized image; (d) segmentation with λ = 0.01 × 2552; (e) segmentation with λ = 0.1 × 2552 (two objects); and (f) segmentationwith λ = 1. × 2552 (one object).

displayed on Fig. 3(c). Note that level sets of area toosmall are suppressed. Figure 3(d)–(f) show how thepenalization parameter influences the segmentation re-sults. The white borders denote the object boundariesresulting from the segmentation.

Figure 4 illustrates how our method selects the num-ber of segments in the “house” image (256 × 256).The first row shows, from left to right, the originalimage, the topographic map when L = 4 and |�min| =0.0005 ×|S| and the image histogram. The four equal-sized intervals are displayed in Fig. 4(c). In this exper-iment, the maximum number of admissible objects isT = 49. The second row shows the piecewise-constantapproximation of the original image when λ = 0.001×2552, λ = 0.005×2552 and λ = 0.05×2552. The thirdrow shows the corresponding boundary sets. Notice,even if the number of objects is low (T = 49), the

Isophotes Selection and Reaction-Diffusion Model 75

Figure 4. Segmentation of the “house” image (L = 4, |�min| = 0.0005 × |S| pixels, T = 49). (a) Original image; (b) topographic map;(c) image histogram; (d) piecewise-constant approximation using λ = 0.001 × 2552 (P = 49); (e) piecewise-constant approximation usingλ = 0.005 × 2552 (P = 29); (f) piecewise-constant approximation using λ = 0.05 × 2552 (P = 17); (g) boundary set using λ = 0.001 × 2552;(h) boundary set using λ = 0.005 × 2552; and (i) boundary set using λ = 0.05 × 2552.

76 Kervrann, Hoebeke and Trubuil

obtained boundaries in Fig. 4 are relatively numer-ous since the background region is a non-connectedset. A coarser segmentation (P = 17) is obtained inFig. 4(f) and (i) by choosing a higher value for λ.It takes 12 seconds (1176 iterations) of computingtime for building the object bank and selecting thebest configuration using the stepwise greedy algorithm.Enumerating all the configurations is infeasible since2T = 5, 631014 iterations!

3.2. Dynamic Segmentation Process

An example of cloud detection is provided in Figs. 5and 6. For this set of parameters L = 4, λ = 0.1×2552

and |�min| = 0.0001 × |S|, the algorithm labeled seas,continents and small clouds as “background”. The sig-nificant clouds are crudely extracted from the 383×260image and labeled as “objects” (first row in Fig. 6).Figure 5 shows, from left to right, the original image,the topographic map and the piecewise constant ap-proximation of the image if all the T = 207 regionsare used. Figure 6 depicts four samples of the dynamicsegmentation process: the piecewise-constant approx-imation and the boundary sets superimposed on theoriginal image are shown at significant iterations ofthe stepwise greedy algorithm. Figure 6 shows the re-sults at the 214th (P = 2), 2053th (P = 11) and thelast 9781th (P = 55) iteration. The algorithm finallyselects P = 55 regions and stops at the 9781th it-eration (=21 s of CPU time), i.e. before the maximaliteration T ×(T +1)

2 = 21528. For this experiment, nu-merical results of the stepwise greedy algorithm areshown in Table 1. Note that the larger significant re-gions are first extracted before examination of smallclouds.

Figure 5. Level lines of a meteorological image. (a) Original image; (b) topographic map (L = 4, |�min| = 0.0001 × |S|); and (c) piecewise-constant approximation (T = 207).

3.3. Exhaustive Searchvs. Deterministic Minimization

To go further in the comparisons, we can focus onthe results we obtain by minimization of (2) usingthe stepwise greedy algorithm (Fig. 2) to those ob-tained by examining all the object configurations. Theperformance of the minimization procedure is demon-strated for a MR image (256 × 228) where the flatblack background has been previously eliminated dur-ing a pre-processing step. The first row in Fig. 7 shows,from left to right, the original image, the image his-togram and the topographic map. In order to run thealgorithms, we experimentally adjust the parameters:L = 8 and |�min| = 0.0025 × |S| pixels. For this setof parameters, the total number of objects is T = 25.The aim is to segment both the corpus callosum andbrain by selecting λ = 0.5 × 2552. In Fig. 7, the sec-ond row shows the convergence to the global minimaof the energy function using both optimization algo-rithms. In this experiment, the suboptimal stepwisegreedy algorithm converges to a local minima corre-sponding to P = 5 objects after 276 < (T +1)×T

2 iter-ations and 6.4 s of CPU time on a 296 Mhz worksta-tion. We need 395 s to evaluate all the 2T = 33554432configurations and obtain the same segmentation re-sult shown in Fig. 7. In addition, the time necessaryfor building the bilevel sets and extracting the con-necting components is about 5 s in most experiments.We are not claiming that this algorithm always con-verges to the global minima since many toy-examplescan be designed to contradict this fact. However, thisheuristic approach seems adequate in most cases andthen allow us to manipulate many objects in naturalimages.

Isophotes Selection and Reaction-Diffusion Model 77

Figure 6. Segmentation of a meteorological image (λ = 0.1 × 2552). Left column (top to bottom): piecewise-constant approximation atiteration 214 (P = 2), 2053 (P = 11) and 9781 (P = 55). Right column (top to bottom): boundary sets superimposed on the original image atiteration 214, 2053 and 9781.

3.4. Mumford and Shah Functional

In these experiments, we compare our segmentationresults to those provided by a region growing methodusing the simplified version of the Mumford and Shahmodel which aims at approximating a given image withpiecewise constant functions (Koepfler et al., 1994;Morel and Solimini, 1994). The prior model ensures

the discontinuity set has a small length (Mumford andShah, 1989). A satisfying reconstruction using one hun-dred regions is shown in Fig. 8(c) and (d). Our approachproduces comparative visual results using 29 regionsand a non-connected background (see Fig. 4(e)). Incomparison, we run the region-growing algorithm toextract only 24 regions. Figure 8(b) displays the crudelypiecewise-constant approximation result.

78 Kervrann, Hoebeke and Trubuil

Table 1. Object selection at each step of the stepwisegreedy algorithm.

P Iteration Eλ

P∑i=1

|�i |

1 5 8804 158642 214 8360 221853 420 8089 276704 662 7983 290515 832 7928 300826 1028 7902 304187 1229 7887 306768 1449 7877 308389 1630 7868 31026

10 1844 7862 3113711 2053 7855 3123012 2258 7850 3129313 2431 7846 3136214 2620 7843 3141515 2817 7839 3150216 3019 7836 3156917 3203 7833 3162918 3401 7830 3167119 3602 7827 3171320 3783 7825 3174621 3977 7822 3178822 4147 7820 3183323 4352 7818 3219424 4530 7816 3222225 4694 7814 3226026 4881 7812 3229327 5069 7811 3232328 5248 7809 3235229 5420 7808 3237830 5610 7807 3240231 5785 7806 3242432 5960 7805 3244233 6140 7804 3245634 6313 7803 3247335 6481 7802 3249136 6661 7801 3250637 6834 7800 3252038 7003 7799 3253139 7170 7798 3254540 7341 7798 3255641 7508 7797 3256742 7674 7796 3257943 7834 7796 3259144 8007 7795 3261645 8167 7794 3262846 8352 7794 3266347 8490 7793 3267448 8651 7793 3268549 8824 7792 3272050 8969 7792 3273151 9142 7791 3274752 9286 7791 3276553 9457 7791 3278254 9614 7790 3279355 9781 7790 32809

4. Reaction-Diffusion Model

The results of Sections 2 and 3 suggest that the forms ofconnected components in noisy images are not gener-ally smooth. This section deals with object boundaries-preserving isophotes regularization. Numerous partialdifferential equations-based algorithms have been pro-posed recently to tackle problems of noise removal,image enhancement and image restoration (Peronaand Malik, 1990; Nordstrom, 1990; Alvarez et al.,1992; Catte et al., 1992; Alvarez et al., 1993; Alvarezand Mazorra, 1994; Ter Haar Romeny, 1994; Rougonand Preteux, 1995; Charbonnier et al., 1997; Malladiand Sethian, 1996; Kornprobst et al., 1997; Weickert,1998a; Caselles et al., 1998). They have been recentlyinjected into the segmentation problem since they canbe cast as the minimization of an energy functional(Morel and Solimini, 1994; You et al., 1996). With thisvariational approach, edges are traditionally unknownand must be detected at the same time as the object isreconstructed and regularized.

The method we present consists in smoothing theinput image to impose some length on level curves ofthe image. We also intend to keep the object bound-aries, encourage intra-region diffusion and find an over-smoothed representation of the background in the inputimage. It is clear that such requirement will be fulfilledwhen we use an a priori knowledge about objects togenerate a geometry-driven reaction-diffusion schemethat reflects the characteristics of a given spatial po-sition in the image. For that reason we have chosenthe anisotropic diffusion of Perona and Malik (1990)which will not blur or move the edges, where an addi-tional reaction term forces the solution to remain closeto the data in objects and to be a constant in the back-ground. Due to this modeling, the coefficient of spa-tial diffusion will be designed to smooth adaptivelythe foreground, background and edges. Moreover, themethod which consists in moving the iso-intensity con-tours in an image under curvature dependent speedlaws is not appropriate: the mean curvature operator(Alvarez et al., 1992) is free parameter and thereforeenables to easily encourage diffusion inside regions anddiscriminate background and foreground regions overtime. The background region cannot be labeled aftereach iteration of this diffusion process since all the iso-intensity contours are similarly smoothed accordingto curvature. Nevertheless, the classification betweenforeground and background can be performed inde-pendently of the diffusion process in a post-processing

Isophotes Selection and Reaction-Diffusion Model 79

Figure 7. Segmentation results of a MR image with λ = 0.5 × 2552. (a) Original image; (b) image histogram; (c) topographic map (L = 8);(d) segmentation map; and (e) boundary set superimposed on the image.

Figure 8. Segmentation of the “house” image using the simplified Mumford and Shah model (Koepfler et al., 1994; Morel and Solimini, 1994).(a) Piecewise-constant approximation using 29 regions; (b) boundary set of 29 regions; (c) piecewise-constant approximation using 100 regions;and (d) boundary set of 100 regions.

step. In our approach, the luminance statistics insidelevel lines are examined after each iteration to estimatethe coefficients of spatial diffusion inside and outsideregions. When the original image is used as the ini-tial condition, this process produces over time, a set ofsmooth regions, a background with a nearly constantluminance separated by well located and sharp edges.

4.1. Related Works

4.1.1. Perona-Malik Equation. Anisotropic diffusionas proposed by Perona and Malik (1990) is an alter-native to the linear scale-space described by the heatequation. Intuitively, the idea of non-hom*ogeneous dif-fusion is to limit the flow of intensity according to local

80 Kervrann, Hoebeke and Trubuil

gradient information. This would result in the blurringof details within an object while keeping the edgessharp, producing a cartoon of the image. Basically,the idea is to evolve from an original image f0(x), afamily of increasingly smooth images f (x, t) derivedfrom the solution of the following partial differentialequation (Perona and Malik, 1990):

∂ f (x, t)

∂t= div (g(|∇ f (x, t)|)∇ f (x, t)),

∂ f (x, t)

∂n

∣∣∣∣∂S

= 0,

f (x, 0) = f0(x)

(13)

where n denotes the normal to the image boundary ∂S,t is the time or evolution parameter, ∇ the gradient op-erator, div the divergence operator and g(·) the conduc-tance modulating term which is a bounded, positive, de-creasing and tending to zero at infinity as a function of|∇ f |. For this reason g(·) is also called edge-stoppingfunction. Perona and Malik suggested an edge-stoppingfunction g(|∇ f (x, t)|) = (1 + |∇ f (x,t)|2

2σ 2 )−1 which de-pends explicitly on the gradient magnitude of the func-tion itself: Eq. (13) has the effect of limiting blurringnear edges and also increasing the steepness or gradientof edges that are sufficiently steep in the original imagef0. In practice, areas in which the gradient magnitudeis lower than the scale parameter σ will be blurredmore strongly than areas with a higher gradient magni-tude. This tends to smooth noise and unwanted textureregions, while preserving the edges between differentregions. To remove undesirable “staircase” artifacts,we can iterate these equations but, although the noiseis eradicated, the edges are softened.

Let us now consider the simplest discrete approxima-tion of this process. We consider discrete times t ∈ N

and a time step �t . The standard anisotropic diffusionequation is known to present potential numerical in-stability. However this numerical instability is not ob-served by adopting the discrete anisotropic diffusionupdate proposed in Perona and Malik (1990):

f (x, t + 1) = f (x, t) + �t

|Nx |×

∑x−∈Nx

g(δ( f, x, x−, t), σ )

× δ( f, x, x−, t) (14)

where Nx represents the first-order spatial neighbor-hood of x , |Nx | is the number of directions in which dif-fusion is computed (typically |Nx | = 4) and σ is a scale

parameter. The directional derivative δ( f, x, x−, t) (notto be confused with gradient operator ∇ f (x, t)) locatedat x can be defined as a simple difference on the East,North, West and South directions:

δ( f, x, x−, t) = f (x−, t) − f (x, t), x− ∈ Nx . (15)

This discretization, also used in our own framework,has been widely referred to as anisotropic diffusion andwe will refer hereafter to this implementation.

4.1.2. Robust Anisotropic Diffusion. It has beenshown that anisotropic diffusion models derive fromfunctionals of the form (Geiger and Yuille, 1991; Youet al., 1996)

E( f ) =∫

Sρ(|∇ f (x)|, σ ) dx . (16)

E( f ) specifies a generalized membrane defined by anadaptive potential ρ(·) related g(·) by

g(|∇ f (x)|, σ ) = ρ ′(|∇ f (x)|, σ )

|∇ f (x)| . (17)

This form of energy can be interpreted in the senseof robust statistics where the diffusion process aimsat estimating a piecewise smooth image from noisyinput image (Black et al., 1998). In this robust statisticalframework, edges are interpreted as outliers and σ as ascale parameter that will adaptively reduce unwantednoise when the local statistics properties of the imageare Gaussian and will preserve edges when the gradientis high.

In Black et al. (1998), it is established that the func-tion g(·) is proportional to the Lorentzian error normused in robust statistics. Nevertheless, experimental re-sults showed that the Tukey’s biweight ρ-error func-tion produces sharper discontinuities than the originalLorentzian function (Black and Rangarajan, 1996). Wehave also used in our anisotropic smoothing context theedge-stopping function derived from Tukey’s biweightρ-error

g(|∇ f (x)|, σ )

(1 − |∇ f (x)|2

σ 2

if |∇ f (x)| ≤ σ,

0 otherwise.

(18)

The scale parameter σ can be directly estimated in arobust way (see Rousseeuw and Leroy, 1987) from theimage gradients. Intuitively σ should characterize thevariance of the majority of the data within a region.

Isophotes Selection and Reaction-Diffusion Model 81

Black et al. (1998) took the region for computing σ =√5 σ to be the entire image where

σ = 1.4826 median f ||∇ f (x)| − median f (|∇ f (x)|)|.(19)

This follows from the fact that the median value ofthe absolute values of a large enough sample of unit-variance normal distributed one-dimensional values is0.6745 = 1/1.4826. This approach works well whenedges are distributed hom*ogeneously across the imagebut this is rarely the case. The authors recently exploredthe computation of this measure in image patches toprevent the amplification of undesirable texture gradi-ents and noise in hom*ogeneous image regions (Blackand Sapiro, 1999). However, information about the im-age structure and the data generation process can bevaluable to improve tasks such as noise reduction andimage segmentation.

4.1.3. Variational Methods for Denoising Images.Coming from the optimization approach for image de-noising, it is also interesting to add a reaction termwhich asserts that we must remain close to the originalimage (Nordstrom, 1990; Morel and Solimini, 1994;Catte et al., 1992; Charbonnier et al., 1997; Kornprobstet al., 1997; Caselles et al., 1998). Nordstrom has pro-posed to modify purely diffusive models by introducinga bias term constraining the solution to remain close tothe initial condition f0. This reaction term avoids theselection of a stopping time and termination of the dif-fusion at a trivial solution, such as a constant image(Nordstrom, 1990).The biased anisotropic diffusion isdefined by

∂ f (x, t)

∂t= [ f0(x) − f (x, t)]

+ α div(g(|∇ f (x, t)|, σ )∇ f (x, t)),

∂ f (x, t)

∂n

∣∣∣∣∂S

= 0,

f (x, 0) = f0(x).

(20)

Obviously, the left-hand side [ f0(x) − ft (x)] enforcesan additional constraint that penalizes deviation fromthe input image. The intermediate solutions of (20) maybe regarded as a descent on this following energy func-tional (Morel and Solimini, 1994):

E( f ) =∫

S( f − f0)2 dx + α

∫Sρ(|∇ f (x)|, σ ) dx .

(21)

The first summand encourages similarity between therestored image and the original one, while the sec-ond summand rewards smoothness. The smoothnessweight α > 0 is called regularization parameter. Choos-ing a non-convex potential does not guarantee well-posedness and stable algorithms, but (20) however con-verges globally to the steady-state of the regularizationprocess. Notice that computing g(·) using some reg-ularized first order derivatives (e.g. Gaussian deriva-tives) leads to a well-posed PDE, even for non-convexρ-function (Catte et al., 1992).

4.2. Reaction-Diffusion Modelfor Background Simplification

In addition to the classical anisotropic diffusion givenin Section 4.1, other diffusion mechanisms may be em-ployed to adaptively filter an image for image segmen-tation purposes. For instance, the idea of Alvarez andEsclarin (1997) is to propose a reaction-diffusion par-tial differential equation which attracts the image overtime to a finite number of specified grey levels. This ap-proach aims at adaptively quantizing the image using avariational method. The philosophy of our approach issimilar to Alvarez and Esclarin (1997) in a certain pointof view but the methodologies are quite separate. Wefollow hereafter the philosophy of Nordstrom (1990)by incorporating a reaction term to penalize adaptivelythe background and foreground regions. The additionalreaction term forces the solution to remain close to thedata in objects and to be a constant in the background.Hence, one can attempt to solve the following set ofmore sophisticated partial differential equations:

∂ f (x, t)

∂t= [ f0(x) − f (x, t)]

+ α� div(g(|∇ f (x, t)|, σ�(t))

× ∇ f (x, t)), if x ∈ S/�,

∂ f (x, t)

∂t= [ f �(t) − f (x, t)]

+ α� div(g(|∇ f (x, t)|, σ�(t))

× ∇ f (x, t)), if x ∈ �,

∂ f (x, t)

∂n

∣∣∣∣∂S

= 0,

∂ f (x, t)

∂n

∣∣∣∣∂�i

= 0,

f (x, 0) = f0(x)

(22)

where � is the background, � = S\� is the unionof objects, that is the foreground, f �(t) is the average

82 Kervrann, Hoebeke and Trubuil

value of f over � at iteration t , α� and α� are two dis-tinct regularization parameters which control the com-petition between the two energy terms and determinethe smoothing, and σ�(t) and σ�(t) are respectively thescale parameters related to the foreground and back-ground at time t . We note n the normal to either theimage boundary ∂S or the object boundary ∂�i . No-tice the background region is a non-connected regionsince it is defined as the complementary subset of theunion of objects: � = S\⋃P

i=1 �i . It is made up of theunion of several connected regions of any size andf (x) ∈ [ fmin, fmax] if x ∈ �. Accordingly, the varianceσ� region may be large. Unlike the background, the ob-jects �i are connected regions and l fmax− fmin

L ≤ f (x) <

(l+1) fmax− fmin

L if x ∈ �i , i = 1, . . . , P and l = 1, . . . , L .These two foreground/background regions must bethen treated and smoothed differently. Finally, each ob-ject could be also separately processed. However, thecomputational load can be prohibitive if the number Pof objects is large due to the robust but computationallyexpensive estimation of variances (see (19) and (23))used in our framework.

The idea behind the model given by (22) is tocombine diffusion for noise filtering and a reactionterm to either remain close to the initial conditionf0 inside objects or remove structures in the back-ground region. A constant value f �(t) for the back-ground is herein desired since it is regarded as a non-informative region for image interpretation. In our case,we choose a non-convex potential g(·) correspond-ing to the Tukey’s biweight norm. Of course, the non-convexity of the Tukey’s biweight function will givean ill-posed problem. However, we use this functionlater in our experiments because non-convex func-tions leads to edge enhancement and give better re-sults in practice. Notice that finite difference discreteschemes are known to compensate the ill-posedness ofthe continuous model, i.e. introduce additional regu-larization. We refer to Weickert (1998b) for a thoroughdiscussion concerning the implementation of discreteschemes.

In this modeling, the choice of σ� and σ� is thenan important issue and the main concern of this sec-tion. The estimation of these scale parameters differsfrom previous definitions (Black et al., 1998; Black andSapiro, 1999). We propose values with an explicit de-pendence on the spatial coordinates and the data func-tion f . First, we interpret the standard deviations σ�(t)and σ�(t) defined by Rousseeuw and Leroy (1987),Black and Rangarajan (1996), and Black and Sapiro

(1999))

σ�(t) = 1.4826 median f ‖∇ f | − median f (|∇ f |)|,if x ∈ � = S\�,

σ�(t) = 1.4826 median f ‖∇ f | − median f (|∇ f |)|,if x ∈ �,

(23)

as the gradient magnitudes at which outliers begin tobe downweighted in the background and foregroundregions. Hence, we choose a value for the scale pa-rameters so that the Tukey influence function beginsrejecting outliers at the values 2σ�(t) = σ�(t)/

√5 and

2σ�(t) = σ�(t)/√

5 (Black et al., 1998). Thus, start-ing form f (x, 0) = f0(x), the first term of (22) diffusesadaptively the image in the background and foregroundregions, while the reaction terms form either a constantbackground or preserve the original image. Since thecomputation involves determining the background andforeground regions for solving (22), a way for effi-ciency is to estimate a segmentation map at each iter-ation t and update the scale parameters {σ�(t), σ�(t)}and the average intensity f �(t). The set of parame-ters {σ�(t), σ�(t), f �(t)} regulates the filtering froman a priori knowledge of the geometry of the image(composed of objects). In the following, we fix regu-larization parameters α� and α�. High values for α�

suggest that a simplified background is desired. Thischoice of these parameters has of course to dependon the image data and the desired application. Thismodeling improved the existing PDE-schemes: stan-dard anisotropic diffusion are inefficient for filteringstructural distortions and partial non-gaussian clutter(Zhu and Mumford, 1997).

4.3. Discrete Formulation and Algorithmic Scheme

In the previous section, we have introduced a robustmethod to incorporate “fresh” information from thedata structure and thus improving the behavior the dif-fusion process in presence of meaningful objects. Inprevious works (Perona and Malik, 1990; Black et al.,1998; Black and Sapiro, 1999; Weickert, 1998b), theregions for computing the scale parameter are eitherthe entire image or rectangular patches. Besides, imagesmoothing and image segmentation are two indepen-dent tasks (Weickert, 1998b). We found that two dis-tinct adaptive scale parameters for the foreground andbackground are a simple way to control the diffusion

Isophotes Selection and Reaction-Diffusion Model 83

rate since we have

lfmax − fmin

L≤ f (x) < (l + 1)

fmax − fmin

Lif x ∈ �i , i = 1, . . . , P,

fmin ≤ f (x) < fmax if x ∈ �.

(24)

The new conductance function depends not only on thelocal behavior of the data function f which is subjectto noise, but also on the particular characteristics of thesystem at every location x of the image.

In what follows, we consider that, on the outsideof �i , f has a constant value corresponding to thevalue on ∂�i , i.e, the value of a level line. Adoptingthe discrete anisotropic diffusion update in Perona andMalik (1990), diffusion may be implemented by

f (x, t + 1) = f (x, t) + �t α�

|Nx |×

∑x− ∈ Nx

g(δ( f, x, x−, t), σ�) δ( f, x, x−, t)

+ �t[ f0(x) − f (x, t)] if x ∈ S/�

f (x, t + 1) = f (x, t) + �t α�

|Nx |×

∑x− ∈ Nx

g(δ( f, x, x−, t), σ�) δ( f, x, x−, t)

+ �t[ f �(t) − f (x, t)] if x ∈ �

(25)

where g(δ( f, x, x−, t), σ ) ∈ [0, 1] given in (18) indi-cates the presence (g(·) close to 0) or absence (g(·) closeto 1) of discontinuities or outliers inside patches (Blacket al., 1998). Traditionally, the number of iterations inthe discrete implementation is set interactively by vi-sually inspecting the results. The discrete scheme is anaive scheme and probably a more performant schemewould be suitable. Nevertheless, at t = 0, a crude seg-mentation is obtained by selection of a subset of levellines. According to (15) and (25), the directional deriva-tive δ( f, x, x−, t) is computed as:

δ( f, x, x−, t)

f (x−, t) − f (x, t)

if x, x− ∈ �i and x− ∈ Nx ,

if x ∈ �i , x− �∈ �i and x− ∈ Nx ,

where Nx is is the first order-spatial neighborhood ofx . This naive discretization enables to diffuse acrossobject boundaries.

The number of iterations of the diffusion process isindirectly controlled by the segmentation map corre-sponding to a partition of the image into P objects �i :identical level lines/isophotes are invariably selected ifthe restored image tends to a piecewise-constant im-age. In our experiments, this can be accomplished withless than tmax = 50 iterations in sufficient precision. Incontrast to other approaches (Perona and Malik, 1990;Black et al., 1998; Black and Sapiro, 1999; Weickert,1998b), the diffusion stopping rule is then defined withrespect to objects that are being selected and the seg-mentation map. This relationship is more underlinedthrough a practical segmentation algorithm describedin Fig. 9.

5. Experimental Results in Noiseand Clutter Removal

In this section, we compare results of different dif-fusion schemes applied to synthetic data, optic andconfocal images. All images were produced using thesame set of parameters. In our experiments, �t = 0.025and α� = α� = 20 show better performance in form-ing piecewise constant regions and simplifying back-ground. These values have been fixed manually fromexperiments. In addition, α� and α� are not adjusted toselect automatically the stopping time of the diffusionprocess but to regularize the solution; in our experi-ments, the stopping time tmax is heuristically estimated.A subjective examination of examples suggests that thechoice of parameters L , |�min| and λ will depend onthe image that is being analyzed.

We point out the simplified Mumford-Shah func-tional yields to a piecewise-constant approxima-tion/segmentation of the original image. A smoothedand diffused version of the original of the image is notavailable if we consider only this functional. In ourapproach, the overall algorithm yields to, on the onehand a piecewise-constant approximation of the im-age by optimizing (2) along with an identification ofthe non-connected background, and on the other hand,an object-based restoration of the original image. Infact, further comparison with the complete Mumford-Shah model for segmentation/restoration (Mumfordand Shah, 1989) are probably more suitable but theoptimization/implementation of this complex model isquite tricky. Finally, we have chosen to reconstruct ob-jects when the preliminary segmentation (i.e. at timet = 0) is clearly unsatisfying (see Fig. 10 vs. Fig. 4),which is a more challenging task.

84 Kervrann, Hoebeke and Trubuil

Figure 9. Object selection and reaction-diffusion algorithm.

Figure 10. Results of different diffusion schemes after 10 iterations applied to a noisy version (gaussian noise of standard deviation 14.1)of the synthetic image “gdr”. (a) Original image; (b) noisy image; (c) robust anisotropic diffusion (σ = 26.52) (Black et al., 1998); and (d)geometry-driven diffusion-reaction.

5.1. Synthetic Data

In order to assess the performance of the proposedand other common diffusion schemes, we have used

a synthetic 256 × 256 image where a ground truth isknown, coming from the GdR ISIS corrupted with awhite gaussian noise. This image shown in Figs. 10and 12 is composed of six objects against a background

Isophotes Selection and Reaction-Diffusion Model 85

and includes two main difficulties: the irregularity ofthe shape on bottom right-hand side and a grey-levelgradation on top right-hand side. Figure 10 showsthe results after 10 iterations of two diffusion pro-cesses. The first row shows from left to right, theoriginal synthetic “gdr” image, the image after addingnoise (standard deviation 14.1), the image processedwith the robust anisotropic diffusion (Black et al.,1998) (the scale parameter σ = 26.52 is robustly es-timated using (19)), and the image processed with ourgeometry-driven reaction-diffusion scheme. The sec-ond row shows the magnified intensity gradient normcorresponding to each of images on the row above.Both diffusion schemes smooth hom*ogeneous regionsbut geometry-based approach has a better enhancingeffect on the object boundaries after 10 iterations aswe can see on the gradient images. Our method lo-cates them with satisfying visual precision and sharp-ness as shown in Fig. 10. In Fig. 12, the segmentationmaps at iteration t = 1, 2, 5, 10 are displayed. On thisfigure, the background and objects are respectively rep-resented by a “white” label and the average values off computed over each object domain (second row inFig. 12). The third row shows object boundaries repre-sented by “white” solid lines. After the first iteration,we see that object surfaces are not adequately regular-ized and regions are not complete (Fig. 12). We retrievethe meaningful objects after 5 iterations using L = 8and λ = 0.001 × 2552 but the linear gradation on topright-hand side is partially merged in the backgroundwhen t > 10. The histogram of the restored image after10 iterations and the evolution of scale parameters σ�

and σ� as a function of the iteration number are plottedin Fig. 11. The computational time is 20 seconds periteration on a 296 Mhz workstation for the completesegmentation-filtering procedure.

5.2. Clutter Removal

In many applications, distortions in images are notgaussian. For example, the tree branches in Fig. 13may be regarded as clutter that cause occlusions of thebuilding. Modeling such clutter is a challenging prob-lem which can be solved by learning Gibbs distribu-tions for each set of objects, i.e. trees and buildings,and maximizing a posterior distribution using a variantof simulated annealing (Langevin equation) (Zhu andMumford, 1997). This situation can be partially han-dled by the reaction-diffusion filter described in theprevious section. We assume herein that clutter corre-

sponds to non-constant areas in images with respectto hom*ogeneous objects. These inhom*ogeneous areasare then labeled as background and can be eliminatedusing the reaction-diffusion process. However, label-ing background seems a challenging problem in imagesegmentation as we have addressed in Section 2.

Applying our approach to the image (248 × 288)in Fig. 13, that is setting L = 4, λ = 0.025 × 2552

and |�min| = 0.001 × |S| pixels, produces visually ap-pealing constant images with straight edges. As wecan observe on Figs. 13 and 15, trees branches arepartially eliminated from the image and boundariesare enhanced. As a comparison, we run the affinemorphological scale-space and the robust anisotropicdiffusion (constant scale parameter σ = 59.67) pro-cesses in Fig. 15(a) and (b). This corresponding resultshave been introduced just to demonstrate the perfor-mance of this method after 50 iterations when appliedon the same data. To go further in the comparisons,the restored image using the Gibbs reaction-diffusionequations (Zhu and Mumford, 1997) is displayed inFig. 15(d). The underlying image formation model isclearly quite different. However, we aims also at dis-criminating two regions/classes (buildings in the fore-ground and trees in the background). It appears that theforcing term in the model (22) acts to better eliminatebranches without degrading building edges. The imageeventually becomes oversmoothed using the two firstdiffusion schemes with continued iterations whereasour reaction-diffusion process is stopped before 50 it-erations in this experiment. The histogram of the re-stored image after 50 iterations and the evolution ofscale parameters σ� and σ� for increasing number ofiterations are plotted in Fig. 14.

5.3. Confocal Microscopy Imagery

Confocal systems offer the chance to image thick bio-logical tissue in 2D + t or 3D dimensions. They oper-ate in the bright-field and fluorescence modes, allowingthe formation of high-resolution images with a depthof focus sufficiently small that all the detail which isimaged appears in focus and the out-of-focus informa-tion is rejected. Some of the current applications inbiological studies are in neuron research. The biolog-ical study aims at examining the effect of intracellu-lar calcium concentration (Ca2) on neurite outgrowthin individual neuronal cells. Indeed, Ca2 is knownto be implicated as an important regulator of neuriteextension.

86 Kervrann, Hoebeke and Trubuil

Figure 11. Left: histogram of the restored image after 10 iterations. Right: evolution of scale parameters σ� and σ�.

Figure 12. Different stages of the object selection and reaction-diffusion process applied to a noisy version (gaussian noise of standard deviation14.1) of the synthetic image “gdr”. The columns, from left to right, show the diffused images after 1, 2, 5 and 10 iterations respectively. The firstrow show the diffused images. The second row shows the piecewise constant approximation of objects where the background is shown using a“white” label. The third row shows the boundary sets superimposed on the input noisy image.

Isophotes Selection and Reaction-Diffusion Model 87

Figure 13. Results of the geometry-driven reaction-diffusion scheme after 50 iterations applied to a optical image showing a building andbranches. (a) Original image; (b) geometry-driven reaction-diffusion; (c) difference image between (a) and (b); (d) topographic map; (e)piecewise-constant approximation; and (f) boundary set.

Figure 14. Left: histogram of the restored image after 50 iterations. Right: evolution of scale parameters σ� and σ�.

88 Kervrann, Hoebeke and Trubuil

Figure 15. Results of different diffusion schemes applied to a confocal image. First row: original image. Second row, left to right: affine mor-phological scale space, robust anisotropic diffusion, geometry-driven reaction-diffusion, Gibbs reaction-diffusion equations (Zhu and Mumford,1997). Third row: magnitude of the intensity gradient corresponding to each image in the row above.

We have tested the proposed algorithm on a 2D con-focal microscopy 115 × 512 image (Fig. 16, courtesyof INSERM 413 IFRMP n◦23 (Rouen, France)) de-picting a neurite in cultured cerebellar granule cells.In Fig. 16(a), high grey-level values correspond to ele-

vated calcium concentration. Accordingly, our goal istwofold: extracting the axon and nucleus in the fore-ground and segmenting these objects of interest intosignificant sub-regions corresponding to different lev-els of calcium concentration. Additionally, undesirable

Isophotes Selection and Reaction-Diffusion Model 89

Figure 16. Results of the geometry-driven reaction-diffusion scheme applied to a confocal image. (a) Original image; (b) geometry-drivenreaction-diffusion; (c) difference image between (a) and (b); (d) topographic map; (e) piecewise-constant approximation; and (f) boundary set.

structures (occluded axons) in the background shouldbe ideally eliminated when the segmentation is per-formed. The boundaries are quite accurately delin-eated in Fig. 16 (L = 4, λ = 0.05 × 2552, |�min| =0.005 × |S|) after 50 iterations. The level lines (of therestored image) corresponding to L = 4 are shownon Fig. 16(d). On this figure, the background and ob-jects are respectively represented by a “white” label(see second row of Fig. 16). The axon and nucleus

Figure 17. Left: histogram of the restored image after 50 iterations. Right: evolution of scale parameters σ� and σ�.

have been roughly extracted and segmented into sub-regions associated to mainly three classes of calciumconcentration.

This application (Figs. 16 and 18) illustrates the ro-bustness of our method with respect to noise and partialclutter. As we mentioned before, this shows that theproposed method performs reasonably well in pres-ence of non-significant geometric structures. Figure 17displays the histogram of the restored image after 10

90 Kervrann, Hoebeke and Trubuil

Figure 18. Results of the different diffusion schemes applied to a confocal image. Results of the different diffusion schemes applied to aconfocal image. (a) Original image; (b) affine morphological scale space; (c) robust anisotropic diffusion (σ = 62.99); and (d) geometry-drivenreaction-diffusion.

iterations and the decreasing of scale parameters σ�

and σ�. In Fig. 18, if we compare the different diffu-sion processes and observe that boundaries are bettermaintained and enhanced using our approach. Noiseand irrelevant features are eliminated.

6. Conclusion and Perspectives

In this paper, we have presented a level line se-lection approach for extracting structures in images.We proved that the minimizer of our segmenta-tion energy can be directly characterized. We intro-duced a geometry-driven reaction diffusion processto preserve sharper boundaries of objects and regu-larize level lines and object surfaces in the image.A total CPU time of a few seconds for segment-ing a 256 × 256 image on a 296 Mhz workstationmakes the method attractive for many time-criticalapplications. The contribution of this approach has

been illustrated on synthetic as well as real-worldimages.

Several promising directions may be explored forcontinued research. In particular, some open math-ematical questions must be addressed. On the onehand, given the image under consideration, mathe-matical conditions and constraints to ensure that nopathological minima exists, must be investigated; wejust mention the Morse theory which can be the rele-vant functional framework as it has been already in-vestigated in image segmentation by Olsen (1996).On the other hand, the underlying global energyfunctional including segmentation and diffusion isnot actually well-defined due to the dependency of{�i , �, f �i , f �, σ�, σ�} and should be examined infuture works. Moreover, the proposed discrete schemeis a naive scheme and probably a more efficient schemewould be suitable. The influence of parameters �t , α�

and α�, manually adjusted in our experiments, must bealso studied. In other respects, we currently examine

Isophotes Selection and Reaction-Diffusion Model 91

other energy functionals which minimizers can be char-acterized as well, and study an entropy-based quanti-zation technique instead of the uniform quantizationused at present to estimate the objects. Finally, an otherdirection for future work is to extend the proposed ap-proach to operate on multi-spectral 3D images. In thissetting, the structure of the algorithm would be largelythe same, although there are a number of points whichwould need to be examined closely.

Appendix A: Proof of Lemma 2

First we consider the usual distance between two closedsubsets A and B is the Hausdorff distance, defined by

d∞(A, B) = max

{supx∈A

d(x, B), supx∈B

d(x, A)

}, (26)

where d(x, A) denotes, as usually, the distance of apoint x to the set A

d(x, A) = infy∈A

|x − y|. (27)

For two sets A and B such that B ⊆ A, denote∫A\B f

def= ∫A f − ∫

B f.Now, let �δ be a variation of a set � such that

� ⊆ �δ , i.e. the Hausdorff distance d∞(�δ, �) ≤ δ.To prove Lemma 2, we assume that, for any connectedperturbation of � such d∞(�δ, �) ≤ δ, two neigh-boring sets � and �′ do not merge into one single set� ∪ �′ and, for any connected perturbation of � suchd∞(�δ, �) ≤ δ, � does not split into two new sets. Thiscorresponds to prohibited topological changes. Noticethe Morse theory can be the relevant functional frame-work to ensure that no topological changes does occuror not Olsen (1996). Without loss of generality, weprove Lemma 2 for one object � and a background �,that is the closure of the complementary set of �. Then,we have

∫�δ\�

1 def= |�δ| − |�| and

(∫�δ

−(∫

�

= 2∫

�

f∫

�δ\�f +

(∫�δ\�

. (28)

The difference between the involved energiesis equal to �Eλ( f, �, �) = Eλ( f, �δ, �δ) − Eλ

( f, �, �) where

Eλ( f, �, �) =∫

�

( f (x) − f �)2 dx

+∫

�

( f (x) − f �)2 dx + λ |�|

Eλ( f, �δ, �δ) =∫

�δ

(f (x) − f �δ

)2dx

+∫

�δ

(f (x) − f �δ

)2dx + λ |�δ|.

(29)

Therefore, �Eλ( f, �, �) = T1 + T2 + T3 + T4 + T5

with

T1 =∫

�δ

f 2 −∫

�

f 2,

T2 = − 1

|�δ|(∫

�δ

+ 1

|�|(∫

�

T3 =∫

S\�δ

f 2 −∫

S\�f 2, (30)

T4 = − 1

|S| − |�δ|(∫

S\�δ

+ 1

|S| − |�|(∫

S\�f

T5 =∫

�δ\�λ.

Denote �|�| = |�δ|− |�|. Using (28), and passing tothe limit �|�| → 0, i.e. |�δ| � |�|, we obtain (higherorder terms are neglected)

T1 = −T3 =∫

�δ\�f 2,

T2 = − 2

|�|∫

�δ\�f∫

�

f − 1

|�|(∫

�δ\�f

+ 1

|�|2∫

�δ\�1

(∫�

, (31)

T4 = 2

|S| − |�|∫

�δ\�f∫

S\�f − 1

|S| − |�|(∫

�δ\�f

− 1

(|S| − |�|)2

∫�δ\�

1(∫

S\�f

and T5 = λ∫�δ\� 1. Define the following image

moments

m0 =∫

�

1, m1 =∫

�

f, M0 =∫

S1, M1 =

∫S

(32)

92 Kervrann, Hoebeke and Trubuil

Using the mean value theorem for double integral,which states that if f is continuous and a connectedsubset E is bounded by a simple curve, then for somepoint x0 in E we have

∫E f (x)d E = f (x0) · |E | where

|E | denotes the area of E , it follows that

�Eλ( f, �, �) =

a0︷︸︸︷[m2

m20

− (M1 − m1)2

(M0 − m0)2+ λ

]

×∫

�δ\�1 +

a1︷︸︸︷[−2m1

m0+ 2(M1 − m1)

M0 − m0

]f (x0)

×∫

�δ\�1 −

m0+ 1

M0 − m0

]f (x0)2

(∫�δ\�

1)2

(33)

Let xb be a fixed point of the boundary ∂�. Choose �δ

such that ∂�δ = ∂� except on a small neighborhood ofxb. The energy having a minimum for �, f (xb) needsto be solution of the following equation

�Eλ( f, �, �)

�|�| = a0 + a1 f (xb) + O(�|�|) = 0.

(34)

By passing to the limit �|�| → 0, we obtain

a0 + a1 f (xb) = 0. (35)

Equation (35) has one single solution. The coefficientsa0 and a1 do depend on neither xb nor f (xb), and a0 �= 0.The function f is continuous on S ⊂ R

2 and ∂� is aconnected curve. Therefore f (xb) is constant when xb

covers ∂�. In other words, ∂� is a level line of f . Thiscompletes the proof.

Acknowledgments

A part of this work has been done using theMegaWave2 image processing environment (GeorgesKOEPFLER—Copyright (C)1993-1199 CMLA, ENSCachan, 94235 Cachan cedex, France—http://www.ceremade.dauphine.fr—All rights reserved—).The authors would like also to thank the GdR ISIS(Information, Signal, Images et viSion), http://www-isis.enst.fr of CNRS for the “gdr” test image.

Note

1. We can reformulate the minimization of∑P

i=1

∫�i

( f (x) −f �i )2dx + ∫

�( f (x) − f �)2dx + λ

∑P

i=1|�i |︸︷︷︸

|S|\|�|

as the minimiza-

tion of∑P

i=1

∫�i

( f (x)− f �i )2 dx +∫�

[( f (x)− f �)2 −λ]dx . λmay be then interpreted as a classical Gaussian threshold chosenaccording to the normal law table for a given percentage.

References

Acton, S. and Mukherjee, D. 2000. Scale space classification usingarea morphology. IEEE Trans. Image Processing, 9(4):623–635.

Alvarez, L. and Esclarin, J. 1997. Image quantization using reaction-diffusion equations. SIAM Journal of Applied Mathematics,57(1):153–175.

Alvarez, L., Gousseau, Y., and Morel, J. 1999. Scales in naturalimages and a consequence on their bounded variation. In Int. Conf.on Scale-Space Theories Comp. Vis., Kerkyra, Greece, pp. 247–258.

Alvarez, L., Guichard, F., Lions, P., and Morel, J. 1993. Axioms andfundamental equations in image processing. Arch. Rational Mech.Aanalysis, 16(9):200–257.

Alvarez, L., Lions, P., and Morel, J. 1992. Image selective smoothingand edge detection by non-linear diffusion—2. SIAM J. NumericalAnalysis, 29(3):845–866.

Alvarez, L. and Mazorra, L. 1994. Signal and image restoration usingshock filters and anisotropic diffusion. SIAM Journal of NumericalAnalysis, 31(2):590–605.

Amadieu, O., Debreuve, E., Barlaud, M., and Aubert, G. 1999. Si-multaneous inward and outward curve evolution. In Int. Conf. onImage Processing, Kobe, Japan.

Beaulieu, J. and Goldberg, M. 1989. Hierarchy in picture segmenta-tion: A stepwise optimization approach. IEEE Trans. Patt. Anal.and Mach. Int., 11(2):150–163.

Black, M. and Rangarajan, A. 1996. On the unification of line pro-cesses, outlier rejection, and robust statistics with applications inearly vision. Int. J. Computer Vision, 19(1):57–92.

Black, M. and Sapiro, G. 1999. Edges as outliers: Smoothing usinglocal image statistics. In Int. Conf. on Scale-Space Theories Comp.Vis., Kerkyra, Greece, pp. 259–270.

Black, M., Sapiro, G., Marimont, D., and Heeger, D. 1998. Robustanisotropic diffusion. IEEE Trans. Image Processing, 7(3):421–432.

Blake, A. and Zisserman, A. 1987. Visual Reconstruction. MIT Press:Cambridge, Mass.

Caselles, V., Coll, B., and Morel, J. 1999. Topographic maps andlocal contrast changes in natural images. Int J. Computer Vision,33(1):5–27.

Caselles, V., Kimmel, R., and Sapiro, G. 1997. Geodesic active con-tours. Int J. Computer Vision, 22(1):61–79.

Caselles, V., Morel, J., Sapiro, G., and Tannenbaum, A. 1998. Spe-cial issue on partial differential equations and geometry-drivendiffusion in image processing and analysis. IEEE Trans. ImageProcessing, 7(3):269–473.

Catte, F., Lions, P., Morel, J., and T., C. 1992. Image selective smooth-ing and edge detection by non-linear diffusion. SIAM J. NumericalAnalysis, 29(1):182–193.

Isophotes Selection and Reaction-Diffusion Model 93

Chan, T. and Vese, L. 1999. Active contour model without edges. InInt. Conf. on Scale-Space Theories Comp. Vis., Kerkyra, Greece,pp. 141–151.

Charbonnier, P., Blanc-Feraud, L., Aubert, G., and Barlaud, M. 1997.Deterministic edge-preserving regularization in computed imag-ing. IEEE Trans. Image Processing, 6(2):298–311.

Chesnaud, C., Refregier, P., and Boulet, V. 2000. Statistical re-gion snake-based segmentation adapted to different physicalmodels. IEEE Trans. Patt. Anal. and Mach. Int., 21(11):1145–1157.

Chu, C. and Agarwal, J. 1993. The integration of image segmentationmaps using region and edge information. IEEE Trans. Patt. Anal.and Mach. Int., 15:1242–1252.

Cohen, L. 1991. On active contour models and balloons. CVGIP:Image Understanding, 53(2):211–218.

Cohen, L. 1996. Deformable curves and surfaces in image analysis.In Int. Conf. Curves and Surfaces, Chamonix, France.

Darrell, T. and Pentland, A. 1995. Cooperative robust estimationusing layers of support. IEEE Trans. Patt. Anal. and Mach. Int.,17(5):474–487.

Deriche, R. 1987. Using canny’s criteria to derive a recursively imple-mented optimal edge detector. Int J. Computer Vision, 2(1):167–187.

Figueiredo, M., Leitao, J., and Jain, A. 2000. Unsupervised contourrepresentation and estimation using b-plines and mininim descrip-tion length criterion. IEEE Trans. Image Processing, 6(9):1075–1087.

Froment, J. 2000. Perceptible level lines and isoperimetric ratio. InInt. Conf. on Image Processing, Vancouver, Canada.

Geiger, D. and Yuille, A. 1991. A common framework for imagesegmentation. Int. J. Computer Vision, 6(3):1186–1191.

Geman, S. and Geman, D. 1984. Stochastic relaxation, gibbs distri-butions, and the bayesian restoration of images. IEEE Trans. Patt.Anal. and Mach. Int., 6(6):721–741.

Grenander, U. and Miller, M. 1994. Representations of knowledge incomplex systems. J. Royal Statistical Society, Series B, 56(4):549–603.

Ishikawa, H. and Geiger, D. 1998. Segmentation by grouping junc-tions. In Comp. Vis. and Patt. Rec., Santa Barbara, CA.

Istas, J. 1997. Statistics of Processes and Signal-Image Segmentation.University of Paris VII.

Jermyn, I. and Ishikawa, H. 1999. Globally optimal regions andboundaries. In Int. Conf. on Comp. Vis., Kerkyra, Greece, pp. 904–910.

Kass, M., Witkin, A., and Terzopoulos, D. 1987. Snakes: Activecontour models. Int. J. Computer Vision, 12(1):321–331.

Kichenesamy, S., Kumar, A., Olver, P., and Yezzi, A. 1996. Confor-mal curvature flows: From transition to active contours. Archivefor Rational Mechanics and Anal., 134:275–301.

Koepfler, G., Lopez, C., and Morel, J. 1994. A multiscale algorithmfor image segmentation by variational method. SIAM J. NumericalAnalysis, 31(1):282–299.

Kornprobst, P., Deriche, R., and Aubert, G. 1997. Image coupling,restoration and enhancement via pde’s. In Int. Conf. on ImageProcessing, Santa Barbara, USA, pp. 458–461.

Leclerc, Y. 1989. Constructing simple stable descriptions for imagepartitioning. Int. J. Computer Vision, 3:73–102.

Malladi, M. and Sethian, J. 1996. Image processing: Flows undermin/max curvature and mean curvature. Graphical Models andImage Processing, 58(2):127–141.

Matheron, G. 1975. Random Sets and Integral Geometry. John Wiley:New York.

Merriman, B., Bebce, J., and Osher, S. 1994. Motion of multiplejunctions: A level set approach. J. Comput. Physics, 112(2):334–363.

Monasse, P. and Guichard, F. 2000. Fast computation of a contrast-invariant image representation. IEEE Trans. Image Processing,9(5):860–872.

Morel, J. and Solimini, S. 1994. Variational Methods in Image Seg-mentation. Birkhauser.

Mumford, D. 1994. The bayesian rationale for energy functionals.In Geometry-Driven Diffusion in Computer Vision, Bart Romeny(Ed.). Kluwer Academic: Dordrecht, pp. 141–153.

Mumford, D. and Shah, J. 1989. Optimal approximations by piece-wise smooth functions and variational problems. Com. on Pureand applied Mathematics, 42(5):577–685.

Nitzberg, N. and Mumford, D. 1990. The 2.1d sketch. In Int. Conf.on Comp. Vis., Osaka, Japan, pp. 138–144.

Nordstrom, K. 1990. Biased anisotropic diffusion: A unified regular-ization and diffusion approach to edge detection. Image and VisionComputing, 8(4):318–327.

Olsen, O. 1996. Multi-scale segmentation of grey-scale images. InTechnical Report DIKU Nr. 96/30, ISSN 0107-8283, Departmentof Computer Science, University of Copenhagen.

Osher, S. and Sethian, J. 1988. Fronts propagating with curvaturedependent speed: Algorithms based on the hamilton-jacobi for-mulation. J. Computational Physics, 79:12–49.

Paragios, N. and Deriche, R. 2000. Coupled geodesic active regionsfor image segmentation: A level set approach. In Euro. Conf. onComp. Vis., Dublin, Ireland.

Pavlidis, T. and Liow, Y. 1990. Integrating region growing andedge detection. IEEE Trans. Patt. Anal. and Mach. Int., 12:225–233.

Perona, P. and Malik, J. 1990. Scale-space and edge detection us-ing anisotropic diffusion. IEEE Trans. Patt. Anal. and Mach. Int.,12(7):629–639.

Rougon, N. and Preteux, F. 1995. Controlled anisotropic diffusion.In Proc. SPIE Conf. on Nonlinear Image Processing VI, San Jose,CA, vol. 2424, pp. 329–340.

Rougon, N. and Preteux, F. 1998. Directional adaptive deformablemodels for segmentation. J. of Electronic Imaging, 7(1):231–256.

Rousseeuw, P. and Leroy, A. 1987. Robust Regression and OutlierDetection. Wiley Series in Probability and Mathematical Statistics,Wiley.

Salembier, P. and Serra, J. 1995. Flat zones filtering, connected oper-ators, and filters by reconstruction. IEEE Trans. Image Processing,4(8):1153–1160.

Samson, C., Blanc-Feraud, L., Aubert, G., and Zerubia, J. 1999. Alevel set model for image classification. In Int. Conf. on Scale-Space Theories Comp. Vis., Kerkyra, Greece, pp. 306–317.

Schnorr, C. 1998. A study of a convex variational diffusion approachfor image segmentation and feature extraction. J. Math. Imagingand Vision, 3(8):271–292.

Sethian, J. 1996. Level Sets Methods: Evolving Interfaces in Geom-etry, Fluid Mechanics, Computer Vision, and Material Science.Cambridges University Press: Cambridge.

Ter Haar Romeny, B. 1994. Geometry-Driven Diffusion in Com-puter Vision. Kluwer Academic Publishers: Dordrecht, TheNetherlands.

94 Kervrann, Hoebeke and Trubuil

Vincent, L. and Soille, P. 1991. Watershed in digital spaces: Anefficient algorithm based on immersion simulations. IEEE Trans.Patt. Anal. and Mach. Int., 13(6):583–598.

Wang, J. 1998. Stochastic relaxation on partitions with connectedcomponents and its application to image segmentation. IEEETrans. Patt. Anal. and Mach. Int., 20(6):619–636.

Weickert, J. 1998a. Anisotropic Diffusion in Image Processing.Teubner-Verlag: Stuttgart, Germany.

Weickert, J. 1998b. Fast segmentation methods based on partial dif-ferential equations and the watershed transformation. In Muster-erkennung 1998, P. Levi, R.-J. Ahlers, F. May, and M. Schwanz(Eds.). Springer: Berlin, pp. 93–100.

Yezzi, A., Tsai, A., and Willsky, A. 1999. A statistical approach to

snakes for bimodal and trimodal imagery. In Int. Conf. on Comp.Vis., Kerkyra, Greece, pp. 898–903.

You, Y., Xu, W., Tannenbaum, A., and Kaveh, M. 1996. Behavioralanalysis of anisotropic diffusion in image processing. IEEE Trans.Image Processing, 5:1539–1553.

Younes, L. 2000. Calibrating parameters of cost functionals. In Euro.Conf. on Comp. Vis., Dublin, Ireland.

Zhu, S. and Mumford, D. 1997. Prior learning and gibbs reaction-diffusion. IEEE Trans. Patt. Anal. and Mach. Int., 19(11):1236–1250.

Zhu, S. and Yuille, A. 1996. Region competition: Unifying snakes,region growing, and bayes/mdl for multiband image segmentation.IEEE Trans. Patt. Anal. and Mach. Int., 18(9):884–900.

(PDF) Isophotes Selection and Reaction-Diffusion Model for Object Boundaries Estimation - DOKUMEN.TIPS (2024)

References