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


\begin{figure}[!ht]
\centerline{\hbox{ \epsfxsize=12cm \epsfbox{DDIR/dynsh.ps} }}
 \caption[Gauss convolution+threshold=median]{
\label{dynamicshapeexperiment} ``Gauss convolution+threshold=median''.  A convolution with a gaussian followed by a
threshold yields a weighted median filter. The alternate iterations on a binary shape of the linear smoothing and of a
$1/2$ threshold emulate the mean curvature motion. Left : original picture, then, result of the convolution with a
small gaussian kernel, in continuation threshold at the mid-level of the blurred image, right : result after few
iterations of this alternate process. }
\end{figure}

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}[!t]\label{fatteningeffect}
\setkeys{Gin}{width=12cm}
\includegraphics{FDIR/extrema.eps}
 \caption[Extrema killer implementation]{\label{extremakillerimplementation}    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{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.


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