%gmChap02.tex

\graphicspath{{./figures/}}

\chapter{Linear Scale Space and Edge Detection\label{ChapLinearScaleSpaceEdgeDetection}}


 The general analysis framework in  which an image is associated
with smoothed versions of itself at several scales is called \textit{scale space}.  Following the results of Chapter
\ref{ChapCOIFTTHE}, a linear scale space must be performed by applying the heat equation to the  image.  The main aim
of this smoothing is to find out {\it edges} in the image.  We shall first explain this doctrine. In the second
section, we discuss experiments and several serious objections to such an image representation.

\section{The edge detection doctrine\label{sectionTEDD}}

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 defines 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$ (see Exercise
\ref{calculsHaralick}). 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  $Du\not=0$ and $\Delta u$ changes sign. These points are called zero-crossings of the Laplacian, or simply zero-crossings.
\item (Optional) Eliminate the zero-crossings where the gradient is below some prefixed threshold.
\item track back from large scales to fine scales the ``main edges" detected at large scales.
\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)$.  The backtracking step across scales is the same as for Marr--Hildreth.
\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 Proposition \ref{proposition1.8}, $u(t,\cdot)=G_t*u_0$ is a $C^\infty$ function
for each $t>0$ if $u_0\in {\cal F}$. 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 the mentioned proposition. 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=12cm} \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=12cm} \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.7. Note that such an edge detection theory depends on no fewer than two parameters
that must be fixed by the user: smoothing scale and gradient threshold . }
\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=12cm} \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=12cm} \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 local 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 \ref{Chap9}. }
\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.  We display the level lines of the image at
levels multiples of 16. }
\end{figure}


\subsubsection*{Contrast invariance\label{section2.3}}

The use of contrast-invariant operators can solve some of the technical problems associated with linear smoothing and
other linear image operators. 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*}

\begin{exercise} Prove this last relation. Prove that if $g(s)=as+b$ then $g(u)$ satisfies the heat equation if $u$ does.
\finex
\end{exercise}

\section{Exercises}

\begin{exercise}\label{calculsHaralick} Define  an edge point $\fxx$
in a smooth image $u$ as a point $\fxx$ at which $g(t)$ attains a maximum,  where
$$g(t)=|Du\left(\fxx+t\frac{Du(\fxx)}{|Du(\fxx)|}\right)|.$$
 Prove by differentiating $g(t)$ that edge points satisfy $D^2u(\fxx)(Du(\fxx),Du(\fxx))=0$   \finex
\end{exercise}
\begin{exercise}
Construct simple functions $u$, $g$, and $k$ such that $g(k*u)\not=k*(g(u))$. \finex
\end{exercise}

\begin{exercise}
Consider the 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}$. Prove that 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}.
\]
Interpret the local behavior of the equation as a heat equation or a reverse heat equation
  according to the
size of $|Du|$ compared to $\lambda^{-1}$. \finex
 \end{exercise}

\section{Comments and references}


\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*{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 extensively 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}.


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