%======================================================================%
% Filename: Article.tex
% Project: Insight Journal Article - Minimal path
% Author: Dan Mueller - d.mueller[at]qut.edu.au, dan.muel[at]gmail.com
% Date: $Date$
% Revision: $Revision$
%----------------------------------------------------------------------%
% Description:
% This artcile is for the Insight Journal.
% It describes minimal path extraction using FastMarching action maps.
%======================================================================%
\documentclass{InsightArticle}

%\usepackage{courier}                 %Font
%\usepackage{charter}                 %Font
%\usepackage{times}                   %Font
%\usepackage{newcent}                 %Font
%\usepackage{pslatex}                 %Font
%\usepackage{lmodern}                 %Font
%\usepackage{palatino}                %Font
%\usepackage{bookman}                 %Font
\usepackage{amsfonts}                 %Maths fonts
\usepackage{amssymb}                  %Maths symbols
\usepackage{listing}                  %Used for list of listings
\usepackage{listings}                 %List program code source
\usepackage{subfigure}                %For subfigures
\usepackage{color}                    %Allow color
\usepackage{fancyhdr}                 %Use headers and footers
\usepackage{graphicx}                 %Graphics
\usepackage{ccaption}
\usepackage[dvipdfmx=true,
            bookmarks=true,
            bookmarksopen=true,
            bookmarksopenlevel=2,
            bookmarksnumbered=true,
            backref=section,
            colorlinks=true,
            linkcolor={blue},
            citecolor={blue},
            urlcolor={blue}]{hyperref}

%Set default font
%\renewcommand{\familydefault}{\rmdefault}  %Roman
\renewcommand{\familydefault}{\sfdefault}   %Sans serif
%\renewcommand{\familydefault}{\ttdefault}  %Typewriter

%Keywords
\newcommand{\keywords}[1]{\noindent\textbf{Keywords:}\it{\ #1}}

% TOC depth
% \setcounter{tocdepth}{2}

%Stuff for images
\graphicspath{{images/}}

%Create indent for description
\newenvironment{indentdescription}
               {\list{}{\leftmargin \leftmargin}%
                \item\relax}
               {\endlist}
               
%Allow symbol footnotes
\long\def\symbolfootnote[#1]#2{\begingroup%
\def\thefootnote{\fnsymbol{footnote}}\footnote[#1]{#2}\endgroup}

% Center last line of captions
%\captionstyle{\centerlastline}

%Allow code listings
\input{listcode}                 

%Create conditional to include or exclude figures 
%(saves time during debug compilation).
\makeatletter
\newif\if@includefigs@
\@includefigs@true

%======================================================================%
%                   F r o n t   M a t t e r                            % 
%======================================================================%
% The title should be descriptive enough for people to be able to find
% the relevant document.
\title{Fast Marching Minimal Path Extraction in ITK}

% Increment the release number whenever significant changes are made.
% The author and/or editor can define 'significant' however they like.
\release{1.0}

% Name and affiliation
\author{Dan Mueller$^{\scriptscriptstyle\textsf{1}}$}
\authoraddress{$^{\textsf{1}}$Queensland University of Technology, Brisbane, Australia}

% Set custom header style
\pagestyle{fancy}                               
\fancyhead{}            % clear all header fields
\fancyhead[L]{\textsf{Fast Marching Minimal Path Extraction in ITK}}
\fancyhead[R]{\textsf{\thepage}}
\fancyfoot{}            % clear all footer fields

% Document
\begin{document}

% Show the title
\maketitle

% Make the title page header/footer empty
\thispagestyle{empty}

%======================================================================%
%                       A b s t r a c t                                % 
%======================================================================%
\begin{abstract}
\noindent
This paper describes the ITK implementation of a minimal path extraction
framework based on Fast Marching arrival functions.
The method requires the user to provide three inputs:
1. a meaningful speed function to generate an arrival function,
2. path information in the form of start, end, and way-points
(which the path must pass near),
and 3. an optimizer which steps along the resultant arrival function
perpendicular to the Fast Marching front.
A number of perspectives for choosing speed functions and optimizers
are given, as well as examples using synthetic and real images.
%
\newline
\newline
\keywords{minimal path, centerline, vessel segmentation, ITK}

\end{abstract}

%======================================================================%
%                     M a i n   M a t t e r                            % 
%======================================================================%
\tableofcontents

%----------------------------------------------------------------------%
% Section
%----------------------------------------------------------------------%
\newpage
\section{Background}
%
Minimal path extraction is useful for a range of application domains
including medical image analysis, robot navigation, and 
artificial intelligence (eg. emulating human behaviour in computer games).
A number of optimisation approaches have been employed within these
problem domains, such as 
tracking methods \cite{Aylward2002a,Wink2000a},
minimal cost methods (ie. Dijkstra's algorithm) \cite{Wink2004a},
and
A* (pronounced A star) \cite{WikipediaAStar2007}.

Sethian \cite[pp. 284–-312]{Sethian1999a} and Andrews et al. \cite{Andrews2007a}
explored the use of Fast Marching arrival functions to extract minimal paths.
This method relies on the fact that the gradient of the Fast Marching arrival
function $T$ is orthogonal to the wave front.
Furthermore, $T$ has only one local minimum, with is guaranteed
to be the global minimum (see \cite[pp.17]{Deschamps2001a}).
Therefore the minimal path can be extracted by back-propagating
from a given seed (corresponding to the end point of the desired path)
to the starting point implicitly embedded in $T$
(see \figurename~\ref{fig:Background:Synthetic:01}).
The back-propagation can be accomplished using the existing ITK optimizer
framework (see Section~\ref{sec:Implementation:Optimizer}).

The advantages of Fast Marching minimal path extraction are three-fold:
%
\begin{enumerate}
	\item Fast Marching is computationally efficient with $O\left( n\: \log\: n \right)$ complexity.
  \item The resultant path is continuous and directed (compare with skeletonization methods
        which return a set of unordered discrete points).
  \item To the best of this author's knowledge, the generic method is not covered by any patents
        (however, please note that the framework \emph{can} be used to implement
        vessel extraction methods which are covered by patents,
        as demonstrated in Section~\ref{sec:Implementation:SpeedFunction}).
\end{enumerate}

\if@includefigs@
\begin{figure}[htb]
    \centering
    \subfigure[Path]{\label{fig:Background:Synthetic:01:Speed}
    \includegraphics[width=0.4\textwidth]{Synthetic-01-Image.png}}
    \hspace{2mm}
    \subfigure[Arrival]{\label{fig:Background:Synthetic:01:Arrival}
    \includegraphics[width=0.4\textwidth]{Synthetic-01-Arrival.png}}
    \caption{\label{fig:Background:Synthetic:01}%
             The minimal path is extracted by back-propagating from an end point
             to a start point embedded in the arrival function.
             Because the Fast Marching front moves orthogonal to itself,
             tracking perpendicular to the arrival function is guaranteed to
             follow the geodesic (minimal) path.}
\end{figure}
\fi

%----------------------------------------------------------------------%
% Section
%----------------------------------------------------------------------%
\section{Implementation}
%
%----------------------------------------------------------------------%
\subsection{Overview}
%
This project implements a number of auxiliary functions,
however the main filter expected to be employed by the user is
\code{itk::SpeedFunctionToPathFilter}.
This filter is a subclass of \code{itk::ArrivalFunctionToPathFilter},
which provides the basic functionality to convert an arrival function to a path.
The main filter provides the added functionality of computing
appropriate arrival functions from a given speed function,
which are then feed to the subclass for processing.
\figurename~\ref{fig:Implementation:Framework} depicts an overview
of \code{itk::SpeedFunctionToPathFilter}.

\code{itk::SpeedFunctionToPathFilter}
expects path information (consisting of start, end, and optional way-points)
and a speed function (a real-valued image in the range $[0,1]$).
Starting with the first way-point (or end point if there are no way-points)
a front is propagated in the typical Fast Marching manner
(each point from which a front is propagated is called a \emph{Trial} point).
The front terminates when both the previous and next points in the path have been reached
(this is accomplished using the \code{itk::FastMarchingUpwindGradientImageFilter}).
The resultant arrival function is viewed as a cost function using a new class
called \code{itk::SingleImageCostFunction},
and is optimized using an appropriate \code{itk::SingleValuedNonLinearOptimizer}.
The optimizer is initialised with the next point in the path list,
and at each step the cost function (ie. arrival function) is minimised to step
closer to the local (and global) minimum (ie. the \emph{Trial} point).

\code{itk::SpeedFunctionToPathFilter} monitors the optimizer \code{Iteration}
event, and saves each position to the current path.
This means the filter only supports step-based optimizers,
such as \code{itk::GradientDescentOptimizer}.
Once the minimum of the current arrival function has been reached,
the process is repeated with the next point in the list of path information.
Once the path list has been exhausted, the process is repeated for each
path (start, end, way-points) specified by the user.
An extracted path is stored in a separate filter output
(ie. specifying multiple path information objects generates multiple outputs).

At the moment paths are represented using \code{itk::PolyLineParametricPath}, 
which stores each point as an \code{itk::ContinuousIndex}.
This situation is perhaps not ideal because the path points should really
be stored in physical space (ie. using \code{itk::Point}s).
Additionally, the \code{itk::Path} framework does not
support input/output operations (the user is required to implement reading
and writing to files manually).

\if@includefigs@
\begin{figure}[!b]
    \centering
    \includegraphics[width=0.70\textwidth]{Path-Framework-01.eps}
    \caption{\label{fig:Implementation:Framework}%
    An overview of the internal workings of \code{itk::SpeedFunctionToPathFilter}.}
\end{figure}
\fi

%----------------------------------------------------------------------%
\subsection{Path Information}\label{sec:Implementation:PathInformation}
%
Fast Marching minimal path extraction is a semi-autonomous segmentation
method --- the user is required to provide start and end points.
This is accomplished using an internal class \code{PathInfo}
(see \code{itkSpeedFunctionToPathFilter.h}).
The following code fragment demonstrates how the set the start and
end points, and how to add a single way-point which the path
must pass near.

\listcpluspluspsnip
\lstinputlisting[float, escapechar=\|,
				         label=lst:PathInfo,
				         linerange={27-29,31-31,35-36,73-85}]
                 {../Source/examples.cxx}

Multiple \code{PathInfo} objects can be added to the 
\code{itk::SpeedFunctionToPathFilter} instance using the
\code{AddPathInfo()} method.
Previously added path information objects added to the filter can be cleared using the
\code{ClearPathInfo()} method.

Because the user can not be expected to exactly specify points on the path,
the \code{TerminationValue} parameter provides a configurable buffer
which protects against the path needing to directly pass through each point.
The optimizer is terminated when the current arrival value is less than 
\code{TerminationValue};
the smaller the value, the closer the path will get to each user specified point.
The default value is 1.0
(note the parameter is specified in arrival time units, not physical spacing).
To prevent oscillations,
it is recommended that the specified optimizer has a small step size when
\code{TerminationValue} is small.

%----------------------------------------------------------------------%
\subsection{Speed Function}\label{sec:Implementation:SpeedFunction}
%
Choosing an appropriate speed function is the most difficult part
of the entire process (similar to choosing a speed function when performing
active contour segmentation).
Recall that the speed function must be a real-valued image
(ie. \code{float} or \code{double}) in the range $[0,1]$.
For computational efficiency it is desirable to have the speed function
as close to $1.0$ near the desired path, and $0.0$ elsewhere.
In some situations this may not be possible,
in which case there will be a computational penalty because the Fast
Marching front will visit most (if not all) of the speed function.

TODO: Discuss efficient method for generating speed function for vessel segmentation.

When using the framework to extract the medial axis of a vessel,
it may be necessary to use a method to center the path.
Deschamps \cite{Deschamps2001a} proposed one such method which
firstly used the method to extract a rough centerline,
performed segmentation using the rough centerline,
then computed a speed function using the signed distance transform of the
resultant binary volume, and repeated the minimal path extraction.
This method is depicted in \figurename~\ref{fig:Speed:RealDSA01}.
It should be noted that this method (ie. extracting centered paths using a distance
transform) may be covered by a 
\href{http://math.lbl.gov/~deschamp/html/publications.html#SECT3}{patent}.
As can be seen, because the minimal path extraction framework described in this
paper is capable of implementing a range of methods,
care must be taken to ensure the framework is not used to implement a
patented method.

\if@includefigs@
\begin{figure}[htb]
    \centering
    \subfigure[Uncentered Speed]{
    \includegraphics[width=0.35\textwidth]{Real-DSA-01-Uncentered.png}}
    \hspace{1mm}
    \subfigure[Centered Speed]{
    \includegraphics[width=0.35\textwidth]{Real-DSA-01-Centered.png}}
    \subfigure[Overlay]{
    \includegraphics[width=0.7225\textwidth]{Real-DSA-01-Both.eps}}
    \caption{\label{fig:Speed:RealDSA01}%
    A digital subtraction angiography (DSA) image of the cerebral vessels.
    The resultant paths have been laid over the original image:
    the minimal path (green) and centered path (yellow).
    Original image courtesy of
    \href{http://upload.wikimedia.org/wikipedia/en/6/63/Cerebral_Angiogram_Lateral.jpg}{Wikipedia}.}
\end{figure}
\fi

%----------------------------------------------------------------------%
\subsection{Optimizer}\label{sec:Implementation:Optimizer}
%
At present there are three step-wise optimizers which function together
with the given minimal path framework:
\code{GradientDescentOptimizer},
\code{RegularStepGradientDescentOptimizer}, and
\code{IterateNeighborhoodOptimizer}.
\code{IterateNeighborhoodOptimizer} is a new optimizer which operates
by iterating the neighbors of the current position, and moving to
the minimum value.

The gradient descent methods are suitable for minimizing continuous
cost functions (ie. functions without discontinuities).
In comparison, the neighborhood method is suitable for minimizing cost functions
with large discontinuities
(ie. when pixels visited by the arrival function are directly adjacent to non-visited pixels).
The gradient descent optimizers are not suitable for discontinuous cost functions
because the step vector, which is computed using the gradient, is skewed in regions
direct adjacently to non-visited pixels.
Examples of continuous and discontinuous cost
functions are depicted in \figurename~\ref{fig:Optimizers:Continuous}.
It is anticipated that for most medical imaging applications
--- in which the speed function is continuous --- the gradient descent
optimizers will be most suitable.

\if@includefigs@
\begin{figure}[htb]
    \centering
    \subfigure[Continuous Arrival]{
    \includegraphics[width=0.45\textwidth]{Synthetic-02-Arrival.png}}
    \hspace{1mm}
    \subfigure[Paths]{
    \includegraphics[width=0.45\textwidth]{Synthetic-02-Optimizers.png}}
    \subfigure[Discontinuous Arrival]{
    \includegraphics[width=0.45\textwidth]{Synthetic-03-Arrival.png}}
    \hspace{1mm}
    \subfigure[Paths]{
    \includegraphics[width=0.45\textwidth]{Synthetic-03-Paths.png}}    
    \caption{\label{fig:Optimizers:Continuous}%
    The top images (a and b) depict a continuous arrival function
    and extracted paths using the gradient descent (blue) and
    iterate neighborhood (red) optimizers for a synthetic dataset
    (\code{Data/Synthetic-02}).
    The bottom image (c) depicts a discontinuous arrival function
    (note that visited pixels are directly adjacent to non-visited pixels).    
    The maze image (\code{Data/Synthetic-03}) is courtesy of
    \href{http://en.wikipedia.org/wiki/Image:Cg_pp_maze.png}{Wikipedia}.}
\end{figure}
\fi

%----------------------------------------------------------------------%
% Section
%----------------------------------------------------------------------%
\section{Examples}
The code from these examples can be found in \code{Source/examples.cxx}.

\subsection{Using the Gradient Descent Optimizer}
%
The first step for using the proposed framework is to setup the
relevant type definitions:
\listcpluspluspsnip
\lstinputlisting[escapechar=\|,
				         linerange={28-29,31-31,35-36}]
                 {../Source/examples.cxx}
Next create a cost function:
\lstinputlisting[escapechar=\|,
				         linerange={58-59}]
                 {../Source/examples.cxx}
Now we create and setup the desired optimizer,
noting that the \code{LearningRate} parameter will default to \code{1.0}:
\lstinputlisting[escapechar=\|,
				         linerange={63-65}]
                 {../Source/examples.cxx}
Finally, we plug the cost function and optimizer into a newly created
path filter object:
\lstinputlisting[escapechar=\|,
				         linerange={68-72}]
                 {../Source/examples.cxx}
The path information is specified as described in
Section~\ref{sec:Implementation:PathInformation} and the filter
updated to compute the path.

\subsection{Using the Regular Step Gradient Descent Optimizer}
%
The regular step gradient descent optimizer gives the user more control
over the step size, via the \code{SetMaximumStepLength},
\code{SetMinimumStepLength}, and \code{SetRelaxationFactor} methods.
In the following example, the step length is configured to be roughly \code{0.5}
but can vary as small as \code{0.1}:
\lstinputlisting[escapechar=\|,
				         linerange={167-172}]
                 {../Source/examples.cxx}

%----------------------------------------------------------------------%
% Conclusion
%----------------------------------------------------------------------%
\section{Conclusion}
%
This paper has described the implementation of Fast Marching minimal
path extraction for ITK.
The framework allows users to specify a speed function,
from which a Fast Marching arrival function is computed
and the geodesic (minimal) path between a number of points is extracted.
The proposed implementation uses the existing optimizer framework within ITK
to perform the minimal path extraction.
For suggestions or bugs, feel free to contact 
the author\footnote{Corresponding author: Dan Mueller: %
\href{mailto:d.mueller@qut.edu.au}{d.mueller@qut.edu.au} or
\href{mailto:dan.muel@gmail.com}{dan.muel@gmail.com}.
}.

%======================================================================%
%                    B i b l i o g r a p h y                           % 
%======================================================================%
\newpage
\bibliographystyle{plain}
\bibliography{article}

\end{document}
