7 avril pour jury

[these_gilles.git] / THESE / Chapters / chapter3 / chapter3.tex~

\section{Présentation de l'algorithme}
La principale difficulté soulevée par l'emploi d'algorithmes de type \textit{snake} orientés contour est le choix de la fonction d'énergie externe et la détermination de la nature des images auxquelles elle convient. 
Dans l'approche orientée régions, les deux régions que sont l'extérieur et l'intérieur du contour (cas mono cible) sont prises en compte dans l'estimation de la forme du contour ;  cela permet d'extraire des formes dans des images où les contours de la cible sont mal définis, en raison d'un fort niveau de bruit par exemple.
Les algorithmes découlant de cette approche n'ont fait l'objet, à notre connaissance, d'aucune parallèlisation sur GPU, malgré le grand intérêt qu'elles revêtent dans l'interprétation d'images fortement bruitées ( RADAR, médicales,\dots ) et le besoin d'en réduire les temps d'exécution suffisamment  pour permettre l'interactivité. 
Nous proposons dans la suite de ce chapitre de détailler tout d'abord l'algorithme séquentiel que nous avons pris comme référence, puis d'en présenter la version parallèle pour GPU que nous en avons conçu.
L'algorithme a été décrit et proposé initialement en 1999 par Chesnaud \textit{et al.} dans \cite{ChesnaudRB99}. L'implémentation que les auteurs ont développé a continué d'être améliorée jusqu'à aujourd'hui et est employée comme brique élémentaire dans des algorithmes plus complexes. La version qui sert de référence ici est une implémentation séquentielle optimisée qui met aussi à profit les capacités de parallélisme des CPU actuels en employant le jeu d'instruction SSE2 des microprocesseurs. La description que nous en faisons dans les lignes qui suivent est très largement inspirée de \cite{ChesnaudRB99} à la différence que nous n'implémentons pas les critères de régularisation du contour ni de minimisation de la longueur de description pour nous focaliser sur la déformation du contour et sa convergence. 

\subsection{Formulation}
À l'intérieur de l'image observée $\bar{v}$, soient $\bar{t}$ le vecteur composé par les niveaux de gris des $N_t$ pixels de la région cible $\Omega_t$ et $\bar{b}$ celui des $N_b$ pixels du fond $\Omega_b$. Les vecteurs $\bar{t}$ et $\bar{b}$ sont supposés non corrélés et sont caractérisés par leurs densités de probabilité (PDF) respectives $p^{\Theta_t}$ et $p^{\Theta_b}$ ; $\Theta_t$ et $\Theta_b$ étant les vecteurs des paramètres de leurs PDF. Dans le cas gaussien que nous supposerons ici, $\Theta = (\mu, \sigma)$ où $\mu$ est la moyenne et $\sigma^2$ est la variance.
On note $\Gamma$ le contour de la région cible ($\Gamma \in \Omega_t$), que l'on suppose continu en connexité à 8 voisins. 

Le but de la segmentation est alors de déterminer la géométrie de $\Gamma$ qui maximise un critère de vraisemblance généralisée (GL).
La vraisemblance sur l'ensemble de l'image, \textit{ie.} la région $\Omega$ est donnée par

\begin{equation}
P\left(\bar{v} | \Omega_t, \Omega_b, \Theta_t, \Theta_b\right) = P\left(\bar{v} | \Omega_t, \Theta_t\right)P\left(\bar{v}|\Omega_b, \Theta_b\right)   
\end{equation}
soit en développant 
\begin{equation}
P\left(\bar{v} | \Omega_t, \Omega_b, \Theta_t, \Theta_b\right) = \prod_{x_k\in\Omega_t}p^{\Theta_t}\left(v_k, \Theta_t\right)\prod_{x_k\in\Omega_b}p^{\Theta_b}\left(v_k, \Theta_b\right)  
\label{eq-lhprod1}
\end{equation}

Dans le cas gaussien, la PDF étant de la forme 
\begin{equation}
p(\alpha) = \frac{1}{\sigma\sqrt{2\pi}}\mathrm{e}^{-\frac{\left(\alpha - \mu\right)^2}{2\sigma^2}}  
\label{eq-pdfgauss}
\end{equation}

La substitution de \eqref{eq-pdfgauss} dans \eqref{eq-lhprod1}, suivie du logarithme, permet d'obtenir l'expression de la \textit{log-vraisemblance}
\begin{equation}
  \label{eq-gl1}
  -N_tln\left(\sqrt{2\pi}\right)-N_tln\left(\sigma_t\right)-\frac{1}{2\sigma_t^2}\sum_{x_k\in\Omega_t}\left(v_k-\mu_t\right)^2
-N_bln\left(\sqrt{2\pi}\right)-N_bln\left(\sigma_b\right)-\frac{1}{2\sigma_b^2}\sum_{x_k\in\Omega_b}\left(v_k-\mu_b\right)^2
\end{equation}
dans laquelle les vecteurs $\Theta_t$ et $\Theta_b$ sont estimés suivant la méthode du  maximum de vraisemblance, qui donne l'expression générique suivante pour l'estimée de $\Theta_t$, notée $\widehat{\Theta_t}$ et transposable à l'identique pour $\Theta_b$ 
\begin{equation}
  \label{eq-teta}
  \widehat{\Theta_t} \left(
\begin{array}{l}
\widehat{\mu_t} = \frac{1}{N_t} \displaystyle\sum_{x_k\in \Omega_t} v_k \\
\widehat{\sigma^2_t} = \frac{1}{N_t} \displaystyle\sum_{x_k\in \Omega_t} \left(v_k - \widehat{\mu_t}\right)^2 \\
  \end{array}
\right.
\end{equation}

En intégrant \eqref{eq-teta} dans \eqref{eq-gl1}, il reste, à une constante près, le critère de vraisemblance généralisée suivant, noté GL, que l'on cherche à  optimiser en déterminant la géométrie de $\Gamma$ qui en maximise la valeur et épousera alors au mieux la forme du contour de la cible.
\begin{equation}
  \label{eq-gl}
  GL=\frac{1}{2}\left(N_tln\left(\widehat{\sigma_t^2}\right)+N_bln\left(\widehat{\sigma_b^2}\right)\right)
\end{equation}

\subsection{Optimisation des calculs}
La maximisation de GL est effectuée en employant une technique itérative où sa valeur doit être calculée à chaque déformation du contour $\Gamma$.
Si l'on se reporte à l'équation \eqref{eq-teta}, on voit que l'obtention de la valeur de GL nécessite, à chaque évaluation d'une géométrie donnée de $\Gamma$, le calcul des sommes  
\begin{align}
\label{eq-sommes1}
  S_v(\Omega_t) &= \sum_{x_k\in \Omega_t} v_k & S_{v^2}(\Omega_t) &= \displaystyle\sum_{x_k\in \Omega_t} v_k^2 \nonumber\\
  S_v(\Omega_b) &= \sum_{x_k\in \Omega_b} v_k & S_{v^2}(\Omega_b) &= \displaystyle\sum_{x_k\in \Omega_b} v_k^2 \nonumber\\
\end{align}

Considèrons la région cible $\Omega_t$, les pixels de coordonées $(i,j)$ qui la composent, et généralisons l'écriture des sommes de \eqref{eq-sommes1} en 
\begin{equation}
  \label{eq-sommes-gene}
  S_f(\Omega_t) = \sum_{i=i_{min}}^{i=i_{max}}\sum_{j=j_{min}(i)}^{j=j_{max}(i)}g\left(v(i,j)\right)
\end{equation}
où $f$ représente la fonction de valeurs de niveaux de gris à sommer.

En posant 
\begin{equation}
  \label{eq-cumuls1}
  T_g(y,\tau) = \sum_{j=0}^{\tau}g\left(v(y,j)\right)
\end{equation}
L'équation \eqref{eq-sommes-gene} devient 
\begin{equation}
  \label{eq-somme-cumuls2}
  S_f(\Omega_t) = \sum_{i=i_{min}}^{i=i_{max}}\left[T_g(i,j_{max}(i))-T_g(i,j_{min}(i)-1)\right]
\end{equation}
qui représente  une sommation sur le contour $\Gamma$ que l'on peut écrire 
\begin{equation}
  \label{eq-somme-contour}
  S_f(\Omega_t) = \sum_{(i,j)\in \Gamma}C(i,j)\gamma(i,j)
\end{equation}
où $C(i,j)$ est un coefficient lié à la direction du contour au point $(i,j)$ et $\gamma(i,j)$ prend sa valeur selon les règles suivantes

\begin{equation}
  \label{eq-coefC}
  \gamma(i,j)=
  \begin{cases}
    T(i,j)  & \text{si $C(i,j)= 1 $}\\
    T(i,j-1)& \text{si $C(i,j)= -1$}\\
    0       & \text{si $C(i,j)= 0$}\\
  \end{cases}
\end{equation}

La valeur de $C(i,j)$ est déterminée pour chaque pixel d'indice $l$ du contour en considérant les pixels d'indices $l-1$ et $l+1$ qui définissent les deux vecteurs $f_{in}$ et $f_{out}$ et leur code selon le codage de Freeman, comme l'illustre la figure \ref{fig-freeman}. La table \ref{tab-freeman} donne les valeurs de $C(i,j)$ selon les valeurs des codes de Freeman des vecteurs $f_{in}$ et $f_{out}$.
\begin{figure}[htb]
  \centering
  \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/codage-freeman.png}
  \caption{À gauche : détermination des vecteurs $f_{in}$ et $f_{out}$. À droite : code de Freeman d'un vecteur en fonction de sa direction, l'origine étant supposée au pixel central, en noir. }
  \label{fig-freeman}
\end{figure}

\begin{table}[htb]
  \centering
\begin{tabular}[htb]{ccccccccc}
      \toprule
      &\multicolumn{8}{c}{$f_{out}$}\\
      \cmidrule(r){2-9}
      $f_{in}$& 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 \\
      \midrule
      0     &0&0&0&0&0&-1&-1&-1 \\
      1     &1&1&1&1&1&0&0&0\\
      2     &1&1&1&1&1&0&0&0\\
      3     &1&1&1&1&1&0&0&0\\
      4     &0&0&0&0&0&-1&-1&-1\\
      5     &0&0&0&0&0&-1&-1&-1\\
      6     &0&0&0&0&0&-1&-1&-1\\
      7     &0&0&0&0&0&-1&-1&-1\\
      \bottomrule
\end{tabular}
   \caption{Valeur du coefficient $C(i,j)$ en fonction des valeurs des codes de Freeman des vecteurs $f_{in}$ et $f_{out}$.}
      \label{tab-freeman}
\end{table}

L'intérêt de cette transformation est majeur :
\begin{itemize}
\item La sommation en deux dimensions sur la région $\Omega_t$ est ainsi réduite à une sommation à une dimension sur le contour $\Gamma$.
\item Les valeurs $T_g(i,j)$ peuvent être calculées préalablement à la phase de segmentation proprement dite. Pour le cas gaussien qui nous intéresse, cela revient à pré-calculer les trois images \textit{cumulées} $S_1$, $S_x$ et $S_{x^2}$ définies par
  \begin{alignat}{4}
    \label{eq-img-cumul}
    S_1(i,j) &= \sum_{x=0}^jx & \quad \text{,}\quad S_x(i,j) &= \sum_{x=0}^jv(i,x) & \quad \text{et}&\quad & S_{x^2}(i,j) &= \sum_{x=0}^jv(i,x)^2 
  \end{alignat}
\item Les valeurs du coefficient $C(i,j)$ se calculent très facilement durant la génération du contour $\Gamma$.
\end{itemize}


Par ailleurs, le choix d'un contour polygonal permet également d'améliorer l'efficacité de l'algorithme car lors de la phase de segmentation, le déplacement d'un sommet du polygone n'influe que sur les pixels des deux segments qui s'y rapportent, réduisant ainsi la quantité de calculs à effectuer à chaque nouvelle déformation du contour.  

\begin{upminfo}
  L'approche décrite dans ce chapitre n'est valide que si les segments formant le polygone du contour ne se croisent pas. Il est donc nécessaire, lors de la convergence de la segmentation, d'empêcher les croisements de segments. Une solution simple a été proposée dans \cite{ChesnaudRB99} et nous l'avons parallélisée dans le cadre des travaux présentés ici. 
\end{upminfo}

\subsection{Implémentation séquentielle}
Un des inconvénients des algorithmes de type \textit{snake} est l'influence du contour initial sur la convergence de la segmentation. Pour pallier simplement ce défaut, une technique progressive est adoptée, en initialisant le contour avec peu de sommets (4) puis en augmentant le nombre au fur et à mesure de la convergence. L'algorithme \ref{algo-snake-cpu1} décrit macroscopiquement la solution mise \oe uvre tandis que l'algorithme \ref{algo-snake-cpu2} en présente les détails.

\begin{algorithm}
\label{algo-snake-cpu}
\caption{Principe mis en \oe uvre pour la convergence du snake polygonal}
  Calculer les images cumulées\;
  Initialiser le contour avec 4 sommets\;
  \Repeter(\tcc*[f]{niveau contour}){aucun sommet ne peut être ajouté}{
      \Repeter(\tcc*[f]{niveau sommet}){aucun sommet ne peut être déplacé}{
        Déplacer chaque sommet autour de sa position actuelle.\;
        Déplacer le sommet vers la position induisant le meilleur GL\;
     }
      Ajouter un sommet au milieu de chaque \textit{grand} segment\;
  }
\end{algorithm}

\begin{algorithm}[h]
\caption{Détail de l'implémentation du snake polygonal} 
\label{cpualgo}
   Lire l'image $\bar{v}$\;
   Calculer les images cumulées $S_1$, $S_x$ $S_{x^2}$ \nllabel{li-img-cumul}\tcc*[r]{en parallèle via SSE2} 
   $n \leftarrow 0$ \tcc*[r]{indice de boucle niveau contour}
   $N_n \leftarrow 4$ \tcc*[r]{nombre de n\oe uds}
   $\Gamma \leftarrow \{\Gamma_0,\Gamma_1,\Gamma_2,\Gamma_3\} $\;
   $d \leftarrow d_{max}$ \tcc*[r]{pas de déplacement des n\oe uds}
   $l_{min} = 32$ \tcc*[r]{longueur mini des segments sécables}
   $\Gamma_i \leftarrow \Gamma_0$ \tcc*[r]{sommet courant}
   $GL_{ref} \leftarrow GL(\Gamma, N_n, \bar{v}, S_y, S_{y^2})$ \tcc*[r]{voir à partir de ligne 18 pour le détail}
  \Repeter(\tcc*[f]{niveau contour}){$N_{add}=0$}{
    $N_{add}\leftarrow 0$\;
    \Repeter(\tcc*[f]{niveau n\oe ud}){$N_{move}=0$}{
      $N_{move}\leftarrow 0$\;
      \Pour{$i=0$ à $i=N_n-1$}{
        Calculer les positions $\{\Gamma_i^0, \dots, \Gamma_i^7\}$ \tcc*[r]{les 8 voisins de $\Gamma_i$ }
        \Pour{$w=0$ à $w=7$}{
          Soustraire à $GL_{ref}$ la contribution des segments $\Gamma_{i-1}\Gamma_i$ et $\Gamma_{i}\Gamma_{i+1}$\;
          Discrétiser les segments $\Gamma_{i-1}\Gamma_i^w$ et $\Gamma_{i}^w\Gamma_{i+1}$\nllabel{li-bresen}\;
          Lire dans $S_1$, $S_x$ et $S_{x^2}$ les contributions des pixels de $\Gamma_{i-1}\Gamma_i^w$ et $\Gamma_{i}^w\Gamma_{i+1}$\nllabel{li-contrib-seg-deb}\;
          Calculer les directions et lire les codes de Freeman \;
          Calculer $GL_w$ incluant les contributions de $\Gamma_{i-1}\Gamma_i^w$ et $\Gamma_{i}^w\Gamma_{i+1}$ \nllabel{li-contrib-seg-fin}\;
          \lSi{$GL_w > GL_{ref}$}{
            $GL_{ref} \leftarrow GL_w$\;
            $\Gamma_i \leftarrow \Gamma_i^w$\;
            $N_{move} \leftarrow N_{move}+1$\;
          }
        }
      }
     $l \leftarrow l+1$\;
   }
   \PourCh{segment $\Gamma_i\Gamma_{i+1}$}{
     \Si{$\|\Gamma_i\Gamma_{i+1}\| > l_{min}$}{
       Ajouter un n\oe ud au milieu de $\Gamma_i$ et $\Gamma_{i+1}$\;
       $N_{add} \leftarrow N_{add}+1$\;
     }
   }
   $N_n \leftarrow N_n + N_{add}$\;
   \lSi{$d > 1$}{ $d \leftarrow d/2$ } \lSinon{ $d \leftarrow 1$ }\;
   $GL_{ref} \leftarrow GL(\Gamma, N_n, \bar{v}, S_y, S_{y^2})$ \;
  }
\end{algorithm}

Les différentes sommations nécessaires au calcul de la valeur du critère $GL$ sont effectuées en parallèle à l'aide du jeu d'instructions SSE2. La taille des registres utilisables est de 128 bits et permet ainsi de traiter des images de 4096$\times$4096 pixels dont les niveaux de gris sont codés sur 16 bits. Cela ne laisse toutefois que 12 bits pour le codage des segments du contour et limite ainsi leur longueur à 4096 pixels. L'organisation d'un registre SSE 128 bits est donc la suivante :
\begin{itemize}
\item 24 bits pour les sommes à opérandes dans $S_1$
\item 24+16 = 40 bits pour les sommes à opérandes dans $S_x$
\item 24+32 = 60 bits pour les sommes à opérandes dans $S_x^2$
\end{itemize}
Soit un total de 124 bits, qui peuvent donc être représentés par un registre SSE2.

\subsection{Performances}
Les images de 1024$^2$ pixels de la figure \ref{fig-snakecpu-cochon512} montrent l'évolution du contour lors de la segmentation d'une image photographique prise en faible éclairement et bruitée artificiellement par un bruit gaussien d'écart type 25. Les paramètres de la séquence sont fixés empiriquement aux valeurs $d_{max}=16, l_{min}=8$.
Les temps d'exécution indiqués sont mesurés sur Intel Xeon E5530-2.4GHz with 12Go RAM et sont les valeurs moyennes obtenues pour 10 exécutions.
\begin{figure}
  \centering
  \subfigure[Initialisation : 4 n\oe uds]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im000.png}}
  \subfigure[Itération 1 : 8 n\oe uds 3~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im001.png}}
  \subfigure[Itération 2 : 16 n\oe uds 1~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im002.png}}
  \subfigure[Itération 3, 32 n\oe uds 1~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im003.png}}\\
  \subfigure[Itération 7 : 223 n\oe uds 3~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im007.png}}
  \subfigure[Itération 10 : 244 n\oe uds 3~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im010.png}}
  \subfigure[Itération 13 : 256 n\oe uds 3~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im013.png}}
  \subfigure[Itération 14 : 256 n\oe uds 3~ms]{\includegraphics[height=3.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im014.png}}
  \caption{Évolution du contour lors de la segmentation d'une image de 512$^2$ pixels. La convergence est obtenue à l'itération 14 après 44~ms pour un total de  256 n\oe uds.}
 \label{fig-snakecpu-cochon512}
\end{figure}

La dépendance vis à vis du contour initial qui est un des principaux soucis liés au snake est ici fortement relativisée. La figure \ref{fig-snakecpu-compinit} montre le contour final segmentant l'image de test de la figure \ref{fig-snakecpu-cochon512} à partir d'un état initial très éloigné du précédent et \textit{a priori} très défavorable compte tenu du fait qu'il est loin de la cible et sans intersection avec elle. Toutefois, le contour final obtenu est très proche de celui obtenu à partir d'un état initial englobant la cible, malgré un n\oe ud qui s'est ``accroché'' au bord de l'image. La convergence est également plus longue à obtenir dans ce cas avec Un total de 17 itérations en 87~ms et 273 n\oe uds. 

\begin{figure}
  \centering
  \subfigure[Initialisation 2 ]{\includegraphics[height=4.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/snakecpu-init2.png}}
  \subfigure[Contour final 2 : 273 n\oe uds 87~ms]{\includegraphics[height=4.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/snakecpu-compinit2.png}}
  \subfigure[Contour final 1 : 256 n\oe uds 44~ms]{\includegraphics[height=4.5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/im014.png}}
  \caption{Influence du contour initial sur la segmentation. Le contour final 1 est celui de la figure \ref{fig-snakecpu-cochon512}.}
  \label{fig-snakecpu-compinit}
\end{figure}

La dimension de l'image à traiter a également un effet sur le résultat et naturellement sur le temps de calcul. Si l'on conserve les mêmes paramètres d'optimisation que pour la segmentation de l'image 512$^2$ pixels et un contour initial dont les cotés sont à une distance des bords équivalente à 10\% des cotés de l'image, le résultat sur une image identique de 4000$^2$ pixels est  obtenu en 1.3~s avec 1246 n\oe uds ; il est reproduit  à la figure \ref{fig-snakecpu-cochon4ka}. Le nombre de pixels appartenant à la région cible est tel que l'amplitude des déplacements autorisés pour chaque n\oe ud ($d$) peut se révéler trop faible vis à vis du seuil d'acceptation des mouvements. On observe que les zones à gradient élevé ne posent pas de problème et sont détourées de la même manière, alors que dans le bas de l'image où figure une zone de gradient faible (ombre), la cible se trouve maintenant quelque peu surévaluée en surface là ou elle était plutôt sous évaluée dans l'image en 512$^2$ pixels. 
On parvient à un résultat très proche beaucoup plus rapidement en adaptant les paramètres à la taille de l'image, comme le montre par exemple la segmentation de la figure \ref{fig-snakecpu-cochon4kb}, effectuée avec $d_{max}=128$ et $l_{min}=32$ et qui converge vers un contour de 447 n\oe uds en moins de 0.7~s.
Au delà des 16 millions de pixels (4000$^2$ pixels), l'implémentation séquentielle est toujours possible mais doit se priver des instructions SSE. Nous avons, avec leur accord, adapté le code des auteurs en ce sens et réalisé les mesures pour des tailles allant jusqu'à 150~MP. La table \ref{tab-snakecpu-speed-size} en synthétise les résultats en distinguant chaque fois le temps pris par les pré-calculs et celui nécessaire à la convergence de la segmentation. 

\begin{figure}
  \centering
  \subfigure[$d_{max}=16$ et $l_{min}=8$, 1246 n\oe uds en 1.3~s]{\label{fig-snakecpu-cochon4ka}\includegraphics[height=5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/snakecpu-cochon4k.png}}
  \subfigure[$d_{max}=128$ et $l_{min}=32$, 447 n\oe uds en 0.7~s]{\label{fig-snakecpu-cochon4kb}\includegraphics[height=5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/snakecpu-cochon4k-128-32.png}}
  \caption{Segmentation de l'image de test en 4000$^2$ pixels.}
  \label{fig-snakecpu-cochon4k}
\end{figure}


\begin{table}
  \centering
  \begin{tabular}{rrrr}
      \toprule
      &\multicolumn{3}{c}{Taille de l'image}\\
      &\multicolumn{3}{c}{(millions de pixels)}\\
      \cmidrule(r){2-4}
      & 15 & 100 & 150 \\
      \midrule
      Pré-calculs &0,13&0,91&1,4\\
      Segmentation&0,46&3,17&4,3\\
      {\bf Total} &0,51&4,08&5,7\\
      \bottomrule
\end{tabular}
   \caption{Performances (en secondes) de la segmentation par snake polygonal sur CPU en fonction de la taille de l'image à traiter. Le temps sont obtenus avec la même image de test dilatée et bruitée et un contour initial carré dont la distance aux bords est proportionnelle à la taille de l'image.}
      \label{tab-snakecpu-speed-size}
\end{table}


Enfin, il faut aussi considérer les tailles relatives de la cible et de l'image. Ainsi, si on fait l'hypothèse d'une cible de petite taille ``noyée'' dans une image de grandes dimensions, les résultats de la segmentation seront impactés en raison, cette fois, d'une moindre adaptation à la cible lors des toutes premières itérations, les plus grossières, où le nombre de n\oe uds et réduit et le pas de déplacement potentiellement grand vis à vis de la cible. Ce cas de figure est illustré par la segmentation reproduite à la figure \ref{fig-snakecpu-cochon4kc3} et qui met en évidence une qualité moindre par la confusion des zones les plus sombres de la cible avec le fond.

\begin{figure}[h]
  \centering
  \includegraphics[height=5cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/cochon4kc-128-8.png}
  \caption{Segmentation de l'image de test en 4000$^2$ pixels avec une cible de petite taille. Le contour initial est celui utilisé à la figure \ref{fig-snakecpu-cochon4k}.}
  \label{fig-snakecpu-cochon4kc3}
\end{figure}
 


\section{Implémentation parallèle GPU du snake polygonal}
L'analyse de l'exécution du programme séquentiel révèle la prépondérance des blocs fonctionnels suivants, dans l'ordre d'importance, qui occupent à eux seuls plus de 80\% du temps total d'exécution :
\begin{itemize}
 \item Le calcul de la contribution des segments (lignes \ref{li-contrib-seg-deb} à \ref{li-contrib-seg-fin} dans l'algorithme \ref{cpualgo}) 
  \item La génération des trois images cumulées, avant le début des itérations (ligne \ref{li-img-cumul}).
  \item La discrétisation des segments définis par les coordonnées de leurs extrémités (ligne \ref{li-bresen}).
\end{itemize}
Cette proportion est globalement conservée lorsque la taille de l'image à traiter varie, comme le montre le graphique de la figure \ref{fig-snakecpu-chronos1} 

\begin{figure}
  \centering
  \includegraphics[width=0.5\linewidth, height=0.25\linewidth]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter3/img/cpu_profil.png}
\label{fig-snakecpu-chronos1}
  \caption{Évolution du coût relatif des trois fonctions les plus consommatrices en temps de calcul en fonction de la taille de l'image à traiter.}
\end{figure} 

Si l'effort de parallélisation porte essentiellement sur ces fonctions coûteuses, l'ensemble du traitement est réalisé sur le GPU afin de réduire autant que possible les transferts entre le GPU et le système hôte qui, selon le volume concerné, sont susceptibles de grever considérablement la performance globale. L'hôte ne conserve que l'initiative du transfert initial et le contrôle de la boucle principale, ne nécessitant l'échange que d'un seul octet à chaque itération (le nombre de nouveau n\oe uds $N_{add}$).

Les traitements étant totalement indépendants, nos traitons séparément la parallélisation des pré-calculs et celle de la segmentation.

\subsection{pré-calculs des images cumulées}   
Pour réduire la quantité de mémoire requise, nous avons choisi de ne pas générer l'image $S_1$ mais plutôt d'en calculer les valeurs à la volée. L'expression en est simple et le temps pris par les opération élémentaires qu'elle met en jeu est largement compensé par le gain obtenu en économisant les accès mémoire qui auraient été nécessaires, ce qui n'est pas le cas des deux autres images $S_x$ et $S_x^2$ dont le calcul est quant à lui réalisé en appliquant la méthode des \textit{prefixsums} inclusifs décrite dans \cite{BlellochTR90}. 

La prise en compte de grandes images nécessite cependant de les traiter par tranches successives calibrées en fonction de la contrainte limitative qui peut-être :
\begin{itemize}
\item le nombre maximal de threads dans une grille, pour le cas où l'on adopte la règle ``1 pixel par thread''.
\item la quantité de mémoire partagée par blocs de threads, si l'on traite plus d'un pixel par thread.
\end{itemize}

La solution qui apporte les meilleures performances est celle du traitement 1 pixel par thread, associé à un découpage de l'image en tranche équilibrées. Le seuil au delà duquel ce découpage devient nécessaire se situe à 33 millions de pixels pour C1060 et 66 millions pour C2050.