This chapter has been written in collaboration with Jean-Michel Morel.

\section[D\'etection de bords et segmentation]{Edge detection and segmentation.}

\subsection[Vers une th\'eorie g\'en\'erale]{Towards a general theory.}

A good explanation of why the edge detection problem has become such an
important, and even obsessing problem, can be found in Lowe~\cite{Low1} :
\begin{quote} 
\em Edge detection plays an important role in many computer vision systems (and
apparently in biological vision) by identifying points of intensity
discontinuity in an image.  The locations of these intensity discontinuities
usually reflect underlying discontinuities in the geometry of surface
reflectance of a scene and thereby discount the effects of varying
illumination and imaging parameters. For this reason, edges are proved to
be one of the most reliable low-level features for bridging the gap between
image intensities and scene properties.
\end{quote}
A striking confirmation of the importance of edge detection in vision (which
even leads the perceptual system to ``see'' illusory contours ) was found by
R. von der Heydt and his lab.~\cite{HePB}. The responses of single cells in visual areas
are codified by describing the visual field : the area within which moving or
stationary bars and edges produces activity.  They found, however, that
many cells in visual area V2 responded when no actual stimulus was present
in their visual field, but when edges outside this field produced illusory
contours that cross this field.

As stated by Zucker, some of the earliest attempts to build computer
vision systems limited their input data to TV-images of block-world scenes
(e.g. Guzman~\cite{Guz}), let us quote Zucker~\cite{Zu2} :
\begin{quote}
\em The belief was widespread that a suitable local edge
operator could be found that would give a response which could be tracked
into a line drawing of the scene.  This line drawing would then serve as the
low-level segmentation upon which subsequent processes would operate. 
Now, the problem with local edge operators is that they respond not only to
the presence of edges but also to various noise configurations.  This
ambiguity led to substantial problems in the design of algorithms for
joining these edge responses into the elusive perfect line drawing.
\end{quote}
Early attempts to solve the edge detection by more global methods are Martelli~ \cite{Mart} in 1972 and Montaneri~\cite{Mont} in 1971.  \\
The ``low level'' analysis of  images presented therefore already the two
basic opposite requirements of the so called segmentation problem :  
\begin{itemize}
\item One wishes to smooth the homogeneous regions of the picture with two
scopes : noise elimination and image interpretation.  
\item On the other side, one wants to keep the accurate location of  the
boundaries of these regions.  Those boundaries are called ``step edges'' (see~\cite{Marr2}). 
\end{itemize}
This problematic has been raised  more recently by Haralick and Shapiro~\cite{HarS}: 
\begin{quote}
\em What should a good image segmentation be?  Regions of an image
segmentation should be uniform and homogeneous with respect to some
characteristic such as gray tone or texture.  Region interiors should be
simple and without many small holes.  Adjacent regions of a segmentation
should have significantly different values with respect to the
characteristic on which they are uniform.  Boundaries of each segment
should be simple, not ragged, and must be spatially accurate.

Achieving all these desired properties is difficult because strictly uniform
and homogeneous regions are typically full of small holes and have ragged
boundaries.  Insisting that adjacent regions have large differences in values
can cause regions to merge and boundaries to be lost. 
\end{quote}
These authors conclude with the following assertion.
\begin{quote}
\em As there is no theory of clustering, there is no theory of image
segmentation. Image segmentation techniques are basically ad hoc and
differ precisely in the way they emphasize one or more of the desired
properties and in the way they balance and compromise one desired property
against another.
\end{quote}
We cannot agree completely with this pessimistic conclusion.  Indeed, the
obvious fact that the segmentation criteria can vary a lot only indicates
that a general theory of segmentation should not focus primarily on the
details of these criteria.  If we trust the other assertions of these authors,
we shall be led to the following questions
\begin{description}
\item[ a)] Have the segmentation criteria a generic form, which could be therefore studied in whole generality ?
\item[ b)] Is there a computational theory to find automatically a compromise
between those competing criteria ?
\end{description}
These questions might have no good answer in the tremendously general
framework of data analysis.  However, the segmentation problem is much
more specific (two dimensional, geometrically invariant, ...).  As we shall
see, the questions a) and b) can now be adressed in such a way that :
\begin{itemize}
\item a common variational formulation can be given to most segmentation and
edge detection problems;
\item most segmentation and edge detection methods (and even the most
heuristic) can be interpreted as attempts to solve this common variational
formulation.
\end{itemize}
This point of view is certainly now well admitted in the vision community, and was implicit in Zucker's following remark~\cite{Zu2}, stated many years ago :
\begin{quote}
\em The separation of edge and region-based approaches expresses more of a
difference in initial orientation than a difference in practice.
\end{quote}
\medskip

Now, what is a variational formulation ?  It consists in summarizing all criteria concerning a set of edges $K$ in a single functional $E(K)$ with real value. In other terms, with any set of edges or ``segmentation'' $K$, is associated a
value $E(K)$ which states how ``good'' the segmentation is.  A classical
convention is to define $E$ in such a way that $E(K)$ is the better the lower it is. Is that a strong assumption to associate such a functional with a
segmentation method ?  Certainly not. We shall give in this chapter four
kinds of arguments in favour of an explicit variational formulation.
\medskip

The first one is the complexity of the heuristic methods dealing with image
segmentation.  Let us take as a characteristic example Perkins' method for
image segmentation~\cite{Per}.\\
The 9 successive steps of the proposed algorithm are  : 
\begin{description}
\item[A)] Compute ``edge regions'' by a threshold on the gradient.  An edge point is a point with high enough gradient.
\item[B)] Thin the gradient regions in order that they become onedimensional.
\item[C)] Expand the edges thus obtained into new edge regions, in order to close the small gaps between edges.
\item[D)] Determine the different regions by a connectivity algorithm.
\item[E)] Eliminate the small regions.
\item[F)] Shrink again the edge regions.
\item[G)] Add new edge pixels to close the newly appearing gaps between regions.
\item[H)] Eliminate small edge regions.
\item[I)] Calculate properties of uniform intensity regions.
\end{description}

All of the tools used in the algorithm are certainly necessary for a ``good''
segmentation and are indeed present in most algorithms.  Recent works on
image segmentation have emphasized the need for the combination of
numerous tools to obtain a sound segmentation (Gagalowicz and Monga~\cite{GagM},
Pavlidis and Liow~\cite{PavL}, Prager~\cite{Pra}).  Now, the number of parameters of such methods might be reduced as much as possible.  Only a rigorous axiomatic of the segmentation problem can provide this reduction.
\medskip

The second argument for a variational formalization of the segmentation
problem is an a priori argument  based on an obvious comparison principle. 
Any reasonable segmentation method must obey a comparison criterion. 
Indeed, given two segmentations, it must be able to state which is better. 
In other terms, a segmentation method implicitly orders all segmentations. 
Now,  every order can be quantified by a functional $E(K)$ such that $E(K) > E(K')$ if $K'$ is better than $K$, and it is no strong assumption that $E(K)$ should be a real number.  
\medskip

The third argument is factual (a posteriori).  We shall verify in this chapter
that most practical segmentation methods lead to the minimization an
explicitable (and simple) functional.  Moreover, we shall explain in section \ref{sec_translate} a simple method for ``translating'' any heuristical method into a
variational one and we shall apply it to several examples.
\medskip

The fourth argument is axiomatic and consists in deducing the variational
formulation from a classical axiomatization of image processing : the so
called ``multiscale analysis''.  Multiscale analysis introduces in image
processing a basic parameter, the ``scale'', which we shall denote by $\lambda$. There are of course several ways of defining the scale.
\medskip

Roughly speaking, the definition of scale which we adopt can be interpreted 
to be the minimal size of the details which will be kept in the analysed
image.   Now, the precise meaning of the scale depends of the kind of
multiscale analysis in consideration.  Therefore, a common notation across
several theories and variants corresponds to the fact that in the axioms to
follow, the scale has the same role.

\subsection[Principes de l'analyse multi\'echelle]{Principles of multiscale analysis.}

We define a picture as a ``grey level'' or ``brightness'' real function $u_0(x, y)$ defined at each point (``pixel'') of a domain $\Omega$, which is generally the plane, a square or a rectangle.

Multiscale analysis is concerned with the generation, from a single initial
picture $u_0$, of a sequence of simplified pictures $u_{\lambda}(x, y)$, where $u_{\lambda}(x, y)$ appears to be a rougher and sketchier version of $u_0$ as $\lambda$ increases.  In $u_{\lambda}$, the  details and features like edges are kept if their ``scale'' exceeds $\lambda$.

Denote by ${\cal S}_\lambda$ the map : $u_0 \mapsto S_{\lambda}(u_0)= u_{\lambda}$.   In the image processing litterature, $S_{\lambda}(u_0)$ is either a picture, denoted $u_{\lambda}$, or more generally a pair $(K_{\lambda}, u_{\lambda})$, where $K_{\lambda}$ is a set of ``boundaries'' or ``edges at scale $\lambda$''.

The basic (and rather obvious) properties which any multiscale analysis must verify are :
\begin{description}
\label{P_MSA}
\item[1)] Fidelity : $u_\lambda \rightarrow u_0 \mbox { as } \lambda\rightarrow 0$.
\item[2)] Causality : $S_\lambda(u_0)$ only depends on $S_{\lambda'}(u_0)$  if $\lambda > \lambda'$.
\item[3)] Euclidian invariance : If  $A$ is an isometry, then $S_\lambda(u_0\circ A) = (S_\lambda(u_0))\circ A $.
\item[4)] (In case of boundary detection) Strong causality :   $K_\lambda\subset K_{\lambda'}$  if  $\lambda > \lambda'$.
\end{description}

Since $S_\lambda(u_0)$ is obtained from $u_0$ by a ``simplification'', as axioms 1) and 2) imply, it is rather intuitive that this loss of information should be
quantified by some energy functional $E(u_\lambda)$ such that $E(u_\lambda)$ decreases as $\lambda$ increases.  Moreover, the euclidian invariance of this functional, coming from axiom 3) makes it likely to be a Lebesgue or Hausdorff integral of some bidimensional and/or monodimensional terms (in case of boundary
detection) depending on $u_\lambda$ and $K_\lambda$. We shall verify that, in practical approaches, it is always the case.  In the next sections, we examine a range of fundamental examples of multiscale analysis and discuss their variational formulation.  


\section[Filtrage multi\'echelle lin\'eaire et non lin\'eaire]{Linear and nonlinear multiscale filtering.}

\subsection[Le mod\`ele lin\'eaire]{The linear model. \newline (Hildreth, Marr, Witkin, Koenderink)}

As recalled by Zucker~\cite{Zu2}, the edge detection theory of Hildreth and Marr is
embedded in one of the most complete attempts to implement a
general-purpose low-level vision system.   Marr conjectured that the
purpose of low-level vision is to construct a Primal Sketch, i.e. a symbolic
description of the intensity changes in an image.   His system begins by
convolving local edge and bar masks over the image in several orientations
and sizes.   The responses of these masks are then parsed into a set of
assertions dealing with  the presence of various edge, line, and blob-like
structures.  Finally, a number of grouping operations attempt to eliminate
noise responses and to cluster these initial symbols into assertions about
longer lines and edges.

Let us develop with more details how edges are found, according to the
original Hildreth and Marr theory~\cite{RosT,Wi2,YuP1}.  The edges are defined as the curves where the gradient of the smoothed picture has a maximum.  (The set
of ``edges'' is therefore contained in the set of the points where the
laplacian of the smoothed signal changes sign).

This theory has been improved by Witkin~\cite{Wi2}, Koenderink~\cite{Koen}, and Canny \cite{Can}.  The low pass filtering is generally made by convolution with gaussians  of increasing variance. One understands easily the necessity of a previous low pass filtering : if the signal is noisy, the gradient will have a  lot of irrelevant maxima which must be eliminated.  Of course, strong oscillations
can be due to different causes, for instance to the presence of ``textures''.  

 In a recent paper, Torre and Poggio~\cite{TorP}  compare the properties of
several integrodifferential operators which have been proposed as edge
detectors.  One of the main results is a rigorous justification of filtering
before differentiation. It justifies the use of suitable derivatives of
Gaussian-like filters in edge detection.  Indeed, they don't create edges as
the accuracy of the filter decreases and they ensure with high probability
that zero-crossing edges form closed contours.

Koenderink~\cite{Koen} noticed that the convolution of the signal with gaussians at each scale was equivalent to the solution of the heat equation with the
signal as initial datum.  Denoting by $u_0$ this datum, the ``scale space''
analysis associated with $u_0$ consists in solving the system
\begin{equation}
\label{chal_Koen}
\left\{
\begin{array}{rcl}
    \Disp{\DP{u(x,y,t)}{t}} & = & \Delta u(x,y,t) \\
			\\
                 u(x,y,0)   & = & u_0(x,y)
\end{array}
\right.
\end{equation}
The solution of this equation for an initial datum with bounded quadratic
norm is $u(x, y, t) = G_t\star u_0$  where $G_\sigma(x, y) = \Dfrac{1}{4\pi\sigma} exp\left(-\Dfrac{x^2+y^2}{4\sigma}\right)$ is the Gauss function.
\medskip

The sequence $u_\lambda = u(., ., \lambda^2)$ defines clearly a multiscale analysis verifying the properties 1 (fidelity), 2 (causality) and 3 (invariance) listed in section~\ref{P_MSA}.  Now, many attempts have been made to obtain, or restore a strong causality property based on this multiscale analysis (Witkin, Mallat and Zhong, aso.).  Hildreth and Marr define the ``edges'' as follows.

\subsubsection*{Restoration of the strong causality.}

We shall say that $(x, y)$ is an edge point for the ``scale'' $\sqrt{t}$ if $\Delta u(x, y, t)$ changes sign and $|\nabla u(x, y, t)|$ is ``big''.  Of course, this last condition introduces some a priori defined threshold.  Unfortunately, it is well known (and it is enough to look at the ``edges'' found by this method \cite{Marr3} to see) that the edges at low scales give an inexact account of the boundaries which, according to our perception, should be considered as correct.  This is still true for the low pass filtering of Canny~\cite{Can} which is generally used as the best linear filter for white noise elimination and edge detection.  On the opposite side, if one makes a sharp low filtering, with
small variance, all the edges will keep their correct location.  Now, the
``main'' edges will be embedded in a crowd of ``spurious'' edges due to noise,
texture, and so on. The theory of Witkin~\cite{Wi2} proposes therefore to identify the main edges at a low scale, and then to ``follow them backward'' by making
the scale decreasing again.  This method could theoretically give the exact
location of all main edges.  However, its implementation is rather heavy
from the computational viewpoint and unstable, because of  the following of
edges across scales and the multiple thresholdings involved in the edge
detection at each scale. 

\subsubsection*{Variational formulation of the Hildreth-Marr theory.}

As we announced in the last section, there is an energy naturally associated
with the Hildreth-Marr multiscale analysis.  It evaluates the loss of
information on $u$ as $t$ increases : 
$$E(u) = \int_\Omega |\nabla u|^2 d\!x .$$
Indeed, the heat equation is the steepest descent method associated with
this energy.  If one looks at a discretization of the evolution equation, it can
also take a variational form, since it is $u_t-u_0 = t\Delta u$, which is the solution of a variational problem.  Roughly speaking, to solve numerically the heat equation is therefore equivalent to minimize the following functional iteratively, with
$t=\lambda^2$ small :
$$E_\lambda(u) =\lambda^2 \int_\Omega|\nabla u|^2 d\!x  \;   + \; \int_\Omega (u-u_0)^2 d\!x .$$
This formulation has its own interest : it  makes the scale parameter appear
explicitly as the multiplier of the smoothing term.  The solution $u_\lambda$ of this problem can be defined as the ``analysis of $u$ at scale $\lambda$'' as well as the solution of the heat equation $u(x, y, \lambda^2)$.  Notice, however, that the solutions $v_\lambda$ of the minimizing problem associated with $E_\lambda$ do not verify the causality property.  

\subsection[La th\'eorie de  Malik et Perona]{The Malik and Perona theory.}

The  above mentioned attempts to restore the strong causality are, as we
shall see, quite close to an important improvement of the edge detection
theory proposed by Perona and Malik~\cite{MaliP1}.

Their main idea is to introduce a part of the edge detection step into the
filtering itself, allowing an interaction between scales from the beginning
of the algorithm. \\
They propose to replace the heat equation by a nonlinear equation :
\begin{equation}
\label{chal_MalPer}
\left\{
\begin{array}{rcl}
    \Disp{\DP{u}{t}} & = &  {\rm div}(f(|\nabla u|)\nabla u) \\
			\\
                 u(O)   & = & u_0
\end{array}
\right.
\end{equation}
In this equation, $f$ is a smooth nonincreasing function with $f(0) = 1$, $f(x)\geq 0$ and $f(x)$ tending to zero at infinity.  The idea is that the smoothing process obtained by the equation is conditional: if $\nabla u(x, y)$ is big, then the diffusion will be low  and therefore the exact localization of the edges
will be almost kept.  If $\nabla u(x, y)$ is small, then the diffusion will tend to smooth still more around $(x, y)$.  Thus the choice of $f$ corresponds to a sort of thresholding which has to be compared to the thresholding of $|\nabla u|$ used in the final step of the classical theory.  Since this thresholding introduced anyway a nonlinear device, it was natural to  use it earlier in the method, in the smoothing process itself.  The experimental results obtained by Malik
and Perona are perceptually impressive and show that an edge detector
based on this theory gives edges which remain much more stable across the
scales, making therefore the backward following of edges across scales
unnecessary.  (This property has been seeked by many researchers in the two
last decades).

\subsubsection*{Variational formulation of the Malik and Perona model.}

As in the preceding section let us look for the energy associated with the Malik and Perona model. Set $\tilde{f}(a)=f(\sqrt{a})$ . Equation~\Pref{chal_MalPer} can now be written :
$$\left\{
\begin{array}{rcl}
    \Disp{\DP{u}{t}} & = &  {\rm div}(\tilde{f}(|\nabla u|^2)\nabla u) \\
			\\
                 u(O)   & = & u_0
\end{array}
\right.$$
Now set $\tilde{F}(t)=\Disp{\int_0^t}\tilde{f}(s)d\!s$ and define an energy $E$ :
$$E(u)=\int_\Omega\tilde{F}(|\nabla u|^2)d\!x.$$
Suppose that $u$ is a minimum of $E$. Calculating the first variation of $E$ we have for every test function $\varphi$ :
$$\begin{array}{cll}
\Disp{\int_\Omega}\tilde{f}(|\nabla u|^2)\nabla u\nabla\varphi\, d\!x & = & 0    \\
\\
-\Disp{\int_\Omega}\varphi\;{\rm div}\left(\tilde{f}(|\nabla u|^2)\nabla u\right)\,d\!x & = & 0    \\
\\
-\Disp{\int_\Omega}\varphi\;{\rm div}\left(f(|\nabla u|)\nabla u\right)\,d\!x & = & 0   
\end{array}$$
So equation~\Pref{chal_MalPer} is a steepest descent method for the energy defined above and we can associate the following variational multiscale analysis
$$E_\lambda(u)=\lambda^2\int_\Omega\tilde{F}(|\nabla u|^2)d\!x + \int_\Omega(u-u_0)^2d\!x.$$

\subsubsection*{Discussion.}

However, the Malik and Perona model has several serious practical and
theoretical difficulties which have been improved in a recent work~\cite{CatCLM}.  The first difficulty is a straightforward objection  that Malik and Perona
themselves acknowledged.  Assume that the signal is noisy, with white
noise for instance.  Then the noise introduces very large, in theory
unbounded, oscillations of the gradient $\nabla u$.  Thus, the conditional smoothing introduced by the model will not give good results, since all these noise
edges will be kept !  

The second difficulty arises from the equation itself: among the functions $f$
which Malik and Perona consider advisable, one finds functions of the type
$f(s) = e^{-s}$ or $f(s) = (1+s^2)^{-1}$ for which no correct theory of equation~\Pref{chal_MalPer} is available.  Indeed, in order to obtain both existence and uniqueness of the solutions, $f$ must verify that $sf(s)$ is nondecreasing.  If this condition is not verified, one can observe for some functions $f$ with $sf(s)$ nonincreasing a nondeterministic and therefore unstable process : the same picture can in theory be the initial condition of solutions divergent in time (see \cite{HolN} for simple and explicit examples, and also \cite{Dia}). In practice, that means that very close pictures could produce divergent solutions and therefore different edges.   However, it was reasonable to try to define a model where the
function $f$, which is a sort of thresholding,  could be  quickly decreasing, as
$e^{-s}$ for instance.  

In \cite{CatCLM}, a synthesis of Malik and Perona's ideas is proposed, which avoids the above mentioned difficulties : it is robust in the presence of noise and
consistent from the formal viewpoint mentioned above. 

Define the ``selective smoothing of $u_0$ at scale $\sqrt{t}$ based on estimates at the scale $\sigma$'' as the function $u(x, y, t)$ verifying
\begin{equation}
\label{chal_CaCoMo}
\left\{
\begin{array}{rcl}
    \Disp{\DP{u}{t}} & = &  {\rm div}(f(|\nabla G_\sigma\star u|)\nabla u) \\
			\\
                 u(O)   & = & u_0
\end{array}
\right.
\end{equation}
where $G_\sigma(x, y) = \Dfrac{1}{C\sigma} \exp(-\Dfrac{x^2+y^2}{4\sigma})$.  

As we have seen, $G(x, y, t) = G_t(x, y)$ is nothing but the fundamental solution of the heat equation.  Therefore the term $(\nabla G_\sigma*u)(x, y, t)$  which appears inside the divergence term of~\Pref{chal_CaCoMo} is simply the gradient of the solution at time $\sigma$ of the heat equation with $u(x, y, t)$ as initial datum.   Therefore, it appears to be an estimate of the gradient of u at point $(x, y)$, obtained by the classical Marr-Hildreth-Witkin theory recalled above.  Thus, the modification of the model of Malik and Perona is only to replace the gradient $|\nabla u|$ by its estimate $|\nabla G_\sigma\star u|$. As it is proved in~\cite{CatCLM}, this slight change of the model is enough to avoid both inconsistencies of the Malik and Perona model. The equation  diffuses at a point with more or less strength, according to the estimate of the gradient.  This estimate is achieved by the new term.  This term, like in Witkin's theory, serves to recognize the location of the main edges.  This information is then used in the equation to avoid too much diffusion at these locations.

However, this last model keeps some of the drawbacks of the previous
models: first it has no clear geometric interpretation, because the term
inside the divergence is hybrid and combines the estimate of the gradient
and the gradient.  Moreover, even if existence and uniqueness are proved, the
stability of the model as the scale parameter $\sigma$ tends to zero is generally not true, because the limit model can ``invert'' locally the heat equation. As
we shall see, these drawbacks are related to the excessive generality of the
Perona and Malik model.  

As commented in~\cite{CatCLM}, the function $G$ to be considered  can be any ``low pass filter'', or, to use the calculus terminology , any smoothing kernel.  However, in order to preserve the notion of scale in the gradient estimate, it is
convenient that this kernel depends on a scale parameter. A good and
classical example is, as mentioned above, the gaussian. It is important to
keep this particular case in mind. Indeed, a question which arises
immediately in the consideration of model~\Pref{chal_CaCoMo} is what time is best for ``stopping'' the evolution of the signal $u(x, y, t)$.  Now, one may appeal to  the Witkin model to answer this question : according to this model, time is interpreted as a ``scale factor''.  (More precisely, the solution $u(x, y,t)$ at time $t$ corresponds to a spatial scale $\sqrt t$.  Indeed, roughly speaking, $u(x, y, t)$ appears in the Witkin model as a smoothed version of $u_0$
 obtained by convolving it with a filter of spatial width $\sqrt t$).
 Thus in the model~\Pref{chal_CaCoMo} it is coherent to correlate the stopping time $t$  and the time introduced via the estimator $G_\sigma$.  One should therefore choose a stopping time $t$ of the order of $\sigma$.  Then the scale above which the signal is smoothed in regular zones of the image will be of the order of  $\sqrt t$.

On the parts of the picture where edges are present, the situation is different. Since the scope of the equation is to delay diffusion in these zones, the scale at which  edge information is lost will depend on the shape of the thresholding function $f$. Therefore, even if we wrote above that one should stop the equation at a time of the order of $\sigma$, this is rather a lower estimate and there is no inconsistency to look at what happens to the signal $u(x, y, t)$ for the time  parameter greater than $\sigma$.

In experimentations, it might be convenient to play with two parameters which are anyway  at hand in any edge detection model : the scale parameter (spatial width of the filtering) on one side, $\lambda = \sqrt\sigma$, and, on the other side, the threshold parameter for edges which is implicit in the shape of $f$ : if $f$ is near 1 on some interval containing zero and decreases briskly at the end $b$ of this interval, then $b$ is the threshold for edges : where $|\nabla u|$ is greater than $b$, the edges will remain and where it is smaller, they will disappear.  Thus it is clear that if this model is used as a preliminar step for an edge detection device, one might use the same threshold on the gradient for keeping the edges as the one implicit in $f$, i.e. $b$.  


\subsection[Diffusion anisotropique et ``mean curvature motion'']{Anisotropic diffusion and ``mean curvature motion''.}

The intuitive idea of ``edges'' is that they are generally piecewise smooth . 
Therefore, it seems natural to modify the diffusion operator so that it
diffuses more in the direction parallel to the edge and less in the
perpendicular one. In the extreme case, one could think of a diffusion
which is made only  in the direction of the edge : such a diffusion would
apparently keep exactly the location and sharpness of the edge, while
smoothing the picture on both sides on this edge.  Let us first consider this
limit case in a linear framework.  It is not difficult to see that the
diffusion equation which does not diffuse at all in the direction of the
gradient $\nabla u$ can be written as  :
\begin{equation}
\label{chal_mcm}
    \DP{u}{t}  =   \Delta u - \frac{D^2u(\nabla u,\nabla u)}{|\nabla u|^2}.
\end{equation}
 
The first term, the Laplacian, is the same as in scale space theory and the
second is an ``inhibition'' of the diffusion  in the direction of the gradient. 
Let us denote by $\xi$ the coordinate associated with the direction orthogonal
to $|\nabla u|$.  Therefore a formulation of the preceding equation with respect to this new coordinate is
$$\DP{u}{t}  = u_{\xi\xi}\;,$$
where, of course, $\xi$ depends on $\nabla u$.\\
In a quasi ``divergence form'', the equation can be also written as
$$ \DP{u}{t} = |\nabla u| {\rm div}\left(\frac{\nabla u}{|\nabla u|}\right)$$ 
and in a more litteral formulation
$$\DP{u}{t}  = \frac{u_y^2u_{xx} - 2u_xu_y u_{xy} + u_x^2u_{yy}}{u_x^2+u_y^2}.$$
This equation, which models the ``mean curvature motion'', has recently received a lot of attention because of its
geometrical interpretation : indeed, at least formally (see S. Osher and J.
Sethian \cite{OshS}, L. C. Evans and J. Spruck \cite{EvS}) the level sets of the solution move in the normal direction with a speed proportional to their mean
curvature.  

A theory of weak solutions of this equation, based upon the so-called
viscosity solution theory, has been proposed by Y. G. Chen, Y. Giga and S. Goto
\cite{ChGG},  Y. Giga, S. Goto, H. Ishii and M.H. Sato \cite{GiGIS} and H.M. Soner \cite{Son}.


\subsection[Restauration d'images par ``mean curvature motion'']{Image restoration by mean curvature motion.}

Let us now see which equation is given by the combination of the degenerate
diffusion with the above mentioned improvements of the Malik and Perona
theory : Alvarez, Lions and Morel~\cite{AlvLM} propose and study a class of nonlinear parabolic differential equations for image processing based on the above
explained ``mean curvature motion''
\begin{equation}
\label{chal_AlLiMo}
\left\{
\begin{array}{rcl}
   \Disp{\DP{u}{t}} & = &  f(|\nabla G_\sigma\star u|)\,|\nabla u|\,{\rm div}\left(\Dfrac{\nabla u}{|\nabla u|}\right) \\
			\\
         u(x,y,0)   & = & u_0(x,y)
\end{array}
\right.
\end{equation}
Roughly speaking, the interpretation of the terms of the equation are as
follows :
\begin{description}
\item[a)] The term $|\nabla u|\,{\rm div}\Disp{\left(\frac{\nabla u}{|\nabla u|}\right)} = \Delta u - \Dfrac{D^2u (\nabla u, \nabla u)}{|\nabla u|^2}$  represents a degenerate diffusion term, which diffuses $u$ in the direction orthogonal to its gradient $\nabla u$ and does not diffuse at all in the direction of $\nabla u$.  The aim of this term is to make $u$ smooth on both sides of an ``edge'' with a minimal smoothing of the edge itself.  (An edge is defined as a line along which the gradient is ``large'').
\item[b)] The term $f(|G\star\nabla u|)$ is used for the ``enhancement'' of the edges.  Indeed, it controls the speed of the diffusion : if $\nabla u$ has a small mean in a neighbourhood of a point $(x, y)$ this point is considered as the interior point of a smooth region of the image and the diffusion is therefore strong. If $\nabla u$ has a large mean value on the neighbourhood of $(x, y)$, this point is considered as an edge point and the diffusion spead is lowered, since $f(s)$ is small for large $s$.
\end{description}
Thus, the proposed model is a selective smoothing of the image, where the
``edges'' are relatively enhanced.  It preserves the strong causality as much
as possible.

As the preceding ones, this model has the minimal number of parameters
required by any image processing model : 
\begin{itemize}
\item A ``contrast'' function, represented by $f$ and which allows to decide wether a detail is sharp enough to be kept. This parameter corresponds to the
threshold on $\nabla u$ in the Hildreth-Marr theory.
\item A ``scale parameter'', given by the variance of $G$, which fixes the minimal size of the kept details in the processed picture.
\end{itemize}
Among the PDE models, this model seems to be at the frontier of the stable
models which preserve uniqueness of the solution and decrease $L^\infty$ norms of the image, being therefore associated to efficient numerical schemes~\cite{AlvLM}.  


\subsubsection*{Variational interpretation of the ``mean curvature motion'' model.}

Note that if $f(s) = \Dfrac{1}{s}$ and $\sigma=0$, then we obtain as a limit case of this new model the equation
\begin{equation}\label{eq_Cas.part}
\DP{u}{t} = {\rm div}\left(\frac{\nabla u}{|\nabla u|}\right)
\end{equation}
which is a particular case of the Malik and Perona model.  This limit case is particularly important  because it explains the geometric behaviour of the solution, which might be hidden in our general formulation.  Indeed the equation~\Pref{eq_Cas.part} corresponds to the descent method ($\partial u/\partial t = -\nabla E(u)$) associated with the following energy functional :
$$E(u) = \int_\Omega |\nabla u(x, y)| dx\!dy.$$
As we did for the Hildreth-Marr theory, we can associate with this
evolution equation a multiscale energy.  This energy has to be minimized
when the evolution equation is discretized in time : 
$$E_\lambda(u) = \lambda^2\int_\Omega |\nabla u|dx \;  + \; \int_\Omega(u-u_0)^2dx.$$ 
How can we interpret this energy ?  A particular case for $u$, which is also
the ideal case for the detection of ``step edges'' is when $u_0 =\chi_A$ is the
characteristic function of some set $A$ with smooth boundary.   Then $E(u_0)$ is
nothing but the perimeter of $A$ and the evolution equation can therefore be
interpreted as a smoothing of the boundary of this characteristic function
which tends to preserve the ``edges''.


\section[Croissance de r\'egions et de bords]{Region and edge growing methods.}

\subsection[Formalisation des m\'ethodes de croissance de r\'egions]{Formalization of region growing methods. \newline(Brice and Fennema, Pavlidis, \ldots)}

The general idea of these methods is to create a partition of the picture
into homogeneous regions by any kind of process which starts with small
regions and then grows them according to some homogeneity criterion.   One
of the first attempts to do that is the ``regional neighbour search'' of Muerle
and Allen \cite{MueA}.

Let us adopt the terms of a classification of these methods proposed by
Zucker~\cite{Zu1} and Rosenfeld and Kak~\cite{RosK}.  They distinguish between
\begin{itemize}
\item Region aggregation (which is the properly called ``region growing'') which
begins with a ``seed region'' and makes the region grow by aggregation of
pixels.  
\item Region growing : ``{\em start with a set  of highly uniform regions, e.g.,
connected components of constant grey level, or, more generally regions
extracted by some preliminary segmentation process.  Growth then takes
place by starting with one of these regions and merging neighbouring
regions with it, one at a time.}'' (\cite[p.141]{RosK}). It is important to
notice that the regions growed from seed regions may overlap, in contrast
with the definition of segmentation  as a partition given by Horowitz and
Pavlidis.
\item Partitionning :  ``{\em Many of the methods of region growing can also be used to partition the picture into regions.   More generally, picture partitions can be constructed by starting with an initial partition and allowing both merging and splitting of regions.  The initial partition may be trivial (e.g., the entire picture is the sole region; each pixel or $k\times k$ block of pixels is a region; the connected components of constant grey level are the regions ; and so on \ldots), or it can be the result of a previous segmentation process.}'' (\cite[p.145]{RosK}).
\end{itemize}

Pavlidis~\cite{Pav1} is the first author to notice that the segmentation problem
implies the introduction of a smooth approximate version of the image on
each region.  The main point of his approach is to define  a family $\cal F$ of
approximating functions, and to try to find regions $O_i$ and functions $u_i$ of the family for which $\|u_i-u_0\|$ is small enough, where the norm is an appropriate one for functions defined on $O_i$.  If for instance the family $\cal F$ is the family of constant functions, then the regions which are seeked are regions where $u_0$ is as constant as possible.  The number of regions will be decreased by the merging of regions whose constants are close.  

Horowitz and Pavlidis~\cite{HoroP2} contains probably the first general definition of region growing, which is generalized into a ``split and merge'' paradigm:
\begin{quote}
 \em Arbitrary initial partition in pyramidal data structure; split regions with
large approximation errors ; and merge adjacent regions with similar
approximation ; group similar squares into irregular regions. \\
 The notion of splitting is now algorithmically joined with the original merging
approach.  The overall orientation is functional approximation, in which
regions are described in terms of an approximating function.  The
closedness of the approximation is evaluated by an error norm.
\end{quote}
The ``split and merge principle'', developed by Pavlidis and Horowitz, is to merge adjacent regions having similar approximations and to split those regions
that have large error norms. One of the possible strategies is to split the
whole picture into four nonoverlapping squares, and split again those whose
dissimilarity measure is to high, and so on.  This induces a so called ``pyramidal structure'' on the picture, where at the bottom are the single pixels as
elementary regions, and at the top the whole picture.  At each intermediate
level, the picture is divided into squares of the same size and each level has
four times less squares than the preceding.    Thus the ``split and merge''
consists in moving up and down in this pyramid, according to the merging or
splitting operations, and thus decrease the approximation norm, until a
minimum is attained. Notice that the height in the pyramid is nothing but a
scale parameter.

The final result of a region growing method is always a partition of the
image, but the ``growing'' presentation already implies the use of a scale
parameter.  Typically, this parameter is initially very small and increases
during the region growing process.  Thus all the methods quoted above are
multiscale and verify immediately the principles 1, 3 and 4 stated on page~\pageref{P_MSA} (except the quad-tree pyramid). They can all be summarized (or \underline{improved}) in the following generic algorithm :
\begin{quote}
\em
{\bf Generic multiscale region growing algorithm.}\\
\label{G_MSRGA}Fix an initial scale $\lambda_0$ very small.  Start with an initial partition of pixels or very homogeneous regions. 
\begin{description}
\item[a)] Merge all pairs of regions whose merging ``improves'' the segmentation. Increase the scale parameter $\lambda$. Iterate.  
\item[b)] If convenient, decrease $\lambda$ and split the regions whose splitting might improve the segmentation. Go to a).
\end{description}
\end{quote}

In many methods, the scale parameter is not explicit, but appears as a
threshold parameter which evolves during the region growing process.  Now,
as in the methods which we presented in the preceding sections of this
chapter, one can distinguish between the multiscale presentations and the
energetic presentation.  For instance, Horowitz and  Pavlidis directly define
a segmentation as an optimal partition :
\begin{quote}
\em If $u_0(x,y)$ is the brightness function of a picture defined on a domain $D$, then we may attempt to divide $D$ into the minimum  number of regions where $u_0(x,y)$ satisfies certain constraints (e.g. is approximately constant).
\end{quote}
Of course, this definition depends on a hidden scale parameter, which is implicit in the constraints. 

\subsubsection*{Variational formulation of the region growing methods.}

We shall follow the remark made in section~\ref{P_MSA}, according to which there is a ``natural'' energy associated with any multiscale analysis : this energy
measures the amount of information contained in the smoothed version of
$u_0$.  In fact, a region growing method provides two results at each scale :
the boundaries of the obtained regions, whose union is denoted by $K_\lambda$, and the smoothed version $u_\lambda$ of $u$ in each region.  The amount of information remaining in the segmentation has therefore two terms : one measuring the
amount of boundaries $K_\lambda$ and their smoothness, and one term measuring the smoothness of $u_\lambda$ on each region.

The simplest energy associated with a segmentation $(u,K)$ obtained by
region growing is therefore
$$E(u, K) = \mu\,{\rm length}(K) + \int_{\Omega\setminus K} |\nabla u|^2dx.$$
If we follow the formal approach developed in the preceeding sections, we can associate with the multiscale region growing scheme a discretized  scheme
where the scale appears explicitly.  We are led to define a variational
multiscale energy by :
$$E(u, K) = \lambda^2\int_{\Omega\setminus K}|\nabla u|^2dx    +  \mu\lambda^2\int_K d\sigma   +  \int_{\Omega\setminus K}(u-u_0)^2dx.$$
This energy is known as the Mumford and Shah functional (see~\cite{MumS1,MumS2}). When $u$ is imposed to be constant on each region (the first integral term is not enough to ensure this constancy) we are led to the functional 
$$E(u, K) =  \lambda^2\int_K d\sigma    +  \int_{\Omega\setminus K}(u-u_0)^2dx.$$
This functional is the simplest possible among those associated to the region
growing methods.  It will therefore be studied with much detail in the next
chapter.  The Mumford and Shah formulation is a simple instance of a
general variational formulation of image processing problems due to
Terzopoulos~\cite{Te3}.  Following ideas of Tikhonov and Arsenin~\cite{TiA} for the solution of ``ill-posed'' problems, this author has proposed a generic variational formulation for most vision problems.
\smallskip
\newline
In the energy associated with each problem, one has two terms in competition : a ``regularity'' term and an approximation term.  The first term penalizes unsmooth estimates of the picture and the second estimates which are not close enough to the picture. (The same idea has been developped by Kass, Witkin and Terzopoulos~\cite{KasWT} for global edge detection.  See next section).


\subsection[M\'ethodes de croissance de bords]{Edge growing methods.}

 A good justification of multiscale edge growing can be found in Lowe~\cite{Low2} :
\begin{quote}
\em While edge detection is an important first step for many vision systems,
the linked lists of edge points produced by most existing edge detectors
lack the higher level of curve description needed for many visual tasks.  For
example, they do not specify the tangent direction or curvature of an edge or
the location of tangent discontinuities. (...) The multiscale analysis of
curves is complementary to any multiscale detection of the original edge
points.
\end{quote}
The most obvious method for transforming local edge information into a
global one is the "edge linking" as decribed in Rosenfeld and Kak~\cite[Vol.2]{RosK} :
\begin{quote}
\em When large pieces of edges or curves have been extracted from a picture, 
they can be further linked into global curves.  The criteria for linking two
curves into a large curve might depend on their strength (i.e. the average
strengths of the original responses), lengths, and mutual alignment or ``good
continuation''.  Of course, if global curves of specific shapes are desired,
this can also be taken into account in defining the linking criteria.
\end{quote}
Another good description of edge linking, quoted below, is given by Weiss and Boldt~\cite{WeiB} who base their edge growing method on Hildreth-Marr's edge detection technique (see also Bajcsy and  Tavakoli~\cite{BajT}).
\begin{quote}
\em Zero crossing points of the Laplacian of intensity images and the gradients
at those points are used to extract initial information about local edges. 
These edges are input to a \underline{hierarchical} linking and merging 
algorithm. Edges are linked on both geometric and intrinsic properties, e.g. if their endpoints are close and their contrasts are similar. In the merging
process, a sequence of linked edges is replaced by a longer straight line if
the approximation is sufficiently good.
\end{quote}
We underlined ``hierarchical'' in this quotation because this term means once
again that edge linking methods are multiscale.  As in the preceding section
on region growing, we now give a generic definition of edge growing and a
formula for the natural associated energy.

The initial datum for this algorithm is a set of edges $K_0$.  Each edge $e$
 of $K_0$ is defined as a curve with a length  and  a ``strength'' which is the mean value of $|\nabla u_0|$, or $\Disp{\left|\DP{u_0}{n}\right|}$.  If the edge is smooth, its curvature may be also considered.  
\begin{quote}
\em
{\bf Generic multiscale edge linking algorithm.}
\label{G_MSEGA}
\begin{description}
\item[a)] Begin with a set $K_0$ of ``very small'' edges and $\lambda$ very small.
\item[b)] Link every pair of edges for which the linking is convenient for the scale depending criterion and smooth the linked curve.  We obtain $K_\lambda$.
\item[c)] Increase $\lambda$ and go back to b).
\end{description}
\end{quote}

\subsubsection*{Variational formulation for the edge growing methods.}

We shall proceed as in the preceding section by first identifying the ``loss of
information'' in the edge growing process, and then giving the associated
multiscale functional.  In the case of region growing, the increasing
regularity of $u$ in each region was measured by some integral of its
gradient.  It is natural to take as a one dimensional equivalent an integral
depending on  length and curvature (which are the lowest order
geometrically intrinsic measures on a curve).  The loss of information in the
region growing was also given by the length of the boundaries which were
removed.  Our one dimensional equivalent is of course the cardinality of the
tips of the curves (the zero-dimensional measure of $\partial K_\lambda$).  Therefore, the ``natural'' functional associated with the edge growing methods is
$$E(K_\lambda) =\int_{K_\lambda}\left(1 + {\rm curv}(\sigma)^2\right) d\!\sigma   +   \mu\,{\rm Card}(\partial K_\lambda).$$
This is the ``smoothing'' term of the method.  By analogy with the other
multiscale methods, we need to define a simple ``fidelity'' term, which
measures how good the approximation of $K_0$ by $K_\lambda$ is.  The simplest term of the kind is $\rm{length}(K_\lambda\setminus K_0)$.  However, if the method tends to smooth the edges, this term is no good measure of the proximity of $K_\lambda$ to $K_0$ because $K_\lambda$ will move by the smoothing operation. Therefore, we can simply use the same edge criterion for $K_\lambda$ than the one which was used for $K_0$.  Since $K_0$ is found by an edge detector, which consists in detecting maxima of the gradient of $u_0$, we can simply measure the fidelity of $K_\lambda$ to edges by $-\Disp{\int_{K_\lambda}}|\nabla u(\sigma)|^2 d\!\sigma$, or, more accurately, $-\Disp{\int_{K_\lambda}|\DP{u}{n}}(\sigma)| d\!\sigma$, which measures the ``strength'' of the edge.
\medskip

Thus the most simple and generic multiscale variational energy associated
with edge growing methods is (we omit the obvious normalisation
constants):
\begin{eqnarray*}
E_\lambda(K_\lambda) &=&\lambda\left(\int_{K_\lambda}\left(1+{\rm curv}(\sigma)^2\right) d\!\sigma  + \mu\,{\rm Card}(\partial K_\lambda)\right) \\
& & \mbox{\hspace{0.5cm}}-\int_{K_\lambda}|\nabla u(\sigma)|^2 d\!\sigma -	\int_{K_\lambda}\left|\DP{u}{n}(\sigma)\right| d\!\sigma   +  {\rm length}(K_\lambda\setminus K_0).
\end{eqnarray*}
Not all terms in the preceding functional are necessary.  However, one needs
at least one term in each part of the functional (the ``causality'' part, which
has the coefficient $\lambda$, and the ``fidelity'' part).  For instance, the simplest edge linking algorithm is defined by the minimizing of 
$$E_\lambda(K_\lambda) = \lambda\,{\rm Card}(\partial K_\lambda)  + {\rm length}(K_\lambda\setminus K_0).$$
The classical ``snake method'', found by Kass, Witkin and Terzopoulos~\cite{KasWT} (and quite close to the original idea of Montaneri~\cite{Mont} and Martelli~\cite{Mart}) keeps only the following terms 
$$E_\lambda(K_\lambda)=\lambda \int_{K_\lambda}\left(1+{\rm curv}(\sigma)^2 \right) d\!\sigma  - \int_{K_\lambda}|\nabla u(\sigma)|^2 d\!\sigma.$$
We shall list other examples below.



\subsection[M\'ethodes hybrides]{Hybrid methods.}

They consist in combining both edge and region growing methods.  In most
cases, either they first run an edge growing and then a region growing, or
they run first a region growing and then an edge growing.  Of course, the
energy associated with these methods will be simply the sum of the
energies proposed in the last two subsections.  Let us consider two typical
cases.
\medskip

R.M. Haralick~\cite{Har2} has been one of the first authors to propose an hybrid
method. His approach combines edge detection and pure region growing:\\
  An edge map is computed by some linear edge operator, labeling thus each pixel as edge or non-edge.  Then neighbouring pixels, neither of which are edges, are linked to be in the same region.  The results of this technique completely rely on the quality of the edge operator.  Haralick detects edges by making at each
pixel a  cubic spline of the image on a neighbourhood of the pixel.  The
subsequent region growing may create new edges, but it keeps all edges
found in the first stage of the algorithm.
\medskip

In  a recent paper of T. Pavlidis and Y.T. Liow~\cite{PavL}\label{ref_PavL}, one can find the converse
method :     
\begin{quote}
\em We present a method that combines region growing and edge detection for
image segmentation.  We start with a split-and-merge algorithm where the
parameters have been set up so that an oversegmented image results.  Then
region boundaries are eliminated or modified on the basis of criteria that
integrate contrast with boundary smoothness, variation of the image
gradient along the boundary, and a criterion that penalizes for the presence
of  artefacts reflecting the data structure used during segmentation (quad-tree in our case).
\end{quote}
Pavlidis and Liow remark that combining two techniques for segmentation
produces much better results than what would have been obtained by either
technique alone.
\medskip

It is important to notice that the general hybrid algorithm is nothing but a
slightly more general split and merge method.  Indeed, the main new tool
which is allowed in the process is what can be called ``weak splitting''. 
Weak splitting consists in determining by some edge detector an ``edge map''. 
This edge map does not necessarily define different regions.   Now, these
edges influence the region growing by favouring the regions whose
boundaries contain edges.  Thus, the general hybrid algorithm will contain,
in arbritary order  all tools of the above analysed methods.
\begin{quote}
\em
{\bf General hybrid split and merge algorithm.}
\label{G_HSMA}
\begin{description}
\item[ a)] Find an edge map by some edge detector (or define it to be the initial grid).
\item[ b)] Grow the edges.
\item[ c)] Grow the regions.
\item[ d)] Split again the regions.
\item[ e)] Smooth the edges.
\end{description}
\end{quote}
One can find recent hybrid methods in R. Bajcsy, M. Mintz and E. Liebman~\cite{BajLL} and L. Cohen, L. Vinet, P. T. Sander and A. Gagalowicz~\cite{CohVSG}, who propose a hybrid linkage region growing method which works simultaneously on both pictures of a stereo pair, making therefore stereo matching, edge detection and region growing collaborate. In Laurent Vinet's thesis one can find more details about this technique, see~\cite[part II]{Vin}. Part I of this thesis consists of a summary and tests of region-growing devices with edge-detectors and color segmentation.

\subsection[Comment traduire un algorithme heuristique en th\'eorie variationnelle ?]{How to translate a heuristic algorithm into a variational theory ?}
\label{sec_translate}

The reader might ask why such a translation is necessary.  The above
algorithmic description yields an obvious answer : if you want to define
heuristically each one of the five steps of the preceding algorithm, you will
have to find heuristically and independently a lot of parameters (generally
thresholds fixing when to merge, when to grow, aso \ldots). As we shall see, the
advantage of an a priori  variational formulation is to reduce the number of
parameters and thresholds and to fix them automatically.  Indeed, they can
be deduced from the weights given to the different terms of the functional
to be minimized.  We shall illustrate the ``translation'' method with the most
characteristic algorithms mentionned above.  The main tool of our
translation method is the following proposition, which explains why any
split and merge method is variational.
\medskip

\begin{Proposition} Consider a merging method consisting in merging iteratively
the regions of a segmentation $(O_i)$.  Assume that the criterion for merging
$O_i$ and $O_j$  is  $a( O_i\cup O_j )-a(O_i)-a(O_j) <0$ for some function $a(.)$ depending on $(O_i)$. Then the method minimizes the energy
$$E((O_i)) = \sum_i a(O_i).$$  
\end{Proposition}

\Proof Let $(O_{i'})$ be the segmentation obtained after the merging of $O_i$  and $O_j$.  Then $E((O_{i'}))-E((O_i)) =  a(O_i\cup O_j)-a(O_i) -a(O_j) < 0$.  Thus the merging process decreases $E$.\Demf

\begin{Remark} Of particular interest is the case where several procedures are
used in the segmentation process.  For instance the merging process may
depend on two different conditions, $a(i,j)> 0$ or $b(i,j)>0$.  Then they must be both added in the energy, with a free parameter which fixes their respective
weights : $\alpha a(i,j) + b(i,j)$.  Of course a splitting  condition $c(i,j) >0$ must be put into into the energy with the opposite sign, and so on\ldots
\end{Remark}
\medskip

\subsubsection*{Example 1 : Pavlidis' first algorithm (1972).}
\addcontentsline{toc}{subsection}{\mbox{\hspace{1.4cm}}Exemple 1 : Le premier algorithme de Pavlidis}

The first algorithm of Pavlidis~\cite{Pav1} depends on a ``merging'' threshold $\lambda$:  
\begin{enumerate}
\item Initialization : take all pixels as regions.  
\item For every pair of regions $O_i$ and $O_j$ such that ${\rm Var}(g,O_i\cup O_j) < \lambda$, merge $O_i$ and $O_j$.
\end{enumerate}
The algorithm is no direct application of our preceding proposition, but
corresponds to the \underline{constrained optimization} problem
$$\min_{{\rm Var}(g,O_i)<\lambda}{\rm Card}((O_i)_i),$$
where the minimum is taken on all partitions $(O_i)_i$ of $\Omega$ into connected components $O_i$.  This formulation has the drawback of giving no control on
the smoothness of the boundaries.

\subsubsection*{Example 2 : The Beaulieu and Goldberg clustering algorithm (1989).}
\addcontentsline{toc}{subsection}{\mbox{\hspace{1.4cm}}Exemple 2 : L'algorithme de Beaulieu et Goldberg}

Let us consider now a much more recent example of region growing with a
quite close formulation.  We shall see that the improvement of the J.M.
Beaulieu and M. Goldberg~\cite{BeaG} algorithm with respect to the preceding algorithm
is to avoid the thresholding effect due to the constraints. Let us quote
Beaulieu and Goldberg's description of their algorithm.
\begin{quote}
\em A hierarchical segmentation algorithm inspired from hierarchical data
clustering and based upon stepwise optimization is now presented. (...)\\
In a merging scheme, hierarchical clustering starts with $N$ clusters corresponding to each of the $N$ data points, and sequentially reduces the number of clusters by merging.  At each iteration, the similarity measures $d(O_i, O_j)$ are calculated for all clusters pairs $(O_i,O_j)$, and the clusters of the pair that minimizes the measure are merged.  This merging is repeated sequentially until the required number of clusters is obtained.
\end{quote}
The merging criterion is the same as in numerical taxonomy : The regions which
are merged are those for which ${\rm Var}(O_i\cup O_j) - {\rm Var}(O_i) - {\rm Var}(O_j)$ is minimal, where ${\rm Var}(O)$ denotes the variance of the grey level function in the region $O$. (Or, more generally, the quadratic distance of the grey level in the region to a polynomial spline of a fixed order).  

The subjacent functional $E$ which is minimized by this process is easy to
extract. Indeed, Beaulieu and Goldberg give us the amount of energy which 
has to be ``paid'' at each merging.  They make clear that this amount, ${\rm Var}(O_i\cup O_j) - {\rm Var}(O_i)- {\rm Var}(O_j)$, must be minimal. Define, following the above proposition, the energy of a partition as the sum of the variances of the regions, that is
$$E((O_i)_i) = \sum_i  {\rm Var}(O_i).$$
Denote by $O$ the  partition before merging of $O_i$  and $O_j$ and $O'$ the partition after the merging. Then $E(O') - E(O) = {\rm Var}(O_i\cup O_j) - {\rm Var}(O_i) - {\rm Var}(O_j)$.  Thus the Beaulieu and Goldberg algorithm tends to minimize $E(O)$ under the constraint that ${\rm Card}((O_i)_i)$ is fixed. It solves the problem
$$\min_{{\rm Card}((O_i)_i)=c}\sum_i {\rm Var}(O_i),$$
which is the exact converse of the above mentionned Pavlidis method.

Thus the method of Beaulieu and Goldberg appears to be a path to obtain a
partition with small enough cardinality and variance.  It is therefore natural
to consider it as a method minimizing the energy
$$E((O_i)_i) = \sum_i  {\rm Var}(O_i)  + \lambda\,{\rm Card}((O_i)_i).$$
Instead of fixing a priori the cardinality of the segmentation, which is a
highly arbitrary choice, it seems more natural to weight the cost of the
cardinal by a scale parameter $\lambda$. We can reinterpret the Beaulieu and
Goldberg algorithm in this context : it is simply a merging method applied to
the preceding energy, with a different stopping time however.  Indeed, the
merging method stops as for any $i$ and $j$, ${\rm Var}(O_i\cup O_j) - {\rm Var}(O_i) - {\rm Var}(O_j) > \lambda$.
 
Conversely, does this variational method allows to find a solution with
fixed cardinality ? Yes, since it is enough to make the algorithm multiscale. 
One starts with $\lambda$ small.  Then the cardinality of the initial solution will be high. The scale $\lambda$ is increased slowly in the merging process until the seeked cardinality is attained.

Notice that the method has the same drawback as the above mentionned
Pavlidis method : no control on the amount or smoothness of the boundaries.


\subsubsection*{Example 3 : The Brice and Fennema ``multiregional heuristics'' (1970).}
\addcontentsline{toc}{subsection}{\mbox{\hspace{1.4cm}}Exemple 3  : L'algorithme de Brice et Fennema}

The method (see~\cite{BriF}) consists in forming first connected components of equal intensity (region growing) ; refining then with ``phagocyte'' heuristic; refining with ``weakness'' heuristic.  The phagocyte heuristic  consists in evaluating the ``weak part'' of the common boundary of two regions $O_i$ and $O_j$, defined as the part across which the jump of the gray level is small. If this weak part is large with respect to the lengths of the perimeters of both regions, then the regions are merged.  The weakness heuristic is a variant, consisting in estimating the ratio between the length of the ``weak'' part of $\partial(O_i,O_j)$ and the length of $\partial(O_i,O_j)$.  This method depends on two thresholds, one fixing what is a ``small'' jump of the gray level, the second fixing which ratio of ``weak'' boundary is enough to merge two regions.

Let us give a translation of the ``weakness heuristic''.  The merging
condition of this heuristic is
$$\begin{array}{ll}
\mbox{If}  & \beta\,{\rm length}(\partial(O_i,O_j)) - {\rm length} \left\{ (x,y)\in\partial(O_i,O_j) \Disp{\left/\, \left|\DP{g}{n}\right|\right.} < \delta \right\} < 0,    \\
\\
\mbox{then} & \mbox{merge $O_i$ and $O_j$}. 
\end{array}$$
Let us follow again the above proposition and define an energy functional by
simply summing  the opposites of these merging decision coefficients over
all pairs of neighbouring regions $O_i$ and $O_j$. Denoting by $K$ the union of all boundaries of the $O_i$, we obtain 
\begin{eqnarray*}
E((O_i)_i)&=&-\sum_{i,j} \beta\,{\rm length}(\partial(O_i,O_j)) + {\rm length}\left\{(x,y)\in\partial(O_i,O_j)\Disp{\left/\, \left|\DP{g}{n}\right|\right.} < \delta\right\}         \\
\\
          &=& -2\beta\,{\rm length}(K) + 2\,{\rm length}\left\{(x,y)\in K \Disp{\left/\, \left|\DP{g}{n}\right|\right.} < \delta \right\} \\
\\
          &=& 2(1-\beta)\,{\rm length}(K) - 2\,{\rm length}\left\{(x,y)\in K \Disp{\left/\, \left|\DP{g}{n}\right|\right.} > \delta\right\}.
\end{eqnarray*}
This achieves the translation of the Brice and Fennema method as a
variational method.  Indeed, according to the merging condition, two regions
will be merged if and only if the energy $E$ decreases by the merging.

Let us now look whether this energy can be rationalized. We mean by
rationalization every simplification of the algorithm which will avoid
thresholding effects.  The second term of the merging criterion is
essentially defined to measure the contrast of the boundary.  This can be
done without using any threshold $\delta$ by simply setting
$$E(K) = \lambda\,{\rm length}(K) - \int_K \left|\DP{g}{n}\right|d\!\sigma.$$



\section[Exemples et variantes]{Examples and variants.}

\subsection[Le mod\`ele 2.1-D de Nitzberg et Mumford]{The 2.1-D sketch of Nitzberg and Mumford.}

The aim of Nitzberg and Mumford~\cite{NiM} is to propose a variant of the Mumford and Shah segmentation model allowing regions to overlap.  If a region occludes
another, this can be interpreted as a ``depth'' information : one of the objects
is in front and the other in back.  One good indication for occlusion is given
by the so called ``T-junctions'' : That means that when an object partially
occludes another, its edges will cut the edges of the object in back,
creating thus a junction of two edges with the shape of a ``T''.  The Gestalt
psychology work of Kanisza~\cite{Kani} has shown how humans make such reconstructions unconsciously and automatically using T-junctions and other singularities in the image edges. The T-junctions seem to be one of the most relevant informations for ordering objects in space.  Indeed, K. Nakayama~\cite{NakSS} has explored human perception of data in which various clues for depth (T-junctions, stereo, motion, etc.) conflict.  His results show in many cases that T-junctions override other clues and are always powerful organizing forces in an image.
\medskip

The data structure proposed by Nitzberg and Mumford consists in regions $R_i$
in the domain $D$ of the image corresponding to the distinct objects $O_i$ in
front of the lens where the $R_i$ are partially ordered by occlusion between
the objects $O_i$ and where the parts of the objects $O_i$ behind other $O_j$s are reconstructed, so that the $R_i$s overlap.\\
To compute this data structure, they seek a decomposition of $D$ into
overlapping ordered regions $\{R_i\}$ such that 
\begin{description}
\item[a)] The intensity of the image on the visible part of each $R_i$ is as constant as possible. 
\item[b)] The edges $\partial R_i$ are as short and straight as possible. In particular, invisible edges are reconstructed by the family of curves, called {\em elastica}, that minimizes a sum of length and the square of curvature.
\end{description}
The paper of Nitzberg and Mumford describes an algorithm that constructs a \mbox{2.1-D} sketch of an image from prescribed edges, and gives results for
several example images, including some well-known ``optical illusions''.

The segmentation energy is therefore the following variant of the first
energy proposed by Mumford and Shah~\cite{MumS2} :
$$E_{2.1}(\{R_i\},<) = \sum_{i=1}^n  \left[  \mu^2\int_{R'_i} (g-m_i)^2d\!x + \epsilon\int_{R_i}\; d\!x + \int_{\partial R_i\setminus\partial D}  \phi(\kappa)d\!\sigma \right],$$
where :
\begin{quote}
\begin{description}
\item[] $R_i$ is a family of regions of $D$ such that $\Disp{\bigcup_i}R_i=D$;
\item[] `` $<$ '' is a partial ordering on $R_i$, that represents occlusion :\\
 $R_i<R_j$ if $R_i$ is occluded by $R_j$;
\item[] $R'_i$ is the non occluded or ``visible'' part of $R_i$, defined by\\
 $R'_i= R_i\cup\{R_j \mbox{ such that } R_i<R_j\}$;
\item[] $m_i$ is the mean value of $g$ on $R'_i$;
\item[] $\kappa$ is the curvature on $\partial R_i$, and $\phi$ is defined by
$$\phi(\kappa)=
\left\{
\begin{array}{lcl}
\nu + \alpha \kappa^2 & \mbox{ for } & |\kappa|<\Dfrac{\beta}{\alpha}  \\
\\
\nu + \beta \kappa    & \mbox{ for } & |\kappa|\geq\Dfrac{\beta}{\alpha}
\end{array}
\right.
$$
with  real parameters $\alpha$ and $\beta$.
\end{description}
\end{quote}
The first term of the functional $E_{2.1}$ minimizes the variance of the visible parts.  The second term corresponds to the assumption that the occluded parts have a minimal surface.  However, the coefficient $\epsilon$ is assumed to be small and this criterion is therefore used only in ambiguous cases.  Finally, the third term asks that the boundaries of regions be short and not too curvy. Let us quote the authors to explain the definition of $\phi(.)$ :
\begin{quote}
\em The purpose of $\phi(.)$ changing from a quadratic to a linear function at $\Dfrac{\beta}{\alpha}$ is to allow us to extend the definition of the second term to a contour which is smooth except at a finite number of infinite-curvature points, or ``corners''.
\end{quote}



\subsection[Le mod\`ele de Sha'ashua et Ullman]{The ``structural saliency'' of Sha'ashua and Ullman.}

Sha'ashua and Ullman's~\cite{Sha'U} goal is to construct a saliency map which is a representation of the image emphasizing salient locations.  They seek to
associate, therefore, a measure of saliency denoted by the function $\Phi(x)$ to
each location $x$ in the image.  A property that seems to play a role in
structural saliency is the combination of length and smoothness measured
at a particular scale.  That is, a measure of saliency would favour long
smooth curves, where the smoothness of a curve is related to its curvature
or its curvature variation.\\
They therefore face the following problems : 
\begin{enumerate}
\item Defining an appropriate measure $\Phi(x)$ that, when applied to a point along
a given curve, will increase when the curve increases in length and
smoothness.  
\item A selection problem.  The measure $\Phi(x)$ depends on the curve passing
through $x$.  Since the curves we are considering are either continuous or
separated by any number of gaps, there will usually be many possible curves
to consider. Their approach to this problem will be to select the curve that
maximizes $\Phi(x)$ over all curves passing through $x$.
\end{enumerate}

Let us describe in a continuous framework how $\Phi$ is defined.   We assume
that an ``edge map'' $\cal E$ is given, that is, a set with bounded Hausdorff measure
which can be (in a computational context) assumed to be a finite union of
curves.  Fix a reference length $\lambda$ for the curves along which the ``saliency'' of $x$ will be computed.  Given a curve $c(s)$, parametrized with its arc length $s$ in $[0,\lambda]$ starting at $x$ (i.e. $c(0) = x$), we define the saliency measure of $c$ as
$${\cal F}(c) = \int_{c([0,\lambda])\setminus {\cal E}}  \exp\left\{ -\int_{c([0,\sigma])} {\rm curv}(s)^2 d\!s  -  \alpha \MH(c([0,\sigma])\setminus {\cal E})\right\} d\!\sigma ,$$
where $\alpha$ is a positive attenuation coefficient, $\MH$ denotes the Hausdorff length, ${\rm curv}(s)$ is the curvature of the curve $c(s)$ at point $s$.\\
The saliency measure of $x$ is the upper bound of ${\cal F}(c)$, taken over all smooth curves $c$ starting at $x$ and with length $\lambda$. In other terms,
$$\Phi_\lambda (x) = \sup_{\{c\,/\,\MH(c)=\lambda\}} {\cal F}(c).$$
The attenuation coefficient $\alpha$ penalizes ``gaps'' in the edge map along the curve $c$.  If for instance $c$ is entirely contained in $\cal E$, then there is no
attenuation of the integrand.  It is interesting to notice that, in a discrete
framework, the computation of the structural saliency at each point can be
made quick and parallel by using classical algorithms of minimal path
search~\cite{Mars}.  In contrast with many problems in computer vision, there is
therefore no combinatorial blow up in the computation of the saliency map.


\subsection[La m\'ethode de ``longueur de description minimale'' d'Yvon Leclerc]{The ``minimum description length'' method of Yvon Leclerc.}

Yvon Leclerc~\cite{Lec} proposes to interpret the segmentation problem as a
partitioning problem, where the criterion for prefering a partition to
another should be its ``description length''.  Of course, the measure of the
description length must be done according to an a priori fixed language. 
The Minimum Description Length method (MDL) thus defined is related to the
classical stochastic approach (MAP : Maximum a posteriori probability), but
it requires only intuitive information theory concepts. We shall not
describe the details of the formalisation, but give directly the functional to
which Yvan Leclerc is lead :
\medskip
\newline
a) Piecewise constant model for a picture :
$$E(u, K) =  \frac{1}{\sigma^2} \sum_{i\in I}(u(i) - g(i))^2  +  \beta\sum_{i\in I} \sum_{j\in N(i)} (1- \delta(u_i-u_j)),$$
where $g$ is the picture defined on the set of pixels $I$, $\sigma$ is the variance of the white noise added to the picture, $\delta$ is the Kronecker symbol, $N(i)$ is the set of the neighbouring pixels of $i$, $\beta$ is inversely proportional to the average length of the boundary of a region.
\medskip
\newline
b) Piecewise polynomial model for a picture : In this case, the first term
remains the same, but is augmented with the number $n(i)$ of nonzero Taylor
coefficients at $i$, multiplied by the number $\alpha$ of bits used for coding such a coefficient.  The ``length'' term is  replaced by a term $m(i, j)$ proportional to the number of jumps of the Taylor coefficients of $u$ between $i$ and its neighbours.  Thus the functional becomes 
$$E(u, K) = \frac{1}{\sigma^2} \sum_{i\in I}(u(i) - g(i))^2  + \alpha\sum_{i\in I}n(i) + \beta\sum_{i\in I} \sum_{j\in N(i)} m(i,j).$$
However, the coefficient $\alpha$ has to be normalized in order to be inversely
proportional to the average surface of a region, and $\beta$ is normalized in order to be inversely proportional to the average length of a region boundary.
\medskip
\newline
Further improvements :\\
Leclerc makes the interesting remark that the noise
variance is neither necessarily known, nor necessarily constant on each
region.  Therefore, the first version of the functional (the piecewise
constant case) should be replaced by
$$\sum_r \sum_{i\in O_r}\left[ \frac{1}{\sigma_r^2}(u_r - g(i))^2 \right] + \sum_r \log \sigma_r + \sum_r \MH(\partial O_r),$$
where the $O_r$ denote the regions where $u$ is constant, $\partial O_r$ their boundaries, $\MH$ the length.  The second term measures the bit length of the description of the $\sigma_r$'s.
\medskip
\newline
Numerical simulations for the Leclerc model : \\
It is a homotopy method, quite close to Blake and Zissermann's ``Graduate
Non Convexity Algorithm''~\cite{BlakZ}.  It consists in introducing a smoothing parameter $s$ and replacing $\delta(u_i-u_j)$ by $\exp\left(-\Dfrac{1}{s^2\sigma_r^2}(u_i-u_j)^2\right)$.  Then a local minimum of the new functional $E_s(u, K)$ is obtained by a simple descent method for a first large value of $s$. The obtained solution is used as a starting point for the descent method applied with a smaller value of $s$. This process is iterated while $s$ tends to zero. Then the point between pixels $i$ and $j$ is a boundary point provided $\delta(u_i-u_j)$ is large enough.
\medskip
\newline
Comments : \\
It is clear that the Leclerc model is a variant of the Mumford
and Shah general model, where piecewise analytic functions are replaced by
piecewise polynomials.  However, the Mumford and Shah model allows ``open''
boundaries to arise, that is, boundaries which do not separate two different
regions. Now, the piecewise constant -or piecewise polynomial with an
upper bound on the degree- model is a partitionning model.  Therefore, it
should be expected that an adapted algorithm computes a partition.  The
Leclerc simulations (like those of Blake and Zissermann and Richardson~\cite{Ric}) do not yield such a partition, but only broken boundaries which give not enough information for properly partitionning the image.  Another drawback of the
Leclerc model is the presence of ``hidden parameters'' inside the weights of
the different terms : indeed, the parameters $\alpha$, $\beta$, $\gamma$ contain a priori estimates on the average surface of the regions and the average length of their boundaries.  So the method should require that these parameter
estimates are confirmed a posteriori.  Or, if it is not the case, they should
also be considered as unknowns of the problem. 


\subsection[Segmentation d'images de profondeur (``range-images'')]{Segmentation of range images.} 

A. Leonardis, A. Gupta and R. Bajcsy~\cite{LeoGB} present an algorithm which aim it is to divide the image into regions well described by a polynomial spline.  The ideas of this paper  are quite similar to P.J. Besl and R. Jain~\cite{BeslJ}, but go further in the analysis of the model selection for the overlapping regions.

The first step of the algorithm is a region growing.  From seed regions
located on a grid-like pattern, regions are growed.  The growing of a region
stops when the quadratic distance between the polynomial model and the
image becomes greater than some a priori fixed threshold.  The growing
procedure outputs many different regions of which many are partially or
totally overlapped.  The criteria for selecting the models are : first that
their number is as small as possible, second that the size of each model is
as large as possible, and third that the error measure between the original 
data and the recovered models is small.  These aims are attained by
minimizing boolean objective functions of the models.
\medskip 

Denote by $m_i$ the boolean variable which is 1 if the model $i$ is used and 0 if not. For $M$ models the model selection functional to be minimized is
$$E(m_1,\ldots,m_i,\ldots, m_M) = \sum_{i,j\in\{1,\ldots,M\}} c_{ij}\,m_im_j.$$
The efficiency $c_{ii}$ of the $i$-th model is measured by 
$$c_{ii} = K_1n_i - K_2\xi_i-K_3N_i,$$ 
where the $K_j$ are suitable positive constants, $n_i$ is the number of pixels explained by the $i$-th model, $\xi_i$ is the quadratic description error of the model $i$ in the $i$-th region, and $N_i$ is the number of nonzero coefficients in the polynomial description of the model $i$. This last term is also present in the ``minimal length  description'' proposed by Y. Leclerc.  

Then the interaction coefficient $c_{ij}$ which gives an account of  what is
gained when both models $i$ and $j$ are simultaneously selected, is given by
$$c_{ij} = -K_1 n_{ij}  + K_2 \xi_{ij},$$
where $n_{ij}$ is the area of the subregion common to both models, and $\xi_{ij}$ is an error correction term  on the union of both regions $i$ and $j$.  This last term is computed by assuming that on the intersection of both regions, the best fitting model is chosen.

The highly nonconvex functional $E$, analogous to those defined by Blake and
Zissermann, Mumford and Shah, Leclerc, \ldots, is minimized by a similar device:
a continuation method where a scale parameter is introduced and a serie of
more and more nonconvex problems is solved by a descent method which
takes at step $n$ the solution at step $n-1$ as starting point.  After this model selection has been achieved, the regions, which may overlap, need
refinement of their spatial extent.  This is done by some local
convexification of the regions, which tends to avoid cusps on the region
boundaries and to control the curvature of the boundaries of the regions.


\subsection[Algorithmes de minimisation de l'\'energie de Mumford et Shah]{Algorithms for minimizing the Mumford and Shah energy.}

Blake and Zisserman have devoted a whole book~\cite{BlakZ} to the problem of minimizing the Mumford and Shah energy ( the weak membrane model ) . We shall limit ourselves to give here an accurate idea of the ``Graduated Non Convexity'' Algorithm.  There is no restriction in working in the case of a one-dimensional datum instead of a picture ; the principles of the algorithm being quite the same.  Let us consider the Mumford and Shah
energy for analysing a one dimensional datum $d(x)$ defined on a real interval
$I$ :
$$E(u, K) = \int_I\, [(u(x)-d(x))^2 + u'(x)^2 ]\,d\!x  + {\rm Card}(K).$$
(For sake of simplicity we omit the scale parameters).  In order to
discretize this energy, one defines a grid of points on $I$ : $i$ = 1, 2, \ldots, $n$ . \\
Define
$$
g^0(t)= 
\left\{
\begin{array}{ccl}
 t^2 & &\mbox{ if } t\leq 1 \\
 1   & &\mbox{ if } t\geq 1
\end{array}
\right.
$$
Then a discrete version of $E$ is
$$E(u) = \sum_{i=1}^n(u_i-d_i)^2 + g^0(u_i-u_{i-1}),$$
and $K$ is defined as the set of the $i$ for which $|u_i-u_{i-1}| >1$. This energy is non convex and can have a lot of local minimizers.  The idea of Blake and
Zisserman is to convexify $E$, in order to get a global minimizer.  Then this
minimizer is used as a starting point for a less convex version of $E$.  Thus a
series of more and more nonconvex versions of $E$ is built which are successively minimized.  Of course, nothing ensures that this method will provide a global minimum of $E$.  It is only a sound heuristic.
\smallskip
\newline
Define
$$
g^p(t)=
\left\{
\begin{array}{ccccccc}
      t^2                &\mbox{ , if } &            &    &|t|&<&\Dfrac{1}{r}\\
 \\
1 - \Dfrac{(|t|-r)^2}{4p}&\mbox{ , if } &\Dfrac{1}{r}&\leq&|t|&<&  r        \\
 \\
       1                 &\mbox{ , if } &     r      &\leq&|t|
\end{array}
\right.
$$
where $r = \sqrt{4p+1}$.
\medskip
\newline
Set the smoothed versions of $E$ to be
$$E^p(u) = \sum_{i=1}^n(u_i-d_i)^2 + g^p(u_i-u_{i-1}),$$
Then it is shown that if $p\geq 1$, $E^p$ is convex.  Thus a global minimizer can be easily found for $E^p$ if $p =1$.  Note that $E = E^0$. So the progressive
nonconvexity method consits in letting $p$ decrease from 1 to 0.
\medskip
\newline

In the first part of his thesis, N.K. Nordstr\"{o}m~\cite{Nor} studies the problem of defining procedures for minimizing the Mumford and Shah functional, with a term added, which is the number of curves of the segmentation.  This addition, also studied in Richardson's thesis (see below), makes  very easy the existence results for the minimization problem.  Indeed, take a minimizing sequence $(K_n, u_n)$ of segmentations.  One can assume that $K_n$ is made of a fixed number of curves.  Since the total length of these curves is therefore bounded from above,  each one of them can be assumed to converge by a trivial application of the Ascoli-Arzela theorem.  The length being lower semicontinuous, it is clear that the limit segmentation $K$ has a length which does not exceed the $\inf$ limit of the lengths of the minimizing sequence.  The fact that the $u$ associated with $K$ has less energy than the energies of the $u_n$ is an easy consequence of the lower semicontinuity of the $H^1$ Sobolev norm on any open set of $\R^2$. However, Nordstr\"{o}m does not mention or use this fact, but tends to define much more restrictive classes of segmentations. Roughly speaking, he imposes that each region (or connected component of $\Omega\setminus K$ ) of the segmentation is an open set with Lipschitz controlled boundary. The assumption is strong enough to ensure that each region will remain without cusps and without endfree edges.\\ 
The proposed algorithm however is general and admits a priori any kind of
configuration for edges.  Its steps are :
\begin{description}
\item[a)] Initialize the segmentation using some edge detector.\\
\item[b)] Evolve the segmentation according to the steepest descent direction
suggested by the Euler equations.  Since the proposed Euler equations do not
consider the behavior of endfree edges, these are treated heuristically.
\end{description}
Nordstr\"{o}m argues that edges with free endpoints are rare, and that
their value for image segmentation purposes, one of the prime motivations
for edge detection, is limited, as they do not give rise to any new
components of the continuity set.
\medskip
\newline

Tom Richardson~\cite{Ric} studies several properties of the Mumford and Shah
functional and presents experiments about its minimization in his thesis. A first result is on a stronger version of the functional, where the number of connected components of the boundaries is assumed to be bounded : a simple proof of
the lower semicontinuity of the functional is then given.  This is a slightly
more general result than the lower semicontinuity result under the
assumption that $K$ is made of a finite number of curves.  Another interesting
development is made on the recovery properties of the minimization method
: assume that the image completely fits  with the a priori model of the
functional.   In other terms, assume that the picture is piecewise constant,
with a finite amount of length for the boundaries.   Then Richardson proves
that the minimisation process finds all these boundaries as the scale
parameter tends to zero. He uses finite elements to discretize the Mumford and Shah functional with the $\Gamma$-convergence approach of Ambrosio and Tortorelli \cite{AmbT} .


\subsection[Variantes de la diffusion anisotropique]{Variants of the anisotropic diffusion}

Osher and Rudin's theory \cite{OshR} tries to get as close as possible to the inverse heat equation by defining some conservative scheme like 
$$\DP{u}{t} = -|u_x|F(u_{xx}) \mbox{ \hspace{0.5cm} with initial value $u_0$},$$
where $F$ is a function such that $sF(s)\geq 0$.  The big advantage of this method is to have a scheme which lets the image develop true edges, that is shock lines along which $u(x, y, t)$ becomes discontinuous in $(x, y)$.  
\medskip

Nordstr\"{o}m~\cite{Nor} proposes a new presentation of Perona and Malik's theory (\cite{MaliP1}) which relates it to the well known variational global edge detection methods of Kass, Witkin, Terzopoulos~\cite{KasWT},  Mumford and Shah~\cite{MumS1}.  He introduces a new term in the equation~\Pref{chal_MalPer} which forces $u(x, y, t)$ to remain close to $u_0$.  Because of the forcing term $u-u_0$, the new equation 
$$\DP{u}{t}- {\rm div}(f(\nabla u)\nabla u) = u_0-u,$$
has a priori the advantage of having a nontrivial steady state, eliminating
therefore the problem of choosing a stopping time.
\medskip

The approach of Nitzberg and Shiota~\cite{NiS}  is related to adaptive filtering
methods which where introduced for TV images denoising by Graham~\cite{Gra}.  
The idea is to blurr selectively and anisotropically the signal with ``oblong''
gaussians.  The gaussian used for blurring at a point $(x, y)$  depends on the
intensity and direction of the gradient in the neighbourhood.  Roughly speaking, the blurring will be faster in the direction orthogonal to the gradient.  Therefore, the signal will be smoothed on both sides of an edge, but the edge is conserved.  Since a corner is the crossing point of two edges, there will be two directions of ``nondiffusion'' instead of one, and therefore corners are well conserved by this method.  Moreover, Nitzberg and Shiota prove by a scaling argument that as the size of the Gaussians tends to zero, their diffusion method tends to some ill posed partial differential equation analogous to Malik and Perona's one.


\subsection[M\'ethodes d'histogrammes]{Histogram methods.}

Histogram methods have appeared early in image processing (see Prewitt and 
Mendelsohn in 1966 \cite{PreM}).  A brief account of histogram methods can be found in
Zucker~\cite{Zu2} :
\begin{quote}
\em A histogram, in its simplest form,  is a plot of the number of
picture points at each possible intensity level.  Assume that the different
objects in a picture are homogeneous and have different grey levels.  Then
the grey level histogram will have several modes or peaks.  It is a simple
matter to separate these peaks with  valley-centered thresholds, to label
all the points in a given mode identically, and to map these labelled points
back into their pictorial locations.   The final appropriate pictorial regions
can be obtained by joining neighbouring points which have identical mode
labels.
\end{quote}

Weszka~\cite{Wesz} gives a complete study of the particular case where only two
different modes are to be separated (for instance a dark object on bright
background).  Then the method is equivalent to the choice of a threshold
separating the histogram in two peaks. 
\medskip

A refinement of this methods involving not only the grey level, but also the
gradient has been proposed by Panda and Rosenfeld~\cite{PanR}.  They describe an approach for segmenting a white blob against a dark background:
\begin{quote}
\em Consider the histogram of gray levels for all pixels which have small gradients. If a pixel has a small gradient, then it is not likely for it to be an edge. If it not an edge, then it is either a dark background pixel or  a bright blob pixel.  Hence, the histogram of all pixels having small gradients will be bimodal and for pixels with small gradients, the valley between the two modes of
the histogram is an appropriate threshold point.  Next consider the histogram of grey levels for all pixels which have high gradients. If a pixel has a high gradient, then it is likely for it to be an edge. If it is an edge separating a bright blob against a dark background, and if the separation boundary is not sharp but somewhat diffused, then the histogram will be unimodal,  the mean being a good threshold separating the dark background pixels from the bright blob pixels.
\end{quote}
Thus Panda and Rosenfeld suggest determining two thresholds; one for low gradient pixels and one for high gradient pixels. By this means, they perform the clustering in the two dimensional measurement space consisting of gray tone intensity and gradient.
\medskip

Another refinement of the method, quite analogous to Pavlidis' split and
merge has been experimented by Ohlander, Price and Reddy~\cite{OhPR} : it
consists in making recursive histogram segmentation from ``coarse to fine''. 
At each step of the proposed algorithm,  a feature histogram of the still non
segmented part of the picture is made. If this histogram contains several
significant peaks, the best one is selected.  The points whose features
belong to the range of values of the peak are then selected and the region
that they form is smoothed, its small holes eliminated.  This region is saved
(but can be checked again for further segmentation) and the process
iterated.  Another integration of histogram and region growing methods can
be found in Beveridge and al.~\cite{BevGKHR}.
\medskip

Further attempts have been made for the histogram analysis of multispectral images (Goldberg and S. Shlien~\cite{GoS}). Now, for multiband images, the histogram becomes multidimensional and therefore enormous.  Thus only very rough methods can be used to identify the main clusters.
\medskip

The histogram based segmentation methods are very promising methods,
whose correct theory is still to build.  They have the big advantage of
making globas statistics on the image, where most segmentation methods
remain local.  Among their drawbacks, one can list the following remarks.
\begin{itemize}
\item The histogram methods do not provide any control on the amount and
smoothness of the boundaries of the obtained regions.
\item Since the grey level histograms give an information which is invariant
under any rearrangement of the pixels, they provide no information on the
size or connectivity of the labeled ``regions'' : For instance, a black blob on
bright background can be rearranged into a checkerboard like picture
without changing the histogram.  Therefore, most authors have found it
convenient to reintroduce some local spatial information by localizing the
histograms, or by introducing some gradient information.  Now, local
histograms tend to contain an information which is quite close to that
contained in the successive derivatives of the grey level function.  That is for example  the case in the Ohlander and al. method.  
\item The Panda and Rosenfeld method tends to reintroduce explicitly the gradient in the histogram, which becomes two dimensional.  Then the segmentation problem is led back to a new segmentation problem, this time in the measurement
space!  Indeed, the homogeneous regions will be identified by segmenting
the low gradient region of the two dimensional histogram.\\
To summarize, most histogram-based methods translate the image segmentation problem into a histogram segmentation problem.  If this histogram is one
dimensional, there is some reduction of the difficulty.  The approach is
more delicate if the image is defined by several channels, since the
complexity of the clustering problem increases quickly with the dimension.  
\item Since no control is made on the boundary smoothness by this methods, a shrink operation and a fill operation have to be made on the segmentation
to eliminate the spurious boundaries.  See the survey of Haralick and
Shapiro~\cite{HarS} for the obtained pictures, before and after this complementary operation.  Thus a complete and useable segmentation process by the
histogram techniques needs two steps : one for identifying the homogeneous
regions, and one successive to create acceptable boundaries and eliminate
the ``holes'' and spurious regions.  
\end{itemize}


\subsection[La th\'eorie de Mallat et Zhong]{Mallat and Zhong theory : edges as a complete representation.}

In a recent paper, S.G. Mallat and S. Zhong~\cite{MallZ} follow David Marr's
suggestion that the edge representation at several dyadic scales should be a
complete representation of a picture.  They propose an algorithm, based on
some stability properties of the wavelet transform, which indeed
reconstructs the picture from the edges at five scales.
 
An interesting and suggestive fact of this algorithm is that if one removes the ``spurious edges'' (e.g. the edges which have low gradient or cannot be prolongated), then the reconstruction is still almost perfect.  This justifies a posteriori the importance given to edge detection and enhancement in the litterature.  An important application of this algorithm is selective
smoothing.  Indeed, these authors are able to clean up the picture by
eliminating the ``spurious'' edges (i.e. the edges which are not stable across
several scales) and thereafter reconstructing the picture.
\medskip

Another work worth mentionning on this theme is by Hummel and Moniot \cite{HumM}, where another algorithm for image reconstruction from edges is proposed.  They give a mathematical proof of this reconstruction property, but only for continuous
scale edges. Their proof generalizes Logan's theorem, which was proved for
one dimensional signals.
