This chapter has been written in collaboration with T. Cohignac, F. Eve, F. Guichard, and J.M. Morel and has been submitted to the International Journal of Computer Vision.

\skipaline

\skipaline
 

{\small {\noindent \bf Abstract.}
 The fundamental equation of image processing
 $$\frac{\partial u}{\partial t} = \vert Du \vert (curv(u))^{\frac{1}{3}} $$
has been discovered recently with existence proofs in the sense of viscosity 
solutions [AGLM-1-2-3]. Another approach, by modelling curve evolution has been introduced
in [SaTa] and yields the intrinsic equation $\frac{\partial C}{\partial t}=C_{ss}$
where $s$ is the "affine length"  along the curve.
The second equation describes the motion of the level curves of a solution $u$ 
of the first equation; therefore both models are equivalent. We introduce a third 
reformulation for this equation, coming from the "mathematical morphology" set
analysis. This diversity of formulations leads to a rich discussion about 
numerical schemes: classical finite differences, curve evolution schemes, and
set evolution schemes. We compare the numerical results and discuss the
new kind of problems arising in their numerical analysis: Indeed, the fundamental
equation has been created for achieving multiscale shape analysis. Therefore, geometric, affine and
scale invariance properties should be very accurately satisfied. We
introduce two new curve evolution schemes for the equation which
achieve all invariance properties. We also state (and justify) several convergence
conjectures for the new introduced schemes and discuss the experimental results
obtained on several shapes. They lead to a crossvalidation of two very different schemes. 
}

\newpage


\section [ Introduction.] { Introduction.}


In the following, we consider an image $u_0(x)$ defined on the plane and with
real values representing the brightness.  We also consider smooth Jordan curves
 $C(s)$ where s denotes either the euclidean parameter, defined up to a shift
by the relation $\vert C_s\vert = 1$, or the affine parameter, defined by $\vert det(C_s, C_{ss})\vert = 1$ ([SaTa]).
We shall discuss some numerical schemes for five equations in $\real^2$ related with image processing.
 
 The most important equation to be considered is the fundamental equation of image processing (the {\it Affine Morphological Scale Space (AMSS)}),
\fb& 
\frac{\partial u}{\partial t} = \vert Du \vert (curv(u))^{\frac{1}{3}} ,\qquad u(0)=u_0
\label{m2cm} \fe&



It can be interpreted as the motion of the level curves of u by the intrinsic 
equation 
\fb& 
  \frac {\partial C}{\partial t} = C_{ss} \fe& 
where $s$ is the affine parameter. Indeed, this last equation can be rewritten [SaTa2] as

\fb& 
 \frac {\partial C}{\partial t} = curv(C(s))^{1/3} \vec n(s) \fe&

where $\vec n(s)$ is the normal at $s$ to the curve $C(s)$ and $curv(C(s))$ is the curvature of $C$ at $s$.

\skipaline

In order to make the discussion more fruitful, 
we shall also include in the discussion the euclidian analogs of these equations,
that is 
\fb&  
  \frac{\partial u}{\partial t}= \vert Du \vert curv(u) ,\qquad u(0)=u_0 \label{eq:mcm} \fe& 

\fb& 
  \frac {\partial C}{\partial t} = C_{ss} = curv(C(s)) \vec n(s) \fe& 

 where s is the euclidian parameter. Indeed, at a point $x=C(s)$ where $u$ is smooth and $Du\neq 0$, one has 
$$curv(u).\frac {Du}{\vert Du \vert} = curv(C(s)) \vec n(s) = C_{ss}$$ where $C$ 
denotes the level curve of $u$ passing by $x$.
(What $s$ is, length or affine parameter, will be clear in context). The relation between the equations (4), (5) has been discovered by Barles, Osher and Sethian ([Bar], [OsSe]).

\skipaline

Equations (2) and (5) have been solved, with existence, regularity 
and uniqueness results when the initial Jordan curve is smooth 
in [Grayson, Gage-Hamilton, Sapiro-Tannenbaum].  A general theory has been
made for equations (1) to (5) by using the viscosity solution framework [AGLM-1-2-3].  It
also yields existence and uniqueness results and implies the equivalence of both
approaches ([EvSp], [CGG] ...).  Therefore, (at least) two kinds of schemes
are available for the applications in image processing :{\it we can either identify
a picture with the set of its level curves and then evolve them jointly by using
the curve evolution framework.  We also can evolve the image itself by using the
PDE formulation}.  We shall introduce a third numerical approach, inspired from the 
{\it Mathematical Morphology} [Serra].  This approach consists in considering the level 
sets of the image $X_a = \{ u(x) \ge a\} $ and to handle them by operations which are called
"openings", "closings", "median filters", etc.... We shall use recent works by Evans-Spruck, Merriman-Bence-Osher, Barles-Georgelin, Evans, Lasry-Morel which prove that the iteration of 
these filters tends theoretically towards equation (1) or (4).
Before starting with the discussion, let us mention that many ideas for the schemes which we discuss
here come from the image processing litterature, but they were in many cases 
not correctly implemented or interpreted {\it because the subjacent PDE models were not 
identified}.  We shall make this assertion more precise in the next sections.

\skipaline

In section $2$, we shall start with an explanation of the use of the 
fundamental equation in image processing.  This explanation is necessary, 
because it leads to ask for new kinds of (long time) invariance properties 
when defining numerical schemes. These invariance properties are a key to 
applications, e.g. in shape recognition. Therefore, we shall recall an axiomatic
characterization of equation (2) which was developed in [LoMo], simplifying 
the founding approach of [AGLM1].  This axiomatization leads us to a list of 
invariance requirements for numerical schemes which will be the basis of our
discussion through the paper.

In section $3$, we present and discuss classical finite difference schemes, close to the Osher-Sethian formulation.

In section $4$, we relate the fundamental equation to a new "affine morphological" set 
theory and derive new schemes for these equations.

In section $5$, we discuss curve evolution schemes and show why they can be understood
as special implementations of morphological schemes.

Section $6$ is devoted to the experimental comparison of the numerical schemes and
 to the precise statement and justification of several conjectures about their convergence.

\skipaline

\noindent {\bf Aknowledgements:\rm} \rm The last author thanks Jean-Michel Lasry and
Stanley Osher for valuable conversations. Part of this work was made while the last
author was visiting Cognitech-Inc.

\newpage 

 
\begin{center}

{\bf \Large PLAN}

\end{center}

\skipaline


\skipaline


{\noindent \bf 1. Introduction.

\skipaline


\noindent2. The four principles of Multiscale Shape Analysis and the fundamental equation.

\skipaline


2.1 The fundamental equation of Shape analysis.

\skipaline


2.2 The fundamental equation of Image Processing. 

\skipaline


2.3 Invariance and stability requirements for numerical schemes of the fundamental equation.

\skipaline


\noindent 3. Quasilinear finite difference Schemes.

\skipaline


3.1 Luis Alvarez's scheme.

\skipaline


3.2  Luis Alvarez, Frederic Guichard's scheme.

\skipaline


\noindent 4. Morphpological (set evolution) schemes.

\skipaline

4.1 A theoretical scheme.

\skipaline

4.2 Iterated median filters and curvature motion.

\skipaline



\noindent 5. Purely morphological (curve evolution) schemes.

\skipaline

5.1 Introduction.

\skipaline

5.2 Quasilinear intrinsic diffusion schemes.

\skipaline

5.3 A geometrical curve evolution scheme.

\skipaline

\noindent 6. Experimental results.

\skipaline

\noindent 7. Conclusion and conjectures.


\skipaline


\noindent Appendix
  


}

\newpage





\section [The four principles of Multiscale Shape Analysis and the fundamental equation.]{The four principles of Multiscale Shape Analysis and the fundamental equation.}

Let us start with a list of four principles of shape analysis which lead to the
fundamental equation (For more details, see [AGLM], [LoMo]). These principles will
be essential for establishing requirements for numerical schemes and in order to
justify the conjectures.

\skipaline

Define a {\it shape or silhouette} as a closed set X whose boundary is a Jordan curve of $\real ^2$. 
 We denote by $T_t(X)$  the {\it shape analysed at scale} $t$. $X$ is identified with its
characteristic function $X(x)=1$ if $x \in X$ and $0$ else.  
We call multiscale analysis any family operators $(T_t)_{t \ge 0}$ acting on shapes
and we set $X(t)=T_t(X)$. 
  
\skipaline

\rm We first consider a principle which is called in vision theory
  the {\it pyramidal architecture} property.
  According to this assumption, $T_t$ can be computed from $T_s$  for any
  $s\le t$ and $T_0$ is of course the identity. 
 This is natural, since a coarser analysis of the original picture is 
likely to be deduced from a finer one without any dependence upon the
 original picture.  Of course the finest picture analysis is the identity.
   In order to formalize the pyramidal architecture we include $T_t = T_{t,0}$ 
into a family of {\it transition operators} $T_{t,s}$ 
 indexed by $0 \le s \le t$  and satisfying 

\skipaline

\noindent {\bf [Pyramidal architecture] \rm}	    	$T_{t+s} = T_{t+s,s}T_s$ for all $0 \le s$, $ 0 \le t$


\skipaline



  The shape analysis process is thought of as a flow of the shape 
through different filters, each associated with a scale $t$. If we discretize
 such a process, it yields a serie of operators
 $S_n = T_{(n+1)h,nh} $ applied successively.   


\skipaline



\rm Assume that $X$ and $Y$ are two silhouettes and that for some
 $x \in \partial Y$ and some $r>0$,
 one has $X \cap B(x,2r) \subset Y\cap B(x,2r)$.
  Assume further that the inclusion is strict in the sense that
 $\partial X$ and $\partial Y$ only meet eventually at $x$.
  Then we shall say that {\it the shape} $X$ {\it is included in shape} $Y$ {\it around} $x$.


\skipaline



\noindent {\bf [Shape local inclusion] \rm} If $X$ is included in $Y$ around $x$,
 then for $h$ small enough, $T_{t+h,t}(X) \cap B(x,r) \subset T_{t+h,t}(Y)\cap B(x,r)$.

\skipaline



This last axiom implies that the value of $T_{t+h,t}(X)$ for $h$ small,
 at any point $x$, is determined by the behaviour of $X$ near $x$.
  We are allowed to take $r$ infinite.
   Therefore the shape local inclusion principle also implies that if
 a shape is globally contained in another, this order is preserved for every scale [MaMo].
 Notice that the combination of both pyramidal and shape local inclusion principles
is generally called in Image Processing {\it causality} 
[Koe], [Wit], [Kim], [SaTa], [KTZ].
Indeed, the causality assumption in Image Processing states that when the
scale increases, no new structure appears in the image. Hummel [Hum] noticed
that the causality assumption can be reduced to maximum principle.


\skipaline



Both preceding principles allow, as we shall see,
 to localize the shape analysis in space and time and we therefore only need
 to state what the multiscale analysis makes of very simple shapes in order
 to specify it.  So we add a basic "principle" which will state what happens
 to disks.   As we shall see, disks are somehow the "basis" on which every
 shape can be decomposed.


\skipaline





\noindent {\bf [Basic principle] \rm} Let $D=D(x,1/r)$ be a disk with curvature $1/r$ and center $x$.
  Then $T_{t+h,t}(D)$ is a disk with radius $\rho (t, h,1/r ) $ and center $x$.
  Moreover, the function $h \rightarrow \rho(t,h,1/r) $ is differentiable with
 respect to $h$ at $h=0$ and the differential is continuous with respect to $1/r$.

\skipaline

The "basic principle" implies that the multiscale analysis behaves
 in a smooth and isotropic way. 



In the following, we set $$g(t, 1/r) = \frac{\partial \rho}{\partial h}(t, 0, 1/r)$$
\rm Note that $g(t,s)$ is defined for $t \ge 0$ and $s \in \real $, now to define
$g(t,0)$ we must assume that $lim \; g(t,1/r)$ exists when $r$ tends to $+\infty $ or $-\infty $ .
 The radius $r$ may be positive or negative, according to the orientation of the normal $n(x)$.
 In the case where the curve is a Jordan curve enclosing a set $X$,
 we take as a {\it convention that} $n(x)$ is {\it pointing outside} $X$ and {\it the curvature
 is negative if} $X$ {\it is convex at} $x$, {\it positive else}.
It may seem natural to assume therefore that $g(t,\kappa )$ is odd with respect to $\kappa $.
Indeed, this corresponds to the {\bf assumption that a black disk on grey background and
a white disk behave in a same way.} 
 

\skipaline




By the shape local inclusion principle, it is easily seen that the function $g(t, \kappa) = \frac{\partial \rho}{\partial h}(t,0,\kappa)$
 is nondecreasing with respect to $\kappa$.


\skipaline

    


The last shape preserving principle is the affine invariance of shape analysis.
Set $AX(x) = X(Ax)$ for any linear map $A$. The affine invariance can be stated as:


\skipaline



\noindent {\bf [Affine inv]\rm}. \rm	 There exists a $C^1$ function 
$t'(t,A) \ge 0$ defined for any A and $t \ge 0$, such that
 $AT_{t'(t,A),t'(s,A)}=T_{t,s}A$.

It was proved in [AGLM] that we can always take $t'(t,A)=t.\vert det A\vert ^{1/2}$ {\bf up to a rescaling} (See Appendix).

\skipaline

\rm This relation means that the result of the multiscale analysis $T_t$ is
 independent of the size and position in space of the analyzed features. 
This is very important in the world where we live,
 since the same plane object $X$ can be seen at very different distances and
 orientations and therefore with any affine change $X(x) \rightarrow X(Ax)$.
 A linear map corresponds to the anamorphosis of a plane image when it is presented to
the eye at an arbitrary large distance and with an arbitrary orientation in space. 

\skipaline

\subsection [The fundamental equation of Shape Analysis.]{ The fundamental equation of Shape Analysis.}

 \rm  When a point $x$ belongs to an evolving curve, we denote by
 $\dot {x}$ the time derivative of $x$,
 which is a vector of $\real ^2$.  By $curv(x)$ we denote the curvature
 of a curve which is $C^2$ at $x$. Recall that the curvature is defined
 as the inverse of the radius of the osculatory circle to the curve at $x$.
  The curvature is zero if the radius is infinite.
  By $\vec n(x)$ we denote the vector normal to the curve and pointing 
towards the concavity of the curve (These definitions for curves are compatible with the ones given for sets.)
  
\skipaline

{\noindent \bf Theorem 1.} [LoMo]{\it (i) Under the  three principles pyramidal, local shape inclusion and
 "basic", the multiscale analysis of shapes is governed by a 
 curvature motion equation
 \fb& \label{(Eq1)} \; \dot{x} = g(t, curv(x))n(x) \fe&
 where $g$ is a nondecreasing function.

\noindent (ii) If we assume [Affine inv.], 
then the equation of the multiscale analysis is, {\bf up to a rescaling} :
\fb& \label{(Eq2)} \; \dot{x} = \gamma (t.curv(x))n(x) \fe&
where $\gamma$ is defined as following:
$\gamma (x)=a.x^{1/3}$ if $x\ge 0$ and $\gamma (x)=b.x^{1/3}$ if $x\le 0$     
and $a$, $b$ are two nonnegative values.

\noindent (iii) If we add that $T_t(X^c)=T_t(X)^c$ [Reverse contrast invariance] 
then the function $g$ in (i) is odd and $a=b$ in (ii).
 }
  

\skipaline



\rm By the expression "governed by", we mean that the equation must be verified 
at any point $(t,x)$ of $X(t)$ where the silhouette is smooth enough to give a classical
 sense to both terms of the equation.


\skipaline



\noindent {\bf Remark: \rm}
It has been proved [GaHa], [Gra], [Ang] that the equation (\ref{(Eq1)})
 with $g(t,s)=s$ has smooth solutions.
This equation has been introduced in image processing by [KTZ],
 [MaMo] and [AGLM1,2,3] and [Yui] in different contexts.
 An early version of an algorithm leading to equation (\ref{(Eq1)}) is in [KoDo]. 



The equation (\ref{(Eq2)}) has been first introduced and axiomatically
justified in [AGLM1,2,3].
It has been independently proposed in [SaTa1].  
Existence and regularity of the solution are proved in [SaTa2].



\skipaline




\subsection [ The fundamental equation of Image Processing.] { The fundamental equation of Image Processing.} 

Since the multiscale analysis verifies the obvious inclusion principle that
 if $ A \subset B$ then $T_{t+h,t}(A) \subset T_{t+h,t}(B)$,
 then we can, as well known in the "mathematical morphology" school [Ser], [Mara],
 associate with a picture $u$ the set of its level sets
 $$X_au =\{x, u(x) \ge a\}$$  
Then, assuming that each level set is a union of silhouettes, 
we can simply define $T_t(u)$ as

\skipaline



\noindent {\bf [morphological principle]} For any a, t, h and u: $X_aT_{t+h,t}(u) = T_{t+h,t}(X_au)$  

\skipaline

It is easily seen [LoMo] that if we define a multiscale analysis of images
 from the preceding analysis of shapes by using [morphological principle], it is governed by 
\fb& \label{(Eq3)} \; \frac{\partial u}{\partial t} = g(t, curv(u))\vert Du \vert \fe&



Conversely, this equation has a unique solution in the sense of
 "viscosity solutions" [AGLM]. {\it So a consequence of Theorem 1 above is that
there is a single image multiscale analysis, satisfying [pyramidal], [local
shape inclusion], [basic principle], [affine invariance] and the [morphological principle]}.
Therefore, we call the associated equation the fundamental equation of image processing
$$\frac{\partial u}{\partial t} = \vert Du \vert (t.curv(u))^{\frac{1}{3}} ,\qquad u(0)=u_0$$



\skipaline

{\noindent \bf Remark :} It is easy to obtain the equations (1) to (5) by rescaling the last equations ($t\to \phi (t)$).


\subsection [Invariance and stability requirements for numerical
 schemes of the fundamental equation.] {Invariance and stability requirements for numerical
 schemes of the fundamental equation.}


In order to compare schemes, we shall fix in this section a list of formal requirements which 
derive from the above developed axiomatic analysis.
We have first the properties which imply convergence of the schemes toward
the equations, and second the invariance properties of the equations themselves.
All schemes will be defined as the iteration of a nonlinear discrete filter, which we call $T$.
We set $u^n=T^nu_0$, where $u_0$ is the (discretized) initial image. All schemes
will depend on a scale step $\Delta t$ and a space step $\Delta x$.

\skipaline

We shall consider three different discrete representations of $u^n$ :

\begin{enumerate}
\item The classical discrete representation on a grid $u^n_{i,j}=u^n(i,j)$, with $1\le i\le N$ $1\le j\le N$ .

\item "The morphological representation" of u by a bunch of level sets $X_a=\{ u(x)\ge a\} $,
where $a$ admits, for technological reasons, $255$ values.

\item The "purely morphological representation" of u by its level curves, a 
sequence $C_k$ of Jordan curves discretized by affine finite elements.
(Each curve is given by a sequence of vertices, with real components).

\end{enumerate}

Let $H\la$, $R$, and $A$ the operators defined by $ H\la u(x) = u(\lambda x) $, 
$ R u(x) = u(R x) $
where $R$ is an isometry of $\real^2$, $ A u(x) = u(A x)$  where $A$ is a linear application of $\real^2$.
Since the schemes to be considered will be very different in structure, we cannot 
state all following properties in a completely formal way. What they mean in each
particular case will be clear in context.

\vspace{2mm}

\begin{itemize}

\item \CONS : We shall say that a scheme is consistent if
 the discrete operator $T$ approximating the second member of (1) or (2) 
tends to the differential one,
 when the steps $\Delta x$ and $\Delta t$ tend
 to 0 in a suitable way.

\item \CONV : $\forall t \; T^n u \rightarrow T_t u \; a.e. $
when $\Delta t $ and $\Delta x$ tend to 0 in a suitable way, and 
$n \Delta t \rightarrow t$.

\item \STAB : Let c,d be real numbers. If c$\leq u^{n}(x) \leq$d
then c$\leq u^{n+1}(x) \leq$d.

\item \MONO : If $u^{n}(x) \ge v^{n}(x)$ for all x in $\real^2$, then $u^{n+1}(x) \ge v^{n+1}(x)$.

\item \Morp : For all nondecreasing functions h : $T h(u) = h(T u)$. This means
that {\it the level sets} are handled independently.

\item \MORP : There is no interaction between the connected components of level sets.
(This is a more restrictive property than \Morp and means that {\it the level curves} are
handled independently).

\item \ROTA : $ T R_\theta = R_\theta T $.

\item \ECHE : $\forall l, n \in \entier \; \exists n'(n,l),\; T^n H_l = H_l T^{n'} $.

\item \AFFI : $\forall n \in \entier \; \exists n'(n,A),\; T^nA = AT^{n'} $
(for any matrix $A$ with entire coefficients (only for equation \ref{m2cm}.))

\item \QUAS : Let T be a discrete operator approximating a differential operator of u on a grid. (For instance, the second derivative in a direction). T is quasilinear if it can be written as a linear combination of the $u^n_{ij}$, where i,j belong to a fixed stencil. The coefficients of this linear combination may vary. 

We add this (not very formal) notion because quasilinear schemes are easy to implement.


\end{itemize}

 In the next sections, we define schemes and discuss which of the above properties is satisfied.

\skipaline

\section [Quasilinear finite difference schemes.]{ Quasilinear finite difference schemes.}

These schemes are based on a diffusive interpretation of the equation (4).
Indeed, if $\xi$ is a unit vector, such that $\nabla u . \xi = 0$,
and $u_{\xi\xi}$ the second derivative of u in the direction $\xi$,
we have :
$$ (4) \equi& \ut& = u_{\xi\xi} $$
This formulation yields an anisotropic diffusion equation. Anisotropic, 
because we diffuse only in one direction $\xi$, depending on the gradient.
(compare with the isotropic diffusion ($\ut& = \Delta u$) i.e. the
heat equation).

Equation (1) can be rewritten as 
$ \ut& = ( \nabla u )^{\frac{2}{3}} (u_{\xi\xi})^{\frac{1}{3}} $, 
where $( \nabla u )^{\frac{2}{3}}$ is understood as the speed of the diffusion.


The schemes which we wish to discuss compute $u_{\xi\xi}$ by a quasilinear formula in the sense defined above.

\skipaline



\subsection [ Luis Alvarez's scheme.]{ Luis Alvarez's scheme.}

The numerical scheme of Alvarez uses the main directions of the grid,
 given by the vectors
 $$(1,0),(2,1),(1,1),(1,2),(-1,0),(-2,1),(-1,1),(-1,2)$$
At every point and every time there will be diffusion in only one of these directions.
So, the algorithm is:

Step 1 : We search the nearest main direction $\eta$ of $\xi$.

Step 2 : We assume that $u_{\xi\xi}$ is equal to $u_{\eta\eta}$. And we discretize $u_{\eta\eta}$ by  
$$u_{\eta\eta}^{n+1}=\frac{u_{i-1,j}^{n+1} - 2 u_{i,j}^{n+1} + u_{i+1,j}^{n+1}}{(\Delta x)^2}, \hbox{ if }
\; \; \eta = \vec{i} $$
$$ u_{\eta\eta}^{n+1}=\frac{u_{i-1,j-1}^{n+1} -2 u_{i,j}^{n+1} + u_{i+1,j+1}^{n+1}}{2(\Delta x)^2}, \hbox{ if }
\; \; \eta = \vec{i}+\vec{j} $$
and so on... Then define an implicit scheme by

$$ u^{n+1} - u^{n} = \Delta t u_{\eta\eta}^{n+1}$$

Step 3 : Solve the implicit scheme and return to step 1.

\skipaline

This scheme satisfies \ROTA for rotations of radius $\frac{\pi}{4}$, and symmetries with respect to main axis.
\STAB is also true. But it does not satisfy \CONS nor \Morp! We remark experimentally that 
the scheme is robust (\STAB). But, the code has not enough 
geometrical invariants. For instance, a disk becomes after some time
 a polygon... This can be checked on the analysis of the spiral (fig. {\it Spiral.img})
So, {\it it cannot be used for sharp shape analysis because it creates corners which
do not exist in the original shape}. We need more rotational invariance and this 
will be achieved by the next scheme.

\skipaline


\subsection [Luis Alvarez, Frederic Guichard's scheme.] { Luis Alvarez, Frederic Guichard's scheme. }

In the last subsection, we tried to use an adaptative diffusion (in order to keep \STAB),
 now we seek for a good 
estimate of the curvature and gradient direction. We want to have more
isometrical invariants.
 Let us begin by computing a discretization of $curv(u)$ as
 rotation invariant as possible. We write curv(u) as a linear combination of the values of u
in the 9 points.

$$ (u_{\xi \xi})_{i,j}= \sum_{-1\leq l\leq 1} \sum_{-1\leq k\leq 1}
 \alpha_{l,k}^{i,j}(\theta_{i,j}) u_{i,j} + o(\Delta x)^2 $$
where $\theta_{i,j}$ is the angle of the vector $\xi$.
At each point (i,j),we would like to define the weights $\alpha_{l,k}^{i,j}$, in order to have :

\skipaline

(i)  the preceding equality.

(ii)  the obvious geometrical invariants of the grid.

(iii) A relation which implies \STAB
 $$ \sum_{-1\leq l\leq 1} \sum_{-1\leq k\leq 1}
 \vert \alpha_{l,k}^{i,j}(\theta_{i,j}) \vert - 2*\alpha_{0,0}^{i,j}(\theta_{i,j}) \leq 0$$

Now, if we have (i) and (ii) , we cannot have (iii) for all angles $\theta$. Thus, we shall not
have \STAB. In order to have so much stability as possible, we change therefore (iii) into
 
 $$ \sum_{-1\leq l\leq 1} \sum_{-1\leq k\leq 1}
 \vert \alpha_{l,k}^{i,j}(\theta_{i,j}) \vert - 2*\alpha_{0,0}^{i,j}(\theta_{i,j})
\qquad \hbox{as small as possible.} $$ 

At last the functions $\alpha_{l,k}$ must be as smooth as possible 
for stability.

Properties (i), and (ii) imply that 
$$ \left\{ 
\begin{array}{l} 
 \alpha_{1,0} = \alpha_{-1,0} = \lambda_1 (\theta) \\
 \alpha_{0,1} = \alpha_{0,-1} = \lambda_2 (\theta) \\
 \alpha_{1,1} = \alpha_{-1,-1} = \lambda_3 (\theta) \\
 \alpha_{-1,1} = \alpha_{1,-1} = \lambda_4 (\theta) \\
 \alpha_{0,0} = -4 \lambda_0 (\theta)  \end{array} \right. $$

$$ \hbox{where } \qquad \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. $$


and

$ \lambda_0(\theta+\frac{\pi}{2})= \lambda_0(\theta) $
$\lambda_0(\theta+\frac{\pi}{4})= \lambda_0(\frac{\pi}{4} - \theta) $
$ \lambda_0(0) = 0.5$  $\lambda_0(\pi/4) = 0.25$.

We choose for $\lambda_0$ the simplest trigonometric polynomial which verifies these 
equalities, and property (iii) as much as possible (see Appendix). So we have 
$ \lambda_0(\theta)= 0.5 - \cos^2(\theta) + \cos^4(\theta) $ and we deduce the value of 
 $\lambda_1(\theta)$, $\lambda_2(\theta)$, ...

\skipaline

We can show that it is equivalent to say that the computation of $u_{\xi\xi}$ is exact
for a set of polynomials of degree as high as possible. We can include all polynomials of degree 3 plus the polynomial : $(x^2+y^2)(x \sin\theta - y \cos\theta ) ^2$. (See Appendix). 

\skipaline

Our algorithm is following:


\begin {itemize}

\item We estimate the gradient direction $\theta$, and $\xi$ the orthogonal direction
using the 9 points of the stencil 3x3.

\item Computation of $curv(u)$ as above.

\item We proceed, in case of equation

\end {itemize}

$$ \begin{array}{l} \hbox{(\ref{eq:mcm}) : } \qquad u^{n+1}(x) = u^{n}(x) + \Delta t u^n_{\xi\xi}(x)
\\ \hbox{(\ref{m2cm}) : } \qquad u^{n+1}(x) =u^{n}(x) +
 \Delta t (Du)^{2/3} (u^n_{\xi\xi}(x))^{1/3} \end{array} $$

In the case of (\ref{eq:mcm}), we must choose $\Delta t \leq 0.5 (\Delta x)^2$
 in order to have some stability.
Whereas for (\ref{m2cm}), the maximum value of $\Delta t$ cannot be imposed 
theoretically. This scheme verifies \CONS, because it is exact at order 2.
All the geometrical properties (\ROTA, \ECHE, ...) are true for
all the polynomials mentionned above. For other functions, we have these properties only
 for the affine transforms which leave the grid invariant.

\underbar{But, the scheme is not morphological.} 
And, it does not satisfy \MONO nor \STAB. However, it is experimentally  
quite stable.
In fact, we do not expect to have \MONO and \CONS in such a quasilinear scheme
with a fixed stencil. If you impose both, as noticed in [ALM] and proved by [CrLi],
 you do not limit the size of the stencil.



\skipaline


\input{cor2.tex}


