# Bayesian image segmentation through level lines selection

Charles Kervrann

<figure><img src="image_1.png" /></figure>

To cite this version:

Charles Kervrann. Bayesian image segmentation through level lines selection. Image Analysis & Stereology, 2001, 20, pp.251-256. 10.5566/ias.v20.p163-168. hal-02676725

HAL Id: hal-02676725

https://hal.inrae.fr/hal-02676725

Submitted on 31 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
BAYESIAN IMAGE SEGMENTATION THROUGH
LEVEL LINES SELECTION

CHARLES KERVRANN

INRA - Biométrie, Domaine de Vilvert, 78352 Jouy-en-Josas, France
e-mail: ck@jouy.inra.fr
*(Accepted October 16, 2001)*

ABSTRACT

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.

Keywords: area distribution, connected components, energy minimization, image segmentation, level curves.

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 ad-
hoc methods (Paragios and Deriche, 2000). In practical
imaging, these methods still suffer from the problem
of initialization of curves (Zhu and Yuille, 1996), off-
line 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 $\mathbb{R}^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 \to \mathbb{R}^+$ 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 \in S$ is assigned a grey value $f(x)$. According to Matheron (1975), we interpret the image $f$ as a family of sets defined by $L_\tau(f) = \{x \in S : f(x) \ge \tau\}$, $\tau \in \mathbb{R}^+$. Each level set $L_\tau(f)$ is assumed to be of finite perimeter. Therefore, $f$ will belong to the bounded variation (BV) space as shown by Alvarez et al. (1999).

Let $\{\Omega_i \subset S\}$ be a set of disjoint and non-empty image domains or objects, and $\{\partial\Omega_i\}$ their boundaries. A partition of the image domain $S$ consists in finding a set $\{\Omega_i\}_{i=1}^P$ and a background $\bar{\Omega}$ defined as the complementary subset of the union of objects $\bar{\Omega} = S \setminus \cup_{i=1}^P \Omega_i$, $\Omega_i \cap_{i \neq j} \Omega_j = \emptyset$ and $\Omega_i \cap \bar{\Omega} = \emptyset$. We assume that the observed image $f$ has been produced by the model $f = f_{\text{true}} + \varepsilon$, where $\varepsilon$ is a zero-mean Gaussian white noise: $\varepsilon(x) \stackrel{iid}{\sim} \mathcal{N}(0, \sigma^2)$, $x \in S$. The true image $f_{\text{true}}(x) = \sum_{i=1}^P \bar{f}_{\Omega_i} 1_{x \in \Omega_i} + \bar{f}_{\bar{\Omega}} 1_{x \in \bar{\Omega}}$ is supposed piecewise constant, where $\bar{f}_{\Omega_i}$ and $\bar{f}_{\bar{\Omega}}$ denote respectively the unknown average values of $f$ over $\Omega_i$ and $\bar{\Omega}$, and $1_{x \in E}$ is the set indicator function of the set $E$. The variance $\sigma^2$ is assumed to be known and constant over the entire image. So, the likelihood for the data $f$ given $\{\Omega_1, \dots, \Omega_P\}$ is specified by

$$p(f | \Omega_1, \dots, \Omega_P) \propto \exp -\frac{1}{2\sigma^2} \left\{ \sum_{i=1}^{P} \int_{\Omega_i} (f(x) - \bar{f}_{\Omega_i})^2 dx + \int_{\bar{\Omega}} (f(x) - \bar{f}_{\bar{\Omega}})^2 dx \right\}. \quad (1)$$

We seek a partition of the rectangle $S$ into a finite set of objects $\Omega_i$, each of which corresponds to a part of the image where $f$ is constant. Given the objects $\{\Omega_i\}$, the unknown background is explicitly determined as the complementary subset of the union of estimated objects. Therefore, we define the following collection $\mathcal{C}_P$ of $P \ge 0$ admissible, closed and connected objects: $\mathcal{C}_P = \{\{\Omega_1, \dots, \Omega_P\} \subset S ; S \setminus \bar{\Omega} = \cup_{i=1}^P \Omega_i ; \Omega_i \cap_{1 \le i \ne j \le P} \Omega_j = \emptyset\}$. When $P=0$, there is no object in the image. Following the Bayesian approach, we use some functional of the posterior distribution $p(\Omega_1, \dots, \Omega_P | f) \propto p(f | \Omega_1, \dots, \Omega_P) \pi(\Omega_1, \dots, \Omega_P)$. The likelihood $p(f | \Omega_1, \dots, \Omega_P)$ is given by (1) and $\pi(\Omega_1, \dots, \Omega_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 $\{\Omega_1, \dots, \Omega_P\}$. We define a density that penalizes the area $|\Omega_i|$ of objects. Additionally, the variables $\{|\Omega_i|\}$ may be considered as independent random variables with density $g(|\Omega_i|)$. Hence, the prior distribution is of the form $\pi(\Omega_1, \dots, \Omega_P) = Z_p^{-1} \prod_{i=1}^P g(|\Omega_i|)^a$ where $Z_p$ is a normalization constant and $a$ a real positive value. The density $g(|\Omega_i|)$ is chosen to be a non-negative monotonically decreasing function of the object area $|\Omega_i|$. For instance, Alvarez et al. (1999) have empirically observed that the area distribution of homogeneous parts in images follows a power law $\beta|\Omega_i|^{-\gamma}$. In what follows, we shall consider this model for the density $g(|\Omega_i|)$. There are other possible choices of $g(|\Omega_i|)$: the case of $g(|\Omega_i|) \propto \exp\{-\beta|\Omega_i|^{\gamma}\}$ related 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 $p(\Omega_1, \dots, \Omega_P | f)$. 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 function $E_\lambda(f, \Omega_1, \dots, \Omega_P)$ defined as

$$E_{\lambda}(f, \Omega_{1}, \dots, \Omega_{P}) = E_{p}(\Omega_{1}, \dots, \Omega_{P}) + \lambda E_{d}(f, \Omega_{1}, \dots, \Omega_{P}), \quad (2)$$

where

$$E_p(\Omega_1, \dots, \Omega_p) = \sum_{i=1}^{P} (\gamma \log(\Omega_i) - A)$$

is the penalty functional,

$$E_d(f, \Omega_1, \dots, \Omega_P) = \sum_{i=1}^{P} \int_{\Omega_i} (f(x) - \bar{f}_{\Omega_i})^2 dx + \int_{\bar{\Omega}} (f(x) - \bar{f}_{\bar{\Omega}})^2 dx$$

the data model, $\lambda = 2a\sigma^2 > 0$ the regularization parameter and $A = \log(\beta)$. The penalty functional tends to regulate the emergence of objects $\Omega_i$ in the image. The regularization parameter $\lambda$ can be then interpreted as a scale parameter that only tunes the number of regions (Morel and Solimini, 1994; Kervrann et al., 2000). If $\lambda = 0$, each point is potentially a region and $\bar{\Omega} = \emptyset$; the global minimum coincides with zero and this segmentation is called the “trivial segmentation” (Morel and Solimini, 1994).

By using classical arguments on lower semi-continuous functionals on the BV space, we assume here the existence of minimizers of $E_\lambda(f, \Omega_1, \dots, \Omega_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)

$$ (\hat{\Omega}_1, \dots, \hat{\Omega}_{\hat{P}}) = \underset{\{\Omega_1, \dots, \Omega_p\} \in \mathcal{C}_p}{\operatorname{argmin}} \operatorname{argmin}_{0 \le P \le T} \operatorname{argmin}_{\Omega_i} E_{\lambda}(f, \Omega_1, \dots, \Omega_p), \quad (3) $$

where $\mathcal{C}_p \subseteq \mathcal{C}_T, \forall P \le T$, and $T$ is the maximum number of admissible objects registered in a bank $\mathcal{C}_T$. We recall that $\hat{\Omega} = S \setminus \cup_{i=1}^{P} \hat{\Omega}_i$ is the complementary subset of estimated objects $\{\hat{\Omega}_1, \dots, \hat{\Omega}_{\hat{P}}\}$. A direct minimization with respect to all unknown domains $\Omega_i$ and parameters $\hat{f}_{\Omega_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_{\lambda}(f, \Omega_1, \dots, \Omega_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 $\partial\hat{\Omega}_i$ of each $\hat{\Omega}_i$ is a boundary of a connected component of a level set of $f$.

**Proof of Lemma 1.** Let $\Omega_\delta$ be a variation of a set $\Omega$, i.e. the Hausdorff distance $d_\infty(\Omega_\delta, \Omega) \le \delta$. To prove Lemma 1, we assume that, for any connected perturbation of $\Omega$ such that $d_\infty(\Omega_\delta, \Omega) \le \delta$, two neighboring sets $\Omega$ and $\Omega'$ do not merge into one single set $\Omega \cup \Omega'$ and, for any connected perturbation of $\Omega$ such $d_\infty(\Omega_\delta, \Omega) \le \delta$, $\Omega$ does not split into two new sets. This corresponds to prohibited topological changes. Without loss of generality, we prove Lemma 1 for one object $\Omega$ and a background $\bar{\Omega}$, that is the closure of the complementary set of $\Omega$. For two sets $A$ and $B$ such that $B \subseteq A$, we denote $\int_{A \setminus B} f \stackrel{\text{def}}{=} \int_A f - \int_B f$. Then, we have

$$ \Delta|\Omega| = \int_{\Omega_\delta \setminus \Omega} 1 \stackrel{\text{def}}{=} |\Omega_\delta| - |\Omega| \quad \text{and} \\ \left(\int_{\Omega_\delta} f\right)^2 - \left(\int_{\Omega} f\right)^2 = 2 \int_{\Omega} f \int_{\Omega_\delta \setminus \Omega} f + \left(\int_{\Omega_\delta \setminus \Omega} f\right)^2. \quad (4) $$

We define $\Delta E_\lambda(f, \Omega) = E_\lambda(f, \Omega_\delta) - E_\lambda(f, \Omega) = T_1 + T_2 + T_3 + T_4 + T_5$, where

$$ T_1 = \int_{\Omega_\delta} f^2 - \int_{\Omega} f^2, $$

$$ T_2 = \frac{1}{|\Omega|} \int_{\Omega_\delta \setminus \Omega} \lambda \gamma, $$

$$ T_3 = -\frac{1}{|\Omega_\delta|} \left( \int_{\Omega_\delta} f \right)^2 + \frac{1}{|\Omega|} \left( \int_{\Omega} f \right)^2, $$

$$ T_4 = \int_{S \setminus \Omega_\delta} f^2 - \int_{S \setminus \Omega} f^2, $$

$$ T_5 = -\frac{1}{|S| - |\Omega_\delta|} \left( \int_{S \setminus \Omega_\delta} f \right)^2 + \frac{1}{|S| - |\Omega|} \left( \int_{S \setminus \Omega} f \right)^2. \quad (5) $$

By definition we have $T_1 + T_4 = 0$. Passing to the limit $\Delta|\Omega| \to 0$, we obtain expressions for $T_2, T_3$ and $T_5$ (higher order terms are neglected)

$$ T_2 = \frac{\lambda \gamma}{|\Omega|} \int_{\Omega_\delta \setminus \Omega} 1, $$

$$ T_3 = -\frac{2}{|\Omega|} \int_{\Omega_\delta \setminus \Omega} f \int_{\Omega} f - \frac{1}{|\Omega|} \left( \int_{\Omega_\delta \setminus \Omega} f \right)^2 + \frac{1}{|\Omega|^2} \int_{\Omega_\delta \setminus \Omega} 1 \left( \int_{\Omega} f \right)^2, $$

$$ T_5 = \frac{1}{|S| - |\Omega|} \left\{ 2 \int_{\Omega_\delta \setminus \Omega} f \int_{S \setminus \Omega} f - \left( \int_{\Omega_\delta \setminus \Omega} f \right)^2 - \frac{1}{|S| - |\Omega|} \int_{\Omega_\delta \setminus \Omega} 1 \left( \int_{S \setminus \Omega} f \right)^2 \right\}. \quad (6) $$

We define the image moments $m_0 = \int_{\Omega} 1, m_1 = \int_{\Omega} f, K_0 = \int_{S} 1, K_1 = \int_{S} f$. 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 $\int_E f(x) dE = f(x_0) \cdot |E|$ where $|E|$ denotes the area of $E$, it follows that

$$ \Delta E_{\lambda}(f, \Omega) = \left\{ \underbrace{\left[ \frac{m_1^2}{m_0^2} - \frac{(K_1 - m_1)^2}{(K_0 - m_0)^2} + \frac{\lambda \gamma}{m_0} \right]}_{M_1} + \underbrace{\left[ \frac{2(K_1 - m_1)}{K_0 - m_0} - \frac{2m_1}{m_0} \right]}_{M_0} f(x_0) - \frac{K_0}{m_0(K_0 - m_0)} f(x_0)^2 \int_{\Omega_\delta \setminus \Omega} 1 \right\} \int_{\Omega_\delta \setminus \Omega} 1. \quad (7) $$

Let $x_b$ be a fixed point of the border $\partial\Omega$. Choose $\Omega_\delta$ such that $\partial\Omega_\delta = \partial\Omega$ except on a small neighborhood of $x_b$. The energy having a minimum for $\Omega$, $f(x_b)$ needs to be solution of the following equation

$$ \frac{\Delta E_{\lambda}(f, \Omega)}{\Delta|\Omega|} = [M_0 + M_1 f(x_b)] + O(\Delta|\Omega|) = 0. \quad (8) $$
By passing to the limit $\Delta|\Omega| \to 0$, we obtain $M_0 + M_1 f(x_b) = 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 \neq 0$. The function $f$ is continuous and $\partial\Omega$ is a connected curve. Therefore $f(x_b)$ is constant when $x_b$ covers $\partial\Omega$. 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$, $\lambda$, $|\Omega_{\min}|$ be the input parameters set by the user.

**Bilevel set construction.** The first step completes a quantization of the function $f \in [f_{\min}, f_{\max}]$ in $K = \{4, \dots, 8\}$ non-equal-sized and non-overlapping intervals $[t_{l-1}, t_l)$, $l = \{1, \dots, K\}$. Given this set of intervals estimated using the maximum entropy sum method (Kapur et al., 1985), let $b_l$ be the bilevel set image with $b_l(x) = 1$ if $f(x) \in [t_{l-1}, t_l)$ and $b_l(x) = 0$ 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 $\Omega_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 $|\Omega_{\min}|$. 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 $\mathcal{C}_T$ of admissible objects $\{\Omega_1, \dots, \Omega_T\}$ with $|\Omega_i| \ge |\Omega_{\min}|$.

**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 $\mathcal{C}_T$. The background $\overline{\Omega}$ 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 $\{\int_{\Omega_i} (f(x) - \overline{f}_{\Omega_i})^2 dx\}$ are computed once and stored on a RAM memory. The energy term $\int_{\overline{\Omega}} (f(x) - \overline{f}_{\overline{\Omega}})^2 dx$ is efficiently updated for each configuration. For a fixed bank $\mathcal{C}_T = \{\Omega_1, \dots, \Omega_T\}$, one way to choose the optimal set of objects $\{\hat{\Omega}_1, \dots, \hat{\Omega}_P\}$, $\hat{P} \le T$, is to search for all possible combinations of $P$ objects and compute the corresponding energy $E_{\lambda}(f, \Omega_1, \dots, \Omega_P)$. Enumerating 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_{\lambda}(f, \Omega_1, \dots, \Omega_P)$. We start from $P=0$ and introduce one object $\Omega_j$ at a time. At the first step, we compute the $T$ energies with one single object $\Omega_j$ at once against the complementary subset $\overline{\Omega} = S \setminus \cup_{j \neq i=1}^T \Omega_i$. Let $\hat{\Omega}_1$ be the estimated object that best lowers $E_{\lambda}(f, \Omega_1, \dots, \Omega_P)$. This object is stored on a RAM memory as an object of the optimal configuration. It is removed from the initial bank $\mathcal{C}_T$. At any step of the algorithm, a new object is chosen to maximally decrease the energy $E_{\lambda}(f, \Omega_1, \dots, \Omega_P)$. Suppose that at the $P$-th step, $\hat{P}$ and $\overline{\Omega}$ are not known but we have estimated $P$ objects $\{\hat{\Omega}_1, \dots, \hat{\Omega}_P\}$ and a current background $\overline{\Omega} = S \setminus \{\hat{\Omega}_1, \dots, \hat{\Omega}_P\}$. Let $E_{\lambda}(f, \hat{\Omega}_1, \dots, \hat{\Omega}_P)$ be the current computed energy. Then at the $(P+1)$-th step, we choose the object $\Omega_j \in \mathcal{C}_T \setminus \{\hat{\Omega}_1, \dots, \hat{\Omega}_P\}$ which has the maximal difference, i.e.

$$ \hat{\Omega}_{P+1} = \arg \max_{\Omega_j \in \mathcal{C}_T \setminus \{\hat{\Omega}_1, \dots, \hat{\Omega}_P\}} \left( E_{\lambda}(f, \hat{\Omega}_1, \dots, \hat{\Omega}_P) - E_{\lambda}(f, \hat{\Omega}_1, \dots, \hat{\Omega}_P, \Omega_j) \right). \quad (9) $$
The algorithm stops at the *P*-th step when the addition of any object does not decrease $E_{\lambda}(f, \Omega_1, \dots, \Omega_p)$. This means that the optimal number of objects is $\hat{P} = P$ and the remaining objects of the bank are a part of the estimated background $\hat{\Omega} = S \setminus \{\hat{\Omega}_1, \dots, \hat{\Omega}_{\hat{P}}\}$. This algorithm selects a suboptimal configuration of objects corresponding to a local minimum. Using this algorithm, at most $(T \times (T + 1))/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 $|\Omega_i| < [0.0001, 0.001] \times |S|$ are discarded. To estimate $A$ and $\gamma$, we consider the sets of observations $\{\log(|\Omega_i|), \log(g(|\Omega_i|)), 1 \le i \le T\}$. We perform a linear regression on this set so as to find the straight line (in the log-log coordinates) $\log(g(|\Omega_i|)) = A - \gamma\log(|\Omega_i|)$ closest to the data in the least squares sense (Alvarez et al., 1999). The choice of the hyperparameter $\lambda$ determines mostly

the properties of the segmentation result. If $f$ is a function from $S$ to $[0, 255]$, a default choice for the hyperparameter is $\lambda \in [0.1, 1.] \times 255^2$. Fig. 1a shows an aerial 256 x 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$, $|\Omega_{\min}| = 0.00025 \times |S|$ and $\lambda = 0.25 \times 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 = 3.727$ and $\gamma = 1.486$ by linear regression (Fig. 1c). In that case, the residual sum of squares is 2.007 and $17 \le |\Omega_i| \le 2.886 \cdot 10^4$ pixels. It takes about 15 seconds ($25095 < (T \times (T+1))/2 = 42486$ iterations) of computing time for building $C_T$ and selecting the best configuration shown in Fig. 1e ($P = 105$ objects) using the stepwise greedy algorithm ($\lambda = 0.25 \times 255^2$). Enumerating all the configurations is infeasible since $2^T = 3.98 \cdot 10^{87}$ iterations! The non-connected background is labeled in “white” in Fig. 1e and the objects are filled with their mean gray values $\{\bar{f}_{\Omega_i}\}$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Segmentation of a (256 x 256) satellite image ($K = 8$, $|\Omega_{\min}| = 0.00025 \times |S|$, $\lambda = 0.25 \times 255^2$).</figcaption></figure>

Finally, the performance of the segmentation procedure is demonstrated for a (181 x 217) MR image shown in Fig. 2 and a microscopic medical breast (256 x 256) image shown in Fig. 3. In Fig. 2, the dark background has been previously eliminated before processing. Segmentation of both grey and white matter is achieved using the set of parameters $K = 8$, $\lambda = 0.5 \times 255^2$ and $|\Omega_{\min}| = 0.0001 \times |S|$. Fig. 3 exhibits an inflammatory carcinoma with metastases. Fig. 3b shows the segmentation results when $K = 8$, $|\Omega_{\min}| = 0.001 \times |S|$ and $\lambda = 0.5 \times 255^2$.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Segmentation of a MR image.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 3. *Segmentation of a microscopic image.*</figcaption></figure>

## 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.

## REFERENCES

Alvarez L, Gousseau Y, Morel JM (1999). Scales in natural images and a consequence on their bounded variation. In: Nielsen M, Johansen P, Olsen OF, Weickert J, eds. Scale-Space Theories in Computer Vision. Berlin, Heidelberg: Springer, 247-58.

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

Caselles V, Coll B, Morel JM (1999). Topographic maps and local contrast changes in natural images. Int J Comput Vision 33(1):5-27.

Chan T, Vese L (1999). Active contour model without edges. In: Nielsen M, Johansen P, Olsen OF, Weickert J, eds. Scale-Space Theories in Computer Vision. Berlin, Heidelberg: Springer, 141-51.

Geman S, Geman D (1984). Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE T Pattern Anal 6(6):721-41.

Grenander U, Miller MI (1994). Representations of knowledge in complex systems. J Roy Stat Soc B Met 56(4):549-603.

Kapur JN, Sahoo PK, Wong AKC (1985). A new method for gray-level picture thresholding using the entropy of the histogram. Comput Vision Graph 29:273-85.

Kass M, Witkin A, Terzopoulos D (1987). Snakes: active contour models. Int J Comput Vision 12(1):321-31.

Kervrann C, Hoebeke M, Trubuil A (2000). Level lines as global minimizers of energy functionals in image segmentation. In: Vernon D, ed. Computer Vision – ECCV 2000, 6th European Conference on Computer Vision, Dublin, Ireland, July 26 – July 1, 2000, Proceedings, Part II. Berlin, Heidelberg: Springer, 241-56.

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

Møller J, Waagepertersen RP (1998). Markov connected component fields. Adv Appl Probab 30: 1-35.

Morel JM, Solimini S (1994). Variational methods in image segmentation. Boston: Birkhauser.

Mumford D, Shah J (1989). Optimal approximations by piecewise smooth functions and variational problems. Commun Pur Appl Math 42(5):577-685.

Paragios N, Deriche R (2000). Coupled geodesic active regions for image segmentation: a level set approach. In: Vernon D, ed. Computer Vision – ECCV 2000, 6th European Conference on Computer Vision, Dublin, Ireland, July 26 – July 1, 2000, Proceedings, Part II. Berlin, Heidelberg: Springer, 224-40.

Salembier P, Serra J (1995). Flat zones filtering, connected operators, and filters by reconstruction. IEEE T Image Process 4(8):1153-60.

Vincent L, Soille P (1991). Watershed in digital spaces: an efficient algorithm based on immersion simulations. IEEE T Pattern Anal 13(6):583-98.

Zhu SC, Yuille A (1996). Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation. IEEE T Pattern Anal 18(9):884-900.