%gmChap02.tex

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

\section{Applications of the Heat Equation to Image Analysis}


\subsection{Convergence of iterated filters to the heat equation}

We wish to prove  that iterated linear smoothing is asymptotic to
applying the heat equation, or equivalently, to convolution with
the Gaussian at different scales.

\begin{definition}
We say that a real-valued kernel $g\in L^1(\RR^N)$ is Laplacian consistent if it satisfies the following moment conditions:
\begin{enumerate}%[$(i)$]\setlength\itemsep{2pt}
\item $\int_{\RR^N} g(\fxx)\,\rmd\fxx=1.$
\item For $i=1,2,\ldots,N$, $\int_{\RR^N} x_ig(\fxx)\,\rmd\fxx=0.$
\item For each pair $i,j=1,2,\ldots,N$, $i\neq j$,
$\int_{\RR^N}x_i x_jg(\fxx)\,\rmd\fxx=0.$
\item For $i=1,2,\ldots,N$, $\int_{\RR^N} x_i^2g(\fxx)\,\rmd\fxx=\sigma,$ where $\sigma>0.$
\item $\int_{\RR^N}|\fxx|^3|g(\fxx)|\,\rmd\fxx<+\infty.$
\end{enumerate}
\end{definition}
Note that we do not assume that $g\geq0$; in fact, many important filters used in signal and image processing are not positive. However, condition $(i)$ implies that $g$ is ``on average'' positive. More precisely, if we consider $g*u(\fxx)$ to be a weighted average of the values of $u(\fxx)$, then condition $(i)$ says that the sum of the weights is one.

We say that a function $g$ is \textit{radial} if $g(\fxx)=g(|\fxx|)$, $\fxx\in\RR^N$. This is equivalent to saying that $g$ is invariant under all rotations around the origin in $\RR^N$.

\begin{exercise}
Assume $g\in L^1(\RR^N)$ is radial with finite second moments, $\int_{\RR^N}|\fxx|^2|g(\fxx)|\,\rmd\fxx<+\infty,$ such that $\int_{\RR^N}\fxx_i^2g(\fxx)\,\rmd\fxx>0.$ Show that $g$ satisfies conditions $(ii)$ and $(iii)$ and that, for suitably chosen $a,b\in\RR$, the rescaled function $\fxx\mapsto ag(\fxx/b)$ satisfies conditions $(i)$ and $(iv)$, where $\sigma$ can be taken to be an arbitrary positive number. \finex
\end{exercise}

Exercise 2.1 says that any radial function $g\in L^1(\RR^N)$ can be rescaled to be Laplacian consistent if it decays fast enough at infinity and if $\int_{\RR^N}\fxx_i^2g(\fxx)\,\rmd\fxx$ and $\int_{\RR^N}g(\fxx)\,\rmd\fxx$ have the same sign. (The significance of this last condition will become clear from the proof of Theorem 2.2.)

We consider rescalings of a kernel $g$ defined by
\begin{equation}\label{gh}
g_h(\fxx)=\frac{1}{h^{{N}/{2}}}\,g\Big(\frac{\fxx}{h^{{1}/{2}}}\Big)
\end{equation}
for $h>0$ (see Figure \ref{fig21}). Notice that this rescaling differs slightly from the one used in section 1.2. We have used the factor $h^{1/2}$ here because it agrees with the factor $t^{1/2}$ in the Gaussian. We denote the convolution of $g$ with itself $n$ times by $g^{n*}$. The main result of this section concerns the behavior of $g_h^{n*}$ as $n\to+\infty$ and $h\to0$.

\begin{exercise}Prove the following two statements:
\begin{enumerate}[(i)]\setlength\itemsep{2pt}
\item $g_h$ is Laplacian consistent if and only it $g$ is Laplacian consistent.
\item If $g\in L^1(\RR^N)$, then $(g_h)^{n*}=(g^{n*})_h$. \finex
\end{enumerate}
\end{exercise}

Our first result concerns the behavior of $g_h$ as $h\to0$. This will establish a more general and precise form of equation \eqref{Mhx}.

\begin{figure}[!t]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.01.1.eps}
}
\caption[Kernel rescaling]{\label{fig21}The rescalings $g_t(\fxx)= ({1}/{t^2}) g({\fxx}/{t})$ of a kernel for t=4, 3, and 2.}
\end{figure}

\begin{theorem} If $g$ is Laplacian consistent, then for every $u\in L^{\infty}_C\cap C^3$,
\begin{equation}\label{LC}
g_h*u(\fxx)-u(\fxx)=h\frac{\sigma}{2}\Delta u(\fxx)+\varepsilon(h,\fxx)
\end{equation}
where $|\varepsilon(h,\fxx)|\leq Ch^{3/2}$.
\end{theorem}
\begin{proof}We use condition $(i)$, the definition of $g_h$, and rescaling inside the integral to see that
\begin{align*}
g_h*u(\fxx)-u(\fxx)=&\int_{\RR^N}\frac{1}{h^{{N}/{2}}}\,g\Big(\frac{\fxx}{h^{{1}/{2}}}\Big)\big(u(\fxx-\fxy)-u(\fxx)\big)\,\rmd\fxy\\
=&\int_{\RR^N}g(\fxz)\big(u(\fxx-h^{1/{2}}\fxz)-u(\fxx)\big)\,\rmd\fxz.
\end{align*}
Using Taylor's formula with the Lagrange remainder, we have
\[
u(\fxx-h^{1/{2}}\fxz)-u(\fxx)=-h^{1/{2}}Du(\fxx)\cdot\fxz+\frac{h}{2}D^2u(\fxx)(\fxz,\fxz)-\frac{1}{6}h^{3/{2}}D^3u(\fxx-h^{1/{2}}\theta\fxz)(\fxz,\fxz,\fxz),
\]
where $\theta=\theta(\fxx, \fxz,h)\in [0,1]$. By condition $(ii)$, $\int_{\RR^N}g(\fxz)Du(\fxx)\cdot\fxz\,\rmd\fxz=0$; by conditions $(iii)$ and $(iv)$, $\int_{\RR^N}g(\fxz)D^2u(\fxx)(\fxz,\fxz)\,\rmd\fxz=\sigma\Delta u(\fxx)$. Thus,
\[
g_h*u(\fxx)-u(\fxx)=h\frac{\sigma}{2}\Delta u(\fxx)-\frac{1}{6}h^{3/{2}}\int_{\RR^N}g(\fxz)D^3u(\fxx-h^{1/{2}}\theta\fxz)(\fxz,\fxz,\fxz)\,\rmd\fxz.
\]
We denote the error term by $\varepsilon(h,\fxx)$. Then we have the following estimate:
\begin{align*}
|\varepsilon(h,\fxx)|\leq&\frac{1}{6}h^{{2}/{3}}\int_{\RR^N}|g(\fxz)D^3u(\fxx-h^{1/{2}}\theta\fxz)(\fxz,\fxz,\fxz)|\,\rmd\fxz\\
\leq&\frac{1}{6}h^{3/{2}}N^{3/{2}}\sup_{\alpha,\fxx}|\partial^{\alpha}u(\fxx)|\int_{\RR^N}|\fxz|^3|g(\fxz)|\,\rmd\fxz,
\end{align*}
where the supremum is taken over all vectors $\alpha=(\alpha_1,\alpha_2,\cdots,\alpha_N)$, $\alpha_j\in\{1,2,3\}$, such that $|\alpha|=3$ and over all $\fxx\in C$. This proves the theorem with
\[
 C=\frac{1}{6}N^{3/{2}}\sup_{\alpha,\fxx}|\partial^{\alpha}u(\fxx)|\int_{\RR^N}|\fxz|^3|g(\fxz)|\,\rmd\fxz.\qedhere
\]
\end{proof}

Theorem 2.2 shows a direct relation between smoothing with a Laplacian-consistent kernel and the heat equation. It also shows why we require $\sigma$ to be positive: If it is not positive, the kernel is associated with the inverse heat equation. The next theorem is the asymptotic result that generalizes the example in section 1.1. This result is illustrated in Figure \ref{22}.

\begin{theorem}Let $g$ be a nonnegative Laplacian-consistent kernel with $\sigma=2$ and define $g_h$ by \eqref{gh}. Write $T_hu_0=g_h*u_0$ for $u_0\in L^1_C$, and let $u(t,\cdot)=G_t*u_0$ be the solution of the heat equation \eqref{heatsol}. Then, for each $t>0$,
\begin{equation}\label{37}\addtolength\belowdisplayskip{-12pt}
T_h^nu_0\to u(t,\cdot)\,\text{ in }\, L^1(C)\,\text{ as }\, n\to+\infty\,\text{ and }\, nh\to t.
\end{equation}
\end{theorem}

\begin{proof}The first order of business is to say precisely what is meant by the asymptotic limit \eqref{37}: Given $t>0$ and given $\varepsilon>0$, there exists an $n_0=n_0(t,\varepsilon)$ and a $\delta=\delta(t,\varepsilon)$ such that $\|T_h^nu_0-u(t,\cdot)\|_{L^1(C)}\leq\varepsilon$ if $n>n_0$ and $|nh-t|\leq\delta$. This is what we must prove. It says that $T_h^nu_0\to u(t,\cdot)$ in $L^1(C)$ as $n\to+\infty$ if, at the same time, $h\to0$ in a certain way, which is that $h$ must be close to the asymptote $t/n$. We will first prove the result by requiring that $h$ is actually on the asymptote, which means that $h=t/n$. We will then show that the result is true when $h$ is suitably close to $t/n$.

\begin{figure}[!t]
\setkeys{Gin}{width=10cm}
\centerline{\includegraphics{fig02.02.1.eps}}
\caption[Iterated linear smoothing converges to the heat equation]{\label{22}Iterated linear smoothing converges to the heat equation. In this
experiment with one-dimensional functions, it can be appreciated how
fast an iterated convolution of a positive kernel converges to a
Gaussian. On the left are displayed nine iterations of the convolution of the characteristic function of an interval with itself, with appropriate rescalings. On the right, the same experiment is repeated with a much more irregular kernel. The convergence is almost as fast as the first case.}
\end{figure}


We begin with comments about the notation. By Exercise 2.1, $(T_h)^n=(T^n)_h$, so there is no ambiguity in writing $T^n_h$. We will be applying $T^n_h$ to the solution $u$ of the heat equation, which is $C$-periodic and $C^{\infty}$ on $(0,+\infty)\times\RR^N.$ In this situation, $t$ is considered to be a parameter, and we write $T^n_hu(t,\fxx)$ as shorthand for $T^n_hu(t,\cdot)(\fxx).$ Throughout the proof, we will be dealing with error terms that we write as $O(h^r)$. These terms invariably depend on $h$, $t$, and $\fxx$. However, in all cases, given a closed interval $[t_1,t_2]\subset(0,+\infty)$, there will be a constant $C$ such that $|O(h^r)|\leq Ch^r$ uniformly for $t\in[t_1,t_2]$ and $\fxx\in C$. Finally, keep in mind that everything in sight is $C$-periodic.

We wish to fix the interval $[t_1,t_2]$, but since this depends on the point $t$ in \eqref{37} and on $\varepsilon$, we must first choose these numbers. Thus, choose $\tau>0$ and keep it fixed. This will be the ``$t$'' in \eqref{37}. Next, choose $\varepsilon>0$. Here are the conditions we wish $t_1$ and $t_2$ to satisfy:
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item $t_1$ is small enough so $\|u(t_1,\cdot)-u_0\|_{L^1(C)}<\varepsilon$. (This is possible by Theorem 1.10.)
\item $t_1$ is small enough so $\|u(t_1+\tau,\cdot)-u(\tau,\cdot)\|_{L^1(C)}<\varepsilon$. (Again, by Theorem 1.10.)
\item $t_2$ is large enough so $t_1+\tau<t_2$.
\end{enumerate}

There is no problem meeting these conditions, so we fix the interval $[t_1,t_2]\subset(0,+\infty)$ and use Theorem 2.2 to write
\begin{equation}\label{38}
T_hu(t,\fxx)-u(t,\fxx)=h\Delta u(t,\fxx)+O(h^{{3}/{2}}),
\end{equation}
where $t\in[t_1,t_2]$. That the error function is bounded uniformly  by $Ch^{3/2}$ on the compact set $[t_1,t_2]\times C$ follows from the fact that $\sup_{\alpha,t,\fxx}|\partial^{\alpha}u(t,\fxx)|$ is finite for $(t,\fxx)\in[t_1,t_2]\times C$ (see the proof of Theorem 2.2). Since $u$ is a solution of the heat equation, we also have
\begin{equation}\label{39}
u(t+h,\fxx)-u(t,\fxx)=h\Delta u(t,\fxx)+O(h^2).
\end{equation}
This time the error term is bounded uniformly by $Ch^2$ on $[t_1,t_2]\times C$ because $u$ is $C^{\infty}$ on $(0,+\infty)\times\RR^N.$ By subtracting \eqref{39} from \eqref{38} we see that
\begin{equation}\label{310}
T_hu(t,\fxx)=u(t+h,\fxx)+O(h^{{3}/{2}}).
\end{equation}
This shows that applying $T_h$ to a solution of the heat equation at time $t$ advances the solution to time $t+h$, plus an error term.

So far we have not used the assumption that $g$ is nonnegative. Thus, \eqref{310} is true for any Laplacian-consistent kernel $g$ with $\sigma=2$. However, we now wish to apply the linear operator $T_h$ to both sides of equation \eqref{310}, and in doing so we do not want the error term to increase. Since $g\geq0$, this is not a problem:
\[
|T_hO(h^{{3}/{2}})|\leq \int_{\RR^N}|O(h^{{3}/{2}})|g_h(\fxx-\fxy)\,\rmd\fxy\leq \int_{\RR^N}Ch^{{3}/{2}}g_h(\fxx-\fxy)\,\rmd\fxy=Ch^{{3}/{2}}.
\]
With this in hand, we can apply $T_h$ to both sides of \eqref{310} and obtain
\begin{equation}\label{311}
T^2_hu(t,\fxx)=T_hu(t+h,\fxx)+O(h^{{3}/{2}}).
\end{equation}
If we write equation \eqref{310} with $t+h$ in place of $t$ and substitute the expression for $T_hu(t+h,\fxx)$ in equation \eqref{311}, we have
\begin{equation}
T^2_hu(t,\fxx)=u(t+2h,\fxx)+2O(h^{{3}/{2}}).
\end{equation}
We can iterate this process and get the relation
\begin{equation}
T^n_hu(t,\fxx)=u(t+nh,\fxx)+nO(h^{{3}/{2}})
\end{equation}
with the same constant $C$ in the estimate $|O(h^{{3}/{2}})|\leq Ch^{{3}/{2}}$ as long as $t+nh\in[t_1,t_2]$. To ensure that this happens, we take $t=t_1$ and $h=\tau/n$. Then
\begin{equation}\label{315}
T^n_hu(t_1,\fxx)=u(t_1+\tau,\fxx)+O\Big(\Big(\frac{\tau}{n}\Big)^{{1}/{2}}\Big)
\end{equation}
and
\begin{equation}
\lim_{n\to+\infty\atop nh=\tau}T^n_hu(t_1,\fxx)=u(t_1+\tau,\fxx),
\end{equation}
where the convergence is uniform for  $\fxx\in C$. If we could take $t_1=0$, this would end the proof. This is not possible because all of the $O$ terms were based on a fixed interval $[t_1,t_2]$. However, we have taken $t_1$  small enough to finish the proof .

By Proposition 1.6,
\[
\|g_h*v\|_{L^1(C)}\leq\|g_h\|_{L^1(\RR^N)}\|v\|_{L^1(C)}
\]
for $v\in L^1_C$. Since $g\geq0$ and $\int_{\RR^N}g(\fxx)\,\rmd\fxx=1$, $\|g_h\|_{L^1(\RR^N)}=1$, and thus
\[
\|g_h^{n*}*v\|_{L^1(C)}\leq\|v\|_{L^1(C)}.
\]
If we take $v=u(t_1,\cdot)-u_0$, then this inequality and condition (1) imply that
\begin{equation}\label{317}
\|T_h^{n}u(t_1,\cdot)-T_h^{n}u_0\|_{L^1(C)}<\varepsilon.
\end{equation}

If we write \eqref{315} as
\[
T^n_hu(t_1,\fxx)-u(t_1+\tau,\fxx)=O\Big(\Big(\frac{\tau}{n}\Big)^{{1}/{2}}\Big)
\]
and integrate both sides over $C$, then we also have
\begin{equation}\label{318}
\|T^n_hu(t_1,\cdot)-u(t_1+\tau,\cdot)\|_{L^1(C)}=O\Big(\Big(\frac{\tau}{n}\Big)^{{1}/{2}}\Big)<\varepsilon
\end{equation}
for all sufficiently large $n$, say $n>n_0$. Inequalities \eqref{317} and \eqref{318} imply that
\begin{equation}\label{319}
\|T_h^{n}u_0-u(t_1+\tau,\cdot)\|_{L^1(C)}<2\varepsilon.
\end{equation}
This inequality and condition (2) show that
\begin{equation}\label{215}
\|T_h^{n}u_0-u(\tau,\cdot)\|_{L^1(C)}<3\varepsilon
\end{equation}
for $n>n_0$ and $h=\tau/n$. This proves the theorem in the case $h=\tau/n$.

It is a simple matter to obtain the more general result. Again, by Theorem 1.10, there is a $\delta=\delta(\tau,\varepsilon)$ such that $|nh-\tau|<\delta$ implies that $\|u(nh,\cdot)-u(\tau,\cdot)\|_{L^1(C)}<\varepsilon$ and that $nh\in[t_1,t_2]$ (by condition (3)). Combining this with \eqref{215} shows that
\[
\|T_h^{n}u_0-u(nh,\cdot)\|_{L^1(C)}<4\varepsilon
\]
if $n>n_0$ and $|nh-\tau|<\delta$, and this completes the proof.\hfill
\end{proof}

We will denote the characteristic function of a set $A\subset\RR^N$ by $\bmone_A$. Thus, $\bmone_A(\fxx)=1$ if $\fxx\in A$ and $\bmone_A(\fxx)=0$ otherwise.


\begin{exercise}The kernel in the example in section 1.1 was $g=(1/\pi)\bmone_{D(0,1)}$, where $D(0,1)$ is the disk of radius one centered at zero. In this case, $g$ is a radial function and it is clearly Laplacian consistent. For $N=2$, let $A=[-1/2,1/2]\times[-1/2,1/2]$. Then $g=\bmone_A$ is not radial. Show that it is, however, Laplacian consistent. If we take $B=[-1,1]\times[-1/2,1/2]$, then $g=(1/2)\bmone_B$ is no longer Laplacian consistent because it does not satisfy condition $(iv)$. Show that this kernel does, however, satisfy a relation similar to \eqref{LC}. \finex
\end{exercise}

\section{Directional averages and directional heat equations}

The purpose of this section is to show how to approximate directional second derivatives using local averages. Consider a $C^2$ function from $\RR^N$ into $\RR$ and a vector $\fxz\in\RR^N$ with $|\fxz|=1$. We wish to compute the mean value of $u$ along a segment of the line through $\fxx$ parallel to the vector $\fxz$. To do this, we define the operator $T^{\fxz}_h$, $h\in[-1,1]$, by
\[
T^{\fxz}_hu(\fxx)=\frac{1}{2h}\int_{-h}^h u(\fxx+s\fxz)\,\rmd s.
\]
This operator is the directional counterpart of the isotropic operator $M_h$ defined by equation (1.1). We use Taylor's formula to expand $u$ at the point $\fxx$ along the line through $\fxx$ parallel to the vector $\fxz$:
\begin{equation}\label{linez}
u(\fxx+s\fxz)=u(\fxx)+sDu(\fxx)\cdot\fxz+\frac{s^2}{2}D^2u(\fxx)(\fxz,\fxz)+o(s^2).
\end{equation}
By averaging both sides of \eqref{linez} over $s\in[-h,h]$, we obtain the next result.
\begin{proposition}
\[
T^{\fxz}_hu(\fxx)=u(\fxx)+\frac{h^2}{6}D^2u(\fxx)(\fxz,\fxz)+o(h^2).
\]
\end{proposition}
Proposition 2.4 is similar to equation (1.6) and to Theorem 2.2, and it suggests that iterations of the operator $T^{\fxz}_h$ are associated with the directional heat equation
\begin{equation}\label{dirheateqn1}
\frac{\partial u}{\partial t}(t,\fxx)=\frac{1}{6}D^2u(t,\fxx)(\fxz,\fxz)
\end{equation}
in the same way that the iterations of the operator $T_h$ in Theorem 2.3 are associated with the ordinary heat equation. If $\fxz$ is fixed, then the operator $T^{\fxz}_h$ and equation \eqref{dirheateqn1} act on $u$ along each line in $\RR^N$ parallel to $\fxz$ separately; there is no ``cross talk'' between lines. Exercise 2.4 formalizes and clarifies these comments.

\begin{exercise}
Let $\fxz$ be a fixed vector in $\RR^N$ with $|\fxz|=1$ and let $u_0$ be in $L^1_C$. Define the kernel $g$ by $g(s)=\frac{1}{2}{\bmone}_{[-1,1]}(s)$.
\begin{enumerate}[(i)]\setlength\itemsep{2pt}
\item Show that $g$ satisfies the conditions of Definition 2.1 with $N=1$. Compute the variance $\sigma$ of $g$.
\item Show that
\[
u(t,\fxx)=\int_\RR u_0(\fxx+s\fxz)G_t(s)\,\rmd s
\]
is a solution of the directional heat equation
\begin{equation}\label{dirheateqn2}
\frac{\partial u}{\partial t}(t,\fxx)=D^2u(t,\fxx)(\fxz,\fxz),
\end{equation}
in the sense that $u$ satisfies \eqref{dirheateqn2} and
\[
\int_C |u(t,\fxx)-u_0(\fxx)|\,\rmd\fxx \to 0\;\text{ as }\; t\to 0.
\]
Give an example to show that $u(t,\cdot)$ is not necessarily $C^2$. This being the case, how does one interpret the right-hand side of \eqref{dirheateqn2}?
\item Let $g_h(s)=(6h)^{-1/2}g(s/(6h)^{-1/2})$ and $T_hu(\fxx)=\int_{\RR}u(\fxx+s\fxz)g_h(s)\,\rmd s$. By applying Theorem 2.3 for $N=1$, show that, for each $t>0$,
\begin{equation}
T_h^nu_0\to u(t,\cdot)\,\text{ in }\, L^1(C)\,\text{ as }\, n\to+\infty\,\text{ and }\, nh\to t.\, \finex
\end{equation}
\end{enumerate}
\end{exercise}

Exercise 2.4 has been established under the assumption that $\fxz$ is fixed. However, Proposition 2.4 is true when $\fxz$ is a function of $\fxx$. This means that we are able to approximate the directional second derivative by taking directional averages when $\fxz$ varies from point to point. The main choices considered in the book are $\fxz={Du}/{|Du|}$ and $\fxz={Du^{\perp}}/{|Du|}$, where $Du=(u_x,u_y)$ and $Du^{\perp}=(-u_y,u_x)$. Then by Proposition 2.4 we have the following limiting relations:
\begin{enumerate}[$\bullet$]\setlength\itemsep{2pt}
\item Average in the direction of the gradient. By choosing $z={Du}/{|Du|}$,
\[
\frac{1}{|Du|^2}D^2u(Du,Du)=6\lim_{h\to0}\frac{T^{Du/|Du|}_hu-u}{h^2}.
\]
We will interpret this differential operator as Haralick's edge detector in section 2.3.
\item Average in the direction orthogonal to the gradient. By choosing $z={Du}/{|Du^\perp|}$,
\[
\frac{1}{|Du|^2}D^2u(Du^\perp,Du^\perp)=6\lim_{h\to0}\frac{T^{Du^\perp/|Du|}_hu-u}{h^2}.
\]
This differential operator appears as the second term of the curvature equation.
\end{enumerate}
Although we have not written them as such, the limits are pointwise in both cases.

This section has been about approximating directional second derivatives by computing directional averages. Exercise 2.5 provides another method for computing second derivatives.


\begin{exercise} Suppose that $u\in C^2(\RR)$. Assuming that $u'(x)\not=0$, show that
\begin{equation}\label{ex251}
u''(x)=\lim_{h\to0}\frac{1}{h^2}\Big(\max_{s\in[-h,h]}u(x+s)+\min_{s\in[-h,h]}u(x+s)-2u(x)\Big).
\end{equation}
What is the value of the right-hand side of \eqref{ex251} if $u'(x)=0$?

Now consider $u\in C^2(\RR^2)$. We wish to establish an algorithm similar to \eqref{ex251} to compute the second derivative of $u$ in the direction of the gradient $Du=(u_x,u_y)$. For this to make sense, we must assume that $Du(\fxx)\not=0$. With these assumptions, we know from \eqref{ex251} that
\begin{equation}\label{ex252}
u_{\xi\xi}(\fxx)=\frac{\partial^2 v}{\partial \xi^2}(\fxx,0)=\lim_{h\to0}\frac{1}{h^2}\Big(\max_{s\in[-h,h]}u(\fxx+s\fxz)+\min_{s\in[-h,h]}u(\fxx+s\fxz)-2u(\fxx)\Big),
\end{equation}
where $v(\fxx,\xi)=u(\fxx+\xi\fxz)$ and $\fxz=Du/|Du|$. The second part of the exercise is to prove that, in fact,
\begin{equation}\label{ex253}
u_{\xi\xi}(\fxx)=\lim_{h\to0}\frac{1}{h^2}\Big(\max_{\fxy\in D(0,h)}u(\fxx+\fxy)+\min_{\fxy\in D(0,h)}u(\fxx+\fxy)-2u(\fxx)\Big),
\end{equation}
where $D(0,h)$ is the disk of radius $h$ centered at the origin. Intuitively, \eqref{ex253} follows from \eqref{ex252} because the gradient indicates the direction of maximal change in $u(\fxx)$, so in the limit as $h\to0$, taking max and min in the direction of the gradient is equivalent to taking max and min in the disk. The point of the exercise is to formalize this. \finex
\end{exercise}

\section{Edge detection and linear scale space}

The general analysis framework in which an image is associated with smoothed versions of itself at several scales is called \textit{scale space}.
We are first going to discuss some of the ways linear smoothing has been used in image analysis. We will then indicate the problems associated with the linear theory.

\subsection{The edge detection doctrine}

One of the uses of linear theory in two dimensions is \textit{edge detection}. The assumption of the edge detection doctrine is that relevant information is contained in the traces produced in an image by the apparent contours of physical objects. If a black object is photographed against a white background, then one expects the silhouette of the object in the image to be bounded by a closed curve across which the light intensity $u_0$ varies strongly. We call this curve an \textit{edge}. At first glance, it would seem that this edge could be detected by computing the gradient $Du_0$, since at a point $\fxx$ on the edge, $|Du_0(\fxx)|$ should be large and $Du(\fxx)$ should point in a direction normal to the boundary of the silhouette. It would therefore appear that finding edges amounts to computing the gradient of $u_0$ and determining the points where the gradient is large. This conclusion is unrealistic for two reasons:
\begin{enumerate}[(a)]\setlength\itemsep{2pt}
\item There may be many points where the gradient is large due to small oscillations in the image that are not related to real objects. Recall that digital images are always noisy, and thus there is no reason to assume the existence or computability of a gradient.
\item The points where the gradient exceeds a given threshold are likely to form regions and not curves.
\end{enumerate}

As we emphasized in the Introduction, objection (a) is dealt with by smoothing the image. We associate with the image $u_0$ smoothed versions $u(t,\cdot)$, where the scale parameter $t$ indicates the amount of smoothing. In the classical linear theory, this smoothing is done by convolving $u_0$ with the Gaussian $G_t$.

One way that objection (b) has been approached is by redefining edge points. Instead of just saying an edge point is a point $\fxx$ where $|Du_0(\fxx)|$ exceeds a threshold, one requires the gradient to satisfy a maximal property. We illustrate this in one dimension. Suppose that $u\in C^2(\RR)$ and consider the points where $|u'(x)|$ attains a local maximum. At some of these points, the second derivative $u''$ changes sign, that is, ${\rm{sign}}(u''(x-h))\not={\rm{sign}}(u''(x+h))$ for sufficiently small $h$. These are the points where $u''$ \textit{crosses zero}, and they are taken to be the edge points. Note that this criterion avoids classifying a point $x$ as an edge point if the gradient is constant in an interval around $x$. Marr and Hildreth generalized this idea to two dimensions by replacing $u''$ with the Laplacian $\Delta u$, which is the only isotropic linear differential operator of order two that generalizes $u''$ \cite{Marr:1980:TED}. Haralick's edge detector is different but in the same spirit \cite{Har2}. Haralick gives up linearity and define edge points as those points where the gradient has a local maximum in the direction of the gradient. In other words, an edge point $\fxx$ satisfies $g'(0)=0$, where $g(t)=|Du(\fxx+tDu(\fxx)|/|Du(\fxx)|$. This implies that $D^2u(\fxx)(Du(\fxx),Du(\fxx))=0$.


\begin{exercise}
Check that $g'(0)=0$ implies that $D^2u(\fxx)(Du(\fxx),Du(\fxx))=0$. \finex
\end{exercise}

We are now going to state these two algorithms formally. They are illustrated in Figures \ref{fig24} and \ref{fig26}, respectively.


\begin{algorithm}[Edge detection: Marr--Hildreth zero-crossings]\leavevmode
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item Create the multiscale images $u(t,\cdot)=G_t*u_0$ for increasing values of $t$.
\item At each scale $t$, compute all the points where $\Delta u$ changes sign and $Du\not=0$. These points are called zero-crossings of the Laplacian, or simply zero-crossings.
\item (Optional) Eliminate the zero-crossings where the gradient is small.
\end{enumerate}
\end{algorithm}

\begin{algorithm}[Edge detection: The Haralick--Canny edge detector]\leavevmode
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item As before, create the multiscale images $u(t,\cdot)=G_t*u_0$ for increasing values of $t$.
\item At each scale $t$, find all points $\fxx$ where $Du(\fxx)\not=0$ and $D^2u(\fxx)(\fxz,\fxz)$ crosses zero, $\fxz=Du/|Du|$. At such points, the function $s\mapsto u(\fxx+s\fxz)$ changes from concave to convex, or conversely, as $s$ passes through zero.
\item At each scale $t$, fix a threshold $\theta(t)$ and retain as \textit{edge points at scale} $t$ only those points found above that satisfy $|Du(\fxx)|>\theta(t)$.
\end{enumerate}
\end{algorithm}


\noindent In practice, edges are computed for a finite number of dyadic scales, $t=2^n$, $n\in\ZZ$.

\subsection{Discussion and critique}

The Haralick--Canny edge detector is generally preferred for its accuracy to the Marr--Hildreth algorithm. Their use and characteristics are, however, essentially the same. There are also many variations---attempted improvements---of the algorithms we have described, and the following discussion adapts easily to these related edge detection schemes. The first thing to notice is that, by Theorem 1.10, $u(t,\cdot)=G_t*u_0$ is a $C^\infty$ function for each $t>0$ if $u_0\in L^1_C$. Thus we can indeed compute second order differential operators applied to $u(t,\cdot)=G_t*u_0$, $t>0$. In the case of linear operators like the Laplacian or the gradient, the task is facilitated by the formula proved in Proposition 1.8. For example, we have $\Delta u(t,\fxx)=\Delta(G_t*u_0)(\fxx)=(\Delta G_t)*u_0(\fxx)$, where in dimension two (Figure \ref{fig23}),
\[
\Delta G_t(\fxx)=\frac{|\fxx|^2-4t}{16\pi t^3}{\mathrm e}^{-{|\fxx|^2}/{4t}}.
\]

\begin{figure}[!t]
\setkeys{Gin}{width=8cm}
\centerline{
\includegraphics{fig02.03.1.eps}
}
\caption[The Laplacian of the Gaussian]{\label{fig23}A three-dimensional representation of the Laplacian of the Gaussian. This convolution kernel, which is a wavelet, is used to estimate the Laplacian of an image at different scales of linear smoothing.}
\end{figure}

\vspace{1\baselineskip}

\begin{figure}[!h]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.04.1.eps}
}
\caption[Zero-crossing of the Laplacian at different scales]{\label{fig24}
Zero-crossings of the Laplacian at different scales. This figure illustrates
the original scale space theory as developed by David Marr \cite{Marr:1982:V}. To extract more global structure, the image is convolved with Gaussians whose variances are powers of two. One computes the Laplacian of the smoothed image and displays the lines along which this Laplacian changes sign: the zero-crossings of the Laplacian. According to Marr, these zero-crossings represent the ``raw primal sketch" of the image, or the essential information on which further vision algorithms should be based. Above, left to right: the results of smoothing and the associated Gaussian kernels at scales 1, 2, and 4. Below, left to right: the zero-crossings of the Laplacian and the corresponding kernels, which are the Laplacians of the Gaussians used above.}
\end{figure}

In the same way, Haralick's edge detector makes sense because $u$ is $C^\infty$ at all points where $Du(\fxx)\not=0$. If $Du(\fxx)=0$, then $\fxx$ cannot be an edge point, since $u$ is ``flat'' there. Thus, thanks to the filtering, there is no theoretical problem with computing edge points. There are, however, practical objections to these methods, which we will now discuss.


\subsubsection*{Linear scale space}

The first serious problems are associated with the addition of an extra dimension: Having many images $u(t,\cdot)$ at different scales $t$ confounds our understanding of the image and adds to the cost of computation. We no longer have an absolute definition of an edge. We can only speak of edges at a certain scale. Conceivably, a way around this problem would be to track edges across scales. In fact, it has been observed in experiments that the ``main edges'' persist under convolution as $t$ increases, but they lose much of their spatial accuracy. On the other hand, filtering with a sharp low-pass filter, that is, with $t$ small, keeps these edges in their proper  positions, but eventually, as $t$ becomes very small, even these main edges can be lost in the crowd of spurious edge signals due to noise and texture. The scale space theory of Witkin proposes to identify the main edges at some scale $t$ and then to track them backward as $t$ decreases \cite{Witkin:1983:SSF}. In theory, it would seem that this method could give an accurate location of the main edges. In practice, any implementation of these ideas is computationally costly due to the problems involved with multiple thresholdings and following edges across scales. In fact, tracking edges across scales is incompatible with having thresholds for the gradients, since such thresholds may remove edges at certain scales and not at others. The conclusion is that one should trace all zero-crossings across scales without considering whether they are true edges or not. This makes matching edges across scales very difficult. For example, experiments show that zero-crossings of sharp edges that are sparse at small scales are no longer sparse at large scales. (Figure \ref{fig25} shows how zero-crossings can be created by linear smoothing.) The Haralick--Canny detector suffers from the same problems, as is well demonstrated by experiments.

\begin{figure}[!t]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.06.1.eps}
}
\caption[Canny's edge detector]{\label{fig26}Canny's edge detector. These images illustrate the Canny edge detector. Left column: result of the Canny filter without the threshold on the gradient. Middle column: result with a visually ``optimal'' scale and an image-dependent threshold (from top to bottom: 15, 0.5, 0.6). Right column: result with a fixed gradient threshold equal to 0.5. Note that this edge detection ``theory'' depends on no fewer than three parameters that must be fixed by the user.
}
\end{figure}


Other problems with linear scale space are illustrated in Figures \ref{fig28} and \ref{fig27}. Figure \ref{fig28} illustrates how linear smoothing can create new gray levels and new extrema. Figure \ref{fig27} shows that linear scale space does not maintain the inclusion between objects. The shape inclusion principal will be discussed in Chapter 21.

We must conclude that the work on linear edge detection has been an attempt to build a theory that has not succeeded. After more than thirty years of activity, it has become clear that no robust technology can be based on these ideas. Since edge detection algorithms depend on multiple thresholds on the gradient, followed by ``filling-the-holes'' algorithms, there can be no scientific agreement on the identification of edge points in a given image. In short, the problems associated with linear smoothing followed by edge detection have not been resolved by the idea of chasing edges across scales.


\begin{figure}[!t]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.05.1.eps}
}
\caption[Zero-crossings of the Laplacian of a synthetic image]{\label{fig25}Zero-crossings of the Laplacian of a synthetic image. Left to right: the original image; the image linearly smoothed by convolution with
a Gaussian; the sign of the Laplacian of the filtered image (the gray color corresponds to values close to 0, black to clear-cut negative values, white to clear-cut positive values); the zero-crossings of the Laplacian. This experiment clearly shows a drawback of the Laplacian as edge detector.
}
\end{figure}

\begin{figure}[!ht]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.08.1.eps}
}
\caption[The heat equation creates structure.]{\label{fig28}The heat equation creates structure. This experiment shows that linear scale space can create new structures and thus increase the complexity of an image. Left to right: The original synthetic image (a) contains three gray levels. The black disk is a regional and absolute minimum. The ``white'' ring around the black disk is a regional and absolute maximum. The outer gray ring has a gray value between the other two and is a regional minimum. The second image (b) shows what happens when (a) is smoothed with the heat equation: New extrema have appeared. Image (c) illustrates the action on (a) of a contrast-invariant local filter, the iterated median filter, which is introduced in Chapter 9.
}
\end{figure}

\begin{figure}[!ht]
\setkeys{Gin}{width=10cm}
\centerline{
\includegraphics{fig02.07.1.eps}
}
%\centerline{
%\includegraphics{fig02.07.2.eps}
%}
\caption[Violation of the shape inclusion principle by the linear scale space]{\label{fig27}Violation of the inclusion by the linear scale space.
Top, left: an image that contains a black disk enclosed by a white disk. Top, right: At a certain scale, the black and white circles mix together. Bottom, left: The boundaries of the two circles. Bottom, right: After smoothing with a certain value of $t$, the inclusion that existed for very small $t$ in no longer preserved. The topology of the situation has changed.
}
\end{figure}


\subsubsection*{Contrast invariance}

As already mentioned in the Introduction, a central theme of the book is that the use of contrast-invariant operators will solve most of the technical problems associated with linear smoothing and other linear image operators. The development of these ideas starts in Chapter 3.

Recall from section I.3 that an (image) operator $u\mapsto Tu$ is contrast invariant if $T$ commutes with all nondecreasing functions $g$, that is, if
\begin{equation}\label{contrast}
g(Tu)=T(g(u)).
\end{equation}
If image analysis is to be robust, it must be invariant under changes in lighting that produce contrast changes. It must also be invariant under the nonlinear response of the sensors used to capture an image. These, and perhaps other, contrast changes are modeled by $g$. If $g$ is strictly increasing, then relation \eqref{contrast} ensures that the filtered image $Tu=g^{-1}(T(g(u)))$ does not depend on $g$. A problem with linear theory is that linear smoothing, that is, convolution, is not generally contrast invariant:
\[
g(k*u)\not=k*(g(u)).
\]
In the same way, the operator $T_t$ that maps $u_0$ into the solution of the heat equation, $u(t,\cdot)$ is not generally contrast invariant. In fact, if $g$ is $C^2$, then

\begin{align*}
\frac{\partial(g(u))}{\partial t}&=g'(u)\frac{\partial u}{\partial t}\\[-2ex]
\intertext{and}\\[-5ex]
\Delta(g(u))&=g'(u)\Delta u + g''(u)|Du|^2.
\end{align*}
Of course, if $g$ is linear, then $g(u)$ satisfies the heat equation if $u$ does.

\begin{exercise}
Construct simple functions $u$, $g$, and $k$ such that $g(k*u)\not=k*(g(u))$. \finex
\end{exercise}

\section{Dynamic shape}

In 1986, Koendernik and van Doorn defined a \textit{shape} in $\RR^N$ to be a closed subset $X$ of $\RR^N$ \cite{Koenderink:1986:DS}. They then proposed to smooth the shape by applying the heat equation $\partial u/\partial t-\Delta u=0$ directly to $\bmone_X$, the characteristic function of $X$. Of course, the solution $G_t*\bmone_X$ is not a characteristic function. The authors defined the evolved shape at scale $t$ to be
\[
X_t=\{\fxx\mid u(t,\fxx)\geq 1/{2}\}.
\]
The value $1/2$ is chosen so the following simple requirement is satisfied: Suppose that $X$ is the half-plane $X=\{(x,y)\mid (x,y)\in\RR^2, x\geq0\}$. The requirement is that
\[
X=X_t=\{(x,y)\mid G_t*\bmone_X(x,y)\geq\lambda\},
\]
and this is true only if $\lambda=1/2$. There are at least two problems with dynamic shape evolution for image analysis. The first concerns nonlocal interactions, as illustrated in Figure \ref{fig210}. Here we have two disks that are near one another. The evolution of the union of both disks, considered as a single shape, is quite different from the evolution of the disks separately. A related problem, also illustrated in Figure \ref{fig210}, is the creation of singularities. Note how a singularity in orientation and the curvature of the boundary of the shape develops at the point where the two disks touch. Figure \ref{fig29} further illustrates the problems associated with the dynamic shape method.

\begin{figure}[!t]
\setkeys{Gin}{width=3cm}
\centerline{
\includegraphics{fig02.10.1.eps}\hspace{-10pt}
\includegraphics{fig02.10.2.eps}\hspace{-10pt}
\includegraphics{fig02.10.3.eps}\hspace{-10pt}
\includegraphics{fig02.10.4.eps}
}
\caption[Dynamic shape and nonlocal interactions.]{\label{fig210} Nonlocal interactions in the dynamic shape method. Left to right: Two close disks interact as the scale increases. This creates a new, qualitatively different, shape. The change of topology, at the scale where the two disks merge into one shape, also entails the appearance of a singularity (a cusp) on the shape(s) boundaries.}
\end{figure}










\section{Curve evolution using the heat equation}

We consider shapes in $\RR^2$ whose boundaries can be represented by a finite number of simple closed rectifiable Jordan curves. Thus, each curve  we consider can be represented by a continuous mapping $f:[0,1]\to\RR^2$ such that $f$ is one-to-one on $(0,1)$ and $f(0)=f(1)$, and each curve has a finite length. We also assume that these curves do not intersect each other. We will focus on smoothing one of these Jordan curves, which we call $C_0$. We assume that $C_0$ is parameterized by $s\in[0,L]$, where $L$ is the length of the curve. Thus, $C_0$ is represented as $\fxx_0(s)=(x(s),y(s))$, where $s$ is the length of the curve between $\fxx_0(0)$ and $\fxx_0(s)$.

\begin{figure}[!t]
\setkeys{Gin}{width=4cm}
\centerline{
\includegraphics{fig02.09.1.eps}\hspace{-4pt}
\includegraphics{fig02.09.2.eps}\hspace{-4pt}
\includegraphics{fig02.09.3.eps}
}\vspace{-20pt}
\centerline{
\includegraphics{fig02.09.4.eps}\hspace{-4pt}
\includegraphics{fig02.09.5.eps}\hspace{-4pt}
\includegraphics{fig02.09.6.eps}
}
\caption[Dynamic shape, nonlocal behavior of shapes.]{\label{fig29}  Nonlocal behavior of shapes with the dynamic shape method. This image displays the smoothing of two irregular shapes by the dynamic shape method (Koenderink--van Doorn). Top left: initial image, made of two irregular shapes. From left to right, top to bottom: dynamic shape smoothing with increasing Gaussian variance. Notice how the shapes merge more and more. We do not have a separate analysis of each shape but rather a ``joint analysis'' of the two shapes. The way the shapes merge is of course sensitive to the initial distance between the shapes. Compare with Figure \ref{fig212}.
}
\end{figure}

At first glance, it might seem reasonable to smooth $C_0$ by smoothing the coordinate functions $x$ and $y$ separately. If this is done linearly, we have seen from Theorem 2.3 that the process is asymptotic to smoothing with the heat equation. Thus, one is led naturally to consider the vector heat equation
\begin{equation}\label{vheat}
\frac{\partial \fxx}{\partial t}(t,s)=\frac{\partial^2 \fxx}{\partial s^2}(t,s)
\end{equation}
with initial condition $\fxx(0,s)=\fxx_0(s)$. If $\fxx(t,s)=(x(t,s),y(t,s))$ is the solution of \eqref{vheat}, then we know from Proposition 1.9 that
\begin{align*}
\inf_{s\in[0,L]}x_0(s)&\leq x(t,s)\leq\sup_{s\in[0,L]}x_0(s),\\
\inf_{s\in[0,L]}y_0(s)&\leq y(t,s)\leq\sup_{s\in[0,L]}y_0(s),
\end{align*}
for $s\in[0,L]$ and $t\in[0,+\infty)$. Thus, the evolved curves $C_t$ remain in the rectangle that held $C_0$. Also, we know from Theorem 1.10 that the coordinate functions $x(t,\cdot)$ and $y(t,\cdot)$ are $C^{\infty}$ for $t>0$. There are, however, at least two reasons that argue against smoothing curves this way:
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item When $t>0$, $s$ is no longer a length parameter for the evolved curve $C_t$.
\item Although $x(t,\cdot)$ and $y(t,\cdot)$ are $C^{\infty}$ for $t>0$, this does not imply that the curves $C_t$ have similar smoothness properties. In fact, it can be seen from Figure \ref{fig211} that it is possible for an evolved curve to cross itself and it is possible for it to develop singularities.
\end{enumerate}

\begin{figure}[!t]
\setkeys{Gin}{width=6cm}
\centerline{
\includegraphics{fig02.11.1.eps}
\includegraphics{fig02.11.2.eps}
}
\centerline{
\includegraphics{fig02.11.3.eps}
\includegraphics{fig02.11.4.eps}
}
\caption[Curve evolution by the heat equation.]{\label{fig211} Curve evolution by the heat equation. The coordinates of the curves are parameterized by the arc length and then smoothed as real functions of the length using the heat equation. From A to D: the coordinates are smoothed with an increasing scale. Each coordinate function therefore is $C^\infty$;  the evolving curve can, however, develop self-crossings (as in C) or singularities (as in D).}
\end{figure}

In image processing, we say that a process that introduces new features, such as described in item (2) above, is not \textit{causal}. (This informal definition should not be confused with the use of ``causality,'' as it is used, for example, when speaking about filters: A filter $F$ is said to be causal, or realizable, if the equality of two signals $s_0$ and $s_1$ up to time $t_0$ implies that $Fs_0(t)=Fs_1(t)$ for the same period.)

\begin{exercise}
Construct a $C^{\infty}$ mapping $f:[0,1]\to\RR^2$ such that the image of $[0,1]$ is a square. This shows that a curve can have a $C^{\infty}$ parameterization without being smooth. \finex
\end{exercise}



\section{Restoring locality and causality}

The main objective of this book is to redefine the smoothing processes so they are local and do not create new singularities. This can be done, as we prove in the following chapters, by alternating a small-scale linear convolution with a natural renormalization process. We end this section by introducing two algorithms as a preview to these later developments.

\subsection{Localizing the dynamic shape method}

In the case of dynamic shape analysis, we define an alternate dynamic shape algorithm as follows:
\begin{algorithm}[The Merriman--Bence--Osher algorithm]\leavevmode
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item Convolve the characteristic function of the initial shape $X_0$ with $G_h$, where $h$ is small.
\item Define $X_1=\{\fxx\mid G_h*\bmone_{X_0}\geq 1/2\}$.
\item Set $X_0=X_1$ and go back to (1).
\end{enumerate}
\end{algorithm}
This is an iterated dynamic shape algorithm. The dynamic shape method itself  is an example of a \textit{median filter}, which will be defined in Chapter 9. The Merriman--Bence--Osher algorithm is thus an \textit{iterated median filter} (see Figure \ref{fig212}). We will see in Chapters 10 and 11 that median filters have asymptotic properties that are similar to those expressed in Theorem 2.3. In the case of median filters, the associated partial differential equation will be a curvature motion equation (defined in Chapter 5).

\begin{figure}[!t]
\setkeys{Gin}{width=4cm}
\centerline{
\includegraphics{fig02.12.1.eps}\hspace{-4pt}
\includegraphics{fig02.12.2.eps}\hspace{-4pt}
\includegraphics{fig02.12.3.eps}
}\vspace{-20pt}
\centerline{
\includegraphics{fig02.12.4.eps}\hspace{-4pt}
\includegraphics{fig02.12.5.eps}\hspace{-4pt}
\includegraphics{fig02.12.6.eps}
}
\caption[Localized, iterated version of the dynamic shape.]{\label{fig212} The Bence--Merriman--Osher shape smoothing method is a localized and iterated version of the dynamic shape method. A convolution of the binary image with small-sized Gaussians is alternated with mid-level thresholding.  It uses the same initial data (top, left) as in Figure \ref{fig29}.  From left to right, top to bottom: smoothing with increasing scales. Notice that the shapes remain separate. In fact, their is no interaction between the evolving shapes. Each one evolves as if the other did not exist.
}
\end{figure}

\subsection{Renormalized heat equation for curves}

In 1992, Mackworth and Mokhtarian noticed the loss of causality when the heat equation was applied to curves \cite{Mackworth:1992:TMC}. Their method to restore causality looks, at least formally, like the remedy given for the nonlocalization of the dynamic shape method. Instead of applying the heat equation for relatively long times (or, equivalently, convolving the curve $\fxx$ with the Gaussian $G_t$ for large $t$), they use the following algorithm:
\begin{algorithm}[Renormalized heat equation for curves]\leavevmode
\begin{enumerate}[(1)]\setlength\itemsep{2pt}
\item Convolve the initial curve $\fxx_0$, parameterized by its length parameter $s_0\in[0,L_0]$, with the Gaussian $G_h$, where $h$ is small.
\item Let $L_n$ denote the length of the curve $\fxx_n$ obtained after $n$ iterations and let $s_n$ denote its length parameter. For $n\geq1$, write $\tilde\fxx_{n+1}(s_n)=G_h*\fxx_n(s_n)$. Then reparameterize $\tilde\fxx_{n+1}$ by its length parameter $s_{n+1}\in[0,L_{n+1}]$, and denote it by $\fxx_{n+1}$.
\item Iterate.
\end{enumerate}
\end{algorithm}

This algorithm is illustrated in Figure \ref{fig213}; it should be compared with Figure \ref{fig211}.

\begin{theorem} Let $\fxx$ be a $C^2$ curve parameterized by its length parameter $s\in[0,L]$. Then for small $h$,
\begin{equation}\label{225}
G_h*\fxx(s)-\fxx(s)=h\frac{\partial^2 \fxx}{\partial s^2} +o(h).
\end{equation}
\end{theorem}
\begin{exercise}
Prove Theorem 2.5. If $\fxx$ is a $C^3$ function of $s$, then the result follows directly from Theorem 2.2. The result holds, however, for a $C^2$ curve. \finex
\end{exercise}

\begin{figure}[!t]
\setkeys{Gin}{width=6cm}
\centerline{
\includegraphics{fig02.13.1.eps}
\includegraphics{fig02.13.2.eps}
}
\centerline{
\includegraphics{fig02.13.3.eps}
\includegraphics{fig02.13.4.eps}
}
\caption[Curve evolution by the renormalized heat equation.]{\label{fig213} Curve evolution by the renormalized heat equation (Mackworth--Mokhtarian). After each smoothing step, the coordinates of the curve are reparameterized by the arc length of the smoothed curve. From A to D: the curve is smoothed with an increasing scale. Note that, in contrast with the linear heat equation (Figure \ref{fig211}), the evolving curve shows no singularities and does not cross itself.}
\end{figure}



In view of \eqref{225} and what we have seen regarding asymptotic limits in Theorem 2.3 and Exercise 2.4, it is reasonable to conjecture that, in the asymptotic limit, Algorithm 2.4 will yield the solution of following evolution equation:
\begin{equation}\label{226}
\frac{\partial \fxx}{\partial t}=\frac{\partial^2 \fxx}{\partial s^2},
\end{equation}
where $\fxx_0=\fxx(0,\cdot)$. It is important to note that \eqref{226} is \textit{not} the heat equation \eqref{vheat}. Indeed, from Algorithm 2.4 we see that $s$ must denote the length parameter of the evolved curve $\fxx(t,\cdot)$ at time $t$. Since ${\partial^2 \fxx}/{\partial s^2}$ is the curvature vector $\text{Curv}(\fxx)$, this is easier to remember if we write \eqref{226} as
\[
\frac{\partial \fxx}{\partial t}=\text{Curv}(\fxx).
\]
We will study this nonlinear curve evolution equation in Chapter 5.

\section{Comments and references}

\paragraph*{Smoothing and the Laplacian.}
One of the first tools proposed in the early days of image processing in the 1960s came, not surprisingly, directly from signal processing. The idea was to restore an image by averaging the gray levels locally (see, for example, \cite{Golay:1972:SDL} and \cite{Helstrom:1967:IRM}). The observation that the difference between an image and its local average is proportional to the Laplacian of the image has proved to be one of the most fruitful contributions to image processing. As noted in the Introduction, this method for deblurring an image was introduced by Kovasznay and Joseph in 1955 \cite{Kovasznay:1955}, and it was studied and optimized by Gabor in 1965 \cite{Gabor:1965} (information taken from \cite{Lindenbaum:1994:GCI}). (See also \cite{Hummel:1983:GBH} and \cite{Hummel:1987:DGB}.) Burt and Adelson based their Laplacian pyramid algorithm on this idea, and this was one of the results that led to multiresolution analysis and wavelets \cite{Burt:1983:LPC}.

\paragraph*{Asymptotics.}
Our proof that iterated and rescaled convolutions of a Laplacian-consistent kernel tend asymptotically to the Gaussian is a version of the De Moivre--Laplace formula, or the central limit theorem, adapted to image processing \cite{Breiman:1968:P}. This result is particularly relevant to image analysis, since it implies that iterated linear smoothing leads inevitably to convolution with the Gaussian, or equivalently, to the application of the heat equation. We will prove this again in Chapter 19 in the context of the axiomatic development of scale space. We do not wish to imply, however, that the Gaussian is the only important kernel for image processing. The Gaussian plays a significant role in our form of image analysis, but there are other kernels that, because of their spectral and algebraic properties, have equally important roles in other aspects of signal and image processing. This is particularly true for wavelet theory.

\paragraph*{Directional diffusion.}
Directional diffusion has a long history that began when Hubel and Wiesel showed the existence of direction-sensitive cells in the visual areas of the neocortex \cite{Hubel}. There has been an explosion of publication on directional linear filters, beginning, for example, with influential papers such as that by Daugman \cite{Daugman:1985:URR}. We note again that Gabor's contribution to directional filtering is described in \cite{Lindenbaum:1994:GCI}.

\paragraph*{Gaussian smoothing and edge detection.}
The use of Gaussian filtering in image analysis is so pervasive that it is impossible to point to a ``first paper.'' It is, however, safe to say that David Marr's famous book, \textit{Vision} \cite{Marr:1982:V}, and the original paper by Hildreth and Marr \cite{Marr:1980:TED} have had an immeasurable impact on edge detection and image processing in general. The term ``edge detection'' appeared as early as 1959 in connection with television transmission \cite{Julesz:1959:MCT}. The idea that the computation of derivatives of an image necessitates a previous smoothing has been aggressively developed by the Dutch school of image analysis \cite{Blom:1991:SDP,Florack:1994:LSS}. See also the books by Florack \cite{Florack:1997:IS}, Lindeberg \cite{Lindeberg:1994:SSTa}, and Romeny \cite{Romeny:1998:FEV}, and the paper \cite{Falzon:1994:SAD}. Haralick's edge detector \cite{Har2}, as implemented by Canny \cite{Canny:1986:CAE}, is probably the best known image analysis operator. A year after Canny's 1986 paper, Deriche published a recursive implementation of Canny's criteria for edge detection \cite{Deriche:1987:UCC}.

\paragraph*{Scale space.}
The term ``scale space'' was introduced by Witkin in 1983. He suggested  tracking the zero-crossings of the Laplacian of the smoothed image across scales \cite{Witkin:1983:SSF}. Yuille and Poggio proved that these zero-crossings can be tracked for one-dimensional signals \cite{Yuille:1985:FTZ}. Hummel and Moniot \cite{Hummel:1986:RBZa,Hummel:1989:RZC} and Yuille and Poggio \cite{Yuille:1986:STZ} analyzed the conjectures of Marr and Witkin according to which an image is completely recoverable from its zero-crossings at different scales. Mallat formulated Marr's conjecture as an algorithm in the context of wavelet analysis. He replaced the Gaussian with a two-dimensional cubic spline, and he used both the zero-crossings of the smoothed images and the nonzero values of the gradients at these points to reconstruct the image. This algorithm works well in practice, and the  conjecture was that these zero-crossings and the values of the gradients determined the image. A counterexample given by Meyer shows that this is not the case. Perfect reconstruction is possible in the one-dimensional case for signals with compact support if the smoothing kernel is the Tukey window, $k(x)=1+\cos x$ for $|x|\leq\pi$ and zero elsewhere. An account of the Mallat conjecture and these examples can be found in \cite{JMR}. Koenderink presents a general and insightful theory of image scale space in \cite{Koenderink:1984:SI}.

\paragraph*{Dynamic shape, curve evolution, and restoring causality.} Our account of the dynamic shape method is based on the well-known paper by Koenderink and van Doorn in which they introduced this notion \cite{Koenderink:1986:DS}. The curve evolution by the heat equation if from the first 1986 version of curve analysis proposed by Mackworth and Mokhtarian \cite{Mackworth:1986:SBD}. See also the paper by Horn and Weldon \cite{Horn:1986:FCC}. The algorithm that restores causality and locality to the dynamic shape method was discovered by Merriman, Bence, and Osher, who devised this algorithm for a totally different reason: They were looking for a cleaver numerical implementation of the mean curvature equation \cite{Merriman:1992:DGM}. There were model errors in the 1986 paper \cite{Mackworth:1986:SBD} that were corrected by the authors in their 1992 paper \cite{Mackworth:1992:TMC}. There, they also proposed the correct intrinsic equation. However, this 1992 paper contains several inexact statements about the properties of the intrinsic equation. The correct theorems and proofs can be found in a paper by Grayson written in 1987 \cite{Grayson:1987:HESb}.

\paragraph*{Topological change under smoothing.} We have included several figures that illustrate how essential topological properties of an image change when the image is smoothed with the Gaussian. Damon has made a complete analysis of the topological behavior of critical points of an image under Gaussian smoothing \cite{Damon98}.





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