%TEMPLATE FOR PRODUCING IEEE-FORMAT ARTICLE FOR ICRA98 USING LaTeX.
%written by Matthew Ward, CS Department, Worcester Polytechnic Institute.
%use at your own risk.  Complaints to /dev/null.
%make two column with no page numbering, default is 10 point
%This template can also be downloaded from the Conf.'s WWW site at:
%     Europe              http://www.laas.fr/icra-98/
%     Japan               http://www.iis.u-tokyo.ac.jp/icra-98/
%     United States       http://www.ecn.purdue.edu/icra-98/
\documentstyle[twocolumn]{article}
\pagestyle{empty}

%set dimensions of columns, gap between columns, and space between paragraphs
\setlength{\textheight}{8.75in}
\setlength{\columnsep}{2.0pc}
\setlength{\textwidth}{6.8in}
\setlength{\footheight}{0.0in}
\setlength{\topmargin}{0.25in}
\setlength{\headheight}{0.0in}
\setlength{\headsep}{0.0in}
\setlength{\oddsidemargin}{-.19in}
\setlength{\parindent}{1pc}

%I copied stuff out of art10.sty and modified them to conform to IEEE format
\makeatletter
%as LaTeX considers descenders in its calculation of interline spacing,
%to get 12 point spacing for normalsize text, must set it to 10 points
\def\@normalsize{\@setsize\normalsize{12pt}\xpt\@xpt
\abovedisplayskip 10pt plus2pt minus5pt\belowdisplayskip \abovedisplayskip
\abovedisplayshortskip \z@ plus3pt\belowdisplayshortskip 6pt plus3pt
minus3pt\let\@listi\@listI} 

%need an 11 pt font size for subsection and abstract headings
\def\subsize{\@setsize\subsize{12pt}\xipt\@xipt}

%make section titles bold and 12 point, 2 blank lines before, 1 after
\def\section{\@startsection {section}{1}{\z@}{24pt plus 2pt minus 2pt}
{12pt plus 2pt minus 2pt}{\large\bf}}

%make subsection titles bold and 11 point, 1 blank line before, 1 after
\def\subsection{\@startsection {subsection}{2}{\z@}{12pt plus 2pt minus 2pt}
{12pt plus 2pt minus 2pt}{\subsize\bf}}
\makeatother

\begin{document}
%don't want date printed
\date{}

%make title bold and 14 pt font (LaTeX default is non-bold, 16 pt)
\title{\Large\bf Snakes, Shapes and Gradient Vector Flow}

%for single author (just remove % characters)
\author{Junaed Sattar \\
  School of Computer Science \\
  McGill University \\
  Montreal QC Canada H3A 2A7}
 
\maketitle
%I don't know why I have to reset thispagesyle, but otherwise get page numbers

\thispagestyle{empty}

\subsection*{\centering Abstract}
%IEEE allows italicized abstract
{\em
A novel method for finding active contours, or snakes as developed by Xu and Prince \cite{paper:one} is presented in this paper. The approach uses a regularization based technique and calculus of variations to find what the authors call a Gradient Vector Field or GVF in binary-values or grayscale images. The GVF is in turn applied to 'pull' the snake towards the required feature. The approach presented here differs from other snake algorithms in its ability to extend into object concavities and its robust initialization technique. Although their algorithm works better than existing active contour algorithms, it suffers from computational complexity and associated costs in execution, resulting in slow execution time.
%end italics mode
}
\input{epsf}
\section{Introduction}
Active contours or snakes are of great importance in computer vision. Applications like object tracking, shape modelling, segmentation and edge detection rely heavily on finding the contour of the object in question. Although, problems relating to initialization and convergence in object concavities limit their applications. Xu and Prince \cite{paper:one} present a new approach of snake modelling which they claim removes the abovementioned difficulties in active contour generation. This report discusses their algorithm in detail for binary-valued and grayscale images (section 2) and also discusses some results obtained by implementing and testing their algorithm (section 3). 

\section{Description of the algorithm}
This section discusses the algorithm by Xu and Prince \cite{paper:one} for active contour modelling in details. The outline of the algorithm is as follows. The image on which the snake is to be computed is taken and is converted into a gray-level or binary-valued representation. If the objects of interest are present like line-like structures in the image, the input image can be directly used as the edge map. If the objects are present as homogeous regions whose boundaries separate the regions from the background of different intensity value, an edge map has to be computed. The edge map can be computed using any type of edge detector in the vision or image processing literature. The edge map is normalized to have all edge intensities fall between 0 and 1, if required This normalized edge map is input into the GVF solver. This will produce the GVF field of the edge map. The active contours can be formed by following the direction of the gradient field vectors over a certain number of iterations until a statistical equilibrium is reached.

\subsection{Active Contours}
A traditional snake or active contour is a curve $\textbf{x(s)} =[ x(s),y(s) ], s \in [0,1]$. that moves through the spatial domain of an image to minimize the energy functional 
\begin{equation}
E=\int_0^1 {1 \over 2}(\alpha|x'(s)|^2 + \beta|x''(s)|^2)+E_{ext} (x(s))ds
\end{equation}

Here $\alpha$ and $\beta$ are the weighting parameters that control the snake's tension and rigidity, respectively. $x'(s)$ and $x''(s)$ are the first and second derivatives of $x(s)$ with respect to $s$. The external energy function $E_{ext}$ is derived in such a way from the image that it takes on smaller values at the feature of interests; e.g. at boundary positions. An example of $E_{ext}$ for a gray level image $I(x,y)$ leading to step edges could be
\begin{equation}
E_{ext}(x,y) = -|\nabla (G_\sigma (x,y) * I(x,y)) |
\end{equation}
where $G_\sigma (x,y)$ is a two-dimensional Gaussian function with standard deviation $\sigma$ and $\nabla$ is the gradient operator. Larger values of $\sigma$ will cause the boundaries to become blurry. Even though, large values of $\sigma$ are sometime required to increase the capture range of the active contour, as can been seen from the results(section 3). A snake that minimizes $E$ must satisfy the Euler equation 

\begin{equation}
\alpha x''(s) - \beta x'''(s) - \nabla E_{ext} = 0
\end{equation}

This is equivalent to a force balance equation

\begin{equation}
F_{int}+F_{ext}^{(p)}=0
\end{equation}

where $F_{int}=\alpha x''(s) - \beta x'''(s)$ and $F_{ext}}^{(p)} = -\nabla E_{ext}$. The internal force $F_{int}$ discourages stretching and bending whereas the external potential force $F_{ext}^{(p)}$ pulls the snake towards the desired image edges. To solve for (3), the snake is made dynamic by treating $x$ as a function of time $t$ as well as $s$. Then the partial derivative of $x$ with respect to $t$ is set as follows:

\begin{equation}
x_t (s,t) =  \alpha x''(s,t) - \beta x'''(s,t) - \nabla E_{ext}
\end{equation}
When statistical equilibrium for $x(s,t)$ is reached, the term $x_t (s,t)$ vanishes and a solution to equation (3) is obtained.

\subsection{Problems with Snake Formation}
Applications of snakes in computer vision and image processing has had their limitations. There are two main reasons for that. One, the poor convergence property of the active contours. Specifically, concavities in the object of interest are rarely covered; i.e. the snake does not extend to the concave regions of the object. An example of that is shown in figure 1. The snake was obtained by following the general algorithm outlined in the previous subsection.
\begin{figure}
\epsfysize=2in
\epsfxsize=3in
\epsffile{fig1.eps}
\caption{Traditional snake and force field after 500 iterations}
\label{Figure 1}
\end{figure}
As can be seen, the steady state active contour fails to reach inside the concave region of the U-shaped object. The second problem with snakes is the limited capture range, which is related to the initialization of the snake around the object of interest. The magnitude of external forces die out far away from the boundaries. This phenomenon might fail convergence of the snake. Increasing the value of $\sigma$ of the Gaussian might improve the range, but at higher values, the edge will get blurry and eventually the concavitiy will be obliterated. Cohen and Cohen \cite{paper:two} proposed a distance potential force model to get around this problem. The general idea behind this model is to have large external forces far away from the boundaries of the object, thus increasing the capture range of the snake. Even then, snakes based on this model fails to converge to concavities, due to horizontal forces pulling the snake apart but not downward into the concavities. The technique presented by Xu and Prince \cite{paper:one} addresses these issues and presents a new formulation for active contour modelling. 

\subsection{Gradient Vector Fields}
The poor convergence of the snake can be shown to happen due to convergence of the solution to a local minimum. The authors present a solution to the problem by replacing the standard external force $F_{ext}$ in the force balance equation (4) with a static external force which does not change with time or depend on the position of the snake itself. This new static external external force field $F_{ext}^{(g)}=\textbf{v}(x,y)$ is called the Gradient Vector Field or GVF. Replacing the external potential force $-\nabla E_{ext}$ in (5) with \textbf{v} yields the following equation:

\begin{equation}
x_t (s,t) =  \alpha x''(s,t) - \beta x'''(s,t) + \textbf{v} 
\end{equation}

The parametric curve solving the above dynamic equation is termed as a \emph{GVF snake}. Standard numerical methods can be used to solve for this equation and yield the GVF snake. 

\subsection{GVF Snake Formation}
The process starts by calculating the edge map of the given image, using any edge finding algorithm from the image processing literature. The edge map has three imporant features realting to snake formation. One, the gradient of this edge map has vectors pointing towards the edges which is a desirable property for snakes. Two, these vectors have large magnitude at the vicinity of the edges. Three, in homogenous regions (regions with little variation in image intensity) $\nabla f$ is almost zero. The second and third features can be problematic when constructing an active contour. To keep the first feature and nullify the effect of the latter two, the gradient map is extended farther away from the edges and into homogenous regions using a computational diffusion process. 

The gradient vector flow field is defined as the vector field $\textbf{v}(x,y)=(u(x,y), v(x,y))$ that minimizes the following enerygy functional

\begin{equation}
\varepsilon = \int \int \mu( u_x^2 + u_y^2 + v_x^2 + v_y^2 ) + |\nabla f|^2|\textbf{v} - \nabla f|^2 dxdy
\end{equation}

As can be seen, this is an example of variational formulation of regularization. The parameter $\mu$ is a regularizing parameter which adjusts the tradeoff between the first and second terms of the integrand and is set according to the level of noise present in the image. Also, where the value of the edge gradient is small, energy is dominated by the sum of the partial derivatives of the gradient field. When the gardient is large, the second term dominates. In this case, setting $\textbf{v}=\nabla f$ minimizes the energy. Overall, this formulation transforms the gradient vector flow field by keeping it equal to the edge gradient at the boundaries; it also keeps $\textbf{v}$ slowly varying at the homogenous regions of the image. Using the calculus of variations, it can be shown that the GVF field can be found by solving the pair of Euler equations stated below:
\begin{equation}
\mu \nabla ^2 u -(u-f_x)(f_x^2 + f_y^2) = 0
\end{equation}
\begin{equation}
\mu \nabla ^2 v -(v-f_y)(f_x^2 + f_y^2) = 0
\end{equation}

Here, $\nabla ^2$ is the Laplacian operator. Equations (8) and (9) can be solved numerically by treating $u$ and $v$ as a function of time. The resulting equations are: 
\begin{eqnarray}
u_t (x,y,t)&=&\mu \nabla ^2 u_t (x,y,t)-(u(x,y,t)- f_x(x,y))\nonumber\\
&&\cdot (f_x(x,y)^2 + f_y(x,y)^2)
\end{eqnarray}

\begin{eqnarray}
v_t (x,y,t)&=&\mu \nabla ^2 v_t (x,y,t)-(v(x,y,t) - f_y(x,y))\nonumber\\
&&\cdot (f_x(x,y)^2 + f_y(x,y)^2)
\end{eqnarray}
The steady-state solution of equation (10) and (11) yields the solution to the Euler equations (8) and (9). An iterative solution can be set up for solving the equations above.

\section{Experiments and Results}
The algorithm for snake generation and active contour modelling was implemented in Matlab and tested on images obtained from the author's websites and also on images obtained by myself. The images contained both binary valued and grayscale images. All edge detection was performed using the Canny edge detection algorithm.

\subsection{Convergence to Boundary Concavity}
An implementation of the GVF snake algorithm was first run on a binary valued 64x64 pixel image containing a U-shaped object as shown in figure 2. The gradient maps, the GVF of the flow vector and the final GVF snake is also shown in the figure. 
\begin{figure}
\epsfysize=3in
\epsfxsize=3in
\epsffile{u1.eps}
\caption{GVF Snake Formation}
\label{Figure 2}
\end{figure}
As can be seen from the outputs, the algorithm correctly converges to the object boundary, even in the concavity in the given shape. The algoritm was next run on a binary valued image with dimensions of 64x64 pixels containing a non-convex polygon. The results can be seen in figure 3. Initialization of the snake, the gradient maps and the gradient vector fields for the given curve is also shown. The snake succeeded in converging to the boundary positions in this example too.
The next run of the implementation was performed on a gray level image. The edge map was first computed using the Canny edge operator as mentioned before. No threshholding was performed on the image to binarize the edge map. The results are shown in figure 4. The results show that the snake was able to latch on to the correct object boundaries on this test image as well. 
\begin{figure}
\epsfysize=3in
\epsfxsize=3in
\epsffile{fig3.eps}
\caption{GVF Snake Formation}
\label{Figure 3}
\end{figure}


\subsection{Insensitivity to Initialization}
The GVF implementation was tested alongside an implementation of a traditional snake algorithm. The traditional snake algorithm, as can be seen from figure 5, has a limited capture range and is also very sensitive to the way the snake contour is initialized. Particularly, in cases where the snake is initialized accross the object, part of the snake falls inside and part outside the object of interest. In such cases, the algorithm gets confused about the direction along which it should progress the snake, and converges to a highly undesirable configuration.
\begin{figure}
\epsfysize=3in
\epsfxsize=3in
\epsffile{fig4.eps}
\caption{Graylevel GVF Snake with across-the-object initialization}
\label{Figure 4}
\end{figure}

\begin{figure}
\epsfysize=2in
\epsfxsize=3.5in
\epsffile{trad1.eps}
\caption{Traditional Snake with across-the-object initialization}
\label{Figure 5}
\end{figure}

The GVF snake, on the other hand, follows the gradient vector flow field to reach a point of steady state configuration. The flow vectors both inside and outside of the object point to the object boundaries. The GVF snake thus eventually converges correctly along the object boundary. The tests show that GVF has insensitivity to initialization as well as ability to progress into bounary concavities.

\begin{figure}
\epsfysize=3in
\epsfxsize=3in
\epsffile{fig8.eps}
\caption{GVF Snake on Noisy Image}
\label{Figure 6}
\end{figure}

\subsection{Noise Sensitivity}
The binary valued images on which the snake algorithm was tested came out to be quite robust upto a certain level of noise. After being corrupted by Gaussian white noise of zero mean and a small variance, the GVF snake algorithm was run on the images. The result can be shown in figure 6. The grayscale images took a larger value of sigma to reduce the effect of noise on the snake algorithm.  

\subsection{Performance}
The single biggest drawback of the GVF snake algorithm was the speed of exeucution. The algorithm was implemented in Matlab for Linux (kernel 2.4.22) running on an AMP Athlon XP 1.4 GHz processor with 512 megabytes of RAM. For the binary-valued images (i.e. the image containing the U-shaped object) the algorithm took somewhere around the 50 second mark to complete the execution. The standard snake algorithm, on the other hand took around 11 seconds, although it failed to correctly converge to the boundary concavities. For grayscale images, the time taken was a lot more (almost around 90 seconds). The slowness could be due to the lack of optimization in the Matlab execution engine.Even then, this drawback makes the algorithm unsuitable for real time applications like tracking, inspite of its robustness and felxibility. The authors nevertheless claim that the speed of the execution could be improved by using native optimized C code and by using an algorithm optimization method as the multigrid method.

\subsection{Higher-Dimensional GVF}
The GVF snake algorithm can be extended to higher dimension and applied to 3-D objects. The paper presents an example of computing the contour of a 3-D object by using the GVF snake, although this particular case was not tested in the implementation. 

\section{Discussion}
The gradient vector flow based active contour generation algorithm by Xu and Prince \cite{paper:one} was presented. The algorithm, along with a traditional snake generation algorithm was implemented in Matlab and tested on various images, both grayscale and binary valued. It was found that the algorithm succeeds in convergin the active contour to boundary concavitie in both types of images, even with the presence of noise. The drawback of the methods is its execution speed. Inspite of its robustness to initialization and increased capture range, the algorithm takes a long time to converge to object contours. Also, the value of the regulrization parameter $\mu$ and the snake parameters $\alpha$ and $\beta$ are set manually. The paper does not adress the issue of setting these parameters according to a particular image. 
The technique presented by Xu and Prince is attractive in terms of accuracy and robustness. It has found applications in medical imaging and offline tracking. With these issues in mind, the approach may not be suitable for real time tracking or contour modelling in unconstrained (e.g. excessively noisy) environments. 

\begin{thebibliography}{9}
\bibitem{paper:one}
Chengyang Xu and Jerry L. Prince.
``Snakes, Shapes and Gradient Vector Flow.''
{\em IEEE Transactions on Image Prccessing}, Vol. 7, No. 3, March 1998.

\bibitem{paper:two}
L. D. Cohen and I. Cohen.
``Finite-element methods for active contour models and balloons for 2-D and 3-D images.''
{\em IEEE Transactions on Pattern Analysis and Machine Intelligence}, 15(11):1131-1147, November 1993.

\bibitem{paper:three}
Chengyang Xu and Jerry L. Prince. 
``Gradient Vector Flow: A New External Force for Snakes.''
{\em Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition}, 66-71, March 1997.

\bibitem{paper:four}
B. K. P. Horn and B. G. Schunck.
``Determining Optical Flow.''
{\em Artificial Intelligence}, 17:185-203, 1981.

\end{thebibliography}

\end{document}

