Bayesian image segmentation through level lines selection

Bayesian statistical theory is a convenient way of taking a priori information into consideration when inference is made from images. In Bayesian image segmentation, the a priori distribution should capture the knowledge about objects. Taking inspiration from (Alvarez et al. , 1999), we design a prior density that penalizes the area of homogeneous parts in images. The segmentation problem is further formulated as the estimation of the set of curves that maximizes the posterior distribution. In this paper, we explore a posterior distribution model for which its maximal mode is given by a subset of level curves, that is the boundaries of image level sets. For the completeness of the paper, we present a stepwise greedy algorithm for computing partitions with connected components.


INTRODUCTION
Image segmentation and object boundaries estimation are among the most challenging and fundamental addressed problems in image analysis (Mumford and Shah, 1989;Vincent and Soille, 1991;Morel and Solimini, 1994).Segmentation can be achieved by minimizing an energy model designed in conjunction with Bayes's theorem as shown by Mumford and Shah (1989) and Zhu and Yuille (1996).Indeed, it is straightforward to transfer a Bayesian criterion into an energy minimization criterion (Morel and Solimini, 1994;Zhu and Yuille, 1996).Thereby, the discrete (Geman and Geman, 1984;Blake and Zisserman, 1987) or continuous (Mumford and Shah, 1989;Morel and Solimini, 1994) energy functional is traditionally comprised of two terms: the first term is a fidelity term describing the interaction between the observed data and the model (data model) and the second is a regularity term (prior model).
In Bayesian image segmentation, the prior model should capture the knowledge about objects.In particular, the incorporation of prior information about the outline of objects can be applied.There has been a growing interest in this field, particularly along the guidelines of Grenander's general pattern theory using deformable templates (Grenander and Miller, 1994).In a separate way, Zhu and Yuille (1996) attempted to unify snakes (Kass et al., 1987) and region growing methods within a general energy/Bayes framework.Both approaches estimate the curves that maximally separate unknown statistics inside and outside the curves.Finally, the designed energy functionals are complex, for which it is difficult to specify global minimizers corresponding to "best" segmentations.The maximum a posteriori (MAP) estimate is generally determined by prohibitive stochastic search procedures (Grenander and Miller, 1994) or other variants of steepest ascent algorithms (Zhu and Yuille, 1996).With these tools, additional a priori knowledge may be specified to ease the segmentation task: statistics inside region boundaries are assumed to be known (Grenander and Miller, 1994;Chan and Vese, 1999) or estimated using adhoc methods (Paragios and Deriche, 2000).In practical imaging, these methods still suffer from the problem of initialization of curves (Zhu and Yuille, 1996), offline estimation of the mixture model of Gaussians approximating the probability density function of the image (Paragios and Deriche, 2000), or selection of hyperparameters weighting the contribution of energy terms (Zhu and Yuille, 1996;Chan and Vese, 1999;Paragios and Deriche, 2000).
In this paper, we address these problems and follow the Bayesian approach for recovering simply connected objects in the plane.The prior model focuses on how the area and number of objects can be varied in images.Unlike other approaches (Grenander and Miller, 1994;Zhu and Yuille, 1996), we shall see that maximizing the posterior distribution is here equivalent to selecting a subset of connected components of image bilevel sets.For the completeness of the paper, we present a stepwise greedy algorithm for computing partitions with connected components.Finally, we illustrate this approach with some experiments on satellite and medical images.

THE BAYESIAN FRAMEWORK
Let S be an open subset (rectangle) of 2 and f a grey-scale image treated as a function defined on S. Below we will work in the continuous setup, where S is a subset of a Euclidian space and f : S ¡ £¢ represents the observed data function.We use the terminology "site" or "pixel" to denote a point of the image, even in the continuous case.Each point x ¤ S is assigned a grey value f ¥ x¦ .According to Matheron (1975), we interpret the image f as a family of sets defined by assumed to be of finite perimeter.Therefore, f will belong to the bounded variation (BV ) space as shown by Alvarez et al. (1999).
Let © Ω i S be a set of disjoint and non-empty image domains or objects, and © ∂ Ω i their boundaries.
A partition of the image domain S consists in finding a set © Ω i P i 1 and a background Ω defined as the complementary subset of the union of objects Ω § S P i 1 Ω i , Ω i i j Ω j § / 0 and Ω i i Ω § / 0. We assume that the observed image f has been produced by the model f § f true ε, where ε is a zeromean Gaussian white noise: supposed piecewise constant, where f Ω i and f Ω denote respectively the unknown average values of f over Ω i and Ω, and 1 x & E is the set indicator function of the set E. The variance σ 2 is assumed to be known and constant over the entire image.So, the likelihood for the data f given © Ω 1 $ (' (' (' )$ Ω P is specified by We seek a partition of the rectangle S into a finite set of objects Ω i , each of which corresponds to a part of the image where f is constant.Given the objects © Ω i , the unknown background is explicitly determined as the complementary subset of the union of estimated objects.Therefore, we define the following collection F P of P 0 admissible, closed and connected objects: F P § G© H© Ω 1 $ (I (I (I )$ Ω P S ; S Ω § P i 1 Ω i ; Ω i 1P i jP P Ω j § / 0 QI When P § 0, there is no object in the image.Following the Bayesian approach, we use some functional of the posterior distribution The likelihood p¥ f R Ω 1 $ (' (' (' S$ Ω P ¦ is given by ( 1) and π ¥ Ω 1 $ (' (' (' )$ Ω P ¦ is the prior distribution of objects.
The posterior distribution is used in a further inferential issue concerning the objects within the Bayesian paradigm.The a priori distribution should capture the knowledge about © Ω 1 $ (' (' (' S$ Ω P .We define a density that penalizes the area R Ω i R of objects.Additionally, the variables © R Ω i R may be considered as independent random variables with a normalization constant and a a real positive value.
The density g¥ R Ω i R ¦ is chosen to be a non-negative monotically decreasing function of the object area R Ω i R .
For instance, Alvarez et al. (1999) have empirically observed that the area distribution of homogeneous parts in images follows a power law β R In what follows, we shall consider this model for the to the Markov connected components fields, has been already discussed by Kervrann et al. (2000), Alvarez et al. (1999) and Møller and Waagepetersen (1998).

BAYESIAN INFERENCE
All kinds of inference are made from Finding the maximum a posteriori (MAP) estimate is herein our choice of inference.As a consequence, the MAP estimation of objects is equivalent to the minimization of a global energy where is the penalty functional, the data model, λ § 2a σ 2 d 0 the regularization parameter and A § log¥ β ¦ .The penalty functional tends to regulate the emergence of objects Ω i in the image.The regularization parameter λ can be then interpreted as a scale parameter that only tunes the number of regions (Morel and Solimini, 1994;Kervrann et al., 2000).If λ § 0, each point is potentially a region and Ω § / 0 ; the global minimum coincides with zero and this segmentation is called the "trivial segmentation" (Morel and Solimini, 1994).
By using classical arguments on lower semicontinuous functionals on the BV space, we assume here the existence of minimizers of E λ ¥ f $ Ω 1 $ (' (' (' )$ Ω P ¦ among functions of sets finite perimeter (or of bounded variation) (Morel and Solimini, 1994;Zhu and Yuille, 1996).Our MAP estimator is defined by (when it exists) where F P s F T $ ut P v T , and T is the maximum number of admissible objects registered in a bank F T .
We recall that e Ω § S w f prove Lemma 1, we assume that, for any connected perturbation of Ω such that d ∞ ¥ Ω δ $ Ω¦ yv δ , two neighboring sets Ω and Ω do not merge into one single set Ω Ω and, for any connected perturbation of Ω such d ∞ ¥ Ω δ $ Ω¦ v δ , Ω does not split into two new sets.This corresponds to prohibited topological changes.Without loss of generality, we prove Lemma 1 for one object Ω and a background Ω, that is the closure of the complementary set of Ω.For two sets A and B such that B s A, we denote A B f def § A f V B f .Then, we have where By definition we have T 1 T 4 § 0. Passing to the limit ∆ R Ω R ¡ 0, we obtain expressions for T 2 $ T 3 and T 5 (higher order terms are neglected) We define the image moments . Using the mean value theorem for a double integral, which states that if f is continuous and a connected subset E is bounded by a simple curve, then for some point x 0 in E we have Let x b be a fixed point of the border ∂ Ω.
Choose Ω δ such that ∂ Ω δ § ∂ Ω except on a small neighborhood of x b .The energy having a minimum for Ω, f ¥ x b ¦ needs to be solution of the following equation By passing to the limit ∆ R Ω R ¡ 0, we obtain M 0 M 1 f ¥ x b ¦ o § 0. This equation has a unique solution.
The coefficients M 0 and M 1 depend on neither x b nor f ¥ x b ¦ , and M 0 p § 0. The function f is continuous and ∂ Ω is a connected curve.Therefore f ¥ x b ¦ is constant when x b covers ∂ Ω.This completes the proof.

A STEPWISE GREEDY ALGORITHM FOR IMAGE SEGMENTATION
To implement our level set image segmentation based on energy minimization, a four step method is used.The proposed algorithm is not a region growing algorithm (Morel and Solimini, 1994) since all objects are built once and for all.It differs from the watershed approach since regions that emerge from the watershed segmentation are not necessarily connected components within the image level sets (Vincent and Soille, 1991).Let K, λ , R Ω min R be the input parameters set by the user.
Bilevel set construction.The first step completes a quantization of the function f ¤ ql f min $ f max n in K § r© 4$ (' (' (' X$ 8 non-equal-sized and non-overlapping intervals l t l U 1 $ t l ¦ , l § © 1$ (' (' (' X$ K .Given this set of intervals estimated using the maximum entropy sum method (Kapur et al., 1985), let b l be the bilevel set otherwise.The connected components of bilevel sets can be characterized by their surrounding curves, that is the level lines.If we map these level lines for a given set of K levels, we get a segmentation of the image also called topographic map (Caselles et al., 1999).More generally, one can consider a segmentation achieved using only some connected components of bilevel sets, which is the philosophy of our approach.The most perceptible level lines could be determined by the detection of T-junctions of level lines (Caselles et al., 1999).Instead, we use herein a simpler criterion where perceptually significant level lines are the bilevel sets boundaries of an quantized image by using K quantizers and an entropy method.The entropy method due to Kapur et al. (1985) chooses the thresholds © t l to be the values at which the information is maximum.
Object extraction.A crude way to build pixels sets corresponding to objects is to proceed to a connected components labeling of bilevel set images © b l and to associate each label with an object Ω i .
Though this process may work in the noise-free case, in general we would also need some smoothing effect of the connected components labeling.So, we consider a size-oriented morphological operator acting on sets that consists in keeping all connected components of the output of area larger than a limit R Ω min R .This connected operator in mathematical morphology will never introduce new features or edges and boundaries of remain connected components are preserved (Salembier and Serra, 1995).The list of connected components then forms the bank F Configuration determination.The connected components are then combined during the third step to form object configurations.For instance, these configurations can be built by enumeration of all possible object combinations, i.e. 2 T configurations.Each configuration is made of a subset of connected components taken in the bank F T .The background Ω corresponds to the complementary set of objects selected for each configuration.
Energy computation and object configuration selection.Energy calculations take the image intensities of the original (not quantized) image to establish piecewise-constant approximation errors.

Energies of the form
computed once and stored on a RAM memory.

The energy term
is efficiently updated for each configuration.For a fixed bank

F
T § t© Ω 1 $ (' (' (' X$ Ω T , one way to choose the optimal set of of objects © e Ω 1 $ (' (' (' X$ e Ωf P , e P v T , is to search for all possible combinations of P objects and compute the corresponding energy all possible sets of objects in the object bank and comparing their energies is computationally too expensive if T is large.Instead of a such brute force search, we propose a stepwise greedy algorithm for minimizing E λ ¥ f $ Ω 1 $ (' (' (' S$ Ω P ¦ .We start from P § 0 and introduce one object Ω j at a time.At the first step, we compute the T energies with one single object Ω j at once against the complementary subset Ω 1 be the estimated object that best lowers E λ ¥ f $ Ω 1 $ (' (' (' X$ Ω P ¦ .This object is stored on a RAM memory as an object of the optimal configuration.It is removed from the initial bank F T .At any step of the algorithm, a new object is chosen to maximally decrease the energy E λ ¥ f $ Ω 1 $ (' (' (' )$ Ω P ¦ .Suppose that at the P-th step, e P and e Ω are not known but we have estimated P objects © xe Ω 1 $ (' (' (' )$ )e Ω P and a current background Ω § S E© e Ω 1 $ (' (' (' X$ e Ω P .Let E λ ¥ f $ )e Ω 1 $ (' (' (' S$ )e Ω P ¦ be the current computed energy.Then at the ¥ P 1¦ -th step, we choose the object Ω j ¤ F T v© xe Ω 1 $ (' (' (' )$ )e Ω P which has the maximal difference, i.e.
e Ω P ¢ 1 § arg max The algorithm stops at the P-th step when the addition of any object does not decrease E λ ¥ f $ Ω 1 $ (' (' (' X$ Ω P ¦ .
This means that the optimal number of objects is e P § P and the remaining objects of the bank are a part of the estimated background e Ω § S © e Ω 1 $ (' (' (' S$ e Ωf P .This algorithm selects a suboptimal configuration of objects corresponding to a local minimum.Using this algorithm, at most ¥ T } ~¥ T 1¦ X¦ ) 2 object configurations are examined.

EXPERIMENTAL RESULTS
Experiments were conducted on satellite and medical images.The number of bilevel sets K was set fairly low (K § 4 or K § 8) to obtain large regions and to improve robustness to noise and artifacts in the image.Regions with areas R Ω i R l 0I 0001$ 0I 001n } R S R are discarded.To estimate A and γ, we consider the We perform a linear regression on this set so as to find the straight line (in the log-log coordinates) the data in the least squares sense (Alvarez et al., 1999).The choice of the hyperparameter λ determines mostly the properties of the segmentation result.If f is a function from S to l 0$ 255n , a default choice for the hyperparameter is λ ¤ ~l 0I 1$ 1I n a} 255 2 .Fig. 1a shows an aerial 256 } 256 image (in the visual spectrum) depicting the region of Saint-Louis during the rising of the Mississippi and Missouri rivers in July 1993.We are interested in extracting the rivers and a background corresponding to textured urban areas.Fig. 1 shows the segmentation results when K § 8, R Ω min R § 0I 00025 } R S R and λ § 0I 25 } 255 2 .In this experiment, the maximum number of significant components is T § 291.The corresponding topographic map is shown in Fig. 1d.The image histogram has been quantized with K § 8 quantizers and an entropic method (Fig. 1b).We estimated the values of parameters A § 3I 727 and γ § 1I 486 by linear regression (Fig. 1c).In that case, the residual sum of squares is 2I 007 and 17 v R 3 exhibits an inflammatory carcinoma with metastases.

CONCLUSION AND PERSPECTIVES
In this paper, we have presented a Bayesian approach for extracting structures in images.Although our work is related to morphological approaches based on connected operators (Salembier and Serra, 1995), it is an independent approach since we seek minimizers of a global objective functional.In addition, we proved that our MAP estimator can be determined by selecting a subset of image level lines.A total CPU time of a few seconds on a 296 MHz workstation makes the method attractive for many time-critical applications.In terms of future directions for research, we propose to create a non-linear scale-space by successive applications of an area morphology operator to select most meaningful regions in the image.
the complementary subset of estimated objects © xe Ω 1 $ (I (I (I S$ )e Ωf P .A direct minimization with respect to all unknown domains Ω i and parameters f Ω i is a very intricate problem, even if T is low since objects are not designated.In what follows (Lemma 1), we prove that the object boundaries that minimize E λ ¥ f $ Ω 1 $ (I (I (I X$ Ω P ¦ are level lines of the function f , which makes the problem tractable.LEMMA 1 If there exist minimizers and no pathological minimum exists, then the energy minimizing set of curves is a subset of level lines of f , i.e. the border ∂ e Ω i of each e Ω i is a boundary of a connected component of a level set of f .Proof of Lemma 1.Let Ω δ be a variation of a set Ω, i.e. the Hausdorff distance d