

\chapter{Introduction}

This book adresses a possible theory  of image {\it low level} analysis.  Image
``low  level'' analysis  aims  at  extracting   reliable,  local  geometric
informations  from  a digital  image.  Such  informations  are often  called  ``features" and they are used in  order to compare an image to other images. For
instance,  these features can  be used  for motion  estimation, or  to retrieve
shapes, or  to build the still hypothetic  ``high level"  vision. The observed
image is the result of a smoothing of the original photon flux and is therefore
continuous. It is  nontheless well admitted that the subjacent  ``real image",
namely the  focused photon flux,  is either a  measure or, for  more optimistic
authors, a function which presents  strong discontinuities. Rudin and de Giorgi
proposed independently in 1984 the space BV of functions with bounded variation
as  the right  function  space for  ``real"  images.   More recently  (1999),
however, Gousseau  and Alvarez used a  statistical device on  digital images to
estimate how their real subjacent  images oscillate. They deduced, by geometric
measure arguments,  that the ``real" physical  images have in  fact unbounded
variation.  We may therefore accept the idea that the subjacent high resolution
image  behaves in  a strongly  oscillatory  way.  Although  the digital  images
present a averaging of this  oscillatory phenomenon, common sense tells us that
they must  have anyway strong discontinuities at  transitions between different
observed objects,  i.e. on  the apparent contours  of physical objects.  The BV
space looked  at first well adapted  to that aim because  it contains functions
having step discontinuities.
\newline
One of the  goals of image analysis has ever been  to find such discontinuities
in an  image. This  search is called  ``edge  detection" because  early vision
research played with  images of cubes. Along the edges of  the cubes, the light
intensity behaved, in a first approximation, as a step function. Unfortunately,
the early research in vision led to the sad discovery that one could find edges
``everywhere"  in a  digital image  (Marr, Vision),  due to  the oscillations
remaining  in  the  digital  image   after  the  digitization  problem.   As  a
consequence, the image  analysis process was conceived as  a smoothing process,
permitting to  decluster the  true ``edges" from  the inherent noise.   As in
Distribution theory, a smoothing was necessary before computing any derivative.
This is why the heat equation was  proposed and a new doctrine proposed : the
``scale space". Scale space  means that, instead  of talking of features  of an
image at a given  location, we talk of them at a given  location and at a given
scale.  The scale quantifies  the amount  of smoothing  performed on  the image
before computing the feature.  We shall therefore see in experiments ``edges at
scale 4" and ``edges at scale 7" as different outcomes of an edge detector.
\newline
\newline
Which kind of smoothing should be performed ? Three terms associated with image analysis operators arise here, to which we will give a more and more precise meaning.
\newline
The first one, ``locality", is related to the occlusion problem :  Most optical  images are made of a superposition of different objects partly hiding each other.   It is plain that we must avoid mixing them in the analysis, as would do e.g.  a wide convolution. Thus, the analysis must be made as local as possible. As we shall see, the heat equation is the worst candidate to the task, since it makes a wide-range m\'{e}lange of grey levels.
\newline
The second key word is ``iteration".  Indeed, we shall see that it is generally better, from the locality viewpoint, to iterate a very local smoothing operator than to apply it directly at a large scale.  This is precisely not true for the heat equation !  Iterating the convolution of small gaussian kernels is stricly the same as convolving directly the image with a big gaussian.  Now, iteration  of very local filters will bring a significant improvement for some of the most relevant nonlinear filters which we shall consider, namely  the median filter and the affine erosion-dilations.  At this point, it must be immediately announced that the combination of smoothing,  locality and iteration implies that we are talking about parabolic partial differential equations.   This announcement is heuristic and the object of the book is precisely to formalize the necessity and the role of several P.D.E.'s in image analysis.
\newline
Our last key word is ``invariance".  The invariance requirements play a central role in image analysis because the objects to be recognized have to be recognized under varying conditions of illumination (contrast invariance) and from different points of view (projective invariance).    Contrast invariance is one of the key requirements of a famous image analysis theory, the Mathematical Morphology (Matheron, Serra).  This theory proposed a list of contrast invariant image analysis operators (dilations, erosions, median filters, openings, closings,...)  We shall involve this theory, as we shall attempt to localize as much as possible the ``morphomath" operators to extract their behaviour at small scale, and then iterate them.  As an outcome, we shall prove that several geometric partial differential equations, namely the curvature motions, can be considered as the common asymptotic denominators to many ``morphomath" operators.  These P.D.E.'s permit therefore to fuse the Scale Space doctrine and Mathematical Morphology. In particular, affine invariant morphomath operators, which looked  unpractical, turn out to yield in their local iterated version a very affordable P.D.E., the so called ``affine morphological scale space" (A.M.S.S.).
\newline In the next section, we shall make a survey of most  P.D.E.'s which have been proposed for image analysis (Section 1) and thereafter give a detailed mathematical account of how this book is organized.
\newline\newline
We would like to end this short foreword with a light warning to the reader.  We do not claim that what will be developped here is a necessary future for image analysis.  What the mathematical analysis can provide is more hypothetical, as noticed Von Neumann and Morgenstern  in their book on game theory.  We would say, by imitation : ``\underline{If} Image Analysis requires a smoothing theory, then here is how it should be done, and here is the proof that there is no other way to do it."  This statement does not exclude the possibility of other theories, based on different principles, or even the impossibility of making any theory.
\newline
We have tried to prove every single mathematical statements, assuming only a two or three years mathematical training.  Thus, most of the P.D.E's adressed, and all of the relevant ones, will be proved existence, uniqueness and given invariant monotone approximation schemes.
This has been technically possible by combining tools from the recent and remarkbly simple theory of viscosity solutions on the one side, and of the Matheron formalism for monotone set and function operators.  Thus, the really necessary mathematical knowledge amounts to elementary differential calculus, linear algebra, and some notions of integration in the chapter on the heat equation. We nether put mathematical statements {\it pour l'amour de l'art} ; all of them are directed at proving the rightfulness of a model, of its properties and of the associated numerical schemes.  Many experiments are spread out in the text, with detailed comments. They can provide an independent, parallel, reading to the mainstream of the text.  There are no direct industrial applications treated. Now, many of the image representation tools commented in the text are being used in the public software MegaWave (http://www.cmla.ens-cachan.fr) developped jointly by several university research groups in France and Spain. They are also extensively used in the industrial video and image processing software Investigator, sold by Cognitech, Inc., Pasadena, California.

 \section {The when, why and how of P.D.E. models in image processing}

Image processing is divided in three parts, corresponding to as many different goals.  The first one derives from the discrete nature of images and the search of their minimal  representation in terms of digital memory.  This discipline is called {\it image compression} (see Figure \ref{IMGPREFACE1}). \PPRECOMPRESSION  The second goal is the {\it restoration} of a better version of an image, given a generation model with noise and blur, or other perturbations.  This is illustrated in Figure \ref{IMGPREFACE2}. \PPREDEBRUITE  The image on the right is apparently destroyed: more than 75\% of its pixels have been put at a random value.  We can nontheless restore it significantly: here is, on the right,  such a restored version.
 The third goal  is analysis , which means in Greek  ``breaking into parts".  Look at the level curve of Figure \ref{IMGPREFACE3}, extracted from a hand image: it is full with a mix of details and noise. What if we ask for a sketchy version, where, however, all essential details are kept ?  The curve on the right is such a sketchy version, where most  of  the spurious details have disappeared,  but the main structures are maintained.  This is what we shall call  {\it image analysis }. The aim is not denoising or compression : it is to construct an invariant code putting in evidence the main parts (here, for instance, the fingers) and permitting a fast recognition in a large database of shapes.
\PPREANALYSE
\newline\newline
The heat equation arises naturally in the image generation process.  Indeed, according to Shannon's theory, an image can be correctly represented by a discrete set of values, the ``samples", only if it has been previously smoothed. This is illustrated in Figure \ref{IMGPREFACE4} : Let us call the original baby image  ``Victor".  If we attempt to reduce the size of Victor by a mere subsampling, that is by taking a point of each sixteen, we obtain a new and smaller image, in which the subsampling has created  new and unstable patterns : see how new stripes have been created, with a frequency an direction which has nothing to do with the original !  If, instead of being steady, the camera moved, those newly created patterns would move and flicker in a totally uncontrolled way.  This kind of  moving pattern appears often  in recent commercial DVD's. They have simply been subssampled against the Shannon rule.  Let us now comment briefly how the subsampling should be done.  According to Shannon's theory, a previous smoothing must to be done before the subsampling.  We start with $u_0$, the original image.  Then a blur  kernel $k$ is applied, i.e. we convolve $u_0$ with $k$ to obtain a new image $k\ast u_0$.  A subsequent subsampling is thereafter possible, where the distance between samples is related to the band-width of the blur kernel by the Nyquist rule. Stability of the image representation is maintained.
\newline \PSHANNON
This simple remark, that smoothing is a necessary part of image formation, leads us to our first PDE's. Gabor remarked in 1960 that the difference between the original and the blurred image is roughly proportional to its laplacian.  In order to formalize this remark, we have to notice that  $k$ is spatially concentrated, and that we may introduce a scale parameter for $k$, namely $k_h(\fxx)=\frac 1{h}k(\frac  {\fxx}{h^\frac 12})$. Then
$$\frac {u_0\ast k_h(\fxx) -u_0(\fxx)}h  \to\Delta u_0(\fxx),$$
so that when $h$ gets smaller, the blur process looks more and more like the heat equation
$$\frac{\partial u}{\partial t}=\Delta u, \quad u(0)=u_0.$$
\PVICBLUR
Conversely, Gabor deduced that we can, in some extent, deblur an image by reversing time in the heat equation :
$$\frac{\partial u}{\partial t}= -\Delta u, u(0)=u_{observed}.$$
Numerically, this amounts to substracting its laplacian from the observed image :
$$u_{restored}=u_{observed}- h\Delta u_{observed}.$$
This operation can be repeated several times with some small values of $h$, until it... blows up.  Indeed, the reverse heat equation is extremely ill-posed. All the same, this Gabor method is efficient and can be applied with some success to most digital images obtained from an optical device.  Let us examine what happens with Victor (Figure \ref{IMGPREFACE6}). \PVICDEBLUR
We see that the method yields some improvement at the beginning and then blows up.  We can also simulate a blur on Victor and try to go back : again, the process blows up but yields a significant improvement at some scale.
\newline\newline
\PVICBLURDEBLUR

We therefore see two directions.  One is to improve, to stabilize, the reverse heat equation. We shall see that this is doable by nonlinear models.  The second direction is to go on with the heat equation : we can numerically simulate a further blurring of the image.  Why should we do so ?  Because, first, this leads to the wavelet theory and its applications to optimal multiscale sampling and compression. Second, iterated linear and nonlinear smoothing (that is, nonlinear PDE's) will be relevant to our main goal : image analysis.

We can indeed improve the time-reverse heat equation.  The first example, due to Rudin and Osher in 87 and 92  proposes an pseudoinverse, where the propagation term |Du| is tuned by the sign of the laplacian.
$$\frac {\partial u}{\partial t} = - sign{\Delta u} |Du|.$$
The equation is called ``shock filter".
As we shall see, this equation propagates, with a constant speed, the level lines of the image in the same direction as the reverse heat equation would do.  It therefore enhances the image.   The equation is more or less equivalent to a good old nonlinear filter due to Kramer in the seventies.  Kramer's filter can be interpreted as a partial differential equation, by the same kind of heuristic arguments which Gabor developed to derive the heat equation.  This equation is
$$\frac {\partial u}{\partial t} = - sign{D^2u(Du, Du)} |Du|.$$
Thus, the laplacian is replaced by a directional second derivative of the image,
$D^2u(Du, Du)$. We shall later on interpret this differential operator as an ``edge detector".  Kramer's equation yields a slightly better version of shock filter as is illustrated in Figure \ref{IMGPREFACE10}. \PVICKRAMER
Both deblurring  equations work... to some extent. They experimentally do not blow up and attain steady states !    The third deblurring method we can mention here is, to our knowledge, the best version.  It poses the deblurring problem as an inverse problem.  Given the observed image $u_0$, we try to find a restored version $u$ such that $k\ast u$ is as close as possible to $u_0$ and the oscillation of $u$ is nontheless bounded :
$$u_{restored}=Argmin (\int |Du| + \lambda (k\ast u -u_0)^2.$$
The parameter $\lambda$ tunes the oscillation we allow for the restored version.  If $\lambda$ is large, the restored version will satisfy accurately the equation
$k\ast u=u_0$, but may be very oscillatory. If instead $\lambda $ is small, we get a smooth but unaccurate solution.  This parameter can be computed in principle as a Lagrange multiplier.
 The obtained restoration can be  remarkable. We display the best result we can obtain with the blurred Victor in Figure \ref{IMGPREFACE10}-right. This total variation restoration method also has fast wavelet packets versions. It recently won a benchmark in satellite image deblurring organised by the French Space Agency (CNES).

The original remark of Gabor, about image generation being related to the Laplacian of the image, leads to the wavelet theory as well.  Here is how it works : if we convolve the image with some smoothing kernel and thereafter make the difference, we obtain a new image, actually a laplacian, which turns out to be faded with respect to the original.  In Figure \ref{IMGPREFACE8},  the  last image on the right shows  in black the values  of this laplacian image of Victor which differ significantly from zero : In most natural images, as here,  this representation is sparse and adapted to compression.    This is why one of the first wavelet representations, due to Burt and Adelson in 83, was called ``Laplacian pyramid".  It boils down to the  iteration of a convolution followed by subsampling.   We only keep the differences between images smoothed at different scales, i.e. their laplacians.  The objective is a compressed representation, but to the price of a loss of invariance due to the multiscale subsampling. \PWAVELETS

In image analysis,  the heat equation has had a very different use: Marr, Hildreth, Canny, Witkin, Koenderink proposed in the eighties to analyse an image by applying the heat equation.  As Rudin and Florack noticed, this is related to distribution theory.  Indeed, details of the images, like boundaries, corners and other singularities cannot be computed without some previous smoothing because they are derivatives of a nonsmooth function.  And this smoothing has to be multiscale because the image  is multiscale !  The heat equation is easily proved to be the only good candidate to the task if image analysis has to be linear. What derivatives should be computed in an image~?
The early research in computer vision proposed  ``edge detection" as a main tool : it is assumed that the apparent contours  of the objects and  also the boundaries of the facets of objects,  result in step discontinuities in the image, while, inside those boundaries, the image oscillates mildly. The apparent contour points, or ``edges points" will be computed as points where the gradient is is some sense largest.  Two ways to do so : Hildreth and Marr proposed the points where $\Delta u$ crosses zero.  A significant improvement was done by Canny, who proposed to compute the points where $Du$ is maximal  on the gradient lines.  Such points satisfy $D^2u(Du, Du)=0$.
Figure \ref{IMGPREFACE9} displays what happens when we smooth the image with the heat equation and compute the points where $D^2u(Du, Du)=0$ and $|Du|$ is large enough.  At first, everything in the image is boundary : the image, being a very oscillatory function, has inflexion points everywhere~!  After some evolution of the heat equation, we can see  what happens : we are able to extract some structure. \PBOUNDHEAT

If the heat equation is, under sound invariance requirements,  the only good linear smoother, there are instead many nonlinear ways to smooth an image.   The first one was proposed by Perona and Malik in 87.  The idea is roughly to smooth out what has to be smoothed, the irrelevant, homogeneous, regions and to enhance instead the boundaries.  Thus, the diffusion should look like the heat equation when $|Du|$  is small and an inverse heat equation should  instead be applied when $|Du|$ is large.  Here is the  equation in divergence form.
$$\frac {\partial u}{\partial t} = div (g( |Du|^2)),$$
where $g(s) = \frac 1 {1+\lambda s^2}$ decreases when $s$ increases. It is easily checked that we have a diffusion equation when $|Du|\leq \lambda$ and
an inverse diffusion equation when $|Du|>\lambda$. In order to do so,
we rewrite the equation in the following way.  We consider the second derivative of $u$ in the direction of $Du$,
$$u_{\eta\eta}= D^2u(\frac {Du}{|Du|}, \frac {Du}{|Du|})$$ and the second derivative in the orthogonal direction,
$$u_{\xi\xi}= D^2u(\frac {Du\perp}{|Du|}, \frac {Du^\perp}{|Du|}),$$
where $Du=(u_x, u_y)$ and $Du^\perp = (-u_y, u_x).$  The laplacian can be  rewritten in the intrinsic coordinates $(\xi, \eta)$ as
$\Delta u = u_{\xi\xi}+u_{\eta\eta}.$
The Perona-Malik equation rewrites
$$\frac {\partial u}{\partial t} = \frac{u_{\xi\xi}}{1+\lambda^2|Du|^2}
+\frac{(1-\lambda^2|Du|^2)u_{\eta\eta}}{(1+\lambda^2|Du|^2)^2}.
$$ \PBOUNDPERONA
So the first term always appears as a one-dimensional heat equation in the direction orthogonal to the gradient, tuned by the size of the gradient though. The second term can be a directional heat equation, or reverse heat equation in the direction of the gradient. So we indeed mix in this model the  heat equation and  the reverse heat equation !
We compare in Figure \ref{IMGPREFACE11} the  Perona-Malik with the classical heat equation in terms of accuracy on the boundaries obtained by Canny's edge detector : at a comparable scale of smoothing, we clearly gain some accuracy in the boundaries and get rid of more ``spurious" boundaries. The representation is both more sparse and more accurate.
\newline\newline

Now, this ambitious model attempts to put in a single operator two  very different goals which we already mentionned, namely restoration and analysis.  This has a cost : the model contains a ``constrast threshold"  which can only be fixed manually. Mathematical existence and uniqueness are not guaranteed, despite some attempts by Kichenassamy and Weickert.  Let us  summarize the involved parameters : we need to fix both $\lambda$ and the smoothing scale(s) $t$ and the threshold on the gradient in Canny's detector as well. We obviously must take the same gradient threshold $\lambda$ in Canny's detector and in the Perona-Malik equation. All the same, we have a two parameters game : how will this be dealt with in automatic image analysis ?  This question seems to have no general answer for the time being. An interesting attempt based on statistical arguments is made, however, in Black, Sapiro.
\newline\newline

If any nonlinear diffusion can be an image analysis model, why not trying them all ? This is exactly what has happened in the past ten years.  We can claim that almost all possible nonlinear parabolic equations have been proposed. The logic in this proliferation of models is this : each attempt  fixes one intrinsic diffusion direction and tunes the diffusion by the size of the gradient or by the size of a nonlinear estimate of the gradient.  Sometimes, the proposed models are even systems of PDE's, but in order not to blow up this preface, we shall focus on the simplest proposed examples.
We can start with Rudin-Osher-Fatemi's model, which consists, for the smoothing term, of minimizing the total variation of $u$.  The gradient descent for $\int |Du|$ writes
$$\frac {\partial u}{\partial t} = div (\frac {Du}{|Du|})=\frac 1 {|Du|} u_{\xi\xi}.$$
Written in that way, the method appears as diffusion in the direction orthogonal to the gradient, tuned by the magnitude of the gradient.  Caselles, ? ? and ? ? proved  that this equation is indeed well posed in the space of bounded variation.
A variant was proposed independently by Alvarez and al.,
$$\frac {\partial u}{\partial t} =\frac {|Du|}{|k\ast Du|} div (\frac {Du}{|Du|})=\frac 1 {|k\ast Du|} u_{\xi\xi},$$
where the tuning of the gradient is nonlocal.  Kimia, Tannenbaum and Zucker proposed, endowed in a more general shape analysis framework, the simplest equation of the list,
$$\frac {\partial u}{\partial t} =|Du|div (\frac {Du}{|Du|})= D^2u(\frac {Du\perp}{|Du|}, \frac {Du^\perp}{|Du|})=u_{\xi\xi}.$$
This equation had been proposed some time before in another context by Sethian as a tool for front propagation
algorithms. This equation, which we call in continuation ``curvature equation", is a ``pure" diffusion in the direction
orthogonal to the gradient. The Weickert equation is a variant of the curvature equation, with  nonlocal estimate of
the direction orthogonal to the gradient : the diffusion direction $d =SEigen (k\ast (Du\otimes Du))$ is computed as
the eigenvector of the least eigenvalue of $k\ast (Du \otimes Du)$ : if the convolution kernel is removed, this
eigenvector simply is $Du^\perp$. The three mentioned models can be interpreted as diffusions in a direction orthogonal
to (an estimate of) the gradient, tuned by the magnitude of the gradient (Figure \ref{IMGPREFACE13}).
\PFIRSTJUNGLEDIFFUSION Other diffusions have been considered as well : For interpolation goals, Caselles et al.
proposed a diffusion which may be interpreted as the strongest possible image smoothing,
$$\frac {\partial u}{\partial t} = D^2u(Du, Du).$$
This equation is not used as the other ones as a preprocessing of the image, but a way to interpolate between the level lines an image with sparse level lines (Figure \ref{IMGPREFACE14}). \PTHIRDJUNGLEDIFFUSION
Zhong and Carmona proposed a diffusion in the direction $d=SEigen(D^2u)$ of the eigenvector with least eigenvalue of $D^2u$ (Figure \ref{IMGPREFACE12}). \PSECONDJUNGLEDIFFUSION
Sochen, Kimmel and Malladi propose instead a nondegenerate diffusion, associated with a  minimal surface variational formulation : their idea was to make a gradient descent for the area of the graph of $u$,  $\int \sqrt{1+|Du|^2}$, which leads to  the diffusion equation (Figure \ref{IMGPREFACE12}).
$$\frac {\partial u}{\partial t} = div (\frac {Du}{\sqrt{1+|Du|^2}}).$$


Among the mentionned models, only the curvature motion was explicitly proposed by Kimia, Tannenbaum and Zucker as shape analysis tool.  We shall now explain why.
\newline
In order to do so, we have to give a definition of  image analysis.  There might be as many ways to define this discipline as they are applicational goals involving digital images.  Now, the range of applications is as wide as the human activity, since most of the scientific and technical human activity, including even sound analysis (visual sonagrams), involves the perceptual analysis of images.   Fortunately, we have at hand a mathematical shortcut to avoid an endless list of partial and specific requirements : This shortcut, well known in Mechanics, consists of stating  invariance requirements.  Invariance requirements will be a short list and they  will, as we shall see, give a possible classification of  models and point out the ones which are the most adequate for all purposes image analysis tools.  The first invariance requirement is the Wertheimer principle according to which visual perception (and therefore, may we add, image analysis) should be independent of the image contrast.  We state this in the following way :
\newline
\newline
{\bf Contrast invariant classes} \newline
$u$ and $v$ are said to be (perceptually) equivalent if there is a continuous increasing function $g$ such that $v=g(u)$.
\newline
\newline
{\bf Contrast invariance requirement} : An image analysis operator $T$ must act directly on the equivalence class.  As a consequence, we may ask that
$T(g(u))=g(Tu)$, i.e. a commutation of the image analysis operator with contrast changes.
\newline
\ newline

The contrast invariance requirement rules out the heat equation and all models stated before, except the curvature motion.
Contrast  invariance led Matheron in seventy-five to reduce image analysis to a set analysis, namely the analysis of level sets.  We call  upper level set with level $\lambda$ of an image $u$ the set
$$\levelset_\lambda u= \{\fxx, u(\fxx)\geq \lambda\}.$$
We can define in exactly the same way the lower level sets, by changing
``$ \geq$'' into ``$\leq$''.  The main point to retain here is the global invariance of level sets under contrast changes, namely, if $g$ is a continuous increasing contrast change, then,
$$\levelset_{g(\lambda)}g(u)=\levelset_\lambda u.$$
According to Mathematical Morphology, this image analysis doctrine founded by Matheron, all of the image shape information is therefore contained in the level sets : it can be proved that an image can be reconstructed, up to a contrast change, from its set of level sets (Figure \ref{IMGPREFACE16} : an image and some of its level sets). \PPRELEVELSET
\newline

The contrast invariance requirement yields powerful and simple denoising operators as the so called ``Extrema killer" defined by Vincent and Serra in 1993.  This image operator simply removes all connected components of upper and lower level sets with area smaller than some fixed scale.  This is not a PDE, actually it's much simpler !  Now, its effect is amazingly good for impulse noise i.e. local destructions of the image, spots.  In Figure \ref{IMGPREFACE15}, we see a image degraded up to 75\%. Below, its restoration by the extrema killer. Left, result of the same operator applied to the original. \PPREFEXTKILL
\newline
\newline
Caselles and Coll localized farther in 1996  the contrast invariance requirement in image analysis.  They proposed as the main object of analysis the level lines of the images, that is, the boundaries of level sets.  This proposition makes sense for a digital image, which is assumed to be a sampling of  a very smooth function as the result of the optical smoothing.  We can therefore define the level lines if, e.g., the interpolated image is $C^1$ as is guaranteed by the canonical Shannon interpolation. There may be other interpolation methods, and even interpolations into a discontinuous functions : this is the case if, e.g., we consider the digital image as constant on each pixel.  We must then for each interpolation method make clear how the level lines are computed and what their structure is.  Two properties are desirable : that the level curves indeed are curves in some affordable sense (Jordan rectifiable curves) and that they are nested, i.e. never cross, so that they make an inclusion tree.
 A study of Kronrod (1950) shows that if the function $u$ is continuous, then the isolevels sets $\{\fxx, u(\fxx)=\lambda\}$ are nested : they build a tree ordered by inclusion.  Now, these isolevel sets need not be really curves.   Monasse (2000) generalized recently the preceding result to lower semicontinuous or upper semicontinous functions.  His result implies that the simplest,  piecewise constant, interpolation of an image yields a nested set of Jordan curves bounding the pixels.  Thus, we have two good ways to associate with the digital image a set of nested Jordan curves.  We call this set ``topographic map".
\footnote{This point of view also is coherent with the ``BV assumption" which we mentionned at the beginning of the introduction, according to which  the right function space for images should be the space BV of functions with bounded variation. By coarea formula, we can then describe  the image by a bunch of Jordan level curves (see Ambrosio et al.) Now, it is in general false for BV functions that boundaries of lower and upper level sets make a nested set of curves : these curves may cross. }
We display in Figure \ref{IMGPREFACE17} the level lines of a digital image at some fixed level.  By the introduction of the topographic map, the search for  image smoothing, which we had already reduced to set smoothing,  is further reduced to  curve smoothing,  provided of course this smoothing preserves curve  inclusion.  \PPRELEVELLINES
\newline
Chen, Giga and Goto and Alvarez, et al.  proved that, under the usual invariance requirements for image processing, including the contrast invariance, all image multiscale analyses should have the form of a curvature motion, namely
$$ \frac {\partial u}{\partial t} = F(curv (u), t) |Du|,$$
where $F$ is increasing with respect to its first argument.
This equation can be interpreted as this : we consider a point $\fxx$ on a given level curve of  $u(t)$, at time $t$.  We call $n(\fxx)$ the normal vector to the level curve and $curv(\fxx)$ its curvature.  Then the preceding equation turns out to be associated with the curve motion equation
$$\frac {\partial \fxx}{\partial t} = F(curv(\fxx)) n(\fxx),$$
which describes how the point $\fxx$ moves in the direction of the normal.
Not much more can be said at this level of generality on $F$.  Now, two particular cases happen to play a prominent role.  First, the case $F(curv(u), t)=curv(u)$, the so called curvature equation which we already met, and second the case $F(curv(u), t) = curv(u)^\frac 13.$
\newline\newline
This particular form for the curvature dependence, the power one third, permits to get a very relevant additional invariance, the affine invariance.  We would like to have a full projective invariance, but a theorem proved by Alvarez et al.  shows  that this is impossible.  The best we can have is invariance with respect to the so called chinese perspective, which preserves parallelism.   Most of the mentionned equations, particularly when $F$ is a power of curvature,  have a viscosity solution in the sense of Crandall and Lions, as shown by recent works of Ishii and Souganidis.
\newline\newline
As we already mentionned, contrast invariant processing boils down to level set, and finally level curve processing. The above mentionned equations indeed are equivalent to curve evolution models, provided strong existence results are at hand. This is the case for the most important cases, namely the power 1, the so called ``curve shortening" and the power 1/3, known as ``affine shortening".   Grayson proved existence, uniqueness and analycity for the first equation,
$$\frac {\partial \fxx}{\partial t} =curv(\fxx)n(\fxx)$$
 and Angenent, Sapiro and Tannenbaum for the affine shortening
$$\frac {\partial \fxx}{\partial t} = curv(\fxx)^\frac 13 n (\fxx).$$
\PPREAMSS
Those results are very relevant to image analysis as they ensure that the diffusion process indeed reduces the curve to a more and more sketchy version. We check  the  affine invariance in Figure \ref{IMGPREFACE21} .  The numerical test we shall make here is as follows : we apply a very fast and fully affine invariant numerical scheme designed by Lionel Moisan.  In the middle, the initial shape is an affine transform of the first one ; the shape on the right will be an inverse affine transform of the middle shape.   If everything is correct, we can expect that, after processing, the shape on the right will be identical to the shape on the left.  We make the experiment with both the the curve shortening and the affine shortening.   So, it works !
\newline
\newline \PPRELIONEL
Evans-Spruck and Chen-Giga-Goto proved in 1991 that  a continuous function moves by curvature motion if  and only almost all of its level curves move by curve shortening. This yields, in that case, a mathematical justification of the now classical Osher-Sethian numerical method for moving fronts by moving a distance function to the front.    The same result is true for the affine invariant curve evolution.  The Osher-Sethian point of view is just converse to the point of view adopted here : they associate with some curve $C$  or surface its signed distance function $u$ , so that the curve or surface is handled indirectly as the zero isolevel set of $u$. Then $u$ is evolved by, say, the curvature motion with a classical numerical difference scheme.  In that way, the curve evolution is dealt with efficiently and accurately.  From our point view, the image can be viewed as a distance function to all and each of its level sets, since we are interested in all of them.
\newline   \PCOURBUN
We show in Figure \ref{IMGPREFACE19}  an application  of this numerical method, with both curvature and affine invariant curvature motions.  In order to gain visibility, we do not display all level curves, but only for about eighteen levels.  Notice that the aim is not here subsampling ; we keep the same resolution.  It is not either restoration : the processed image is clearly worse than the original.  The aim is invariant simplification leading to shape recognition.

\PCOURBDEUX

Before proceeding to shape recognition, let us mention that a well adapted variant of curvature equation can be used for shape detection. It's a by now famous method of contour detection in an image, initially proposed by Kass, Witkin and Terzopoulos.   This method was very unstable and the winning method turns out to be a variant of curvature motion proposed by Caselles, Catt�, Coll, Dibos and improved simultaneously  by Caselles, Kimmel, Sapiro, and Malladi, Sethian.  Here is how it works.  The user draws roughly what are the contours he wants in the image and the algorithm then finds the best possible contour in terms of some variational criterion.  This turns out to be very useful in medical imaging.   The motion of the contour is a tuned curvature motion which tends to minimize the energy $E$ which we will now explain.
Given an original image $u_0$ containing some circular contours which we wish to approximate, we start with an ``edge map"
$$g(\fxx)=\frac 1{1+|Du_0(\fxx)|^2},$$
that is, a function which vanished on the edges of the image.  The user then points out the contour he is interested in,  by drawing a polygon $\gamma_0$ surrounding roughly the desired contour.  The ``geodesic snake" algorithm then builds a distance function $v_0$ to this initial contour, so that $\gamma_0$ is the zero level set of $v_0$.  The energy to be minimized is
$$E(\gamma)=\int_\gamma g(\fxx(s))ds,$$
where $g$ is the edge map associated with the original image $u_0$ and $s$ denotes the length parameter on  $\gamma$.  The motion of the ``analysing image" $v$ is governed by
$$\frac {\partial v}{\partial t} =g |Dv|curv(v) - Dv.Dg.$$
We display an example in Figure \ref{IMGPREFACE22}.  \PPRESNAKE
\newline
\newline
The main obvious application of invariant PDE's seems to be shape retrieval in large databases.  There are thousands of different definitions of shapes, and of shape recognition algorithms.  Now, the real bottleneck has ever been extraction of the relevant shapes.  The discussion above points to a brute force strategy : all contrast invariant local elements, are are the level lines of the image, are candidates to be  ``shape elements".  Of course, this name of shape element suggests the contours of some object, but there is no way to give a simple geometric definition of objects.  We must give up the hope of jumping from the geometry to the common sense world.  We may instead
 simply ask the question :  given two images, can we retrieve all level lines similar in both ? This would give a factual, a posteriori definition of shapes : they would be defined as pieces of level lines common to two different images, no matter what their relationships to real physical objects are.  \PPRELISANI
Of course, this brute force strategy would be impossible without the  previous invariant filtering (AMSS).   It is instead doable if the level lines have been significantly simplified.  This simplification entails the possibility of compressed invariant encoding. In Figure \ref{IMGPREFACE20}, we present an experiment due to Lisani et al.. Two images of a desk  are taken from different angles, and then, in white,  all of the pieces of level lines in Image 1 and in Image 2 which found a match in the other image.   In continuation, we present some of the matches. We notice that several of these matches are doubled : indeed, there are two similar chairs in each images ! A Gestalt law comes immediately to mind. This law states that human perception  tends to group similar shapes. We now see the numerical necessity of this perceptual grouping : a previous self-matching of each image, with grouping of similar shapes, must be performed before we can compare it to other images !



\section{The mathematical organisation}
The mathematical organisation of this book derives from the above discussion of PDE models.  We chose to focus on image analysis models, so that invariance requirements will lead the mathematical analysis.
Most chapters include numerical experiments and indications on the right way to make them.
\begin{description}
\item {Chapter \ref{chapter1}} We extensively analyse the heat equation, because it is  a useful mathematical and algorithmic tool throughout the book. We first prove existence and uniqueness results for the heat equation and then prove that iterated linear smoothing filters are asymptotically equivalent to the heat equation.  This is a first way to show its uniqueness as a linear smoothing filter.

\item{Chapter \ref{chapter2}} We then use the heat equation as a tool to :
\begin{itemize}
     \item edge detection (Hildreth-Marr, Witkin)
     \item set smoothing (Koenderink-Van Dorn, Bence-Merriman-Osher)
    \item curve smoothing (Mackworth-Mokhtarian)
\end{itemize}

\item {Chapter \ref{chapter3}} We prove one of the main principles of mathematical morphology : a contrast invariant class of images is completely described by its  set of level sets. Thus, contrast invariance reduces image analysis to set analysis.

\item {Chapter \ref{chapter4}} Level lines and level surfaces are introduced and the definitions and main formulae  of curvature of level lines and, in higher dimension, of the principal curvatures of an isolevel surface given and proved. We also give a definition of level lines for discontinuous digital images and adress their visualization.

\item {Chapter \ref{chapter5}} We prove that
monotone image operators that are contrast invariant boil down to set monotone operators. We explain how to define the set operator from the image operator and, conversely, how to construct a contrast invariant  image operator from a set operator.  Commutation of the monotone contrast invariant operators with thresholds is proved.
As a first application, the Vincent-Serra ``extrema killer" is formally defined.

\item{Chapter \ref{chapter6}}  An analytic form is given to monotone contrast invariant  and translation invariant image operators : If $T$ is such an operator, it has a canonical ``sup-inf" form,
$$Tu(\fxx) =\sup_{B\in \setB} \inf_{\fxy\in B}u(\fxx+\fxy),$$
where $\setB $ is a set of subsets of $\real^2$ called ``set of structuring elements".  Conversely, every image operator in the sup-inf form is monotone and contrast invariant.
This theorem is a nonlinear version of the Riesz theorem, according to which a linear, continuous, translation invariant  operator can be associated a convolution kernel $k$, so that
$$Tu(\fxx)=\int k(\fxy)u(\fxx-\fxy) d\fxy.$$
The kernel $k$ is called ``impulse response". In the same way, $\setB$ is an impulse response for the nonlinear operator.

\item {Chapter \ref{chapter7}}
The simplest sup-inf operators of mathematical morphology, the so called ``dilations (only sup)" and ``erosions (only inf)" are presented and analysed.   Their subjacent PDE are identified as propagation equations
$$\frac {\partial u}{\partial t} = c|Du|,$$
where $c=1$ for dilations and $c=-1$ for erosions.

\item{Chapter \ref{chapter8}}  The median filter is one of the most emblematic and efficient contrast invariant monotone operator.  It is defined and its numerical implementation discussed.

\item{Chapter \ref{chapter9}} This chapter focuses on the asymptotic expansion of contrast invariant operators when their scale tends to zero.  We associate with each operator $T$ its scaled version $T_h$ associated with the shrinked set of structuring elements $\setB_h=\{hB, B \in \setB\}.$  It is then proved roughly that $$T_hu(\fxx)-u(\fxx)= ch |Du|+ hH(h curv(u))|Du|,$$
where $c$ is a constant and $H$ a nondecreasing function, both associated with
$T$ by simple formulas.  In the case of the median filter it is proved that $c=0$ and $H(s)= s$.  We therefore have two main classes of  contrast invariant operators  :
\begin {itemize}
\item $c\neq 0$ :the operator has the same asymptotic behaviour as a dilation if $c>0$ and an erosion if $c<0$.
\item $c=0$ : the operator corresponds to a curvature dependent motion.
The overall conclusion of this chapter is that all contrast invariant monotone operators boil down to erosions, dilations, and curvature motions.
\end{itemize}

\item{Chapter \ref{chapter10}} : performs the same analysis in arbitrary dimension.

\item {Chapter \ref{chapter11}} Because of the affine invariance requirement, a special attention is paid to   affine invariant operators : an affine invariant distance of a point to a set is defined.  In continuation, affine dilation and erosion operators are deduced and their set of structuring elements identified.

\item{Chapter \ref{chapter12}} Affine invariant families of structuring elements are not bounded, since by an affine map with determinant 1 we can arbitrarily stretch any set of the family.  It is shown that the asymptotic behaviour of the affine invariant operator is not altered if we impose adequate bounds on the affine structuring elements.

\item {Chapter \ref{chapter13}} It is shown that if $T$ is an affine invariant operator such that $T(-u)=-T(u)$, then  its asymptotic behaviour is
$$T_h u(\fxx)-u(\fxx)= ch^\frac 43 |Du| curv(u)^\frac 13.$$
In summary, all affine invariant contrast invariant operators are asymptotically equivalent to the affine PDE,
$$\frac {\partial u}{\partial t} = |Du| curv(u)^\frac 13.$$

\item{Chapter \ref{chapter14}} : adresses the extension of the preceding theory to the so called ``nonflat" mathematical morphology, i.e. when the operators are no more contrast invariant, but only monotone.  It is shown that an asymptotic expansion of a monotone operator can lead either to any of the Hamilton-Jacobi equations
$$ \frac {\partial u}{\partial t} =H( Du)$$
or to parabolic equations (e.g. the heat equation).  Actually, all  of the intrinsic parabolic equations can derived from the Mathematical Morphology framework !
As applications, we analyse Rudin-Osher and Kramer's shock filters.

\item{Chapter \ref{chapter15}}  The curvature equations are not smooth enough to be given classical solutions : we have a need for a concept of solutions compatible with contrast changes. Now, if $u$ is a solution of the curvature equation and $g$ a continuous increasing nondifferentiable contrast change, $g(u)$ should be another solution and is clearly not more than continuous !
The concept of viscosity solutions in the  framework  of Crandall-Lions, permits to solve this difficulty.  We define this concept of solution and prove the main needed properties.  In particular, we check that the erosions are dilations are viscosity solutions of the propagation equations $$\frac {\partial u}{\partial t} =c |Du| $$. We also prove a theorem, due to Barles and Souganidis,
which permits to prove that the iteration of a contrast invariant operators $T_h$ yields a discrete approximation of the subjacent curvature equation.
We also give a brief account of the uniqueness theory for viscosity solutions.

\item{Chapter \ref{chapter16}} is devoted to the application of the mathematical techniques developped in the former chapter.  We consider the main relevant contrast invariant operators, namely the median filter and the affine invariant erosion-dilation.  We consider adequate rescaled versions of both, which we denote generically by $T_h$.  Then we prove that $(T_h)^nu_0\to u(t)$ when
$nh\to t$ and $u(t)$ is a viscosity solution of the associated PDE,
$\frac {\partial u}{\partial t} = |Du| curv(u)^\frac 13$ with $u(0)=u_0$, when $T $ is an affine contrast invariant operator and $\frac {\partial u}{\partial t} = |Du| curv(u)$ when $T$ is a median filter.   Since the viscosity solution  of these equations has been proved to be unique, we have both constructed this unique solution and proved that the considered iterated filter converges towards this unique solution !

\item {Chapter \ref{chapter17}} is an application of the results and techniques of the former chapters to the classical ``active contour" or ``geodesic snakes" method.  We explain this method, show its structural properties and prove existence, uniqueness and contrast invariance of the  motion of the analysing function $v$.


\item {Chapters \ref{chapter18} and \ref{chapter19} : image scale space theory} introduces a different point of view on image analysis, namely the Scale Space  theory.  A scale space is family of smoothing operators
$T_t$, depending upon a scale parameter $t$, which associate with the original image $u_0$ more and more sketchy images $u(t)$. The original scale space theory proposes to compute $u(t)$ as the solution of the heat equation with $u_0$ as initial condition. Clearly, the viscosity solutions $u(t)$ of curvature equations give other scale space theories, more invariant. Now, aren't we missing other possibilities ?  The axiomatic analysis performed in Chapters 18 and 19 permits first to explain which requirements (locality, invariance, comparison principle...) lead to a parabolic P.D.E.  and  then to classify all possible P.D.E.'s according to their invariance properties.  The affine invariant equations are given in any dimension ; in dimension 2, the AMSS model (power
$1/3$ is proved to be the only possibility.

\item{Chapter \ref{chapter20}} A similar axiomatics analysis is performed on shape Scale Spaces. Here are also stated (without proof for once) the existence, uniqueness and regularity results for the curve evolutions to the powers 1 and $1/3$. This shape analysis axiomatics turns out to be simpler than the general image scale space  axiomatics, because it is based on the shape inclusion principle.
It gives also a formal justification of the algorithm which, instead of moving the image by a  curvature P.D.E., moves all level lines by the corresponding curve evolution equation.  It is shown that the viscosity solution for the first equation indeed is obtained by moving all level lines with the second
\item{Chapter \ref{chapter21}} : Scale Space of movies. Here, a new axiom is introduced, the galilean invariance, which permits to single out a single most invariant P.D.E. for movies.

\item{Chapter \ref{chapter22}}  is devoted to several different numerical strategies with finite difference schemes to implement image evolution by a curvature motion.

%%\item{Chapter 23} gives a formal justification of the algorithm which, instead of moving the image by a  curvature P.D.E., moves all level lines by the corresponding curve evolution equation.  It is shown that the viscosity solution for the first equation indeed is obtained by moving all level lines with the second.
\end{description}



\ifodd\count0
    \clearpage
    \leavevmode
    \thispagestyle{empty}
    \clearpage
\else\relax \fi
