\documentclass[10pt,twocolumn,letterpaper]{article}

\usepackage{cvpr}
\usepackage{times}
\usepackage{epsfig}
\usepackage{amsmath}
\usepackage{amssymb}

% Include other packages here, before hyperref.

% If you comment hyperref and then uncomment it, you should delete
% egpaper.aux before re-running latex.  (Or just hit 'q' on the first latex
% run, let it finish, and you should be clear).
\usepackage[pagebackref=true,breaklinks=true,letterpaper=true,colorlinks,bookmarks=false]{hyperref}


% \cvprfinalcopy % *** Uncomment this line for the final submission

\def\cvprPaperID{193} % *** Enter the CVPR Paper ID here
\def\httilde{\mbox{\tt\raisebox{-.5ex}{\symbol{126}}}}

% Pages are numbered in submission mode, and unnumbered in camera-ready
\ifcvprfinal\pagestyle{empty}\fi
\begin{document}


% Don't want date printed
\date{}
% Make title bold and 14 pt font (Latex default is non-bold, 16pt)
\title
{\Large\bf Level Set Segmentation using Statistical Shape Priors}

% For single author (just remove % characters)
%\author{I. M. Anonymous \\
%  My Department \\
%  My Institute \\
%  My City, STATE, zip}
% For two authors (default example)
%\author{\begin{tabular}[t]{c@{\extracolsep{3em}}c}
%A. El-Baz, A. A. Farag  & G. Gimel'farb\\
%\\
%Computer Vision and  & CITR, Tamaki Campus\\
%Image Processing Laboratory & Dept. of Computer Science\\
%University of Louisville & University of Auckland \\
%Louisville, KY 40292 & Auckland, New Zealand
%\end{tabular}}
\maketitle

%%%%%%%%% ABSTRACT
\begin{abstract}
A novel and robust 3-D segmentation approach is proposed using
level sets based on shape constraints. The approach depends on
both the gray level and shape information. A partial differential
equation (PDE) is developed to control the evolution of the level
sets. The PDE does not include weighting coefficients to be tuned,
overcoming the disadvantages of other PDE approaches. The shape
information is gathered from a set of the signed distance maps
representing the training data as a histogram of the occurrences
of the points inside and outside the object. We use a novel
statistical approach to get a probability density function (pdf)
for the signed distance map of the points inside and outside and
also the distribution of gray level inside and outside the object.
The proposed statistical approach is based on modelling the
empirical density function (normalized histogram of occurrence)
for either the gray level distribution or signed distance map with
a linear combination of Gaussians (LCG) with positive and negative
components. We modify an Expectation-Maximization (EM) algorithm
to deal with the LCGs and also propose a novel EM-based sequential
technique to get a close initial LCG approximation for the
modified EM algorithm to start with. The pdf's of the signed
distance and intensity gray level are embedded in the speed
function of the level set specifying the direction of evolution.
Experimental results show how the approach is accurate in
segmenting different types of 2-D and 3-D data sets including
medical images.
\end{abstract}

\section{Introduction}

Surgical planning, navigation, medical visualization and
diagnostics all benefit from image segmentation; and level sets
segmentation has retained its attention (see
\cite{A3,A4,A6,A15,A1,A2,Han03}) due to their topological
flexibility and independence of the parameterizations of the
evolving contour. However, segmentation process is still a
challenge because of the image noise and inhomogeneities;
therefore segmentation algorithms can not depend only on image
information but also have to exploit the prior knowledge of shapes
and other properties of the structures to be segmented.

The incorporation of shapes and deformable models became popular
with Leventon {\em et al.}\cite{Leventon2000,Leventon00} and Shen
{\em et al.}~\cite{Shen00} by attracting the level set function to
the likely shapes from a training set specified by principal
component analysis (PCA). Recently, the approach of
\cite{Leventon2000} was extended in \cite{Tsai03,Tsai04} and
\cite{Yang04}. In \cite{Tsai03}, shapes are represented with a
linear combination of 2D distance maps where the weight estimates
maximize the distance between the mean gray values inside and
outside the shape. In \cite{Tsai04}, the idea is to use a
deterministic model to represent the desired shape as a linear
combination of the weighted signed 2D distance maps, the weights
being estimated by minimizing a mutual information based cost
function. However, as noted in \cite{Tsai03}, the space of
distance maps is not closed with respect to linear operations, so
that the distance map for zero level shape of a linear combination
of the maps does not necessarily coincide with the latter
combination. As a result, the model may produce an unpredictable
shape.

The training shapes are represented with the like linear
combination of the training signed 2D distance maps also in the
papers of a second research group \cite{Yang04}. Here, all the
shapes are described with a multidimensional Gaussian probability
model of the weights in the linear combination. But apart from a
probabilistic description, the model has similar inconsistencies.

In \cite{A5}, registration is combined with the segmentation
process in an energy minimization problem. The evolving
curve/surface is embedded in a higher dimensional level set
function and registered iteratively with a shape model. In
\cite{Rouss_2004}, a 3D shape-segmentation approach is proposed
where a shape model is built from a set of training shapes using
distance functions. A level set function evolves minimizing the
shape alignment energy and the intensity gray level.

In this paper, a novel and robust segmentation approach is
proposed based on the level set technique with shape constraints.
The segmentation approach depends on both gray level intensity and
shape information; whereas the shape information is gathered from
a set of training shapes. A signed distance function is assigned
to each shape and the signed distance values are collected in the
form of a histogram representing the occurrences of each value.
Probability density functions of the object and background are
formed based on the signed distance value in addition to the
 gray level intensity. These functions are estimated using our new
approach known as the modified expectation maximization (EM) which
estimates the density function by a linear combination of
Gaussians (LCG) with positive and negative components. The
estimated pdf's are used in a variational approach making the
segmentation accurate and fast.

\section{Shape Modelling by Level Sets}\label{sec:Sh_ls}

Shape representation is the main task in shape analysis. The
selection of such representation is very important in several
computer vision and medical applications such as registration and
segmentation. There are several ways described in
\cite{Leventon2000,Seb2001,Siddi1999}. Although some of these ways
are powerful enough to capture local deformations, they require a
large number of parameters to deal with important shape
deformations. So an emerging way to represent shapes is derived
using level sets \cite{Osher_1988}. This representation is
invariant to translation and rotation. Given a curve/surface $V$
that represents boundaries of a certain shape, we can define the
following level set function:
\begin{equation}\label{e:eq_sh}
\phi(x,y,z)= \left \{
\begin{array}{lll}\
0 \mathrm{,}~~~~~~~~~~~~~~~~~~~~~~~~~~(x,y,z) \in V\\
d((x,y,z),V)~~~~~~~ (x,y,z) \in R_{V}\\
-d((x,y,z),V)~~~~~\mathrm{Otherwise}\\
\end{array}
\right.
\end{equation}
where $R_V$ is the region defined by the shape and $d((x,y,z),V)$
is the minimum Euclidean distance between the image location
$(x,y,z)$ and the curve/surface $V$.

Such representation can account for local deformations that are
not visible for iso-contours that are far away from the original
shape, and for geometrical features of the shape that can also be
derived naturally from this representation.

By this representation, we can construct a database of
curves/surfaces and signed distance functions that represent
variations for a certain shape. So, given a set of aligned
curves/surfaces $V_1,...,V_N$, level sets $\phi_1,...,\phi_N$ are
calculated as training data. From this information, we can extract
a histogram of the occurrences of signed distance values which
characterizes the shape and its local variations. Also, a mean
curve/surface $V_M$ is calculated as an average of the
corresponding points of all the training curves/surfaces. Simply,
we can get this curve/surface by picking each point in the first
training curve/surface and get the nearest points in the other
curves/surfaces, then calculate the average as follows:

\begin{equation}\label{e:eq_mean}
V_M(x_j,y_j,z_j)=\frac{1}{N}\sum_{i=1}^N~V_i{(x_i,y_i,z_i)},
\end{equation}
where $(x_i,y_i,z_i)$ are the corresponding points.

\section{Curve/Surface Evolution and Level Sets}\label{sec:ls}

Within the level set formalism \cite{A7,A20}, the evolving
curve/surface is a propagating front embedded as the zero level of
a scalar function $\phi(x,y,z,t)$. The continuous change of
$\phi(x,y,z,t)$ can be described by the partial differential
equation:
\begin{equation}\label{e:eq1}
\frac{\partial\phi(x,y,z,t)}{\partial{t}}+F(x,y,z)|\nabla\phi(x,y,z,t)|=0,
\end{equation}
where $F(x,y,z)$ is a velocity function and $\nabla=
[\frac{\partial}{\partial x},\frac{\partial}{\partial
y},\frac{\partial}{\partial z}]^{T} $. The function
$\phi(x,y,z,t)$ deforms iteratively according to $F(x,y,z)$, and
the position of the 2D/3D front is given at each iteration by
solving the equation $\phi(x,y,z,t)=0$. Practically, instead of
Eq.(\ref{e:eq1}), the value $\phi(x,y,z,t_{n+1})$ at step $n + 1$
is computed from $\phi(x,y,z,t_{n})$ at step $n$ by the relation:

\begin{equation}\label{e:eq2}
\phi(x,y,z,t_{n+1})=\phi(x,y,z,t_{n})-\triangle{t}\cdot
F|\nabla\phi(x,y,z,t_{n})|.
\end{equation}
The design of the velocity function $F(x,y,z)$ plays the major
role in the evolutionary process. Among several formulations
proposed in \cite{A8,A9}, we have chosen the following
formulation:
\begin{equation}\label{e:eq3}
F(x,y,z)=\nu-\epsilon{k(x,y,z)},
\end{equation}
where $\nu = 1$ or $-1$ for the contracting or expanding front
respectively, $\epsilon$ is a smoothing coefficient which is
always small with respect to $1$, and $k(x,y,z)$ is the local
curvature of the front. The latter parameter acts as a
regularization term.

\section{Statistical Model for Density Estimation}\label{sec:Model}

In this paper we introduce a new algorithm called a modified
Expectation-Maximization algorithm that approximates an empirical
probability density function of scalar data with a linear
combination of Gaussians (LCG) with positive and negative
components. Due to both positive and negative components, the LCG
approximates inter-class transitions more accurately than a
conventional mixture of only positive Gaussians.

This approach is suitable for estimating the marginal density for
either the gray level distribution $p_{g}(q)$ or signed distances
$p_{s}(d)$ in each region in the given image. In the following
section we will describe this model for estimating the marginal
density for the gray level distribution $p_{g}(q)$ in each region
and the similar way can be used to estimate the density of the
signed distances $p_{s}(d)$ in the given image.

To identify the model accurately, we approximate the marginal gray
level probability density in each region with a LCG having
$C_{\mathrm{p},i}$ positive and $C_{\mathrm{n},i}$ negative
components:% \cite{goshtasby:oneill:99,poggio:girosi:90}:

\begin{equation}
p_{g}(q|i) = \sum_{r = 1}^{C_{\mathrm{p},i}}
w_{\mathrm{p},i,r}\varphi(q|\theta_{\mathrm{p},i,r}) -\sum_{l =
1}^{C_{\mathrm{n},i}}
w_{\mathrm{n},i,l}\varphi(q|\theta_{\mathrm{n},i,l}); \label{eq:2}
\end{equation}
such that \(\int_{-\infty}^{\infty}p_{g}(q|i)dq = 1\). Here, $q$
is the gray level, and \(\varphi(q|\theta)\) is a Gaussian density
having a shorthand notation \(\theta = (\mu,\sigma^{2})\) for its
mean, $\mu$, and variance, $\sigma^{2}$. In contrast to more
conventional normal mixture models, the components are now both
positive and negative and have only one obvious restriction in
line with Eq.~(\ref{eq:2}): \( \sum_{r
=1}^{C_{\mathrm{p},i}}w_{\mathrm{p},i,r} -\sum_{l
=1}^{C_{\mathrm{n},i}}w_{\mathrm{n},i,l} = 1\). These weights are
not the prior probabilities, and the LCG of Eq.~(\ref{eq:2}) is
considered as a functional form of the approximation of a
probability density depending on parameters $(w,\theta)$ of each
component.

The mixture of $K$ LCGs,  $p(q)=\sum_{i=1}^{K}w_{i}p(q|i)$, has
the same form but a larger number of components, e.g.,
$C_{\mathrm{p}}=\sum_{i=1}^{K}C_{\mathrm{p},i}$ and
$C_{\mathrm{n}}=\sum_{i=1}^{K}C_{\mathrm{n},i}$ if all the values
$\theta_{\mathrm{p},i,r}$ and $\theta_{\mathrm{n},i,l}$ differ for
the individual models:

\begin{equation}
p_{g}(q) = \sum_{r = 1}^{C_{\mathrm{p}}}w_{\mathrm{p},r}
\varphi(q|\theta_{\mathrm{p},r}) - \sum_{l =
1}^{C_{\mathrm{n}}}w_{\mathrm{n},l}
\varphi(q|\theta_{\mathrm{n},l}) \label{eq:4}
\end{equation}
To identify this model in the unsupervised mode, the mixed
empirical distribution of gray levels over the image has to be
first represented by a joint LCG of Eq.~(\ref{eq:4}) and then
partitioned into individual LCG-models for each class $i =
1,\ldots,K$.

Under the fixed number of the positive and negative components,
$C$, the model parameters $\textbf{w} = \{w_{c}; c = 1,\ldots,C\}$
and $\mathbf{\Theta} = \{\theta_{c}: c = 1,\ldots,C \}$ maximizing
the image likelihood can be found using an EM algorithm introduced
in Section~\ref{s:par-alg}. It modifies the conventional EM-scheme
to take account of the components with alternating signs.

The modified EM algorithm is sensitive to both its initial state
specified by the numbers of positive and negative Gaussians, and
the initial parameters (mean and variance) of each component. To
find a close initial LCG-approximation of the empirical
distribution, we develop in Section~\ref{sec:sequential} a
sequential initializing EM-based algorithm.

\subsection{Modified EM Algorithm for LCGs} \label{s:par-alg}

Let $f(q), q\in\mathbf{Q}$ be an empirical relative frequency
distribution representing an unknown probability density function
\(\psi(q)\) such that \( \int_{-\infty}^{\infty}\psi(q)dq \equiv
\sum_{q=0}^{Q}f(q) = 1 \). We assume that \(f(q)\) is approximated
by an LCG $P_{g:\mathbf{w},\mathbf{\Theta}}$ with
\(C_{\mathrm{p}}\) positive and \(C_{\mathrm{n}}\) negative
components \(\varphi(q|\theta)\):
\begin{equation}
p_{g:\mathbf{w},\mathbf{\Theta}}(q) =
\sum\limits_{r=1}^{C_{\mathrm{p}}}
w_{\mathrm{p},r}\varphi(q|\theta_{\mathrm{p},r}) -
\sum\limits_{l=1}^{C_{\mathrm{n}}}
w_{\mathrm{n},l}\varphi(q|\theta_{\mathrm{n},l}) \label{e:lcg}
\end{equation}
In line with Eq.~(\ref{e:lcg}), the positive weights
\(\mathbf{w}\) are restricted as follows:
\begin{equation}
\sum\limits_{r=1}^{C_{\mathrm{p}}}w_{\mathrm{p},r} -
\sum\limits_{l=1}^{C_{\mathrm{n}}}w_{\mathrm{n},l} = 1
\label{e:condition}
\end{equation}
We also assume here that the numbers \(C_{\mathrm{p}}\) and
\(C_{\mathrm{n}}\) of the components of each type are known after
the initialization in Section~\ref{sec:sequential} and do not
change during the EM process. The initialization provides also the
starting parameter values \(\mathbf{w}^{[0]}\) and
\(\mathbf{\Theta}^{[0]}\).

The probability densities form a proper subset of the set of the
LCGs due to the additional restriction
\(p_{\mathbf{w},\mathbf{\Theta}}(q) \ge 0\), which holds
automatically only for probability mixtures with no negative
components.

The LCG that provides a local maximum of the log-likelihood of the
empirical data:
\begin{equation}
L(\mathbf{w},\mathbf{\Theta})= \sum\limits_{q\in\mathbf{Q}}f(q)
\log p_{g:\mathbf{w},\mathbf{\Theta}}(q) \label{e:likelhd}
\end{equation}
can be found using the iterative block relaxation process
extending conventional EM schemes.

Let \( p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q) =
\sum_{r=1}^{C_{\mathrm{p}}}
w_{\mathrm{p},r}^{[m]}\varphi(q|\theta_{\mathrm{p},r}^{[m]}) -
\sum_{l=1}^{C_{\mathrm{n}}}
w_{\mathrm{n},l}^{[m]}\varphi(q|\theta_{\mathrm{n},l}^{[m]}) \) be
the LCG at step, or iteration $m$. Relative contributions of each
data item \(q=0,\ldots,Q\) into each positive and negative
Gaussian at the step $m$ are specified by the following respective
conditional weights
\begin{eqnarray}
\begin{array}{l}
\pi_{\mathrm{p}}^{[m]}(r|q) =
\frac{w_{\mathrm{p},r}^{[m]}\varphi(q|\theta_{\mathrm{p},r}^{[m]})}
{p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q)}; \; \; \;
\pi_{\mathrm{n}}(l|q) =
\frac{w_{\mathrm{n},l}^{[m]}\varphi(q|\theta_{\mathrm{n},l}^{[m]})}
{p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q)}
\label{e:posteriors}
\\
\sum\limits_{r=1}^{C_{\mathrm{p}}}\pi_{\mathrm{p}}^{[m]}(r|q) -
\sum\limits_{l=1}^{C_{\mathrm{n}}}\pi_{\mathrm{n}}^{[m]}(l|q) = 1;
\; q = 0,\ldots,Q \label{e:post:cond}
\end{array}
\end{eqnarray}
Using these weights, the log-likelihood of Eq.~(\ref{e:likelhd})
can be rewritten in the equivalent form:
\begin{eqnarray}
%\begin{array}{l}
L({\mathbf{w}}^{[m]},{\mathbf{\Theta}}^{[m]}) &= %\\
\sum\limits_{q=0}^{Q}f(q)[
\sum\limits_{r=1}^{C_{\mathrm{p}}}\pi_{\mathrm{p}}^{[m]}(r|q) \log
p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q)]\nonumber \\ &-
\sum\limits_{q=0}^{Q}f(q)[
\sum\limits_{l=1}^{C_{\mathrm{n}}}\pi_{\mathrm{n}}^{[m]}(l|q) \log
p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q)]
%\end{array}
 \label{e:lklh:equiv}
\end{eqnarray}
where \(\log p_{g:{\mathbf{w}},{\mathbf{\Theta}}}^{[m]}(q)\) in
the first and the second brackets should be replaced with the
equal terms: \(\log w_{\mathrm{p},r}^{[m]} + \log
\varphi(q|\theta_{\mathrm{p},r}^{[m]}) - \log
\pi_{\mathrm{p}}^{[m]}(r|q) \) and \( \log w_{\mathrm{n},l}^{[m]}
+ \log \varphi(q|\theta_{\mathrm{n},l}^{[m]}) - \log
\pi_{\mathrm{n}}^{[m]}(l|q) \), respectively.

The block relaxation converging to a local maximum of the
likelihood function in Eq.~(\ref{e:lklh:equiv})  repeats
iteratively the following two steps:
\begin{enumerate}
\item E-step \([m +1]\): to find the parameters \(\mathbf{w}^{[m
+1]}\), \(\mathbf{\Theta}^{[m +1]}\) by maximizing
\(L(\mathbf{w},\mathbf{\Theta})\) under the fixed conditional
weights of Eq.~(\ref{e:posteriors}) for the step $m$, and \item
M-step \([m +1]\): to find these latter weights by maximizing
\(L(\mathbf{w},\mathbf{\Theta})\) under the fixed parameters
\(\mathbf{w}^{[m +1]}\), \(\mathbf{\Theta}^{[m +1]}\)
\end{enumerate}
until the changes of the log-likelihood and all the model
parameters become small.

The E-step performs the conditional Lagrange maximization of the
log-likelihood of Eq.~(\ref{e:lklh:equiv}) under the restriction
of Eq.~(\ref{e:condition}) to obtain the following estimates of
the weights:
\[
w_{\mathrm{p},r}^{[m+1]} =
\sum\limits_{q\in\mathbf{Q}}f(q)\pi_{\mathrm{p}}^{[m]}(r|q);\;
w_{\mathrm{n},l}^{[m+1]}
=\sum\limits_{q\in\mathbf{Q}}f(q)\pi_{\mathrm{n}}^{[m]}(l|q)
\]
Then the parameters of each Gaussian are obtained by the
unconditional maximization just as in the conventional EM scheme
(below ``$c$'' stands for ``$p$'' or ``$n$'', respectively):
\[
\begin{array}{l}
\mu_{\mathrm{c},r}^{[m+1]} = \frac{1} {w_{\mathrm{c}, r}^{[m+1]}}
\sum\limits_{q\in\mathbf{Q}} q \cdot f(q) \pi_{\mathrm{c}}^{[m]}(r|q)\\
(\sigma_{\mathrm{c},r}^{[m+1]})^{2} = \frac{1}{w_{\mathrm{c},
r}^{[m+1]}} \sum\limits_{q\in\mathbf{Q}} \left (q -
\mu_{\mathrm{c}, i}^{[m+1]}\right )^{2} \cdot
f(q)\pi_{\mathrm{c}}^{[m]}(r|q)
\end{array}
%\label{e:e-step:gauss}
\]

The M-step performs the conditional Lagrange maximization of the
log-likelihood of Eq.~(\ref{e:lklh:equiv}) under the \(Q+1\)
restrictions of Eq.~(\ref{e:post:cond}), and determines the
conditional weights \(\pi_{\mathrm{p}}^{[m+1]}(r|q)\) and
\(\pi_{\mathrm{n}}^{[m+1]}(l|q)\) of Eq.~(\ref{e:posteriors}) for
all \(r=1,\ldots,C_{\mathrm{p}}\); \(l=1,\ldots,C_{\mathrm{n}}\)
and \(q=0,\ldots,Q\). The modified EM-algorithm is valid until
these weights are strictly positive, and the initial
LCG-approximation should comply to this limitation. The iterations
have to be terminated when the log-likelihood of
Eq.~(\ref{e:lklh:equiv}) begins to decrease.
%Generally, if the initialization is incorrect, this algorithm may
%diverge at the first iteration. Thus the initial LCG has to
%closely approximate the empirical distribution.

\subsection{Sequential EM-Based Initialization}\label{sec:sequential}

We assume that the number of dominant modes is equal to the given
number of classes. To simplify the notation, let the empirical
distribution have only two separate dominant modes representing
the object and the background, respectively. The algorithm we
present below is easily extended to the general case of \(K>2\)
dominant modes. We assume that each dominant mode is roughly
approximated with a single Gaussian and the deviations of the
empirical density from the two-component dominant Gaussian mixture
are described by other components of the LCG in Eq.~(\ref{eq:4}).
Therefore the model has the two dominant positive weights, say,
$w_{\mathrm{p},1}$ and $w_{\mathrm{p},2}$ such that
$w_{\mathrm{p},1} + w_{\mathrm{p},2} = 1$, and a number of
``subordinate" weights of smaller absolute values such that
$\sum_{r=1}^{C_{\mathrm{p}}}w_{\mathrm{p},r} -
\sum_{l=1}^{C_{\mathrm{n}}}w_{\mathrm{n},l} = 0$.

The following sequential algorithm allows for estimating both the
weights and parameters of the individual Gaussians in the latter
LCG model, including the number of the non-dominant components.
\begin{enumerate}
\item Approximate a given empirical distribution \(f(q)\), of gray
levels in the given image, with a dominant mixture $p_{2}(q)$, of
two Gaussians using the conventional EM-algorithm.

\item Find the deviations
$\mathbf{\Delta}=[\Delta(q)=f(q)-p_{2}(q):\ q\in\mathbf{Q}]$
between $f(q)$ and $p_{2}(q)$ and split them into the positive and
negative parts such that \(\delta(q) = \delta_{\mathrm{p}}(q) -
\delta_{\mathrm{n}}(q)\):

\begin{equation}
\begin{array}{lll}
\mathbf{\Delta}_{\mathrm{p}} & = & [ \delta_{\mathrm{p}}(q) =
\max\{
\delta(q),0\}: q\in \mathbf{Q} \}\\
\mathbf{\Delta}_{n} & = & [ \delta_{n}(q) = \max\{-\delta(q),0\}:
q\in \mathbf{Q} \}
\end{array}
\label{eq:8}
\end{equation}

\item Compute the scaling factor for the deviations: \( scale =
\int_{-\infty}^{\infty}\delta_{\mathrm{p}}(q)dq \equiv
\int_{-\infty}^{\infty}\delta_{\mathrm{n}}(q)dq \).

\item If the factor $s$ is less than a given accuracy threshold,
terminate and return the model $p_{g}(q) = p_{2}(q)$.

\item Otherwise consider the scaled-up absolute deviations
$\frac{1}{scale}\mathbf{\Delta}_{\mathrm{p}}$ and
$\frac{1}{scale}\mathbf{\Delta}_{\mathrm{n}}$ as two new
``empirical densities'' and iteratively the conventional
EM-algorithm to find sizes $C_{\mathrm{p}}$ and $C_{\mathrm{n}}$
of the Gaussian mixtures, $p_{p}(q)$ and $p_{n}(q)$, respectively,
approximating the scaled-up deviations. The size of each mixture
corresponds to the minimum of the integral absolute error between
the scaled-up absolute deviation $\mathbf{\Delta}_{\mathrm{p}}\)
(or \(\mathbf{\Delta}_{\mathrm{n}}$) and its model $p_{p}(q)$ (or
$p_{{n}}(q)$). The number of the components is increasing
sequentially by unit step while the error is decreasing. \item
Scale down the subordinate  models $p_{p}(q)$ and $p_{n}(q)$ (i.e.
scale down the weights of their components) and add the scaled
model $p_{p}(q)$ to and subtract the scaled model $p_{n}(q)$ from
the dominant model $p_{2}(q)$ in order to form the desired model
$p_{g}(q)$ of the size \(C=2+C_{\mathrm{p}}+C_{\mathrm{n}}\).
\end{enumerate}
Since the EM algorithm converges to a local maximum of the
likelihood function, it may be repeated several times with
different initial parameter values for choosing the model giving
the best approximation. In principle, this process can be repeated
iteratively in order to approximate more and more closely the
residual absolute deviations between $f(q)$ and \(p_{g}(q)\). But
because each Gaussian in the latter model impacts all the values
\(p(q)\), the iterations should be terminated when the
approximation quality begins to decrease.

The final mixed LCG-model $p_{g}(q)$ has to be split into $K$
LCG-submodels, one per class, by associating each subordinate
component with a particular dominant term in such a way as to
minimize the expected misclassification rate. To illustrate the
association principle, let us consider the bi-modal case with the
two dominant Gaussians having the mean values $\mu_{1}$ and
$\mu_{2}$; $ 0 < \mu_{1} <\mu_{2}< Q$. Let all the subordinate
components be ordered by their mean values, too. Then let those
with the mean values smaller than $\mu_{1}$ and greater than
$\mu_{2}$ relate to the first and second class, respectively. The
components having the mean values in the range $[\mu_{1},
\mu_{2}]$ are associated with the classes by simple thresholding
such that the means below the threshold, \(t\),  belong to the
components associated with the first class. The desired threshold
minimizes the classification error $e(t)$:

\begin{equation}
e(t) = \int\limits_{-\infty}^{t}p_{g}(q|2)dq +
\int\limits_{t}^{\infty}p_{g}(q|1)dq. \label{eq:15}
\end{equation}

\section{Evolutionary Curve/Surface Model}\label{sec:em}

The term $(\nu=\pm~1)$ in Eq.(~\ref{e:eq3}) specifies the
direction of the front propagation. Several approaches were
developed to make all fronts either contracting or expanding (see,
e.g., \cite{A11}) in order to evolve in both directions and avoid
overlaps between the regions.The problem can be reformulated as
classification of each point at the evolving front. If the point
belongs to the associated class (object), the front expands
otherwise it contracts.

\subsection{PDE System}
The classification decision is based on Bayesian
decision~\cite{A24} at voxel $(x,y,z)$ at the front. The term
$(\nu)$ for each point is replaced by the function $\nu(x,y,z)$ so
the velocity function is defined as:
\begin{equation}\label{e:eq11}
F(x,y,z)=\nu(x,y,z)-\epsilon\cdot k(x,y,z).
\end{equation}
where
\begin{equation}\label{e:eq110}
\nu(x,y,z)= \left \{
\begin{array}{ll}\
-1 & \mathrm{if} ~~ p_{g}(q|1)*p_{s}(d|1) > p_{g}(q|2)*p_{s}(d|2)\\
~~~~1     & \mathrm{otherwise}
\end{array}
\right.
\end{equation}

If the voxel $(x,y,z)$ belongs to the object, the front will
expand, otherwise it will contract. Now, we put the
Eq.(~\ref{e:eq1}) in the general form using the derivative of the
Heaviside step function ($\delta_{\alpha}(.)$)(\cite{A18}) as
follows :
\begin{eqnarray}\label{e:eq12}
\frac{\partial\phi(x,y,z,t)}{\partial{t}} &=& \big(\epsilon \cdot
k(x,y,z)-\nu(x,y,z)\big)
\nonumber\\
&\times&\delta_{\alpha(x,y,z)}\big(\phi(x,y,z)\big)\big|\nabla\phi(x,y,z)\big|~~~.
\end{eqnarray}
The function $\delta_{\alpha}(.)$ selects the narrow band points
around the front and the parameter $\alpha$ controls the width of
the narrow band.

%\subsection{Data Consistency Coefficient $h$}
%Data consistency term is always a function of image gradient. But
%this is a disadvantage when we have images with high noise. Our
%approach is similar to the one described in \cite{Baillard:00}. We
%have related the stopping factor $h$ to the posterior probability
%of having a transition between the $Object$ and the $Background$.
%Let $(x)$ be a voxel of the current interface, and $\lambda$ be
%the estimated class $Object$ or $Background$ in this case. The
%posterior probability of $(x)$ being a transition, given
%$\lambda$, is given by:
%\begin{equation}\label{e:hI}
%p_T (x|\lambda)= \left \{
%\begin{array}{lll}\
%p_{s}(d|2)*p_{g}(q|2) & \mathrm{if} & \lambda \in Object\\
%p_{s}(d|1)*p_{g}(q|1) &&\mathrm{otherwise}
%\end{array}
%\right .
%\end{equation}
%
%Where in this case the probability density functions are
%calculated at the point $(x)'$ which is a neighboring voxel of
%$(x)$ located outside the area defined by hypersurface $\phi$.
%Thus if $(x)$ is more likely to be inside the object to be
%segmented, then the posterior transition probability is the
%probability of $(x)'$ to be located outside the object to be
%segmented.
%
%The data consistency term $h(p_T)$ is defined as a decreasing
%function of $p_T (x|\lambda)$. Then $h$ is defined to be as
%follows:
%\begin{equation}\label{e:hI_final}
%h(p_T)=1-H_{\zeta}(p_T),
%\end{equation}
%where $\zeta$ is the prior probability of the object.
%\vspace{-3mm}
\subsection{Registration Step}

To make the proposed approach of shape-based segmentation
invariant for the scaling, rotation, and translation of the
object, the first step of the proposed approach is to align the
image with any image from our aligned database. The deformations
that we use are defined using the Free Form Deformations (FFD) as
explained in ~\cite{Rueckert:99}. The essence of FFD is to deform
an object by manipulating a regular control lattice overlaid on
its volumetric embedding space. One of the main advantages of the
FFD technique is that it imposes implicit smoothness constraints
during deformation, since it guarantees continuity at control
points and continuity elsewhere. Therefore there is no need to
introduce computationally expensive regularization components on
the deformed shapes.

\subsection{Algorithm}

\begin{enumerate}
    \item Estimate the pdf's of the object and background for the intensity
gray level and shape signed distance values using the modified
    EM.
    \item Register the image with any one of the aligned images from the database.
    \item Initialize the level set function.
    \item Mark the points of the narrow band.
    \item For each point calculate the signed distance value ($d$) as the
minimum Euclidean distance from
    $V_M$.
    \item If $p_{s}(d|1)*p_{g}(q|1)>p_{s}(d|2)*p_{g}(q|2)$, then the front
expands at this point, otherwise it
    contracts.
    \item Go to step 4. This process is repeated until the change in the level
set function is not significant.
\end{enumerate}
    %1. Estimate the pdf's of the object and background for the intensity
%gray level and shape signed distance values using the modified
%    EM.\\
%    2. Register the image with any one of the aligned images from the database.\\
%    3. Initialize the level set function.\\
%    4. Mark the points of the narrow band.\\
%    5. For each point calculate the signed distance value ($d$) as the
%minimum Euclidean distance from
%    $V_M$.\\
%    6. If $p_{s}(d|1)*p_{g}(q|1)>p_{s}(d|2)*p_{g}(q|2)$, then the front
%expands at this point, otherwise it
%    contracts.\\
%    7. Go to step 4. This process is repeated until the change in the level
%set function is not significant.

\section{Experimental Results}\label{sec:Model}

We illustrate the performance of the proposed techniques by
applying it on different 2-D and 3-D data sets. The first example
which we show is the segmentation of starfish. The segmentation
separates starfish from the surrounding background so that each
image has only two dominant objects ($K=2$): the darker background
and the brighter starfish. Figure~\ref{fig:star-1}(a) demonstrates
one of the aligned starfish images from our database.
Figure~\ref{fig:star-1}(b) shows one of the starfish images that
we need to segment. Figure~\ref{fig:star-1}(c) shows the result of
the registering the image shown in (a) to image in (b) using the
algorithm proposed in ~\cite{Rueckert:99}.

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=1.0in]{Image_1} &
\includegraphics[width=1.0in]{Image_2} &
\includegraphics[width=1.0in]{Image_3} \\
(a) & (b)& (c) \\
\end{tabular}
\caption{(a) One of the aligned starfish images from the database,
(b) Non-aligned image of starfish, (c) Results of the registration
of (a) and (b)} \label{fig:star-1}
\end{center}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cc}
  \includegraphics[width=1.5in]{2_Dominant})&
  \includegraphics[width=1.5in]{Deviation}\\
  (a) & (b) \\
   \includegraphics[width=1.5in]{Number_Comp}&
  \includegraphics[width=1.5in]{Error_Estimation}\\
  (c) & (d)
\end{tabular}
\end{center}
\caption{(a) Empirical density $f(q)$ approximated with the two
dominant Gaussian components $p_{g:2}(q)$ for the image shown in
Fig.~\ref{fig:star-1}(b), (b) Deviation and absolute deviation
between $f(q)$ and $p_{g:2}(q)$, (c) Estimating the number of
subordinate mixture, (d) Estimated density for the absolute
deviation} \label{star-initial}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=1.5in]{Final_Estimation}&
\includegraphics[width=1.5in]{EM2modifGLGlk_b}\\
(a) & (b) \\
\includegraphics[width=1.5in]{All_Comp}&
\includegraphics[width=1.5in]{Marginal_classes}\\
(c) & (d)
\end{tabular}
\end{center}
\caption{(a) Final density estimation of the bi-modal
distribution, (b) The dynamic changes of the log-likelihood of the
modified EM algorithm, (c) All components of the final LCG, (d)
The marginal density estimation for each class.}
\label{star:final}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=0.6in]{S_Image_0001.eps} &
\includegraphics[width=1.2in]{Final_Estimation_Star_Shape}&
\includegraphics[width=1.2in]{Star_Sea_Shape_Marginal}\\
(a) & (b) & (c)\\
\end{tabular}
\end{center}
\caption{(a) Average signed distance of aligned star images, (b)
Final estimated density, (c) The marginal density estimation for
each class.} \label{shape-star}
\end{figure}

Figure~\ref{star-initial} illustrates the sequential EM-based
initialization, (a) the empirical density of starfish image shown
in Fig.~\ref{fig:star-1}(b) and the initial mixture of two
Gaussians approximating the dominant modes, (b) the deviation
between the empirical density $f(q)$ and the mixture of the two
dominant components $p_{g:2}(q)$, (d) the estimated density of the
scaled deviation using the six Gaussian components which give the
minimum error between the estimated density for the deviation and
the empirical deviation as shown in Fig.~\ref{star-initial}(c).

Figure~\ref{star:final} presents the final LCG-model and its 8
components obtained by the modified algorithm as well as the
successive changes of the log-likelihood. The first 10 iterations
of the refining EM-algorithm increase the log-likelihood of
Eq.~(\ref{e:lklh:equiv}) from $-4$ to $-2.9$; then the modified EM
algorithm terminates since the log-likelihood begins to decrease.
The minimum classification error of 0.00018 between the starfish
and the background for the final LCG-model is obtained with the
threshold $t = 82$ in Eq.~(\ref{eq:15}). In this case the
LCG-components 1--3 and 4--6 correspond to the starfish and the
background, respectively.

Figure~\ref{shape-star}(a) shows the average signed distance of
all aligned segmented starfish. We used the approach as discussed
in the previous sections to estimate the marginal density that
describes the distribution of signed distance inside and outside
the starfish object. The results of this approach are shown in
Fig.~\ref{shape-star}(b) and (c).

%\vspace{-20mm}
\begin{figure}[htb!]
\centering
\begin{tabular}{ccc}
\includegraphics[width=1.0in]{Image_Initial.eps}&
\includegraphics[width=1.0in]{Image_Final.eps}&
\includegraphics[width=1.0in]{Image_Final_T.eps}\\
(a) & (b) & (c)
\end{tabular}
\caption{(a) The initialization of the level sets function to
segment the aligned starfish shown in Fig.~\ref{fig:star-1}(c),
(b) The segmentation of the aligned image, (c) The final
segmentation of starfish after multiplying the aligned starfish
image by the inverse of the transformation function (which we used
in the rigid registration) with error 0.4\% with the ground
truth.} \label{f:star_seg}
\end{figure}

Figure~\ref{f:star_seg} shows the result of our proposed
segmentation approach with error of 0.4\% with the ground truth.
%Also, more results on different real star-fish images are given in
%Fig. \ref{fig:Star-more}.

%\begin{figure}[ht!]
%\centering
%\begin{tabular}{|c|c|c|c|}
%\hline
%\includegraphics[width=1.25in,height=1.25in]{Star_More_001}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_002}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_003}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_004}\\ \hline
%Error 0.1034\% & Error 0.3082\% & Error 0.082\% & Error 0.0363\%\\
%\hline
%\includegraphics[width=1.25in,height=1.25in]{Star_More_005}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_006}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_007}&
%\includegraphics[width=1.25in,height=1.25in]{Star_More_008}\\ \hline
%Error 0.7722\% & Error 0.8375\% & Error 0.9060\% & Error
%1.2551\%\\ \hline
%\end{tabular}
%\caption{More results for starfish with different deformations and
%appearance models} \label{fig:Star-more}
%\end{figure}


In the following part, we test the proposed segmentation approach
on 3D data sets to show that it can work robustly for 3D images as
well. To get an accurate 3D shape model of the lateral ventricles
of the brain, the 3D images were taken from 20 subjects. However,
20 data sets are not enough to get an accurate shape for the
ventricles because the ventricles vibrate during the MRI or CT
scans. Therefore, to cover all the shape variations of the brain
ventricles for each subject, we performed finite element analysis
on the motion of the real brain ventricles.

For finite element analysis of the ventricle shape changes, we
assume the cerebrospinal fluid (CSF) inside lateral ventricles is
isotropic and linear elastic. So, the linear elastic mechanical
model is employed for finite element analysis. The Young's modulus
of CSF is 1000 Pascals and the Poisson ratio is 0.499
(\cite{Hung}). For adults, the pressure of CSF under normal
conditions can range between $50$ and $180$ $mmH_2O$. The median
pressure $115$ $mmH_2O$ of the CSF should be the one most people
have statistically. Therefore, we apply the uniform pressure $115$
$mmH_2O$ over the surfaces of the ventricles for each subject, and
perform finite element analysis to capture the variation of the
ventricles. After finite element analysis, the 3D structure is
re-sliced and 10 states of the ventricles are obtained resulting
in $20$ $subjects \times 10$  $states = 200$ $datasets$ for the
ventricles. The four states of one subject's ventricles are shown
in Fig.~\ref{ventreal:final}. Using the resulting 200 datasets for
the ventricles, we are following our density estimation approach
to estimate the density of the gray level distribution and the
signed distance map inside and outside the ventricles. In this
estimation, we are considering that the brain MRI consists of two
classes: one class is composed of the gray matter, white matter,
fat and bones; and the other class is the CSF of the brain (inside
and outside the ventricles). The results of density estimation
using the proposed approach are shown in Fig.~\ref{Real-ven-d} and
the segmentation results at different signal to ratios (SNR)
(obtained by adding Gaussian noise with different variance) are
shown in Fig.~\ref{vent:final3}.

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cccc}
\includegraphics[width=0.85in]{ventreal0}&
\includegraphics[width=0.75in]{ventreal7}&
\includegraphics[width=0.65in]{ventreal9}&
\includegraphics[width=0.55in]{ventreal11}\\
(a) & (b) & (c) & (d)
\end{tabular}
\end{center}
\caption{The four states of the real ventricles at t = 0 sec, 0.7
sec, 0.9 sec, 1.1 sec (final) from left to right and top to
bottom.} \label{ventreal:final}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=1.5in]{Mixed_Empr_sim_vent_Gray}&
\includegraphics[width=1.5in]{Gray_claases}\\
(a) & (b) \\
\includegraphics[width=1.5in]{1}&
\includegraphics[width=1.5in]{2}\\
(c) & (d)
\end{tabular}
\end{center}
\caption{(a) Final density estimation for the mixed frequency gray
level distribution, (b) The class model of gray level for each
class, (c) Final density estimation for the sign map inside and
outside the 3-D brain ventricle, (d) The final class model of
signed distance map for each class.} \label{Real-ven-d}
\end{figure}

\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cccc}
\includegraphics[width=0.75in]{Real_V_N_S_0}&
\includegraphics[width=0.75in]{Real_V_N_S_2}&
\includegraphics[width=0.75in]{Real_V_N_S_2_2}&
\includegraphics[width=0.75in]{Real_V_N_S_3}\\
(a) & (b) & (c) & (d)
\end{tabular}
\end{center}
\caption{ (a) Results of our segmentation approach at SNR = 28dB,
$error=0.23\%$ (b) SNR = 7.8dB, $error=0.9\%$ (c) SNR = -1.9dB,
$error=4.8\%$. (d) Result of the segmentation errors using the
conventional EM algorithm as a density estimator at SNR = 28dB,
$error=18.8\%$. The errors are shown in red.}\label{vent:final3}
\end{figure}

The hand segmentation of the radiologists may not always be
accurate because of hand shaking. Therefore, to get accurate
evaluation of our proposed approach, we are using a ventricle
phantom that resembles the geometrical structure of the real
ventricles. To use the mould ventricle phantom, it is scanned
using cone-beam CT machine so that the scans can be used for
finite element simulation. Following the same procedure of the
ventricle motion estimation as described above, we captured all
the variations of the motion of the phantom ventricles. Figure
\ref{vent:final-simulated} shows the four states of this ventricle
phantom. The shape changes of the ventricles for the subject and
the ventricle phantom are recorded in short videos that will be
supplemented.
\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cccc}
\includegraphics[width=1.1in]{vent0}&(a)&
\includegraphics[width=1.0in]{vent4}&(b)\\
%(a) & (b) \\
\includegraphics[width=0.9in]{vent7}&(c)&
\includegraphics[width=0.8in]{vent10}&(d)\\
%(c) & (d)\\
\end{tabular}
\end{center}
\caption{The four states of the phantom ventricles at t=0 sec, 0.4
sec, 0.7 sec, 1 sec (final) from left to right and top to bottom.}
\label{vent:final-simulated}
\end{figure}


\begin{figure}[htb!]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=1.50in]{Mixed_Empr_sim_vent}&
\includegraphics[width=1.50in]{Marginal_Empr_sim_vent}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{(a) Final density estimation for the sign map inside and
outside the 3-D ventricle phantom, (b) The final class model of
signed distance map for each class.} \label{Real-ven-phantom}
\end{figure}

For the database that was obtained from the geometrical phantom,
the inverse mapping method was used to get the same gray level
distribution as of the real ventricles; the gray level
distribution of which was shown in Fig.~\ref{Real-ven-d}(b). The
final step of our algorithm is to estimate the pdf that describes
the signed distance map inside and outside the geometrical
phantom. The results of our modified EM algorithm are shown in
Fig.~\ref{Real-ven-phantom}. Figure
\ref{vent:final-phantom-segment} shows the results of our
segmentation for the geometrical phantom at different signal to
noise ratios and the errors are calculated with respect to the
ground truth from the phantom.

\begin{figure}[ht!]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=1.0in]{S_N_1}&
\includegraphics[width=1.0in]{S_N_3}&
\includegraphics[width=1.0in]{S_N_4}\\
(a) & (b) & (c)
\end{tabular}
\end{center}
\caption{Segmentation using our approach at different SNR (a)
$SNR=28dB, error=0.01\%$, (b) $SNR=7.8dB, error=0.8\%$ (c)
$SNR=-1.9dB, error=3.68\%$. } \label{vent:final-phantom-segment}
\end{figure}

\section{Conclusion} We have presented a segmentation approach that depends on both the
intensity gray level and the shape information. Our modified EM is
used to estimate the density distribution of the intensity and
signed distance values. The density distributions are embedded in
the PDE that controls the evolution of the level set function. We
consider the registration between the average shape (2D/3D) and
the object to be segmented as a basic step in our approach. Unlike
the other approaches, our segmentation does not need energy
minimization avoiding the local minimum problem. Different types
of images are used and the results are promising. This technique
is very suitable to segment the anatomical structures that have
noise and inhomogeneity problems.

This algorithm is very robust and accurate, it is invariant to
translation, rotation and scaling. Therefore, the presented
segmentation algorithm is not only useful for medical imaging
society but also for the computer vision applications. Our future
work will include the segmentation of 2D and 3D color objects.


\begin{thebibliography}{99}
\small  % Use 9 point text.
% Use 9 point text.
\bibitem{A3}
V. Caselles, R. Kimmel, and G. Sapiro, ``Geodesic active
contours,'' {\em IJCV}, 22:61--79, 1997.
%
\bibitem{A4}
H-K. Zaho, T. Chan, B. Merriman, and S. Osher, ``A variational
level set approach to multiphase motion,'' {\em J. of
Computational Physics}, 127:179--195, 1996.
%
\bibitem{A6}
R. Goldenberg, R. Kimmel, E. Rivlin, and M. Rudzsky, ``Cortex
Segmentation: A Fast Variational Geometric Approach,'' {\em IEEE
Trans. on Medical Imaging}, Vol. 21, No. 2, pp. 1544--1551,
December 2002.
%
\bibitem{A15}
T. Chan and L. Vese, ``A multiphase level set framework for image
segmentation using the mumford and shah model,'' {\em
International Journal of Computer Vision}, 50(3):271�-293,2002.
%
\bibitem{A1}
S. Osher and N. Paragios, {\em Geometric Level Set Methods in
Imaging, Vision, and Graphics,} Springer, 2003.
%
\bibitem{A2}
D. Terzopoulos, ``Regularization of inverse visual problems
involving discontinuities,'' {\em IEEE Trans. on PAMI},
8(2):413--424, 1986.
%
\bibitem{Han03}
X. Han, C. Xu, and J. L. Prince, ``A Topology Preserving Level Set
Method for Geometric Deformable Models," IEEE Transactions on
PAMI, Vol. 25, No. 6, pp.755-768, June 2003.
%
\bibitem{Leventon2000}%
M. Leventon, ``Statistical Models in Medical Image Analysis,''
{\em Ph.D. dissertation, Massachusetts Inst. Technol., Dept.
Elect. Eng.}, 2000.
%
\bibitem{Leventon00}
M. Leventon, E. Grimson, and O. Faugeras, ``Statistical shape
influence in geodesic active contours", Proc. IEEE Conf. Computer
Vision Pattern Recognition, vol. 1, pp.316-323, 2000.
\bibitem{Shen00}
D. Shen and C. Davatzikos, ``An Adaptive-Focus Deformable Model
Using Statistical and Geometric Information," IEEE Transactions
Pattern Analysis and Machine Intelligence, Vol. 22(8), 2000, pp.
906-913.
\bibitem{Tsai03}
A. Tsai, A. Yezzi, W. Wells, C. Tempany, D. Tucker, A. Fan, E.
Grimson, and A Willsky, ``A Shape-Based Approach to Curve
Evolution for Segmentation of Medical Imagery", IEEE Transactions
on Medical Imaging, Vol. 22, No. 2,  pp.137-154, 2003.
\bibitem{Tsai04}
A. Tsai, W. Wells, C. Tempany, E. Grimson, and A Willsky, ``Mutual
Information in Coupled Multi-Shape Model for Medical Image
Segmentation", Medical Image Analysis, Vol 8, No. 4, pp. 429-445,
2004.
\bibitem{Yang04}
J. Yang and J. Duncan, ``3D image segmentation of deformable
objects with joint shape-intensity prior models using level sets",
Medical Image Analysis, Vol.8, pp. 285-294, 2004.
%
\bibitem{A5}
X. Huang, D. Metaxas, and T. Chen, ``MetaMorphs: Deformable Shape
and Texture Models,'' {\em CVPR} , pp 496--503, Washington DC.,
USA, June 2004.
%
\bibitem{Rouss_2004}%MICCAI 2004
M. Rousson, N. Paragios and R. Deriche, ``Implicit Active Shape
Models for 3D Segmentation in MRI Imaging,'' {\em Medical Image
Computing and Computer Assisted Intervention (MICCAI)}, Part 1, pp
209--216, Saint-Malo, France, September 26-29, 2004.
%
\bibitem{Seb2001}
T. Sebastin, P. Klein, and B. Kimia, ``Recognition of Shapes by
Editting Shock Graphs,''{\em In IEEE International Conference in
Computer Vision}, pp 755--762, 2001.
%
\bibitem{Siddi1999}
K. Siddiqi, A. Shokoufandeh, S. Dickinson, and S. Zucker, ``Shocks
Graphs and Shape Matching,''{\em International Journal of Computer
Vision}, 35:13--32, 1999.
%
\bibitem{Osher_1988}
S. Osher and J. Sethian, ``Fronts Propagating with
Curvature-Dependent Speed: Algorithms Based on the Hamilton-Jacobi
Formulation,''{\em Journal of Computational Physics}, 79:12--49,
1988.
%
\bibitem{A7}
J.A. Sethian, {\em Level Set Methods and Fast Marching Methods,}
Combridge University Press, USA, 1999.
%
\bibitem{A20}
R. Malladi, J. Sethian, and B. Vemuri, ``Shape modeling with front
propagation: A level set approach,'' {\em IEEE Tr. on PAMI},
17(2):158�-175, February, 1995.
%
\bibitem{A8}
J. Gomes and O. Faugeras, ``Reconciling distance functions and
Level-Sets,'', {\em Technical Report 3666}, INRIA, April 1999.
%
\bibitem{A9}
N. Paragios and R. Deriche, ``Unifying boundary and region-based
information for geodesic active tracking,'' {\em CVPR}, Vol. 2, pp
300--305, Fort Collins, Colorado, June 1999.
%
%\bibitem{Baillard:00}
%C. Baillard , C. Barillot , and P. Bouthemy, ``Robust Adaptive
%Segmentation of 3D Medical Images with Level Sets", in {\em
%Technical Report, INRIA, France, November 2000}
%
\bibitem{A11}
X. Zeng, L.H. Staib, R.T. Schultz, H. Tagare, L. Win, and J.S.
Duncan, ``A new approach to 3D sulcal ribbon finding from MR
images,'' {\em MICCAI}, pp 148�-157, sep. 1999.
%
\bibitem{A24}
R. Duda, P. Hart, and D. Stork, ``Pattern Classification,'' {\em
John Wiley and Sons Inc.}, 2001.
%
\bibitem{A18}
C. Samson, L. Blanc-F�raud, G. Aubert, and J. Zerubia,
``Multiphase Evolution and Variational Image Classification,''
{\em Technical Report 3662}, INRIA, France, 1999.
%
\bibitem{Rueckert:99}
D. Rueckert, L. I. Sonoda, C. Hayes, D. L. G. Hill, M. O. Leach,
and D. J. Hawkes, ``Nonrigid Registration Using Free-Form
Deformations: Application to Breast MR Images'', {\em IEEE Tran.
Medical Imaging}, Vol. 18, No. 8, pp. 712--721, 1999
%
\bibitem{Hung}
H. Takizawa, K. Sugiura, M. baba {\em et al.}, ``Analysis of
intracerebral hematoma shapes by numerical computer simulation
using the finite element method,'' Neurol. Med. Chir. Vol. 34, pp.
65-69, 1994
%
\end{thebibliography}
\end{document}
