\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}.
+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}.
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.
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
+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
\[ 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.
+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.
\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]{\label{ny-sap}\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\%.}
+ \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}\label{sec-op-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
+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
\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}
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_{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}= \begin{bmatrix}1&2&1\\2&4&2\\1&2&1\end{bmatrix}\]
+\[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}\]
\begin{figure}
\centering
\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 du 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 \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.
+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.
+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.
+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.
+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.
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.
\begin{figure}
\centering
\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).
+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).
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$.
\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.
+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.
+
+ 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.
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.
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.
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}.
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).
-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.
+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.}.
\subsection{Le filtre médian}\label{sec-median}
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}.
\begin{figure}
\centering
- \includegraphics[width=10cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/shmem_prefetch_zheng2011.png}
+ \includegraphics[width=12cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/shmem_prefetch_zheng2011.png}
\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.}
\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, 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.
+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.
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.
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.
