%gmIntro.tex

\newcommand{\add}{}
\newcommand{\remove}[1]{}
\newcommand{\comment}[1]{}
\newcommand{\makeexercise}[1]{}

 %Si on enlève [1], la commande remove n'agit pas et le texte apparaît
 %même  chose pour comment et makeexercise


\comment{\bf I am putting labels to each chapter and changing
explicit references into labels. This is important since the
chapter numbering will change. Also, this means that we'll have to
revise the  whole book searching for "chapter".}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\chapter*{Introduction}

\addcontentsline{toc}{chapter}{Introduction}
\markboth{Introduction}{Introduction}

\renewcommand{\thesection}{{I}.\arabic{section}}
\renewcommand{\thefigure}{{I}.\arabic{figure}}
\renewcommand{\theequation}{{I}.\arabic{equation}}
 This book addresses the problem of low-level image analysis and,
as such, is a contribution to image processing and imaging
science. \remove{\it If there is to be a comprehensive theory of
image analysis, then the material developed here is a possible
candidate for such a theory.} While the material touches on
several aspects of image analysis---and peripherally on other
parts of image processing---the main subject is image smoothing
using partial differential equations (PDEs). The rational for a
book devoted to smoothing is the assumption that a digital image
must be smoothed before reliable features can be extracted.

The purpose of this introduction  is to establish some of the
language, conventions, and assumptions that are used throughout
the book, to review part of the history of PDEs in image
processing, and to introduce notation and background material.

\section{Images}

Since the objects of our study are ultimately digital images,  we
begin by defining what we mean by ``digital image'' and by
describing some of the ways these images are obtained and some
current assumptions about the ``original images'' from which the
digital images are derived.

Most of the images dealt with will be natural images, that is,
images from nature (people, landscapes, cityscapes, etc.). We
include medical images and astronomical images, and we do not
exclude drawings, paintings, and other man-made images. All of the
images we consider will be grayscale images. Thus, mathematically,
an image is a real-valued function $u$ defined on some subset
$\Omega$ of the plane $\mathbb R^2$. The value $u(\bm x)$,
$\bm{x}=(x, y)\in\Omega$, represents the gray level of the image
at the point $\bm x$. If $u$ is a digital image, then its domain
of definition is a finite grid with evenly spaced points. It is
often square with $2^n\times2^n$ points. The gray levels $u(\bm
x)$ are typically coded with the integers 0--255, where 0
represents black and 255 represents white. If $h$ is the distance
between grid lines, then the squares with sides of length $h$
centered at the points $u(\bm x)$ are called {\it pixels}, where
``pix'' is slang for ``picture'' and ``el'' stands for
``element.''

The mathematical development in this  book proceeds along two
parallel lines. The first is theoretical and deals with images $u$
that belong to function spaces, generally spaces of continuous
functions that are defined on domains of $\mathbb R^2$. The second
line concerns numerical algorithms, and for this the images are
digital images. To understand the relations between the digital
and continuous images, it is useful to consider some examples of
how images are obtained and some of the assumptions we make about
the processes and the images. Perhaps the simplest example is that
of taking a picture of a natural scene with a digital camera. The
scene---call it $S$---is focused at the focal plane of the camera
forming a representation of $S$ that we denote by $u_f$. When we
take the picture, the image $u_f$ is sampled, or captured, by an
array of charged coupled devices (CCDs) producing the digital
image $u_d$. This image, $u_d$, is the only representation of $S$
that is directly available to us; the image $u_f$ is not directly
available to us. Even more elusive is the completely hypothetical
image that we call $u_S$. This is the representation of $S$ that
would be formed at the focal plane of an ideal camera having
perfect optics. A variation on this example is to capture $u_f$ on
film as the image $u_p$. Then $u_p$ can be sampled (scanned) to
produce a digital image $u_d$. For example, before the advent of
CCDs, astronomical images were captured on Schmidt plates. Many of
these plates have been scanned recently, and the digital images
have been made available to astronomers via the Internet.

Aspects of the photographic example could be  recast for medical
imaging. Although photography plays an important role in medicine,
images for diagnostic use are often obtained using other kinds of
radiation. X-rays are perhaps closest to our photographic example.
In this case, there is an image corresponding to $u_p$ that can be
scanned to produce a digital image $u_d$. Other medical imaging
processes, such as scintigraphy and nuclear magnetic resonance,
are more complicated, but these processes yield digital images.
The images examined by the experts are often ``negatives''
produced from an original digital images. Irrespective of the
process, digital images captured by some technology all have one
characteristic in common: They are all noisy.

One way to relate the different representations of $S$, is to write
\[
u_d=Tu_S+n,
\]
where $T$ is a hypothetical operator representing some technology
and $n$ is noise. In the case of photography, we might write this
in two steps,
\begin{equation*}
\begin{cases}
u_f=P*u_S+n_1,\\
u_d=Ru_f+n_2,
\end{cases}
\end{equation*}
where $P$ represents the optics and $R$ represents the sampling.
This is a useful model in optical astronomy, since astronomers
have considerable knowledge about the operators $P$ and $R$ and
about the noises $n_1$ and $n_2$. Similarly, experts in other
technologies know a great deal about the processes and noise
sources. Noise and pixels are illustrated in Figure \ref{biga}

In the photographic example, the image $u_f$ is a smoothed version
of $u_S$. Furthermore, $Ru_f(\bm{x})$ is not exactly $u_f(\bm{x})$
but rather an average of values of $u_f$ in a small neighborhood
of $\bm{x}$, which is to say that the operator $R$ does some
smoothing. Thus, in this example, $u_d$ is sampled from a smoothed
version of $S$. We are going to assume that this is the case for
the digital images considered in the book, except for  digital
images that are artificially generated. This is realistic, since
all of the processes $T$ that we can imagine for capturing images,
smooth the original photon flux. In fact, this is more of an
observation about technology than it is an assumption. We are also
going to assume that, for any technology considered, the sampling
rate used to produce $u_d$ is high enough so that $u_d$ is a
``good'' representation of the smoothed version of $S$, call it
$u_f$, from which it was derived. Here, ``good'' means that the
parallel development in the book mentioned above make sense; it
means that, from a practical point of view, the theoretical
development that uses smooth functions to model the images $u_f$
is indeed related to the algorithmic development that uses the
digital images $u_d$. We will say more about smoothing and
sampling in section {\sc I}.2.

\begin{figure}[!ht]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig00.01.0.eps}
}
\caption[A noisy image]{\label{biga} A noisy image magnified to show the pixels.}
\end{figure}

It is widely assumed that the underlying ``real image''  $u_S$ is
either a measure or, for more optimistic authors, a function that
has strong discontinuities. Rudin in 1987 \cite{Rudin:1987:INA}
and De~Giorgi and Ambrosio in 1988 \cite{DeGiAm} proposed
independently the space $BV(\mathbb R^2)$ of functions with
bounded variation as the correct function space for modeling the
images $u_S$. A function $f$ is in $BV(\mathbb R^2)$ if its
partial derivatives ${\partial f}/{\partial x}$ and ${\partial
f}/{\partial y}$, taken as distributions, are Radon measures with
finite total mass. $BV(\mathbb R^2)$ looked at first well adapted
to modeling digital images because it contains functions having
step discontinuities. In fact, the characteristic functions of
smooth domains in $\mathbb R^2$ belong to $BV(\mathbb R^2)$.
However, in 1999, Alvarez, Gousseau, and Morel used a statistical
device on digital images $u_d$ to estimate how the corresponding
images $u_S$ oscillate \cite{Alvarez:1999:SON}. They deduced by
geometric-measure arguments, that the $u_S$ have, in fact,
unbounded variation. We may therefore accept the idea that these
high-resolution images contain very strong oscillations. Although
the images $u_f$ are smoothed versions of the $u_S$, and hence the
oscillations have been averaged, common sense tells us that they
also have large derivatives at transitions between different
observed objects, that is, on the apparent contours of physical
objects. Furthermore, we expect that these large derivatives
(along with noise) are passed to the digital images $u_d$.

\section{Image processing}

For the convenience of exposition, we divide  image processing
into separate disciplines. These are distinguished not so much by
their techniques, which often overlap, as they are by their goals.
We will briefly describe two of these areas: compression and
restoration. The third area, image analysis, is the main subject
of the book and will be discussed in more detail.

\subsection*{Image compression}
 Compression is based on the
discrete nature of digital images, and it is motivated by economic
necessity: Each form of storage and transmission has an associated
cost, and hence one wishes to represent an image with the least
number of bits that is compatible with end usage. There are two
kinds of compression: lossless compression and lossy compression.
Lossless compression algorithms are used to compress digital files
where the decompressed file must agree bit-by-bit with the
original file. Perhaps the best known example of lossless
compression is the zip format. Lossless algorithms can be used on
any digital file, including digital images. These algorithms take
advantage of the structure of the file itself and have nothing to
do with what the file represents. On the other hand, lossy
compression algorithms take advantage of redundancies in natural
images and subtleties of the human visual system. Done correctly,
one can throw away information contained in an image without
impairing its usefulness. The goal is to develop algorithms that
provide high compression factors without objectionable visible
alterations. Naturally, what is visually objectionable depends on
how the decompressed image is used. This is nicely illustrated
with our photographic example. Suppose that we capture the image
$u_f$ at our camera's highest resolution. If we are going to send
$u_d$ over the Internet to a publisher to be printed in a
high-quality publication, then we want no loss of information and
will probably send the entire file in the zip format. If, however,
we just want the publisher to have a quick look at the image, then
we would probably send $u_d$ compressed as a .jpg file, using the
Joint Photographic Expert Group (JPEG) standard for still image
compression. This kind of compression is illustrated in Figure
\ref{IMGPREFACE1}.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.01.1.eps}
\includegraphics{fig00.01.2.eps}
\includegraphics{fig00.01.3.eps}
\includegraphics{fig00.01.4.eps}
} \caption[Compression]{\label{IMGPREFACE1} Compression. Left to
right: the original image and its increasingly  compressed
versions. The compression factors are roughly 7, 10, and 25. \add{
Up too a 10 factor, alterations are hardly visible}.}
\end{figure}

\subsection*{ Image restoration}

A second area is restoration or denoising.  Restoring digital images is much like restoring dirty or damaged paintings
or photographs. Beginning with a digital image that contains blurs or other perturbations (all of which may be
considered as noise), one wishes to produce a better version of the image; one wishes to enhance aspects of the image
that have been attenuated or degraded. Image restoration plays an important role in law enforcement and legal
proceedings. For example, surveillance cameras generally produce rather poor images that must often be denoised and
enhanced as needed. Image restoration is also important in science. When the Hubble Space Telescope was first launched
in 1990, and until it was repaired in 1993, the images it returned were all blurred due to a spherical aberration in
the telescope's primary mirror. Elaborate (and costly) algorithms were developed to restore these poor images, and
indeed useful images were obtained during this period. Restoration is illustrated in Figure \ref{IMGPREFACE2} with an
artificial example. The image on the left has been ostensibly destroyed by introducing random-valued pixels amounting
to 75\% of the total pixel count. Nevertheless, the image can be significantly restored, and a restored version is
shown on the right, by using a Vincent and Serra operator which we will study in Chapter \ref{nonflat}, the  ``area
opening".


\begin{figure}[!ht]
\setkeys{Gin}{width=6cm} \centerline{
\includegraphics{fig00.02.1.eps}
\includegraphics{fig00.02.2.eps}
} \caption[Denoising]{\label{IMGPREFACE2} Denoising.  Left: an
image with up to 75\% of its pixels contaminated by simulated
noise. Right: a denoised version by the Vincent-Serra algorithm
(area opening). \comment{I doubled the size of these images}}
\end{figure}

\subsection*{Image analysis} A third area of image
processing is low-level image analysis, and since this is the main topic of the book, it is important to explain what
we mean by ``low-level'' and ``analysis.'' ``Analysis'' is widely used in mathematics, with various shades of meaning.
Our use of ``analyze,'' and thus of ``analysis,'' is very close to its common meaning, which is to decompose a whole
into its constituent parts, to study the parts, and to study their relation to the whole. For our purposes, the
constituent parts are, for the most part, the ``edges'' and ``shapes'' in an image. These objects, which are often
called features, are things that we could, for a given image, point to and outline, although for a complex natural
image this would be a tedious process. The goal of image analysis is to create algorithms that do this automatically.

The term ``low-level'' comes from the study of human vision and means extracting reliable, local geometric information
from an image. At the same time, we would like the information to be minimal but rich enough to characterize the image.
The goal here is not compression, although some of the techniques may provide a compressed representation of the image.
Our goal is rather to answer questions like, Does a feature extracted from image A exist in image B?\ We are also
interested in comparing features extracted from an image with features stored in a database. As an example, consider
the level set at the left of Figure \ref{IMGPREFACE3}. It consists of major features (roughly, the seven appendages)
and noise. The noise, which is highly variable, prevents us from comparing the image directly with other images having
similar shapes. Thus we ask for a sketchy version, where, however, all essential features are kept. The images on the
right are such a sketchy versions, where most of the spurious details (or noise) have disappeared, but the main
structures are maintained. These sketchy versions may lead to concise invariant encoding of the shape. Notice how the
number of inflexion points of the shape has decreased in the simplification process. This is an example of what we mean
by image analysis. The aim is not denoising or compression. The aim is to construct an invariant code that puts in
evidence the ``main parts'' of an image (in this case, the appendages) and that facilitates fast recognition in a large
database of shapes.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.03.1.eps}
\includegraphics{fig00.03.2.eps}
\includegraphics{fig00.03.3.eps}
\includegraphics{fig00.03.4.eps}
}
\caption[Analysis of a shape]{\label{IMGPREFACE3} Analysis of a shape. The original scanned shape is on the left. Simplified versions are to the right.  }
\end{figure}

\subsection*{Edge detection and scale space}

Since the earliest work in the 1960s,  one of the goals of image
analysis has been to locate the strong discontinuities in an
image. This search is called {\it edge detection}, and it derives
from early research that involved working with images of cubes.
This seemingly simple goal turned out to be exceedingly difficult.
Here is what David Marr wrote about the problem in the early 1980s
(\cite{Marr:1982:V}, p. 16):

\begin{quotation}
The first great revelation was that  the problems are difficult.
Of course, these days this fact is a commonplace. But in the 1960s
almost no one realized that machine vision was difficult. The
field had to go through the same experience as the machine
translation field did in its fiascoes of the 1950s before it was
at last realized that here were some problems that had to be taken
seriously. The reason for this misconception is that we humans are
ourselves so good at vision. The notion of a feature detector was
well established by Barlow and by Hubel and Wiesel, and the idea
that extracting edges and lines from images might be at all
difficult simply did not occur to those who had not tried to do
it. It turned out to be an elusive problem: Edges that are of
critical importance from a three-dimensional point of view often
cannot be found at all by looking at the intensity changes in an
image. Any kind of textured image gives a multitude of noisy edge
segments; variations in reflectance and illumination cause no end
of trouble; and even if an edge has a clear existence at one
point, it is as likely as not to fade out quite soon, appearing
only in patches along its length in the image. The common and
almost despairing feeling of the early investigators like B.K.P.
Horn and T.O. Binford was that practically anything could happen
in an image and furthermore that practically everything did.
\end{quotation}

The point we wish to emphasize is that  textures and noise (which
are often lumped together in image analysis) produce unwanted
edges. The challenge was to separate the ``true edges'' from the
noise. For example, one did not want to extract all of the small
edges in a textured wall paper; one wanted the outline of the
wall. The response was to blur out the textures and noise in a way
that left the ``true edges'' intact, and then to extract these
features. More formally, image analysis was reformulated as two
processes: smoothing followed by edge detection. At the same time,
a new doctrine, the {\it scale space}, was proposed. Scale space
means that instead of speaking of features of an image at a given
location, we speak of them at a given location and at a given
scale, where the scale quantifies the amount of smoothing
performed on the image before computing the features. We will see
in experiments that ``edges at scale 4" and ``edges at scale 7"
are different outputs of an edge detector.

\subsection*{Three requirements for image smoothing
operators}

We have advertised that this book is about image  analysis, which
we have just defined to be smoothing followed by edge detection,
or feature extraction. In fact, the text focuses on smoothing and
particularly on discussing and answering the question, What kind
of smoothing should be used? To approach this problem, we need to
introduce three concepts associated with image analysis operators.
These concepts will be used to narrow the field of smoothing
operators. We introduce them informally at first; more precise
meanings will follow.

\subsection*{Localization}


The first notion is {\it localization}. Roughly  speaking, to say
that an operator $T$ is localized means it essentially uses
information from a small neighborhood of $\bm{x}$ to compute the
output $Tu(\bm x)$. Recall that the sampling operator $R$ in the
photographic example was well localized. As another example,
consider the classic Gaussian smoothing operators $\mathcal{G}_t$
defined by
\[
\mathcal{G}_tu(\bm{x})=G_t*u(\bm{x})=\int_{\RR^2} G_t(\bm{y})u(\bm{x}-\bm{y})\,\mathrm{d}\bm{y},
\]
where $G_t(\bm{x})=({1}/{4\pi t})\mathrm{e}^{-|\bm{x}|^2/4t}$. If $t>0$ is small, then the Gaussian $G_t$ is well localized around zero and $\mathcal{G}_tu(\bm{x})$ is essentially an average of the values of $u(\bm{x})$ in a small neighborhood of $\bm{x}$. The importance of localization is related to the occlusion problem: Most optical images consist of a superposition of different objects that partially obscure one another. It is clear that we must avoid confusing them in the analysis, as would, for example, $G_t$ if $t$ is large. It is for reasons like this that we want the analysis to be as local as possible.

We will prove in Chapter 1 under rather general conditions that $u(t,\bm{x})=G_t*u_0(\bm{x})$ is the unique solution of the heat equation
\[
\frac{\partial{u}}{\partial t} =\Delta u
\]
with initial  value $u_0$. Thus, we can say that smoothing $u_0$
with the Gaussian $G_t$ is equivalent to applying the heat
equation to $u_0$. We will see that the heat equation is possibly
the worst candidate in our search for the ideal smoothing
operator, since, except for small $t$, it is poorly localized and
produces a very blurred image.

\subsection*{Iteration}
 One might conjecture that a way
around this problem with the heat equation would be to replace $G_t$ with a more suitable positive kernel. This is not
the case, but it does serve to introduce the second concept, which is {\it iteration}. We will show in Chapter 2 that
under reasonable assumptions and appropriate rescalings, iterating a convolution with a positive kernel leads to the
Gaussian, and thus directly back to the heat equation. There is, however, a different point of view that leads to
useful smoothing operators: Instead of looking for a different kernel, look for other PDEs that provide smoothing. This
program leads to a class of nonlinear PDEs, where the Laplacian in the heat equation is replaced by various nonlinear
operators. We will see that for these operators it is generally better, from the localization point of view, to iterate
a well localized operator than to apply it directly at a large scale. This, of course, is just not true for the heat
equation; if you iterate $n$ times the convolution $G_t*u$  you get exactly $G_{nt}*u$. This is a good place to point
out that if we are dealing with smoothing, localization, and iteration, then we are talking about parabolic PDEs. This
announcement is heuristic, and the object of the book is to formalize and to make precise the necessity and the role of
several PDEs in image analysis.

\subsection*{Invariance }
 Our last concept is {\it invariance}. Invariance requirements play a central role in image
analysis because the objects to be recognized must be recognized under varying conditions of illumination (contrast
invariance) and from different points of view (projective invariance). Contrast invariance is one of the central
requirements of the theory of image analysis called mathematical morphology (see, for example, Matheron
\cite{Matheron:1975:RSI} or Serra \cite{Serra:1982:IAM}). This theory involves a number of contrast-invariant image
analysis operators, including dilations, erosions, median filters, openings, and closings. We are going to use this
theory by attempting to localize as much as possible these {\it morphomath} operators to exploit their behavior at
small scales. We will then iterate these operators. This will lead to the proof that several geometric PDEs, namely,
the curvature motions, are asymptotically related to certain   morphomath operators in much the same way that linear
smoothing is related to the heat equation. Thus, through these PDEs, one is able to combine the scale space doctrine
and mathematical morphology. In particular, affine-invariant morphomath operators, which seemed at first to be
computationally impractical, turn out to yield in their local iterated version a very affordable PDE, the so called
{\it affine morphological scale space} (AMSS) equation.

\subsection*{Shannon's sampling theory }
 We mentioned in
section {\sc I}.1 that most of the digital images $u_d$ that come to us in practice have been sampled from a smoothed
version, call it $u_f$, of the ``real image'' $u_S$. This was basically a comment about the technology. Another comment
(or assumption) was that the sampling rate was high enough to capture all of the information in $u_f$ that is needed in
practice. What we mean by this is that the representations of $u_f$ that we reconstruct from $u_d$ show no signs that
$u_f$ was undersampled. This is an empirical statement; we will comment on the theory in a moment, but first we wish to
illustrate in Figure \ref{IMGPREFACE4} what can happen if an image is undersampled.

We call the original image on the left Victor. Notice that
Victor's sweater contains a striped pattern, which has a spatial
frequency that is high relative to other aspects of the picture.
If we attempt to reduce the size of Victor by simply sampling, for
example, by taking one pixel in sixteen in a square pattern, we
obtain a new image (the third panel) in which the sampling has
created new and unstable patterns. Notice how new stripes have
been created with a frequency and direction that has nothing to do
with the original. This is called {\it aliasing}, and it is caused
by high spatial frequencies being projected onto lower
frequencies, which creates new patterns. If this had been a video
instead of being a still photo, these newly created patterns would
move and flicker in a totally uncontrolled way. This kind of
moving pattern often appears in recent commercial DVDs. They have
simply not been sampled at a high enough rate. The second panel in
Figure \ref{IMGPREFACE4} is a version of Victor that has been
smoothed enough so that we no longer see the stripes in the
sweater. This image is sampled the same way---every fourth pixel
horizontally and vertically---and appears in panel four. It is not
a good image, but there are no longer the kinds of artifacts that
appear in the third image. To  compare the images we have
magnified the sampled versions by a factor of four. This example
also shows that simply subsampling an image is a poor way to
compress it.

\begin{figure}[!t]
\setkeys{Gin}{width=3cm} \centerline{
\includegraphics{fig00.04.1.eps}
\includegraphics{fig00.04.2.eps}
\includegraphics{fig00.04.3.eps}
\includegraphics{fig00.04.4.eps}
} \caption[Shannon Theory and Sampling]{\label{IMGPREFACE4}
Shannon theory and sampling.  Left to right: original image;
smoothed image; sampled version of the original image; sampled
version of the smoothed image. \add{This illustrates the famous
Shannon-Nyquist law that an image must be smoothed before sampling
in order to avoid the aliasing artifacts. } }
\end{figure}

This pragmatic discussion and the experiment have their
theoretical counterpart, namely, Shannon's theory of sampling.
Briefly, Shannon's theorem, in the two-dimensional case, states
that for an image to be accurately reconstructed from samples, the
image must be bandlimited, which means that it contains no spatial
frequencies greater than some bound $\lambda$, and the sampling
rate must be higher than a factor of $\lambda$. Some implications
of these statements are that the image $u$ must be infinitely
differentiable, that its domain of definition is all of
$\mathbb{R}^2$, and that there must be an infinite number of
samples to accurately reconstruct $u$. Furthermore, in Shannon's
theory, the image $u$ is reconstructed as an infinite series of
trigonometric functions.  Note that this is very different from
what was done in Figure \ref{IMGPREFACE4}. So what does this have
to do with the problems addressed in this book? What does this
have to do with, say, a hypothesized $u_S$ in $BV(\mathbb{R}^2)$
that is definitely not bandlimited? Our answer, which may smack of
smoke and mirrors, is that we always are working in two parallel
worlds, the theoretical one and the practical one based on
numerical computations, and that these two worlds live together in
harmony at a certain scale. Here is an example of what we mean:
Suppose that $u$ is not a bandlimited image. To sample it properly
we would first have to smooth it with a bandlimited kernel.
Suppose that instead we smooth it with the Gaussian $G_t$, which
is not bandlimited. Theoretically this is wrong, but practically,
the spectrum of $G_t$, which is $G_t$ itself, decays
exponentially. If $|\bm{x}|^2/4t$ is sufficiently large, then
$G_t(\bm{x})$ appears as zero in computations, and thus it is
``essentially'' bandlimited. Arguments like this could be made for
other situations, but the important point for the reader to keep
in mind is that the parallel developments, theory and practice,
make sense in the limit.

In the next section, we present a survey of most of the PDEs that have been proposed for image analysis. This provides an informal account of the mathematics that will be developed in detail in the following chapters.

We wish to end this section with a mild disclaimer, and for this we take a page from {\it Theory of Games and Economic Behavior} by John von Neumann and Oskar Morgenstern where they comment on their theory of a zero-sum two-person game \cite{JvonN47} p.~147:
\begin{quotation}
We are trying to find a satisfactory theory,---at this stage for the zero-sum two-person game. Consequently we are not arguing deductively from the firm basis of an existing theory---which has already stood all reasonable tests---but we are searching for such a theory.\ldots This consists in imagining that we have a satisfactory theory of a certain desired type, trying to picture the consequences of this imaginary intellectual situation, and then drawing conclusions from this as to what the hypothetical theory must be like in detail. If this process is applied successfully, it may narrow the possibilities for the hypothetical theory of the type in question to such an extent that only one possibility is left,---i.e. that the theory is determined, discovered by this device. Of course, it can happen that the application is even more ``successful,'' and that it narrows the possibilities down to nothing---i.e. that it demonstrates that a consistent theory of the kind desired is inconceivable.
\end{quotation}

We take much the same philosophical position, and here is our variation on the von Neumann--Morgenstern statement: We do not suggest that what will be developed here is a necessary future for image analysis. However, 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.

\section{PDEs and image processing}

We have argued that smoothing---suppressing high spatial frequencies---is a necessary part of image processing in at
least two situations: An image needs to be smoothed before features can be extracted, and  images must be smoothed
before they are sampled. We have also mentioned that, while smoothing with the Gaussian is not a good candidate for the
first situation (we will see that it is not contrast invariant, and it is not well localized except for small $t$), it
is not unreasonable to use it numerically in the second situation, since it does a good job of suppressing high
frequencies. These smoothing requirements and the fact that the Gaussian is the fundamental solution of the heat
equation mean that the heat equation appears completely naturally in image processing,  and indeed it is the first PDE
to enter the picture in Chapters 1 and 2. Smoothing with the heat equation is illustrated in Figure \ref{IMGPREFACE5}.
\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.05.1.eps}
\includegraphics{fig00.05.2.eps}
} \caption[Heat equation and smoothing]{\label{IMGPREFACE5} Heat
equation and smoothing. The original image is on the left; the
heat equation has been applied at some scale, and the resulting
blurred image is on the right.}
\end{figure}

There is another path hinted at in section {\sc I}.1 that leads to
the Gaussian and thus to the heat equation. Suppose that $k$ is
any positive kernel such that $k(\bm{x})=k(|\bm{x}|)$ and such
that $k$ is localized in the sense that $k(\bm{x})\to 0$
sufficiently rapidly as $|\bm{x}|\to\infty$. If $k$ is normalized
properly and if we write $k_h(\fxx)=(1/h)k({\fxx}/{h^{1/2}})$,
then
\[
\frac {k_h\ast u_0(\fxx)  -u_0(\fxx)}h  \to\Delta u_0(\fxx)
\]
as $h\to 0$ whenever the image $u_0$ is sufficiently smooth. We write this as
\begin{equation}\label{r1}
k_h\ast u_0(\fxx)  -u_0(\fxx)=h\Delta u_0(\fxx)+o(h).
\end{equation}
Now let $u(t,\bm{x})$ denote the solution of the heat equation
\[
\frac{\partial u}{\partial t}=\Delta u, \quad u(0,\bm{x})=u_0(\bm{x}).
\]
If $u_0$ is sufficiently smooth, then we can write
\begin{equation}\label{r2}
u(t,\bm{x})-u(0,\bm{x})=t\Delta u_0(\fxx)+o(t).
\end{equation}

\subsection*{The reverse heat equation } Equations
\eqref{r1} and \eqref{r2} suggest that blurring $u_0$ with a kernel $k_h$ for small $h$ is equivalent to applying the
heat equation to $u_0$ at some small scale $t$. This is true, and it will be made precise in Chapter 2. These equations
also lead to another idea: We read in the paper \cite{Lindenbaum:1994:GCI} by Lindenbaum, Fischer, and Bruckstein that
Kovasznay and Joseph \cite{Kovasznay:1955} introduced in 1955 the notion that a slightly blurred image could be
deblurred by subtracting a small amount of its Laplacian. Numerically, this amounts to subtracting a fraction $\lambda$
of the Laplacian of the observed image from itself:
\[
u_{\mathrm{restored}}=u_{\mathrm{observed}}- \lambda\Delta u_{\mathrm{observed}}.
\]

Dennis Gabor, who received the Nobel prize in 1971 for his
invention of optical holography, studied this process and
determined that the best value of $\lambda$ was the one that
doubled the steepest slope in the image
\cite{Lindenbaum:1994:GCI}. Empirically, one can start with a
small value of $\lambda$ and repeat the process until a good image
is obtained; with further repetitions the process blows up.
Indeed, this process is just applying the reverse heat equation to
the observed image, and the reverse heat  equation is notoriously
ill-posed. On the other hand, the Kovasznay--Joseph--Gabor method
is efficient for sufficiently small $\lambda$ and can be
successfully applied to most images obtained from optical devices.
This process is illustrated in Figure \ref{IMGPREFACE6}. A few
iterations can enhance the image (second panel), but the inverse
heat equation finally blows up (third panel).

\begin{figure}[!ht]
\setkeys{Gin}{width=4cm} \centerline{
\includegraphics{fig00.06.1.eps}
\includegraphics{fig00.06.2.eps}
\includegraphics{fig00.06.3.eps}
} \caption[Kovasznay--Joseph--Gabor  deblurring
I]{\label{IMGPREFACE6} Kovasznay--Joseph--Gabor deblurring.
\comment{I  enlarge these three images (and the next three ones)
as much as possible in order to make clear the improvement in the
second one} Left to right: original image; three iterations of the
algorithm; ten iterations of the algorithm.\label{vicdeblur}}
\end{figure}

Figure \ref{IMGPREFACE7} shows that same experiment applied to an image of Victor that has been numerically blurred. Again, the process blows up, but it yields a significant improvement at some scales.

\begin{figure}[!ht]
\setkeys{Gin}{width=4cm} \centerline{
\includegraphics{fig00.07.1.eps}
\includegraphics{fig00.07.2.eps}
\includegraphics{fig00.07.3.eps}
}
\caption[Kovasznay--Joseph--Gabor deblurring II]{\label{IMGPREFACE7} Kovasznay--Joseph--Gabor deblurring. This is the same deblurring experiment as in Figure \ref{vicdeblur}, but it is applied to a much more blurred image.}
\end{figure}


We have now seen the heat equation used in two senses, each with a different objective. In both cases, we have noted drawbacks. In the first instance, the heat equation (or Gaussian) was used to smooth an image, but as we have mentioned, this operator is not contrast invariant, and thus is not appropriate for any theory of image analysis that requires contrast-invariant operators. This does not mean that the Gaussian should be dismissed; it only means that it is not appropriate for our version of image analysis. To meet our objectives, we will replace the Laplacian, which is a linear isotropic operator, with nonlinear, nonisotropic smoothing operators. This will bring us to the central theme of the book: appropriate smoothing for a possible theory of image analysis.

In the second instance, the heat equation is run backward (the inverse heat equation) with the objective of restoring a blurred image. As we have seen, this is successful to some extent, but the drawback is that it is an unstable process. The practical problem is more complex than the fact that the inverse heat equation is not well posed. In the absence of noise, the best way to deblurr a slightly blurred image is to use the inverse heat equation. However, in the presence of noise, this isotropic operator acts equally in all direction, and while it enhances the definition of edges, the edges become jagged due to the noise. This observation led Gabor to try to improve matters by using more directional operators in place of the Laplacian. Gabor was concerned with image restoration, but his ideas will appear later in our story in connection with smoothing. (For an account of Gabor's work see \cite{Lindenbaum:1994:GCI}.)

\subsection*{Shock filters}

The objective for running the heat equation backward is image restoration, and although restoration is not the main subject of the book, we are going to pause here to describe two ways to improve the stability of the inverse heat equation. Image restoration is an extremely important area of image processing, and the techniques we describe illustrate another use of PDEs in image processing. There are indeed stable ways to ``reverse'' the heat equation. More precisely, there are ``inverse diffusions'' that deblurr an image and reach a steady state. The first example, due to Rudin in 1987 \cite{Rudin:1987:INA} and Osher and Rudin in 1990 \cite{Stanley:1990:FOI} is  a pseudoinverse for the heat equation, where the propagation term $|Du|=|(u_x, u_y)|$ is controlled by the sign of the Laplacian:
\begin{equation}\label{ro}
\frac {\partial u}{\partial t} = -\mathrm{sign}(\Delta u) |Du|.
\end{equation}
This equation is called a {\it shock filter}. We will see later
that this operator propagates the level lines of an image with a
constant speed and in the same direction as the reverse heat
equation would propagate these lines; hence it acts as a
pseudoinverse for the heat equation. This motion enhances the
definition of the contours and thus sharpens the image. Equation
\eqref{ro} is similar to a classic nonlinear filter introduced by
Kramer in the seventies \cite{Kramer:1975:INL}.  Kramer's filter
can be interpreted in terms of a PDE using the same kinds of
heuristic arguments that have been used to derive the heat
equation. This equation is
\begin{equation}\label{kramer}
\frac {\partial u}{\partial t} = -\mathrm{sign}(D^2u(Du, Du)) |Du|,
\end{equation}
where the Laplacian has been replaced by \begin{equation} D^2u(Du,
Du)=u_{xx}(u_x)^2+2u_{xy}u_xu_y+u_{yy}(u_y)^2.
\label{Haralickintro}
\end{equation}
We will see in Chapter 2 that $D^2u(Du, Du)/|Du|^2$ is the second
derivative of $u$ in the direction of its gradient $Du$,  and we
will interpret the differential operator \eqref{Haralickintro} as
Haralick's edge detector. Kramer's equation yields a slightly
better version of a shock filter. The actions of these filters are
illustrated in Figure \ref{IMGPREFACE10}. The image on the left is
a blurred image of Victor. The next image has been deblurred using
the Rudin--Osher shock filter. This is a pseudoinverse of the heat
equation that attains a steady state. The third image has been
deblurred using Kramer's improved shock filter, which also attains
steady state. The fourth image was deblurred using the
Rudin--Osher--Fatemi restoration scheme, which is described below
\cite{Rudin:1992:NTV}.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm} \centerline{
\includegraphics{fig00.08.1.eps}
\includegraphics{fig00.08.2.eps}\\
\includegraphics{fig00.08.3.eps}
\includegraphics{fig00.08.4.eps}
}
\caption[Deblurring with shock filters and a variational method]{
\label{IMGPREFACE10}
Deblurring with shock filters and a variational method.  Left to right: blurred image; Rudin--Osher shock filter; Kramer's improved shock filter; Rudin--Osher--Fatemi restoration method.}
\end{figure}

The deblurring algorithms \eqref{ro} and \eqref{kramer} work to the extent that, experimentally, they attain steady states and do not blow up. However, a third deblurring method, the Rudin--Osher--Fatemi algorithm, is definitely better. It poses the deblurring problem as an inverse problem. It is very efficient when the observed image $u_0$ is of the form $k*u+n$, where $k$ is known and where the statistics of the noise $n$ are also known. Given the observed image $u_0$, one tries to find a restored version $u$ such that $k\ast u$ is as close as possible to $u_0$ and such that the oscillation of $u$ is nonetheless bounded. This is done by finding $u$ that minimizes the functional
\begin{equation}\label{rof}
\int \big(|Du(\bm{x})| + \lambda (k\ast u(\bm{x}) -u_0(\bm{x}))^2\big)\,\rmd\fxx.
\end{equation}
The parameter $\lambda$ controls the oscillation in the restored
version $u$. If $\lambda$ is large, the restored version will
closely satisfy the equation $k\ast u=u_0$, but it may be very
oscillatory. If instead $\lambda $ is small, the solution is
smooth but inaccurate. This parameter can be computed in principle
as a Lagrange multiplier. The obtained restoration can be
remarkable. The best result we can obtain with the blurred Victor
is shown in the fourth panel of Figure \ref{IMGPREFACE10}.  This
scheme was selected by the French Space Agency (CNES) after a
benchmark for satellite image deblurring, and it is currently
being used by the CNES for satellite image restoration. This total
variation restoration method also has fast wavelet packets
versions.

\subsection*{From the heat equation to wavelets } The observation by Kovasznay, Joseph, and Gabor (and undoubtedly others)
that the difference between a smoothed image and the original image is related to the Laplacian of the original image
is also the departure of one of the paths that lead to wavelet theory. Here, very briefly, is the idea: If we convolve
an image with an appropriate smoothing kernel and then take the difference, we obtain a new image related to the
Laplacian of the original image (see equation \eqref{r1}). This new ``Laplacian image'' turns out to be faded with
respect to the original, and if one retains only the values greater than some threshold, the image is often sparse.
This is illustrated in Figure \ref{IMGPREFACE8}. The last panel on the right shows in black the values of this
Laplacian image of Victor that differ significantly from zero. Here, and in most natural images, this representation is
sparse and thus useful for compression. This experiment simulates the first step of a well-known algorithm due to Burt
and Adelson.

In 1983, Burt and Adelson developed a compression algorithm called the Laplacian pyramid based on this idea
\cite{Burt:1983:LPC}. Their algorithm consists of iterating two operations: a convolution followed by subsampling.
After each convolution, one keeps only the difference $k_n*u_n-u_n$, where $n$ is used here to indicate that each step
takes place at a different scale due to the subsampling. The image is then coded by the (finite) sequence of these
differences. These differences resemble the Laplacian of $u_n$, hence the name ``Laplacian pyramid.'' An important
aspect of this algorithm is that the discrete kernels $k_n$, which are low-pass filters, are all the same kernel $k$;
the index $n$ merely indicates that $k$ is adjusted for the scale of the space where the subsampled image $u_n$ lives.
Ironically, the smoothing function cannot be the Gaussian, since the requirements for reconstructing the image from its
coded version rule out the Gaussian. Burt and Adelson's algorithm turned out to be one of the key steps that led to
multiresolution analyses  and wavelets. Burt and Adelson were interested in compression, and, indeed, the differences
$k_n*u_n-u_n$ tend to be sparse for natural images. On the other hand, we are interested in image analysis, and for us,
the Burt and Adelson algorithm has the drawback that it is not translation invariant or isotropic because of the
multiscale subsampling.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.09.1.eps}
\includegraphics{fig00.09.2.eps}
\includegraphics{fig00.09.3.eps}
\includegraphics{fig00.09.4.eps}
}
\caption[The Laplacian pyramid]{\label{IMGPREFACE8} The Laplacian pyramid of Burt and Adelson. Left to right: the original image; the image blurred by Gaussian convolution; the difference between the original image and the blurred version, which approximates the Laplacian of the original image; the points where this Laplacian image is large.}
\end{figure}

\subsection*{Back to edge detection} Early research in computer vision focused on {\it edge detection} as a main tool for
image representation and analysis. It was assumed that the apparent contours of objects, and also the boundaries of the
facets of objects, produce step discontinuities, while inside these boundaries, the image oscillates mildly. The
apparent contour points, or {\it edges points}, were to be computed as points where the gradient is in some sense
largest. Two ways were proposed to do this: Marr and Hildreth proposed computing the points where $\Delta u$ crosses
zero, the now-famous zero-crossings \cite{Marr:1980:TED}. A significant improvement was made by Harakick  who defined
the boundaries, or edges, of an image as those points where $|Du|$ attains a local maximum along the gradient lines
\cite{Har2}. Two years later, Canny implemented Haralick's detector in an algorithm that consists of Gaussian smoothing
followed by computing the (edge) points where $D^2u(Du, Du)=0$ and $|Du|$ is above some threshold
\cite{Canny:1986:CAE}. We refer to this algorithm as the Haralick--Canny edge detector. The fourth panel in Figure
\ref{IMGPREFACE9} displays what happens when we smooth the image with the Gaussian (the heat equation) and then compute
the points where $D^2u(Du,Du)=0$ and $|Du|$ is above some threshold. If this computation is done on the raw image
(first panel), then "edges" show up everywhere (second panel) because the raw image is a highly oscillatory function
and contains a very dense set of inflexion points. After applying the heat equation and letting it evolve to some scale
(third panel), we see that the Haralick--Canny edge detector is able to extract some meaningful structure.

\begin{figure}[!ht]
\setkeys{Gin}{width=6cm}
\centerline{
\includegraphics{fig00.10.1.eps}
\includegraphics{fig00.10.2.eps}
} \caption[\label{heatequationandCanny}Heat equation  and
Haralick's edge detector]{\label{IMGPREFACE9} Heat equation and
Haralick's edge detector.  Left to right: original image; edge
points found in the original image using Haralick's detector;
blurred image; edges found in the blurred image using the
Haralick--Canny detector. \add{The image ``edges" are singled out
after the image has been smoothed. This  smoothing eliminates tiny
oscillations and maintains the big ones. }}
\end{figure}

\subsection*{The Perona-Malik equation }  Given certain
natural requirements such as isotropy, localization, and scale
invariance, the heat equation is the only good linear smoothing
operator. There are, however, many nonlinear ways to smooth an
image. The first one was proposed by Perona and Malik in 1987
\cite{Perona:1987:SSE, Perona:1990:SSE}. Roughly, the idea is to
smooth what needs to be smoothed, namely, the irrelevant
homogeneous regions, and to enhance the boundaries. With this in
mind, the diffusion should look like the heat equation when $|Du|$
is small, but it should act like the inverse heat equation when
$|Du|$ is large. Here is an example of a Perona--Malik equation in
divergence form:
\begin{equation}\label{PM}
\frac {\partial u}{\partial t} = \mathrm{div}(g(|Du|)Du),
\end{equation}
where $g(s) = 1/(1+\lambda^2 s^2)$. It is easily checked that we have a diffusion equation when $\lambda|Du|\leq1$ and an inverse diffusion equation when $\lambda|Du|>1$. To see this, consider the second derivative of $u$ in the direction of $Du$,
\[
u_{\xi\xi}= D^2u\bigg(\frac {Du}{|Du|}, \frac {Du}{|Du|}\bigg),
\]
and the second derivative of $u$ in the orthogonal direction,
\[
u_{\eta\eta}= D^2u\bigg(\frac {Du^\perp}{|Du|}, \frac {Du^\perp}{|Du|}\bigg),
\]
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 then becomes
\[
\frac {\partial u}{\partial t} = \frac{1}{1+\lambda^2|Du|^2}u_{\eta\eta}
+\frac{1-\lambda^2|Du|^2}{(1+\lambda^2|Du|^2)^2}u_{\xi\xi}.
\]
The first term in this representation always appears as a one-dimensional diffusion in the direction orthogonal to the
gradient, tuned by the size of the gradient. The nature of the second term depends on the value of the gradient; it can
be either diffusion in the direction $Du$ or inverse diffusion  in the same direction. This model indeed mixes the heat
equation and the reverse heat equation. Figure \ref{IMGPREFACE11} is used to compare the Perona--Malik equation with
the classical heat equation (illustrated in  Figure \ref{IMGPREFACE9}) in terms of accuracy of the boundaries obtained
by the Haralick--Canny edge detector (see Chapter \ref{ChapLinearScaleSpaceEdgeDetection}). At a comparable scale of
smoothing, we clearly gain some accuracy in the boundaries and remove more ``spurious" boundaries using this
Perona--Malik equation. The representation is both more sparse and more accurate.

\begin{figure}[!ht]
\setkeys{Gin}{width=6cm}
\centerline{
\includegraphics{fig00.11.1.eps}
\includegraphics{fig00.11.2.eps}
}
\caption[A Perona and Malik equation and edge detection]{\label{IMGPREFACE11} A Perona--Malik equation and edge detection. This is the same experiment as in Figure I.11, but here the Perona--Malik equation is used in place of the heat equation. Notice that the edge map looks slightly better in this case.}
\end{figure}

The ambitious Perona--Malik model  attempts to build into a single
operator the ability to perform two very different tasks, namely,
restoration and analysis. This has its cost: The model contains a
``contrast threshold" $\lambda^{-1}$ that must be set manually,
and although experimental results have been impressive, the
mathematical existence and uniqueness of solutions are not
guaranteed, despite some partial results by Kichenassamy
\cite{Kichenassamy:1997} and Weickert and Benhamouda
\cite{WeickertBen}. There are three parameters involved in the
overall smoothing and edge-detecting scheme: the gradient
threshold $\lambda^{-1}$ in the equation \eqref{PM}, the smoothing
scale(s) $t$ (or the time that equation \eqref{PM} evolves), and
the gradient threshold in the Haralick--Canny detector. We can
use the same gradient threshold in both the Haralick--Canny
detector and the Perona--Malik equation, but this still leaves us
with a two-parameter algorithm. Can these parameters be dealt with
automatically for an image analysis scheme? This question seems to
have no general answer at present. An interesting attempt based on
statistical arguments had been made, however, by Black et al.
\cite{Black:1998}.

\subsection*{A proliferation of PDE's\label{proliferation}}
 If one believes that some nonlinear
diffusion might be a good image analysis model, why not try them
all? This is exactly what has happened during the last ten years.
We can claim with some certainty that almost all possible
nonlinear parabolic equations have been proposed. A few of the
proposed models are even systems of PDEs. The common theme in this
proliferation of models is this: Each attempt fixes one intrinsic
diffusion direction and tunes the diffusion using the size of the
gradient or the value of an estimate of the gradient. To keep the
size of this introduction reasonable, we will focus on a few of
the simplest models.

We begin with the Rudin--Osher--Fatemi model \cite{Rudin:1992:NTV}. In this model  the $BV$ norm of $u$, $\int
|Du(\bm{x})|\,{\rm d}\bm{x}$, is one of the terms in the expression \eqref{rof} that is minimized to obtain a restored
image. It is this term that provides the smoothing. The gradient descent for $\int |Du(\bm{x})|\,{\rm d}\bm{x}$
translates into the equation
\[
\frac {\partial u}{\partial t} = \mathrm{div} \bigg(\frac {Du}{|Du|}\bigg)=\frac 1 {|Du|} u_{\eta\eta}.
\]
Written this way, the method appears as a diffusion in the
direction orthogonal to the gradient, tuned by the size of the
gradient.  Andreu et al. proved that this equation is well posed
in the space $BV$ of functions of bounded variation \cite{ABCM1,
ABCM2}. A variant of this model was proposed independently by
Alvarez, Lions, and Morel \cite{Alvarez:1992:ISS}. In this case,
the relevant equation is
\[
\frac {\partial u}{\partial t} =\frac {1}{|k\ast Du|}\,|Du|\, \mathrm{div}\bigg (\frac {Du}{|Du|}\bigg)=\frac 1 {|k\ast Du|} u_{\eta\eta},
\]
and again the diffusion is in the direction $Du^\perp$ orthogonal to the gradient. Note that the rate of diffusion depends on the average value $k\ast Du$ of the gradient in a neighborhood of $\bm{x}$, whereas the direction of diffusion, $Du^\perp(\bm{x})/|Du(\bm{x})|$, depends on the value of $Du(\bm{x})$ at $\bm{x}$. The kernel $k$ is usually the Gaussian. Kimia, Tannenbaum, and Zucker, working in a more general shape-analysis framework, proposed the simplest equation of our list \cite{Kimia:1995:SSD}:
\begin{equation}\label{curveqn}
\frac {\partial u}{\partial t}  =|Du|\;\mathrm{div}\bigg (\frac
{Du}{|Du|}\bigg)= D^2u\bigg(\frac {Du^\perp}{|Du|}, \frac
{Du^\perp}{|Du|}\bigg)=u_{\eta\eta}.
\end{equation}
This equation had been proposed earlier in another context by Sethian as a tool for front-propagation algorithms
\cite{Sethian:1985:CEF}. This equation is a ``pure" diffusion in the direction orthogonal to the gradient. We call this
equation {\it the curvature equation}; this is to distinguish it from other equations that depend on the curvature of
$u$ in some other way. These latter will be called {\it curvature equations}. When we refer to the action of the
equations, we often write {\it curvature motions} or {\it curvature-dependent motions}. (See Chapters \ref{Chapter4}
and \ref{Chap5}.)

The Weickert equation can be viewed  as a variant of the curvature equation \cite{Weickert:1998:ADI}. It uses a
nonlocal estimate of the direction orthogonal to the gradient for the diffusion direction. This direction is computed
as the direction $v$ of the eigenvector corresponding to the smallest eigenvalue of $k\ast(Du\otimes Du)$, where
$(\bm{y}\otimes\bm{y})(\bm{x})=(\bm{x}\cdot\bm{y})\bm{y}$. Note that if the convolution kernel is removed, then this
eigenvector is simply $Du^\perp$. \add{ So the equation writes
\begin{equation}\label{diffusiondansunedirection}
\frac {\partial u}{\partial t}  = u_{\eta\eta},
\end{equation}
where $\eta$ denotes the coordinate in the direction $v$.
 The three models just described can be
interpreted as diffusions in a direction orthogonal to the
gradient (or an estimate of this direction), tuned by the size of
the gradient. They are illustrated in Figure \ref{IMGPREFACE13}.
(The original image is in the first panel of Figure
\ref{IMGPREFACE12}.)

Carmona and Zhong proposed a diffusion  in the direction of the
eigenvector $w$ corresponding to the smallest eigenvalue of $D^2u$
\cite{Carmona:1998:ASR}. So the equation is again
\ref{diffusiondansunedirection}, but this  time $\eta$ denotes the
coordinate in the direction of $w$.} This is illustrated in panel
three of Figure \ref{IMGPREFACE12}. Sochen, Kimmel, and Malladi
propose instead a nondegenerate diffusion associated with a
minimal surface variational formulation \cite{Sochen:1998:GFL}.
Their idea was to make a gradient descent for the area, $\int
\sqrt{1+|Du(\bm{x})|^2}\,{\mathrm d}\bm{x}$, of the graph of $u$.
This leads to the diffusion equation
\[
\frac {\partial u}{\partial t} = \mathrm{div}\bigg (\frac {Du}{\sqrt{1+|Du|^2}}\bigg).
\]

At points where $Du$ is large this equation behaves like  $\frac {\partial u}{\partial t}= \mathrm{div}\bigg (\frac
{Du}{|Du|}\bigg),$ where we retrieve the Rudin-Osher-Fatemi model of Section \ref{proliferation}.  At points where $Du$
is small we have $\frac {\partial u}{\partial t}=\mathrm{div}(Du)$ which is the heat equation. This equation is
illustrated in panel four of Figure \ref{IMGPREFACE12}. Other diffusions have also been considered. For purposes of
interpolation, Caselles, Morel, and Sbert proposed a diffusion that may be interpreted as the strongest possible image
smoothing \cite{Caselles:1998:AAI},
\[
\frac {\partial u}{\partial t} = D^2u(Du, Du)=|Du|^2u_{\xi\xi}.
\]
This equation is not used for preprocessing the image as the others are; rather, it is a way to interpolate between the
level lines of an image with sparse level lines (Figure \ref{IMGPREFACE14}). Among the models mentioned, only the
curvature motion proposed by Kimia, Tannenbaum, and Zucker was specifically introduced as a shape analysis tool. We are
going to explain this, but to do so we must say more about image analysis.



\begin{figure}[!ht]
\setkeys{Gin}{width=6cm} \centering{
\includegraphics{fig00.14.1.eps}
\includegraphics{fig00.14.2.eps}\\
\vspace{0.2cm}
\includegraphics{fig00.14.3.eps}
\includegraphics{fig00.14.4.eps}
} \caption[Diffusion models I]{\label{IMGPREFACE13}  Diffusion models I. Left to right: Osher, Sethian 1988: the
curvature equation; Rudin, Osher, Fatemi 1992: minimization of the image's total variation; Alvarez, Lions, Morel 1992:
nonlocal variant of the preceding; Weickert 1994: nonlocal variant of the curvature equation.  All of these models
diffuse only in the direction orthogonal to the gradient, using a more or less local estimate of this direction.  This
explains why the results of the filters are so similar. \add{However, the Weickert model captures better the texture
direction.} \comment{These images have to be doubled in size and put on two rows. Otherwise, no difference can be seen
between the three PDE's.}}
\end{figure}


\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.12.1.eps}
\includegraphics{fig00.12.2.eps}
\includegraphics{fig00.12.3.eps}
\includegraphics{fig00.12.4.eps}
} \caption[Diffusion models II]{\label{IMGPREFACE12} Diffusion
models II. Left to right: original image; Perona--Malik equation
1987, \add{creating blurry parts separated by sharp edges};
Carmona, Zhong 1998 which actually blurs the  whole image:
diffusion along the least eigenvector of $D^2u$; Sochen, Kimmel,
Malladi 1998: minimization of the image graph area. \add{This last
equation has effects similar to the Perona-Malik model.} }
\end{figure}

\begin{figure}[!ht]
\setkeys{Gin}{width=4cm} \centerline{
\includegraphics{fig00.13.1.eps}
\includegraphics{fig00.13.2.eps}
\includegraphics{fig00.13.3.eps}
} \caption[Diffusion models III]{\label{IMGPREFACE14} Diffusion
models III.   \comment{These three images must be enlarged to see
something and this is easily fixed letting them occupy the whole
width.} Left to right: original image; quantized image (only 10
levels are kept - 3.32 bits/pixel); the quantized image
reinterpolated using the Caselles--Sbert algorithm 1998. They
apply a diffusion on the quantized image with values on the
remaining level lines taken as boundary conditions.}
\end{figure}

\subsection*{Principles of image analysis }  There are
probably as many ways to approach image analysis as there are uses of digital images, and today the range of
applications covers much of human activity. Most scientific and technical activities, including particularly medicine,
and even sound analysis (visual sonograms), involve the perceptual analysis of images. Our goal is to look for
fundamental principles that underlie most of these applications and to develop algorithms that are widely applicable.
From a less lofty point of view, we wish to examine the collection of existing and potential image operators to
determine which among them fit our vision of image analysis. Instead of examining an endless list of partial and
specific requirements, we rely on a mathematical shortcut, well known in mechanics, that consists of stating a short
list of invariance requirements. These invariance requirements will lead to a classification of models and point out
the ones that are the most suitable as image analysis tools. The first invariance requirement is the Wertheimer
principle according to which visual perception (and therefore, we add, image analysis) should be independent of the
image contrast \cite{Wertheimer:1923:ULG}. We formalize this as follows: \vspace{-8pt} \subparagraph{Contrast-invariant
classes.}Two images $u$ and $v$ are said to be (perceptually) equivalent if there is a continuous increasing function
$g$ such that $v=g(u)$. In this case, $u$ and $v$ are said to belong to the same contrast-invariant class.
(``Increasing'' always means ``strictly increasing.'') \vspace{-8pt} \subparagraph{Contrast invariance requirement.}An
image analysis operator $T$ must act directly on the equivalence class. As a consequence, we ask that $T(g(u))=g(Tu)$,
which means that the image analysis operator commutes with contrast changes.

\medskip
The contrast invariance requirement rules out the heat  equation
and all of the models described above except the curvature motion
\eqref{curveqn}. Contrast invariance led Matheron in 1975 to
formulate image analysis as set analysis, namely, the analysis of
the level sets of an image. The {\it upper level set} of an image
$u$ at level $\lambda$ is the set
\[
\levelset_\lambda u= \{\fxx\mid u(\fxx)\geq \lambda\}.
\]
We define in exactly the same way the {\it lower level sets}  by
changing ``$ \geq$'' into ``$\leq$.''  The main point to retain
here is the global invariance of level sets under contrast
changes. if $g$ is a continuous increasing contrast change, then
\[
\levelset_{g(\lambda)}g(u)=\levelset_\lambda u.
\]

According to mathematical morphology, the image analysis doctrine founded by Matheron and Serra, the essential image
shape information is contained in its level sets. It can be proved (Chapter \ref{Chap3}) that an image can be
reconstructed, up to a contrast change, from its set of level sets \cite{Matheron:1975:RSI}. Figure \ref{IMGPREFACE16}
shows an image and one of its level sets.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.15.1.eps}
\includegraphics{fig00.15.2.eps}
} \caption[level set of an image]{\label{IMGPREFACE16} An image
and  one of its level sets. On the right is level set 140 of the
left image. This experiment illustrates Matheron's thesis that the
main shape information is contained in the level sets of an image.
Level sets are contrast invariant.}
\end{figure}

The contrast invariance requirement  leads to powerful and simple denoising operators like the so-called {\it extrema
killer, or area opening}, (Chapter \ref{Chap6}) defined by Vincent in 1993 \cite{Vincent:1993:GAO}. This image operator
simply removes all connected components of upper and lower level sets with areas smaller than some fixed value. This
operator is not a PDE; actually it's much simpler. Its effect is amazingly good for impulse noise, which includes the
local destruction of the image and spots. The action of the extrema killer is illustrated in Figure \ref{IMGPREFACE15}.
The original image is in the first panel. In the third panel, the image has been degraded by adding ``salt and pepper''
noise to 75\% of the pixels. The next panel shows its restoration using the extrema killer set to remove upper and
lower level sets with areas smaller than 80 pixels. The second panel shows the result of the same operator applied to
the original.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.16.1.eps}
\includegraphics{fig00.16.2.eps}
\includegraphics{fig00.16.3.eps}
\includegraphics{fig00.16.4.eps}
} \caption[Extrema killer]{\label{IMGPREFACE15} The  extrema
killer filter. Left to right: original image; extrema killer
applied with area threshold equal 80 pixels; 75\% salt and pepper
noise added to the original image; the same filter applied.}
\end{figure}

\subsection*{Level lines as a complete contrast invariant
representation}  In 1996, Caselles, Coll, and Morel further localized the contrast invariance requirement in image
analysis. They proposed as the main objects of analysis the {\it level lines} of an image, that is, the boundaries of
its level sets \cite{Caselles:1999}. For this program---and the previous one involving level sets---to make sense, the
levels sets and level lines must have certain topological and analytic properties. Level sets and isolevel sets
$\{\fxx\mid u(\fxx)=\lambda\}$, which we would like to be the ``level lines," can be defined for any image (or
function) $u$, but they will not necessarily be useful for image analysis. In particular, we cannot directly define
useful level sets and level lines for a digital image $u_d$. What is needed is a representation of $u_d$ for which
these concepts make sense. But this is not a problem. By the assumptions of section {\sc I}.1, a digital representation
$u_d$ of a natural image $S$ has been obtained by suitably sampling a smooth  version of $S$, call it $u_f$, and a
smooth approximation of $u_f$ is available to us by interpolation. There are, of course, different interpolation
methods to produce smooth representations of $u_d$. One can also obtain a useful discontinuous representation by
considering the extension of $u_d$ that is constant on each pixel. For an interpolation method to be useful, the level
lines should have certain minimal properties: They should be composed of a finite number of rectifiable Jordan curves,
and they should be nested. This means that they do not cross, and thus that they form a tree by inclusion (Section
\ref{levellinestructure}.)

A study by Kronrod in 1950 shows that if the function $u$ is continuous, then the isolevels sets $\{\fxx\mid
u(\fxx)=\lambda\}$ are nested and thus form a tree when ordered by inclusion \cite{Kronrod}. These isolevel sets are
not necessarily curves; they are curves, however, if $u$ has continuous first derivatives. Monasse proved Kronrod's
result for lower semicontinuous and upper semicontinuous functions in 2000 \cite{MonassePhD} (see also \cite{BaCaMo}).
His result implies that the extension of $u_d$ that is constant on each pixel yields a nested set of Jordan curves
bounding the pixels. Thus we have at least two ways to associate a set of nested Jordan curves with a digital image
$u_d$, depending on how $u_d$ is interpolated. Given an interpolation method, we call this set of nested curves a {\it
topographic map} of the image.\footnote{The use of level lines is also consistent with the ``$BV$ assumption" mentioned
in section {\sf I}.1, according to which  the correct function space for modeling images is the space $BV$ of functions
of bounded variation. In this case, the coarea formula can be used to associate a set of Jordan curves with an image
(see \cite{AMCAMAMO}) It is, however, in general false for $BV$ functions that the boundaries of lower and upper level
sets  form a nested set of curves; these curves may cross \add{(see again \cite{MonassePhD}.)}} By introducing the
topographic map, the search for image smoothing, which had already been reduced to set smoothing, is further reduced to
curve smoothing. Of course, we require that this smoothing preserves curve inclusion. Level lines of an image at a
fixed level are shown in Figure \ref{IMGPREFACE17}.

\begin{figure}[!ht]
\setkeys{Gin}{width=6cm} \centerline{
\includegraphics{fig00.17.1.eps}
\includegraphics{fig00.17.2.eps}
} \caption[Level lines of an image]{\label{IMGPREFACE17}\comment{I
doubled the size of these images because  otherwise the structure
of level lines as curves is not visible} Level lines of an image.
Level lines, defined as the boundaries of level sets, can be
defined to be a nested set of Jordan curves. They provide a
contrast-invariant representation of the image. On the right are
the level lines at level 183 of the left image.}
\end{figure}

\subsection*{Contrast invariant PDE's}
 Chen, Giga, and Goto \cite{Chen:1989:UEV,
Chen:1991:UEV} and Alvarez et al. \cite{Alvarez:1993:AFE} proved
that if one adds contrast invariance to the  usual isotropic
invariance requirement for image processing, then all multiscale
image analyses should have a curvature-dependent motion of the
form
\begin{equation}\label{curve}
 \frac {\partial u}{\partial t} = F(\mathrm{curv}(u), t) |Du|,
\end{equation}
where $F$ is increasing with respect to its first argument
\add{(see chapters \ref{Chap17} and \ref{Chap18})}. This equation
can be interpreted as follows: Consider a point $\fxx$ on a given
level curve $C$ of $u$ at time $t$. Let $n(\fxx)$ denote the unit
vector normal to $C$ at $\fxx$ and let $\mathrm{curv}(\fxx)$
denote its curvature. Then the preceding equation is associated
with the curve motion equation
\[
\frac {\partial \fxx}{\partial t} = F(|{\bm\kappa|}(\fxx), t) n(\fxx)
\]
that describes how the point $\fxx$ moves in the direction of the normal. The formula defining $curv(u)$ at a point
 $\fxx$ is (Chapter \ref{Chapter4})
$$curv(u)(\fxx)=\frac{1}{|Du|^3}D^2u(Du^{\perp},Du^{\perp})(\fxx)=
\frac{u_{xx}u_y^2-2u_{xy}u_xu_y+u_{yy}u_x^2}{(u_x^2+u_y^2)^{3/2}}(\fxx).$$ The curvature vector at a point of a $C^2$
curve is the second derivative for a curve $\fxx(s)$ parameterized by length : ${\bm\kappa}=\rmd^2\fxx/\rmd s^2$. We
refer to Chapter \ref{Chapter4} for the detailed definitions and the links between the curvature vector of a level line
of $u$ and $curv(u)$.
 Not much more can be said at this level of generality about $F$. Two specific cases play prominent roles in this
subject. The first case is $F(\mathrm{curv}(u), t)=\mathrm{curv}(u)$, the curvature equation \eqref{curveqn}. The
second case is $F(\mathrm{curv}(u), t) = (\mathrm{curv}(u))^{1/3}$.

This particular one-third power form for the curvature dependence provides an important additional invariance, namely,
affine invariance. We would like to have complete projective invariance, but a theorem proved by  Alvarez et al. shows
that this is impossible \cite{Alvarez:1993:AFE}  \add{(Chapter \ref{Chap18})}. The best we can have is invariance with
respect to the so-called Chinese perspective, which preserves parallelism. Most of these equations, particularly when
$F$ is a power of the curvature, have a viscosity solution in the sense of Crandall and Lions \cite{Crandall:1992:UGV}.
This was shown in 1995 by Ishii and Souganidis \cite{Ishii:1995:GMN}. We refer to Chapters \ref{Chapviscositysolutions}
and \ref{CurvatureequationsandIterated} for all details.

As we have mentioned, contrast-invariant processing can be reduced to level set processing and, finally, to level curve processing. The equations mentioned above are indeed equivalent to curve evolution models if existence and regularity have been established. These results exist for the most important cases, namely, for $F(\mathrm{curv}(u), t)=\mathrm{curv}(u)$, called {\it curve shortening}, and for $F(\mathrm{curv}(u), t) = (\mathrm{curv}(u))^{1/3}$, known as {\it affine shortening}. Grayson proved existence, uniqueness, and analyticity for the curve shortening equation \cite{Grayson:1987:HESb},
\begin{equation}\label{curveshort}
\frac {\partial \fxx}{\partial t} =\mathrm{curv}(\fxx)n(\fxx),
\end{equation}

\medskip
NE PAS LAISSER COMMME C'EST : $curv$ n'est pas la meme notation
qu'apres et n'est pas meme defini!

\medskip
and Angenent, Sapiro, and Tannenbaum proved the same results for the affine shortening equation \cite{Angenent:1998:AHF},
\begin{equation}\label{affineshort}
\frac {\partial \fxx}{\partial t} = (\mathrm{curv}(\fxx))^\frac13 n (\fxx).
\end{equation}
These results are very important for image analysis because they ensure that the shortening processes do indeed reduce a curve to a more and more sketchy version of itself.


\subsection*{Affine invariance}
 An experimental
verification of affine invariance for affine shortening is illustrated in Figure \ref{IMGPREFACE21}. The numerical
tests were made using a very fast numerical scheme for the affine shortening designed by Lionel Moisan
\cite{Moisan:1998:APC}.  The principle of this algorithm is explained in Chapter \ref{ChapAffineInvariantMathMorpho}.
Unlike many numerical schemes, this one is itself affine invariant. Each of the three panels in Figure
\ref{IMGPREFACE21} contains three shapes. The first panel shows the action of  an affine transformation $A$: Call the
first shape in the first panel $X$; then the second shape is $A(X)$ and the third shape is $A^{-1}A(X)=X$. The second
panel shows that affine shortening, $S$, commutes with $A$: The shapes are, from left to right, $S(X)$, $SA(X)$, and
$A^{-1}SA(X)$. Since this third shape is the same as the first, we see that $A^{-1}SA(X)=S(X)$, or that $SA(X)=AS(X)$.
The third panel shows the same experiment with affine shortening replaced with curve shortening. Since the first and
third shapes are different, this illustrates that $A$ does not commute with curve shortening, and hence that curve
shortening is not affine invariant.

\begin{figure}[!ht]
\setkeys{Gin}{width=4.2cm}
\centerline{
\includegraphics{fig00.19.1.eps}
\includegraphics{fig00.19.2.eps}
\includegraphics{fig00.19.3.eps}
}
\caption[Experimental verification of the affine invariance of the affine shortening (AMSS)]{\label{IMGPREFACE21} Experimental verification of the affine invariance of the affine shortening (AMSS). The first panel contains three shapes, $X$, $A(X)$, and $A^{-1}A(X)$. The second panel contains $S(X)$, $SA(X)$, and $A^{-1}SA(X)$. The congruence of the first and third shapes implies that $S$ and $A$ commute. In the third panel, the same procedure has been applied using equation \eqref{curveshort}. Here the first and third shapes are not congruent, which shows that the curve shortening is not affine invariant, as expected.}
\end{figure}

Evans and Spruck \cite{Evans:1991:MLS} (also \cite{C:1992:MLSa, C:1992:MLSb, C:1995:MLS}) and Chen, Giga, and Goto \cite{Chen:1989:UEV, Chen:1991:UEV} proved in 1991 that a continuous function moves by the curvature motion (equation \eqref{curve} with $F(\mathrm{curv}(u), t)=\mathrm{curv}(u)$) if and only if almost all of its level curves move by curve shortening (equation \eqref{curveshort}). The same result is true for the affine invariant curve evolution (equation \eqref{curve} with $F(\mathrm{curv}(u), t)=(\mathrm{curv}(u))^{1/3}$) and affine shortening (equation \eqref{affineshort}).

In the case of the curvature motion, this result provides a mathematical justification for the now-classic Osher--Sethian numerical method for moving fronts \cite{Osher:1988:FPCb}:\label{OS} They associate with some curve or surface $C$ its signed distance function $u(\fxx)=\pm d(\fxx,C)$, and  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 classic numerical difference scheme. Thus, the evolution of the curve $C$ is dealt with efficiently and accurately as a by-product of the evolution of $u$. The point of view that we adopt is slightly different from that of Osher and Sethian. We view the image as a generalized distance function to each of its level sets, since we are interested in all of them.

We show in Figure \ref{IMGPREFACE19} how the level lines are simplified by  evolving the image numerically using affine invariant curvature motion. For clarity, we display only sixteen levels of level curves. Notice that the aim here is not subsampling; we keep the same resolution. Nor is the aim restoration; the processed image is clearly worse than the original. The aim is invariant simplification leading to shape recognition.

\begin{figure}[!ht]
\setkeys{Gin}{width=6cm} \centering{
\includegraphics{fig00.18.1.eps}
\includegraphics{fig00.18.2.eps}
\\ \vspace{0.2cm}
\includegraphics{fig00.18.3.eps}
\includegraphics{fig00.18.4.eps}
} \caption[The affine and morphological scale
space]{\label{IMGPREFACE19} \comment{Here again, I doubled the
size because this is one of the  few occasions to get a view of
the nested level lines and  the  effect of AMSS} The affine and
morphological scale space (AMSS model). Left to right: original
image; level lines of this image (16 levels only); original image
smoothed using the AMSS equation; level lines of the third image.}
\end{figure}

Figures \ref{curvaturescalespaceexperiment1} and \ref{curvaturescalespaceexperiment2} illustrate the effect of affine
curvature motion on the values of the curvature of an image. In Figure \ref{curvaturescalespaceexperiment1} the sea
bird image has been smoothed by affine curvature motion at calibrated scale 1. In Figure
\ref{curvaturescalespaceexperiment2} the smoothing is stronger at calibrated scale 4. (A calibrated scale $t$ means
that at this scale a disk with radius $t$ disappears.) The absolute values of the curvature of the smoothed images are
shown in the upper-right panels of both figures, with the convention that the  darkest points have the largest
curvature. For clarity, the curvature is shown only at points where the gradient of the image was larger than 6 \add{
in a scale ranging from 0 to 255}. Note how the density of points having large curvature is reduced in the second
figure where the smoothing is stronger. On the other hand, the regions with large curvature are more concentrated with
stronger smoothing. Each degree of smoothing produces a different curvature map of the original image, and thus
curvature motions can be used as a nonlinear means to compute a "multiscale" curvature of the original image. The
bottom two panels of the figures show, from left to right, the positive curvature and the negative curvature.

\begin{figure}[!ht]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig00.20.1.eps}
} \caption[Curvature scale space I]{
\label{curvaturescalespaceexperiment1}    Curvature scale space I.
Top, left to right: original sea bird image smoothed by affine
curvature motion at calibrated scale 1; the absolute value of the
curvature. Bottom, left to right: the positive part of the
curvature; the negative part. Compare with Figure I.22, where the
calibrated smoothing scale is 4. }
\end{figure}

\begin{figure}[!ht]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig00.21.1.eps}
}
\caption[Curvature scale space II]{\label{curvaturescalespaceexperiment2}Curvature scale space II. Top, left to right: original sea bird image  smoothed by affine curvature motion at calibrated scale 4; the absolute value of the curvature. Bottom, left to right: the positive part of the curvature; the negative part.
Compare with Figure I.21, where the calibrated smoothing scale is 1.
}
\end{figure}

\subsection*{The snake method}  Before proceeding to shape
recognition, we mention that a variant of the curvature equation
can be used for shape detection. This is a well-known method of
contour detection, initially proposed by Kass, Witkin, and
Terzopoulos \cite{Kass:1988:SAC}. Their method was very unstable.
A better method is a variant of curvature motion proposed by
Caselles, Catt\'{e}, Coll, and Dibos \cite{Caselles:1993:GMA} and
improved simultaneously by Caselles, Kimmel, and Sapiro
\cite{Caselles:1997:GAC} and Malladi, Sethian, and Vemuri
\cite{Malladi:1995:SMF}. Here is how it works. The user draws
roughly the desired contour in the image, and the algorithm then
finds the best possible contour in terms of some variational
criterion. This method is  very useful in medical imaging. The
motion of the contour is a tuned curvature motion that tends to
minimize an energy function $E$. Given an original image $u_0$
containing some closed contour  that we wish to approximate, we
start with an edge map
\[
g(\fxx)=\frac 1{1+|Du_0(\fxx)|^2},
\]
that is, a function that vanishes on the edges of the image. The user then designates the contour of interest by drawing a polygon $\gamma_0$ roughly following the desired contour. The {\it 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))\,\mathrm{d}s,
\]
where $g$ is the edge map associated with the original image $u_0$ and $s$ denotes the parameter measuring the length along $\gamma$. The motion of the ``analyzing image" $v$ is governed by
\[
\frac {\partial v}{\partial t}(\bm{x},t) =g(\bm{x}) |Dv(\bm{x})|\mathrm{curv}(v)(\bm{x}) - Dv(\bm{x})\cdot Dg(\bm{x}).
\]
This algorithm is illustrated with a medical example in Figure \ref{IMGPREFACE22}.

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig00.22.1.eps}
\includegraphics{fig00.22.2.eps}
\includegraphics{fig00.22.3.eps}
\includegraphics{fig00.22.4.eps}
}
\caption[Active Contour]{\label{IMGPREFACE22}Active contour, or ``snake.'' Left to right: original image; initial contour; evolved distance function; final contour.}
\end{figure}

\subsection*{Shape retrieval} It seems to us that the most
obvious application of invariant PDEs is shape retrieval in large
databases. There are thousands of different definitions of shapes
and a multitude of shape recognition algorithms. The real
bottleneck has always been the ability to extract the relevant
shapes. The discussion above points to a brute force strategy: All
contrast-invariant local elements, or the level lines of the
image, are candidates to be ``shape elements.''  Of course, this
notion 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 the level lines that are similar
in both images? This would give a factual, a posteriori,
definition of shapes. They would be defined as pieces of level
lines common to two different images, irrespective of their
relationships to real physical objects.

Of course, this brute force strategy would be impossible without
the  initial invariant filtering (AMSS). It is doable only 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. \cite{Lisani:2000:AIM}. Two images of a desk
and the backs of chairs, viewed from different angles, are shown
in the first two panels. All of the pieces of level lines in the
two images that found a match in the other image are shown in the
last two panels. Notice that several of these matches are doubled.
Indeed, there are two similar chairs  in each image. This brings
to mind a Gestalt law that states that human perception tends to
group similar shapes. We now see the numerical necessity of this
perceptual grouping: A preliminary self-matching of each image,
with grouping of similar shapes, must be performed before we can
compare it with other images.

This concludes our overview of the use of PDEs in image analysis. The rest of the book is devoted to filling in the
mathematical details that support most of the results mentioned in this introduction. We have tried to prove all of the
mathematical statements, assuming only two or three years of mathematical training at the university level. Thus, for
most of the PDEs addressed, and for all of the relevant ones, we prove the existence and uniqueness of solutions. We
also develop invariant, monotone approximation schemes. This has been technically possible by combining tools from the
recent, and remarkably simple, theory of viscosity solutions with the Matheron formalism for monotone set and function
operators. Thus, the really necessary mathematical knowledge amounts to elementary differential calculus, linear
algebra, and some results from the theory of Lebesgue  integration, which are used in the chapters on the heat
equation. Mathematical statements are not introduced as art for art's sake; all of the results are directed at proving
the correctness of a model, of its properties, or of the associated numerical schemes. Numerical experiments, with
detailed comments, are described throughout the text. They provide an independent development that is parallel to the
central theoretical development. Most image processing algorithms mentioned in the text are accessible in the public
software MegaWave. MegaWave was developed jointly by several university research groups in France, Spain and America,
and it is available at \texttt{http://www.cmla.ens-cachan.fr}.

\begin{figure}[!t]
\setkeys{Gin}{width=12cm} \centering{
\includegraphics{fig00.23.1.eps}\\
\includegraphics{fig00.23.2.eps}
} \caption[Level-lines based shape parser]{\label{IMGPREFACE20}
\comment{Here again, I doubled the size so that the reader can
appreciate something which is probably  the main application  of
the theory.  I also added a comment on the black parts of level
lines which  was missing and explained better the experiment} A
shape parser based on level lines. The two left images are of a
desk and the backs of chairs viewed from different angles. In the
far left panel, one level line has been selected (in white). In
the second panel we show, also in white, all  matching pieces of
level lines. The match is ambiguous, as must be expected when the
same object is repeated in the scene. In the two panels on the
right, we display all the matching pairs of pieces of level lines
(in white).  The  non matching parts of the same level lines are
shown in black.  Usually, recognized shape elements are
 pieces of level lines, seldom  whole level lines. See \cite{} }
\end{figure}




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