final

[these_gilles.git] / DOCS / paper_snake_gpu / snake_gpu.tex~

%% bare_conf.tex
%% V1.3
%% 2007/01/11
%% by Michael Shell
%% See:
%% http://www.michaelshell.org/
%% for current contact information.
%%
%% This is a skeleton file demonstrating the use of IEEEtran.cls
%% (requires IEEEtran.cls version 1.7 or later) with an IEEE conference paper.
%%
%% Support sites:
%% http://www.michaelshell.org/tex/ieeetran/
%% http://www.ctan.org/tex-archive/macros/latex/contrib/IEEEtran/
%% and
%% http://www.ieee.org/

%%*************************************************************************
%% Legal Notice:
%% This code is offered as-is without any warranty either expressed or
%% implied; without even the implied warranty of MERCHANTABILITY or
%% FITNESS FOR A PARTICULAR PURPOSE! 
%% User assumes all risk.
%% In no event shall IEEE or any contributor to this code be liable for
%% any damages or losses, including, but not limited to, incidental,
%% consequential, or any other damages, resulting from the use or misuse
%% of any information contained here.
%%
%% All comments are the opinions of their respective authors and are not
%% necessarily endorsed by the IEEE.
%%
%% This work is distributed under the LaTeX Project Public License (LPPL)
%% ( http://www.latex-project.org/ ) version 1.3, and may be freely used,
%% distributed and modified. A copy of the LPPL, version 1.3, is included
%% in the base LaTeX documentation of all distributions of LaTeX released
%% 2003/12/01 or later.
%% Retain all contribution notices and credits.
%% ** Modified files should be clearly indicated as such, including  **
%% ** renaming them and changing author support contact information. **
%%
%% File list of work: IEEEtran.cls, IEEEtran_HOWTO.pdf, bare_adv.tex,
%%                    bare_conf.tex, bare_jrnl.tex, bare_jrnl_compsoc.tex
%%*************************************************************************

% *** Authors should verify (and, if needed, correct) their LaTeX system  ***
% *** with the testflow diagnostic prior to trusting their LaTeX platform ***
% *** with production work. IEEE's font choices can trigger bugs that do  ***
% *** not appear when using other class files.                            ***
% The testflow support page is at:
% http://www.michaelshell.org/tex/testflow/



% Note that the a4paper option is mainly intended so that authors in
% countries using A4 can easily print to A4 and see how their papers will
% look in print - the typesetting of the document will not typically be
% affected with changes in paper size (but the bottom and side margins will).
% Use the testflow package mentioned above to verify correct handling of
% both paper sizes by the user's LaTeX system.
%
% Also note that the "draftcls" or "draftclsnofoot", not "draft", option
% should be used if it is desired that the figures are to be displayed in
% draft mode.
%
\documentclass[10pt, conference, compsocconf]{IEEEtran}
% Add the compsocconf option for Computer Society conferences.
%
% If IEEEtran.cls has not been installed into the LaTeX system files,
% manually specify the path to it like:
% \documentclass[conference]{../sty/IEEEtran}


%  \usepackage[latin1]{inputenc}
%  \usepackage[cyr]{aeguill}
%  \usepackage[francais]{babel}


% Some very useful LaTeX packages include:
% (uncomment the ones you want to load)


% *** MISC UTILITY PACKAGES ***
%
%\usepackage{ifpdf}
% Heiko Oberdiek's ifpdf.sty is very useful if you need conditional
% compilation based on whether the output is pdf or dvi.
% usage:
% \ifpdf
%   % pdf code
% \else
%   % dvi code
% \fi
% The latest version of ifpdf.sty can be obtained from:
% http://www.ctan.org/tex-archive/macros/latex/contrib/oberdiek/
% Also, note that IEEEtran.cls V1.7 and later provides a builtin
% \ifCLASSINFOpdf conditional that works the same way.
% When switching from latex to pdflatex and vice-versa, the compiler may
% have to be run twice to clear warning/error messages.






% *** CITATION PACKAGES ***
%
\usepackage{cite}
% cite.sty was written by Donald Arseneau
% V1.6 and later of IEEEtran pre-defines the format of the cite.sty package
% \cite{} output to follow that of IEEE. Loading the cite package will
% result in citation numbers being automatically sorted and properly
% "compressed/ranged". e.g., [1], [9], [2], [7], [5], [6] without using
% cite.sty will become [1], [2], [5]--[7], [9] using cite.sty. cite.sty's
% \cite will automatically add leading space, if needed. Use cite.sty's
% noadjust option (cite.sty V3.8 and later) if you want to turn this off.
% cite.sty is already installed on most LaTeX systems. Be sure and use
% version 4.0 (2003-05-27) and later if using hyperref.sty. cite.sty does
% not currently provide for hyperlinked citations.
% The latest version can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/cite/
% The documentation is contained in the cite.sty file itself.






% *** GRAPHICS RELATED PACKAGES ***
%
\ifCLASSINFOpdf
  \usepackage[pdftex]{graphicx,color}
  % declare the path(s) where your graphic files are
  \graphicspath{{img/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  \DeclareGraphicsExtensions{.pdf,.jpeg,.png}
\else
  % or other class option (dvipsone, dvipdf, if not using dvips). graphicx
  % will default to the driver specified in the system graphics.cfg if no
  % driver is specified.
  % \usepackage[dvips]{graphicx}
  % declare the path(s) where your graphic files are
  % \graphicspath{{../eps/}}
  % and their extensions so you won't have to specify these with
  % every instance of \includegraphics
  % \DeclareGraphicsExtensions{.eps}
\fi
% graphicx was written by David Carlisle and Sebastian Rahtz. It is
% required if you want graphics, photos, etc. graphicx.sty is already
% installed on most LaTeX systems. The latest version and documentation can
% be obtained at: 
% http://www.ctan.org/tex-archive/macros/latex/required/graphics/
% Another good source of documentation is "Using Imported Graphics in
% LaTeX2e" by Keith Reckdahl which can be found as epslatex.ps or
% epslatex.pdf at: http://www.ctan.org/tex-archive/info/
%
% latex, and pdflatex in dvi mode, support graphics in encapsulated
% postscript (.eps) format. pdflatex in pdf mode supports graphics
% in .pdf, .jpeg, .png and .mps (metapost) formats. Users should ensure
% that all non-photo figures use a vector format (.eps, .pdf, .mps) and
% not a bitmapped formats (.jpeg, .png). IEEE frowns on bitmapped formats
% which can result in "jaggedy"/blurry rendering of lines and letters as
% well as large increases in file sizes.
%
% You can find documentation about the pdfTeX application at:
% http://www.tug.org/applications/pdftex





% *** MATH PACKAGES ***
%
%\usepackage[cmex10]{amsmath}
% A popular package from the American Mathematical Society that provides
% many useful and powerful commands for dealing with mathematics. If using
% it, be sure to load this package with the cmex10 option to ensure that
% only type 1 fonts will utilized at all point sizes. Without this option,
% it is possible that some math symbols, particularly those within
% footnotes, will be rendered in bitmap form which will result in a
% document that can not be IEEE Xplore compliant!
%
% Also, note that the amsmath package sets \interdisplaylinepenalty to 10000
% thus preventing page breaks from occurring within multiline equations. Use:
%\interdisplaylinepenalty=2500
% after loading amsmath to restore such page breaks as IEEEtran.cls normally
% does. amsmath.sty is already installed on most LaTeX systems. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/amslatex/math/





% *** SPECIALIZED LIST PACKAGES ***
%
\usepackage[ruled,lined,linesnumbered]{algorithm2e}
%\usepackage{algorithmic}
% algorithmic.sty was written by Peter Williams and Rogerio Brito.
% This package provides an algorithmic environment fo describing algorithms.
% You can use the algorithmic environment in-text or within a figure
% environment to provide for a floating algorithm. Do NOT use the algorithm
% floating environment provided by algorithm.sty (by the same authors) or
% algorithm2e.sty (by Christophe Fiorio) as IEEE does not use dedicated
% algorithm float types and packages that provide these will not provide
% correct IEEE style captions. The latest version and documentation of
% algorithmic.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithms/
% There is also a support site at:
% http://algorithms.berlios.de/index.html
% Also of interest may be the (relatively newer and more customizable)
% algorithmicx.sty package by Szasz Janos:
% http://www.ctan.org/tex-archive/macros/latex/contrib/algorithmicx/




% *** ALIGNMENT PACKAGES ***
%
\usepackage{array}
% Frank Mittelbach's and David Carlisle's array.sty patches and improves
% the standard LaTeX2e array and tabular environments to provide better
% appearance and additional user controls. As the default LaTeX2e table
% generation code is lacking to the point of almost being broken with
% respect to the quality of the end results, all users are strongly
% advised to use an enhanced (at the very least that provided by array.sty)
% set of table tools. array.sty is already installed on most systems. The
% latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/required/tools/


\usepackage{mdwmath}
\usepackage{mdwtab}
% Also highly recommended is Mark Wooding's extremely powerful MDW tools,
% especially mdwmath.sty and mdwtab.sty which are used to format equations
% and tables, respectively. The MDWtools set is already installed on most
% LaTeX systems. The lastest version and documentation is available at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/mdwtools/


% IEEEtran contains the IEEEeqnarray family of commands that can be used to
% generate multiline equations as well as matrices, tables, etc., of high
% quality.


%\usepackage{eqparbox}
% Also of notable interest is Scott Pakin's eqparbox package for creating
% (automatically sized) equal width boxes - aka "natural width parboxes".
% Available at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/eqparbox/





% *** SUBFIGURE PACKAGES ***
%\usepackage[tight,footnotesize]{subfigure}
% subfigure.sty was written by Steven Douglas Cochran. This package makes it
% easy to put subfigures in your figures. e.g., "Figure 1a and 1b". For IEEE
% work, it is a good idea to load it with the tight package option to reduce
% the amount of white space around the subfigures. subfigure.sty is already
% installed on most LaTeX systems. The latest version and documentation can
% be obtained at:
% http://www.ctan.org/tex-archive/obsolete/macros/latex/contrib/subfigure/
% subfigure.sty has been superceeded by subfig.sty.



%\usepackage[caption=false]{caption}
%\usepackage[font=footnotesize]{subfig}
% subfig.sty, also written by Steven Douglas Cochran, is the modern
% replacement for subfigure.sty. However, subfig.sty requires and
% automatically loads Axel Sommerfeldt's caption.sty which will override
% IEEEtran.cls handling of captions and this will result in nonIEEE style
% figure/table captions. To prevent this problem, be sure and preload
% caption.sty with its "caption=false" package option. This is will preserve
% IEEEtran.cls handing of captions. Version 1.3 (2005/06/28) and later 
% (recommended due to many improvements over 1.2) of subfig.sty supports
% the caption=false option directly:
\usepackage[caption=false,font=footnotesize]{subfig}
%
% The latest version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/subfig/
% The latest version and documentation of caption.sty can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/caption/




% *** FLOAT PACKAGES ***
%
\usepackage{fixltx2e}
% fixltx2e, the successor to the earlier fix2col.sty, was written by
% Frank Mittelbach and David Carlisle. This package corrects a few problems
% in the LaTeX2e kernel, the most notable of which is that in current
% LaTeX2e releases, the ordering of single and double column floats is not
% guaranteed to be preserved. Thus, an unpatched LaTeX2e can allow a
% single column figure to be placed prior to an earlier double column
% figure. The latest version and documentation can be found at:
% http://www.ctan.org/tex-archive/macros/latex/base/



%\usepackage{stfloats}
% stfloats.sty was written by Sigitas Tolusis. This package gives LaTeX2e
% the ability to do double column floats at the bottom of the page as well
% as the top. (e.g., "\begin{figure*}[!b]" is not normally possible in
% LaTeX2e). It also provides a command:
%\fnbelowfloat
% to enable the placement of footnotes below bottom floats (the standard
% LaTeX2e kernel puts them above bottom floats). This is an invasive package
% which rewrites many portions of the LaTeX2e float routines. It may not work
% with other packages that modify the LaTeX2e float routines. The latest
% version and documentation can be obtained at:
% http://www.ctan.org/tex-archive/macros/latex/contrib/sttools/
% Documentation is contained in the stfloats.sty comments as well as in the
% presfull.pdf file. Do not use the stfloats baselinefloat ability as IEEE
% does not allow \baselineskip to stretch. Authors submitting work to the
% IEEE should note that IEEE rarely uses double column equations and
% that authors should try to avoid such use. Do not be tempted to use the
% cuted.sty or midfloat.sty packages (also by Sigitas Tolusis) as IEEE does
% not format its papers in such ways.



% correct bad hyphenation here
% \hyphenation{op-tical net-works semi-conduc-tor}


\begin{document}
%
% paper title
% can use linebreaks \\ within to get better formatting as desired
\title{GPU implementation of a region based algorithm \\ for large images segmentation}


% author names and affiliations
% use a multiple column layout for up to two different
% affiliations

\author{
\IEEEauthorblockN{Gilles Perrot, St\'{e}phane Domas, Rapha\"{e}l Couturier}
\IEEEauthorblockA{Distributed Numerical Algorithmics team (AND), Laboratoire d'Informatique de Franche-comt\'{e}\\
Rue Engel Gros, 90000 Belfort, France\\
forename.name@univ-fcomte.fr}
}



% use for special paper notices
%\IEEEspecialpapernotice{(Invited Paper)}




% make the title area
\maketitle

\begin{abstract}
Image segmentation is one of the most challenging issues in image computing.
In this work, we focus on region-based active contour techniques (snakes) as they seem to achieve a high level of robustness and fit with a large range of
applications. Some algorithmic optimizations provide significant speedups, but even so, execution times are still non-neglectable  
with the continuing increase of image sizes. Moreover, these algorithms are not well suited for running on multi-core CPU's.
At the same time, recent developments of Graphical Processing Units (GPU) suggest that higher speedups could be obtained 
by use of their specific design. We have managed to adapt a specially efficient snake algorithm that fits recent Nvidia GPU architecture 
and takes advantage of its massive multithreaded execution capabilities. The speedup obtained is most often around 7.
\end{abstract}

\begin{IEEEkeywords}
 GPU; segmentation; snake;
\end{IEEEkeywords}

\section{Introduction}
Segmentation and shape detection are still key issues in image computing. These techniques are used in numerous fields ranging from medical imaging to video tracking, shape recognition or localization.
Since 1988, the active contours (snakes) introduced by and  Kass et al. \cite{KassWT88}, have proved to be efficient and robust, especially against noise, for a wide range of image types. 

The main shortcoming of these algorithms is often their high dependence on the initial shape, though several contributions have lowered this dependency and also brought 
more accurate segmentation of non convex shapes \cite{Ruch01} \cite{XuP98}. 

The information that drives a snake model comes either from the contour itself or from the characteristics of the regions it defines. 
For noisy images, the second option is often more suitable as it takes into account the statistical fluctuations of the pixels. 
One approach \cite{ChesnaudRB99,AllainBG08} proposes a geometric (polygonal) region-based snake driven by the minimization of the stochastic complexity. One significant
advantage is that it runs without any free parameter which can be helpful when dealing with image sequences or slices (3D).

An important issue of image processing, especially segmentation, has always been the computation time of most algorithms. Over the years, the increase of CPU computing capabilities, 
although quite impressive, has not been able to fulfill the combined needs of growing resolution and real-time computation.
Since having been introduced in the early 1980's, the capabilities and speed of graphics accelerators have always been increasing. So much so that the recent GPGPU 
(General Purpose Graphic Processing Units) currently benefit by a massively parallel architecture for general purpose programming, especially when dealing with large matrices 
or vectors. On the other hand, their specific design obviously imposes a number of limitations and constraints. 
Some implementations of parametric snakes have already been tested,  such as \cite{Brunett}. However, a similar solution (computation per small tile)
is not suited for the algorithm we have implemented.

Our goal, in collaboration with the PhyTI team\footnote{Physics and Image Processing Group, Fresnel Institute, Ecole Centrale de Marseille (France)}, was to propose a way to fit their algorithm 
to the Nvidia$^{\textcopyright}$ Tesla GPU architecture.
The remainder of this paper presents the principles of the algorithm and notations in section \ref{secCPUalgooutlines}. In section \ref{secCPUalgodetails}, the details of
the sequential CPU implementation are explained. Section \ref{GPUgeneralites} summarizes Nvidia's GPU 
important characteristics and how to deal with them efficiently. Then sections \ref{GPUimplementation} and \ref{secSpeedups} detail our GPU implementation and timing results. 
Finally, the conclusion of section \ref{secConclusion} evaluates the pros and cons of this implementation and then gives a few direction to be followed in future works.
 


\section{\label{secCPUalgooutlines}Sequential algorithm : outlines} 
The goal of the active contour segmentation (snake) method we studied \cite{Ruch01} is to distinguish, inside an image $I$, a target region $T$ from the background region
$B$. The size of $I$ is L x H pixels of coordinates $(i,j)$ and gray level $z(i,j)$.
We assume that the gray levels of $T$ and $B$ are independent random vectors, each with a distribution $p^{\Omega}$ of its components $(\Omega \in \{T ; B\})$. 
The present implementation uses a Gaussian distribution, but another one can easily be used as Gamma, Poisson,...(Cf. \cite{ChesnaudRB99})\dots 

The \textit{active contour} $S$, which represents the shape of $T$ is chosen as polygonal.
The purpose of the segmentation is then to determine the shape that optimizes a pseudo log-likelihood-based criterion (PLH).
This is done by a very simple iterative process which is initialized with an arbitrary shape, then at each step :
\begin{enumerate}
  \item it modifies the shape
  \item it estimates the parameters of the Gaussian functions for the two regions and evaluates the criterion.
  \item it validates the new shape if the criterion has a better value.
\end{enumerate}
A simplified description of it is given in \emph{Algorithm \ref{cpualgosimple}} which features two nested loops : the main one, on iteration level, is 
responsible for tuning the number of nodes ; the inner one, on step level, takes care of finding the best shape for a given number of nodes.
\emph{Figure \ref{images_algo}} shows intermediate results at iteration level. Sub-figure \emph{\ref{fig:labelinit}} shows the initial rectangular shape, \emph{\ref{fig:labelit1}}
 shows the best four-node shape that ends
the first iteration. Sub-figures \emph{\ref{fig:labelit2}} and \emph{\ref{fig:labelit4}} show the best shape for an eight-node snake (resp. 29-node) 
which occurs at the end of the second  iteration (resp. fourth). 

\begin{algorithm}[h]
\label{cpualgosimple}
\caption{Sequential algorithm : outlines}
\SetNlSty{textbf}{}{:}

  %compute\_cumulated\_images()\;
  begin with a rectangular 4 nodes snake\;
  \Repeat(\tcc*[f]{iteration level}){no more node can be added}{
      \Repeat(\tcc*[f]{step level}){no more node can be moved}{
          Test some other positions for each node, near its current position\;
          Find the best PLH and adjust the node's position\;
      }
  Add a node in the middle of each \emph{long enough} segment\;
  }
\end{algorithm}


\begin{figure}[h]
 \centering
\subfloat[Initial snake ]{\label{fig:labelinit} \includegraphics[width=0.4\linewidth]{./img/cochon_petit_init.jpg}}\qquad
\subfloat[End of first iteration (4 nodes) ]{\label{fig:labelit1} \includegraphics[width=0.4\linewidth]{./img/cochon_petit_it1.jpg}}\\
\subfloat[End of second iteration (8 nodes)]{\label{fig:labelit2} \includegraphics[width=0.4\linewidth]{./img/cochon_petit_it2.jpg}}\qquad
\subfloat[End of fourth iteration (29 nodes)]{\label{fig:labelit4} \includegraphics[width=0.4\linewidth]{./img/cochon_petit_it4.jpg}}
%\subfloat[width=0.4\linewidth]{./img/cochon_b_entier.jpg}
 % cochon_b_entier.jpg: 3960x2970 pixel, 72dpi, 139.70x104.78 cm, bb=0 0 3960 2970
 \caption{segmentation of a noisy image}
  \label{images_algo}
\end{figure}



\section{\label{secCPUalgodetails}Sequential algorithm : details}
\subsection{Criterion}
For $p^{\Omega}$ a Gaussian function, $\Theta_{\Omega}$ ($\Omega \in \{T ; B\}$) has two components, the average value $\mu$ and the deviation $\sigma$ which are estimated by 
$$
\widehat{\Theta_{\Omega}} \left(
\begin{array}{l}
\widehat{\mu} = \frac{1}{N_{\Omega}} \displaystyle\sum_{(i,j)\in \Omega} z(i,j) \\
\widehat{\sigma^2} = \frac{1}{N_{\Omega}} \displaystyle\sum_{(i,j)\in \Omega} z^2(i,j) - \mu^2 \\
  \end{array}
\right.
$$
The likelihood of a region is given by
$$ P[I|S_{n,l}, \Theta_T, \Theta_B] = P(\chi_T | \Theta_T)P(\chi_B | \Theta_B)$$ 
where
$$P(\chi_{\Omega} | \Theta_{\Omega}) = \prod_{(i,j)\in \Omega} p^{\Omega}[z(i,j)] ~~~~(\Omega \in \{T ; B\})$$ 
And then the log-likelihood by 
$$-N_{\Omega}\log\left(\sqrt{2\pi}\right) -N_{\Omega}.log\left(\sigma\right) - \frac{1}{2\sigma^2}\sum_{(i,j)\in \Omega} \left( z(i,j)-\mu \right)^2 $$
Considering the two regions, the criterion to be optimized is then :
$$C = \frac{1}{2}\left( N_B\log\left(\widehat{\sigma_B}^2\right) + N_T\log\left(\widehat{\sigma_T}^2\right)\right)$$ 

\subsection{CPU implementation}
Let $S_{n,l}$ be the snake state at step $l$ of iteration $n$, and $S_{n,l}^i$ the node $i$ of $S_{n,l}$ ($i \in [0;N_n]$).
Each segment of $S_{n,l}$ is considered as an oriented list of discrete points.
Chesnaud \& Refregier \cite{ChesnaudRB99} have shown how to replace the 2 dimensions sums needed to estimate $\Theta_{\Omega}$ by 1 dimension sums along $S_{n,l}$.
However, this approach involves weighing coefficients for every single point of $S_{n,l}$ which leads to compute a pair of transformed images, at the very 
beginning of the process. Such images are called cumulated images and will be used as lookup tables.
Therefore, beyond this point, we will talk about the \emph{contribution} of each point to the 1D sums. By extension, we also talk about the \emph{contribution} of each segment to the 1D sums.

A more detailed description of the sequential algorithm is given by \emph{Algorithm \ref{cpualgo}}. 
The process starts with the computation of cumulated images ; an initialization stage  takes place from line \ref{debinit} to line \ref{fininit}. 
Then we recognize the two nested loops (line \ref{loopnewnodes} and line \ref{loopmovenodes}) and finally the heart of the algorithm stands on line \ref{kernelPLH} which represents
the main part of the calculations to be done :
\begin{enumerate}
 \item compute the various sums without the contributions of both segments connected to current node  $S_{n,l}^i$.
  \item \label{CPUcontrib_segments} compute the contributions of both segments, which requires :
  \begin{itemize}
   \item \label{CPUbresenham} To determine the coordinates of every discrete pixel of both segments connected to $S_{n,l}^{i,w}$.
    \item \label{CPUcontrib_pixels} To compute every pixel contribution.
    \item To sum pixel contributions to obtain segment contributions.
  \end{itemize}
   \item compute the PLH given the contribution of each segment of the tested snake.
\end{enumerate}

\begin{algorithm}[h]
\SetNlSty{textbf}{}{:}
\caption{Sequential simplified algorithm} 
\label{cpualgo}
   read image from HDD\;
   compute\_cumulated\_images()\label{cumuls}\;
   iteration $n \leftarrow 0$\label{debinit}\;
   $N_0 \leftarrow 4$\;
   $S_{n,l} \leftarrow S_{0,0}$\;
   step $d \leftarrow d_{max} = 2^q$\;
   current node $S_{0,0}^i \leftarrow S_{0,0}^0$\;
   $l \leftarrow 0$\;
   compute $PLH_{ref}$, the PLH of $S_{n,0}$\label{fininit}\;
  \Repeat(\tcc*[f]{iteration level}){no new node added}{\label{loopnewnodes}
    \Repeat(\tcc*[f]{step level}){no node move occured}{\label{loopmovenodes}
      \For{$i=0$ to $N_n$}{
        $S_{n,l}^{i,w}$ ($w \in [0;7]$) are the neighbors of $S_{n,l}^i$ by distance $d$\;
        \For{$w=0$ to $7$}{
          compute $PLH_w$ for $S_{n,l}$ when $S_{n,l}^{i,w}$ replaces $S_{n,l}^i$ \label{kernelPLH}\;
          \lIf{$PLH_w$ is better than $PLH_{ref}$}{
            $PLH_{ref} \leftarrow PLH_w$\;
            move node $S_{n,l}^i \leftarrow S_{n,l}^{i,w}$\;
          }
        }
      }
     $l \leftarrow l+1$\;
   }
    add new nodes, $N_n \leftarrow N_n + N_{newnodes}$\;
   \lIf{$d > 1$}{ $d \leftarrow d/2$ } \lElse{ $d=1$ }\;
    $n \leftarrow n+1$\;
    compute $PLH_{ref}$, the PLH of $S_{n,0}$ \;
  }
\end{algorithm}



The profiling results of the CPU implementation shown in \emph{Figure \ref{CPUprofile}} display the relative costs of the most time-consumming functions.
It appears that more than 80\% of the total execution time is always spent by only three functions~:
\begin{itemize}
 \item \texttt{compute\_segment\_contribution()} which is responsible for point \ref{CPUcontrib_segments} above,
  \item \texttt{compute\_cumulated\_images()} which computes the 3 lookup tables at the very beginning,
  \item \texttt{compute\_pixels\_coordinate()} which is called by \texttt{compute\_segment\_contribution()}.
\end{itemize}

\begin{figure}[h]
 \centering
 \includegraphics[width=0.9\linewidth, height=0.5\linewidth]{./img/data_profile_cpu.png}
 \caption{\label{CPUprofile}the three most-consumming functions for various image sizes}
\end{figure}

Measurements have been performed for several image sizes from 15~MPixels (about 3900 x 3900) 
to 144 MPixels (about 12000 x 12000). On the one hand, we can notice that function \texttt{compute\_segment\_contribution()} always lasts more than 45\% of the total running time, and even 
more when the image gets larger. 
On the other hand, the function \texttt{compute\_cumulated\_images()} costs more than 23\%, decreasing with image size, while function \texttt{compute\_pixels\_coordinate()} always takes around 6\%.
It confirms that the need for parallelization resides in line \ref{kernelPLH} and line \ref{cumuls} of Algorithm \ref{cpualgo} as they contain every call to those three functions.

The following sections detail how we managed to implement these time-consumming functions in parallel, but
a brief reminder on GPU's recent architecture is presented first. 
 


\section{\label{GPUgeneralites}NVidia's GPU architecture}
GPUs are multi-core, multi-threaded processors, optimized for highly parallel computation. Their design focuses on SIMT model by devoting
more transistors to data processing rather than data-caching and flow control \cite{CUDAPG}.
 
For example, Figure \ref{GPUC1060} shows a Tesla C1060  with its 4GB of global memory and 30 SM processors, each including  :
\begin{itemize}
 \item 8 Scalar Processors (SP)
  \item a Floating Point Unit (FPU)
  \item a parallel execution unit (SIMT) that runs threads by warps of 32. 
 \item 16KB of shared memory, organized in 16 banks of 32 bits words
\end{itemize}
Nvidia uses a parameter called the \emph{compute capability} of each GPU model. Its value is composed of a major number and a minor number ; for example the C1060 is a sm13 GPU (major=1 minor=3) 
and C2050 is a sm20 GPU.

\begin{figure*}[htbp]
 \centering
 \includegraphics[width=0.7\linewidth]{./img/GPU_block.png} 
 \caption{\label{GPUC1060}schematic diagram of GPU's internal architecture}
\end{figure*}

The recent Fermi cards (eg. C2050,) have improved performances by supplying more shared memory in a 32 banks array, a second execution 
unit and several managing 
capabilities on both the shared memory and level 1 cache memory ( \cite{CUDAPG}, \cite{CUDAFT}, \cite{CUDAFC}. 
However, writing efficient code for such architectures is not obvious, as re-serialization must be avoided as much as possible. Thus, when designing, one must 
keep a few key points in mind :
\begin{itemize}
 \item CUDA model organizes threads by a) threads blocks in which synchronization is possible, b) a grid of blocks with no possible synchronization
      between blocks.
  \item there is no way to know in what order the blocks are to be scheduled during one single kernel execution.
  \item data must be kept in GPU memory, to reduce the overhead due to copying between CPU and GPU memories.
  \item the total amount of threads running the same computation must be maximized.
  \item the number of execution branches inside a block should be reduced as much as possible.
  \item global memory accesses should be coalescent, \emph{ie}. memory accesses done by physically parallel threads (16 at a time) must be consecutive and contained in a 128 Bytes range.
  \item shared memory is organized by 16 x 32 bits wide banks. To avoid bank conflicts, each parallel thread (16 at a time) must access a different bank.
\end{itemize}

All the above charasteristics make it always a quite constrained problem to solve when designing a GPU code. 
%detailler
Moreover, a non suited code would probably run even slower on GPU than on CPU due to the automatic serialization which would be done at run time.
%%
\section{\label{GPUimplementation}GPU implementation}
In the implementation described below, pre-computations and proper segmentation are discussed separately.
To keep data in GPU memory, the whole computation is assigned to the GPU. CPU still hosts :
\begin{itemize}
 \item data reading from HDD
  \item data writing on HDD if needed
  \item main loops control (corresponding to lines \ref{loopnewnodes} and \ref{loopmovenodes} of Algorithm \ref{cpualgo})
\end{itemize}

It must be noticed that controlling these loops is achieved with only a very small amount of data being transferred between host (CPU) and device (GPU),
which does not produce high overhead. \\
Morever, the structures described below need 20 Bytes per pixel of the image to process (plus an offset of about 50~MByte). 
It defines the maximum image size we can accept : approximately 150 M Pixels.

\subsection{Pre-computations}
To replace 2D sums by 1D sums, Chesnaud \textit{et al.} \cite{ChesnaudRB99}  have shown that the three matrices below should be computed :
$$C_1(i,j) = \sum_{k=0}^{k=j} (1+k)$$
$$C_z(i,j) = \sum_{k=0}^{k=j} z(i,k)$$ and
$$C_{z^2}(i,j) = \sum_{k=0}^{k=j} z^2(i,k)$$
Where $z(i,k)$ is the gray level of pixel of coordinate $(i,j)$, so that $C_1$, $C_z$ and $C_{z^2}$ are the same size as image $I$.

\begin{figure*}[htbp]
  \centering
 \resizebox{0.8\linewidth}{0.3\linewidth}{\input{./img/GPUcumuls.pdf_t}}
  \caption{\label{GPUcumuls}\texttt{compute\_blocks\_prefixes()} details.}
\end{figure*}

\medskip 
First, we chose not to generate $C_1(i,j)$, which requires that values should be computed when needed, but saves global memory and does not lead to any overhead.
The computation of $C_{z}$ and $C_{z^2}$ easily decomposes into series of  \emph{inclusive prefixsums} \cite{Harris07}.
However, by keeping the \emph{1 thread per pixel} rule, as the total number of threads that can be run in a grid cannot exceed $2^{25}$ (Cf. \cite{CUDAPG}),
slicing is necessary for images exceeding a size threshold which can vary according to the GPU model (e.g. 33 MPix for sm13 GPU, eg. C1060).
It's quite easy to do, but it leads to a small overhead as the process requires multiple calls to one kernel.
Slicing can be done in two ways :
\begin{itemize}
 \item all slices are of the same size (balanced)
  \item slices fit the maximum size allowed by the GPU, leaving one smaller slice at the end of the process (full-sized).
\end{itemize}
The balanced slice option has proved to run faster.\\
For example : if a given image has 9000 lines and the GPU can process up to 4000 lines at a time, it's faster to run 3 times with 3000 lines rather than twice with 
4000 and once with 1000. 

As the sums in $C_z$ and $C_{z^2}$ are row-wide, it is easy to see that every block-wide sum will be needed before being able to use it in the global sum.
But as mentioned earlier, the scheduling of blocks must be considered as random.
So, in order to ensure synchronizations, each row of the original image is then treated by three different kernels : 
\begin{itemize}
 \item \texttt{compute\_blocks\_prefixes()}.
  \item \texttt{scan\_blocksums()}.
  \item \texttt{add\_sums2prefixes()}.   
\end{itemize}
Figures \ref{GPUcumuls}, \ref{GPUscansomblocs} and \ref{GPUaddsoms2cumuls} show relevant data structures for a given row $i$ of $I$. 
We assume that each thread block runs $bs$ threads in parallel and each row of $C_z$ needs $n$ blocks to cover its $L$ pixels.

Figure \ref{GPUcumuls} shows the details of the process for row $i$ of the original image $I$, already stored in GPU global memory.
Operands are first copied into GPU shared memory for efficiency reasons.
An inclusive prefixsum is then performed inside each independant thread block. 
At this point, only the first shared memory block contains the final values. Its last element contains the sum of all
elements in the corresponding block of $I$.
In order to obtain the right values for the row $i$ of $C_z$, every element value in the other blocks must then be summed with an offset value. 
This offset value is the sum of all element values in every corresponding previous block of row $i$.

As the scheduling of blocks is fully unpredictable, the necessary intermediate results have to be stored in GPU global memory before exiting from kernel.
Each element of the prefixsums in GPU shared memory has been stored in its corresponding position in $C_z$ (GPU global mem),
along with the vector of block sums which will be passed later to the next kernel \texttt{scan\_blocksums()}.

The kernel \texttt{scan\_blocksums()} (Figure \ref{GPUscansomblocs}) only makes an exclusive prefixsum on the vector of block sums described above.
The result is a vector containing, at index $x$, the value to be added to every element of block $x$ in each line of $C_z$.

This summing is done in shared memory by kernel \texttt{add\_sums2prefixes()} as described by Figure \ref{GPUaddsoms2cumuls}.

The values of $C_{z^2}$ are obtained together with those of $C_{z}$ and in exactly the same way.
For publishing reasons, figures do not show the $C_{z^2}$ part of structures. 



\begin{figure*}[htbp]
  \centering
 \resizebox{0.6\linewidth}{0.2\linewidth}{\input{./img/GPUscansomblocs.pdf_t}}
  \caption{\label{GPUscansomblocs}\texttt{scan\_blocksums()} details.}
\end{figure*}

\begin{figure*}[htbp]
  \centering
 \resizebox{0.7\linewidth}{0.4\linewidth}{\input{./img/GPUaddsoms2cumuls.pdf_t}}
  \caption{\label{GPUaddsoms2cumuls}\texttt{add\_sums2prefixes()} details.}
\end{figure*}

With this implementation, speedups are quite significant (Table \ref{tabresults}). Moreover, the larger the image, 
the higher the speedup is, as the step-complexity of the sequential algorithm is of $O(N^2)$ and  $O(N\log(N))$ for the parallel version.
Even higher speedups are achieved by adapting the code to specific-size images, especially when the number of columns is a power of 2. This avoids 
inactive threads in the grid, and thus improves efficiency.
However, on sm13 GPUs, these computations are made with a 2-way bank conflict as sums are based on 64-bit words, thus creating overhead.


\subsection{Segment contributions}
The choice made for this implementation has been to keep the \emph{1 thread per pixel} rule for the main kernels.
Of course, some reduction stages need to override this principle and will be pointed out.

As each of the $N_n$ nodes of the snake $S_{n,l}$ may move to one of the eight neighbor positions as shown in \emph{Figure \ref{GPUtopo}}, 
there is $16 N_n$ segments whose contribution has to be estimated.
The best combination is then chosen to  obtain $S_{n,l+1}$ (Figure \ref{GPUtopo}).
Segment contributions are computed in parallel by kernel \texttt{GPU\_compute\_segments\_contrib()}.

\begin{figure}[h]
  \centering
 \resizebox{0.9\linewidth}{0.81\linewidth}{\input{./img/topologie.pdf_t}}
  \caption{\label{GPUtopo}topology around nodes}
\end{figure}

The grid parameters for this kernel are determined according to the size of the longest segment $npix_{max}$.
If $bs_{max}$ is the maximum theoritical blocksize that a GPU can accept,
\begin{itemize}
 \item the block size $bs$ is taken as 
  \begin{itemize}
   \item $npix_{max}$'s next power of two if \\${npix_{max} \in [33 ; bs_{max} ] }$
    \item 32 if ${npix_{max} < 32 }$
    \item $bs_{max}$ if ${npix_{max} > 256 }$
  \end{itemize}
  \item the number of threads blocks assigned to each segment, $N_{TB} = \frac{npix_{max} + bs -1 }{bs}$ 
\end{itemize}
Our implementation makes intensive use of shared memory and does not allow the use of the maximum theoritical blocksizes 
(512 for sm13, 1024 for sm20, see \cite{CUDAFT} and \cite{CUDAPG}).
Instead we set $bs_{max}^{sm13} = 256$ and $bs_{max}^{sm20} = 512$. 
Anyway, testing has shown that most often, the best value is 256 for both \textit{sm13} and \textit{sm20} GPU's.

\begin{figure*}[htbp]
  \centering
 \resizebox{0.6\linewidth}{0.35\linewidth}{\input{./img/contribs_segments.pdf_t}}
  \caption{\label{contribs_segments}structure for segments contributions computation. Gray symbols help to locate inactive threads as opposed to black 
          ones that figure active threads.}
\end{figure*}

Then \texttt{GPU\_compute\_segments\_contrib()} computes in parallel :
\begin{itemize}
 \item every pixel coordinates for all $16 N_n$ segments. Since the snake is only read in one direction, we have been able 
  to use a very simple parallel algorithm instead of Bresenham's.
It is based on the slope $k$ of each segment~: one pixel per row if $|k|>1$, one pixel per column otherwise.
  \item every pixel contribution by reading the corresponding values in the lookup tables.
  \item every thread blocks sums of individual pixel contributions by running a \textit{reduction} stage for each block.
\end{itemize}
The top line of Figure \ref{contribs_segments} shows the base data structure in GPU shared memory which is relative to one segment.
We concatenate the single segment structure as much as necessary to create a large vector representing every pixel of every test segment.
As each segment has a different size (most often different from any power of two), there is a non-neglectable number of inactive threads scattered in the whole structure.
Two stages are processed separately : one for all even nodes and another one for odd nodes,
as shown in the two bottom lines of Figure \ref{contribs_segments}. 


The process is entirely done in shared memory ; only a small amount of data needs to be stored in global memory for each segment~:
\begin{itemize}
 \item the coordinates of its middle point, in order to be able to add nodes easily if needed.
  \item the coordinates of its first and last two points, to compute the slope at each end of the segment.
\end{itemize}
The five values above are part of the weighing coefficients determination for each segment and node.

 The \texttt{GPU\_sum\_contribs()} takes the blocks sums obtained by \texttt{GPU\_compute\_segments\_contrib()} and computes a second stage parallel summing to provide 
the $16 N_n$ segment contributions.

\subsection{Segments with a slope $k$ such as $|k|\leq1$}
Such a segment is treated with 1 thread per column and consequently, it often has more than one pixel per row as shown by Figure \ref{tripix}.
In an image row, consecutive pixels which belong to the target define an interval which can only have one low and one high ends 
That's why, on each row, we choose to consider only the contributions of the innermost pixels.
This selection is also done inside \texttt{GPU\_compute\_segments\_contrib()} when reading the lookup tables for each pixel contribution. 
We simply set a null contribution for pixels that need to be ignored.
\begin{figure}[h]
  \centering
 \resizebox{0.75\linewidth}{0.35\linewidth}{\input{./img/tripix.pdf_t}}
  \caption{\label{tripix}Zoom on part a of segment with $|k| < 1$,  at pixel level.}
\end{figure}


\subsection{Parameters estimation}
A \texttt{GPU\_compute\_PLH()} kernel computes in parallel :
\begin{itemize}
 \item every $8N_n$ vector of parameters values corresponding to each possible next state of the snake. Summing is done in shared memory but relevant 
  data for these operations are stored in global memory.
  \item every associated pseudo likelihood value.
  \item node substitutions when better PLH have been found and if it does not lead to segments crossing.
\end{itemize}

\subsection{End of segmentation}
Segmentation is considered achieved out when no other node can be added to the snake (Algorithm \ref{gpualgosimple}).
A very simple GPU kernel adds every possible node and returns the number it added.

\begin{algorithm}[h]
\label{gpualgosimple}
\caption{Parralel GPU algorithm : outlines. \texttt{<<<...>>>} indicates a GPU kernel parallel process.}
\SetNlSty{textbf}{}{:}
  load images\;
  transfer image from CPU to GPU\;
  \texttt{<<<}compute the 2 cumulated images\texttt{>>>}\;
  \texttt{<<<}initialize the snake\texttt{>>>}\;
  \Repeat(\tcc*[f]{iteration level}){no more node can be added}{
      \Repeat(\tcc*[f]{step level}){no more node can be moved}{
          \texttt{<<<}find best neighbor snake\texttt{>>>}\;
          \texttt{<<<}adjust node's positions\texttt{>>>}\;
          transfer the number of moves achieved from GPU memory to CPU memory.
      }
  \texttt{<<<}Add nodes\texttt{>>>}\;
  transfert the number of nodes added from GPU memory to CPU memory.
  }
\end{algorithm}

\section{\label{secSpeedups}Speedups}
The CPU (SSE) implementation by N. Bertaux from the PhyTI team, based on \cite{AllainBG08} has been our reference to ensure segmentation's quality and to estimate speedups. 
Results are given in Table \ref{tabresults}. 
CPU timings were measured on an Intel Xeon E5530-2.4GHz with 12Go RAM (LIFC cluster).
GPU timings were obtained on a C2050 GPU with 3GB RAM (adonis-11.grenoble.grid5000.fr).\\
Execution times reported are means on ten executions.
%Measurements on CPU may vary more than on GPU. So CPU results given in \ref{tabresults} are near the fastest values we observed.
The image of figure \ref{fig:labelinit} (scaled down for printing reasons) is a 16-bit gray level photo from PhyTI team,
voluntarily noisy for testing reasons. The contrast has been enhanced for better viewing.

We separately give the timings of pre-computations as they are a very general purpose piece of code.
Segmentations have been performed with strictly the same parameters (initial shape, threshold length).
The neighborhood distance for the first iteration is 32 pixels. It has a slight influence on the 
time process, but it leads to similar speedups values of approximately 7 times faster than CPU.

Though it does not appear in Table \ref{tabresults}, we observed that during segmentation stage, higher speedups are obtained in the very first iterations, when segments are made of a lot of pixels, leading to a higher parallelism ratio.\\
Several parameters prevent from achieving higher speedups~:
\begin{itemize}
 \item accesses in the lookup tables in global memory cannot be coalescent. It would imply that the pixel contributions of a segment are stored in consecutive spaces in $C_z$ and $C_{z^2}$. 
   This is only the case for horizontal segments.
  \item the use of 64-bit words for computations in shared memory often leads to 2-way bank conflicts.
  \item the level of parallelism is not so high, ie. the total number of pixel is not large enough to achieve impressive speedups. For example, on C2050 GPU, a grid can 
  run about 66 million of threads, but a snake in a 10000 x 10000 image would be less than 0.1 million pixel long.
\end{itemize}


% \begin{table}
%  \begin{center}
% \begin{tabular}{|l| r|r r r|}
% \hline
% && CPU & GPU & Speedup\\\cline{3-5}
% Image 15MP    & \bf total     & \bf0.51 s     &  \bf0.06 s & \bf x8.5 \\
%               & pre-comp.     & 0.21 s        &  0.02 s & x10\\
%               & segment.      & 0.34 s        &  0.04 s & x8.5\\\hline
% Image 100MP   & \bf total     & \bf 4.33 s    &  \bf 0.59 s & \bf x7.3\\
%               & pre-comp.     & 1.49 s        &  0.13 s & x11\\
%               & segment.      & 2.84 s        &  0.46 s & x6.1\\\hline
% Image 150Mp   & \bf total     & \bf 26.4 s    &  \bf 0.79 s & \bf x33\\
%               & pre-comp.     & 8.4 s         &  0.20 s & x42\\
%               & segment.      & 18.0 s        &  0.59 s & x30\\\hline
%  \end{tabular}
%  \end{center}
% 
% \caption{\label{tabresults} GPU (C2050, sm20) vs CPU timings.}
% \end{table} 


\begin{table}
 \begin{center}
\begin{tabular}{|l| r|r r r|}
\hline
&& CPU & GPU & Speedup\\\cline{3-5}
Image 15MP      & \bf total     & \bf0.51 s     &  \bf0.06 s & \bf x8.5 \\
                & pre-comp.     & 0.13 s        &  0.02 s & x6.5\\
                & segment.      & 0.46 s        &  0.04 s & x11.5\\\hline
Image 100MP     & \bf total     & \bf 4.08 s    &  \bf 0.59 s & \bf x6.9\\
                & pre-comp.     & 0.91 s        &  0.13 s & x6.9\\
                & segment.      & 3.17 s        &  0.46 s & x6.9\\\hline
Image 150Mp     & \bf total     & \bf 5.7 s     &  \bf 0.79 s & \bf x7.2\\
                & pre-comp.     & 1.4 s         &  0.20 s & x7.0\\
                & segment.      & 4.3 s         &  0.59 s & x7.3\\\hline
 \end{tabular}
 \end{center}

\caption{\label{tabresults} GPU (C2050, sm20) vs CPU timings.}
\end{table} 

\IEEEpeerreviewmaketitle

\vspace{1cm}

\section{\label{secConclusion}Conclusion}
The algorithm we have focused on is not easy to adapt for high speedups on GPGPU, though we managed to make it work faster than on CPU.
The main drawback is clearly its relative low level of parallelism. Nevertheless, we proposed different kernels that allowed us to take advantage of the computation power of  GPUs.
 In future works, we plan  to try and manage to benefit from larger computing grids of thread blocks. Among the possible solutions, we plan to work on:
\begin{itemize}
  \item slicing the image and proceeding the parts in parallel. This is made possible since sm20 GPU provide multi kernel capabilities.
  \item slicing the image and proceeding the parts on two different GPUs, hosted by the same CPU. 
   \item translating the parallelism from pixel level (\emph{1 thread per pixel}) to snake level (\emph{1 thread per snake}), at least during the first iteration, which 
is often the longest lasting one. 
  \item designing an algorithm, in a GPU way of thinking, instead of adapting the existing CPU-designed algorithm to GPU constraints as we did.
\end{itemize}


%%RAPH
%%Est ce qu'on parle du fait qu'on va également réfléchir à repenser l'algo en gpu?


% trigger a \newpage just before the given reference
% number - used to balance the columns on the last page
% adjust value as needed - may need to be readjusted if
% the document is modified later
%\IEEEtriggeratref{8}
% The "triggered" command can be changed if desired:
%\IEEEtriggercmd{\enlargethispage{-5in}}

% references section


\bibliographystyle{IEEEtran}

\bibliography{IEEEabrv,biblio}


% that's all folks
\end{document}