2 La segmentation représente un enjeu particulièrement important dans le domaine du traitement d'image et a ainsi fait l'objet d'abondants travaux et publications touchant les nombreux cas d'analyse dans lesquels une segmentation est utilisée. On peut citer la reconnaissance de formes, la détections et/ou la poursuite de cibles, la cartographie, le diagnostic médical, l'interaction Homme-machine, la discrimination d'arrière plan, etc.
4 On pourrait donner de la segmentation une définition spécifique par type d'usage, mais dans un souci d'unification, nous proposons la formulation générique suivante :\\
5 \og La segmentation consiste à distinguer des zones homogènes au sein d'une image \fg{}.\\
6 Le caractère \textit{homogène} s'entend au sens d'un critère pré-établi, adapté aux contraintes particulières de traitement comme le type de bruit corrompant les images, le modèle d'image ou bien la dimension du signal observé $\bar{v}$ selon que l'image est en couleur ou non. Un tel critère peut ainsi être un simple seuil de niveau de gris ou bien nécessiter de coûteux calculs statistiques dont certains seront détaillés dans les chapitres suivants.
8 Devant la diversité des cas à traiter et des objectifs à atteindre, on sait aujourd'hui qu'à l'instar du filtre unique, la méthode universelle de segmentation n'existe pas et qu'une bonne segmentation est celle qui conduit effectivement à l'extraction des structures pertinentes d'une image selon l'interprétation qui doit en être faite.
10 Les éléments constitutifs de la segmentation sont soit des régions (les segments), soit des contours. Les deux notions sont complémentaires étant donné que les contours délimitent des régions\footnote{Nous excluons ici les contours non-fermés qui peuvent être obtenus à l'aide d'un simple détecteur de bords.}, mais les techniques de calcul basées sur l'un ou l'autre de ces éléments relèvent d'analyses différentes.
12 \section{Les techniques de segmentation orientées régions}
13 Les algorithmes de cette classe s'appuient pour beaucoup sur des techniques de regroupement, ou \textit{clustering}, pour l'identification et le peuplement des régions. Ce lien trouve son origine dans la théorie du \textit{gestalt} \cite{humphrey1924psychology} où l'on considère que la perception conceptuelle s'élabore au travers de regroupements visuels d'éléments.
15 La plupart des approches proposées jusqu'à très récemment consistent à minimiser une fonction d'énergie qui n'a pas de solution formelle et que l'on résout donc à l'aide de techniques numériques souvent itératives.
17 \subsection{Analyse d'histogramme}\label{sec-histo}
18 Les techniques de segmentation les plus simples à mettre en \oe uvre sont dites \og de seuillage \fg{} et basées sur une analyse de l'histogramme des niveaux de gris (ou de couleurs). Elles cherchent à distinguer les différentes classes comme autant d'occurrences représentant des \textit{régions} homogènes.
19 Différents critères peuvent être appliqués pour cette analyse, visant par exemple à maximiser la variance \cite{4310076} ou encore à maximiser le contraste pour déterminer les valeurs pertinentes des seuils.
21 Malgré la multitude de variantes proposées, ces méthodes demeurent peu robustes et présentent l'inconvénient majeur de ne pas garantir la connexité des régions déterminées. On les réserve à des applications très spécifiques où, par exemple, on dispose d'une image de référence dont l'histogramme peut être comparé à celui des images à traiter. C'est le cas de certaines applications de contrôle industriel où la simplicité algorithmique permet de surcroît des implémentations très rapides, voire câblées.
23 Ces techniques peuvent aujourd'hui être considérées comme rudimentaires mais les calculs d'histogrammes et les analyses associées interviennent dans beaucoup d'algorithmes récents parmi les plus évolués et performants.
24 La figure \ref{fig-histo-cochon} illustre le traitement typique de l'histogramme de l'image d'entrée \ref{fig-histo-cochon-a} dans le but de distinguer les deux régions du fond et de la peluche (la cible). La première étape consiste à dresser l'histogramme des niveaux de gris sur tout le domaine de l'image \ref{fig-histo-cochon-b}. Il faut ensuite identifier le seuil de séparation des deux régions supposées, ici, homogènes au sens des valeurs de niveau de gris. Une estimation visuelle peut-être faite, mais on voit immédiatement que même dans une situation aussi claire, le choix du seuil n'est pas évident. Pour un traitement automatique, on peut par exemple proposer la technique itérative présentée par l'Algorithme \ref{algo-histo-cochon} qui conduit à la segmentation de la figure \ref{fig-histo-cochon-c}.
25 Le principe mis en \oe uvre est de fixer arbitrairement une valeur initiale au seuil (ici 128), puis de calculer les \og poids \fg{} respectifs des valeurs situées au-dessus et en-dessous de ce seuil (ligne 6). Cette évaluation de la répartition énergétique permet d'ajuster la valeur du seuil (ligne 8), puis de recommencer jusqu'à convergence, lorsque l'écart entre deux valeurs successives du seuil devient inférieur à une limite fixée à l'avance ($\epsilon$, ligne 9).
27 L'image \ref{fig-histo-cochon-d} est l'image initiale, corrompue par un bruit gaussien de moyenne nulle et d'écart type 25. Les résultats de la segmentation (\ref{fig-histo-cochon-c} et \ref{fig-histo-cochon-f}) de cette image sont clairement insuffisants : les régions identifiées présentent des discontinuités et dans le cas de l'image bruitée, quantité de pixels orphelins subsistent en quantité. Cette technique nécessiterait une étape supplémentaire pour obtenir une segmentation pertinente.
31 \subfigure[Image initiale comportant deux zones : le fond et la peluche (la cible)]{\label{fig-histo-cochon-a} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/cochon256.png}}\quad
32 \subfigure[Histogramme des niveaux de gris]{\label{fig-histo-cochon-b} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/seg_histogramme/histo-cochon256.png}}\quad
33 \subfigure[Image binaire représentant la segmentation. Seuil estimé à 101 après 4 itérations.]{\label{fig-histo-cochon-c} \includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/seg_histogramme/cochon256-seghisto-101-255.png}}\\
34 \subfigure[Image initiale bruitée]{\label{fig-histo-cochon-d} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/cochon256-sig25.png}}\quad
35 \subfigure[Histogramme des niveaux de gris]{\label{fig-histo-cochon-e} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/seg_histogramme/histo-cochon256-sig25.png}}\quad
36 \subfigure[Image binaire représentant la segmentation. Seuil estimé à 99 après 5 itérations.]{\label{fig-histo-cochon-f} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/seg_histogramme/cochon256-sig25-seghisto-99-255.png}}
37 \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par analyse simple d'histogramme. Colonne de gauche : image d'entrée. Colonne centrale : histogramme des niveaux de gris. Colonne de droite : résultat de la segmentation.}
38 \label{fig-histo-cochon}
42 %\SetNlSty{textbf}{}{:}
43 %\SetKwComment{Videcomment}{}{}
44 \caption{Calcul du seuil de séparation des segments de l'histogramme.}
45 \label{algo-histo-cochon}
46 $\overline{h} \leftarrow $ histogramme sur l'image \;
47 $S_{init} \leftarrow 128$ \;
48 $S_k \leftarrow S_{init}$ \;
49 $\epsilon \leftarrow 1$ \;
50 \Repeat{$\|S_k - \frac{1}{2}(\mu_{inf} + \mu_{sup})\| < \epsilon $}{
51 $\mu_{inf}=\displaystyle \frac{\displaystyle\sum_{i<S_k}h_ii}{\displaystyle\sum_{i<S_k}h_i}$ \;
52 $\mu_{sup}=\displaystyle \frac{\displaystyle\sum_{i\geq S_k}h_ii}{\displaystyle\sum_{i\geq S_k}h_i}$ \;
53 $S_k = \frac{1}{2}(\mu_{inf} + \mu_{sup})$ \ ;
57 \subsection{Partitionnement de graphe}
58 Un autre formalisme qui a généré une vaste classe d'algorithmes de segmentation est celui des graphes. Il repose sur l'idée que les régions de l'image sont représentées par les n\oe uds d'un graphe, alors que les liens traduisent les relations de voisinage existant entre les régions, l'idée de base étant d'initialiser le graphe avec un n\oe ud pour chaque pixel. La segmentation est obtenue par partitionnement itératif du graphe, en évaluant les liens et en déterminant ceux à supprimer, et ce jusqu'à convergence.
60 L'essentiel de la problématique réside donc dans la métrique retenue pour évaluer les liens ainsi que dans le critère de sélection et, là encore, la littérature regorge d'une grande variété de propositions.
61 Les premières d'entre elles, qui n'étaient pas spécifiquement dédiées à la segmentation d'images numériques mais au regroupement d'éléments répartis sur un domaine (1D ou 2D), ont été élaborées autour d'une mesure locale des liens basée sur la distance entre les éléments. La réduction du graphe est ensuite effectuée en utilisant un algorithme spécifique, comme le \textit{minimum spanning tree}, dont l'application a été décrite dès 1970 dans \cite{Zahn:1971:GMD:1309266.1309359} et où il s'agit simplement de supprimer les liens \textit{inconsistants}, c'est à dire ceux dont le poids est significativement plus élevé que la moyenne des voisins se trouvant de chaque coté du lien en question.
63 Le principe en a rapidement été étendu aux images numériques en ajoutant l'intensité des pixels au vecteur des paramètres pris en compte dans l'évaluation du poids des liens.
64 D'autres critères de partitionnement ont été élaborés, avec pour ambition de toujours mieux prendre en compte les caractéristiques structurelles globales des images pour parvenir à une segmentation conduisant à une meilleure perception conceptuelle.
65 Le principe général des solutions actuelles repose sur la construction d'une matrice de similarité qui traduit les liens entre les segments et représente le graphe à partitionner.
67 Pour des images en niveaux de gris, l'expression générale des éléments $w_{ij}$ de la matrice de similarité $W$ est :
70 \mathrm{e}^{\|v_i-v_j\|^2/\sigma_v^2}\mathrm{e}^{\|x_i-x_j\|^2/\sigma_x^2} & \text{si $\|x_i-x_j\|<r$}\\
74 On construit également la matrice de connectivité $D$, diagonale et dont les éléments sont :
75 \[d_{i} = \displaystyle\sum_jw_{ij}\]
77 Une famille de méthodes, inspirée par le \textit{graphe optimal} de Wu et Leahy \cite{wu1993optimal}, réalise le partitionnement sur la base des valeurs propres $\lambda_k$ et vecteurs propres $Y_k$ du système
78 \[\left(D-W\right)Y=\lambda DY \]
79 Certains algorithmes proposés plus récemment s'inscrivent dans cette veine \cite{wang2001image,wang2003image,felzenszwalb2004efficient,shi2000normalized}, avec comme principal point faible, une difficulté à trouver un compromis acceptable entre identification de structures globales et préservation des éléments de détails. Cela se traduit dans la pratique par un ensemble de paramètres à régler pour chaque type de segmentation à effectuer.
81 La figure \ref{fig-graph-cochon} montre un exemple de l'application de l'algorithme \textit{normalized cuts} décrit dans \cite{shi2000normalized} et implémenté par Cour, Yu et Shi en 2004. Cette implémentation utilise des valeurs pré-établies des paramètres de calcul de la matrice de similarité produisant de bonnes segmentations de sujets dans les images naturelles, mais requiert de prédéterminer le nombre de segments. Les images de la figure représentent les résultats obtenus avec un nombre de segments variant de 2 à 5 et montrent qu'il est difficile de trouver un compromis acceptable. De plus, les temps d'exécution peuvent devenir très rapidement prohibitifs, même avec des implémentations plus optimisées. Les résultats de la figure \ref{fig-graph-cochon} ont été obtenus en 1,5~s environ (Matlab R2010 sur CPU intel core i5-2520M \@ 2.50GHz - linux 3.2.0)
84 \subfigure[$s = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_2seg.png}}
85 \subfigure[$s = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_3seg.png}}
86 \subfigure[$s = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_4seg.png}}
87 \subfigure[$s = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_5seg.png}}
88 \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par simplification de graphe de type \textit{Normalized cut} pour un nombre $s$ de segments variant de 2 à 5.}
89 \label{fig-graph-cochon}
92 Un autre procédé de partitionnement de graphe, reposant sur le théorème dit du \textit{maximum flow-minimum cut} énoncé par Ford et Fulkerson \cite{ford1955simple} a fait l'objet de beaucoup de travaux, dont les résultats sont comparés dans \cite{boykov2004experimental,chandran2009computational}.
93 Ce procédé est mis en \oe uvre avec de bons résultats dans plusieurs algorithmes, comme le \textit{push-relabel} \cite{cherkassky1997implementing} ou le \textit{pseudoflow} \cite{hochbaum2013simplifications} qui semble aujourd'hui le plus performant.
95 \subsection{kernel-means, mean-shift et apparentés}
96 Parallèlement à la réduction de graphes, d'autres approches ont donné naissance à une multitude de variantes tournées vers la recherche des moindres carrés.
97 Il s'agit simplement de minimiser l'erreur quadratique totale, ce qui peut se résumer, pour une image de $N$ pixels, en la détermination du nombre $K$ de segments $\Omega_i$ et leur contenu, de sorte à minimiser l'expression
98 \[\sum_{i\in[1..K]}\sum_{x_k\in\Omega_i} \left(v_k-\mu_i\right)^2\]
99 où $\mu_i$ représente la valeur affectée au segment $\Omega_i$, i.e la valeur moyenne des observations $v_k$ sur $\Omega_i$, et $\displaystyle{\bigcup_{i\in[1..K]}\Omega_i=\Omega}$.
101 Cette idée est très intuitive et simple, mais n'a pas souvent de solution explicite, d'autant que le nombre des segments est \textit{a priori} inconnu.
102 Dès 1965, Mac Queen a proposé l'appellation k-means pour cette procédure itérative de regroupement \cite{macqueen1967some} qui débute avec $k$ groupes d'un seul pixel\footnote{Dans son article, MacQueen ne parle pas de pixel mais de point. En effet, la méthode décrite ne visait pas à segmenter des images, mais des données de natures diverses.}
103 pris au hasard, puis d'ajouter chaque point au groupe dont la moyenne est la plus proche de la valeur du point à ajouter. La moyenne du groupe nouvellement agrandi doit alors être recalculée avant le prochain ajout.
104 Cette implémentation est extrêmement simple à mettre en \oe uvre \footnote{Même si en 1965, rien n'était simple à programmer} mais elle possède de nombreux défauts dont le principal est qu'elle ne converge pas nécessairement vers le regroupement optimal, même si on connait la \og bonne \fg{} valeur de $k$.
105 Un autre est d'être très dépendant du choix des $k$ éléments initiaux, en nombre et en position.
107 Toutefois, vraisemblablement du fait de sa simplicité d'implémentation et de son temps d'exécution rapide, la communauté scientifique s'est beaucoup penchée sur cette méthode pour en compenser les défauts, jusqu'à en faire une des plus employées, en particulier par les statisticiens.
108 On compte aussi beaucoup de variantes telles les \textit{k-centers} \cite{agarwal2002exact} et les \textit{k-médians} \cite{arora1998approximation} qui n'employent pas la moyenne arithmétique comme expression du \og centre \fg{} d'un segment.
109 Des solutions ont aussi été apportées pour l'estimation de $K$ en employant, par exemple, un critère BIC (Bayesian Information Criterion) pour choisir la meilleure valeur de $K$ dans un intervalle donné \cite{pelleg2000x}.
110 À titre d'illustration et de comparaison, l'image de la peluche a été traitée par une implémentation naïve de l'algorithme original des \textit{k-means} en donnant successivement au nombre de segments les valeurs $K=2$ à $K=5$. Les résultats sont reproduits à la figure \ref{fig-kmeans-cochon} et montrent encore une fois l'influence de $s$ sur la segmentation.
113 \subfigure[$K = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_2seg.png}}
114 \subfigure[$K = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_3seg.png}}
115 \subfigure[$K = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_4seg.png}}
116 \subfigure[$K = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_5seg.png}}
117 \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme \textit{k-means} pour un nombre $K$ de segments variant de 2 à 5. Chaque couleur est associée à un segment. Les couleurs sont choisies pour une meilleure visualisation des différents segments.}
118 \label{fig-kmeans-cochon}
121 Dès 1975, Fukunaga et Hostetler \cite{fukunaga1975estimation} avaient décrit un algorithme générique permettant de déterminer le nombre de segments, ou modes, ainsi que les points, ou pixels, qui les composent. Leur algorithme cherche, pour ce faire, à localiser les $k$ positions où le gradient de densité s'annule.
122 Il utilise à cet effet un voisinage pondéré (ou \textit{kernel}) et détermine le centre de masse des segments en suivant itérativement le gradient de densité dans le voisinage de chaque élément du domaine. Lorsque l'algorithme a convergé, les $k$ segments sont identifiés et contiennent chacun l'ensemble des points qui ont conduit à leurs centres de masse respectifs.
123 Étonnement, malgré ses qualités intrinsèques, cet algorithme du \textit{mean-shift} est resté longtemps sans susciter de grand intérêt, jusqu'à l'étude de Cheng \cite{cheng1995mean} qui en a valorisé les propriétés et établi des liens avec d'autres techniques d'optimisation ou de filtrage, comme la descente/montée de gradient ou le floutage.
125 Comaniciu et Peer en ont alors étendu l'étude et proposé une application à la segmentation utilisant l'espace colorimétrique CIELUV \cite{foley1994introduction} et ont montré qu'elle permettait une meilleure identification des segments de l'image \cite{comaniciu1999mean,comaniciu2002mean}.
126 Une implémentation de la variante proposée par Keselman et Micheli-Tzanakou dans \cite{keselman1998extraction} appliquée à notre image de test, fournit les résultats reproduits à la figure \ref{fig-meanshift-cochon}. Pour se rapprocher des traitements précédents, nous avons identifié, par essais successifs, les tailles de voisinage conduisant à des nombres de segments identiques à ceux des figures précédentes (de 2 à 5), le volume minimal admis pour un segment étant arbitrairement fixé à 100 pixels.
129 \subfigure[$r=100 \Rightarrow K = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r100m100.png}}
130 \subfigure[$r=50 \Rightarrow K = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r50m100.png}}
131 \subfigure[$r=35 \Rightarrow K = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r35m100.png}}
132 \subfigure[$r=25 \Rightarrow K = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r25m100.png}}
133 \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme \textit{mean-shift} pour un rayon de voisinage $r$ de 100, 50, 35 et 25 pixels permettant d'obtenir un nombre $K$ de segments variant respectivement de 2 à 5. Le volume minimal admis pour un segment est fixé à 100 pixels. Chaque couleur est associée à un segment. Les couleurs sont choisies pour une meilleure visualisation des différents segments.}
134 \label{fig-meanshift-cochon}
137 Il est à noter que les segmentations basées sur des algorithmes de \textit{clustering} comme ceux que l'on vient de présenter nécessitent le plus souvent une phase supplémentaire de génération des frontières inter-segments et d'affectation de la valeur de chaque segment aux éléments qui le composent.
138 Par ailleurs, dans les deux cas du \textit{k-means} et du \textit{mean-shift}, chaque itération génère une réduction de la variance (due au moyennage) et on peut donc rapprocher ces techniques de celles de réduction de bruit par minimisation de variance.
140 \section{Les techniques de segmentation par contours actifs, ou snakes}
141 Contrairement aux précédentes techniques et comme leur nom le laisse deviner, les éléments constitutifs de ces méthodes sont cette fois des \textit{contours} et non plus des \textit{régions}. De fait, ils définissent nativement une segmentation de l'image.
142 Le principe général est de superposer une courbe paramétrique $S$ à l'image, le \textit{snake}, puis de lui appliquer des déformations successives destinées à rapprocher le \textit{snake} des contours de l'objet. Les déformations à appliquer sont guidées par l'évaluation d'une fonction d'énergie $E_{snake}$ prenant en compte :
144 \item l'énergie interne $E_{int}$ de la courbe, fonction de son allongement de sa courbure et qui vise à régulariser les courbes obtenues.
145 \item l'énergie externe $E_{ext}$ liée à l'image, fonction de la proximité de la courbe avec les zones de fort gradient et éventuellement une contrainte fixée par l'utilisateur comme des points imposés par exemple.
147 L'expression générique peut alors s'écrire
148 \[E_{snake} = E_{int}+E_{ext}\]
150 \[E_{int} = \displaystyle\sum_{s\in S} \frac{1}{2}\left(\alpha\left|\frac{\partial x_s}{\partial s}\right|^2
151 +\beta \left|\frac{\partial^2x_s}{\partial s^2}\right|\right)ds\]
153 \[E_{ext} = \displaystyle\sum_{s\in S} -\left|\nabla\left[G_{\sigma}(x_s)\ast v_s\right]\right|^2ds\]
155 L'objectif de l'algorithme du \textit{snake} est de trouver une courbe $S$ qui minimise l'énergie totale $E_{snake}$.
156 Ici encore, la résolution du problème revient à minimiser une fonction sous contrainte et les diverses techniques de résolution numérique peuvent s'appliquer comme pour les autres classes d'algorithmes itératifs présentés précédemment. Ici encore, il faut composer avec un nombre de paramètres à régler assez important. Notons également que dans le cas général, les paramètres notés $\alpha$ et $\beta$, que l'on qualifie aussi d'élasticité et de raideur, sont des fonctions de l'abscisse curviligne $s$. La fonction $G_{\sigma}$ est la fonction d'attraction aux forts gradients de l'image.
158 Dans sa version originale proposée par Kass \textit{et al.} en 1988 \cite{KassWT88}, l'algorithme dit du \textit{snake} présente l'intérêt de converger en un nombre d'itérations assez réduit et permet de suivre naturellement une \textit{cible} en mouvement après une convergence initiale à une position donnée, chaque position de convergence fournissant une position initiale pertinente pour la position suivante.
159 Toutefois, il se montre particulièrement sensible à l'état initial de la courbe et requiert souvent de celle-ci qu'elle soit assez proche de l'objet à \og entourer \fg{}, sous peine de se verrouiller dans un minimum local.
160 La sensibilité au bruit n'est pas non plus très bonne du fait de la formulation locale de l'énergie.
161 Les \og concavités \fg{} étroites ou présentant un goulot d'étranglement marqué sont par ailleurs mal délimitées.
162 Enfin, la fonction d'énergie étant calculée sur la longueur totale de la courbe, cela pénalise la bonne identification des structures de petite taille vis à vis de la longueur totale de la courbe.
163 La figure \ref{fig-snake-tradi-cochon} illustre ces défauts en montrant quelques états intermédiaires ainsi que le résultat final d'une segmentation réalisée à partir d'un contour initial circulaire et des paramètres à valeurs constantes et réglés empiriquement, en employant la méthode du \textit{snake} original.
164 On voit que la convergence est assez rapide mais que le contour ainsi déterminé ne \og colle \fg{} pas bien à l'objet que l'on s'attend à isoler.
167 \subfigure[Les états initial et suivant chacune des trois premières itérations]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it3.png}}
168 \subfigure[L'état du contour après la septième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it7.png}}
169 \subfigure[L'état du contour après la dixième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it10.png}}
170 \subfigure[L'état du contour après la centième itération. C'est le contour final.]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_result.png}}
171 \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme dit du \textit{snake}, dans sa version originale. Les paramètres d'élasticité, de raideur et d'attraction ont été fixés respectivement aux valeurs 5, 0.1 et 5. }
172 \label{fig-snake-tradi-cochon}
175 Il est cependant possible de contrôler la finesse de la segmentation mais au prix de temps de calculs qui peuvent devenir très longs.
176 Parmi les variantes élaborées qui tentent de pallier ces défauts, les plus intéressantes sont :
178 \item le \textit{balloon snake}, conçu pour remédier au mauvais suivi des concavités en introduisant une force supplémentaire de pression tendant à \textit{gonfler} le snake jusqu'à ce qu'il rencontre un contour suffisamment marqué. Cela suppose toutefois que l'état initial de la courbe la situe entièrement à l'intérieur de la zone à segmenter. Cette méthode est donc surtout employée dans des applications semi-automatiques où l'utilisateur définit au moins une position et une taille initiales pour la courbe.
179 \item le \textit{snake} GVF (pour Gradient Vector Flow), dont le but est de permettre qu'une initialisation lointaine de la courbe ne pénalise pas la segmentation. Une carte des lignes de gradient est établie sur tout le domaine de l'image et sert à intégrer une force supplémentaire dans l'énergie totale, qui attire la courbe vers la zone de fort gradient.
180 \item les \textit{level-sets}, dont la particularité est de ne pas employer directement une courbe paramétrique plane mais de définir l'évolution des frontières comme l'évolution temporelle de l'ensemble des points d'élévation nulle d'une surface 3D soumise à un champ de force.
181 Les propriétés des contours actifs par \textit{level-sets} se sont révélées intéressantes, en particulier la faculté de se disjoindre ou de fusionner, mais les temps de calcul se sont avérés très pénalisants.
182 Après la formulation initiale de Osher et Sethian en 1988 \cite{osher1988fronts}, plusieurs façons de réduire le coût du calcul ont été formulées, dont les plus importantes restent la technique dite \textit{narrow band} \cite{adalsteinsson1994fast} (bande étroite) qui ne calcule à chaque itération que les points dans une bande étroite autour du plan $z=0$ de l'itération courante et celle du \textit{fast marching} \cite{sethian1996fast} qui s'applique dans le cas particulier d'une évolution monotone des fronts.
183 \item les \textit{snake} orientés régions, qui visent essentiellement à mieux caractériser les zones à segmenter et améliorer la robustesse vis-à-vis du bruit en employant une formulation de l'énergie calculée sur le domaine complet de l'image \cite{cohen1993surface, ronfard1994region}. Les premiers résultats confirment la qualité de cette méthode, mais la nécessité d'effectuer les calculs sur l'image entière générait des temps de traitement prohibitifs jusqu'à ce que les auteurs de \cite{ChesnaudRB99} proposent une amélioration algorithmique exacte permettant à nouveau un calcul en 1D, le long de la courbe, moyennant une simple étape initiale générant un certain nombre d'images intégrales. Des travaux ultérieurs leur ont permis d'étendre la portée de leur algorithme à la segmentation multirégions \cite{GallandBR03,GermainR01}, puis de proposer des variantes non paramétriques ne nécessitant aucune connaissance \textit{a priori} sur le bruit. La plus récente de ces versions présente des temps de traitement très courts, inférieurs à 10~ms pour l'image de la peluche en 256$\times$256 pixels \cite{galland2005minimal,5767240}. La section \ref{snake-formulation} donnera une description détaillée de la variante de \cite{ChesnaudRB99} qui a servi de base à nos travaux.
185 \subsection{Méthodes hybrides}
186 Aujourd'hui, les algorithmes de segmentation les plus performants en termes de qualité emploient des techniques qui tentent de tirer le meilleur parti de plusieurs des méthodes \og historiques \fg{} décrites précédemment.
187 Le meilleur exemple, et le seul que nous citerons, est le détecteur de contour et l'algorithme de segmentation associé proposé par Arbelaez \textit{et al.} en 2010 \cite{arbelaez2011contour}. Les auteurs construisent des histogrammes locaux pour générer une matrice de similitude (\textit{affinity matrix}) et appliquent les techniques liées à la théorie des graphes pour réduire la dimension de l'espace de représentation (calcul des valeurs et vecteurs propres). Ils utilisent ensuite une technique adaptée de l'algorithme intitulé \textit{ligne de partage des eaux} pour en regrouper les segments.
188 Les résultats sont très bons et des implémentations efficaces ont d'ores et déjà été écrites (voir section \ref{sec-seg-hybride}).
190 %peut-être dire deux mots sur le partage des eaux (avec kmeans et meanshift) puisqu'il est employé dans gPb
192 \section{Les implémentations des techniques de segmentation sur GPU}
193 La problématique tant étudiée de la segmentation n'a pas échappé à l'engouement des chercheurs pour les processeurs graphiques modernes. Un certain nombre de travaux proposent ainsi des implémentations GPU plus ou moins directes de méthodes de segmentation tirant parti de l'architecture massivement parallèle de ces matériels.
194 La majorité d'entre elles cherche à répondre à des besoins liés à l'imagerie médicale allant de la simple extraction des contours d'un organe, d'une tumeur, etc., à la mesure de leur volume ; le traitement en 3D n'étant dans ce cas pas un choix mais une obligation, justifiant d'autant plus l'emploi des GPUs.
195 La nature des tissus et les formes à identifier sont extrêmement variées. Ces images sont souvent très bruitées avec des modèles de bruit qui varient selon l'instrumentation employée. Enfin, le diagnostic médical requérant la plus grande précision possible, aucune solution générique satisfaisante de segmentation n'a encore pu émerger dans ce cadre, laissant place à autant d'implémentations adaptées que de besoins médicaux spécifiques.
197 Beaucoup d'algorithmes récents destinés à la segmentation comportent plusieurs phases de calcul et mettent en \oe uvre différents algorithmes réalisant des fonctions élémentaires comme la réduction de bruit ou le calcul d'histogramme.
198 Selon le type de traitement à effectuer sur le GPU, on peut être amené à en concevoir des implémentations parallèles adaptées, ou bien simplement exécuter de multiples instances indépendantes d'une version séquentielle classique du traitement.
199 Dans les deux cas, même si les articles décrivant ces solutions utilisent sans distinction l'expression \og implémentation GPU \fg{}, cela recouvre des réalités et aussi des niveaux de performance souvent très différents.
201 \subsection{Calcul d'histogramme}
202 Étant donné que les segmentations par analyse d'histogramme sont aujourd'hui cantonnées à des applications très particulières et que, dans la pratique, ces traitements sont souvent réalisés par des circuits spécialisés ou programmables de type FPGA, il serait illusoire d'espérer les concurrencer par une solution de type GPU, plus coûteuse, plus volumineuse et vraisemblablement moins robuste.
204 Le calcul d'histogramme est cependant utilisé de manière intensive dans certains algorithmes de haut-niveau, en particulier le \textit{level-set} et le \textit{gPb}. À ce titre, il faut citer les travaux de Fluck \textit{et al.} \cite{fluck2006gpu} qui apportent une réponse efficace au calcul d'histogramme sur GPU, en parvenant à conserver les données dans la mémoire du processeur graphique tout au long de la segmentation par level-set qui était leur motivation initiale \cite{lefohn2003interactive}.
206 Les résultats annoncés ont été obtenus sur un GPU GeForce 7900 et proclament calculer en 1,6~ms les deux histogrammes nécessaires ( 64 classes chacun) sur une image de 256$\times$256 pixels en niveau de gris.
208 \subsection{Partitionnement de graphe}
209 Le domaine du traitement des graphes est très actif et peut fournir des éléments pour la segmentation comme l'implémentation du \textit{minimum spanning tree} décrite dans \cite{Vineet:2009:FMS:1572769.1572796} qui annonce la construction du minimum spanning tree d'un graphe de 5 millions de n\oe uds et 30 millions de liens en moins d'une seconde.
210 En raison de l'indépendance des blocs de threads, la parallélisation GPU des opérations sur les graphes n'est pas simple et peu de travaux font état d'implémentations efficaces mettant en \oe uvre ces techniques : parmi eux, l'implémentation GPU de l'algorithme \textit{push-relabel} qui effectue le partitionnement selon l'approche \textit{min cut/max flow} a donné lieu aux trois versions remarquables que nous détaillons ci-dessous.
212 \item Dans \cite{dixit2005gpu} une approche assez directe est mise en \oe uvre et parvient à \textit{binariser} une image de 1~MP en 29~ms (GeForce 6800GT).
213 \item Les auteurs de \cite{4563095} remarquent qu'après un nombre réduit d'itérations, très peu de n\oe uds changent de segment. En conséquence, certains blocs sont activés sans avoir de traitement à effectuer, ce qui retarde le traitement éventuel des blocs en attente. Pour réduire les effets de ce comportement, un indicateur d'activité est calculé à chaque itération et pour chaque bloc, en se basant sur le nombre de changements de segment qui vient d'y être effectué. À l'itération suivante, seuls les blocs considérés comme \textit{probablement} actifs seront activés, réduisant ainsi la latence globale. Un reparamétrage dynamique du graphe après chaque itération est également effectué selon la méthode décrite par Kohli et Torr \cite{kohli2007dynamic}. Ces optimisations permettent d'atteindre un débit d'environ 30 images de 0.3~MP par seconde sur GTX280, ce qui représente un bond en avant en termes de performance.
214 \item Enfin, Stitch a proposé dans \cite{graphcutscuda} des optimisations plus étroitement liées à l'architecture des GPUs Nvidia en faisant qu'un même thread mette à jour plusieurs liens du graphe et aussi en compactant la représentation mémoire des indicateurs de changement de segment grâce à l'emploi d'un seul bit par lien, ce qui permet de mémoriser l'état de 32 liens dans une seule variable de type entier. Cela a permis aussi d'accélérer la convergence de l'algorithme, comme le montre la courbe de la figure \ref{fig-graphcutscuda} (tirée de \cite{graphcutscuda}), et d'atteindre les 70 images par seconde dans les mêmes conditions que précédemment (sur C1060).
219 \includegraphics[width=\linewidth]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/graphcutscuda_stitch.png}
220 \caption{Évolution du nombre de pixels actifs pour les itérations successives de l'implémentation de l'algorithme push-relabel de \cite{graphcutscuda}. Les petites images montrent la localisation des pixels actifs après chaque itération, en blanc.}
221 \label{fig-graphcutscuda}
224 \subsection{K-means, mean-shift et apparentés}
225 La popularité de l'algorithme des \textit{k-means} a induit des tentatives de portage sur GPU dont \cite{che2008performance} qui a implémenté de manière directe l'étiquetage des éléments ainsi qu'une réduction partielle, par bloc, pour la mise à jour des centres ; la réduction finale étant réalisée par le CPU. Cette solution conduit à un transfert des données à chaque itération et ne permet pas d'atteindre des performances élevées. Le temps annoncé pour l'exécution d'une seule itération sur l'ensemble des 819200 éléments de la base de test KDD-Cup-99 \cite{kddcup99} comportant 23 segments s'élève à 200~ms. Qui plus est, cette durée n'inclut ni la réduction ni les transferts, l'accélération revendiquée semble alors très discutable.
227 Dans \cite{5170921}, l'ensemble des tâches d'étiquetage et de mise à jour des centres est réalisé sur le GPU. L'étape de réorganisation des données est exécutée sur le CPU, mais s'avère moins pénalisante que dans la solution précédente, car elle permet de présenter au GPU des données réorganisées pour l'exécution parallèle de la mise à jour des centres (opération de réduction). Les temps d'exécution par itération sont sensiblement les mêmes que pour \cite{che2008performance} mais ils incluent cette fois l'ensemble des calculs (hors transferts). Les auteurs fournissent aussi des mesures des temps d'exécution à convergence, qui atteignent la vingtaine de secondes pour le même ensemble de test.
229 L'implémentation la plus convaincante de \textit{k-means} reste à notre connaissance celle décrite dans \cite{kmeansgpuopengl} où la totalité du traitement est effectuée sur le GPU, moyennant l'emploi d'une texture par segment de données. Les mesures ont montré que cette multiplication du nombre des textures ne constituait pas un facteur de perte de performance, tout du moins jusqu'aux limites des tests, conduits avec un maximum de 32 segments dans des ensembles de 1 million d'éléments. Sur GPU GeForce 8500GT, les temps d'exécution obtenus dans ces conditions sont de 13,8~ms par itération, avec une dépendance très réduite vis-à-vis du nombre de segments.
231 Un algorithme de \textit{mean-shift} est mis en \oe uvre dans des travaux à orientation non médicale pour la poursuite de cibles dans des séquences vidéo \cite{li2009mean}. L'accélération obtenue par rapport aux implémentations séquentielles existantes n'est que d'un facteur 2. La solution présentée effectue préalablement une réduction de l'espace colorimétrique via un regroupement par la méthode \textit{k-means}, utilisée dans une version séquentielle. Un gain potentiel de performance pourrait être apporté en employant une implémentation GPU du \textit{k-means}, mais serait toutefois limité en raison des itérations nécessaires plus nombreuses pour le traitement \textit{mean-shift}. Par ailleurs, l'implémentation proposée fait un usage intensif de la mémoire partagée et se heurte à sa limite de 16~Ko par bloc, obligeant à réduire la taille des blocs à l'exécution et avec eux, le parallélisme et vraisemblablement aussi la performance de l'application. On peut malgré tout raisonnablement espérer qu'une telle solution présenterait de meilleures performances sur une carte de type Fermi possédant jusqu'à 48~Ko de mémoire partagée par bloc.
233 \textit{Quick shift} est une solution permettant d'obtenir un résultat en une seule passe (sans itérer) par approximation de l'algorithme mean-shift gaussien (dont les masques de pondération sont des gaussiennes). Proposée initialement dans \cite{vedaldi2008quick}, cette technique a été parallélisée sur GPU par ses auteurs et décrite dans \cite{fulkerson2012really}. Les performance y sont obtenues grâce à des approximations, parmi lesquelles la restriction des calculs de pondération à des voisinages de rayon $3\sigma$ (écart type de la gaussienne définissant les coefficients du masque), considérant qu'au delà, les valeurs en sont négligeables.
234 On construit ensuite un arbre des liens entre les pixels, en limitant la recherche à une distance maximale de $\sigma$ et en divisant arbitrairement par 2 la dynamique de l'espace colorimétrique. La segmentation est finalement obtenue par simple partionnnement de l'arbre selon un seuil $\tau$.
236 Pour s'affranchir de la relative petite taille de la mémoire partagée sans pâtir de la grande latence des accès à la mémoire globale du GPU, les auteurs ont choisi d'associer l'image et l'estimation de densité à des textures et ainsi bénéficier du mécanisme de cache 2D.
237 Les expérimentations ont été menées avec différentes valeurs de $\sigma$ et $\tau$ choisies pour les résultats visuels qu'elles induisent et permettent de segmenter une image couleur de 1~MP en environ 1~s avec $\tau=10$ et $\sigma=6$. Toutefois, des valeurs plus petites, requérant moins de calculs, permettent des temps d'exécution beaucoup plus courts. Les courbes présentées permettent d'envisager, pour $\tau=4$ et $\sigma=2$, une réduction par 30 du temps d'exécution, soit environ 33~ms. Une version améliorée récemment, dans laquelle les positions des centres sont stockées en registres, permet selon les auteurs, de diviser encore par 2 les temps d'exécution pour atteindre une segmentation en environ 16,5~ms.
238 La figure \ref{fig-quickshift-yo}, tirée de \cite{fulkerson2012really}, présente quelques segmentations effectuées avec des valeurs différentes, permettant ainsi de juger des effets des variations des paramètres $\tau$ et $\sigma$.
242 \subfigure[Image originale]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-orig.png}}\quad
243 \subfigure[$\tau=10$ et $\sigma=2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s2t10.png}}\quad
244 \subfigure[$\tau=10$ et $\sigma=10$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s10t10.png}}\quad
245 \subfigure[$\tau=20$ et $\sigma=10$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s10t20.png}}\quad
246 \caption{Segmentation d'une image couleur de 512$\times$512 pixels par l'implémentation GPU quick-shift de \cite{fulkerson2012really}.}
247 \label{fig-quickshift-yo}
250 Récemment, Xiao et Liu ont décrit dans \cite{xiao2010efficient} une implémentation de l'algorithme \textit{mean-shift} qui utilise cette fois une construction de \textit{KD-tree} (arbre binaire à K dimensions) pour réduire l'espace colorimétrique et effectuer rapidement les recherches des plus proches voisins. L'ensemble s'exécute sur le GPU et permet ainsi d'obtenir des résultats beaucoup plus probants puisque les auteurs revendiquent une segmentation d'image couleur de 6.6 millions de pixels en 0.2 secondes. Malheureusement, il n'est pas dit combien de segments comprend l'image et il n'est fait référence qu'à une seule image, dont on déduit qu'il s'agit de l'image de la figure \ref{fig-meanshift-castle}, reproduite afin de montrer les différences avec une implémentation standard du \textit{mean-shift}.
254 \subfigure[Image originale]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/castle-meanshift.png}}\quad
255 \subfigure[Image segmentée par mean-shift standard]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/castle-meanshift-std.png}}\quad
256 \subfigure[Image segmentée par mean-shift kd-tree]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/castle-meanshift-kdtree.png}}
257 \caption{Comparaison des segmentations d'une image couleur de 2256$\times$3008 pixels réalisées par \textit{mean-shift} standard et par le \textit{mean-shift kd tree} de \cite{xiao2010efficient}.}
258 \label{fig-meanshift-castle}
261 \subsection{Level set et snakes}
262 Dès 2003, on recense d'importants travaux liés à l'imagerie médicale mettant en \oe uvre des algorithmes \textit{level set} sur GPU. C'est le cas de \cite{lefohn2003inter,lefohn2003interactive} où les auteurs décrivent une solution de visualisation des coupes d'une mesure volumique réalisée par résonance magnétique (IRM) en exploitant pour la première fois le caractère creux du système d'équations à résoudre. Cette variante est nommée \textit{narrow-band} et diffère de la première solution 2D présentée dans \cite{rumpf2001level} qui implémente la version standard des \textit{level set}. En ne transférant au GPU, pour chaque itération, que les petits pavés de données actifs et en les rangeant de manière contiguë en texture pour optimiser les accès en lecture, les auteurs sont parvenu à effectuer, pour des données volumiques de 256$\times$256$\times$175, entre 3,5 et 70 itérations par seconde, qu'il faut comparer aux 50 itérations par seconde en 2D sur image de 128$\times$128 pixels obtenues dans \cite{rumpf2001level}. La limitation principale de cette solution est celle des dimensions maximales admises pour une texture qui était de 2048$\times$2048 pour le GPU ATI Radeon 9800 pro employé (et programmé en openGL, car ni openCL ni CUDA n'étaient encore disponible à l'époque).
264 Les autres solutions GPU proposées depuis sont également basées sur la variante \textit{narrow-band} (bande étroite) des \textit{level-set} \cite{lefohn2005streaming,cates2004gist,jeong2009scalable}, mais seule \cite{jeong2009scalable} s'affranchit des transferts CPU/GPU à chaque itération pour déterminer et transférer les pavés actifs. La solution retenue est d'employer les opérations atomiques pour assurer l'accès exclusif à la liste des pavés en mémoire GPU. Cela permet de descendre à 3~ms par itération pour une image de 512$\times$512 pixels.
266 L'implémentation la plus performante à ce jour est celle décrite dans \cite{Roberts:2010:WGA:1921479.1921499} qui parvient à des itérations dont la durée varie, sur GTX280, de 1,8 à 6,5~ms pour des données volumiques de 256$\times$256$\times$256 pixels, issues d'examens IRM, pour une moyenne de 3,2~ms sur les 2200 itérations de l'exemple fourni (cerveau en 7~s, Figure \ref{fig-l7-brain}). Une optimisation poussée y a été effectuée pour rendre l'algorithme efficace, en particulier grâce à la refonte du code responsable de la détermination des pavés actifs, qui parvient cette fois à déterminer l'ensemble minimal de pavés actifs et à rendre cette détermination efficace sur le GPU en gérant parallèlement plusieurs tampons, chacun associé à une direction particulière en 6-connexité. Une étape de résolution des doublons est ensuite effectuée avant de les compacter de manière contiguë comme cela était déjà fait dans \cite{lefohn2003inter}.Tout cela est réalisé sans recourir à la mémoire partagée qui s'avère complexe, voire impossible, à utiliser efficacement lorsque les éléments à accéder sont très irrégulièrement répartis en mémoire.
268 Ce faisant, les auteurs parviennent à réduire le nombre total cumulé de pavés qu'il est nécessaire de traiter lors des 2200 itérations de la segmentation de l'image d'exemple, avec seulement 294 millions de pavés, alors que l'algorithme \textit{narrow-band} standard devait en traiter 4877 millions. Il est à noter que la durée d'exécution d'une itération dans cette variante dépend plus fortement de la proportion de pavés actifs que pour le \textit{narrow-band} standard. Les deux courbes sont globalement affines et se croisent pour une proportion de pavés actifs proche de 10\%.
269 Si l'on considère que malgré les stratégies adoptées, tenir à jour cette liste de pavés représente encore 77\% du temps de calcul, cela peut représenter une piste pour une optimisation supplémentaire qui, bien qu'apparemment superflue pour l'image du cerveau, pourrait se justifier pour d'autres, comme le suggère le temps de segmentation de 16~s nécessaire pour l'image des reins (Figure \ref{fig-l7-reins}) et de l'aorte, aux dimensions comparables.
273 \subfigure[Cerveau 256$\times$256$\times$256 en 7~s]{\label{fig-l7-brain}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/l7-brain7s.png}}\quad
274 \subfigure[Reins et aorte, 256$\times$256$\times$272 en 16~s]{\label{fig-l7-reins}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/l7-reins16s.png}}
275 \caption{Segmentation d'images issues d'examens IRM par la méthode des level set à bande étroite.}
276 \label{fig-l7-narrow}
279 Les algorithmes de type \textit{snake}, très coûteux en temps de calcul, pouvaient prétendre à bénéficier largement de la technologie des GPU pour améliorer leurs performances, mais seule la variante paramétrique GVF à été véritablement implémentée de manière spécifique et efficace \cite{snakegvf06, bauer2009segmentation, li2011robust, snakegvfopencl12}. Les variantes de type polygonal restent à ce jour sans implémentation GPU, principalement en raison de l'irrégularité des motifs d'accès à la mémoire.
281 Parmi les premières solutions décrites, \cite{snakegvf06} propose une implémentation réalisée en openGL, où les données de gradient sont compactées en texture à quatre canaux (R,V,B,A) de manière à s'affranchir du format 16 bits de la représentation : les deux premiers canaux R et V contiennent les valeurs représentant respectivement les gradients selon $dx$ et $dy$ sous une forme codée selon les valeur des 2 autres canaux B et A.
282 Par ailleurs, une approximation du système linéaire à résoudre est proposée afin de donner une structure bande symétrique à la matrice à inverser, ce qui améliore considérablement l'efficacité des accès aux données au travers du cache.
286 \subfigure[Contour initial]{\label{fig-epaule-init}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/snake-epaule-init-t.png}}\quad
287 \subfigure[Contour final]{\label{fig-epaule-fin}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/snake-epaule-fin-t.png}}
288 \caption{Segmentation d'une image d'épaule en 1024$\times$1024 pixels issue d'un examen IRM par l'implémentation du snake GVF de \cite{snakegvf06}. Le contour est représenté en rouge et son état final est obtenu en 11~s. Le tracé initial du contour a été artificiellement épaissi pour le rendre visible à l'échelle de l'impression.}
292 Les performances annoncées montrent tout d'abord que l'approximation adoptée n'a qu'un impact extrêmement limité sur le résultat de la segmentation avec un écart radial maximal inférieur à 1.3 pixel par rapport au calcul exact effectué sur CPU. Enfin, la segmentation de l'image d'exemple en 1024$\times$1024 pixels s'effectue en un total de 11~s après l'initialisation manuelle reproduite à la figure \ref{fig-snakegvf}. Cela est annoncé comme presque 30 fois plus rapide que l'implémentation CPU de référence, mais demeure beaucoup trop lent pour un usage interactif.
294 Une solution directe employant la transformée de Fourier pour inverser le système à résoudre a été décrite récemment dans \cite{zheng2012fast}, et programmée en employant la bibliothèque openGL. Les exemples fournis montrent des objets segmentés dans des images d'environ 10000 pixels en une durée de l'ordre de la demi seconde.
296 En adaptant sur GPU une variante du snake GVF dite FD-snake (pour \textit{Fourier Descriptors}) \cite{li2011robust} permettant une convergence plus rapide et un calcul parallèle beaucoup plus adapté au GPU, Li \textit{et al.} parviennent à suivre les déformations d'un contour en temps réel dans des images issues d'examens échographiques. Un contour de 100 points peut ainsi converger convenablement en à peine 30~ms. Une contribution supplémentaire de cette implémentation est de permettre une initialisation simplifiée et semi-automatique du contour.
298 La plus aboutie des implémentations actuelles du snake GVF est celle présentée par Smistad \textit{et al.} dans \cite{snakegvfopencl12} où les auteurs ont concentré leur effort sur l'optimisation des accès mémoire lors du calcul du GVF. Ils ont comparé 8 combinaisons possibles impliquant l'emploi des mémoires partagée et de texture ainsi que la représentation des nombres selon le format classique 32 bits ou selon un format compressé sur 16 bits. Il en ressort que l'association la plus performante est celle des textures et du format de données sur 16 bits.
299 Les performances sont alors nettement en hausse avec des segmentations d'images médicales d'IRM de 512$\times$512 pixels effectuées en 41~ms sur Nvidia C2070 et 28~ms sur ATI 5870 (512 itérations). L'implémentation réalisée en openCL permet d'exécuter le code sur les GPUs des deux principaux fabricants.
301 \subsection{Algorithmes hybrides\label{sec-seg-hybride}}
302 Le détecteur de contour \textit{gPb} décrit dans \cite{arbelaez2011contour} et que l'on considère comme la référence actuelle pour la segmentation d'objets et personnages dans des image naturelles, a été implémenté en CUDA par Cantazaro \textit{et al.} et est décrit dans \cite{5459410}. Dans cette implémentation GPU, La qualité des contours extraits est préservée et le temps de traitement est réduit d'un facteur supérieur à 100 par rapport à la version CPU de \cite{arbelaez2011contour} : les contours des images de 0.15~MP de la base de test BSDS \cite{martin2001database} sont ainsi traitées en 2 secondes environ sur GPU C1060.
303 L'apport principal de ces travaux réside dans la solution conçue pour le calcul des histogrammes locaux, qui dans l'algorithme original s'étendaient sur des demi-disques centrés sur chaque pixel. La parallélisation réalisée fait l'approximation de chaque demi-disque en un rectangle de même surface dont un des grands cotés a le centre du disque pour milieu. Les rectangles sont ensuite pivotés par une rotation basée sur la discrétisation de Bresenham \cite{bresenham1965algorithm} pour en aligner les cotés avec les cotés de l'image et pouvoir employer la technique des images cumulées pour calculer rapidement l'histogramme.
304 La figure \ref{fig-gPb} présente quelques résultats d'extraction de contours.
307 \includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/gPb_examples.png}
308 \caption{Extraction de contour par la version GPU de l'algorithme gPb. Les images sont issues de la base BSDS \cite{martin2001database}}
313 La présentation que nous venons de faire des principales techniques de filtrage et de segmentation ainsi que des implémentations sur GPU qui leur ont été consacrées nous ont permis de confirmer le fait que, malgré leur orientation grand public et les langages de haut niveau permettant d'accéder rapidement à la programmation GPU, la parallélisation efficace d'un algorithme séquentiel destiné à s'exécuter sur ces processeurs n'est pas triviale. Le modèle et les contraintes de programmation sont très spécifiques et obtenir un code efficace découle nécessairement d'un compromis qui peut parfois être complexe à mettre au point.
315 Ajoutons que les générations successives de GPU conservent certes des caractéristiques communes mais diffèrent suffisamment quant à la distribution des ressources, rendant toute généralisation vaine et faisant qu'un code optimisé pour un modèle donné peut devenir moins rapide avec un modèle plus récent. Prenons l'exemple du nombre maximal de registres utilisables par thread ; il est de 128 sur GPU C1060 contre seulement de 63 pour un C2070. Un code faisant un usage optimisé des registres sur C1060 pourra s'exécuter plus lentement sur C2070. C'est un cas de figure sur lequel nous reviendrons plus en détail dans le chapitre \ref{ch-median} consacré au filtre médian.
317 Cet état de fait rend les résultats publiés par les chercheurs souvent délicats à interpréter et plus encore à reproduire lorsque l'on souhaite comparer les performances de nos propres codes avec les références du moment, sauf à disposer d'un panel de cartes GPU représentant toutes les évolutions de l'architecture et ce pour au moins les deux grands fabricants de GPUs que sont ATI et Nvidia.
319 Pour aider les développeurs à allouer les ressources de manière optimale, ou tout du moins estimer le degré d'optimisation de leur code à l'aune de la vitesse d'exécution, Nvidia fournit une feuille de calcul appelée \textit{occupancy calculator} dans laquelle on peut entrer les paramètres d'exécution d'un \textit{kernel} parallèle : nombre de registres utilisés par chaque thread, quantité de mémoire partagée, modèle de GPU, dimensions de la grille. Le tableur retourne alors l'indice de charge (l'\textit{occupancy}) qui traduit le rapport, à chaque instant, entre le nombre de warps actifs et le nombre maximal de warps par SM. L'occupancy se traduit donc par un indice compris entre 0 et 100\% et la recherche de performance semble devoir être la recherche de l'occupancy maximale.
321 Toutefois, comme l'a clairement démontré Volkov dans \cite{volkov2010better}, ce paradigme peut aisément être remis en cause et Volkov parvient effectivement à améliorer les performances d'un certain nombre d'exemples génériques dans des conditions de faible valeur d'\textit{occupancy}.
322 Enfin, nous avons pu constater deux grands modèles d'accès aux données : les algorithmes de filtrage usent quasiment tous de la mémoire partagée comme tampon d'accès aux données de l'image en mémoire globale (ou texture) alors que les algorithmes de segmentation performants s'en affranchissent. La raison en est clairement des motifs d'accès très irréguliers et non contigus pour ces derniers, rendant la gestion efficace de la mémoire partagée délicate et potentiellement si coûteuse qu'elle en devient sans intérêt.
323 Les chapitres présentant nos contributions reviendront sur ces aspects en proposant des solutions pour accroître les performances des algorithmes parallèles.
427 % LocalWords: reparamétrage pénalisante partionnnement volumique
428 % LocalWords: compacter
