\newcommand \D[1] {\left[ #1 \right]}
\newcommand \deriva[2]{\frac{\partial #1}{\partial #2}}
\newcommand \bigh{\frac{}{}\!\!}

\chapter[Finite Difference Schemes]{Finite Difference Schemes for Curvature Motions}
We shall consider \label{PAGEORIGINAL}
 the classical discrete representation of an image $u$ on a grid
$u_{i,j}=u(i,j)$, with $1\le i\le N$, $1\le j\le N$ . The image is the union of the squares centered at the points
$(i,j)$.

\section{ Case of Mean curvature motion.}
We start with the ``Mean curvature motion" equation (M.C.M.) given by
$$\frac{\partial u}{\partial t} = \vert D  u\vert curv(u) =
\frac{u_y^2u_{xx}-2u_xu_yu_{xy}+u_x^2u_{yy}}{u_x^2+u_y^2} $$

In order to discretize this equation by finite differences we shall introduce an explicit scheme which uses a fixed
  $3\times 3$ stencil to discretize the differential operators. We denote by $\Delta x=\Delta y$ the pixel
width. From the PDE viewpoint $\Delta x$ is considered as  an infinitesimal length with respect to the image scale.
Thus we shall write formulas containing $o(\Delta x)$. Numerically $\Delta x$ is equal to 1, and the image scale ranges
from 512 to 4096 and more. By the order 1 Taylor formula one can give the following discrete versions of  the first
derivatives $u_x$ and $u_y$ at a point $(i,j)$ of the grid:

$$(u_x)_{i,j} = \frac{2(u_{i+1,j}-u_{i-1,j})+u_{i+1,j+1}-u_{i-
1,j+1}+u_{i+1,j-1}-u_{i-1,j-1}}{8\Delta x} +O(\Delta x^2) ; $$
$$(u_y)_{i,j} = \frac{2(u_{i,j+1}-u_{i,j-1})+u_{i+1,j+1}-u_{i+1,j-
1}+u_{i-1,j+1}-u_{i-1,j-1}}{8\Delta x} +O(\Delta x^2) ;$$
$$ | D  u_{i,j} | = ((u_x)_{i,j}^2+ (u_y)_{i,j}^2)^\frac{1}{2}.$$

\begin{definition}\label{defconsistency}
A discrete scheme approximating a differential operator is said to be consistent if, when the grid mesh $\Delta x$
tends to zero, the discrete scheme tends to the  differential operator.
\end{definition}
Clearly the above discrete versions of the partial derivatives and of the gradient of $u$ are consistent.
 When $| D  u|
\neq 0$, we can denote by $\xi$ the direction orthogonal to the gradient of $u$. It is easily deduced from Definition
\ref{definition4.14} that
$$\vert D u\vert curv(u)=u_{\xi\xi}.$$

\begin{exercise}
Show this formula. \finex
\end{exercise}
Defining $\theta$ as the angle between the $x$ direction and the gradient, we have
$$ \xi = (-\sin \theta , \cos \theta) =
(\frac{-u_y}{\sqrt{u_x^2+u_y^2}},\frac{u_x} {\sqrt{u_x^2+u_y^2}}), \hbox{ and}$$
\begin{equation} u_{\xi\xi} = \sin^2(\theta) u_{xx} - 2
\sin(\theta)\cos(\theta) u_{xy} + \cos^2(\theta)u_{yy}. \label{COS}
\end{equation}
We would like to write $u_{\xi\xi}$ as a linear combination of the values of $u$ on the fixed $3\times 3$ stencil. Of
course, the coefficients of the linear combination will depend on $\xi$.  Since the direction of  $\xi$ is defined
modulo $\pi$, we must assume by symmetry that the coefficients of points symmetrical with respect to the central point
of the stencil are equal (see Figure \ref{STENCIL}.)


\begin{figure}
\centerline{\hbox{ \epsfxsize=3.5cm \epsfbox{FDIR/stencil2.ps} }} \caption[Fig : A $3\times 3$ stencil]{A $3\times 3$
stencil} \label{STENCIL}
\end{figure}


 In order to ensure consistency
with the differential operator $u_{\xi\xi}$, we must find $\lambda_0 , \lambda_1 , \lambda_2 , \lambda_3, \lambda_4$,
such that
$$
(u_{\xi\xi})_{i,j} = \frac{1}{\Delta x^2} ( -4 \lambda_0 u_{i,j} + \lambda_1 (u_{i+1,j} +u_{i-1,j}) + \lambda_2
(u_{i,j+1} +u_{i,j-1}) $$  \begin{equation} \label{KKK}
    \qquad + \lambda_3 (u_{i-1,j-1} +u_{i+1,j+1}) +
    \lambda_4 (u_{i-1,j+1} +u_{i+1,j-1}))
+\varepsilon(\Delta x). \end{equation}
 We write
$$ u_{i+1,j} = u_{i,j} + \Delta x (u_x)_{i,j} +\frac{\Delta
x^2}{2} (u_{xx})_{i,j} + o((\Delta x)^3) ,$$ and the corresponding relations for the other points of the stencil. By
substituting  these relations into (\ref{KKK}) and by using  \eqref{COS} one obtains four links between the five
coefficients, namely

\begin{equation}
\label{LAMBDAS} \left\{
\begin{array}{l}
 \lambda_1(\theta) = 2 \lambda_0(\theta) - \sin^2\theta \\
 \lambda_2(\theta) = 2 \lambda_0(\theta) - \cos^2\theta \\
 \lambda_3(\theta) = - \lambda_0(\theta) + 0.5(\sin \theta \cos
\theta + 1) \\
 \lambda_4(\theta) = - \lambda_0(\theta) + 0.5(-\sin \theta \cos
\theta + 1) \end{array} \right. \end{equation}

\begin{exercise} Prove these four relations.
\finex \end{exercise} Thus, one degree of freedom is left for our coefficients : we can for example choose
$\lambda_0(\theta)$ as we wish. This choice will be driven by stability and geometric invariance requirements. Denoting
by $u^n_{i,j}$ an approximation of  $u(i\Delta x,j\Delta x,n\Delta t)$ we can write our explicit scheme as

\begin{equation} \label{MEANSCHEME} u^{n+1}_{i,j} = u^n_{i,j} +
\Delta t(u^n_{\xi\xi})_{i,j} \end{equation} \noindent Notice that this scheme can be rewritten as $ u^{n+1}_{i,j} =
\sum _{k,l=-1}^{1}\alpha_{k,l}u^n_{i+k,j+l} $ where the $\alpha_{k,l}$ satisfy $\sum_{k,l=-1}^{1}\alpha_{k,l} =1.$ The
following obvious lemma shows a general condition to have
{\it $L^\infty$ stability} in this kind of scheme. %\newtheorem{lemma}{lemma}
\begin{lemma}
  Let a finite difference scheme given by
  $$ T(u)_{i,j} = \sum _{k,l=-1}^{1}\alpha_{k,l}u_{i+k,j+l} $$
  where $\alpha_{k,l}$ satisfy $\sum_{k,l=-1}^{1}\alpha_{k,l} =1 $.
  We say that the scheme is {\rm $L^\infty$-stable} if for all $i,j$, $$\min_{i,j}u(i,j)\leq  T(u)_{i,j}\leq
  \max_{i,j}u(i,j).$$
Then the scheme is
  $L^\infty$ stable if and only if $\alpha_{k,l}\geq0$ for any $k,l$.
\end{lemma}

\begin{proof}
  If $\alpha_{k,l}\geq0$ for any $k,l$, set $min~=~inf_{i,j}\{
   u_{i,j} \}$, $max=sup_{i,j}\{ u_{i,j} \}$ and take a
point $(i,j)$.  Then  the  $L^\infty$ stability follows from the inequality:
   $$min = \sum _{k,l=-1}^{1}\alpha_{k,l}min \leq \sum _{k,l=-
1}^{1}\alpha_{k,l}u_{i+k,j+l} =
   (Tu)_{i,j} \leq \sum _{k,l=-1}^{1}\alpha_{k,l}max = max$$
   On the other hand, if there exists  $\alpha_{k_0,l_0}<0$
then choosing $u$ and $(i,j)$ such that $u_{i+k_0,j+l_0} = min$ and  $u_{i+k,j+l} = max$ for any other $k,l$,
   we obtain
   $$(Tu)_{i,j} = \sum_{k\neq k_0,l\neq l_0}^{1}\alpha_{k,l}max
+ \alpha_{k_0,l_0}min =
   max+\alpha_{k_0,l_0}(min-max) > max,$$
   which means that the $L^\infty$ stability is violated.

\end{proof}

Following this lemma, in order to guarantee the $L^\infty$ stability in the scheme \eqref{MEANSCHEME} we should look
for $\lambda_0$ such that $\lambda_1, \lambda_2,\lambda_3,\lambda_4\geq 0$ and $(1-\frac{4\lambda_0}{\Delta x^2})\geq
0$. Unfortunately the links between these coefficients make it impossible to obtain these relations, except for the
particular values of $\theta =( 0,\frac{\pi}{4}, \frac{\pi}{2},...)$. Indeed, for $\theta$ in $[0,\frac{\pi}{4}]$,
$$ \lambda_1 \geq \lambda_2 \qquad \hbox{and} \qquad \lambda_3 \geq \lambda_4 $$ But
$$ \lambda_2(\theta) \ge 0 \Rightarrow \lambda_0(\theta) \geq
\frac{\cos^2(\theta)}{2} $$
$$ \lambda_4(\theta) \ge 0 \Rightarrow \lambda_0(\theta) \leq
\frac{1 -\sin(\theta)\cos(\theta)}{2}$$
 We cannot find $\lambda_0(\theta)$ satisfying both inequalities, since
$$ \frac{\cos^2(\theta)}{2} \geq \frac{1
-\sin(\theta)\cos(\theta)}{2}$$ If we chose $\lambda_0(\theta) \geq \frac{\cos^2(\theta)}{2}$, $\lambda_4(\theta)$
would be significantly below zero. If we took $\lambda_0(\theta) \leq \frac{1 -\sin(\theta)\cos(\theta)}{2}$,
 $\lambda_2(\theta)$ would be significantly below  zero. Thus we  shall  choose $\lambda_0$ somewhere between
 both functions, so that
$\lambda_2$ and $\lambda_4$ become only slightly negative. (see Figure \ref{LAMBDA0}.)

\vspace{2mm}

\begin{figure}
\begin{center}
\input{FDIR/COURBE.tex}
\caption[Fig : Arbitrary choice for discretizing the MCM]{The middle curve represents the choice of the function
$\lambda_0$ of Formula \ref{CHOICE2}. The upper function represents the smallest possibility for $\lambda_0(\theta)$
securing $\lambda_2 \geq 0$ for all angles and the lower one represents the largest values of $\lambda_0(\theta)$
securing $\lambda_4(\theta) \geq 0$. Thus, it is not possible to satisfy simultaneously both conditions. The
intermediate curve is the simplest trigonometric function which lies between these two bounds.} \label{LAMBDA0}
\end{center}
\end{figure}

In addition, we can try to impose on $\lambda_0$ the following geometrical requirements
\begin{enumerate}
 \renewcommand{\theenumi}{(\roman{enumi})}
  \item Invariance by rotation of angle $\frac{\pi}{2}$
    $$ \lambda_0(\theta +\frac{\pi}{2}) = \lambda_0(\theta)$$
  \item Purely one-dimensional diffusion in the case $\theta = 0,\frac{\pi}{2},...$
    $$ \lambda_0(0) = 0.5 $$
    This condition implies that
$\lambda_2(0)=\lambda_3(0)=\lambda_4(0)=0$\\
  \item Pure one-dimensional diffusion in the case $\theta = \frac{\pi}{4},
\frac{3\pi}{4},...$
    $$ \lambda_0(\frac{\pi}{4}) = 0.25 $$
    This condition implies that
$\lambda_1(\frac{\pi}{4})=\lambda_2(\frac{\pi}{4})
    =\lambda_4(\frac{\pi}{4})=0$\\
  \item  Symmetry with respect to the axes i+j and i-j,
   $$\lambda_0(\frac{\pi}{2}-\theta) = \lambda_0(\theta)$$
\end{enumerate}
We remark that by the above conditions it is enough to define the function  $\lambda_0(\theta)$ in the interval
$[0,\frac{\pi}{4}]$ because it can be extended  by periodicity elsewhere.

Two choices for the  function $\lambda_0(\theta)$ using as basis the trigonometric polynomials were tested. The first
one corresponds to an average of the boundary functions:

\begin{equation} \label{CHOICE1} \lambda_0(\theta) = \frac{\cos^2(\theta) +1
-\sin(\theta)cos(\theta)}{4} \end{equation} As we shall see this choice is well-adapted to the ``affine curvature
motion'' equation. However,  if we extend this function by periodicity, the extended function is not smooth at
$\frac{\pi}{4}$. If we seek for a smooth function for $\lambda_0(\theta)$, we must impose
$\lambda_0'(0)=\lambda_0'(\frac{\pi}{4})=0$. The  trigonometric polynomial with least degree satisfying the above
conditions and lying between both boundary functions is


\begin{equation} \label{CHOICE2} \lambda_0(\theta)) = 0.5 -
\cos^2(\theta)\sin^2(\theta) \end{equation} The formulas of the other $\lambda_i$'s are deduced using (\ref{LAMBDAS}).
For instance with the above choice of $\lambda_0(\theta)$ we have

$$
\left\{
\begin{array}{l}
\lambda_1(\theta) = \cos^2(\theta) (\cos^2(\theta) - \sin^2(\theta));
\\
\lambda_2(\theta) = \sin^2(\theta) (\sin^2(\theta) - \cos^2(\theta));
\\
\lambda_3(\theta) = \cos^2(\theta)\sin^2(\theta) + 0.5
\sin(\theta)\cos(\theta); \\
\lambda_4(\theta) = \cos^2(\theta)\sin^2(\theta) - 0.5
\sin(\theta)\cos(\theta). \\
\end{array}
\right.
$$

\vspace{5mm} When $| D  u| = 0$, the direction of the gradient is unknown. Therefore the diffusion term $u_{\xi\xi}$ is
not defined. We chose to replace this term by half the Laplacian. (The Laplacian is equal to the sum of the two second
derivatives in orthogonal directions, whereas the diffusion term $u_{\xi\xi}$ is the second derivative in just one).
However, other possibilities will be considered in Section \ref{FDSgradientnull}.
 Summarizing, a consistent, almost $L^\infty$ stable finite difference scheme for the mean curvature motion is
  (iterations start with $u^0$ as initial function)

\begin{enumerate} \item If $ \vert D  u\vert \geq T_g $
$$ u^{n+1} = u^{n} + \frac{\Delta t}{\Delta x^2} ( -4
\lambda_0 u_{i,j} + \lambda_1 (u_{i+1,j} +u_{i-1,j}) + \lambda_2 (u_{i,j+1} +u_{i,j-1}) $$
$$   \qquad + \lambda_3 (u_{i-1,j-1} +u_{i+1,j+1}) +\lambda_4 (u_{i-1,j+1}
+u_{i+1,j-1})). $$

\item Otherwise,$$ u^{n+1} = u^{n} + \frac{1}{2}\frac{\Delta t}{\Delta x^2} (  -4 \lambda_0 u_{i,j} + u_{i+1,j}
+u_{i-1,j} + u_{i,j+1} +u_{i,j-1})
$$

\end{enumerate}
Two parameters have to be fixed in the previous algorithm:

\paragraph{$\bullet$ The iteration step scale} $s:=\frac{\Delta t}{\Delta x^2}$ has to be chosen as
large as possible in order to reduce the number of iterations. However, $ 1/2$ is a natural upper bound for $s$.
Indeed, consider the discrete image defined by $u^0_{i,j} = 0$ for all $i,j$, except for $i=j=0$ where $u^0_{0,0}=1$.
Then the second formula yields $u^1_{0,0} = 1 - 2 * s$. If we want $L^\infty$ stability to be ensured we must have
 $u^1(0,0)\geq 0$,  which yields $s \leq 1/2$. Imposing this condition
\begin{equation} \label{stepfixed} \frac{\Delta t}{\Delta x^2} \leq \frac{1}{2}\end{equation}
it is an experimental observation that there is  a (small with respect to 255) $\epsilon>0$ such that  for any $n \in
\natural$ and $(i,j)$,
$$ -\epsilon+\inf_{i,j}\{u_{i,j}^0\} \leq u_{i,j}^n \leq
\sup_{i,j}\{u_{i,j}^0\} + \epsilon.$$

\paragraph{$\bullet$ The threshold on the spatial gradient norm} : $T_g$ has been fixed experimentally to $6$ for $0$ to $255$ images.

\begin{figure}[!ht]
\centerline{\hbox{ \epsfxsize=11cm \epsfbox{DDIR/plancheAMSS.ps} }} \caption[Curvature motion finite difference scheme
and scale calibration]{  \label {curvaturemotiondifferencescheme} Curvature motion finite difference scheme and scale
calibration. Image filtered by curvature motion at scales 1, 2, 3, 4, 5. In order to give a sound numerical meaning to
the scale, a calibration of the numerical scales (number of iterations) is made in such a way that a disk with radius
$t$ shrinks to a point at scale $t$. }
\end{figure}

\begin{figure}[!ht] \centerline{\hbox{ \epsfxsize=11cm \epsfbox{DDIR/plancheAMSS.ps} }} \caption[Curvature
motion finite difference scheme applied on each level set separately (threshold superposition principle)]{ Curvature
motion finite difference scheme applied on each level set separately, at scales 1, 2, 3, 4, 5. The processed image is
then reconstructed by the threshold superposition principle.  In contrast with the same scheme directly applied on the
image, this scheme yields a fully contrast invariant smoothing. However, a comparison with  Figure
\ref{curvaturemotiondifferencescheme} shows that the resulting images are very close to each  other. This shows that
the contrast invariance is almost achieved when applying the finite difference scheme directly on a good quality image.
    As we shall
see in Figure \ref {figure_amss_extrema}, if the initial image is noisy, the difference between both methods can be
huge. \label{curvaturesuperposition}}
\end{figure}

\begin{figure}[!ht] \centerline{\hbox{ \epsfxsize=11cm \epsfbox{DDIR/plancheMED.ps} }} \caption[Curvature motion by
median filtering]{  Iterated median filter with  normalized  scales 1, 2, 3, 4, 5. The scale normalization permits to
compare very different schemes on the same images. Compare with Figure \ref{curvaturesuperposition}. The striking
 similarity of the results anticipates Theorem \ref{theorem10.4}, according to which the application of the median filter is
 equivalent to a mean curvature motion.}\end{figure}



\section{FDS for AMSS}
We will use the ideas developed  in the above section. We rewrite the AMSS equation as

\begin{equation} \label{AFIN} \frac{\partial u}{\partial t} = (\vert D  u\vert^3
curv(u))^{\frac{1}{3}} = (u_y^2u_{xx}-2u_xu_yu_{xy}+u_x^2u_{yy})^\frac{1}{3} \end{equation} We remark that $\vert D
u\vert^3 curv(u) = \vert  D  u\vert^2u_{\xi\xi}$ where $\xi$ corresponds to the direction orthogonal to the gradient.
Therefore, in order to discretize this operator, it is enough to multiply the discretization of $u_{\xi\xi}$ presented
in the above section by $\vert  D  u\vert^2$.  We choose $\lambda_0(\theta)$ given by (\ref{CHOICE1}) because it
corresponds to a trigonometric polynomial of degree two and then multiplying it by $\vert  D  u\vert^2$ the
coefficients $\eta_i=\vert  D  u\vert^2\lambda_i,  i=0,1,2,3,4$, are polynomials of degree two with respect to $u_x$
and $u_y$.  Indeed, we obtain for $ \theta \in [0,\frac{\pi}{4}]$

$$ \label{KKKA}
(\vert  D  u\vert^2u_{\xi\xi})_{i,j} = \frac{1}{\Delta x^2} ( -4 \eta_0 u_{i,j} + \eta_1 (u_{i+1,j} +u_{i-1,j}) +
\eta_2 (u_{i,j+1} +u_{i,j-1})
$$
$$
    \qquad + \eta_3 (u_{i-1,j-1} +u_{i+1,j+1}) +\eta_4 (u_{i-1,j+1}
+u_{i+1,j-1})) +O(\Delta x^2)
$$

where $\eta_0,\eta_1,\eta_2,\eta_3,\eta_4$ are given by

$$
\left\{
\begin{array}{l}
\eta_0 = 0.25( 2u_x^2 + u_y^2-u_xu_y) \\
\eta_1 = 0.5( 2u_x^2 - u_y^2-u_xu_y) \\
\eta_2 = 0.5( u_y^2-u_xu_y) \\
\eta_3 = 0.25( u_y^2+3u_xu_y) \\
\eta_4 = 0.25( u_y^2-u_xu_y) \\
\end{array}
\right.
$$

Finally, the finite difference scheme for the A.M.S.S. equation is

\begin{equation}
\label{AFINSCHEME} u^{n+1}_{i,j} = u^n_{i,j} + \Delta t(\vert  D  u^n \vert ^2u^n_{\xi\xi})_{i,j}^{\frac{1}{3}}
\end{equation}
We have tested this algorithm and we have noticed that in this case the condition for the experimental stability (in
the sense presented in the above subsection) is

$$\frac{\Delta t}{\Delta x^2} \leq \frac{1}{10}.$$

{\bf Remark.} The finite difference schemes  presented above are consistent.
 Contrast invariance can only be  obtained asymptotically by
taking a very small time step $\Delta t$. The
 experimental results presented in Figures \ref{curvaturemotiondifferencescheme} and
\ref{AFINIMAGE} have been obtained by using these schemes with
 $\Delta t = 0.1$ in the case of mean curvature
motion and $\Delta t = 0.01$ in
 the case of affine curvature motion. Indeed, while experimental stability is achieved with $\Delta t \leq 0.1$,
  the experimental affine
invariance  needs $\Delta t < 0.05$ (see Figure \ref{figure_amss_affine_char}.)


\section{\label{Ilmanque}IL MANQUE UNE EXPERIENCE AMSS SUR L'INVARIANCE AFFINE!}

\section{Numerical normalization of scale.}
(or Relation between scale and the number of iterations).


\paragraph{The case of the curvature motion.}
Setting the distance between pixels $\Delta x$ to 1, the scale achieved with $N$ iterations is simply $N\times\Delta
t$. Now, the scale $t$ associated with  the PDE is somewhat arbitrary : It has no geometric meaning. In order to get
it, we need a rescaling $T\to t(T)$ which we will call scale normalization.

A good way to perform this scale normalization is to define the correspondence  $t(T)$ as the time for which a circle
with initial radius $T$ vanishes under curvature motion. Such a circle moves at a speed equal to its curvature, which
is the inverse of its radius. Thus have for a disk with radius $R(t)$
$$\frac{dR(t)}{dt} = - \frac{1}{R(t)}$$
which yields
$$ \frac{1}{2} (R^2(0) - R^2(t)) = t.$$
\begin{exercise} Check this relation!
\finex \end{exercise}
 The disk disappears  when $R(t)= 0$, that is, at scale $T=R^2(0)/2.$ This last relation gives a
scale normalization: In order to arrive at  the normalized scale $T$  (at which any disk with radius less or equal to
$T$ vanishes), we have to evolve the PDE at  $t=N\Delta t = T^2/2$. This fixes  the number of needed iterations as
$$ N= T^2/{2 \Delta t}.$$

\paragraph{The case of AMSS}
We can perform similar calculations. The radius of an evolving disk satisfies
$$\frac{dR(t)}{dt} = - \frac{1}{R(t)^\frac{1}{3}}$$
which yields
$$ \frac{3}{4} (R^\frac{4}{3}(0) - R^\frac{4}{3}(t)) = t$$
The disappearance time is therefore $t= \frac 3 4 R^\frac{4}{3}$. As for the curvature motion, we define the normalized
scale $T$ as the one at  which a disk with radius $T$ vanishes. In order to achieve this scale $T$, the needed number
of iterations is
$$ N = \frac 3 {4 \Delta t} T^\frac{4}{3}. $$
\begin{exercise}
Check the last two formulas! \finex
\end{exercise}


\section{Contrast invariance and the level set extension}

Both schemes (M.C.M and A.M.S.S) presented above are not numerically contrast invariant. We have seen that a contrast
operator cannot create new gray levels (Exercise \ref{necreepasdenouveauxniveaux}.)
 Now, starting with a binary image
$u^0$ and applying a scheme defined by such a formula as
$$ u^{n+1} = u^{n} + \Delta t (...)$$
does not ensure that $u^{n+1}$ will be also a binary image.

A natural idea to overcome this problem is the following. Starting with a binary image (with values 0 and 1): apply the
scheme until the expected scale is achieved, then threshold the obtained image at $\lambda = \frac 12$. This of course
works only for binary images.  However, the level set extension (see Section \ref{section6.3}) gives us the key to
extend this to general images.

The contrast invariance can be fully obtained by first applying the finite difference scheme on each level set
(considered as a binary image) separately. Then by the superposition principle the evolved image is computed from the
evolved level sets. The procedure is  the following~:


\paragraph{Algorithm starting with an image $u_0$ and evolving it to $u(t,\fxx)$ by curvature motion}
\begin{description}

\item{For each $\lambda\in [0, 255]$, in increasing order:}
\begin{itemize}
\item Let $v_\lambda(\fxx)$ be the characteristic function of $\levelset_\lambda u_0$. (This function is equal to $1$
 inside the level set and to $0$ outside.)
\item Apply to $v_\lambda$ the MCM or AMSS FDS-scheme until scale $t$. This yields the images $w_\lambda(t,.)$.
\item Set $u(t,\fxx)=\lambda$ at each point $(t,\fxx)$ where $ w_\lambda(t,\fxx) \geq 0.5.$
\end{itemize}
\end{description}



\section{Problems at extrema \label{FDSgradientnull}}
 For  MCM and AMSS we raised  the question  of  performing numerically the
equation when $| D  u|=0$.
    For MCM the right hand part of the equation is simply not defined. For
AMSS one can set by continuity as $ D  u\to 0$, $(\vert  D  u^n \vert ^2u^n_{\xi\xi})_{i,j}^{\frac{1}{3}}=0$.
    Now, numerically, this would imply that isolated black or white extrema will not evolve by the equation.  We know that
    this is simply wrong, since small sets collapse by curvature motion.

In short, FDS for  MCM and AMSS are not consistent with the equation at extrema. In  Figure \ref{figure_amss_extrema},
we added to an image a strong ``salt and pepper'' noise. More than one fourth of the pixels have been given a uniform
random value in $[0,255]$ and most of them have become local extrema. Not only these values do not evolve but they
contaminate their neighboring pixels. There are easy  ways to avoid this spurious effect~:

\begin{itemize}
\item One can first zoom by 2 the image by duplicating pixels. This, however, multiplies by 16 the number of computations.
\item One can first remove  pixels extrema with diameter $k$ since they  must anyway  disappear by the equation at
normalized scale $\frac k2$.
\item One can use the level set method. This multiplies  the number of computations by the initial number of gray-levels.
\end{itemize}
All of these solutions are efficient, as shown in Figure \ref{figure_amss_extrema}.



\begin{figure}
\begin{center}
 \begin{tabular}{cc}
\epsfxsize=4cm \epsfbox{DDIR/BRUIT40.ps}  & \epsfxsize=4cm
\epsfbox{DDIR/AMSS3.ps} \\
\epsfxsize=4cm \epsfbox{DDIR/ZAMSS2.ps} & \epsfxsize=4cm
\epsfbox{DDIR/ZAMSS3.ps} \\
\epsfxsize=4cm \epsfbox{DDIR/KAMSS2.ps} & \epsfxsize=4cm
\epsfbox{DDIR/KAMSS3.ps} \\
\epsfxsize=4cm \epsfbox{DDIR/OSAMSS2.ps} & \epsfxsize=4cm
\epsfbox{DDIR/OSAMSS3.ps} \\
\end{tabular}
\end{center}
\caption[Various implementations of the curvature motion on a noisy image]{  \label{figure_amss_extrema} Various
implementations of curvature motion on a noisy image. Top left : image with 40$\%$ pixels replaced by a uniform random
value in $[0,255]$. Top right: application of the finite difference scheme (FDS) at normalized scale 3. On the lines 2
to 4, we see various solutions to the disastrous diffusion of extrema. On the left the image is processed at normalized
scale 1 and on the right at normalized scale 3.  Second line:  FDS applied on the image previously zoomed by a factor
2; third line: FDS applied on the image after its extrema have been "killed" (the reference area is given by the area
of the disk vanishing at the desired scale). Fourth line:  FDS applied separately on each level set and application of
the threshold superposition principle. The third scheme offers a good speed-quality compromise.}
\end{figure}


%\begin{figure}  IMAGE ORIG N'EST PAS LA BONNE!
%\begin{center}
%\epsfxsize=7.8cm \epsfbox{DDIR/ORIG.ps} \epsfxsize=3.8cm \epsfbox{DDIR/O3.ps}
%\end{center}
%\caption[Finite difference scheme for AMSS]{  \label{figure_amss_affine_char} Finite difference scheme for AMSS. Left :
%Original noisy image of letters, middle : same image without the extrema with area less than 80 pixels (extrema
%killer).  Right : result of AMSS at calibrated scale 4. }\end{figure}

\section{Conclusion}

We have seen that standard  finite difference schemes are easy to implement but cannot handle properly the invariance
properties satisfied by the equations.


\begin{enumerate}
\item There is no finite difference scheme that insures the monotonicity. This leads to slightly
oscilatory solutions.

\item No full contrast invariance. For instance FDS create new grey levels and blur  edges. Also,
a spurious diffusion  occurs around the image extrema. However this last problem was dealt with efficiently in the
previous section. The full contrast invariance has been restored by the level set extension of the numerical schemes.

\item
The worst drawback of FDS is the lack of Euclidean  or affine invariance which can be only approximately obtained by
grid local schemes. A much more clever strategy to achieve full invariance is to evolve all level curves of the image
and the reconstruct it. This is the aim of Section \ref{sectionlionel}, but we have already seen in Chapter
\ref{chapfourways} how to evolve curves by curvature.

\end{enumerate}

\section{Comments and references}
\medskip
\noindent {\bf Difference schemes for the curvature motion and the AMSS} The presented difference scheme follows mainly
\cite{Guichard:1994:AAM}, improved in
 Alvarez et al. \cite{Alvarez:1999:FCA}.
This scheme is somehow optimal among the rotationally invariant numerical schemes for curvature motion and the AMSS.
Now, this presentation is specific of those two motions, while other many authors have analysed more general nonlinear
anisotropic diffusions in image processing, namely Acton \cite{bb2091}, Kacur and Mikula \cite{Kacur:1995:SND,
kacur1997sad}. Weickert and the \"{U}trecht school \cite{Niessen:1994:NAG, Weickert:1996:NDS, Weickert:1997:PIA,
Weickert:1998:ERSa} address many aspects of implementation of nonlinear scale spaces, namely speed, parallelism and
robustness. Crandall and Lions \cite{Crandall:1996:CDS} also proposed a finite difference scheme for mean curvature
motion, valid in any dimension. Sethian's book \cite{Sethian:1996:LSM} explains how to implement fast the motion of a
curve or surface by the so called "level set method", where a distance function to the curve or surface is evolved.
Dynamic programming allows a fast implementation (the "fast marching method").


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