19 oct

[these_gilles.git] / THESE / Chapters / chapter2 / chapter2.tex

L'étendue des techniques applicables aux images numériques est aujourd'hui si vaste qu'il serait illusoire de chercher à les décrire ici. Ce chapitre présente plus spécifiquement les algorithmes utilisés en présence d'images (fortement) bruitées. Le bruit rend potentiellement délicate l'extraction des informations utiles contenues dans les images perturbées ou en complique l'interprétation, qu'elle soit automatique ou confiée à la vision humaine. 
L'intuition nous incite donc à chercher des méthodes efficaces de pré-traitement pour réduire la puissance du bruit afin de permettre aux traitements de plus haut niveau comme la segmentation, d'opérer ensuite dans de meilleures conditions.           

Toutefois, il faut également considérer que les opérations préalables de réduction de bruit apportent des modifications statistiques aux images et influent donc potentiellement sur les caractéristiques que l'on cherche à mettre en évidence grâce au traitement principal. En ce sens, il peut-être préférable de chercher à employer des algorithmes de haut niveau travaillant directement sur les images bruitées pour minimiser les effets des altérations apportées par les filtres débruiteurs et conserver toute l'information contenue dans les images perturbées.

Enfin, toute opération aussi basique soit elle, prend un certain temps qui réduit potentiellement celui disponible pour le traitement de haut niveau. Lorsque les images à analyser sont de grande taille,  procéder à un débruitage préalable peut s'avérer incompatible avec les contraintes de débit et les temps requis par le traitement de haut niveau.

 Les images auxquelles nous nous intéressons sont généralement les images numériques allant des images naturelles telles que définies par Caselles \cite{Caselles99topographicmaps} aux images d'amplitude issues de l'imagerie radar à ouverture synthétique (ROS ou en anglais SAR) \cite{cutrona1990synthetic}, de l'imagerie médicale à ultrasons (échographie) ou encore biologique dans le cas de la microscopie électronique. 
Ces dispositifs d'acquisition sont naturellement, et par essence, générateurs de bruits divers, inhérents aux technologies mises en \oe uvre au sein de ces systèmes et qui viennent dégrader l'image idéale de la scène que l'on cherche à représenter ou analyser. On sait aujourd'hui caractériser de manière assez précise ces bruits et la section \ref{sec_bruits} en détaille les  origines physiques ainsi que  les propriétés statistiques qui en découlent.
On peut dores et déjà avancer que la connaissance de l'origine d'une image et donc des propriétés des bruits associés qui en corrompent les informations, est un atout permettant de concevoir des techniques de filtrage adaptées à chaque situation. Toutefois, la recherche d'un filtre universel, bien qu'encore illusoire, n'est pas abandonnée, tant les besoins sont nombreux, divers et souvent complexes.    
       
\section{Modèle d'image bruitée}
On considère qu'une image observée, de largeur $L$ pixels et de hauteur $H$ pixels, est un ensemble de $N=HL$ observations sur un domaine $\Omega$ à deux dimensions ($\Omega \subset \mathbb{Z}^2$). À chaque élément de $\Omega$, aussi appelé \textit{pixel}, est associé un indice unique $k \in [\![1;N]\!]$, une position $x_k=(i,j)_k \in\Omega$ et une valeur observée $v_k=v(i,j)_k$.
La valeur observée peut, selon les cas, être de dimension $1$ pour les images représentées en niveaux de gris ou de dimension 3 pour les images couleur représentées au format RVB. Les dimensions supérieures, pour la représentation des images hyperspectrales n'est pas abordé.
L'image observée peut ainsi être considérée comme un vecteur à $N$ éléments $\bar{v}= (v_k)_{k\in [\![1;N]\!]}$.
Les divers traitements appliqués aux images observées ont souvent pour but d'accéder aux informations contenues dans une image sous-jacente, débarrassée de toute perturbation, dont nous faisons l'hypothèse qu'elle partage le même support $\Omega$ et que nous notons $\bar{u}$. L'estimation de $\bar{u}$ réalisée par ces traitements est notée $\widehat{\bar{u}} = (\widehat{u}_k)_{k\in [\![1;N]\!]}$.
Le lien entre $\bar{u}$ et $\bar{v}$ peut être exprimé généralement par la relation $\bar{v}=\bar{u}+\sigma\epsilon$, où $\epsilon \in \mathbb{R}^N$ représente le modèle de perturbation appliquée à $\bar{u}$ et $\sigma$ représente la puissance de cette perturbation qui a mené à l'observation de $\bar{v}$.
Dans le cas général, $\epsilon$ dépend de $\bar{u}$ et est caractérisé par la densité de probabilité (PDF pour probability density function) $p(v|u)$.

\section{Modèles de bruit}\label{sec_bruits}
\subsection{Le bruit gaussien}
Le bruit gaussien est historiquement le plus étudié et celui auquel sont dédiées le plus de techniques de débruitage.
La génération des images numériques au travers les capteurs CMOS et CCD  est le siège de nombreuses perturbations dues à la technologie de fabrication et à la nature du rayonnement dont ils mesurent l'intensité en différents zones de leur surface, appelées \textit{photosites} \cite{mancuso2001introduction,theuwissen2001ccd}.
On distingue en particulier les bruits suivants selon leur origine physique :
\begin{itemize}
\item la non uniformité de réponse des photosites.
\item le bruit de photon
\item le bruit de courant d'obscurité
\item le bruit de lecture
\item le bruit de non uniformité d'amplification des gains des photosites.
\end{itemize}
On trouve des descriptions détaillées des mécanismes concourant à la génération de ces bruits, entre autres dans \cite{healey1994radiometric} et \cite{kodakccd}.  
Dans un certain intervalle usuel d'intensité lumineuse, il est toutefois admis que l'ensemble des ces perturbations peut être représenté par un seul bruit blanc gaussien, de type \textit{additif} (AWGN), dont la densité de probabilité suit une loi normale de moyenne nulle et de variance $\sigma^2$.
On a alors l'expression suivante, où $\sigma >0$ 
\[p(v|u)=\frac{1}{\sqrt{2}\pi\sigma}\mathrm{e}^{-\frac{(v-u)^2}{2\sigma^2}}\]

\subsection{Le speckle}
En imagerie radar, sonar ou médicale, les surfaces que l'on veut observer sont ``éclairées'' par des sources cohérentes. Les propriétés locales de ces surfaces sont  le siège de réflexions multiples qui interfèrent entre elles pour générer un bruit de tavelures, ou speckle, dont l'intensité dépend de l'information contenue dans le signal observé.

Le speckle est ainsi un bruit de type \textit{multiplicatif} qui confère aux observations une très grande variance qui peut-être réduite en moyennant plusieurs  observations, ou vues,  de la même scène. Si $L$ est le nombre de vues, le speckle est traditionnellement modélisé par la PDF suivante :
\[p(v \mid u)=\frac{L^2v^{(L-1)}\mathrm{e}^{-L\frac{v}{u}}}{\Gamma (L)u^L} \]
L'espérance vaut $\mathrm{E}\left[v\right]=u$ et la variance $\sigma^2=\frac{u^2}{L}$ est effectivement inversement proportionnelle à $L$, mais pour le cas mono vue où $L=1$, la variance vaut $u^2$, soit un écart type du signal $v$ égal à sa moyenne.

\subsection{Le bruit ``sel et poivre''}
Le bruit \textit{sel et poivre}, ou bruit \textit{impulsionnel} trouve son origine dans les pixels défectueux des capteurs ou dans les erreurs de transmission. Il tire son nom de l'aspect visuel de la dégradation qu'il produit : des pixels noirs et blancs répartis dans l'image.
Le bruit impulsionnel se caractérise par la probabilité $P$ d'un pixel d'être corrompu. La PDF peut alors être exprimée par parties comme suit, pour le cas d'images en 256 niveaux de gris (8 bits) :

\[p(v \mid u)=
\begin{cases}
\frac{P}{2}+(1-P) & \text{si $v=0$ et $u=0$}\\
\frac{P}{2}+(1-P) & \text{si $v=255$ et $u=255$}\\
\frac{P}{2}       & \text{si $v=0$ et $u \neq 0$}\\
\frac{P}{2}       & \text{si $v=255$ et $u \neq 255$}\\
(1-P)             & \text{si $v=u$ et $u \notin \{0, 255\}$}\\
0                 & sinon
\end{cases}
 \]  

\subsection{Le bruit de Poisson}
Aussi appelé \textit{bruit de grenaille} (shot noise), ce type de bruit est inhérent aux dispositifs de détection des photons. Il devient prépondérant dans des conditions de faible éclairement, lorsque la variabilité naturelle du nombre de photons reçus par un photosite par intervalle d'intégration influe sur les propriétés statistiques du signal.
Le bruit de grenaille est de type multiplicatif et suit une loi de Poisson. La PDF peut s'écrire comme suit :
\[ p(v \mid u)=\mathrm{e}\frac{u^v}{v!}\]

\section{Les techniques de réduction de bruit}
La très grande majorité des algorithmes de réduction de bruit fait l'hypothèse que la perturbation est de type gaussien, même si le développement des systèmes d'imagerie radar et médicale a favorisé l'étude des bruits multiplicatifs du type \textit{speckle} ou \textit{Poisson}.
Un très grand nombre de travaux proposant des méthodes de réduction de ces bruits ont été menés, ainsi que beaucoup d'états de l'art et d'études comparatives de ces diverses techniques, que nous n'avons pas l'ambition d'égaler.

Nous nous focaliserons sur les techniques en lien avec les travaux que nous avons menés et qui ont donné lieu à des implémentations efficaces  susceptibles de fournir des éléments opérationnels rapides pour le pré-traitement des images. 

La figure \ref{fig-ny-noises} montre une image de synthèse issue de la base de test COIL \cite{coil}, supposée sans bruit et qui sera considérée comme référence, ainsi que deux versions bruitées, respectivement avec un bruit gaussien d'écart type 25 et un bruit impulsionnel affectant 25\% des pixels. 
L'indice de qualité le plus employé pour mesurer la similarité entre deux images est le PSNR (pour Peak Signal to Noise Ratio). Il est exprimé en décibels (dB) et se calcule en appliquant la formule  
\[ PSNR = 10log_{10}\left(\frac{D^2}{\displaystyle\frac{1}{N}\sum_{k < N}\left(v_k - u_k\right)^2}\right)\]
si l'on cherche à évaluer le PSNR de l'image observée $\bar{v}$ par rapport à l'image de référence $\bar{u}$. Le nombre $D$ représente la dynamique maximale des images, e.g 255 pour des images en niveaux de gris codés sur 8 bits. 

Cet indicateur seul est cependant insuffisant pour caractériser convenablement la qualité de débruitage d'un filtre, mesure hautement subjective. Un indice global de similarité structurelle (MSSIM pour Mean Structural Similarity Index) a été proposé par Wang \textit{et al.} \cite{Wang04imagequality} et permet, en conjonction avec le PSNR, de garantir une mesure de qualité plus en rapport avec la perception visuelle. Le MSSIM prend ses valeurs dans l'intervalle $[0;1]$ avec une similarité d'autant plus grande que la valeur est proche de 1.  

\begin{figure}
  \centering
  \subfigure[Sans bruit]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/ny256.png}}
  \subfigure[Bruit gaussien $\sigma=25$, PSNR=22.3~dB MSSIM=0.16]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/ny256_gauss25.png}}
  \subfigure[Bruit impulsionnel 25\%, PSNR=9.48~dB MSSIM=0.04]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/ny256_sap25.png}}
  \caption{Images 256$\times$256 en niveau de gris 8 bits utilisées pour l'illustration des propriétés des filtres. a) l'image de référence non bruitée. b) l'image corrompue par un bruit gaussien d'écart type $\sigma=25$. c) l'image corrompue par un bruit impulsionnel à 25\%.}
\label{fig-ny-noises}
\end{figure}

\subsection{Les opérateurs de base}
\subsubsection{Le filtre de convolution}
L'opération la plus employée dans les procédés de traitement d'image est sans doute la convolution. Selon les valeurs affectées aux coefficients du masque, le filtrage par convolution permet de réaliser bon nombre de traitements comme la réduction de bruit par moyennage ou noyau gaussien ou encore la détection de contours. 
Si la fonction définissant le masque de convolution est notée $h$, l'expression générale de la valeur estimée de pixel de coordonnées $(i,j)$ est donnée par
\begin{equation}
\widehat{u}(x, y) = \left(\bar{v} * h\right) = \sum_{(i < H)} \sum_{(j < L)}v(x-j, y-i)h(j,i)
\label{convoDef}
\end{equation}
Dans les applications les plus courantes, $h$ est à support borné et de forme carrée et l'on parle alors de la taille du masque pour évoquer la dimension du support.
 La figure \ref{fig-ny-convo} présente les résultats de la convolution par deux masques débruiteurs \textit{moyenneurs} $h_3$ et $h_5$ de taille différentes, appliqués à l'image corrompue par un bruit gaussien : on voit la diminution des fluctuations mais aussi le flou apporté et qui rend les contours d'autant moins définis que la taille du masque est grande. La troisième image montre l'effet d'un masque gaussien $h_{g3}$. 
Les matrices définissant les masques sont les suivantes :
 
\[h_3=\frac{1}{9}\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix}, h_{25}=\frac{1}{25}\begin{bmatrix}1&1&1&1&1\\1&1&1&1&1\\1&1&1&1&1\\1&1&1&1&1\\1&1&1&1&1\end{bmatrix}, h_{dx}= \begin{bmatrix}1&2&1\\2&4&2\\1&2&1\end{bmatrix}\]  

\begin{figure}
  \centering
  \subfigure[Moyenneur 3$\times$3, PSNR=27.6dB MSSIM=0.34]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/convo/ny256_gauss25_moy3.png}}\quad  
  \subfigure[Moyenneur 5$\times$5, PSNR=27.7dB MSSIM=0.38]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/convo/ny256_gauss25_moy5.png}}\quad
  \subfigure[Filtre gaussien 3$\times$3, PSNR=27.4dB MSSIM=0.33]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/convo/ny256_gauss25_g3.png}}  
\caption{Filtrage par convolution.}
\label{fig-ny-convo}
\end{figure}

\subsubsection{Le filtre médian}
Le filtrage médian \cite{tukey77} est également une opération très employée en prétraitement pour sa simplicité et ses propriétés de préservation des contours alliées à une capacité de réduction de bruit gaussien importante. 
La valeur du niveau de gris de chaque pixel est remplacée par la médiane des niveaux de gris des pixels voisins. Un des intérêts de ce filtre réside dans le fait que la valeur filtrée est une des valeurs du voisinage, contrairement à ce qui se produit lors d'une convolution. Un autre est de bien filtrer les valeurs extrêmes et par conséquent de trouver naturellement son application dans la réduction du bruit impulsionnel.
Toutefois, la non-linéraité de cette technique et sa complexité n'en ont pas fait un filtre très utilisé jusqu'à ce que des implémentation efficaces soient proposées, en particulier le filtre à temps de calcul ``constant'' décrit par Perreault et Hebert \cite{4287006}. Il est à noter que le filtrage médian est souvent appliqué en plusieurs passes de voisinage restreint.
La figure \ref{fig-ny-median} montre la réduction de bruit impulsionnel obtenu grâce au filtre médian, dans trois conditions distinctes : median 3$\times$3 en une ou deux passes, puis médian 5$\times$5.
\begin{figure}
  \centering
  \subfigure[Médian 3$\times$3 une passe, PSNR=26.4~dB MSSIM=0.90]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/median/ny256_sap25_med3.png}}  
  \subfigure[Médian 3$\times$3 deux passes, PSNR=34.4~dB MSSIM=0.98]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/median/ny256_sap25_med3x2.png}}
  \subfigure[Médian 5$\times$5 une passe, PSNR=35.1~dB MSSIM=0.98]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/median/ny256_sap25_med5.png}}  
\caption{Réduction du bruit impulsionnel par filtre médian.}
\label{fig-ny-median}
\end{figure}


\subsubsection{Le filtre bilatéral}
Le filtre bilatéral \cite{710815} est une composition d'opérations que l'on  peut voir comme un  filtre de convolution dont les coefficients ne dépendraient pas uniquement de la position du pixel courant par rapport au pixel central, mais également de la différence de leurs intensités (cas des images en niveaux de gris). 
Si l'on note $\Omega_k$ le voisinage du pixel d'indice $k$, l'expression générale du niveau de gris estimé est donnée par 
\[\widehat{u_k}=\displaystyle\frac{\sum_{p\in \Omega_k}\left(F_S(x_p, x_k)F_I(v_p, v_k)v_p\right)}{\sum_{p\in\Omega_k }\left(F_S(x_p, x_k)F_I(v_p, v_k)\right)} \]
où $F_S$ et $F_I$ sont les fonctions de pondération spatiale et d'intensité. Classiquement, $F_S$ et $F_I$ sont des gaussiennes de moyennes nulles et d'écarts type $\sigma_S$ et $\sigma_I$.
Ce filtre se prête également bien à une utilisation en plusieurs passes sans flouter les contours. Des approximations séparables du filtre bilatéral, comme celle proposée dans \cite{1521458}, permettent d'obtenir des vitesses  d'exécution plus élevées que les versions standard. Une variante à temps de calcul constant à même été proposée en 2008 par Porikli \cite{4587843}.
Ce filtre permet un bon niveau de réduction de bruit gaussien, mais au prix d'un nombre de paramètres plus élevé à régler, ce qu'illustre la figure \ref{fig-ny-bilat} où le filtrage de la même image a été réalisé avec 9 combinaisons de $\sigma_S$ et $\sigma_I$.
\begin{figure}
  \centering
\subfigure[$\sigma_S=1.0$ et $\sigma_I=0.1$, PSNR=25.6~dB MSSIM=0.25]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_1_01.png}}
\subfigure[$\sigma_S=1.0$ et $\sigma_I=0.5$, PSNR=28.0~dB MSSIM=0.36]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_1_05.png}}
\subfigure[$\sigma_S=1.0$ et $\sigma_I=1.0$, PSNR=27.9~dB MSSIM=0.36]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_1_1.png}}\\ 
\subfigure[$\sigma_S=2.0$ et $\sigma_I=0.1$, PSNR=26.7~dB MSSIM=0.29]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_2_01.png}}
\subfigure[$\sigma_S=2.0$ et $\sigma_I=0.5$, PSNR=27.9~dB MSSIM=0.39]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_2_05.png}}
\subfigure[$\sigma_S=2.0$ et $\sigma_I=1.0$, PSNR=27.5~dB MSSIM=0.38]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_2_1.png}}\\  
\subfigure[$\sigma_S=5.0$ et $\sigma_I=0.1$, PSNR=26.8~dB MSSIM=0.29]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_5_01.png}}
\subfigure[$\sigma_S=5.0$ et $\sigma_I=0.5$, PSNR=26.8~dB MSSIM=0.37]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_5_05.png}}
\subfigure[$\sigma_S=5.0$ et $\sigma_I=1.0$, PSNR=25.9~dB MSSIM=0.36]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_5_1.png}}
\caption{Réduction de bruit gaussien par filtrage bilatéral de voisinage 5$\times$5. $\sigma_S$ et $\sigma_I$ sont les écarts type des fonctions gaussiennes de pondération spatiale et d'intensité.}
\label{fig-ny-bilat}
\end{figure}

Il existe beaucoup de  variantes d'algorithmes basés sur des moyennes ou médianes locales efféctuées sur des voisinages de formes diverses, variables et/ou adaptatives afin de sélectionner le plus finement possible les pixels pris en compte dans le calcul de la valeur filtrée. 
Le principal défaut de ces techniques est de générer des aplats dans les zones homogènes et des marches d'escalier dans les zones de transition douce (staircase effect), ces dernières pouvant être considérablement atténuées comme il a été montré dans \cite{BuadesCM06}.
L'un de ces algorithmes tend à utiliser une portion de la ligne de niveau de chaque pixel comme voisinage pour le moyennage. Cette technique a été présentée dans \cite{bertaux2004speckle} et employée pour réduire le bruit de speckle. Nous y reviendrons en détail dans le chapitre \ref{ch-lniv}.    
  

\subsubsection{Les algorithmes de filtrage par dictionnaire}
Ces algorithmes font l'hypothèse qu'il est possible de décrire l'image à débruiter en utilisant une base de fonctions permettant de décomposer l'image en une combinaison linéaire des éléments de cette base. Les bases les plus employées sont les ondelettes \cite{Mallat:2008:WTS:1525499, Daubechies:1992:TLW:130655} ainsi que les fonctions sinusoïdales (DCT \cite{1093941,strang1999discrete}). Les éléments de la base peuvent être prédéterminés ou bien calculés à partir des données de l'image, par exemple en s'appuyant sur une analyse en composantes principales ou après apprentissage \cite{elad2006image}. Le principe du débruitage est de considérer que le bruit est décorellé des fonctions de la base et donc représenté par les petits coefficients de la décomposition, que l'on peut annuler. Diverses politiques de seuillage peuvent alors être appliquées selon le type d'image et le modèle de bruit ayant chacune ses propres avantages et inconvénients. L'intérêt principal de ces méthodes est de bien restituer les transitions rapides (grande énergie), mais elles génèrent en revanche des artefacts dus aux possibles grands coefficients de bruit. 
La figure \ref{fig-ny-dwt} illustre cela en montrant le résultat du débruitage obtenu par décomposition en ondelettes et seuillage ``dur''.
Certains algorithmes récents, en particulier ceux utilisant une base d'ondelettes adaptative, comme dans \cite{elad2006image} sont proches, en terme de qualité, de l'état de l'art du domaine, avec souvent un avantage lié à des vitesses d'exécution assez rapides.

\begin{figure}
  \centering
  \subfigure[$T=20$, PSNR=26.9~dB MSSIM=0.30]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/wave/ny256_gauss25_dwt20.png}}
  \subfigure[$T=35$, PSNR=27.6~dB MSSIM=0.36]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/wave/ny256_gauss25_dwt.png}}
  \subfigure[$T=70$, PSNR=26.7~dB MSSIM=0.37]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/wave/ny256_gauss25_dwt70.png}}
\caption{Filtrage par décomposition en ondelettes et seuillage dur des coefficients inférieurs au seuil $T$.}
\label{fig-ny-dwt}
\end{figure}


\subsection{Les algorithmes de filtrage par patches}
Les techniques de réduction de bruit les plus efficaces sont aujourd'hui celles qui reposent sur les propriétés d'auto-similarité ds images, on les appelles aussi les techniques par patchs. L'idée principale est, comme pour les techniques classiques à base de de voisinage, de rechercher un ensemble de pixels pertinents et comparables afin d'en faire une moyenne. Cependant, dans le cas des techniques à patchs, la recherche de cet ensemble ne se limite pas à un voisinage du pixel central, mais fait l'hypothèse qu'il existe des zones semblables au voisinage du pixel central, réparties dans l'image et pas nécessairement immédiatement contigues.
Le moyennage s'effectue alors sur l'ensemble des ces zones identifiées.
L'algorithme des moyennes non locales (NL-means, \cite{1467423}) fut parmi les premiers de cette lignée à être proposé et bien qu'ayant représenté un progrès notable dans la qualité de débruitage, fut rapidement suivi, en particulier par le BM3D et ses variantes qui représentent actuellement l'état de l'art en terme de qualité de débruitage \cite{Dabov06imagedenoising,Dabov09bm3dimage}.  
 Les différences entre ces algorithmes résident essentiellement dans la méthode de recherche et d'identification des patchs similaires, incluant la possiblité de forme et taille variables. Une telle recherche est d'autant plus coûteuse en temps de calcul qu'elle est effectuée sur une zone étendue autour du patch central et cela représente le principal inconvénient de ces techniques qui peuvent présenter des temps d'exécution prohibitifs dans l'optique d'un traitement en temps réel.
La figure \ref{fig-ny-nlm} montre des résultats de débruitage obtenus par la méthode des NL-means avec plusieurs combinaisons des paramètres de similarité des patchs et de non localité du voisinage, notés $f$ et $t$. La figure \ref{fig-ny-bm3d} montre quant-à elle le résultat du débruitage par BM3D. Les points forts de ces deux techniques sont, comme on le voit, la qualité du débruitage avec pour l'implémentation BM3D l'avantage de ne nécessiter aucun réglage de paramètres.
\begin{figure}
  \centering
  \subfigure[$f=2$ et $t=2$, PSNR=28.5~dB MSSIM=0.37]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/nlmeans/ny256_gauss25_nlm_2_2_25.png}}\quad
  \subfigure[$f=2$ et $t=5$, PSNR=28.6~dB MSSIM=0.38]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/nlmeans/ny256_gauss25_nlm_2_5_25.png}}\quad
\subfigure[$f=5$ et $t=2$, PSNR=29.0~dB MSSIM=0.39]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/nlmeans/ny256_gauss25_nlm_5_2_25.png}}\quad
\subfigure[$f=5$ et $t=5$, PSNR=29.0~dB MSSIM=0.40]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/nlmeans/ny256_gauss25_nlm_5_5_25.png}}
\caption{Filtrage par NL-means pour différentes combinaisons des paramètres de similarité $f$ et de non localité $t$.}
\label{fig-ny-nlm}
\end{figure}
\begin{figure}
  \centering
  \includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/bm3D/ny256_gauss25_bm3D.png}
\caption{Filtrage par BM3D, PSNR=29.3~dB MSSIM=0.41}
\label{fig-ny-bm3d}
\end{figure}

\section{Les implémentations GPU des algorithmes de filtrage}
Le fabricant de processeurs graphiques Nvidia, seul type d'équipements dont nous disposons, fournit des implémentations performantes de certains prétraitements et algorithmes de filtrage. 
C'est le cas des tranformées de fourrier (FFT, DCT), qui sont par exemple utilisées dans l'implémentation d'un algorithme d'\textit{inpainting} \cite{cmla2009Kes}. 

\subsection{Le filtrage par convolution}
C'est aussi vrai pour l'opération de convolution qui a fait l'objet d'une étude et d'une optimisation poussées pour déterminer la combinaison de solutions apportant la plus grande vitesse d'exécution \cite{convolutionsoup}. L'étude a testé 16 versions distinctes, chacune présentant une optimisation particulière quant-à l'organisation de la grille de calcul, aux types de transferts entre l'hôte et le GPU ainsi qu'au types de mémoire employé pour le calcul sur le GPU. 

Les résultats montrent que l'emploi de texture comme mémoire principale pour le stockage des images à traiter apporte un gain d'environ 50\% par rapport à l'utilisation de la mémoire globale. Par ailleurs, les transactions par paquets de 128 bits apportent également une amélioration sensible, ainsi que l'emploi de la mémoire partagée comme zone de travail pour le calcul des valeurs de sortie. Le traitement de référence effectué pour les mesures est la convolution générique (non séparable) d'une image 8 bits de 2048$\times$2048 pixels par un masque de convolution de 5$\times$5 pixels, expression que l'on raccourcira déronavant en \textit{convolution 5$\times$5}.

Le meilleur résultat obtenu dans les conditions détaillées précédemment, sur architecture GT200 (carte GTX280) est de 1.4~ms pour le calcul, ce qui réalise un débit global de 945~MP/s lorsque l'on prend en compte les temps de transfert aller et retour des images (1.5~ms d'après nos mesures).
Nous continuerons d'utiliser cette mesure de débit en \textit{Pixels par seconde} pour toutes les évaluations à venir ; elle permet en particulier de fournir des valeurs de performance indépendantes de la taille des images soumises au traitement.

\subsection{Le filtre médian}
On connait peu de versions GPU du filtre médian, peut-être en raison des implémentations CPU performantes et génériques que l'on a déjà évoquées (voir par exemple \cite{4287006}) et dont le portage sur GPU ne laisse pas entrevoir de potentiel, ou bien reste à inventer. Néanmoins, une bibliothèque commerciale (LibJacket et ArrayFire) en propose une implémentation GPU dont nous avons pu mesurer les performances pour un masque de 3$\times$3 et qui est également prise comme référence par Sanchez \textit{et al.} pour évaluer les performances de leur propre implémentation appelée PCMF \cite{6288187}. 

Sur architecture GT200 (GTX260), les performances maximales de ces deux versions sont obtenues pour un masque de 3$\times$3 pixels avec respectivement 175~MP/s pour libJacket et 60~MP/s pour PCMF. 
Une précédente implémentation avait été réalisée, basée sur l'algorithme BVM décrit dans \cite{5402362}. Elle prouve son efficacité dans l'élimination des artefacts générés par les dispositifs d'imagerie médicale magnétique en 3D \cite{chen09}, mais ne permet pas d'exploiter véritablement le parallélisme des GPU en filtrage d'image en 2D.

\begin{figure}
  \centering
  \subfigure[Sur GPU GTX260. Courbe tirée de \cite{5402362}]{\label{fig-compare-jacket-pcmf1}\includegraphics[width=7cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/compar-median1.png}}\quad
  \subfigure[Sur GPU C2075. Courbe tirée de \cite{sanchez2013highly}]{\label{fig-compare-jacket-pcmf2}\includegraphics[width=7cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/compar-median2.png}}
\caption{Performances relatives des filtres médians implémentés sur GPU dans libJacket/ArrayFire, PCMF et BVM et exécutés sur deux modèle de générations différentes.}
\label{fig-compare-jacket-pcmf}
\end{figure}

La figure \ref{fig-compare-jacket-pcmf1}, tirée de \cite{5402362}, compare ces trois implémentations et montre que le débit permis par la libJacket décroit très vite avec la taille du masque pour passer à 30~MP/s dès la taille 5$\times$5, alors que le PCMF décroit linéairement jusqu'à la taille 11$\times$11 où il permet encore de traiter quelque 40~MP/s. Ceci s'explique simplement par le fait que libJacket utilise un tri simple pour la sélection de la valeur médiane alors que le PCMF exploite les propriétés des histogrammes cumulés et n'est ainsi que très peu dépendant de la taille du masque.
 
Plus récemment, Sanchez \textit{et al.} ont actualisé dans \cite{sanchez2013highly} leurs mesures sur architecture Fermi (GPU C2075) en comparant leur PCMF à la version ré-écrite en C de libJacket, nommée ArrayFire. Les courbes sont celles de la figure \ref{fig-compare-jacket-pcmf2}, où l'on constate que les variations selon la taille du masque demeurent comparables, avec toutefois des valeurs de débit augmentées, avec près de 185~MP/s pour ArrayFire et 82~MP/s pour PCMF. 

Parallèlement, on trouve aussi des implémentations de filtre médian dans des traitements plus complexes comme dans \cite{aldinucci2012parallel} où les auteurs décrivent la plus récente évolution de leur technique itérative de réduction de bruit impulsionnel, sans qu'il soit possible d'évaluer le débit du médian seul. 

Il faut noter enfin que certains codes sont plus performants sur l'ancienne architecture GT200/Tesla que sur la plus récente Fermi ; c'est le cas pour l'implémentation du médian incluse dans la bibliothèque ArrayFire et nous reviendrons sur les raisons de cette perte de performances constatée au passage à une architecture plus récente dans le chapitre consacré à notre implémentation du filtre médian.

\subsection{Le filtre bilatéral}  
Le filtre bilatéral a été plus abordé et un certain nombre de publications font état d'implémentations rapides. 
Une implémentation à temps constant en est proposée par Yang \textit{et al.} \cite{5206542} et s'exécute entre 3.7~ms et 15~ms pour une image de 1~MP. Cela ne constitue pas une référence de vitesse pour les masques de petite taille, mais devient compétitif pour des masque de grande taille (plus de 400 pixels dans le voisinage).
Une autre plus classique, employée dans la génération des images médicales tomographiques, annonce 16~ms pour un masque de 11$\times$11 sur une image de 0.25~MP.
Il demeure souvent  difficile de comparer les implémentations sans disposer des codes sources, en raison de conditions de test très variables, en particulier en ce qui concerne le modèle de GPU et la taille du masque. 
Ceci étant précisé, on peut prendre comme première référence la version proposée par Nvidia dans le SDK CUDA et nommée ``ImageDenoising''. Elle permet d'exécuter sur GPU GTX480 un filtre bilatéral 7$\times$7 sur une image, déjà en mémoire GPU, de 1~MPixels en 0.411~ms, pour un débit global de 133~MP/s.

Dans \cite{zheng2011performance}, les auteurs présentent un cadre général pour optimiser l'accès aux données par les différents kernels en utilisant la mémoire partagée pour les threads d'un même bloc. 
Le principe est de pré-charger les valeurs utiles au bloc de threads dans la mémoire partagée, cela comprend les valeurs (niveaux de gris) des pixels associés aux threads ainsi que le halo correspondant aux voisinages des pixels de la bande périphérique. On appelle communément cet ensemble la \textit{region of interest} ou ROI. La figure \ref{fig-prefetch-zheng} illustre la mise en \oe uvre de cette technique en montrant comment les threads d'un bloc se répartissent les pré-chargements en mémoire partagée des valeurs des pixels de la ROI. La géométrie des blocs de threads est ici choisie carrée, mais elle s'applique aisément à d'autres proportions comme nous le verrons plus loin. Les limites de cette méthode sont 
\begin{itemize}
\item la taille de la mémoire partagée qui doit pouvoir stocker l'ensemble des valeurs des pixels de la ROI, ce qui peut imposer une limite sur la taille des blocs de threads.
\item l'étendue du voisinage qui ne peut être pré-chargé de cette façon (4 pixels par thread) que si la surface de la ROI demeure inférieure à 4 fois le nombre de thread par bloc. 
\end{itemize}

\begin{figure}
  \centering
  \includegraphics[width=10cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/shmem_prefetch_zheng2011.png}
\caption{Illustration pré-chargement en mémoire partagée mise en \oe uvre dans \cite{zheng2011performance} pour l'implémentation, entre autres, du filtre bilatéral. a) en vert le bloc de threads associé aux pixels centraux. b-e) les blocs de pixels successivement pré-chargés en mémoire partagée. f) la configuration finale de la ROI en mémoire partagée.}
\label{fig-prefetch-zheng}
\end{figure}

Cette recette est ensuite appliquée dans l'implémentation d'un filtre bilatéral et d'un filtre à moyennes non locales (NL-means). Concernant le filtre bilatéral, ils pré-calculent aussi les coefficients de la pondération spatiale, alors que ceux de la pondération d'intensité resent calculés à la volée.
Ces deux optimisations permettent un gain de 20\% sur le temps de calcul du filtre bilatéral pour arriver à 0.326~ms dans les mêmes conditions que ci-dessus. Toutefois, le débit global ne gagne que très peu (132~MP/s) en raison de la prépondérance des temps de tranfert annoncés à 7.5~ms pour l'image de 1~MP.

Ce travail d'optimisation ne perd toutefois pas son intérêt, en ce sens où si le filtre fait partie d'une chaîne de traitement entièrement exécutée par le GPU, le transfert des données n'a besoin d'être effectué qu'une seule fois en tout début et en toute fin de traitement.  

Enfin, l'implémentation qui semble à ce jour la plus performante s'attache à réduire les redondances de calculs et parvient à filtrer une image de 9~MP avec un masque de 21$\times$21 en seulement 200~ms, soir un débit de 47~MP/s hors transfers.

\subsection{Les filtres par patches}  
Intuitivement, les algorithmes à base de patches paraissent moins adaptés au parallèlisme des GPU, du fait de la nécessité d'accéder à un voisinage étendu autour de chaque pixel. On recense malgré tout quelques implémentations dont celle présente dans le SDK CUDA qui fait cependant l'hypothèse que les coefficients de pondération spatiale sont localement constants.   
Dans \cite{PALHANOXAVIERDEFONTES}, le modèle de bruit employé vise une adaptation aux images échographiques présentant du bruit proche du speckle. Dans cette implémentation, aucune approximation des coefficients n'est faite, mais la taille maximale du patch est limitée par la quantité de mémoire partagée disponible pour chaque bloc de threads. 
Une version plus récente implémente exactement l'algorithme original \cite{nlmeansgpubelge} en proposant des optimisations algorithmiques exploitant la symétrie des coefficients spatiaux ainsi que l'interprétation du calcul de la similarité comme une convolution séparable, opération aisément parallélisable sur GPU, comme nous le détaillerons plus loin. Les auteurs parviennent ainsi à filtrer des séquences vidéo couleur de dimension 720$\times$480 à plus de 30~fps en améliorant le PSNR de 16~dB (la séquence bruitée présentant un PSNR de 20~dB). 



\section{Les techniques de segmentation}
La segmentation représente également un enjeu important dans le domaine du traitement d'image et à ce titre a 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 diagnostique médical, l'interaction Homme-machine, la discrimination d'arrière plan, etc.

On pourrait donner de la segmentation une définition spécifique par type d'usage, mais dans un souci d'unification, on propose la formulation générique suivante :
``La segmentation consiste à distinguer les zones homogènes au sein d'une image''.
Dans cette définition, 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.

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.

Les éléments constitutifs de la segmentation sont soit des régions, soit des contours. Les deux notions sont complémentaires étant donné que les contours délimitent des régions, mais les techniques de calcul basés sur l'un ou l'autre de ces éléments relèvent d'abords différents.

Les algorithmes de segmentation orientés régions 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 psychologie du \textit{gestalt} \cite{humphrey1924psychology} où l'on considère que la perception conceptuelle s'élabore au travers de regroupements visuel d'éléments.

Généralement, 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.   

\subsection{Analyse d'histogramme}\label{sec-histo}
Les techniques les plus simples à mettre en \oe uvre en segmentation sont les techniques de seuillage, basées sur une analyse de l'histogramme des niveaux de gris (ou de couleurs) et cherchant à en distinguer les différentes classes comme autant d'occurrences représentant des \textit{régions} homogènes.
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. 

Malgré la multitude de variantes proposées, ces méthodes demeurent toutefois 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 application de contrôle industriel où la simplicité algorithmique permet de surcroît des implémentations très rapides, voire câblées.

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. 
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 du cochon (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}. 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 le segment de la cible comporte des discontinuités et dans le cas de l'image bruitée,  des pixels orphelins épars demeurent en quantité. Cette technique nécessiterait une étape supplémentaire pour disposer d'une segmentation pertinente.

\begin{figure}
  \centering
  \subfigure[Image initiale comportant deux zones : le fond et le cochon (la cible)]{\label{fig-histo-cochon-a} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/cochon256.png}}\quad
  \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
  \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}}\\
\subfigure[Image initiale bruitée]{\label{fig-histo-cochon-d} \includegraphics[height=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/cochon256-sig25.png}}\quad
  \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
  \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}}
  \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.}
\label{fig-histo-cochon}
\end{figure}
 
\begin{algorithm}
  %\SetNlSty{textbf}{}{:}
  %\SetKwComment{Videcomment}{}{}
\caption{Calcul du seuil de séparation des segments de l'histogramme.}   
\label{algo-histo-cochon}
$\overline{h} \leftarrow $ histogramme sur l'image \;
$S_{init} \leftarrow 128$ \;
$S_k \leftarrow S_{init}$ \;
$\epsilon \leftarrow 1$ \;
\Repeat{$\|S_k - \frac{1}{2}(\mu_{inf} + \mu_{sup})\| < \epsilon $}{
  $\mu_{inf}=\displaystyle \frac{\displaystyle\sum_{i<S_k}h_ii}{\displaystyle\sum_{i<S_k}h_i}$ \;
  $\mu_{sup}=\displaystyle \frac{\displaystyle\sum_{i\geq S_k}h_ii}{\displaystyle\sum_{i\geq S_k}h_i}$ \;
  $S_k = \frac{1}{2}(\mu_{inf} + \mu_{sup})$ \ ;
} 
\end{algorithm}

\subsection{Partitionnement de graphe}
Un autre formalisme qui a généré une vaste classe d'algorithmes de segmentation est celui des graphes et repose sur l'idée que les régions de l'image sont représentées par les n\oe uds du graphe, alors que les liens traduisent les relations de voisinage existant entre les régions.
L'idée de base est 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.

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.
Nous pouvons retenir que 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.

L'extension a rapidement été faite 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.
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 prétendre à une segmentation qui conduise à une meilleure perception conceptuelle.
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.
Pour des images en niveaux de gris, l'expression générale des éléments $w_{ij}$ de la matrice de similarité $W$ est :
\[w_{ij} = 
\begin{cases}
\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$}\\
0 & \text{sinon}
\end{cases}
\]
On construit également la matrice de connectivité $D$, diagonale et dont les éléments sont :
\[d_{i} = \displaystyle\sum_jw_{ij}\]

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 
\[\left(D-W)\right)Y=\lambda DY \]
Certains algorithmes proposés plus récemment s'inscrivent dans cette veine \cite{wang2001image,wang2003image,felzenszwalb2004efficient,shi2000normalized}. Le principal point faible de ces techniques réside essentiellement dans la 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.

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 d'objets et/ou personnes dans les images naturelles, mais requiert de prédéterminer le nombre de segments à obtenir. Les images de la figure représentent les résultats obtenus avec un nombre de segments variant de 2 à 5 et montrent qu'il difficile de trouver un compromis acceptable. Enfin, les temps d'exécutions peuvent devenir très rapidement prohibitifs, même avec des implémentations plus optimisées. Pour information, 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) 
\begin{figure}
  \centering
  \subfigure[$s = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_2seg.png}}
  \subfigure[$s = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_3seg.png}}
  \subfigure[$s = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_4seg.png}}
  \subfigure[$s = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/graphe/cochon128_ncuts_5seg.png}}
  \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.}
\label{fig-graph-cochon}
\end{figure}

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. Des comparaison en sont rapportée dans \cite{boykov2004experimental,chandran2009computational}. 
Plusieurs algorithmes mettent en \oe uvre ce procédé avec de bons résultats, comme la méthode du \textit{push-relabel} \cite{cherkassky1997implementing} ou le \textit{pseudoflow} \cite{hochbaum2013simplifications} qui semble aujourd'hui le plus peformant.

\subsection{kernel-means, mean-shift et apparentés}
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. 
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 $C$ de segments $\Omega_i$ et leur contenu, de sorte à minimiser l'expression 
\[\sum_{i\in[1..C]}\sum_{x_k\in\Omega_i} \left(v_k-\mu_i\right)^2\]  
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..C]}\Omega_i=\Omega}$ 

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.
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.}
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.
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éfaut dont le principal est qu'elle ne converge pas nécessairement vers le regroupement optimal, même si on connait la ``bonne'' valeur de $k$. 
Un autre est d'être très dépendant du choix des $k$ éléments initiaux, en nombre et en position.

Toutefois, vraisemblablement du fait de sa simplicité d'implémentation et de temps d'exécution rapides, 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.
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 ``centre'' d'un segment. 
Des solutions ont aussi été apportées pour l'estimation de $k$ en employant, par exemple, un critère de vraisemblance pour choisir la meilleure valeur de $k$ dans un intervalle donné \cite{pelleg2000x}.
À titre d'illustration et de comparaison, l'image du cochon 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 $s=2$ à $s=5$. Les résultats sont reproduits à la figure \ref{fig-kmeans-cochon} et montrent encore une fois l'influence de $s$ sur la segmentation.
\begin{figure}
  \centering
  \subfigure[$s = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_2seg.png}}
  \subfigure[$s = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_3seg.png}}
  \subfigure[$s = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_4seg.png}}
  \subfigure[$s = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/kmeans/cochon128_kmeans_5seg.png}}
  \caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme \textit{k-means} pour un nombre $s$ 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.}
\label{fig-kmeans-cochon}
\end{figure}

Un algorithme initiallement proposé en 1975 par Fukunaga et Hostetler \cite{fukunaga1975estimation} permet de manière plus générique de déterminer le nombre de segments, ou modes, ainsi que les points, ou pixels, qui les composent. Il cherche pour ce faire à localiser les $k$ positions ou le gradient de densité s'annule. 
Il utilisé un voisinage pondére (ou \textit{kernel}) et détermine le centre de masse des segments en suivant itérativement le gradient de densité dans le voisinage autour de chaque élément du domaine. Lorsque l'algorithme à convergé, les $k$ segments sont identifiés et continennent chacun l'ensemble des points qui ont conduit à leur centre de masse respectif.
É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 demontré les propriétés et établi les lien avec d'autres techniques d'optimisation commme la descente/montée de gradient ou de filtrage commme le floutage.
Comaniciu et Peer ont alors étendu l'étude et proposé une application à la segmentation en utilisant l'espace colorimétrique CIELUV \cite{foley1994introduction} et montré qu'elle permettait une meilleure identification des modes de l'image \cite{comaniciu1999mean,comaniciu2002mean}.
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 nombre de segments identiques à ceux des figures précedentes (de 2 à 5). Le volume minimal admis pour un segment à été arbitrairement fixé à 100 pixels. 
\begin{figure}
  \centering
  \subfigure[$r=100 \Rightarrow s = 2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r100m100.png}}
  \subfigure[$r=50 \Rightarrow s = 3$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r50m100.png}}
\subfigure[$r=35 \Rightarrow s = 4$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r35m100.png}}
  \subfigure[$r=25 \Rightarrow s = 5$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/meanshift/cochon128_meanshift_r25m100.png}}  
  \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 $s$ 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.}
\label{fig-meanshift-cochon}
\end{figure}

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. 
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.

\subsection{Les contours actifs, ou \textit{snakes}}
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.
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 :
\begin{itemize}
\item l'énergie interne $E_{int}$ de la courbe, fonction de son allongement de sa courbure.
\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.
\end{itemize}
L'expression générique peut alors s'écrire 
\[E_{snake} = E_{int}+E_{ext}\]
où 
\[E_{int} = \displaystyle\sum_{s\in S} \frac{1}{2}\left(\alpha\left|\frac{\partial x_s}{\partial s}\right|^2
+\beta \left|\frac{\partial^2x_s}{\partial s^2}\right|\right)ds\]
et 
\[E_{ext} = \displaystyle\sum_{s\in S} -\left|\nabla\left[G_{\sigma}(x_s)\ast v_s\right]\right|^2ds\]

L'idée générale de l'algorithme du \textit{snake} est de trouver une courbe $S$ qui minimise l'énergie totale $E_{snake}$. 
Ici encore, la résolution du problème revient donc à 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, avec ici encore, 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 aussi des fonctions de l'abscisse curviligne $s$. La fonction $G_{\sigma}$ est la fonction d'attraction aux forts gradients de l'image. 

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 un \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.
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 à ``entourer'', sous peine de se verrouiller dans un minimum local. 
La sensibilité au bruit n'est pas non plus très bonne du fait de la formulation locale de l'énergie.  
Les ``concavités'' étroites ou présentant un goulot d'étranglement marqué sont par ailleurs mal délimitées.
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.
La figure \ref{fig-snake-tradi-cochon} illustre ces défauts en montrant quelques états intérmé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 snake original.
On voit que la convergence est assez rapide mais que le contour ainsi détérminé ne ``colle'' pas bien à l'objet que l'on s'attend à isoler.
\begin{figure}
  \centering
\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}}
\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}}
\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}}
\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}}   
\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. }
\label{fig-snake-tradi-cochon}
\end{figure} 

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.
Parmi les variantes élaborées qui tentent de pallier ces défauts, les plus intéressantes sont :
\begin{itemize}
\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 et est surtout employé dans des applications semi-automatiques où l'utilisateur définit au moins une position et une taille initiales pour la courbe. 
\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.
\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'une surface 3D soumise à un champ de force, tels que leur élévation soit constamment nulle. 
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 très pénalisants.
Après la formulation initiale de Osher et Sethian en 1988 \cite{osher1988fronts}, plusieurs façon de réduire le coût du calcul ont été formulées, dont les plus importantes restent les techniques dites \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 \textit{fast marching} \cite{sethian1996fast} qui s'applique dans le cas particulier d'une évolution monotone des fronts.  
\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 Bertaux \textit{et al.} 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 \cite{ChesnaudRB99,GallandBR03,GermainR01}. La section \ref{sec-contrib-snake} qui introduit notre contribution à cette technique en donnera une description détaillée. 
\end{itemize}
 
\subsection{Méthodes hybrides}
Aujourd'hui, les algorithmes de segmentation les plus performants en terme de qualité emploient des techniques qui tentent de tirer le meilleur parti de plusieurs des méthodes ``historiques'' décrites précédemment.
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}. Il compose avec la constructions d'histogrammes locaux pour générer une matrice de similitude (affinity matrix) et appliquer 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). Il utilise ensuite une technique adaptée de \textit{ligne de partage des eaux} (que l'on aurait rangée avec les mean-shift) pour regrouper les segments. 
Les résultats sont très bons et des implémentations efficaces ont dores et déjà été écrites (voir section \ref{sec_ea_gpu}). 
%TODO 
%peut-être dire deux mots sur le partage des eaux (avec kmeans et meanshift) puisqu'il est employé dans gPb

\section{Les implémentations GPU des techniques de segmentation}

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 massivememnt parallèle de ces matériels.
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 GPU.
 La natures des tissus et les formes à identifier sont extrêmement variées. Les images sont souvent très bruitées et les modèles de bruit divers selon l'instrumentation employée. Enfin, le diagnostique médical requerant 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 besoin médical spécifique.

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 de la réduction de bruit ou du calcul d'histogramme.
 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 indépendemment, pour chaque pixel par exemple, de multiples instances d'une version séquentielle classique du traitement.
Dans les deux cas, on lira ``implémentation GPU'', mais cela recouvrira des réalités et parfois aussi des niveaux de performance très différents.

\subsection{Calcul d'histogramme}
Comme il a été dit au paragraphe \ref{sec-histo}, les segmentations par analyse d'histogramme sont aujourd'hui cantonnées à des applications très particulières et leurs implémentations GPU ne font pas l'objet de recherches, d'autant que dans la pratique, ces traitements sont souvent réalisés par des circuits spécialisés ou programmables de type FPGA et qu'il serait illusoire d'espérer les concurrencer par une solution de type gpu, plus coûteuse, plus volumineuse et vraisemblablement moins robuste.

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 le GPU leur permetttant de conserver les données dans la mémoire du processeur graphique tout au long de l'exécution de la segmentation par level-set qui leur a servi de motivation \cite{lefohn2003interactive}. 

Les résultats annoncés ont été obtenus sur un GPU GeForce 7900 et font état du calcul des deux histogrammes nécessaires ( 64 classes chacun) sur une image de 256$\times$256 pixels en niveau de gris en 1.6~ms.

\subsection{Partitionnement de graphe}
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. 
La parallèlisation GPU des opérations sur les graphes n'est pas simple en raison de l'indépendance des blocs de threads. Peu de travaux font encore état d'implémentations efficaces mettant en \oe uvre ces techniques.
On ne recense que quelques propositions GPU de l'algorithme \textit{push-relabel} pour le partitionnement selon l'approche \textit{min cut/max flow} dont on ne retient que les trois remarquables détaillée ci-dessous. 

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). 

Les auteurs de \cite{4563095} remarquent qu'après un nombre réduit d'itérations, très peu de n\oe ud se voient changer de segment. En conséquence, certains blocs de traitement sont activés alors qu'ils n'ont effectivement pas de traitement à effectuer et retardent ainsi les traitements éventuels des blocs en attente. Pour réduire les effet de ce comportement, un indicateur d'activité est calculé à chaque itération et pour chaque bloc, en se basant 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 éffectué 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 terme de performance. 

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 des indicateurs de changement de segment par 32 par l'emploi d'un seul bit par lien. Cela a permis d'accélérer la convergence de l'algorithme, comme la 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ême conditions que précédemment (sur C1060).
Il faut noter aussi que sur C1060, l'implémentation décrite dans \cite{4563095} est moins performante, avec 17~fps, que sur la carte GTX280.

\begin{figure}
  \centering
  \includegraphics[width=12cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/graphcutscuda_stitch.png}
\caption{Évolution du nombre de pixels actifs pour les itération 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.}
\label{fig-graphcutscuda}
\end{figure}

\subsection{K-means, mean-shift et apparentés}
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'etiquetage 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. La mesure de performance a été faite avec la base de test KDD-Cup-99 \cite{kddcup99}, comportant 23 segments. Le temps annoncé pour l'exécution d'une seule itération sur un ensemble de 819200 éléments est de 200~ms. Toutefois, cette durée n'inclue pas la réduction ni les transferts et l'accélération revendiquée semble alors très discutable.

Dans \cite{5170921}, l'ensemble des tâches d'étiquetage et de mise à jour des centres est réalisé sur le GPU. Une étape de réorganisation des données est encore exécutée sur le CPU, mais elle s'avère moins pénalisante que la solution présentée précédemment, puisqu'elle permet de présenter au GPU des données permettant d'optimiser l'exécution parallèle de l'étape de réduction suivante (mise à jour des centres). 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 cette fois des mesures des temps d'exécution à convergence, qui atteignent la vingtaine de secondes pour le même ensemble de test.

La plus convaincante des implémentations de \textit{k-means} reste à notre connaissance celle décrite dans \cite{kmeansgpuopengl} et 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.

Des travaux à orientation non médicale mettent en \oe uvre sur GPU un algorithme de \textit{mean-shift} pour la poursuite de cibles dans des séquences vidéo \cite{li2009mean}. L'accélération otenue 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 vraisemblement aussi la performance de l'application. On peut malgré tout raisonnablement espérer qu'une telle solution présenterait des performances meilleures sur une carte de type Fermi possédant jusqu'à 48~Ko de mémoire partagée par bloc.

\textit{Quick shift}, une approximation de l'algorithme mean-shift gaussien, c'est à dire utilisant des masques de pondération gaussiens, permettant d'obtenir un résultat en une seule passe (sans itérer) et proposée initiallement dans \cite{vedaldi2008quick} a été parallélisée sur GPU par ses auteurs et décrite dans \cite{fulkerson2012really}. La recherche de performance se traduit par des approximations, en particulier on restreint les 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égligeable.
Ensuite on construit un arbre des liens entre les pixels, mais on limite la recherche à une distance maximale de $\sigma$. Par ailleurs, on diminue arbitrairement la dynamique de l'espace colométrique par 2. Enfin, la segmentation est obtenu par simple partionnnement de l'arbre selon un seuil $\tau$.
Pour s'affranchir de la relative petite taille de la mémoire partagée sans devoir pâtir de la grande latence des accès à la mémoire globale de GPU, les auteurs ont ici choisi d'associer l'image et l'estimation de densité à des textures et ainsi bénéficier du mécanisme de cache.
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, 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.
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$.

\begin{figure}
  \centering
\subfigure[Image originale]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-orig.png}}\quad
\subfigure[$\tau=10$ et $\sigma=2$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s2t10.png}}\quad
\subfigure[$\tau=10$ et $\sigma=10$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s10t10.png}}\quad
\subfigure[$\tau=20$ et $\sigma=10$]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/quick-shift-yo-s10t20.png}}\quad
\caption{Segmentation d'une image couleur de 512$\times$512 pixels par l'implémentation GPU quick-shift de \cite{fulkerson2012really}.}
\label{fig-quickshift-yo}
\end{figure}

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 reproduite à la figure  \ref{fig-meanshift-castle} afin de montrer les différences avec une implémentation standard du \textit{mean-shift}.

\begin{figure}
  \centering
\subfigure[Image originale]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/castle-meanshift.png}}\quad
\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
\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}}
\caption{Segmentation d'une image couleur de 2256$\times$3008 pixels.}
\label{fig-meanshift-castle}
\end{figure}

\subsection{Snakes et Level set}
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és par résonnance magnétique (IRM) en exploitant pour la première fois le caractère creux du système d'équations à résoudre, \textit{i.e.} variante narrow-band, contrairement à la première solution 2D présentée dans \cite{rumpf2001level} qui implémente la version standard. En ne transférant au GPU, pour chaque itération, que les petits pavés de données actifs et en les  rangeant alors de manière contigue en texture pour optimiser les accès en lecture, les auteurs sont ainsi parvenu à effectuer, pour des données volumiques de 256$\times$256$\times$175, entre 3.5 et 70 itérations par seconde, à comparer aux 50 itérations par seconde en 2D sur image de 128$^2$ pixels otenues dans \cite{rumpf2001level}. La limitation principale de cettesolution est celle des dimensions maximales admises pour une texture qui était de 2048$^2$ pour le GPU ATI Radeon 9800 pro employé (et programmé en openGL, car ni openCL ni CUDA n'étaient encore disponible à l'époque).
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$^2$ pixels.

La plus performante des implémentations à 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$^3$ pixels issues d'examen 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 au travers de la refonte du code responsable de la détermination des pavés actifs. Il 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 contigue 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 impossile à utiliser efficacement lorsque les éléments à accéder sont très irrégulièrement répartis en mémoire. 

Ce faisant, le nombre cumulé total de pavés ainsi traités lors des 2200 itérations de la segmentation der l'image d'exemple s'élève à 294 millions à comparer aux 4877 millions traités par l'algorithme \textit{narrow-band} standard. 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 \textit{narrow-band} standard. Les deux courbes sont globalement affines et se croisent pour une proportion de pavés actifs proche de 10\%.
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 ne semble pas su justifier avec l'image et l'initialisation dont les performances sont détaillées, mais qui pourrait l'être dans d'autres conditions, comme peut le suggérer 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.

\begin{figure}
  \centering
\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
\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}}
\caption{Segmentation d'images issues d'examens IRM par la méthode des level set à bande étroite.}
\label{fig-l7-narrow}
\end{figure}

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 à véritablement été implémentée de manière spécifique et efficace \cite{snakegvf06, bauer2009segmentation, li2011robust, snakegvfopencl12}. Les variantes de type géométrique, principalement en raison de l'irrégularité des motifs d'accès à la mémoire, restent à ce jour sans implémentation GPU.
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 RVBA de manière à s'affranchir du format 16 bits de la représentation : les deux premiers canaux R et V contiennent les valeursreprésentant respectivement le gradients selon $dx$ et $dy$ sous une forme codée par la valeurs des 2 autres canaux. 
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.

\begin{figure}
  \centering
  \subfigure[Contour initial]{\label{fig-epaule-init}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/snake-epaule-init.png}}\quad
  \subfigure[Contour final]{\label{fig-epaule-fin}\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/snake-epaule-fin.png}}
\caption{Segmentation d'une image d'épaule en 1024$^2$ pixels issue d'un examen IRM par l'implémentation du snake GVF de \cite{snakegvf06}. Le contour est représenté en rougeet le contour final est obtenu en 11~s. }
\label{fig-snakegvf}
\end{figure}

Les performances annoncées montrent tout d'abord que l'approximation adoptée n'a qu'un impact extrêmement limité sur le résulat 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$^2$ 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.

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.   

En adaptant sur GPU une variante dite FD-snake \cite{li2011robust} du snake GVF (pour Fourier Descriptors) permettant une convergence plus rapide et un calcul parallèle beaucoup plus adapté au GPU, Li \textit{et al.} parviennent quant à eux à suivre les déformations d'un contour en temps réel dans des images issues d'examens échographique ; Un contour de 100 points pouvant 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. 

La plus aboutie des implémentations actuelles du snake GVF est enfin celle présentée par Smistad \textit{et al.} dans \cite{snakegvfopencl12} et 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.
Les performances sont alors nettement en hausse avec des segmentations d'images médicales d'IRM de 512$^2$ 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 openGL permet d'exécuter le code sur les GPU des deux principaux fabricants.   

\subsection{Algorithmes hybrides}
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 semgentation d'objets et personnages dans des image naturelles, à été implémenté en CUDA par Catanzaro \textit{et al.} et est décrit dans \cite{5459410}. La qualité des contours extraits y est préservée et le temps de traitement y est réduit d'un facteur supérieur à 100 : 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.
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 à 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.   
La figure \ref{fig-gPb} présente quelques résultats d'extraction de contours.
\begin{figure}
  \centering
\includegraphics[height=4cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/gPb_examples.png}
\caption{Extraction de contour par la version GPU de l'algorithme gPb. Les images sont issues de la base BSDS  \cite{martin2001database}}
\label{fig-gPb}
\end{figure}

\section{Conclusion}
La présentation que nous venons de faire des principales techniques de filtrage et de ségmentation ainsi que des implémentations sur GPU qui leur ont été consacrées nous ont permis de mettre une évidence en lumière : malgré leur orientation grand public et les langages de huat 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 leur sont spécifiques et obtenir un code rapide découle nécessairement d'un compromis qui peut parfois être complexe à affiner. 

Ajoutons que les générations de GPU qui se succèdent conservent certes des caractéristiques communes mais diffèrent suffisemment quant-à la distribution des ressources, rendant toute généralité 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 consacré au filtre médian.

Cet état de fait rend les résultats publiés par les chercheurs souvent délicats à intérpré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.   

Pour aider les développeurs à allouer les ressources de manière optimale, ou tout du moins estimer le dégré 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 ont peut entrer les paramètres d'exécution d'un \textit{kernel} parallèle : nmobre 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'occupancy) qui traduit le rapport, à chaque instant, entre le nombre de warps actifs et le nombre maximal de warps par processeur (SM = Streaming Multiprocessor). 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.

Toutefois, comme l'a clairement demontré Volkov dans \cite{volkov2010better}, ce paradigme peut aisément être remis en cause et Volkov parvient effectivement à améliorer les peformances d'un certain nombre d'exemples génériques dans des conditions de faible valeur d'occupancy. 
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 devienne sans intérêt.
Les chapitres suivants présentant nos contributions reviendront sur ces aspects en proposant des solutions pour accroître la performance des algorithmes parallélisés.