%gmIntro.tex

\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. 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 use 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.

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.}
\end{figure}

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

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\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.}
\end{figure}

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 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. Briefly, and naively, this is the background: The image $u_f$ that arrives at the retina is sampled by the photoreceptors of the retina to produce a digital image $u_d$. The bandwidth of the optic nerve is too small to transmit all of the information in $u_d$ to the higher centers of the brain. To circumvent this problem, it is hypothesized that ``low-level'' or ``early vision'' processing takes place in the retinal area and that the output of this processing is sent to the higher centers of the brain where our perception of the physical world is formed. Thus our use of ``low-level analysis'' 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}

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.

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.

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.

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 the convolution $G_t*u$ $n$ times 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.

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.

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.}
\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 form 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 some 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. 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}

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=3cm}
\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. 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=3cm}
\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}.)

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
\[
D^2u(Du, Du)=u_{xx}(u_x)^2+2u_{xy}u_xu_y+u_{yy}(u_y)^2.
\]
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{kramer} 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 recently won a competition organized by the French Space Agency (CNES) to establish 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.

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 contrast invariant 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}

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.}
\end{figure}

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 diffusion  in the direction $-Du$. 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. 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 attempts 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 wish to 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}.

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, $\int |Du(\bm{x})|\,{\rm d}\bm{x}$, which is the $BV$ norm of $u$, 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 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}.

The Weickert equation can be viewed 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 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$.

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 corresponding to the smallest eigenvalue of $D^2u$ \cite{Carmona:1998:ASR}.  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).
\]
It 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=3cm}
\centerline{
\includegraphics{fig00.14.1.eps}
\includegraphics{fig00.14.2.eps}
\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.}
\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; Carmona, Zhong 1998: diffusion along the least eigenvector of $D^2u$; Sochen, Kimmel, Malladi 1998: minimization of the image graph area.}
\end{figure}

\begin{figure}[!ht]
\setkeys{Gin}{width=3cm}
\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.  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 as boundary conditions.}
\end{figure}


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. (Note that contrast changes $g$ are not limited to increasing functions.)

\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. Thus, 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, all of an image's shape information is contained in its level sets. It can be proved 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} 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 extreme 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}

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.

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.} 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=3cm}
\centerline{
\includegraphics{fig00.17.1.eps}
\includegraphics{fig00.17.2.eps}
}
\caption[Level lines of an image]{\label{IMGPREFACE17} 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}

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. 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(\mathrm{curv}(\fxx), t) n(\fxx)
\]
that describes how the point $\fxx$ moves in the direction of the normal.
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}. 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}.

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



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}. 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=3cm}
\centerline{
\includegraphics{fig00.18.1.eps}
\includegraphics{fig00.18.2.eps}
\includegraphics{fig00.18.3.eps}
\includegraphics{fig00.18.4.eps}
}
\caption[The affine and morphological scale space]{\label{IMGPREFACE19} 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. (As will be explained later, 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. 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}

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}


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 images. 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. No direct industrial applications are treated. However, many of the image processing tools mentioned in the text are currently being used in the public software MegaWave. MegaWave was developed jointly by several university research groups in France and Spain, and it is available at \texttt{http://www.cmla.ens-cachan.fr}. These tools are also extensively used in the commercial video and image processing software Investigator, which is sold by Cognitech, Inc., Pasadena, California.

\begin{figure}[!t]
\setkeys{Gin}{width=6cm}
\centerline{
\includegraphics{fig00.23.1.eps}
\includegraphics{fig00.23.2.eps}
}
\caption[Level-lines based shape parser]{\label{IMGPREFACE20} 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 put in white. In the second panel we show, also in white, the matching 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 level lines in white.}
\end{figure}

\section{Notation and background material}



$\mathbb R^N$ denotes the real $N$-dimensional Euclidian space. If ${\bf x} \in\mathbb R^N$ and $N>2$, we write ${\bm x}=(x_1,x_2,\ldots,x_N)$; if $N=2$, we usually write ${\bm x}=(x,y)$. For ${\bm x},\, {\bm y}\in\mathbb R^N$, we denote their scalar product by ${\bm x}\cdot{\bm y}=x_1y_1+x_2y_2+\cdots+x_Ny_N$ and write $$|{\bm x}|=({\bm x}\cdot{\bm x})^{1/2}=(x_1^2+x_2^2+\cdots+x_N^2)^{1/2}.$$


Let $\Omega$ be an open set in $\mathbb R^N$, and let $n\in\mathbb N$ be a fixed integer. $C^n(\Omega)$ denotes the set of real-valued functions $f:\Omega\to\mathbb R$ that have continuous derivatives of all orders up to and including $n$. $f\in C^{\infty}(\Omega)$ means that $f$ has continuous derivatives of all orders; $f\in C(\Omega)=C^0(\Omega)$ means that $f$ is continuous on $\Omega$. We will often write ``$f$ is $C^n$'' as shorthand for $f\in C^n(\Omega)$, and we often omit the domain $\Omega$ if there is no chance of confusion.

We use multi-indices of the form $\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_N)\in{\mathbb N}^N$ as shorthand in several cases. For ${\bm x}\in{\mathbb R}^N$, we write ${\bm x}^{\alpha}$ and $|{\bm x}|^{\alpha}$ for $x_1^{\alpha_1}x_2^{\alpha_2}\cdots x_N^{\alpha_N}$ and $|x_1|^{\alpha_1}|x_2|^{\alpha_2}\cdots |x_N|^{\alpha_N}$, respectively. For $f\in C^n(\Omega)$, we abbreviate the partial derivatives of $f$ by writing
\[
\partial^\alpha f= \frac{\partial^{|\alpha|}f}{\partial x_1^{\alpha_1}\partial x_2^{\alpha_2}\cdots\partial x_N^{\alpha_N}},
\]
where $|\alpha|=\alpha_1+\alpha_2+\cdots+\alpha_N$ and $|\alpha|\leq n$.

We also write the partial derivatives of $f({\bm x})=f(x_1,x_2,\ldots,x_N)$ as $f_i=\partial f/\partial x_i$, $f_{ij}=\partial^2 f/\partial x_i\partial x_j$, and so on. In the two-dimensional case $f({\bm x})=f(x,y)$, we usually write $\partial f/\partial x=f_x$, $\partial f/\partial y=f_y$, $\partial^2 f/\partial x\partial y =f_{xy}$, and so on.

The gradient of $f$ is denoted by $Df$. Thus, if $f({\bm x})=f(x_1,x_2,\ldots,x_N)$,
\[
Df=(f_1,f_2,\ldots,f_N),
\]
and
\[
Df=(f_x,f_y)
\]
in case $N=2$. The Laplacian of $f$ is denoted by $\Delta f$. Thus $\Delta f=f_{11}+f_{22}+\cdots+f_{NN}$ in general, and $\Delta f=f_{xx}+f_{yy}$ if $N=2$.

We will often use the symbols $O$, $o$, and $\varepsilon$. They are defined as follows. We assume that $h$ is a real variable that tends to a limit $h_0$ that can be finite or infinite. We assume that $g$ is a positive function of $h$ and that $f$ is any other function of $h$. Then $f=O(g)$ means that there is a constant $C>0$ such that $|f(h)|<Cg(h)$ for all values of $h$. The expression $f=o(g)$ means that $f(h)/g(h)\to 0$ as $h\to h_0$. We occasionally will use $\varepsilon$ to denote a function of $h$ that tends to zero as $h\to 0$. Thus, $f(h)=o(h)$ can be written equivalently as $f(h)=h\varepsilon(h)$.

\subsection*{Taylor's formula}

\hspace*{0pt}{\indent} An $N$-dimensional form of Taylor's formula is used several times in the book. We will first state it and then explain the notation. Assume that $f\in C^n(\Omega)$ for some open set $\Omega\in\mathbb R^N$, that ${\bm x}, {\bm y}\in\Omega$, and that the segment joining ${\bm x}$ and ${\bm x}+{\bm y}$ is also in $\Omega$. Then
\[
f({\bm x}+{\bm y})=f({\bm x})+\frac{1}{1!}Df({\bm x}){\bm y}^{(1)}+\frac{1}{2!}D^2f({\bm x}){\bm y}^{(2)}+\cdots+\frac{1}{n!}D^n f({\bm x}){\bm y}^{(n)}+o(|{\bm y}|^n).
\]

This has been written compactly to resemble the one-dimensional case, but the price to be paid is to explain the meaning of $D^p f({\bm x}){\bm y}^{(p)}$. We have already seen special cases of this expression in section {\sc I}.3, for example, $D^2u(Du,Du)$ in equation ({\sc I}.4). The expression $D^p f({\bm x}){\bm y}^{(p)}$ is itself a special case of a more general multilinear form, but this generality is of no interest to us, so we will go directly to the expansion $D^p f({\bm x}){\bm y}^{(p)}$.
\[
D^pf({\bm x}){\bm y}^{(p)}=D^p f({\bm x})(\underbrace{\bm y,\bm y,\ldots,\bm y}_{\text{$p$ terms}})=\sum_{(i_1,i_2,\ldots,i_p)}\frac{\partial^p f}{\partial x_{i_1}\partial x_{i_2}\cdots\partial x_{i_p}}(\fxx)y_{i_1}y_{i_2}\cdots y_{i_p},
\]
where the sum is taken over all $N^p$ different vectors $(i_1,i_2,\ldots,i_p)$, $i_j=1,2,\ldots,N$. Notice that $Df({\bm x}){\bm y}^{(1)}$ is just $\sum_{j=1}^N f_{j}y_j=Df({\bm x})\cdot{\bm y}$, which is how we usually write it.

\subsection*{The implicit function theorem}

\hspace*{0pt}{\indent}Consider a real-valued $C^1$ function $f$ defined on an open set $\Omega$ in $\RR^N$. For ease of notation we write $\fxz=(\fxx,y)$, where $\fxx=(x_1,\ldots,x_{N-1})$ and $y=x_N$. Assume that $f(\fxz_0)=0$ for a point $\fxz_0\in\Omega$ and that $f_y(\fxx_0)\ne0$. Then there is a neighborhood $M=M(\fxx_0)$ and a neighborhood $N=N(y_0)$ such that for every $\fxx\in M$ there is exactly one $y\in N$ such that $f(\fxx,y)=0$. The function $y=\varphi(\fxx)$ is $C^1$ on $M$ and $y_0=\varphi(\fxx_0)$. Furthermore, if $f\in C^n(\Omega)$, then $\varphi\in C^n(M)$.

\subsection*{Lebesgue integration}

\hspace*{0pt}{\indent}The Lebesgue integral, which first appeared in 1901 and is thus over a hundred years old, has become the workhorse of analysis. It plays a central role in chapters 1 and 2 and appears briefly in other parts of the book. One does not need a profound understanding of abstract measure theory and integration to follow the arguments. One should, however, be familiar with a few key results and be comfortable with the basic manipulations of the integral. With this in mind, we restate some of these fundamentals.

The functions and sets in this book are always measurable. Thus we dispense with phrases like ``let $f$ be a measurable function.'' A function $f$ defined on a subset $A$ of $\mathbb R^N$ is integrable, if
\[
\int_A|f({\bm x})|\,\rm d\bm x<+\infty.
\]
The Banach space of all integrable function defined on $A$ is denoted as usual by $L^1(A)$; we write $\|f\|_{L^1(A)}=\int_A|f({\bm x})|\,\rm d\bm x$ to denote the norm of $f$ in $L^1(A)$. The most important applications in the book are the two cases $A=\mathbb R^N$ and $A=[-1,1]^N$. Here are two results that we use in chapters 1 and 2. We state them not in the most  general form, but rather in the simplest form suitable for our work.

\subsection*{A density theorem for $L^1(\mathbb R^N)$}

\hspace*{0pt}{\indent}If $f$ is in $L^1(\mathbb R^N)$, then there exists a sequence of continuous functions $\{g_n\}$, each of which has compact support, such that $g_n\to f$ in  $L^1(\mathbb R^N)$, that is, $\|g_n-f\|\to0$ as $n\to+\infty$. This result is true for $L^1([-1,1]^N)$, in which case the $g_n$ are continuous on $[-1,1]^N$.

\subsection*{Fubini's theorem}

\hspace*{0pt}{\indent}Suppose that $f$ is a measurable function defined on $A\times B\in\mathbb R^N\times\RR^N$. Fubini's theorem states that
\[
\int_{A\times B}|f(\fxz)|\,{\rm d}{\bm z}=\int_{A}\int_{B}|f(\fxx,\fxy)|\,{\rm d}\fxx\,{\rm d}\fxy=\int_{B}\int_{A}|f(\fxx,\fxy)|\,{\rm d}\fxy\,{\rm d}\fxx,
\]
where we have written $\fxz=(\fxx,\fxy)$.
It further states, that if any one of the integrals is finite, then
\[
\int_{A\times B}f(\fxz)\,{\rm d}{\bm z}=\int_{B}\int_{A}f(\fxx,\fxy)\,{\rm d}\fxx\,{\rm d}\fxy=\int_{A}\int_{B}f(\fxx,\fxy)\,{\rm d}\fxy\,{\rm d}\fxx.
\]

\subsection*{Lebesgue's dominated convergence theorem}

\hspace*{0pt}{\indent}If a sequence of functions $\{f_n\}$ is such that $f_n(\fxx)\to f(\fxx)$ for almost every $\fxx\in\RR^N$ as $n\to+\infty$, and if there is an integrable function $g$ such that $|f_n(\fxx)|\leq g(\fxx)$ almost everywhere, then
\[
\int_{\RR^N}f_n(\fxx)\,\rmd\fxx\to\int_{\RR^N}f(\fxx)\,\rmd\fxx.
\]
We do not use this result itself, but we do use the following result, which is a direct consequence of the dominated convergence theorem.

\subsection*{Interchanging differentiation and integration}

\hspace*{0pt}{\indent}Suppose that a function $f$ defined on $(t_0,t_1)\times\RR^N$, where $(t_0,t_1)$ is any interval of $\RR$, is such that $t\mapsto f(t,\fxx)$ is continuously differentiable (for almost every $\fxx\in\RR^N$) on some interval $[a,b]\subset (t_0,t_1)$. If there exists an integrable function $g$ such that for all $t\in[a,b]$
\[
\bigg|\frac{\partial f}{\partial t}(t,\fxx)\bigg|\leq g(\fxx) \quad\text{almost everywhere},
\]
then the integral $I(t)=\int_{\RR^N}f(t,\fxx)\,\rmd\fxx$ is differentiable for $t\in(a,b)$ and
\[
\frac{dI}{dt}(t)=\int_{\RR^N}\frac{\partial f}{\partial t}(t,\fxx)\,\rmd\fxx.
\]
A brief but comprehensive discussion of the Lebesgue integral can be found in the classic textbook by Walter Rudin \cite{Rudin:1966:RCA}.

\vspace{5pt}

\begin{center}
[I wish to add the essential facts about symmetric matrices.]
\end{center}

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