11 oct

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

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

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

-%TODO
-% dire aussi que le prétraitement, ça prend du temps. C'est évident mais c'est mieux en le disant
- 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 isues de l'imagerie radar à ouverture synthétique (ROS ou en anglais SAR) \cite{cutrona1990synthetic}, de l'imagerie médicale à ultrasons (echographie) ou encore biologique dans le cas de la microscopie électronique. 
-Ces dispositifs d'acquisition sont naturellement, et par essence, générateurs de bruits divers, inhérents aux thechnologies mises en \oe uvre au sein de ces systèmes et qui viennent dégrader l'image idéale de la scène que l'on cherche à représenter ou analyser. On sait aujourd'hui caractériser de manière assez précise ces bruits et la section \ref{sec_bruits} en détaille les  origines physiques ainsi que  les propriétés statistiques qui en découlent.
+
+Enfin, toute opération aussi basique soit elle, prend un certain temps qui réduit potentiellement celui disponible pour le traitement de haut niveau. Lorsque les images à analyser sont de grande taille,  procéder à un débruitage préalable peut s'avérer incompatible avec les contraintes de débit et les temps requis par le traitement de haut niveau.
+
+ Les images auxquelles nous nous intéressons sont généralement les images numériques allant des images naturelles telles que définies par Caselles \cite{Caselles99topographicmaps} aux images d'amplitude issues de l'imagerie radar à ouverture synthétique (ROS ou en anglais SAR) \cite{cutrona1990synthetic}, de l'imagerie médicale à ultrasons (échographie) ou encore biologique dans le cas de la microscopie électronique. 
+Ces dispositifs d'acquisition sont naturellement, et par essence, générateurs de bruits divers, inhérents aux technologies mises en \oe uvre au sein de ces systèmes et qui viennent dégrader l'image idéale de la scène que l'on cherche à représenter ou analyser. On sait aujourd'hui caractériser de manière assez précise ces bruits et la section \ref{sec_bruits} en détaille les  origines physiques ainsi que  les propriétés statistiques qui en découlent.

 On peut dores et déjà avancer que la connaissance de l'origine d'une image et donc des propriétés des bruits associés qui en corrompent les informations, est un atout permettant de concevoir des techniques de filtrage adaptées à chaque situation. Toutefois, la recherche d'un filtre universel, bien qu'encore illusoire, n'est pas abandonnée, tant les besoins sont nombreux, divers et souvent complexes.    
        
 \section{Modèle d'image bruitée}
 On peut dores et déjà avancer que la connaissance de l'origine d'une image et donc des propriétés des bruits associés qui en corrompent les informations, est un atout permettant de concevoir des techniques de filtrage adaptées à chaque situation. Toutefois, la recherche d'un filtre universel, bien qu'encore illusoire, n'est pas abandonnée, tant les besoins sont nombreux, divers et souvent complexes.    
        
 \section{Modèle d'image bruitée}


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

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

 
 La figure \ref{fig-ny-noises} montre une image de synthèse issue de la base de test COIL \cite{coil}, supposée sans bruit et qui sera considérée comme référence, ainsi que deux versions bruitées, respectivement avec un bruit gaussien d'écart type 25 et un bruit impulsionnel affectant 25\% des pixels. 
 L'indice de qualité le plus employé pour mesurer la similarité entre deux images est le PSNR (pour Peak Signal to Noise Ratio). Il est exprimé en décibels (dB) et se calcule en appliquant la formule  
 
 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  


@@ -202,9 +203,17 @@ On connait peu de versions GPU du filtre médian, peut-être en raison des impl
 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.
 
 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.
 

-La figure \ref{fig-compare-jacket-pcmf}, 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.
+\begin{figure}
+  \centering
+  \subfigure[Sur GPU GTX260. Courbe tirée de \cite{5402362}]{\label{fig-compare-jacket-pcmf1}\includegraphics[width=7cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/compar-median1.png}}\quad
+  \subfigure[Sur GPU C2075. Courbe tirée de \cite{sanchez2013highly}]{\label{fig-compare-jacket-pcmf2}\includegraphics[width=7cm]{/home/zulu/Documents/these_gilles/THESE/Chapters/chapter2/img/compar-median2.png}}
+\caption{Performances relatives des filtres médians implémentés sur GPU dans libJacket/ArrayFire, PCMF et BVM et exécutés sur deux modèle de générations différentes.}
+\label{fig-compare-jacket-pcmf}
+\end{figure}
+
+La figure \ref{fig-compare-jacket-pcmf1}, tirée de \cite{5402362}, compare ces trois implémentations et montre que le débit permis par la libJacket décroit très vite avec la taille du masque pour passer à 30~MP/s dès la taille 5$\times$5, alors que le PCMF décroit linéairement jusqu'à la taille 11$\times$11 où il permet encore de traiter quelque 40~MP/s. Ceci s'explique simplement par le fait que libJacket utilise un tri simple pour la sélection de la valeur médiane alors que le PCMF exploite les propriétés des histogrammes cumulés et n'est ainsi que très peu dépendant de la taille du masque.

  
  

-Plus récemment, Sanchez \textit{et al.} ont actualisé 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-arrayfire-pcmf}, 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. 
+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. 
 
 
 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. 
 


@@ -410,7 +419,7 @@ On voit que la convergence est assez rapide mais que le contour ainsi détérmin
 \subfigure[L'état  du contour après la septième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it7.png}}
 \subfigure[L'état du contour après la dixième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it10.png}}
 \subfigure[L'état du contour après la centième itération. C'est le contour final.]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_result.png}}   
 \subfigure[L'état  du contour après la septième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it7.png}}
 \subfigure[L'état du contour après la dixième itération]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_it10.png}}
 \subfigure[L'état du contour après la centième itération. C'est le contour final.]{\includegraphics[width=3cm]{/home/zulu/Documents/these_gilles/THESE/codes/snake/cochon128_tradi_snake_result.png}}   

-\caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme dit du \textit{snake}, dans sa version originale. Les paramètres d'élastictié, de raideur et d'attraction ont été fixés respectivement aux valeurs 5, 0.1 et 5. }
+\caption{Segmentation d'une image en niveaux de gris de 128 $\times$ 128 pixels par algorithme dit du \textit{snake}, dans sa version originale. Les paramètres d'élasticité, de raideur et d'attraction ont été fixés respectivement aux valeurs 5, 0.1 et 5. }

 \label{fig-snake-tradi-cochon}
 \end{figure} 
 
 \label{fig-snake-tradi-cochon}
 \end{figure} 
 


@@ -525,6 +534,14 @@ Les algorithmes de type \textit{snake}, très coûteux en temps de calcul, pouva
 Parmi les premières solutions décrites, \cite{snakegvf06} propose une implémentation réalisée en openGL, où les données de gradient sont compactées en texture RVBA de manière à s'affranchir du format 16 bits de la représentation : les deux premiers canaux R et V contiennent les valeursreprésentant respectivement le gradients selon $dx$ et $dy$ sous une forme codée par la valeurs des 2 autres canaux. 
 Par ailleurs, une approximation du système linéaire à résoudre est proposée afin de donner une structure bande symétrique à la matrice à inverser, ce qui améliore considérablement l'efficacité des accès aux données au travers du cache.
 
 Parmi les premières solutions décrites, \cite{snakegvf06} propose une implémentation réalisée en openGL, où les données de gradient sont compactées en texture RVBA de manière à s'affranchir du format 16 bits de la représentation : les deux premiers canaux R et V contiennent les valeursreprésentant respectivement le gradients selon $dx$ et $dy$ sous une forme codée par la valeurs des 2 autres canaux. 
 Par ailleurs, une approximation du système linéaire à résoudre est proposée afin de donner une structure bande symétrique à la matrice à inverser, ce qui améliore considérablement l'efficacité des accès aux données au travers du cache.
 

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

 Les performances annoncées montrent tout d'abord que l'approximation adoptée n'a qu'un impact extrêmement limité sur le résulat de la segmentation avec un écart radial maximal inférieur à 1.3 pixel par rapport au calcul exact effectué sur CPU. Enfin, la segmentation de l'image d'exemple en 1024$^2$ pixels s'effectue en un total de 11~s après l'initialisation manuelle reproduite à la figure \ref{fig-snakegvf}. Cela est annoncé comme presque 30 fois plus rapide que l'implémentation CPU de référence, mais demeure beaucoup trop lent pour un usage interactif.
 
 Une solution directe employant la transformée de fourier pour inverser le système à résoudre a été décrite récemment dans  \cite{zheng2012fast}et programmée en employant la bibliothèque openGL. Les exemples fournis montrent des objets segmentés dans des images d'environ 10000 pixels en une durée de l'ordre de la demi seconde.   
 Les performances annoncées montrent tout d'abord que l'approximation adoptée n'a qu'un impact extrêmement limité sur le résulat de la segmentation avec un écart radial maximal inférieur à 1.3 pixel par rapport au calcul exact effectué sur CPU. Enfin, la segmentation de l'image d'exemple en 1024$^2$ pixels s'effectue en un total de 11~s après l'initialisation manuelle reproduite à la figure \ref{fig-snakegvf}. Cela est annoncé comme presque 30 fois plus rapide que l'implémentation CPU de référence, mais demeure beaucoup trop lent pour un usage interactif.
 
 Une solution directe employant la transformée de fourier pour inverser le système à résoudre a été décrite récemment dans  \cite{zheng2012fast}et programmée en employant la bibliothèque openGL. Les exemples fournis montrent des objets segmentés dans des images d'environ 10000 pixels en une durée de l'ordre de la demi seconde.   


@@ -545,39 +562,18 @@ La figure \ref{fig-gPb} présente quelques résultats d'extraction de contours.
 \label{fig-gPb}
 \end{figure}
 
 \label{fig-gPb}
 \end{figure}
 

-     
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+\section{Conclusion}
+La présentation que nous venons de faire des principales techniques de filtrage et de ségmentation ainsi que des implémentations sur GPU qui leur ont été consacrées nous ont permis de mettre une évidence en lumière : malgré leur orientation grand public et les langages de huat niveau permettant d'accéder rapidement à la programmation GPU, la parallélisation efficace d'un algorithme séquentiel destiné à s'exécuter sur ces processeurs n'est pas triviale. Le modèle et les contraintes de programmation leur sont spécifiques et obtenir un code rapide découle nécessairement d'un compromis qui peut parfois être complexe à affiner. 

 
 

+Ajoutons que les générations de GPU qui se succèdent conservent certes des caractéristiques communes mais diffèrent suffisemment quant-à la distribution des ressources, rendant toute généralité vaine et faisant qu'un code optimisé pour un modèle donné peut devenir moins rapide avec un modèle plus récent. Prenons l'exemple du nombre maximal de registres utilisables par thread ; il est de 128 sur GPU C1060 contre seulement de 63 pour un C2070. Un code faisant un usage optimisé des registres sur C1060 pourra s'exécuter plus lentement sur C2070. C'est un cas de figure sur lequel nous reviendrons plus en détail dans le chapitre consacré au filtre médian.

 
 

+Cet état de fait rend les résultats publiés par les chercheurs souvent délicats à intérpréter et plus encore à reproduire lorsque l'on souhaite comparer les performances de nos propres codes avec les références du moment, sauf à disposer d'un panel de cartes GPU représentant toutes les évolutions de l'architecture et ce pour au moins les deux grands fabricants de GPUs que sont ATI et Nvidia.   

 
 

+Pour aider les développeurs à allouer les ressources de manière optimale, ou tout du moins estimer le dégré d'optimisation de leur code à l'aune de la vitesse d'exécution, Nvidia fournit une feuille de calcul appelée \textit{occupancy calculator} dans laquelle ont peut entrer les paramètres d'exécution d'un \textit{kernel} parallèle : nmobre de registres utilisés par chaque thread, quantité de mémoire partagée, modèle de GPU, dimensions de la grille. Le tableur retourne alors l'indice de charge (l'occupancy) qui traduit le rapport, à chaque instant, entre le nombre de warps actifs et le nombre maximal de warps par processeur (SM = Streaming Multiprocessor). L'occupancy se traduit donc par un indice compris entre 0 et 100\% et la recherche de performance semble devoir être la recherche de l'occupancy maximale.

 
 

+Toutefois, comme l'a clairement demontré Volkov dans \cite{volkov2010better}, ce paradigme peut aisément être remis en cause et Volkov parvient effectivement à améliorer les peformances d'un certain nombre d'exemples génériques dans des conditions de faible valeur d'occupancy. 
+Enfin, nous avons pu constater deux grands modèles d'accès aux données : les algorithmes de filtrage usent quasiment tous de la mémoire partagée comme tampon d'accès aux données de l'image en mémoire globale (ou texture) alors que les algorithmes de segmentation performants s'en affranchissent. La raison en est clairement des motifs d'accès très irréguliers et non contigus pour ces derniers, rendant la gestion efficace de la mémoire partagée délicate et potentiellement si coûteuse qu'elle en devienne sans intérêt.
+Les chapitres suivants présentant nos contributions reviendront sur ces aspects en proposant des solutions pour accroître la performance des algorithmes parallélisés.