X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/these_gilles.git/blobdiff_plain/5d5e214d54eab1411ef86a32e11fccfaa53dbc45..10d54068846e7aee58e98dc76fa92f6f3a5c957a:/THESE/Chapters/chapter2/chapter2.tex?ds=sidebyside diff --git a/THESE/Chapters/chapter2/chapter2.tex b/THESE/Chapters/chapter2/chapter2.tex index ebdbd55..43ceac1 100644 --- a/THESE/Chapters/chapter2/chapter2.tex +++ b/THESE/Chapters/chapter2/chapter2.tex @@ -1,687 +1,13 @@ -L'étendue des techniques applicables aux images numériques est aujourd'hui si vaste qu'il serait illusoire de chercher à les décrire ici de manière exhaustive. Ce chapitre présente plus spécifiquement les algorithmes utilisés en présence d'images (fortement) bruitées. +L'étendue des techniques applicables aux images numériques est aujourd'hui si vaste qu'il serait illusoire de chercher à les décrire ici de manière exhaustive. Ce chapitre présente plus spécifiquement les algorithmes utilisés en présence d'images (fortement) bruitées, c'est-à-dire présentant une altération de la réalité \og absolue \fg{} qu'elles représentent. Le bruit rend potentiellement délicate l'extraction des informations utiles contenues dans les images perturbées ou en complique l'interprétation, automatisée ou humaine. L'intuition incite donc à chercher des méthodes efficaces de pré-traitement réduisant la puissance du bruit et permettant ainsi aux traitements de plus haut niveau (comme la segmentation), d'opérer dans de meilleures conditions. -Toutefois, il faut également considérer que les opérations préalables de réduction de bruit génèrent des modifications statistiques et peuvent altérer 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 en préserver toute l'information. +Toutefois, il faut également considérer que les opérations préalables de réduction de bruit génèrent des modifications statistiques et peuvent altérer 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 en préserver toute l'information, ce qui est le cadre de notre contribution portant sur un algorithme de segmentation par contour actif polygonal (dit \textit{snake}, voir chapitre \ref{ch-snake}). De plus, toute opération supplémentaire si basique soit elle, réduit le temps de traitement disponible pour l'opération de haut niveau. En effet, 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. -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 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 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. +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}, ou de l'imagerie médicale à ultrasons (échographie). +Ces dispositifs d'acquisition sont, par essence, générateurs de bruits divers, inhérents aux technologies mises en \oe uvre 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 d'ores 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. Quant à la recherche d'un filtre universel, bien qu'encore illusoire, elle 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 de ces perturbations peut être représenté par un seul bruit blanc gaussien, de type \textit{additif} (AWGN, Additive White Gaussian Noise), 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, laquelle peut être réduite, pour une scène donnée, par moyennage de plusieurs observations, ou vues. -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. 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 -\[ 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}\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 -\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} - -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. - - -\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 ``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 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. -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 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. -Le moyennage s'effectue alors sur l'ensemble de ces zones identifiées. -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 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 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. -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. -\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 sur 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 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. - -\subsection{Le filtrage par convolution} -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 entre l'hôte et le GPU ainsi qu'aux 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 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 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}\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}. - -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èles 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 \ref{ch-median} consacré à notre implémentation du filtre médian. - -\subsection{Le filtre bilatéral}\label{sec-bilateral} - -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 \textit{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 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, ils 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. Toutefois, le débit global ne gagne que très peu (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. - -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. - -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 diagnostic 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 des zones homogènes au sein d'une image'', -où 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ées 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 théorie 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 peu robustes et présentent l'inconvénient majeur de ne pas garantir la connexité des régions déterminées. On les réserve à des applications très spécifiques où, par exemple, on dispose d'une image de référence dont l'histogramme peut être comparé à celui des images à traiter. C'est le cas de certaines applications de contrôle industriel où la simplicité algorithmique permet de surcroît des implémentations très rapides, voire câblées. - -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 : les régions identifiées présentent des discontinuités et dans le cas de l'image bruitée, quantité de pixels orphelins subsistent en quantité. Cette technique nécessiterait une étape supplémentaire pour obtenir une segmentation pertinente. - -\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