\graphicspath{{./figures/}}


\chapter{A snake from A to Z... \label{chapter17}}


\section{An active contour model}
Boundaries of  objects perceived against a different background induce some discontinuities in the gray level,
resulting in high gradient. Let us call {\it contour} a closed Jordan curve located mostly at high gradient points. The
aim of the active contour methods is to find such contours, starting from an initial curve, usually sketched by hand.
The curve moves numerically from its original location until it reaches a position where it maximizes the image
contrast. We do not call as usual the image $u_0$ because our main focus is not on the ``analyzed'' image but on the
``analyzing image'' describing the snake motion. This last image will be called $u(t,\fxx)$ and the analyzed image will
be called $I(\fxx)$. The boundary detection problem can be formulated as an optimization problem. We shall treat it in
2D but the 3D case has exactly the same formalism. Of course, then, ``curve'' has to be replaced by ``surface". Let us
choose a function $g$ from $\real^2$ into $\real$ representing for each point $\fxx$ a penalty for the curve to pass by
the point $\fxx$. Ideally, $g$ has to be chosen small when the magnitude of the image gradient is large. We shall set
for example
 \begin{equation}\label{defdeg}g(\fxx) = \sigma^2/(\sigma^2+|DI(\fxx)|^2)\end{equation}
 on the image domain, where  $\sigma^2$ is the estimated
variance of the noise and texture around the object. For convenience we shall extend $I$ and  $g$ outside this domain
to $\RR^2$, but assume that
 $g(\fxx)$ is zero for $\fxx\geq R$ large enough.
 Thus, given an image $I$, and an initial
curve $C_0=(\fxx_0(s)), s\in [0, L(C_0)]$, we want to find a curve $C=(\fxx(s)), s\in [0, L(C)]$ that minimizes the
 energy \fb& E(C) = \int^L_0 g(\fxx(s))ds \label{energysnake} \fe& around $\fxx_0$, where $s$ is an  arc
length parameter and $L=L(C)$ the length of the curve $C$.

In the following, we will assume $g$ to be twice differentiable with respect to $\fxx$. In order to achieve this in
practice, $I$ is previously  smoothed by the heat equation and $DI$ is therefore replaced in \eqref{defdeg} by   $G\ast
DI=D(G\ast I)$ where $G$ is a gaussian with small variance. In that way $g$ becomes $C^\infty$ (see Proposition
\ref{proposition1.8}.)

Let us refresh some differential notation on curves. We denote by $\fxx(z)$ a parameterization of a curve $C$ on a
fixed interval $[0,1]$.  Recall that $\fxtau(z)=\frac{\fxx'(z)}{|\fxx'(z)|}$ is the tangent unit vector to the curve
and $\fxn(z)=\fxtau(z)^\perp$ the unit normal.  Notice that if $\fxv$ is a vector, then we can decompose it on the
mobile frame $(\fxtau(z), \fxn(z))$ as
\begin{equation}\label{decompositiononmobileframe}
\fxv=(\fxv.\fxtau)\fxtau + (\fxv.\fxn)\fxn.
\end{equation}
Calling $s(z)$ a length parameter on the curve, defined up to a constant by $s'(z)=|\fxx'(z)|$,  one has
$\fxtau(s)=\frac{\fxx'(z)}{|\fxx'(z)|}$ and

\begin{equation}\frac{\partial \fxtau}{\partial s} = \bm{\kappa}(\fxx(z)),\label{Curvaturevectordefinition}\end{equation}
 which is the curvature vector. Thus, differentiating the tangent vector with respect to $z$
yields \begin{equation}\label{differentiatingthetangentvector}\left(\frac{\fxx'(z)}{|\fxx'(z)|}\right)'=\frac{\partial
\fxtau}{\partial z}=\frac{\partial \fxtau}{\partial s}\frac{\partial s}{\partial z}=
\bm{\kappa}(\fxx(z))|\fxx'(z)|.\end{equation}



\begin{proposition}
Let $C(t)=\fxx(t,s)$  be a curve resulting from the gradient descent of the energy (\ref{energysnake}), starting from
$C(0)$. Assume that $C(t)$ is $C^2$. Then $C(t)$ satisfies the following equation \fb& \frac{\partial\fxx}{\partial t}
= -(Dg(\fxx).\fxn) \fxn)
 + g(\fxx) \bm{\kappa}(\fxx)  \label{equationcurvesnake}\fe&
$$ \hbox{ and } \fxx(0,s) = \fxx_0(s)   \hbox{ that is } (C(0) = C_0)$$
\end{proposition}

\begin{proof} We shall first change the parameterization of the curve $C$ so that its length is no longer a parameter of the
energy. We parameterize the curve with $z \in [0,1]$. We have $ds= |\fxx'(z)|dz$, where $'$ denotes the derivative with
respect to $z$. Thus
$$ E(C) = \int_0^1 g(\fxx(z)) |\fxx'(z)|dz $$
Consider any $C^2$  perturbation of the curve $\fxx(z)$,  which we call $dC$ and denote its parameterization by
$d\fxx(z)$. By $|dC|$ we mean the $C^2$ sup-norm of $d\fxx(z)$ on $[0,1]$. By an easy differentiation,
\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
  E(C+dC) -E(C) &=& \int_0^1 Dg(\fxx(z)).d\fxx(z) |\fxx'(z)|dz\\
  &+&  \int_0^1 g(\fxx(z)) \frac{\fxx'(z)}{|\fxx'(z)|}.d\fxx'(z) dz +o(|dC|).
\end{eqnarray*}
Integrating by parts the last integral, we therefore have by \eqref{differentiatingthetangentvector},
\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
   E(C+dC) -E(C)&=& \int_0^1 Dg(\fxx(z)).d\fxx(z) |\fxx'(z)|dz \\
 &-& \!\!\int_0^1 (Dg(\fxx(z)).\fxx'(z)) \frac{\fxx'(z)}{|\fxx'(z)|}.d\fxx(z) dz  \\
   &-&   \int_0^1 g(\fxx(z)) \bm{\kappa}(\fxx(z)).d\fxx(z) |\fxx'(z)|dz +o(|dC|).
\end{eqnarray*}
Using \eqref{decompositiononmobileframe} with $\fxv=Dg(\fxx(z))$, the two first integrals can be merged and we obtain
 $$
E(C+dC) -E(C)=$$
$$ \int_0^1 (Dg(\fxx(z)).\fxn)\fxn .d\fxx(z)|\fxx'(z)|dz - \int_0^1 g(\fxx(z)) \bm{\kappa}(\fxx(z)).d\fxx(z)|\fxx'(z)| dz +o(|dC|).$$
Let us denote the intrinsic scalar product between two vectorial functions $f$ and $h$ defined on the curve $\fxx(z)$
by
$$\langle f.g\rangle = \int_0^1 f(\fxx(z)) g(\fxx(z)) |\fxx'(z)| dz$$
Thus, $$E(C+dC) -E(C) = \langle dC  ,  \left(  (Dg(\fxx(z)).\fxn) \fxn
 - g(\fxx(z)) \bm{\kappa}\fxx(z))  \right) \rangle +o(|dC|)
$$
We therefore have
$$ \nabla E(C) = (Dg(\fxx(z)).\fxn) \fxn
- g(\fxx(z)) \bm{\kappa}(\fxx(z)). $$ As a consequence the gradient descent for a curve $C=(\fxx(z))$ following the
steepest gradient descent  can be described by the equation $\frac{\partial C}{\partial t}=-\nabla E(C),$ that is
$$
 \frac{\partial\fxx(t,z)}{\partial t} = -(Dg(\fxx(t,z)).\fxn) \fxn
 + g(\fxx(z)) \bm{\kappa}(\fxx(z)).
$$
From this equation we can deduce the  normal motion
$$
 \frac{\partial\fxx(t,z)}{\partial t} = -(Dg(\fxx(t,z)).\fxn) \fxn
 + g(\fxx(z)) (\bm{\kappa}(\fxx(z)).\fxn)\fxn.
$$
This last evolution is obtained from the former one by projecting $\frac{\partial\fxx(t,z)}{\partial t}$ on the normal
line $\RR\fxn(t, z)$. Indeed, in order to describe the geometric evolution of a curve we only need to give the motion
of each one of its points in the  direction normal to its tangent. A last simplification is now obtained when for each
$t$, we choose to reparametrize the curve by an arc-length parameter $s$. In such a case
$\bm{\kappa}(\fxx(t,s))=\frac{\partial^2\fxx}{\partial s^2}(t,s)$ is normal to the curve. Thus we can simply write
$$
 \frac{\partial\fxx(t,s)}{\partial t} = -(Dg(\fxx(t,s)).\fxn(t,s)) \fxn(t,s)
 + g(\fxx(t,s)) \bm{\kappa}(\fxx(t,s))=$$ $$-(Dg(\fxx(t,s)).\fxn(t,s)) \fxn(t,s)
 + g(\fxx(t,s))\frac{\partial^2\fxx}{\partial s^2}(t,s),
$$
or, if we omit the variables:
$$
 \frac{\partial\fxx}{\partial t} = -(Dg(\fxx).\fxn) \fxn
 + g(\fxx) \bm{\kappa}(\fxx)=-(Dg(\fxx).\fxn) \fxn
 + g(\fxx)\frac{\partial^2\fxx}{\partial s^2}
$$
\endproof

Unfortunately, we cannot be sure that such an evolution yields a regular curve for all $t$. In fact,
 it is in general
false, since topological changes for the curve can occur, which imply the appearance of infinite curvatures (see Figure
\ref{snaketwoballs}.)


By a straightforward adaptation of the proof of Proposition \ref{corollary5.9}  one immediately obtains a formal link
between the snake curve motion and an image motion.

\begin{proposition}\label{snakeequationsequivalence}
Assume that a function $(t,\fxx)\mapsto u(t,\fxx)$ is $C^2$ in a neighborhood of $(t_0,\fxx_0)$ and that
$Du(t_0,\fxx_0)\not=0$. Then $u$ satisfies the snake equation
\begin{equation}\label{snakeequationfunction}
\frac{\partial u}{\partial t}(t,\fxx)= g(\fxx){\rm curv}(u)(t,\fxx)|Du|(t,\fxx) +Dg(\fxx).Du(t,\fxx)
\end{equation}
in a neighborhood of $(t_0,\fxx_0)$ if and only if the normal flow (Definition \ref{defnormalflow}) of the level lines
of $u$ passing in this neighborhood satisfies the intrinsic snake equation,
\begin{equation} \label{snakeequationcurve}
\frac{\partial\fxx}{\partial t}(t,\fxy) =  g(\fxy) \bm{\kappa}(t,\fxy)-(Dg(\fxy).\fxn(t,\fxy)) \fxn(t,\fxy)),
\end{equation}
where $\bm{\kappa}(t, \fxx(t,\fxy))$ denotes the curvature vector of the level line of $u(t)$ passing by $\fxx(t,\fxy)$
and $\fxn(t,\fxy)$ one of its unit normals.
\end{proposition}

\begin{exercise} By imitating the proof of Corollary \ref{corollary5.9}, prove Proposition
\ref{snakeequationsequivalence}.

\finex
\end{exercise}



\section{Study of the snake equation}

We study in this section the  equation \eqref{snakeequationfunction}, which we can abbreviate as \fb& \frac{\partial
u}{\partial t} = g |Du| curv(u) + Dg . Du \label{imagesnakeequation} \fe& and shall call the {\it snake equation}.


\paragraph{Admissibility of the equation and uniqueness of solutions.}
Let us set
$$F(A,p,\fxx) = g(\fxx) A(p^\perp, p^\perp) +  Dg(\fxx).p$$
Equation \eqref{imagesnakeequation} can be obviously written as
$$ \frac{\partial u}{\partial t} = F(D^2u, Du, \fxx)$$
It is easy  checked that $F$ is admissible (see Definition \ref{definition15.1}).  As a consequence Theorem
\ref{theorem15.14} ensures uniqueness of viscosity solutions of  the equation \eqref{imagesnakeequation} for any
Lipschitz initial condition $u_0$.

\begin{exercise}
Check that $F(A,p,\fxx)$ is admissible.
 \finex \end{exercise}

\paragraph{Existence of solutions by approximation.}

Let us now construct an approximation scheme to the solution of Equation (\ref{imagesnakeequation}). It is possible to
construct a family of structuring elements having as asymptotic behavior the right hand term of the equation
(\ref{imagesnakeequation}). However, this term being a sum of two simple operators, it is simpler to associate with
each one of these operators a simple family of structuring elements and to alternate their corresponding filters. Note
that, due to the presence of  $g(\fxx)$ and $Dg(\fxx)$, the equation is not invariant by translation. As a consequence,
the families of structuring elements will depend on $\fxx$. This situation is new. We will need conditions ensuring
that $\inf$-$\sup$ operators with space-varying sets of structuring elements preserve Lipschitz constants.

\begin{lemma}\label{spacevaryingstructuringelement} Let $B\subset B(0,1)$ and  $B(\fxx)=g(\fxx)B$ be
 a space varying structuring element such that
 $g(\fxx)$ is $M$-Lipschitz. Consider the
associated space-varying dilation $Tu(\fxx)=\sup_{\fxz\in B(\fxx)}u(\fxx+\fxz).$  Then if $u(\fxx)$ is a $L$-Lipschitz
function, $Tu$ is a $LM$-Lipschitz function. The same result works with an $\inf$ instead of a  $\sup$.
\end{lemma}


\begin{proof}For every $\fxz\in B$, $\fxx$ and $\fxy$ one has
$$u(\fxx+g(\fxx)\fxz)   \leq u(\fxx+g(\fxy)\fxz) + L |(g(\fxx)-g(\fxy))\fxz|\leq u(\fxx+g(\fxy)\fxz)+ LM|\fxx-\fxy|.$$
Taking the $\sup$ on $B$ on both sides,

$$Tu(\fxx)=\sup_{\fxz\in B} u(\fxx+g(\fxx)\fxz)\leq Tu(\fxy)+ LM|\fxx-\fxy|.$$
Using this and the analogous inequality interchanging $\fxx$  and  $\fxy$ one gets the announced result.
\end{proof}

\begin{corollary} \label{Lipschitzconstantpreservation} Let $\setB$ be a family of structuring elements $B$ such that $B\subset B(0,1)$ and $g(\fxx)$ a
$M$-Lipschitz function. Let
$$Tu(\fxx)=\inf_{B\in \setB}\sup_{\fxy\in g(\fxx)B}u(\fxx+\fxy).$$
Then if $u$ is $L$-Lipschitz, $Tu$ is LM-Lipschitz.  The same result is true with a $\sup\inf$ instead of an
$\inf\sup$.
\end{corollary}

\begin{proof} This follows from Lemma \ref{spacevaryingstructuringelement} and the fact that an arbitrary infimum of
$L$-Lipschitz functions also is Lipschitz with the same constant (see exercise \ref{infdeLipschitz}.)
\end{proof}

\begin{exercise} \label{infdeLipschitz} Prove that if $(u_i)_{i\in I}$ is a family of $L$-Lipschitz functions such that $|u_i(0)|\leq C$ is
bounded independently of $i$, then
 then $u(\fxx)=\inf_{i\in I}u_i(\fxx)$ and $v(\fxx)=\sup_{i\in I}v(\fxx)$ also are $L$-Lipschitz.

\finex
\end{exercise}
Let us now define the space-varying structuring elements naturally associated with the snake equation.

 \subparagraph{Approximation of $-Dg Du$.} We consider the family made of a single element:
$$\setB_h(\fxx) = \{ \{ h  Dg(\fxx) \} \}$$
By Taylor formula we then have, for each point where $u$ is $C^2$:
$$ (S_h u)(\fxx) = \inf_{B \in \setB_h(\fxx)} \sup_{\fxy \in B} u(\fxx+\fxy) =
u(\fxx +h Dg) $$ $$= u(\fxx) + h Dg(\fxx).Du(\fxx) + O(\fxx, h^2),$$ where $O(\fxx, h^2)$ converges uniformly on every
compact set $K$ where $u$ is $C^2$. We can rewrite the last relation
\begin{equation}
 (S_h u)(\fxx) - u(\fxx) =  h
Dg(\fxx).Du(\fxx) + O(\fxx, h^2). \label{uniformconsistencySh}\end{equation}

\subparagraph{Approximation of $g |Du| curv(u)$.} We consider the structuring elements of the median filter (See
Chapter \ref{Chap9}):
$$\setB'_h(\fxx) = \{ B \mid B \subset B(0,\sqrt{6g(\fxx)h}) \hbox{ and } meas(B)
 \geq 3\pi g(\fxx)h)\}.$$
Set $$ (S'_h u)(\fxx) =  \sup_{B \in \setB'_h(\fxx)} \inf_{\fxy \in B} u(\fxx+\fxy).$$ Thanks to Theorem
\ref{theorem10.4}, we have for any   $u:\RR^2\to\RR$ which is $C^2$:

\begin{enumerate}[$(i)$]\setlength\itemsep{2pt}
\item On every compact set $K\subset\{\fxx\mid Du(\fxx)\not=0\}$,
\begin{equation}\label{consistencySprime}
S'_h u(\fxx)=u(\fxx)+ h g(\fxx) |Du(\fxx)|{\rm curv}(u)(\fxx)+O(\fxx,h^\frac 32),
\end{equation}
where $|O(\fxx,h^\frac 32)|\leq C_Kh^\frac 32$ for some constant $C_K$ that depends only on $u$ and $K$.
\item On every compact set $K$ in $\RR^2$,
\begin{equation}
|S'_h u(\fxx)-u(\fxx)|\leq C_Kh \label{uniformconsistencyShprime}
 \end{equation}
 where the constant $C_K$ depends only on
$u$ and $K$.
\end{enumerate}


\subparagraph{Alternating the two filters.} We now consider $T_h = S_h S'_h$, the alternate filter whose iteration
should  mimic the snake equation. Using \eqref{uniformconsistencySh}, \eqref{consistencySprime} and
\eqref{uniformconsistencyShprime}, Lemma \ref{alternatefiltersuniformconsistency} ensures that for any compact set $K$
where $|Du(\fxx)|\neq 0$,
$$ (T_h u) (\fxx) - u(\fxx) = h ( g(\fxx) Du(\fxx) curv(u)(\fxx) + Dg(\fxx) Du(\fxx) ) + O(\fxx, h^\frac 32) $$
and for $\fxx$ in any compact set  $K$:
$$ (T_h u) (\fxx) -u(\fxx) =  O(\fxx, h), $$
where in both cases the convergence of $O(\fxx, h)$ is uniform on $K$.  As a consequence the filter $T_h$ is uniformly
consistent (see Definition \ref{consistencywithequationgenerale}) with the PDE  \eqref{imagesnakeequation}.

\begin{exercise} Check  carefully that in the above argument, Lemma
\ref{alternatefiltersuniformconsistency} applies. \finex
\end{exercise}




\subparagraph{Construction of the approximate solutions.}
We consider a $L$-Lipschitz initial function $u_0$. We then
define
 $u_h(t, \fxx)$ for every $h>0$ by
$$
\forall n\in \natural,\quad u_h(nh, \fxx) = (T_h^nu_0)(\fxx).
$$
From now on we shall assume that

\begin{equation}\label{assumptionsong}  Dg(\fxx) \hbox{ and } \sqrt{6g(\fxx)} \hbox{ are 1-Lipschitz and bounded by 1}.
\end{equation}

\medskip  \noindent By Corollary \ref{Lipschitzconstantpreservation}, the 1-Lipschitz assumptions on $g$
and $\sqrt{6g(\fxx)}$ ensure that $(T_h)^nu_0$ is $L$-Lipschitz for every $n$. All of these bounds can be achieved by
simple scaling, namely multiplying $g$ by a small enough constant.

\subparagraph{Uniform continuity in $t$ of the approximate solutions.}

\begin{lemma} \label{uniformcontinuityofsnakeapproximation}
The approximate solutions of the snake equation $u_h(nh,\fxx)=(T_h^nu_0)(\fxx)$ are uniformly continuous in $t$. More
precisely:
$$ \forall t, \forall n\mid n h \leq t, \qquad -(L \sqrt{2t}+Lt) \leq T^n_h u -u \leq +(L\sqrt{2t}+Lt).$$
\end{lemma}

\begin{proof}
Let us bound the operator $T_h$ by two isotropic and translation invariant operators. For every  L-Lipschitz
function $u$, one has
$$u(\fxx) - L h |Dg|  \leq (S_{h} u)(\fxx) = u(\fxx +h Dg) \leq u(\fxx) + L h |Dg|$$
$$ u(\fxx) - Lh \leq  (S_{h} u)(\fxx)\leq u(\fxx)+ Lh$$
Then due to the fact that $S'_{h}(u+c) = S'_{h}(u)+c$ for any constant $c$, we also have \fb& (S'_h  u) -Lh \leq T_h u
= S'_{h} S_{h} u \leq (S'_{h} u) + Lh  \label{snake3}\fe&

Let us consider $v(\fxy) = L|\fxx-\fxy|$. The family $\setB'_{\sqrt{6 g(\fxx) h}}$ of structuring elements of the
filter $S'_h$ is made of the subsets of the disk of center $0$ and radius $\sqrt{6g(\fxx)h} \leq  \sqrt{h}$. It is easy
to check that for any $B \in \setB'_{\sqrt{6 g(\fxx) h}}$ satisfying $|B|\geq 3\pi g(\fxx)h$, there exists $B' \in
\setB'_{\sqrt{h}}$ satisfying $|B|\geq \frac{\pi h}2$, such that
$$\inf_{\fxy \in B} v(\fxx-\fxy) = \inf_{\fxy \in B'} v(\fxx-\fxy).$$
Thus,

$$(S'_h v)(\fxx) = \sup_{B \in \setB'_{\sqrt{6 g(\fxx) h}}}
\inf_{\fxy \in B} v(\fxx-\fxy) \leq \sup_{B \in \setB'_{\sqrt{h}}} \inf_{\fxy \in B'} v(\fxx-\fxy).$$ This yields

\fb& \forall \fxx, \;(S'_h v)(\fxx) \leq (M_{h} v)(\fxx) \label{snake1} \fe& where $M_h$ denotes the  median filter on
the ball $B(0, \sqrt{h})$, as defined in Chapter \ref{Chap9}. Similarly, for $w(\fxy) = -L|\fxx-\fxy|$, we have \fb&
\forall \fxx, \;(M_{h} w)(\fxx) \leq (S'_h w)(\fxx) \label{snake2}\fe& We deduce from (\ref{snake1}), (\ref{snake2})
and (\ref{snake3}) the  inequalities  \fb& (T_h v)(\fxx) \leq (M_{6   h}v)(\fxx) + Lh \qquad (M_{6 h}w)(\fxx) - Lh \leq
(T_h w)(\fxx) \fe& By monotonicity of $T_h$ and of the median operator, we thus have for all $n \in \entier$, and for
all $\fxx \in \real^2$,

\fb& (T_h^n v)(\fxx) \leq (M^n_{h}v)(\fxx) + nLh \qquad (M^n_{h}w)(\fxx) -n Lh \leq (T^n_h w)(\fxx) \label{snake4}\fe&
Now, since $u_0$ is $L$-Lipschitz, one has
$$  u(\fxx) - L|\fxx-\fxy| \leq u(\fxy) \leq u(\fxx) - L|\fxx-\fxy|.$$
Thus
$$ u(\fxx) + (T^n_h w)(\fxx) \leq (T^n_h u)(\fxx) \leq u(\fxx) +  (T^n_h v)(\fxx)$$
Using (\ref{snake4}), we obtain
$$ (M^n_{h} w)(\fxx)  -nLh\leq (T^n_h u)(\fxx) -u(\fxx) \leq  (M^n_{h} v)(\fxx)+nLh$$


Lemma \ref{controledumedianitere} tells us that for $h$ small enough and for $nh\leq t$, one has
$$(M^n_{h} v)(\fxx) \leq L \sqrt{2t}.$$
 An analogous inequality obviously holds for $w$, so that for $h$ small enough
$$ \forall t, \forall n; n h \leq t \qquad -(L \sqrt{2t}+Lt) \leq T^n_h u -u \leq +(L\sqrt{2t}+Lt)$$
\end{proof}

\begin{exercise} Give all details for the proof of the property used in the above proof:
For any $B \in \setB'_{\sqrt{6 g(\fxx) h}}$ satisfying $|B|\geq 3\pi g(\fxx)h$, there exists $B' \in \setB'_{\sqrt{h}}$
satisfying $|B|\geq \frac{\pi h}2$, such that $\inf_{\fxy \in B} v(\fxx-\fxy) = \inf_{\fxy \in B'} v(\fxx-\fxy).$
\finex \end{exercise}
\subparagraph{Convergence of the approximate solutions}

 \begin{theorem}
Let  $g$ be a $C^2$  function which is zero outside a ball $B(0,R)$ and satisfies the bounds \eqref{assumptionsong}.
Then for every Lipschitz function $u_0\in \mathcal F$, there exists a unique viscosity solution $u(t,\fxx)$ of the
snake equation
$$\frac{\partial u}{\partial t} = F(D^2u, Du, \fxx) = g |Du| curv(u) + Dg . Du  \qquad u(0,\fxx) = u_0(\fxx).
$$
In addition, $u(t,\fxx)$ is Lipschitz in $\fxx$ and holderian in $t$ and when $h$ tends to 0 and $nh \rightarrow t$,
$(T^n_h u_0)(\fxx)$ converges towards $u(t,\fxx)$ uniformly on compact sets of $\RR^+\times \RR^2$.
\end{theorem}

\begin{proof}
The operator $T_h$ is monotone and local (and therefore satisfies the uniform local comparison principle.) It is
uniformly consistent with the PDE (\ref{imagesnakeequation}) and commutes with the addition of constants. By Corollary
\ref{Lipschitzconstantpreservation} and  Lemma \ref{uniformcontinuityofsnakeapproximation}, its associated approximate
solutions $h \rightarrow u_h(t, \fxx)$ are $L$-Lipschitz in $\fxx$ and uniformly Holderian in $t$ for any initial
Lipschitz function $u_0$. Thus, using Ascoli-Arzela theorem, there is a sub-sequence of the sequence $h \rightarrow
u_h$ which is uniformly converging on every compact set towards a function $u(t,\fxx)$. By Proposition
(\ref{limitedeuhestsolution}), this implies that $u$ is a viscosity solution of (\ref{imagesnakeequation}). In other
words, we get the existence of a viscosity solution for any initial Lipschitz function $u_0$. Since this solution is
unique, all subsequences of $u_h$ converge to the same function $u$ and therefore the whole sequence $u_h$ converges to
$u$.
\end{proof}




\section{Back to shape evolution}

 We now consider the operator $T_t$ which associates to any Lipschitz function $u_0$ in $\mathcal F$ the unique
viscosity solution  $u(t,.)$ of (\ref{imagesnakeequation}) with initial condition $u_0$. $T_t$ is clearly a monotone
operator as limit of monotone operators.  According to Proposition \ref{proposition16.2}, $T_t$ is also contrast
invariant.
 Let us check that it also is standard monotone. We have assumed that $g$ is
zero outside a ball $B(0,R)$. It is a straightforward deduction that $T_h^nu_0(\fxx)=u_0(\fxx)$ for every $\fxx$ with
norm larger than $R$. Since $T_h^nu_0$ converges uniformly to $T_tu_0$ on compact sets as $nh\to t$ and $n\to \infty$,
we still have $T_tu_0(\fxx)=u_0(\fxx)$ outside $B(0,R)$ and therefore $T_tu_0(\infty)=u_0(\infty)$. Thus $T_tu_0$
belongs to $\mathcal F$ and  $T_t$ is standard monotone.

\begin{exercise}
Check in detail the above two statements, that $T_t$ is monotone and contrast invariant. The second statement can be
proven as indicated by using Proposition \ref{proposition16.2}, but also directly by using the contrast invariance of
the iterated operators $T_h^n$ which converge to $T_t$ as $nh\to t$.   \finex
\end{exercise}

\begin{proposition} By direct application of the level set extension theorem \ref{evansspruckextension},
the monotone and contrast invariant image operator $T_t$, defined for any initial Lipschitz function $u_0$ in $\mathcal
F$, defines a unique set operator $\mathcal T_t$ on defined on the set $\mathcal L$ of the compact sets of $S_N$.  Then
$\mathcal T_t$ is monotone,
    $T_t$ and $\mathcal T_t$
   satisfy the commutation with thresholds
$ \mathcal T_t( \levelset_{\lambda}u)=\levelset_{\lambda}(T_tu) $ for all $\lambda\in\RR$,  $T_t$ is the stack filter
associated with $\mathcal T_t$ and $\mathcal T_t$ is upper semicontinuous on $\mathcal L$. In addition, since $T_t$ is
standard, so is $\mathcal T_t$.
\end{proposition}


\begin{exercise}
Theorem \ref{evansspruckextension} applies to an operator $T$ defined on $\mathcal F$. Now, we have defined $T_t$ on
the Lipschitz functions of $\mathcal F$ only. In order to show that this is not a problem, prove first that any
function $u_0$ in $\mathcal F$ can be approximated uniformly by a sequence of functions $u_n$ which are $C^1$ and
 Lipschitz. Then show that that $T_tu_n$ is a Cauchy sequence for the uniform convergence and conclude that
$T_t$ can be extended into a contrast invariant standard monotone operator on all of $\mathcal F$. \finex\end{exercise}

 \noindent {\bf The snake algorithm}

Let us now see how we can the  above results to define a curve evolution. Consider a closed curve $C=\fxx(s)$
surrounding a compact set $K$ of $\real^2$. We define the generalized ``curve'' evolution of $C$ by the following
algorithm:
\begin{description}
\item{Step 1}  Construct a Lipschitz function $u_0$ so that:
\begin{itemize}
\item $\levelset_0 u_0 = K$
\item $u_0$ is Lipschitz.
\end{itemize}
Such a function $u$ can be obtained by considering the signed distance function to the set $K$ defined by $u(\fxx)=
dist(\fxx,K^c)$ if $\fxx\in K$ and $u(\fxx)=\max(-dist(\fxx, K), -1)$ if $\fxx\in K^c$.
\item{Step 2} Compute the viscosity solution $u(t,\fxx)$ of equation (\ref{imagesnakeequation}) with initial condition
 $u_0$. We know one way to do it, by computing $u_h(nh, \fxx)=(T_h^nu_0)(\fxx)$.
\item{Step 3} Set $K(t) = \levelset_0 u(t,.)$ and $C(t) = K(t) \cap \bar{K(t)^c}$.
 ($C(t)$ is the topological boundary of $K(t)$.)
\end{description}
According to all preceding considerations, this algorithm defines for any curve $C$ and for any $t \geq 0$ a unique set
of points $C(t)$. The evolution of $C(t)$ is independent from the choice of the initial function $u_0$ and corresponds
to a generalization of the curve evolution PDE(\ref{equationcurvesnake}). The preceding algorithm satisfies the {\it
shape inclusion principle}: If initially $C_1$ and $C_2$ are two curves such that $C_1$ surrounds $C_2$, which means
$K_2\subset K_1$ then $C_1(t)$ surrounds $C_2(t)$ : $K_2(t)\subset K_1(t)$.

However, $C(t)$ is not necessarily a curve of $\real^2$. It simply is the boundary of a set $K(t)$. It is therefore
difficult to check if the initial geodesic snake energy estimated on $C(t)$ is decreasing. This problem is open! We
have defined a very robust and weak evolution of curves and the initial $C$ needs not be even a curve. Any set of
curves or more generally any closed set can be taken as initial datum. The very good point of this generalized curve
evolution is that it allows $C(t)$ to  break into pieces surrounding separated shapes, as illustrated in Figure
\ref{snaketwoballs}.

\section{Implementation by finite difference scheme}
As usual, the approximation of the process by iterated $\inf\sup$ filter is not quite satisfactory, because these
filters fail to be consistent with the equation at small scales, as pointed out for the median filter. Before proposing
another way to simulate the  snake equation,  some heuristic comments  of the snake equation behavior will be useful.
$$
\frac{\partial u}{\partial t} = g |Du| curv(u) +Dg . Du
$$
The first term is the well known mean curvature motion. As we have seen, it tends to shrink the level lines towards
points. The speed of this motion is related to the amplitude of $g$. On an edge, $g$ is small,  but not zero. Thus, the
motion is slowed down, but does not stop.

The second term is the erosion term. It tends to move the level lines of $u$ downwards for $g$, that is, towards the
edges of $I$ (see figure \ref{figureconvectionsnake}), creating therefore shocks for $u$ around edges since level lines
of $u$ converge to them on both sides. Contrarily to the first term, this term is not active on flat regions for $g$.
Even worse, due to noise, little gradients for  $I$ will induce a non negligible variation of amplitude of $g$,
resulting in non negligible $Dg$ term with random direction. In others words, on flat regions, one can expect to
observe random perturbations of the shape of the evolving contour.




\begin{figure}[!ht]
\centerline{\hbox{ \epsfxsize=10cm \epsfbox{snake.eps} }} \caption[Convection term of the active contour
equation]{\label{snakefigure}  Convection term of the active contour equation. The convection term of the active
contour equation tends to create around minima of $g$. Indeed, the level lines of $u$ are moved in the direction
opposite to the gradient of $g$.\label{figureconvectionsnake}}
\end{figure}

More precisely, we have:
\begin{itemize}
\item{\bf near an edge:}
There is a risk that the curvature term pushes the level line over the edge, since it tends to shrink the curve. The
second term instead moves the level line towards the edge. Thus the  effects of both terms can  be opposed.
Modifications should be made in the equation to ensure that the second term always  wins.

\item{\bf far from an edge:}
The first term moves the level line fast, since $g(\fxx)$ is high. The second term attracts the level line towards tiny
edges  of the image, thus creating little shocks. Here again, modifications should be made in the  equation to let the
second term win.

\end{itemize}

Even if the equation does not show any  parameter weighting the two terms, a  weight between them is in fact hidden in
the choice of the function $g$. Finding a function $g$ that makes the correct equilibrium both near an edge and far
from an edge is somehow complex and case sensitive.

Assuming such an equilibrium could be found, relying on the single (weighted) mean curvature motion to shrink the level
line is not a good idea. Indeed, a weighted mean curvature motion will never help transforming a convex level line into
a general non-convex one. If we start with a circle, it is impossible to recover (e.g.) a star, as illustrated in
Figure \ref{snakestarfigure}.

\begin{figure}
\centerline{\hbox{ \epsfxsize=10cm \epsfbox{star.eps} }} \caption[A difficulty : the local minima of the active contour
energy]{\label{snakestarfigure} A difficulty : the local minima of the active contour energy. Left image: assume that
$g$ is  null on the shape (drawn in bold) and that the initial contour is the line \#0. From the initial line to the
contour of the shape, the intermediate state \#1 consisting of the convex hull of the polygon shows a smaller energy
than the intermediate state \#2 (drawn on the right). This illustrates the difficulties arising with the snake equation
when we wish to land the active contour onto concave parts of the desired contour.
 }
\end{figure}


To cope with this problem, we shall  add  an extra term to the equation. This term is a classical erosion with the same
weight as  the curvature motion. In the experiments we therefore considered the modified snake equation
\begin{equation}
\frac{\partial u}{\partial t} = g |Du| (curv(u)-1) + Dg . Du.
\label{modifiedsnakeequation}
\end{equation}
 This new term is added under the assumption that the initial contour $C(0)$
has been drawn on the outside of the shape and that the initial function $u_0$ is positive inside the contour. Then the
term $-1$ tends to erode the contour. If instead the initial contour was drawn inside the shape, then the new term
should be $+1$. By an obvious adaptation of the proof of proposition \ref{???????} one sees that the associated curve
evolution is
\begin{equation} \label{modifiedsnakeequationcurve}
\frac{\partial\fxx}{\partial t}(t,s) =  g(\fxx(t,s)) \bm{\kappa}(\fxx(t,s))-(Dg(\fxx(t,s)).\fxn(\fxx(t,s)))
.\fxn(\fxx(t,s))).
\end{equation}

  Returning to the snake equation, our finite difference scheme will the following:

\begin{eqnarray*}
% \nonumber to remove numbering (before each equation)
  u^{n+1}(\fxx) = && u^{n}(\fxx) + dt ( g(\fxx) (E_1(u^n)(\fxx)+M_{\sqrt{6}}(u^n)(\fxx) -2u^n(\fxx))) \\
   & & + dt( u^n(\fxx + Dg(\fxx)) - u^n(\fxx)),
\end{eqnarray*}
where $E_1u(\fxx)=\inf_{\fxy\in B(0,1)}u(\fxx+\fxy)$ denotes the erosion by the unit ball $B(0,1)$ and $M_{\sqrt{6}}$
is the median filter on the ball $B(0,\sqrt{6}$. The scheme is not contrast invariant, but is maximum decreasing and
minimum increasing, provided $|2g(\fxx)+1|\leq 1$.

\begin{exercise}
Prove the last statement, that the above scheme is maximum decreasing and minimum increasing, provided
$|2g(\fxx)+1|\leq 1$.  Prove that the  scheme is consistent, in the sense that if the pixel size $h$ tends to zero then
the second member of the equation tends to the second member of the modified snake equation
\ref{modifiedsnakeequation}.  \finex
\end{exercise}


Figure \ref{birdsnake} illustrates the extraction of the bird shape on a textured background.  This  experiment
illustrates well the complexity of the figure-background problem: the shape of the bird body has a quickly changing
color from the white head to its dark tail. The background being uniform grey, there is no unique level line
surrounding the whole shape. In fact, the gradient of the contour of the bird vanishes at many points. The fact that we
``see'' this contour is a classical illusion, called {\it subjective contour}. To uncover the illusion, the reader
should scan small parts of the shape contour by using a white sheet with a small hole. Then he or she will realize that
the contour seen globally has no complete numerical local evidence. This  observation implies that no classical edge
detection device would  give out the whole contour. This can be checked by applying a Canny edge detector to the shape.
To some extent, the snake method manages instead to surround the body shape. All the same, there is a risk that,
because of the erosion term, the active contour goes through the subjective contour. If instead the curvature term is
too strong, it can stop the contour before it reaches a concave corner of the shape. These facts explain the obvious
inaccuracy and irregularity of the found contour.

\begin{figure}
\centerline{\hbox{
\begin{tabular}{ccc}
\epsfxsize=4cm \epsfbox{bird.ps} &\epsfxsize=4cm \epsfbox{label.ps} &\epsfxsize=4cm \epsfbox{snake_bird.ps}
\end{tabular}
}} \caption[Silhouette  of a bird by active contour.]{ Silhouette of a bird by active contour. Left: original image,
middle: initial contour, right: final contour (steady state of the snake equation).  The contours of the bird body are
partly subjective.  The snake evolution manages to some extent to find them, but tends to indent the subjective
contours and to round to the concave corners of the shape. \label{birdsnake}}
\end{figure}


\begin{figure}
\centerline{\hbox{
\begin{tabular}{ccc}
\epsfxsize=4cm \epsfbox{boules.ps} &\epsfxsize=4cm \epsfbox{label.ps} &\epsfxsize=4cm
\epsfbox{inter1.ps} \\
\epsfxsize=4cm \epsfbox{inter2.ps} &\epsfxsize=4cm \epsfbox{inter3.ps} &\epsfxsize=4cm \epsfbox{inter4.ps}
\end{tabular}
}} \caption[Active contour with topological change.]{ Active contour with topological change. Top, left: original
image, middle: initial contour, right: intermediate state. Down, left and middle: successive intermediate states,
down-right: final contour (steady state).  This experiment shows that the level lines of $u$ which bound the evolving
contour cannot be classical solutions of the modified curve evolution equation \eqref{modifiedsnakeequationcurve}.
Indeed, the motion generates singularities when the contour splits. The generalized evolution provided by the viscosity
solution of the snake equation yields more flexibility and allows an evolving curve to split. Original image is "Vue
d'esprit 3", by courtesy of {\it e-on software}.\label{snaketwoballs}}
\end{figure}

\section{Exercises}

\begin{exercise}
In the whole exercise $u(t,\fxx)$  and $\fxx(t)$ are supposed as smooth as needed to  make the computations. Our  aim
is to interpret the equation \begin{equation}\frac{\partial u}{\partial t}=Dg.Du\label{propagationpure}\end{equation}
as a motion of the level lines of $u$ towards the minima of $g$. Let us consider a point $\fxx(t)$ on a level line of
$u(t,\fxx)$ with level $\lambda$ and denote by $\fxx'(t)$ the motion vector of $\fxx(t)$ in the direction normal to the
level line. We know that $\fxx(t)$ obeys the normal flow equation \eqref{eq66x}.     Deduce from this equation and
\eqref{propagationpure} that $\fxx(t)$ moves downwards in the landscape given by $g$, that is $g(\fxx(t))$ is a
non-increasing function of $t$. \finex
\end{exercise}

\begin{exercise}
Construction of another inf-sup scheme converging towards the viscosity solution of the equation
(\ref{imagesnakeequation}). Consider the  family of structuring elements
$$
\setB_h(\fxx) = \{ B \mid B \subset B(\fxx + hDg(\fxx), \sqrt{6g(\fxx)h}) \hbox{ and } meas(B) \geq 3\pi g(\fxx) h\}.
$$
and the operator
$$ T_h u(\fxx) = \inf_{B \in \setB_h} \sup_{\fxy \in B} u(\fxx + \fxy) $$

\begin{enumerate}
\item Interpret the operator $T_h$ as a shifted median filter.
\item Show that $T_h$ is uniformly consistent with equation (\ref{imagesnakeequation}).
\item Show that the iteration of $T_h$ converges towards a viscosity solution of  (\ref{imagesnakeequation}).
\end{enumerate}
\finex
\end{exercise}

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