1 \section{Les techniques de réduction de bruit}
2 La très grande majorité des algorithmes de réduction de bruit fait l'hypothèse que la perturbation est de type additif 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}.
3 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. Aussi nous focaliserons nous 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.
5 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.
6 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 est défini par
7 \[ PSNR = 10log_{10}\left(\frac{D^2}{\displaystyle\frac{1}{N}\sum_{k < N}\left(v_k - u_k\right)^2}\right)\]
8 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.
10 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 fournir 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.
14 \subfigure[Sans bruit]{\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/ny256.png}}
15 \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}}
16 \subfigure[Bruit impulsionnel 25\%, PSNR=9.48~dB MSSIM=0.04]{\label{ny-sap}\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/ny256_sap25.png}}
17 \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\%.}
21 \subsection{Les opérateurs de base}\label{sec-op-base}
22 \subsubsection{Le filtre de convolution}
23 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.
24 Si la fonction définissant le masque de convolution est notée $h$, l'expression générale de la valeur de sortie estimée au pixel de coordonnées $(i,j)$ est donnée par
26 \widehat{u}(x, y) = \left(\bar{v} * h\right) = \sum_{(i < H)} \sum_{(j < L)}v(x-j, y-i)h(j,i)
29 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.
30 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}$.
31 Les matrices définissant les masques sont les suivantes :
33 \[h_3=\frac{1}{9}\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix}, h_{5}=\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_{g3}= \frac{1}{16}\begin{bmatrix}1&2&1\\2&4&2\\1&2&1\end{bmatrix}\]
37 \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
38 \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
39 \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}}
40 \caption{Filtrage par convolution.}
44 Lorsque la matrice $h$ définissant le masque peut s'écrire comme le produit de deux vecteurs à une dimension $h_v$ et $h_h$, on parle alors de convolution séparable, qui peut être effectuée en deux étapes distinctes de convolution à une dimension, l'une verticale, l'autre horizontale.
47 \subsubsection{Le filtre médian}
48 Le filtrage médian \cite{tukey77} est également une opération très employée en pré-traitement pour sa simplicité coneptuelle et ses propriétés de préservation des contours alliées à une capacité de réduction du bruit gaussien importante.
49 La valeur du niveau de gris de chaque pixel est remplacée par la médiane des niveaux de gris des pixels dans un voisinage donné. 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.
50 Le filtre médian est aussi un estimateur robuste de la valeur moyenne permettant de bien s'affranchir des perturbations présentant, dans le voisinage, des valeurs très éloignées de la moyenne. Par conséquent, il trouve aussi son application dans la réduction du bruit impulsionnel.
51 Toutefois, la non-régularité et l'importance du temps de calcul de cette technique 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 \og constant \fg{} 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, car cela permet d'obtenir des niveaux de débruitage semblables à ceux permis par des voisinages plus grands, mais avec des temps de calcul réduits.
52 La figure \ref{fig-ny-median} montre la réduction de bruit impulsionnel obtenu sur l'image \ref{ny-sap} grâce au filtre médian, dans trois conditions distinctes : médian 3$\times$3 en une ou deux passes, puis médian 5$\times$5.
55 \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}}
56 \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}}
57 \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}}
58 \caption{Réduction du bruit impulsionnel par filtre médian.}
63 \subsubsection{Le filtre bilatéral}
64 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 dépendent de la position du pixel considéré par rapport au pixel central, mais dépendent également de la différence entre leurs deux niveaux de gris (cas des images en niveaux de gris).
65 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
66 \[\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)} \]
67 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$.
68 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 standards. Une variante à temps de calcul constant a même été proposée en 2008 par Porikli \cite{4587843}.
69 Ce filtre atteint 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$. On y constate que la qualité s'obtient par compromis entre la valeur de $\sigma_S$ et celle de $\sigma_I$ et que le \og meilleur \fg{} résultat parmi nos 9 combinaisons est à la figure \ref{subfig-bruit-bilat}.
72 \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}}
73 \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}}
74 \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}}\\
75 \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}}
76 \subfigure[$\sigma_S=2.0$ et $\sigma_I=0.5$, PSNR=27.9~dB MSSIM=0.39]{\label{subfig-bruit-bilat}\includegraphics[width=4cm]{/home/zulu/Documents/these_gilles/THESE/codes/bilat/ny_gauss25_bilat_2_05.png}}
77 \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}}\\
78 \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}}
79 \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}}
80 \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}}
81 \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é.}
85 Il existe beaucoup de variantes d'algorithmes basés sur des moyennes ou médianes locales effectué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.
86 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}.
88 L'un de ces algorithmes cherche dans l'image bruitée, une portion de la ligne de niveau de chaque pixel et l'utilise ensuite 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}.
91 \subsubsection{Les algorithmes de filtrage par dictionnaire}
92 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}. On considère que le bruit, blanc, est décorellé des fonctions de la base, son énergie est ainsi répartie uniformément sur tout l'espace de décomposition, contrairement au signal informatif. Le principe du filtrage est alors de supprimer les composantes à faible énergie, représentées par les petits coefficients de la décomposition.
94 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 peuvent en revanche générer des artefacts, par exemple lorsque la PDF du bruit partage avec le signal une composante d'énergie élevée.
95 La figure \ref{fig-ny-dwt} illustre cela en montrant le résultat du débruitage obtenu par décomposition en ondelettes et seuillage \og dur \fg{} : lorsque la valeur du seuil croît, des aplats apparaissent et sont visuellement génants bien que les valeurs des indices de qualité n'aient pas diminué significativement.
96 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.
100 \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}}
101 \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}}
102 \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}}
103 \caption{Filtrage par décomposition en ondelettes et seuillage dur des coefficients inférieurs au seuil $T$.}
108 \subsection{Les algorithmes de filtrage par patches}
109 Les techniques de réduction de bruit les plus efficaces sont aujourd'hui celles qui reposent sur les propriétés d'auto-similarité des images, on les appelle aussi les techniques par patches. L'idée principale est, comme pour les techniques classiques à base de voisinage, de rechercher un ensemble de pixels pertinents et comparables afin d'en faire une moyenne. Cependant, dans le cas des techniques à patches, 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 contiguës.
110 Le moyennage s'effectue alors sur l'ensemble de ces zones identifiées.
111 L'algorithme des moyennes non locales (NL-means, \cite{1467423}), parmi les premiers de cette lignée à être proposé, 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 ce qui se fait de mieux en terme de qualité de débruitage \cite{Dabov06imagedenoising,Dabov09bm3dimage}.
113 Les différences entre ces algorithmes résident essentiellement dans la méthode de recherche et d'identification des patches 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.
114 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 patches 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.
117 \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
118 \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
119 \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
120 \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}}
121 \caption{Filtrage par NL-means pour différentes combinaisons des paramètres de similarité $f$ et de non localité $t$.}
126 \includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/bm3D/ny256_gauss25_bm3D.png}
127 \caption{Filtrage par BM3D, PSNR=29.3~dB MSSIM=0.41}
131 \section{Les implémentations sur GPU des algorithmes de filtrage\label{sec-filtresgpu}}
132 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.
133 C'est le cas des transformées de fourrier (FFT, DCT), qui sont par exemple utilisées dans l'implémentation d'un algorithme d'\textit{inpainting} \cite{cmla2009Kes} et de l'opération de convolution dont l'étude est présentée ci-dessous.
135 \subsection{Le filtrage par convolution}
136 L'opération de convolution a fait l'objet d'une étude et d'une optimisation poussées par le fabricant Nvidia 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 de données entre l'hôte et le GPU ainsi qu'aux types de mémoire employés pour le calcul sur le GPU.
138 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 dorénavant en \textit{convolution 5$\times$5}.
140 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 représente 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).
141 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épendamment de la taille des images soumises au traitement\footnote{Dans la pratique, le débit dépend peu de la taille de l'image, pour les images suffisamment grandes que nous traitons. Le débit croit légèrement avec la taille de l'image en raison des transferts plus efficaces.}.
143 \subsection{Le filtre médian}\label{sec-median}
144 On connaît 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}.
146 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 (transferts de données compris).
147 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.
151 \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
152 \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}}
153 \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èles de générations différentes.}
154 \label{fig-compare-jacket-pcmf}
157 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.
159 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.
161 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.
163 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 \ref{ch-median}.
165 \subsection{Le filtre bilatéral}\label{sec-bilateral}
166 Le filtre bilatéral a été plus étudié et un certain nombre de publications font état d'implémentations rapides.
167 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 masques de grande taille (plus de 400 pixels dans le voisinage).
168 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.
169 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.
170 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 \og ImageDenoising \fg. 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~MP en 0,411~ms, pour un débit global de 133~MP/s.
172 Dans \cite{zheng2011performance}, les auteurs présentent un cadre général pour optimiser l'accès aux données par les différents \textit{kernels} en utilisant la mémoire partagée pour les threads d'un même bloc.
173 Le principe est de pré-charger les valeurs utiles au bloc de threads dans la mémoire partagée, c'est-à-dire 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. Dans cet exemple, chaque thread (en vert) gère le chargement de 4 pixels (blocs en gris) en mémoire partagée. On parvient ainsi à éviter les conflits de banques de mémoire partagée et à répartir la charge des threads de sorte à générer le moins de perte de performances. 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 :
175 \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.
176 \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 threads par bloc.
181 \includegraphics[width=12cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/shmem_prefetch_zheng2011.png}
182 \caption{Illustration du pré-chargement en mémoire partagée mis 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.}
183 \label{fig-prefetch-zheng}
186 Cette méthode empirique de pré-chargement 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, les auteurs pré-calculent aussi les coefficients de la pondération spatiale, alors que ceux de la pondération d'intensité restent calculés à la volée.
187 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 (bilatéral 7$\times$7 sur image 1~MP). Toutefois, le débit global ne progresse pas (132~MP/s) en raison de la prépondérance des temps de transfert annoncés à 7,5~ms pour l'image de 1~MP, alors qu'ils n'étaient que de 7,1~ms dans les conditions décrites par Nvidia.
189 Ce travail d'optimisation ne perd toutefois pas son intérêt dans la mesure 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.
191 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 transferts.
193 \subsection{Les filtres par patches}
194 Intuitivement, les algorithmes à base de patches paraissent moins adaptés au parallélisme des GPUs, du fait de la nécessité d'accéder à un voisinage étendu autour de chaque pixel. On en recense tout de même quelques implémentations, comme celle présente dans le SDK CUDA qui fait l'hypothèse que les coefficients de pondération spatiale sont localement constants.
195 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.
196 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).
