%gmChap11.tex.

\graphicspath{{./figures/}}

\chapter[Asymptotic Behavior in $\mathbb R^N$]{Asymptotic
Behavior in Dimension N\label{ChapasymptoticdimensionN}}

We are going to generalize to $N$ dimensions the asymptotic results of Chapter \ref{Chap10}. Our aim is to show that
the action of any local $\SMTCII$ operator, when properly scaled, is a motion of the $N$-dimensional image that is
controlled by its principal curvatures.  In particular, we will relate the median filter to the mean curvature of the
level surface.
\section{Asymptotic behavior theorem in $\RR^N$}

 Let $u:\RR^N\to\RR$ be $C^3$ and assume
that $Du(\fxx)\not=0$. Then we  denote the vector whose terms are
the $N-1$ principal curvatures of the level surface $\{\fxy\mid
u(\fxy)=u(\fxx)\}$ that passes through $\fxx$ by ${\bm
\kappa}(u)(\fxx)={\bm \kappa}(u)=(\kappa_2,\ldots,\kappa_N)$. The
terms $\kappa_i(u)(\fxx)|Du(\fxx)|$ are then the eigenvalues of
the restriction of $D^2u(\fxx)$ to $Du(\fxx)^\perp$.  (See
Definition \ref{definition4.19}.) For $\fxx\in\RR^N$, we write
$\fxx=(x,y_2,\ldots,y_N)=(x,\fxy)$, $\fxy\in\RR^{N-1}$ and in the
same way $\fxs=(s_2,\ldots,s_N)$.

\begin{theorem}\label{theorem11.1} Let $T$ be a local $\SMTCII$ operator.
Define
\begin{equation}\label{eq121}
H(\fxs)=T[x+s_2y_2^2+\cdots+s_Ny_N^2](0).
\end{equation} Then for every $C^3$
function $u:\RR^N\to\RR$,
\begin{enumerate}[$(i)$]\setlength\itemsep{2pt}

\item $T_hu(\fxx)=u(\fxx)+hH(0)|Du(\fxx)|+O(\fxx, h^2); $

\item If $H(0)=0$, then on every compact set $K$ contained in
$\{\fxx\mid Du(\fxx)\not=0\}$
\[
T_hu(\fxx)=u(\fxx)+hH\Bigg(h\frac1{2} {\bm
\kappa}(u)(\fxx)\Bigg)|Du(\fxx)|+O(\fxx, h^3)
\]
where $|O(\fxx, h^3)|\leq C_Kh^3;$

\item If $H(0)=0$, then on every compact set $K\subset\RR^N$,
\[
|T_hu(\fxx)-u(\fxx)|\leq C_Kh^2,
\]
 where  $C_K$ denotes some constant that depends only on $u$ and $K$.
\end{enumerate}
\end{theorem}

\begin{proof} The proof is the same as the proof of Theorems \ref{theorem10.1} and
\ref{theorem10.2}. We simply have to relate the notation used for
the $N$-dimensional case to that used in the two-dimensional case.
We begin by assuming that $Du(\fxx)\not=0$. We then establish the
local coordinate system at $\fxx$ defined by ${\bm
i}_1=Du(\fxx)/|Du(\fxx)$ and ${\bm i}_2,\ldots,{\bm i}_N$, where
${\bm i}_2,\ldots,{\bm i}_N$ are the eigenvectors of the
restriction of $D^2u(\fxx)$ to the hyperplane $Du(\fxx)^\perp$.
Then in a neighborhood of $\fxx$ we can expand $u$ as follows:
\begin{equation}\label{eq123}
u(\fxx+\fxy)=u(\fxx)+px+ax^2+b_2y_2^2+\cdots+b_Ny_N^2+(\fxc\cdot\fxy)x+R(\fxx,\fxy),
\end{equation}
where $\fxy=x{\bm i}_1+y_2{\bm i}_2+\cdots+y_N{\bm i}_N$, $p=|Du(\fxx)|>0$, and for $j=2,\ldots,N$,
\begin{equation}
\begin{aligned}
a&=\frac{1}{2}\frac{\partial^2 u}{\partial x^2}(\fxx)=\frac{1}{2}D^2u(\fxx)({\bm i}_1,{\bm i}_1),\\
b_j&=\frac{1}{2}\frac{\partial^2 u}{\partial y_j^2}(\fxx)=\frac{1}{2}D^2u(\fxx)({\bm i}_j,{\bm i}_j),\\
c_j&=\frac{\partial^2u}{\partial x\partial y_j}(\fxx)=D^2u(\fxx)({\bm i}_1,{\bm i}_j).
\end{aligned}
\end{equation}
We can also write $b_j$ as
\begin{equation}\label{eq125}
b_j=\frac1{2}|Du(\fxx)|\kappa_j(u)(\fxx).
\end{equation}

For the proof of $(i)$, we write
$u(\fxx+\fxy)=u(\fxx)+px+O(\fxx,|\fxy|^2)$  and just follow the
steps of the proof of Theorem \ref{theorem11.1}. The proof of
$(ii)$ and $(iii)$ follows, step by step, the proof of Theorem
\ref{theorem10.2}. We need only make the following
identifications: $cxy\leftrightarrow (\fxc\cdot\fxy)x$,
$by^2\leftrightarrow b_2y_2^2+\cdots+b_Ny_N^2$, and ${\rm
curv}(u)\leftrightarrow\kappa(u)$.
\end{proof}

\begin{figure}[!t]
\setkeys{Gin}{width=12cm} \centerline{
\includegraphics{fig11.01.1.eps}
} \caption[Three-dimensional median filter]{Three-dimensional
median filter. The original three-dimensional image (not shown) is
of 20 slices of a vertebra. Three successive slices are displayed
in the left column. The next column shows their level lines
(multiples of 20). The third column shows these three slices after
one iteration of the median filter based on the three-dimensional
ball of radius two. The resulting level lines are shown in the
last column. \label{figure11.1} }
\end{figure}

\section{Asymptotic behavior of median filters in $\RR^N$}
The action of median filtering in three dimensions is illustrated
in Figures \ref{figure11.1} and \ref{figure11.2}. The median
filters we consider  will be defined in terms of a continuous
weight function $k:\RR^N\to[0,+\infty)$ that is radial,
$k(\fxx)=k(|\fxx|)$, and that is normalized,
$\int_{\RR^N}k(\fxx)\,\rmd\fxx=1$. Recall that, by definition,
\[
|B|_k=\int_B k(\fxx)\,\rmd\fxx.
\]
We also assume that $k$ is nonseparable, which is the case if
$\{\fxx\mid k(\fxx)>0\}$ is connected. Then by Proposition
\ref{egalitedemedianmoinsetplus}, $\fcnM_ku=\fcnM_k^-u$ and  the
 median operator can defined by
\begin{equation}\label{eq126}
\fcnM_ku(\fxx)=\inf_{|B|_k=1/2}\sup_{\fxy\in\fxx+B}u(\fxy).
\end{equation}


Define the scaled weight function $k_h$,  $0<h\leq1$, by
$k_h(\fxx)=h^{-N}k(\fxx/h)$. Then a change of variable shows that
$|B|_k=1/2$ if and only if $|hB|_{k_h}=1/2$, and this implies that
$(\fcnM_k)_h=\fcnM_{k_h}$ (see Exercise \ref{exercise10.5}). Since
we consider only one weight function at a time, there should be no
confusion if we write $\fcnM_h$ for the scaled operator.

We analyzed the asymptotic behavior of a  median filter in $\RR^2$
whose weight function was the characteristic function of the unit
disk in Chapter \ref{Chap10}. This proof can be generalized to
$\RR^N$ by taking $k$ to be the normalized characteristic function
of the unit ball. We will go in a different direction by taking
smooth weight functions. Our analysis will not be as general as
possible because this would be needlessly complicated. The $k$ we
consider will be smooth ($C^\infty$) and have compact support.
This means that the considered median filters are local.  Thus,
the results of Theorem \ref{theorem11.1} apply, provided we get an
estimate near  0 of the structure function  $H$  of the  median
filter.

\begin{lemma}\label{lemma11.2} Let $k$ be a nonnegative radial
function belonging to the Schwartz class $\mathcal S$. Assume
 that $\int_{\RR^N} k(\fxx)\,\rmd\fxx=1$ and that the support of $k$
  is connected in $\RR^N$. Then the structure function of $\fcnM_k$
  $H(h\fxb)=\fcnM_k[x+h(b_2y_2^2+\cdots+b_Ny_N^2)](0)$ can be expressed as
\[
H(h\fxb)=hc_k\Big(\sum_{j=2}^Nb_j\Big)+O(h^2),
\]
where
\[
c_k=\frac{\int_{\RR^{N-1}}y_2^2k(\fxy)\,\rmd\fxy}{\int_{\RR^{N-1}}k(\fxy)\,\rmd\fxy},
\]
$\fxy=(y_2,\ldots,y_N)$, and $\fxb=(b_2,\ldots,b_N)$.
\end{lemma}

\begin{proof}Before beginning the proof, note that  we have not
assumed that $k$ has compact support, so the result applies to the
Gaussian, for example.

We will use the abbreviation
$\fxb(\fxy,\fxy)=b_2y_2^2+\cdots+b_Ny_N^2$, since $\fxb$ is, in
fact, a diagonal matrix. Our proof is based on an analysis of the
function $f(\lambda,
h)=|\levelset_\lambda(x+h\fxb(\fxy,\fxy))|_k$. Since
$\levelset_\lambda(x+h\fxb(\fxy,\fxy))=\{(x,\fxy)\mid
x+h\fxb(\fxy,\fxy)\geq\lambda\}$, we can express $f$ as an
integral,
\[
f(\lambda,h)=\int_{\RR^{N-1}}\int_{\lambda-h\fxb(\fxy,\fxy)}^\infty k(x,\fxy)\,\rmd x\,\rmd\fxy.
\]
It follows from the assumption  that $k$ is in the Schwartz class
that $f:\RR^2\to\RR$ is bounded and $C^\infty$. Also, for every
$h\in\RR$, $\lim_{\lambda\to-\infty}f(\lambda,h)=1$ and
$\lim_{\lambda\to+\infty}f(\lambda,h)=0$. Thus, for every
$h\in\RR$, there is at least one $\lambda$ such that
$f(\lambda,h)=1/2$. In fact, there is only one such $\lambda$;
this is a consequence of the assumption that the $k$ is continuous
and that its support is connected, which implies that it is
nonseparable (see Exercise \ref{exercise9.2}). To see that
$\lambda$ is unique, assume that there are $\lambda<\lambda'$ such
that $f(\lambda,h)=1/2$ and $f(\lambda',h)=1/2$. Then the two sets
$\{(x,\fxy)\mid x+h\fxb(\fxy,\fxy)\geq\lambda'\}$ and
$\{(x,\fxy)\mid x+h\fxb(\fxy,\fxy)\leq\lambda\}$ both have
$k$-measure $1/2$, but their intersection is empty. This
contradicts the fact that $k$ is nonseparable. This means that the
relation $f(\lambda,h)=1/2$ defines implicitly a well-defined
function $h\mapsto\lambda(h)$.

Recall that $\fcnM_k$ was originally defined in terms of the
superposition formula
\[
\fcnM_ku(\fxx)=\sup\{\lambda\mid\fxx\in\setM_k\levelset_\lambda u\}.
\]
This translates for our case into the relation
\[
\fcnM_k[x+h\fxb(\fxy,\fxy)](0)=\sup\{\lambda\mid0\in\setM_k\levelset_\lambda [x+h\fxb(\fxy,\fxy)]\}=\lambda(h)
\]
because  $0\in\setM_k\levelset_\lambda [x+h\fxb(\fxy,\fxy)]$ if
and only if $|\levelset_\lambda[x+h\fxb(\fxy,\fxy)]|_k\geq1/2$.

We are interested in the behavior  of $h\mapsto\lambda(h)$ near
the origin. The first thing to note is that $\lambda(0)=0$. To see
this, write
\[
f(\lambda(0),0)=\int_{\RR^{N-1}}\int_{\lambda(0)}^\infty k(x,\fxy)\,\rmd x\,\rmd\fxy=\frac1{2}.
\]
Since $k$ is radial, the value  $\lambda=0$ solves the equation
$\int_{\RR^{N-1}}\int_{\lambda}^\infty k(x,\fxy)\,\rmd
x\,\rmd\fxy=1/2$. We have just shown that this equation has a
unique solution, so $\lambda(0)=0$.

\begin{figure}[!t]
\setkeys{Gin}{width=12cm} \centerline{
\includegraphics{fig11.02.1.eps}
} \caption[Median filtering of a three-dimensional image]{Median
filtering of a three-dimensional image. The first image is a
representation of the horizontal slices of a three-dimensional
level surface of the three-dimensional image of a vertebra. Right
to left, top to bottom: 1, 2, 5, 10, 20, 30, 60, 100 iterations of
a three-dimensional median filter based on a ball with radius
three. This scheme is a possible implementation of the mean
curvature motion, originally proposed as  such by Merriman, Bence
and Osher.\label{figure11.2}}
\end{figure}

Now consider the first partial derivatives of $f$:
\begin{align}
\frac{\partial f}{\partial\lambda}(\lambda,h)&=-\int_{\RR^{N-1}}k\Big(\big((\lambda-h\fxb(\fxy,\fxy))^2+\fxy\cdot\fxy\big)^{1/2}\Big)\,\rmd\fxy.\label{eq11a}\\
\frac{\partial f}{\partial h}(\lambda,h)&=\int_{\RR^{N-1}}\fxb(\fxy,\fxy)k\Big(\big((\lambda-h\fxb(\fxy,\fxy))^2+\fxy\cdot\fxy\big)^{1/2}\Big)\,\rmd\fxy.\label{eq11b}
\end{align}
These functions are $C^\infty$ because  $k$ is in the Schwartz
class; also, $(\partial f/\partial\lambda)(0,0)\not=0$. Then by
the implicit function theorem, we know that the function
$h\mapsto\lambda(h)$ that satisfies $f(\lambda(h),h)=1/2$ is also
$C^\infty$ and that
\[
\lambda'(h)\frac{\partial f}{\partial\lambda}(\lambda(h),h)+\frac{\partial f}{\partial h}(\lambda(h),h)=0.
\]
Thus, for small $h$,
\[
\lambda'(h)=-\frac{\frac{\partial f}{\partial h}(\lambda(h),h)}{\frac{\partial f}{\partial\lambda}(\lambda(h),h)},
\]
and, using equations \eqref{eq11a} and\eqref{eq11b}, we see that
\[
\lambda'(0)=\frac{\int_{\RR^{N-1}}\fxb(\fxy,\fxy)k\big((\fxy\cdot\fxy)^{1/2}\big)\,\rmd\fxy}{\int_{\RR^{N-1}}k\big((\fxy\cdot\fxy)^{1/2}\big)\,\rmd\fxy}.
\]
Now expand $\lambda$ for small $h$:
\[
\lambda(h)=\lambda(0)+\lambda'(0)h+O(h^2).
\]
Since
$\int_{\RR^{N-1}}\fxb(\fxy,\fxy)k\big((\fxy\cdot\fxy)^{1/2}\big)\,\rmd\fxy=\big(\sum_{j=2}^{N-1}b_j\big)\int_{\RR^{N-1}}y_2^2
k\big((\fxy\cdot\fxy)^{1/2}\big)\,\rmd\fxy$,
$H(h\fxb)=\lambda(h)$, and $\lambda(0)=0$, this proves the lemma.
\end{proof}

\begin{theorem}\label{theorem11.3}Let $k$ be a nonnegative radial function
belonging to the Schwartz class $\mathcal S$. Assume that
$\int_{\RR^N} k(\fxy)\,\rmd\fxy=1$ and that the support of $k$ is
compact and connected. Then for every $C^3$ function
$u:\RR^N\to\RR$:
\begin{enumerate}[$(i)$]\setlength\itemsep{2pt}
\item On every compact set $K\subset\{\fxx\mid Du(\fxx)\not=0\}$,
\[
\fcnM_hu(\fxx)=u(\fxx)+h^2\frac1{2}c_k\Big(\sum_{i=2}^N\kappa_i(u)(\fxx)\Big)|Du(\fxx)|+O(\fxx,h^3),
\]
where $|O(\fxx,h^3)|\leq C_Kh^3$ for some constant that depends
only on $u$ and $K$.

\item On every compact set $K\subset\RR^N$,
$|\fcnM_hu(\fxx)-u(\fxx)|\leq C_Kh^2$ for some constant $C_K$ that
depends only on $u$ and $K$.
\end{enumerate}
\end{theorem}

\begin{proof}   Theorem \ref{theorem11.1} is directly applicable.
We know from Lemma \ref{lemma11.2} that $H(0)=0$, so we can read
$(ii)$ directly from Theorem \ref{theorem11.1}$(iii)$. By Lemma
\ref{lemma11.2},
\[
H(h{\bm\kappa}(u))=hc_k\Big(\sum_{i=2}^N\kappa_i(u)|Du|\Big)+O(h^2).
\]
From this and Theorem \ref{theorem11.1}$(ii)$, we get
\[
\fcnM_hu(\fxx)=u(\fxx)+h^2\frac1{2}c_k\Big(\sum_{i=2}^N\kappa_i(u)(\fxx)\Big)|Du(\fxx)|+O(\fxx,h^3),
\]
and we know that the estimate is uniform on any compact set
$K\subset\{\fxx\mid Du(\fxx)\not=0\}$.
\end{proof}


\section{Exercises : other motions by  the principal curvatures}

This section contains several applications of Theorem
\ref{theorem11.1} in three dimensions. A level surface of a $C^3$
function in three dimensions has two principal curvatures, and
this provides an extra degree of freedom for constructing
contrast-invariant operators based on curvature motion. We develop
the applications in three exercises. For each case, we will assume
that the principal curvatures $\kappa_1$ and $\kappa_2$ are
ordered so that $\kappa_1\leq\kappa_2$. In each example, the set
of structuring elements $\mathcal B$ is constructed from a single
set $B$ in $\RR^2$ by rotating $B$ in all possible ways, that is,
$\mathcal B=\{RB\mid B\in\RR^2, R\in SO(3)\}$. For each example we
write
\[
SI_hu(\fxx)=\sup_{B\in \mathcal B}\inf_{\fxy\in\fxx+hB}u(\fxy)
\]
and
\[
IS_hu(\fxx)=\inf_{B\in \mathcal B}\sup_{\fxy\in\fxx+hB}u(\fxy),
\]
where $0<h\leq1$.

\begin{exercise}\label{exercise11.2}Let $B$ be a segment of length 2 centered at the
origin. Our aim is to show that
\begin{align*}
IS_hu&=u+h^2\frac1{2}\kappa_1^+(u)|Du|+O(h^3),\\
SI_hu&=u+h^2\frac1{2}\kappa_2^-(u)|Du|+O(h^3).
\end{align*}
This implies
\[
IS_hu +SI_hu=u+h^2\frac1{2}({\rm sign}(\kappa_1(u))+{\rm
sign}(\kappa_2(u)))\min(|\kappa_1(u)|,|\kappa_2(u)|)+O(h^3).
\]
\begin{enumerate}[(i)]\setlength\itemsep{2pt}
\item The first step is to compute $H(h\fxb)$.  One way to do this
is to write $x=r\sin\phi$, $y_2=r\cos\phi\cos\theta$,
$y_3=r\cos\phi\cos\theta$, and use an argument similar to that
given in the proof of Theorem \ref{theorem10.5} to show that, for
a fixed $\theta$ and small $h$, the ``inf-sup'' of
\[
r\sin\phi+hb_2r^2\cos^2\phi\cos^2\theta+hb_3r^2\cos^2\phi\cos^2\theta
\]
always occurs at $\phi=0$. Then $H(h\fxb)=hH(\fxb)$ and
\[
H(\fxb)=\inf_{B\in\mathcal B}\sup_{\fxy\in B}\,(b_2y_2^2+b_3y_3^2)=\inf_{\theta}\sup_{0\leq r\leq1}r^2(b_2\cos^2\theta+b_3\sin^2\theta).
\]
Deduce that $b_2<0$ or $b_3<0$ implies $H(\fxb)=0$ and that $0\leq
b_2\leq b_3$ implies $H(\fxb)=b_2$. \item Since $H(0)=0$, deduce
from Theorem \ref{theorem11.1} that
\begin{equation}\label{eq1213}
IS_hu(\fxx)=u(\fxx)+h^2\frac1{2}\kappa_1^+(u)(\fxx)|Du(\fxx)|+O(h^3).
\end{equation}
\end{enumerate}
\end{exercise}

\begin{exercise}\label{exercise11.3}Let $B$ be the union of two symmetric points $(1,0,0)$ and $(-1,0,0)$. Use the techniques of Exercise 11.2 to show that
\begin{align*}
IS_hu&=u+h^2\frac1{2}\min\{\kappa_1(u),\kappa_2(u)\}|Du|+O(h^3);\\
SI_hu&=u+h^2\frac1{2}\max\{\kappa_1(u),\kappa_2(u)\}|Du|+O(h^3);\\
IS_h+SI_hu&=u+h^2\frac1{2}(\kappa_1(u)+\kappa_2(u))|Du|+O(h^3).\\
\end{align*}
The last formula shows that the operator  $IS_h+SI_h$ involves the
mean curvature of $u$ at $\fxx$. \finex
\end{exercise}

\begin{exercise}\label{exercise11.4}Let $B$ consist of two orthogonal segments of length two centered at the origin.
\begin{enumerate}[(i)]\setlength\itemsep{2pt}
\item Show that
\begin{align*}
IS_hu&=u+h^2\frac1{2}\Big(\frac{\kappa_1(u)+\kappa_2(u)}{2}\Big)^+|Du|+O(h^3);\\
SI_hu&=u+h^2\frac1{2}\Big(\frac{\kappa_1(u)+\kappa_2(u)}{2}\Big)^-|Du|+O(h^3).
\end{align*}
\item Show that you can get the mean curvature by simply taking $B$ to be the four endpoints of the orthogonal segments. Check that another possibility for obtaining the mean curvature is to alternate these operators or to add them. \finex
\end{enumerate}
\end{exercise}

\section{Comments and references}
The references for this chapter are essentially the same as those
for  Chapter \ref{Chap10}. The main theorem on the asymptotic
behavior of morphological filters was first stated and proved in
\cite{Guichard:1995:PDE} and \cite{Guichard:1996:PDE}. The
examples developed in Exercises \ref{exercise11.2},
\ref{exercise11.3}, and \ref{exercise11.4} have not been published
elsewhere. The consistency of Gaussian smoothing followed by
thresholding and mean-curvature motion was proved in increasing
mathematical sophistication and generality by Merriman, Bence, and
Osher \cite{Merriman:1992:DGM}, Mascarenhas
\cite{Mascarenhas:1992:DGM}, Barles and Georgelin
\cite{Barles:1995:SPC}, and Evans \cite{Evans:1993:CAM}. Our
presentation is slightly more general than the ones cited because
we allow any nonnegative weight function in the Schwartz class.
The most general result was given by Ishii in
\cite{Ishii:1994:GBM}.
