
\chapter{Numerical Implementations. \label{chapter22}}
\thispagestyle{empty}

\section{Digital level sets and curves.}

Depending on the expected precision of the level set or level curve representation, there exists different
strategies concerning their extraction. We address here the simplest, by considering a pixel precision. For a higher precision, a model of interpolation between  pixel values can be chosen. Candidates of such models stand e.g. in Shannon interpolaton, Lipschitz interpolation, bi-linear interpolation, etc...

We consider the digital image as a array which associates to each pixel a gray level value.
 
{\bf The level sets.}
The upper or lower level set of level $\lambda$ can be simply obtained by thresholding the image at the level $\lambda$. The upper level set is the set of the pixels that have a gray level higher or equal to $\lambda$. (And conversely for the lower).

{\bf The Connected components of level sets}
To extract the connected components of the level sets, we have first to define what is meant by connected. As well known, there exist two kinds of connectivity in a square-based grid : The 4-connectivity and the 8-connectivity (see figure \ref{figureconnectivity}).


\begin{figure}
\begin{center}
\epsfxsize=6cm
\epsfbox{\FDIR/connectivity.eps}
\end{center}
\caption[Connectivity on a grid]{Connectivity on a grid. Left, 4-connectivity that is each pixel is connected to the ones directly up, left, right and down. Right, 8-connectivity. \label{figureconnectivity}}
\end{figure}

One can arbitrary choose one or the other, being aware that this choice violates the symmetry $u \rightarrow -u$.

Choosing the 8-connectivity for the upper level sets, means in fact considering the digital image as an
upper semi-continuous function. This function is constructed by putting the gray level value of each pixel to all the point enclosed in it and by setting the pixel boundaries to the sup of the adjacent pixels values.
Similarly, choosing the 4-connectivity for the upper level sets, yields to a lower semi-continuous representation of the image.

Once the choice of connectivity is done, extracting the connected composants of the level sets is straitghforward.

%%{\bf The level lines representation.}
%%The level lines are defined as the connected components of the boundaries of upper the level sets. 

%%Now, since boundaries of pixel sets are passing through the pixel borders, it is not easy to represent them using the pixel grid. Take example of a single pixel set. How to represent its boundary with pixels~: the pixel itself or its four neighbors ?


%%One simple way to represent the level curves with pixels and without ambiguity is to first zoom in the image with a factor 2. That is to duplicate each pixel in the X and then in the Y direction.
%%As shown in Figure ??, the single pixel set looks surrounded by the level curve drawn itself with pixels. 
{\bf The level lines extraction.}
The level lines are defined as the connected components of the boundaries of upper the level sets. 
Assuming that the connected components of the level sets have been extracted, one has then just to follow their boundaries.

Note that the boundaries of connected components of a level sets are of two kinds. There is in one hand the exterior envelop, in the other hand there are the boundaries that envelop the holes of the level set (if it has some).

Assume that we follow one of boundary (exterior or a hole) of the set, in a way so that the set is located on the left when facing the path direction. The path can be followed without
ambiguity until a junction. In case of 4-connectivity (resp. 8), the most left (resp. right) direction of the junction should be chosen.

\section{Median Filter}

{\bf General median filter}

Many different digital implementations of the median filter have been proposed. If we forget problem of time computation, the simplest implementation is as follow:

Given a mask (support of the measure) with a odd number of pixels: $2N+1$,  the median value at pixel $p$ is the ``$N$'' largest value among the pixel in the mask centered at $p$. Note that partial sorting, using e.g. quick-sort, can provide such a value. More clever ways to extract the median value, with a limited number of comparaisons, have been proposed depending on the chosen mask.

{\bf median filter on a disk}

The problem is to approximate adequatly a disk by a mask made of pixels. One strategy is to use a weighted median filter. We define the weight of a pixel as the area of the intersection of the pixel itseld and the disk. Then the median
value is the smallest value such that the sum of the weight of the pixels with smaller value passes 50\% of the total weight.

\PDYNAMSH

In case of binary images, (or shapes), as remarked in Chapters \ref{chapter1} and \ref{chapter2}, the weigthed median filter is equivalent to a linear convolution with the weighted masked followed by a threshold. Figure \ref{dynamicshapeexperiment} illustrated this fact with a Gaussian weighted median filter.


{\bf Implementing the median filter using the mean curvature motion...}

In any cases, the median filter suffers two drawbacks:

First, it is too local to catch small curvature. Indeed, it is impossible to make the distinction between a straight line and a slightly curved line if one looks locally and with a pixel precision.

Second, it generates quantized evolutions due to the quantization of space (pixels) and of the gray levels.

As we have seen, iterated median filter converges towards viscosity solution of
the curvature motion. One can therefore consider to use the curvature motion equation as an implementation of an ``ideal'' median filter.

\section{Extrema Killer : Maxima and Minima Killer.}

The ``Extrema Killer" has been formaly defined in section \ref{extremakillersection}. Its formal definition
involve an infinite set of structuring elements, which makes its formal definition not in practical use.
We have shown that it effect is just to remove the small (in area) connected components of the upper level sets for the maxima killer ( and lower level sets for the minima killer).

Therefore an idea to implement the maxima killer it is then to start at local maxima and to go down until the area is achieved. Going down means to add adjacent pixels to a stack starting a the local maximum. The stack represents then a region of pixels which could be interpreted as a local ``mountain" if the gray levels are interpreted as altitudes...
Now, problem occurs when the stack meets an saddle points. An saddle point can be where two ``mountains" meets. And the stack of the two ``mountains" should be added, which would mean that the stack corresponding to each local maxima should be computed in parallel.

A strategy to avoid that is then to stop the increase of the stack when it meets an saddle point, that is in other words, when an adjacent pixel to the stack has a larger value than the min value of the pixels in the stack. The pixels in the stack are then set to this min value, correponding to the value of the saddle point. As consequence the saddle point is no more a saddle point, and these pixels will continue to be treated later, see Figure \ref{extremakillerimplementation}

\begin{figure}
\begin{center}
\epsfxsize=14cm
\epsfbox{\FDIR/extrema.eps}
\end{center}
\caption[Extrema killer implementation]{\label{extremakillerimplementation} \Pstyle  Extrema killer implementation. This figure illustrates the implementation of the maxima killer for one dimensional image. (As described in the text, the algorithm is similar in dimension two). The size asked is displayed in bottom. The algorithm starts with a local maxima, a). And, it goes down in gray level until either the size is achieved or a saddle point is touched. In both case, the pixels are set to the achieved gray level value, b). In the shown case, the area is not achieved, but it will be so, when the next maxima will be treated, c) and d). The minima killer works symmetrically.}
\end{figure}


An implementation of the maxima killer is then as follow:
\begin{enumerate}
\item Choose a connectivity 4 or 8 and an area $A$ (number of pixels).
\item Start from a local maximum pixel. (The value of the pixel is larger or equal to its neighbors, and at least strictly larger than one of them).
Let $\lambda$ be its gray level.
Initialize a stack with this pixel.
\item Among the neighbor pixels of all pixels of the stack, choose the one that has the largest value. Add it to the stack.
\item If this pixel has a value less than $\lambda$, set $\lambda$ to its value, else goto 6
\item If the stack is smaller than $A$ goto 3.
\item Set all pixels of the stack to the value $\lambda$. Empty the stack.
\item If it remains a local maxima pixel not treated, find it and goto 2. 
\end{enumerate}

The symmetrical implementation will perform the minima killer. The extrema killer made of an maxima and a minima killer defines in fact two possible filters
depending on which of the maxima or the minima killer is done first.
A symmetrical version of the extrema killer can be found in \cite{Masnou:1998:IRI,MOGu}




\section{Finite Difference Scheme (FDS) for the Curvature Motion and AMSS.}
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)$, and
the brightness in each square is constant~: $u_{i,j}$. Each one
of the squares is called {\it pixel} (for ``picture element").

\subsection{ Case of Mean curvature motion.}
We start with the ``Mean curvature motion" equation (M.C.M.) given by 
$$\frac{\partial u}{\partial t} = \vert\nabla 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 stencil  of $3x3$ points to 
discretize the differential operators. 
For simplicity, we assume that the spatial
increment  $\Delta x$ is the same in the x-axis and the y-axis. We
approach the first derivatives $u_x$ and $u_y$ in a point 
$(i,j)$ of the lattice by using the following linear scheme:

$$(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}}{4\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}}{4\Delta x} +O(\Delta x^2)$$
$$ |\nabla u_{i,j} | = ((u_x)_{i,j}^2+ (u_y)_{i,j}^2)^\frac{1}{2}$$
When $|\nabla u| \neq 0$, we can denote by $\xi$ the direction orthogonal to the gradient of
$u$, one easily sees that $\vert\nabla u\vert curv(u)$ is equal to
$u_{\xi\xi}$. Defining $\theta$ as the angle between the horizontal axe 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 want to write $u_{\xi\xi}$ as a linear
combination of the values of $u$ on the fixed stencil 3$\ast$3.
Of course, the coefficients of the linear combination may depend
on $\xi$.  Because the direction of the gradient (and then
$\xi$) is defined modulo $\pi$, by symmetry we must assume that
the coefficients of   points symmetrical with respect to the
central point are the same (see figure \ref{STENCIL}). 

\Pstencil


In order to be consistant, 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}))
+O{\Delta x^2}
\fe&
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 same
relation for the other points of the stencil. By feeding
(\ref{KKK})  with these relations and by using relation
(\ref{COS}), we obtain four relations  between our five
coefficients

\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}
There remains  one degree of freedom for our coefficients given
by the choice of $\lambda_0(\theta)$. We shall choose 
$\lambda_0(\theta)$ following the stability and geometric
invariance criteria. 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
Note 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
{\bf L$^\infty$ stability} in this kind of schemes: %\newtheorem{lemma}{lemma}
\begin{lem}
  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 $. 
Then the scheme satisfies 
  {\bf L$^\infty$ stability} if and only if $\alpha_{k,l}\geq0$ for any $k,l$.
\end{lem}

\paragraph{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  {\bf 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$$
   And therefore {\bf L$^\infty$ stability} is violated.
   
\endproof

Following this lemma, in order to have {\bf L$^\infty$ stability} in the scheme
(\ref{MEANSCHEME}) we must seek 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, 
because of the relations between our coefficients, it is
impossible to obtain these relations,  except for particular
values of $\theta =( 0,\frac{\pi}{4}, \frac{\pi}{2},...)$. Indeed, We 
remark
that 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}$$

So, we cannot find $\lambda_0(\theta)$ satisfying both inequalities, 
since
$$ \frac{\cos^2(\theta)}{2} \geq \frac{1 
-\sin(\theta)\cos(\theta)}{2}$$
Then, if we choose $\lambda_0(\theta) \geq 
\frac{\cos^2(\theta)}{2}$ we have $\lambda_4(\theta)$ very 
negative. If we take $\lambda_0(\theta) \leq 
\frac{1 -\sin(\theta)\cos(\theta)}{2}$ we  obtain
$\lambda_2(\theta)$ very negative. We prefer to choose
$\lambda_0$ between both functions, and then to have
$\lambda_2$ and $\lambda_4$ negative, but slightly. (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 achieving for all angles $\lambda_2 \geq0$, the lower one represent the largest possibility for achieving $\lambda_4 \geq0$. Therefore, it is not possible to satisfy both conditions. We then have chosen the simplest trigonometric function which is in between these two constraints.}
\label{LAMBDA0}
\end{center}
\end{figure}

On the other hand, we impose to $\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 Pure 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 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.

We have tested two choices for the  function
$\lambda_0(\theta)$ using as basis the trigonometric
polynomials. 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 simplest
polynomial,  of degree as small as possible, satisfying the
above conditions, and between both boundary functions is


\begin{equation} \label{CHOICE2} \lambda_0(\theta)) = 0.5 -  
\cos^2(\theta)\sin^2(\theta) \end{equation}


We deduce the other $\lambda$'s values 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 $|\nabla u| = 0$, the direction of the gradient is not known. Therefore the diffusion term $u_{\xi\xi}$ is not defined. We choosed to replace this term by half the laplacian. (The laplacian corresponds to the sum of the two second derivative in orthogonal directions, whereas the diffusion term $u_{\xi\xi}$ is the second derivaty in one). However others possibilities will be considered in section \ref{FDSgradientnull}.
 So the FDS scheme for the mean curvature motion is, (iterations start with $u^0$ as initial function) 

Where $ \vert\nabla u\vert \leq 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 + \eta_3 (u_{i-1,j-1} +u_{i+1,j+1}) +\eta_4 (u_{i-1,j+1} 
+u_{i+1,j-1})) $$
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})
$$

\vspace{5mm}

\paragraph{The iteration step scale} ($\frac{\Delta t}{\Delta x^2}$) has to be chosen as
large as possible in order to reduce the number of iterations. However, there is a natural upper bound of $\frac{1}{2}$. Indeed, denoting $s$ this step, and
considering the following image $u^0_{i,j} = 0$ for all $i,j$, except for $i=j=0$ where we set $u^0_{0,0}=1$. Then second formula yields $u^1_{0,0} = 1 - 2 * s$. The point (0,0) should not become smaller than its neighbors, and so 
$u^1_{0,0} \geq 0$ which yields $s \geq 1/2$.

Experimentaly we have noticed that if we impose 
\fb& \label{stepfixed} \frac{\Delta t}{\Delta x^2} \leq \frac{1}{2},\fe&
then the algorithm has a good behaviour and remains stable
in the sense that there exists experimentally a (small with
respect to 255) $\epsilon>0$ such that  for any $n \in \entier$
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{The threshold on the spatial gradient norm} : $T_g$ has been fixed experimentally to $6$ for $0$ to $255$ images.


\PAMSS

\POSAMSS

\PMED



\subsection{Case of the AMSS model.} 
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\nabla 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\nabla u\vert^3 curv(u) =
\vert \nabla 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 \nabla 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 \nabla
u\vert^2$ the coefficients $\eta_i=\vert \nabla
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 \nabla 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 \nabla 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 and we conjecture the convergence.
 Contrast invariance are obtained asymptotically by
taking a little time step $\Delta t$. The 
 experimental results presented in figures \ref{MEANIMAGE} and 
\ref{AFINIMAGE}
have been obtained by using these schemes with 
 $\Delta x =1$ and $\Delta t = 0.1$ in the case of mean curvature 
motion and $\Delta t = 0.01$ in 
 the case of affine curvature motion. One has to take $\Delta
t$ that small because unless experimental stability is
achieved with $\Delta t \leq 0.1$, the experimental affine
invariance experimentally needs $\Delta t < 0.05$ (see \ref{figure_amss_affine_char}).

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


\paragraph{The case of the curvature motion.}
Setting the inter-distance between pixel $\Delta x$ to 1, the scale achieved
with $N$ iterations is simply $N / \Delta t$. Now, the scaling is arbitrary.

A good way to normalize the scale is to define the scale by the radius of
a disapearing circle. The boundary of such a circle move at a speed equal to the curvature that is the inverse of the radius. We have for a disk of radius $R(t)$
$$\frac{dR(t)}{dt} = - \frac{1}{R(t)}$$
which yields
$$ \frac{1}{2} (R^2(0) - R^2(t)) = t$$
the disk disappears in scale, when $R(t)= 0$ that is at scale $t=R^2(0)/2.$
This last relation gives the scale normalization:
In order to achieve the normalized scale $T$  (at which any disks (or any shapes included in) of radius less or equal to $T$ has disapeared), we have the equation scale $t = T^2/2$, and the number of iterations needed equal to 
$$ 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 disappearing time is therefore $t= \frac 3 4 R^\frac{4}{3}$.
As for the curvature motion, we define the normalized scale $T$ at which any disks of radius less or equal to $T$ has
disappeared. In order to achieve scale $T$, the number of iteration needed is
$$ N = \frac 3 {4 \Delta t} T^\frac{4}{3} $$


\subsection{The Evans Spruck extension and contrast invariance.}

Both schemes (M.C.M and A.M.S.S) presented above are only asymptotically contrast invariant. But, numerically they are not. 
Indeed, we have seen that a contrast operator can not create new gray level. Now, starting with a binary image $u^0$ and applying a scheme defined by such formula
$$ u^{n+1} = u^{n} + \Delta t (...)$$
we can not be sure that $u^{n+1}$ is also a binary image.

A natural idea to overcome this problem is the following. Starting with a binary image: apply the scheme until the expected scale is achieved, then
binarize the obtained image (just thresholding). This of course works only for binary images, however the Evans Spruck extension (see section \ref{TheEvansSpruckSection}) gives us the key to extend this to general images.

The contrast invariance can be fully obtained, by applying the process of each level set separatly. The procedure is then the following~:


We start with a image $u_0$ and we construct its version at scale $t$ : $u(t,\fxx)$.
\begin{description}
\item{Initialization~:} set $u(t,\fxx) = 0, \forall \fxx, t$

\item{For each $\lambda$, in increasing order}
\begin{itemize}
\item Set $v(\fxx)$ the caracteristic function of the level set $\lambda$ of $u_0$. (that is equal to $1$ inside the level set, $0$ outside.)
\item Apply on $v$ the MCM or AMSS FDS-scheme until scale $t$. This yields the images $w(t,.)$.
\item Set $u(t,\fxx)=\lambda$ at each point $(t,\fxx)$ where $ w(t,\fxx) \geq 0.5$
\end{itemize}
\end{description}
with respect to the simple FDS scheme, this clearly multiplies the number of needed computation, since the AMSS or MCM have to be solved for each level ! However, The resulting images are more closed to the theoritical solutions of these equations.

Note that explained in chapter \ref{CIOMM}, the extension is not unique. The preceding scheme corresponds to the u.s.c. extension. One could also use
the l.s.c extension for which the last step has to be changed into
\begin{itemize}
\item Set $u(t,\fxx)=\lambda$ at each point $(t,\fxx)$ where $ w(t,\fxx) > 0.5$.
\end{itemize}

Of course, these two extensions differ. This can be easily seen in case of image that differs when considered as u.s.c. and l.s.c. An example, is the chessboard image (see figure{\ref{damieretextensions}}). The image data does not say anything about the value of the function at the borders and corners. If the image is considered as u.s.c. then the borders values is white, conversely it would be black when considered as l.s.c. In pratice, the FDS does not care about the borders and corners values, since it is to rough. 

However, making the Evans Spruck extension of the FDS scheme in fact implies an implicit choice for these values. 

\PDAMIERAMSS

\subsection{Problems at extrema. \label{FDSgradientnull}}
 For the M.C.M or the A.M.S.S. occurs the problem to define numerically the
equation when $|\nabla u|=0$.
	For the first the right part of the equation is simply not defined, for the second by continuity, one could set $\vert \nabla u^n \vert ^2u^n_{\xi\xi})_{i,j}^{\frac{1}{3}}=0$. Now, numerically, this would implies that isolated black or white pixels will not change !

As a consequence, the MCM and AMSS finite difference schemes have a slightly bad behaviour at extrema. In the figure \ref{figure_amss_extrema}, we added to the image a strong salt and pepper noise. On a such image more than a quarter of the pixels are local extrema. The extra-diffusion produced on these points by the laplacian (as suggested in the preceding sections) gives strange structures. We show some ways to avoid this spurious behavior~: 

\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 extrema that should anyway theoritically disappear.

\item One can use the Evans-Spruck extensions. This multiplies then the number of computations by the initial number of gray-level. (e.g. 256).
\end{itemize}


\PMAXAMSS

\PAMSSAFFINE

\subsection{Conclusion on finite difference schemes.}

We have seen that standard  finite difference schemes can not 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 It is not fully contrast invariant. We see appearance of a slight blur around edges. And, 
spurious diffusion around the extrema. Extrema problem can be handle by treating them in a different
way : that is flattening them immediately by applying the extrema killer first.

In order to be fully contrast invariant, the only way we know is then
by applying the equation to the characteristic functions of all the level sets. 
This multiplies the total computation by the number of different gray level in the original 
image (often 256 !).

Now, as shown in figure \ref{figure_amss_extrema}, standard scheme with combination with the extrema killer
might give results ``good enough'' for some applications.

\item
The worse drawback is in fact the affine invariance, (or the Euclidean invariance for the curvature motion). 
Since the scheme works on a grid, motion of the level curves is quantized to square steps.
As consequence : we can not garanty the affine invariance and a motion proportional to the curvature power one third.

The only way to cope this problem is to go out of the grid, which is the aim of the section \ref{sectionlionel}.

\end{enumerate}


%%\section{Basic Matheron Operators.}

\section{Curve evolution.}

Given a list of points $\fxx(n)$ approximating in a polygonal way a curve, we would like
to define the curvature evolution of it. As we have seen, the curvature evolution is equivalent to an intrinsic heat equation on the curve itself.

As seen, in chapter \ref{chapter1}, the convolution by a Gaussian can be well approximated by the iteration of a symmetric and positive kernel. Denoting by $s$ the arc length. We choose the simple following kernel: $k(0) = 0.5$, $k(-D)=0.25$, $k(D) = 0.25$ and $k(s)=0$ everywhere else. $D$ is a distance parameter controlling the precision of the scheme.
The smaller it is the better, but higher number of iterations is then needed.

We set:
$$\fxx(t+1,.) = k \times \fxx(t,.)$$
That is, the point number $n$ is changed by
$$\fxx(t+1,n) = \fxx(t,n) + 0.25 * D *  (\frac{(\fxx(t,n+1) - \fxx(t,n) )}{|\fxx(t,n+1) - \fxx(t,n)|}  + \frac{( \fxx(t,n-1) -\fxx(t,n) )}{|\fxx(t,n-1) - \fxx(t,n)|})  $$
Numerical stability problems might happen when two consecutive points have a small distance between them (typically smaller than $D$). One way is to keep one point every two that are too close to each other.

This scheme has been used to generate Figure \ref{MACKW2L}.

To be noted that this simple scheme is very rough. For better precision,  one can
still implement the classical convolution by the gaussian kernel and iterate it with re-parametrization of the curve.
This scheme can be extended to the affine shortening by replacing the arc length by the affine arc length. 
Obtaining that way an monotonous and affine invariant scheme seems to be quite difficult.

\section{Affine Plane Curve Evolution Scheme. \label{sectionlionel}}

The idea here is to go with the curve evolution instead of going with the image evolution.  
In a sense, we go in the converse way than Osher and Sethian idea.
Indeed, they proposed in order to simulate the evolution of a curve by the curvature
motion to make evolving its distance function by the image version of the equation.
Here we do the other way around~: we extract all level lines of the image. We
then compute their evolutions by the affine shortening. And, at last reconstruct
an image out of them. This image is then the AMSS evolution of the original.

To be noted that since the curves are extracted from an digital image, they
are initialy embeded into the image grid. We can therefore represent each curve
by a polygone.
Secondly, the affine shortening will be executed by alternating an affine erosion
and an affine dilation.


\subsubsection{A fast algorithm}

In general, the affine erosion of $X$ is not simple to compute, because
it can be strongly non local. However, if $X$ is convex, then it has
been shown in \cite{Moisan:1998:APC} that it can be exactly computed in linear time.
In practice, $c$ will be a polygon and the exact affine erosion of $X$
---whose boundary is made of straight segments and pieces of hyperbolae---
is not really needed~; numerically, a good approximation  
by a new polygon is enough. Now the point is that we can approximate 
the combination of an affine erosion plus an affine dilation 
of $X$ by computing the affine erosion of each {\em convex component}
of $c$, provided that the erosion/dilation area is small enough.

\bigskip

The algorithm consists
in the iteration of a four-steps process~:

\begin{enumerate}
\item {\bf Break the curve into convex components.}
This operation permits to apply the affine erosion to convex pieces of
curves, which is much faster (the complexity is linear) and can be
done simply in a discrete way. The main point is to take into account
the finite precision of the computer in order to avoid spurious (small and
almost straight) convex components only due to numerical artifacts. 
\item {\bf Sample each component.}
At this stage, points are removed or added in order to guarantee an optimal
representation of the curve that is preserved by step 3.
\item {\bf Apply discrete afffine erosion to each component.}
\item {\bf Concatenate the pieces of curves obtained at step 3.}
This way, we obtain a new closed
curve on which the whole process can be applied again.
\end{enumerate}

\bigskip


The curve has to be broken at points where the sign of the determinant 
$$d_i = \D{P_{i-1}P_i,P_iP_{i+1}}$$
changes. Numerically, we use the formula
\begin{equation}
d_i = (x_i-x_{i-1})(y_{i+1}-y_i) - (y_i-y_{i-1})(x_{i+1}-x_i)
\label{di.num}
\end{equation}
Since we are interested in the sign of $d_i$, we must be careful
because the finite numerical precision of the computer can make
this sign wrong. Let us introduce the relative precision 
of the computer
\begin{equation}
\varepsilon_0 = \max\{x>0, \; (1.0\oplus x)\ominus 1.0 = 0.0 \}.
\label{defprec.equ}
\end{equation}
In this definition, $\oplus$ (resp. $\ominus$) represent the computer addition
(resp.  substraction), which is not associative.
When computing $d_i$ using (\ref{di.num}), the computer gives a result
$\tilde{d}_i$ such that $|d_i-\tilde{d}_i| \leq e_i$,
with \begin{eqnarray*}
e_i = \varepsilon_0 \left( \bigh \right. & \; & 
|x_i-x_{i-1}|(|y_{i+1}|+|y_i|) 
+ (|x_i|+|x_{i-1}|)|y_{i+1}-y_i| \\
&+& \left. |y_i-y_{i-1}|(|x_{i+1}|+|x_i|)
+ (|y_i|+|y_{i-1}|)|x_{i+1}-x_i| \;\;\;\bigh\right).
\end{eqnarray*}
In practice, we take $\varepsilon_0$ a little bit larger than its 
theoretical value to overcome other possible errors (in particular, errors
in the computation of $e_i$). For four-bytes C {\it float} numbers, we use
$\varepsilon_0 = 10^{-7}$, whereas the theoretical value (that can be
checked experimentally using (\ref{defprec.equ})) is 
$\varepsilon_0 = 2^{-24} \simeq 5.96 \;10^{-8}$.
For eight-bytes C {\it double} numbers, the correct value would be
$\varepsilon_0 = 2^{-53} \simeq 1.11 \;10^{-16}$

\bigskip

The algorithm that breaks the polygonal curve into convex components 
consists in the iteration of the following decision rule~:

\begin{enumerate}
\item If $|\tilde{d}_i| \leq e_i$, 
then remove $P_i$ (which means that to new polygon to be considered from
this point is $P_0P_1...P_{i-1}P_{i+1}...P_{n-1}$)
\item If $|\tilde{d}_{i+1}| \leq e_{i+1}$, 
then remove $P_{i+1}$ 
\item If $\tilde{d}_i$ and $\tilde{d}_{i+1}$ have opposite signs, 
then the middle of $P_i, P_{i+1}$
is an inflexion point where the curve must be broken
\item If $\tilde{d}_i$ and $\tilde{d}_{i+1}$ have the same sign, 
then increment $i$
\end{enumerate}  

This operation is performed until the whole curve has been visited. The
result is a chained (looping) list of convex pieces of curves.


\bigskip

\noindent\noindent
$\bullet$ {\bf Sampling}

At this stage, we add or remove points from each
polygonal curve in order to ensure that the Euclidean distance between 
two successive points lies between $\varepsilon$ and $2\varepsilon$
($\varepsilon$ being the absolute
space precision parameter of the algorithm).


\bigskip

\noindent\noindent
$\bullet$ {\bf Discrete affine erosion}

This is the main step of the algorithm~: compute quickly
an approximation of the affine erosion of scale $\sigma$ 
of the whole curve. 

\bigskip

The first step consists in the calculus of the ``area'' $A_j$ 
of each convex component
${\cal C}^j = P^j_0P^j_1...P^j_{n-1}$, given by
$$A_j = \frac 1 2 \sum_{i=1}^{n-2}\D{P^j_0P^j_i,P^j_0P^j_{i+1}}.$$
Then, the effective area used to compute the affine erosion is
$$\sigma_e = \max\left\{\frac{\sigma}{8},\min_j A_j\right\}.$$
We restrict the erosion area to $\sigma_e$ (which is
less than $\sigma$ in general) because the simplified algorithm
for affine erosion (based on the breaking of the initial curve
into convex components) may give a bad estimation of the continuous 
affine erosion+dilation when the area of one component is less than
the erosion parameter. The term $\sigma/8$ is rather
arbitrary and guarantees an upper bound to the number of iterations
required to achieve the final scale.

\bigskip

Once $\sigma_e$ is computed, 
the discrete erosion of each component is defined as the succession
of each middle point of each segment $[AB]$ such that
\begin{enumerate}
\item $A$ and $B$ lie on the polygonal curve
\item $A$ or $B$ is a vertex of the polygonal curve
\item the area enclosed by $[AB]$ and the polygonal curve is equal 
to $\sigma_e$
\end{enumerate}
These points are easily computed by keeping in memory and updating 
the points $A$ and $B$ of the curve plus the associated chord area.

\bigskip

Notice that if the convex component is not closed (which is the case
if the initial curve is not convex), its endpoints are kept.

\PLIONEL1

\bigskip

\noindent\noindent
$\bullet$ {\bf Iteration of the process}



To iterate the process, we use the fact that if $E_{\sigma}$ denotes the
affine erosion plus dilation
operator of area $\sigma$, and $h=(h_i)$ is a subdivision 
of the interval $[0,H]$ with $H= T/\omega$ and 
$\omega = \frac{1}{2}\left(\frac{3}{2}\right)^{2/3}$, then as we are going
to show further,
$$E_{(h_1-h_0)^{3/2}}\circ E_{(h_2-h_1)^{3/2}} \circ ... 
\circ E_{(h_n-h_{n-1})^{3/2}} \;\;\left(\bigh c_0 \right) \;\;
\longrightarrow \;c_T$$
as $|h| = \max_i h_{i+1}-h_i \rightarrow 0$, where $c_T$ is
the affine shortening of $c_0$ described above by (\ref{ass.equ}).

\bigskip

\noindent\noindent
$\bullet$ {\bf Comments}

The algorithm takes a curve (closed or not) as input, and
produces an output curve 
representing the affine shortening of the input curve
(it can be empty if the curve has disappeared) . The parameters are
\begin{itemize}
\item $T$, the scale to which the input curve must be smoothed
\item $\varepsilon_r$, the relative spacial 
precision at which the curve must be 
numerically represented (between $10^{-5}$ and $10^{-2}$ when using four
bytes C {\it float} numbers).
\item $n$, the minimum number of iterations required to compute the
affine shortening (it seems that $n \simeq  5$ is a good choice). 
From $n$, the erosion area $\sigma$ used in step 3 is computed
with the formula 
$$\sigma^{2/3} = \frac{\alpha \cdot T^{4/3}}{n}.$$
Notice that thanks to the $\sigma/8$ lower bound for $\sigma_e$, the effective
number of iterations cannot exceed $4n$.
\item $R$, the radius of a disk containing the input curve, used to obtain
homogeneous results when processing simultaneously several curves. 
The absolute precision $\varepsilon$ used at step 2 is defined by
$\varepsilon = R \varepsilon_r$.
\end{itemize}


The algorithm has linear complexity in time and memory, 
and its stability is ensured by the
fact that each new curve is obtained as the set of the middle points
of some particular chords of the initial curve, defined themselves
by an integration process (an area computation). Hence, no derivation or
curvature computation appears in the algorithm.

\newpage
\zbibnote

\noindent
{\bf Planar curve motion by morphological schemes}
\newline
The fully invariant affine curve evolution geometric algorithm which we 
presented  was found by 
Moisan \cite{Moisan:1998:APC}. Its implementation for {\it all} level lines
of an image was realized in Koepfler \cite{Koepfler:1999:GMR}.
Cao and Moisan \cite{CaoMoisan:2000:CDE} have generalized 
this curve evolution approach to curvature 
motions  at arbitrary speed of the curvature. They succeeded in numerically 
moving curves at velocities proportional to the power 10 of curvature.
Lisani et al. \cite{Lisani:2000:AIM} use the geometric scheme for fast recognition
and image comparison algorithms.


\medskip
\noindent
{\bf Representation of an image by its level lines}
\newline
This is a very recent way of representing an image probably because it is
or was computationaly expensive.  Caselles et al. 
\cite{Caselles:1996:KP, Caselles:1999:TML} have introduced this representation
as very adequate for shape and occlusion analysis. Monasse and Guichard 
propose a way to compute in a computationaly fast and efficient way the set
of level lines organized as a tree. We did not present here their 
Fast Level Set Transform  \cite{Monasse:2000:FCC, MnsseGchrd3} but have 
described their implementation of 
the extrema killer and their comments on 4- and 8-connectivity.
  

\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, Kacur:1997:SAD}.
Weickert and the \"{U}trecht school \cite{Niessen:1994:NAG, Weickert:1996:NDS, 
Weickert:1997:PIA, Weickert:1998:ERSa} adress many aspects of 
implementation of nonlinear scale spaces, namely speed, paralelism 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"). 

\medskip
\noindent 
{\bf Morphological schemes for the mean curvature motion }
\newline

The first selfconscious implementation of the curvature motion as a sup-inf
operator was made by Catt\'{e} et al. \cite{Catte:1994:MSM}. Now, Merriman,
Bence and Osher \cite{Merriman:1992:DGM} discovered 
that the gaussian weighted iterated median filter 
implements the mean curvature motion. Since the median filter
 is  usually implemented by superposition principle as a sup-inf operator,
we might say that most implementations of the iterated median filter are
numerical schemes for the mean curvature motion.  Those schemes are 
more than inaccurate, though, since only integer speed are allowed. 
Due to this inaccuracy, and as shown in 
Gan and Mao  \cite{Gan:1991:TCT}, discrete median filters are nilpotent
and attain fixed points. 



