X-Git-Url: https://bilbo.iut-bm.univ-fcomte.fr/and/gitweb/these_gilles.git/blobdiff_plain/5d5e214d54eab1411ef86a32e11fccfaa53dbc45..1171799649e99aa6b7222c9f180de7523e5e7da4:/THESE/Chapters/chapter4/chapter4.tex?ds=inline diff --git a/THESE/Chapters/chapter4/chapter4.tex b/THESE/Chapters/chapter4/chapter4.tex index 55bc050..22657b8 100644 --- a/THESE/Chapters/chapter4/chapter4.tex +++ b/THESE/Chapters/chapter4/chapter4.tex @@ -1,22 +1,20 @@ +\section{Introduction} +Le concept de ligne de niveau dans les images a été introduit dès 1975 par Matheron \cite{matheron75}, puis Caselles \textit{et al.} \cite{caselles97} l'ont exploité et proposé le cadre définissant les \textit{images naturelles} comme les scènes photographiées, en intérieur ou en extérieur, à l'aide d'un appareil standard. Ces images vérifient alors l'hypothèse de gradient à valeurs bornées et peuvent être décomposées en un ensemble de lignes de niveaux. +Bertaux \textit{et al.} ont plus récemment proposé un algorithme de réduction du \textit{speckle} dans les images éclairées en lumière cohérente en introduisant, pour les pixels de l'image observée, une contrainte d'appartenance aux lignes de niveaux du modèle d'image non bruitée \cite{bertaux2004speckle}. L'image observée étant perturbée, les lignes de niveaux ne sont pas accessibles et il s'agit donc d'en estimer, localement par morceaux, la valeur et la forme, en se basant sur un modèle de forme pré établi. +Pour un pixel dont on cherche à estimer la valeur du niveau de gris, la contrainte d'appartenance à une ligne de niveau demeure locale, avec cependant un voisinage de forme et de taille (en nombre de pixels) variables en fonction des propriétés de l'image bruitée dans la zone concernée. +Ce voisinage, dont la forme, l'étendue et le niveau de gris sont déterminés par maximum de vraisemblance, appelé une \textit{isoline}, est une estimation locale de la ligne de niveau à laquelle appartient le pixel concerné. +Cette technique a montré qu'elle permettait de réduire très significativement le niveau de bruit tout en préservant les contours des objets. Elle s'est en revanche averée gourmande en ressources, ce qui a initialement conduit les auteurs à réduire la résolution de calcul des \textit{isolines} par application d'un maillage sur l'image bruitée. +Malgré cela, les temps de calcul demeuraient prohibitifs, avec une image de 2 millions de pixels traitée en 1 minute par un PIII-1GHz. +Comme nous l'avons déjà évoqué, l'amélioration des performances des microprocesseurs permet aujourd'hui de réduire assez considérablement ce temps de calcul. Cependant, la résolution des images à traiter à crû dans des proportions comparables, laissant les termes du compromis qualité/performance inchangés. \section{Présentation de l'algorithme} - -Bertaux \textit{et al.} ont proposé un algorithme de réduction du \textit{speckle} dans les images éclairées en lumière cohérente en introduisant un terme d'attache aux données lié aux lignes de niveaux du modèle d'image non bruitée \cite{bertaux2004speckle}. -La contrainte demeure locale vis à vis du pixel dont on cherche à estimer la valeur de niveau de gris, avec cependant un voisinage de forme et de taille (en nombre de pixels) variables en fonction des propriétés de l'image bruitée dans la zone concernée. -Ce voisinage, dont la forme, l'étendue et le niveau de gris sont déterminés par maximum de vraisemblance est appelé une \textit{isoline} et est une estimation locale de la ligne de niveau à laquelle appartient le pixel concerné. -Cette technique a montré qu'elle permettait de réduire très significativement le niveau de bruit tout en préservant les contours des objets. Elle s'est en revanche averée gourmande en ressources, ce qui a initialement conduit les auteurs a réduire la résolution de calcul des \textit{isolines} par application d'un maillage sur l'image corrompue. -Malgré cela, les temps de calcul demeurent prohibitfs, avec une image de 2 millions de pixels traitée en 1 minute par un PIII-1GHz. -Comme nous l'avons déjà évoqué, l'amélioration des performances des microprocesseurs permet aujourd'hui de réduire assez considérablement ce temps de calcul, mais la résoltion des images à traiter à elle aussi crû dans des proportions comparables laissant les termes du compromis qualité/performances inchangés. - \subsection{Formulation} Les \textit{isolines} sont des lignes brisées composées d'un ou plusieurs segments et construites par allongements successifs. Le niveau de gris affecté en sortie au pixel considéré est la valeur moyenne des niveaux de gris des pixels appartenant à l'\textit{isoline}. -Les segments sont de longueur $n$ fixe mais paramétrable et selectionnés parmi 32 motifs prédéterminés et mémorisés dans une table de référence notée $P_{n-1}$ dont un extrait est reproduit à la figure \ref{fig-lniv-p5q1} avec les motifs des segments correspondants. Tous les motifs sont composés du même nombre $n-1$ de pixels. -Pour chaque pixel de l'image d'entrée, le premier segment est choisi comme celui présentant la meilleure vraisemblance parmi les 32 possibles. Le choix d'intégrer ou non d'autres segments à l'\textit{isoline} et la sélection des segments à intégrer sont efféctués par évaluation d'un critère de vraisemblance généralisée dont l'obtention est détaillée dans la suite. +Les segments sont de longueur $n$ fixe mais paramétrable et leur \og forme \fg{} est sélectionnée parmi 32 motifs prédéterminés et mémorisés dans une table de référence notée $P_{n-1}$ dont un extrait est reproduit à la figure \ref{fig-lniv-p5q1} avec les motifs des segments correspondants. Tous les motifs sont composés du même nombre $a=n-1$ de pixels. +Pour chaque pixel de l'image d'entrée, le premier segment est choisi comme celui présentant la meilleure vraisemblance parmi les 32 possibles. Un test statistique GLRT (\textit{Generalized Likelihood Ratio Test}) détermine la pertinence d'intégrer ou non un segment supplémentaire à l'\textit{isoline} et le segment effectivement sélectionné pour l'allongement est le meilleur au sens du maximum de vraisemblance, parmi tous les segements satisfaisant au test GLRT. \begin{figure}[h] -\subfigure[Les 8 premières lignes de la table $P_5$ dont les éléments sont les positions relatives des pixels de chaque motif par rapport au pixel central.]{\includegraphics[height=4cm]{Chapters/chapter4/img/P5Q1.jpg}}\quad -\subfigure[Motifs des 8 premier segments candidats pour. Les pixels noirs representent le pixel traité (ou pixel central), qui n'appartient pas au motif. Les pixels gris sont ceux qui constituent le motif.]{ -$ +\subfigure[Les 8 premières lignes de la table $P_5$. Les éléments sont les positions relatives des pixels de chaque motif par rapport au pixel central.]{$ P_5 = \begin{bmatrix} (0,1)&(0,2)&(0,3)&(0,4)&(0,5)\\ @@ -30,53 +28,56 @@ P_5 = \ldots&\ldots&\ldots&\ldots&\ldots\\ \end{bmatrix} $ -} -\caption{\label{fig-lniv-p5q1}Détail des motifs et de leur représentation interne, pour la taille $d=5$. } +}\quad +\subfigure[Motifs des 8 premiers segments associés aux 8 premières lignes de $P_5$. Les pixels noirs représentent le pixel traité (ou pixel central), qui n'appartient pas au motif. Les pixels gris sont ceux qui constituent le motif.]{\includegraphics[height=4cm]{Chapters/chapter4/img/P5Q1.jpg}} +\caption{\label{fig-lniv-p5q1}Détail des motifs et de leur représentation interne, pour la taille $a=5$. } \end{figure} -\subsubsection{Isolines à un seul segment} +\subsubsection{Détermination du premier segment} -Pour chacun des pixels $(i,j)$ de l'image corrompue, on calcule la vraisemblance associée à chaque segment candidat de la table $P_{n-1}$ dans la région carrée $\omega$ centrée en $(i,j)$ et de coté $2n-1$. La région $\omega$ est l'union des deux sous-régions $S^n$ et $\overline{S^n}$ telles que $S^n$ décrit le segment candidat à évaluer comme un ensemble de $n$ pixels de coordonnées $(i_q,j_q)$ où $q\in [0..n[$. -La figure \ref{fig-lniv-regions} montre cette répartition pour $d=5$ et l'un des motifs ($p_{5,3}$). +Nous avons arbitrairement choisi un modèle simple où, au sein de la région carrée $\omega$ centrée en $(i,j)$ et de côté $2n-1$, les pixels du segment de ligne de niveau composent la région $S^n$ et possèdent la même valeur de niveu de gris. Les pixels de $\omega$ n'appartenant pas à $S^n$ composent la région $\overline{S^n}$ et aucun modèle ne leur est attaché. La région $\omega$ est ainsi l'union des des sous régions $S^n$ et $\overline{S^n}$ et le segment à évaluer est l'ensemble des $(a+1)$ pixels de coordonnées $(i_q,j_q)$ où $q\in [0..n[$. + +Pour chacun des pixels $(i,j)$ de l'image observée, on calcule la vraisemblance associée à chaque segment candidat de la table $P_{n-1}$ dans la région $\omega$. +La figure \ref{fig-lniv-regions} montre cette répartition pour $a=5$ et le motifs $p_{5,3}$. \begin{figure}[h] \center \includegraphics[height=3cm]{Chapters/chapter4/img/illustration_mv.jpg} -\caption{\label{fig-lniv-regions}. Exemple de la répartition des pixels dans la région $\omega$ pour le calcul de la vraisemblance, pour $n=6$.} +\caption{\label{fig-lniv-regions}Exemple de la répartition des pixels dans la région $\omega$ pour le calcul de la vraisemblance, pour $n=6$ ($a=5$).} \end{figure} -La densité de probabilité des valeurs des niveaux de gris des pixels de $S^d$ est supposée gaussienne de paramètres $\mu_{S^n}$ (moyenne) et $\sigma$ (écart type) inconnus. Les moyennes des niveaux de gris de chaque pixel sur $S^n$ sont inconnues et supposées indépendantes. +La densité de probabilité des valeurs des niveaux de gris des pixels de $S^n$ est supposée gaussienne de paramètres $\mu_{S^n}$ (moyenne) et $\sigma$ (écart type) inconnus. Les moyennes des niveaux de gris de chaque pixel sur $S^n$ sont inconnues et supposées indépendantes. Soit $Z$ l'ensemble des niveaux de gris des pixels de $\omega$ et $\{\mu_{ij}\}_{\overline{S^n}}$ l'ensemble des valeurs moyennes des pixels de $\overline{S^n}$. On peut écrire la probabilité \[ -P[Z|S^n, \mu_{S^n}, \{\mu_{ij}\}_{S^n}, \sigma] +P[Z|S^n, \mu_{S^n}, \{\mu_{ij}\}_{\{overline{S^n}}, \sigma] \] qui se développe comme suit, en distinguant les contributions de $S^n$ et $\overline{S^n}$ \begin{eqnarray} -\displaystyle \prod_{(i,j)\in S^n}{P\left[z(i,j) | \mu_{S^n}, \sigma \right]} \displaystyle\prod_{(i,j)\in \overline{S^n}}{P\left[z(i,j) | \left\{\mu_{ij}\right\}_{\overline{S^n}}, \sigma \right]} +\displaystyle \prod_{(i,j)\in S^n}{P\left[v(i,j) | \mu_{S^n}, \sigma \right]} \displaystyle\prod_{(i,j)\in \overline{S^n}}{P\left[v(i,j) | \left\{\mu_{ij}\right\}_{\overline{S^n}}, \sigma \right]} \label{LL2} \end{eqnarray} Nous cherchons alors à déterminer l'ensemble $S^n$ qui maximise la valeur de l'expression \eqref{LL2} ci dessus. -Or, sur $S^n$, les niveaux de gris $z(i,j)$ peuvent aussi être pris comme les estimations $\widehat{\mu_{ij}}$ des moyennes $\mu_{ij}$. -Le sesond terme de l'expression \eqref{LL2} devient donc +Or, d'après notre modèle, sur $\overline{S^n}$, les niveaux de gris $z(i,j)$ sont aussi les estimations $\widehat{\mu_{ij}}$ de leurs moyennes $\mu_{ij}$. +Le second terme de l'expression \eqref{LL2} devient donc \begin{eqnarray} -\displaystyle\prod_{(i,j)\in \overline{S^n}}{P\left[z(i,j) | \left\{\widehat{\mu_{ij}}\right\}_{\overline{S^n}}, \sigma \right]=1} +\displaystyle\prod_{(i,j)\in \overline{S^n}}{P\left[v(i,j) | \left\{\widehat{\mu_{ij}}\right\}_{\overline{S^n}}, \sigma \right]=1.} \end{eqnarray} Il reste alors le premier terme de \eqref{LL2}, soit l'expression de la vraisemblance généralisée donnée par : \begin{eqnarray} -\displaystyle \prod_{(i,j)\in S^n}{P\left[z(i,j) | \mu_{S^n}, \sigma \right]} +\displaystyle \prod_{(i,j)\in S^n}{P\left[v(i,j) | \mu_{S^n}, \sigma \right]} \label{GL} \end{eqnarray} -Que l'on peut développer, en remplaçant la probabilité par son expression. Il vient : +On peut développer cette expression en remplaçant la probabilité, gaussienne, par son expression. Il vient : \begin{eqnarray} \displaystyle \prod_{(i,j)\in S^n}\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{\left(z(i,j)-\mu_{S^n}\right)^2}{2\sigma^2}} \label{GL2} \end{eqnarray} -En prenant le logarithme et en developpant, on obtient alors l'expression suivante, dite \textit{log-vraisemblance} +En prenant le logarithme et en développant, on obtient alors l'expression suivante de la \textit{log-vraisemblance} \begin{eqnarray} \displaystyle -\frac{n}{2}log\left(2\pi\right) - \frac{n}{2}log\left(\sigma^2\right) - \frac{n}{2} \label{LL1} @@ -85,8 +86,8 @@ où le vecteur des paramètres $(\mu_{S^n}, \sigma)$ est lui même obtenu par es $$ \left( \begin{array}{l} -\widehat{\mu_{S^n}} = \displaystyle\frac{1}{n} \sum_{(i,j)\in S^n} z(i,j) \\ -\widehat{\sigma^2} = \displaystyle\frac{1}{n} \sum_{(i,j)\in S^n} \left(z(i,j) - \widehat{\mu_{S^n}}\right)^2 \\ +\widehat{\mu_{S^n}} = \displaystyle\frac{1}{n} \sum_{(i,j)\in S^n} v(i,j) \\ +\widehat{\sigma^2} = \displaystyle\frac{1}{n} \sum_{(i,j)\in S^n} \left(v(i,j) - \widehat{\mu_{S^n}}\right)^2 \\ \end{array} \right. $$ @@ -94,25 +95,25 @@ $$ Le motif retenu pour le segment est celui qui maximise l'expression de \eqref{LL1}. \subsubsection{Isolines composées de plusieurs segments - critère d'allongement} - -L'objectif poursuivi en cherchant à étendre la portée des isolines est d'améliorer la force du filtrage en intégrant plus de valeurs de niveaux de gris dans le calcul de la moyenne qui deviendra la valeur de sortie filtrée. -Pour cela nous permettons à chaque \textit{isoline}, comportant initialement un seul segment, d'être prolongée par d'autres segments, chaque allongement faisant l'objet d'une validation selon un critère de vraisemblance généralisée. -L'évaluation de l'ensemble des \textit{isolines} pouvant être construite sur ce modèle présente un coût prohibitif en temps de calcul et l'idée en a donc été abandonnée. -À la place, nous effectuons une sélection à chaque étape d'allongement : on évalue l'ensemble des 32 allongements possibles et si au moins un des motifs est accepté, on retient l'\textit{isoline} ayant la meilleure vraisemblance. Ce processus est répeté tant qu'au moins un motif représente un allongement valide. +Pour un bruit indépendant, la variance varie de manière inversement proportionnelle à la racine du nombre d'échantillons. +L'objectif est alors d'améliorer la force du filtrage en intégrant plus de valeurs de niveaux de gris dans le calcul de la moyenne qui deviendra la valeur de sortie filtrée. +Pour cela nous permettons à chaque isoline, comportant initialement un seul segment, d'être prolongée par d'autres segments, chaque allongement faisant l'objet d'une validation par un test (GLRT) issu de la théorie de la détection statistique (voir \cite{van2004detection}). +L'évaluation de l'ensemble des isolines pouvant être construite sur ce modèle présente un coût prohibitif en temps de calcul et l'idée en a donc été abandonnée. +À la place, nous effectuons une sélection à chaque étape d'allongement : on évalue l'ensemble des 32 allongements possibles et si au moins un des motifs est accepté, on retient l'\textit{isoline} ayant la meilleure vraisemblance. Ce processus est répété tant qu'au moins un motif représente un allongement valide. \begin{figure}[h] \center -\includegraphics[height=4cm]{Chapters/chapter4/img/exemple_extension_1.jpg} -\caption{\label{fig-lniv-allongement}Allongement du segment $S^n$. Deux candidats $S^{p'}$ et $S^{p''}$ sont évalués au travers du critère GLRT de l'équation \eqref{GLRT} que seul $S^{p''}$ s'avère satisfaire. a) Représentation dans le plan de l'image. b) Évolution des niveaux de gris en fonction de la position des pixels dans les lignes brisées ainsi formées.} +\includegraphics[height=5cm]{Chapters/chapter4/img/exemple_extension_1.jpg} +\caption{\label{fig-lniv-allongement}Allongement du segment $S^n$. Deux candidats $S^{p'}$ et $S^{p''}$ sont évalués au travers du test GLRT de l'équation \eqref{GLRT} que seul $S^{p''}$ s'avère satisfaire. a) Représentation dans le plan de l'image. b) Évolution des niveaux de gris en fonction de la position des pixels dans les lignes brisées ainsi formées.} \end{figure} -Soit $S^n$ une \textit{isoline} précédemment validée et $S^p$ un segment connecté à $S^n$ de telle sorte qu'il représente un potentiel allongement de $S^n$. Une situation de cette nature est représentée à la figure \ref{fig-lniv-allongement} avec un premier segment valide et deux candidats $S^{p'}$ et $S^{p''}$. -À gauche de la figure est reproduite une petite zone d'image réelle, suffisamment grossie pour permettre de bien individualiser les pixels et à laquelle ont été superposés les trois segments en question. Les deux relevés de la partie droite montrent quant à eux l'évolution des valeurs des niveaux de gris des pixels des deux \textit{isolines} possibles que sont $S^nS^{p'}$ et $S^nS^{p''}$. On a également identifié les différents sous-ensembles de pixels $S^n$, $S^{p'}$ et $S^{p''}$ ainsi que les valeurs moyennes de chacun. +Soit $S^n$ une isoline précédemment validée et $S^p$ un segment connecté à $S^n$ de telle sorte qu'il représente un allongement potentiel de $S^n$. Une situation de cette nature est représentée à la figure \ref{fig-lniv-allongement} avec un premier segment valide et deux candidats $S^{p'}$ et $S^{p''}$. +À gauche de la figure est reproduite une petite zone d'image réelle, suffisamment grossie pour permettre de bien individualiser les pixels et à laquelle ont été superposés les trois segments en question. Les deux relevés de la partie droite montrent quant à eux l'évolution des valeurs des niveaux de gris des pixels des deux isolines possibles que sont $S^nS^{p'}$ et $S^nS^{p''}$. On a également identifié les différents sous-ensembles de pixels $S^n$, $S^{p'}$ et $S^{p''}$ ainsi que les valeurs moyennes de chacun. À la lecture de ces deux représentations, on peut aisément imaginer que $S^{p'}$ ne soit pas retenu comme extension valide de $S^n$, au contraire de $S^{p''}$. -Pour formaliser ce que notre intuition semble nous dicter dans l'exemple précedent, nous comparons les log-vraisemblances des deux situations suivantes : +Pour formaliser, nous comparons les log-vraisemblances des deux situations suivantes : \begin{enumerate} -\item Le segment $S^p$ {\bf est} une extension valide pour $S^n$. Ils forment donc une \textit{isoline} $S^nS^p$.\\ +\item Le segment $S^p$ {\bf est} une extension valide pour $S^n$. Ils forment donc tous deux une isoline $S^nS^p$.\\ Dans ce cas et par hypothèse, la valeur moyenne des niveaux de gris est définie sur $S^nS^p$ et vaut $\mu_{S^nS^p}$. D'après \eqref{LL2}, la log-vraisemblance est alors donnée par \begin{eqnarray} \displaystyle -\frac{(n+p)}{2}\left(log\left(2\pi\right)+1\right) - \frac{(n+p)}{2}log\left(\widehat{\sigma_1}^2\right) @@ -128,30 +129,33 @@ où \hspace{3cm}$\widehat{\sigma_1^2} = \displaystyle\frac{1}{n+p} \sum_{(i,j)\i où \hspace{3cm}$\widehat{\sigma_2^2} = \displaystyle\frac{1}{n+p} \left( \sum_{(i,j)\in S^n} \left(v(i,j) - \widehat{\mu_{S^n}}\right)^2 + \sum_{(i,j)\in S^p} \left(v(i,j) - \widehat{\mu_{S^p}}\right)^2\right) $. \end{enumerate} -La différence entre \eqref{LLNP} et \eqref{LLNP2} nous donne l'expression du critère GLRT (Generalized Likelihood Ratio Test) +La différence entre \eqref{LLNP} et \eqref{LLNP2} nous donne l'expression du test GLRT \begin{eqnarray} T(S^n, S^p, T_{max}) = T_{max}- (n+p)\left[log\left(\widehat{\sigma_1}^2\right) - log\left(\widehat{\sigma_2}^2\right) \right] \label{GLRT} \end{eqnarray} -où $T_{max}$ est un seuil arbitrairement fixé.\\ -Un allongement de $S^n$ par $S^p$ est validé si $T(S^n, S^p, T_{max}) > 0$. +Dans cette étude, la valeur du seuil de détection $T_{max}$ est déterminée empiriquement, de sorte à produire des résultats visuels et chiffrés satisfaisants. Une détermination rigoureuse de la valeur adéquate de $T_{max}$ est envisageable, en fonction de la probabilité de fausse alarme visée, par le calcul de la PDF de $T$, mais dépasse le cadre de nos travaux. L'allongement de $S^n$ par $S^p$ est finalement validé si $T(S^n, S^p, T_{max}) > 0$. -\section{Implémentation parallèle} -Les \textit{isolines} sont construites segment après segment. On leur permet ainsi de suivre des formes courbes. La validité d'un segment et son éventuelle sélection sont soumises au critère décrit dans le paragraphe précédent. -Une limitation du nombre de segments candidats nous est également apparue pertinente, car elle permet d'apporter une réponse aux points suivants : -\begin{enumerate} -\item la sélection du premier segment est cruciale mais il n'est pas prouvé que l'\textit{isoline} la meilleure ait pour premier segment celui qui a été effectivement selectionné en premier. Une telle erreur sur la direction primaire peut s'avérer très pénalisante pour la qualité du traitement. C'est pourquoi nous conduisons en parallèle les allongements des 32 \textit{isolines}, chacune ayant l'un des motifs permis (voir figure \ref{fig-lniv-p5q1}) comme premier segment. -\item évaluer systématiquement les 32 motifs pour chaque extension peut alors rendre l'algorithme très coûteux. En effet, si $q$ est le nombre de segments maximum autorisés pour une \textit{isoline}, le nombre d'évaluations à effectuer se monte à $32^q$ par pixel. Cela représente un total de $\mathbf{3,5.10^{13}}$ évaluations pour des \textit{isolines} de 5 segments dans une image de 1024$^2$ pixels. -\item permettre à tout allongement de se faire dans chacune des 32 directions risque de générer des \textit{isolines} oscillant entre les deux extrémités de l'un de ses segments, ou bien s'enroulant sur elles-même au delà du simple rebouclage. . -\item une ligne de niveau ne peut pas se couper, donc une \textit{isoline} ne peut être composée de segments qui se croisent. -\end{enumerate} - -Les contraintes des points 3 et 4 ci-dessus nous ont conduit à limiter la déviation angulaire pouvant résulter de toute procédure d'allongement. Nous notons $\Delta d_{max}$ l'écart maximal toléré entre les indices des motifs de deux segments successifs. -Le choix de la valeur de $\Delta d_{max}$ adaptée depend de la taille des segments ainsi que du nombre maximal de segments que peut comporter une \textit{isoline}. -L'autre conséquence de cette limitation est la diminution du nombre total d'évaluations nécessaires. Si $\Delta d_{max} = 2$, le nombre d'évaluations effectuées dans l'exemple du point 2 passe ainsi de $\mathbf{3,5.10^{13}}$ à $1024^2\times 32\times 5^{q-1} = \mathbf{2,0.10^{10}}$ soit 1500 fois moins. +\section{Modélisation des isolines pour l'implémentation parallèle sur GPU} +Les isolines sont construites segment après segment afin de pouvoir approcher des formes courbes. La validité d'un segment et son éventuelle sélection sont soumises au tests statistiques décrits dans le paragraphe précédent. +La solution la plus évidente, que nous qualifierons de \textit{globale}, consiste à soumettre l'ensemble des isolines possibles aux tests statistiques pour sélectionner la meilleure. Toutefois, évaluer systématiquement les 32 motifs pour chaque extension peut rendre l'algorithme très coûteux. En effet, si $q$ est le nombre de segments maximum autorisés pour une isoline, le nombre d'évaluations à effectuer se monte à $32^q$ par pixel. Cela représente par exemple un total de $\mathbf{3,5.10^{13}}$ évaluations pour des isolines de $q=5$ segments dans une image de 1024$\times$1024 pixels. + +De plus, cette solution génère potentiellement des artefacts pour le cas où des isolines s'enrouleraient sur elles-mêmes au delà du simple rebouclage, ou bien feraient des aller-retours entre les deux extrémités d'un segment, ou encore présenteraient des segments qui se croisent. + +Pour toutes ces raisons, nous avons limité le nombre de segments candidats à chaque étape d'allongement, mais n'effectuons pas de sélection directe du premier segment. Nous conduisons en parallèle les allongements de 32 isolines, chacune ayant l'un des motifs de $P_a$ comme premier segment (voir figure \ref{fig-lniv-p5q1}). +Cela permet d'éviter une erreur sur la direction primaire potentiellement très pénalisante pour la qualité du traitement, tout en conservant une quantité d'évaluations réalisable dans un temps assez court. +% \begin{enumerate} +% \item la sélection du premier segment est cruciale mais il n'est pas prouvé que la meilleure isoline soit celle qui a pour premier segment celui qui a été effectivement sélectionné en premier. Une telle erreur sur la direction primaire peut s'avérer très pénalisante pour la qualité du traitement. C'est pourquoi +% \item +% \item permettre à tout allongement de se faire dans chacune des 32 directions risque de générer des isolines oscillant entre les deux extrémités de l'un de ses segments, ou bien s'enroulant sur elles-même au delà du simple rebouclage. +% \item une ligne de niveau ne peut pas se couper, donc une isoline ne peut pas être composée de segments qui se croisent. +% \end{enumerate} +La solution retenue consiste à limiter la déviation angulaire pouvant résulter de toute procédure d'allongement. Nous notons $\Delta d_{max}$ l'écart maximal toléré entre les indices des motifs de deux segments successifs. +Le choix d'une valeur de $\Delta d_{max}$ adaptée dépend de la taille des segments ainsi que du nombre maximal de segments que peut comporter une isoline. +L'autre conséquence de cette limitation est la diminution du nombre total d'évaluations nécessaires. Si $\Delta d_{max} = 2$, le nombre d'évaluations effectuées dans l'exemple du point 2 passe ainsi à $1024^2\times 32\times 5^{q-1} = \mathbf{2,0.10^{10}}$ soit 1500 fois moins (avec $q=5$). \subsection{Isolines évaluées semi-globalement} -La première implémentation proposée et notée PI-LD consiste donc à conduire l'allongement des 32 \textit{isolines} candidates à leur terme, puis de sélectionner la plus vraisemblable des \textit{isolines} parmi celles qui partagent la plus grande longueur. L'exemple de la figure \ref{fig-lniv-pild} illustre ce processus pour la sélection du deuxième segment d'une isoline avec $d=5$ et $\Delta d_{max}=2$. +La première implémentation proposée et notée PI-LD (\textit{Poly Isolines with Limited Deviation}), consiste donc à conduire l'allongement des 32 isolines candidates à leur terme, puis de sélectionner la plus vraisemblable parmi celles qui partagent la plus grande longueur. L'exemple de la figure \ref{fig-lniv-pild} illustre ce processus pour la sélection du deuxième segment d'une isoline avec $a=5$ et $\Delta d_{max}=2$. \begin{figure}[h] \centering @@ -161,11 +165,11 @@ La première implémentation proposée et notée PI-LD consiste donc à conduire \subfigure[Troisième segment évalué, associé au motif $p_{5,2}$.]{\label{pild:sub3} \includegraphics{Chapters/chapter4/img/PI-LD_detail_sub3.jpg}}\quad \subfigure[Quatrième segment évalué, associé au motif $p_{5,3}$.]{\label{pild:sub4} \includegraphics{Chapters/chapter4/img/PI-LD_detail_sub4.jpg}}\quad \subfigure[Cinquième segment évalué, associé au motif $p_{5,4}$.]{\label{pild:sub5} \includegraphics{Chapters/chapter4/img/PI-LD_detail_sub5.jpg}} -\caption{Processus de sélection lors de l'allongement d'une \textit{isoline} comportant initialement deux segment $s_1$ et $s_2$. Dans cet exemple $d=5$ et $\Delta d_{max}=2$. Chaque segment évalué est soumis au critère GLRT. Si au moins un des segments présente un test GLRT positif, alors l'allongement est réalisé avec le segment qui forme l'\textit{isoline} la plus vraisemblable.} +\caption{Processus de sélection lors de l'allongement d'une isoline comportant initialement deux segments $s_1$ et $s_2$. Dans cet exemple $a=5$ et $\Delta d_{max}=2$. Chaque segment évalué est soumis au test GLRT. Si au moins un des segments présente un test GLRT positif, alors l'allongement est réalisé avec le segment qui forme l'isoline la plus vraisemblable.} \label{fig-lniv-pild} \end{figure} -La rapidité de cette implémentation est très supérieure à celle des algorithmes \textit{état de l'art} comme BM3D (\cite{}), la qualité du débruitage étant tout de même moindre. Le tableau \ref{tab-lniv-pild} rassemble les performances comparées de nos implémentations et de celle du BM3D en y ajoutant comme référence de vitesse d'exécution un simple filtre moyenneur. Les mesures ont été réalisées sur l'ensemble des images de la base de test de S. Lansel (université de Berkeley), devenue entre temps indisponible au téléchargement, mais qui représente toujours une base de référence pour comparer des implémentations d'algorithmes de débruitage. Les images en sont reproduites à la figure \ref{fig-lniv-imgslansel}. +La rapidité de cette implémentation est très supérieure à celle des algorithmes \textit{état de l'art} comme BM3D (\cite{Dabov06imagedenoising}), la qualité du débruitage étant tout de même moindre. Le tableau \ref{tab-lniv-results} rassemble les performances comparées de nos implémentations et de celle du BM3D en y ajoutant comme référence de vitesse d'exécution un simple filtre moyenneur. Les mesures ont été réalisées sur l'ensemble des images de la base de test de S. Lansel (université de Berkeley), devenue entre temps indisponible au téléchargement, mais qui représente toujours une base de référence pour comparer des implémentations d'algorithmes de débruitage. Les images en sont reproduites à la figure \ref{fig-lniv-imgslansel}. \begin{figure}[ht] \centering @@ -187,16 +191,16 @@ La rapidité de cette implémentation est très supérieure à celle des algorit \end{figure} -L'adaptation au fonctionnement du GPU n'est pas non plus optimale du fait de la nécessité de réaliser deux différents types de validations à chaque étape d'allongement : un test GLRT et une minimisation de log-vraisemblance. Cela induit de nombreuses branches divergentes à l'exécution du kernel principal qui sont sérialisée par le GPU et peuvent causer une perte de performance considérable. +L'adaptation de ce modèle au fonctionnement du GPU n'est pas non plus optimale du fait de la nécessité de réaliser, à chaque étape d'allongement, deux différents types de validation : un test GLRT et une minimisation de log-vraisemblance. Cela induit de nombreuses branches d'exécution divergentes dans le kernel principal, qui sont sérialisées par le GPU et causent une perte de performance considérable. -Une analyse plus poussée des \textit{isolines} construites nous montre que la proportion est relativement faible d'\textit{isolines} optimales dont le premier segment n'aurait pas été celui sélectionné en appliquant notre algorithme. -L'exemple représentatif de la figure \ref{fig-lniv-histo-singe} montre l'histogramme des différences constatée pour l'image du singe. Les autres images de l'ensemble de test fournissent des histogrammes très semblables qui sont reproduits en petit format à la figure \ref{fig-lniv-histo-autres}. +Une analyse plus poussée des isolines construites nous montre qu'il y a une proportion relativement faible d'isolines déterminées par la méthode \textit{globale} dont le premier segment s'écarte notablement de celui sélectionné en l'absence d'allongement, c'est-à-dire par PI-LD avec $q=1$. +L'exemple représentatif de la figure \ref{fig-lniv-histo-singe} montre l'histogramme des écarts angulaires constatés entre les deux estimateurs, pour l'image du singe en version bruitée ($\sigma=25$). Les autres images de l'ensemble de test fournissent des histogrammes très semblables qui sont reproduits en petit format à la figure \ref{fig-lniv-histo-autres}. On y observe que pour environ 60\% des pixels de l'image, il y a correspondance des directions et que pour 80\% des pixels, l'écart angulaire reste inférieur à 22,5 degrés (soit 2 indices des motifs). \begin{figure}[h] \centering \includegraphics[height=4.5cm]{Chapters/chapter4/img/histodir.png} -\caption{Histogramme des écarts angulaires entre la direction primaire de l'\textit{isoline} optimale et celle de l'\textit{isoline} sélectionnée, pour l'image du singe (Mandrill). -Pour la très grande majorité des pixels, le mode de sélection de l'\textit{isoline} ne génére pas d'erreur sur la direction du premier segment.} +\caption{Histogramme des écarts angulaires entre la direction primaire de l'isoline optimale et celle du segment sélectionné par PI-LD avec $q=1$ (sans allongement), pour l'image du singe (Mandrill). +Pour la très grande majorité des pixels, l'écart est nul.} \label{fig-lniv-histo-singe} \end{figure} @@ -214,21 +218,21 @@ Pour la très grande majorité des pixels, le mode de sélection de l'\textit{is \subfigure[Peppers]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_peppers_25_l3.png}} \subfigure[Stream]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_stream_25_l3.png}} \subfigure[Zelda]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_zelda_25_l3.png}} -\caption{Histogrammes des écarts angulaires entre la direction primaire de l'\textit{isoline} optimale et celle de l'\textit{isoline} sélectionnée, pour les images de l'ensemble de test de S. Lansel. La répartition des erreurs est semblable dans toutes ces images naturelles.} +\caption{Histogrammes des écarts angulaires entre la direction primaire de l'isoline optimale et celle de l'isoline sélectionnée, pour les images de l'ensemble de test de S. Lansel. La répartition des erreurs est semblable dans toutes ces images, mais également dans toute image naturelle.} \label{fig-lniv-histo-autres} \end{figure} -On observe également que les pixels pour lesquels la sélection du premier segment n'est pas robuste sont situés dans les zones de l'image ne contenant pas de forts gradients de niveaux de gris. +On observe également que les pixels pour lesquels la sélection du premier segment n'est pas robuste sont situés dans les zones de l'image ne contenant pas de forts gradients de niveaux de gris, ce qui est cohérent avec la difficulté d'identifier une direction privilégiée dans ces régions. -\subsection{Isolines à segments pre-évalués - modèle PI-PD}\label{subsection-pipd-intro} -Les observations précédentes nous indiquent que, dans les zones où la sélection du premier segment est robuste, il n'est pas nécessaire de conduire l'étape de sélection consécutive à chaque allongement. En effet, dans ce cas, étendre une \textit{isoline} de $t$ segments se terminant au point final $(i_t, j_t)$ du segment $s_t$ revient à sélectionner le premier segment de l'\textit{isoline} débutant en $(i_t, j_t)$. +\subsection{Isolines à segments pré-évalués - modèle PI-PD\label{subsection-pipd-intro}} +Les observations précédentes nous indiquent que, dans les zones où la sélection du premier segment est robuste, il n'est pas nécessaire de conduire l'étape de sélection consécutive à chaque allongement. Sous cette hypothèse, étendre une isoline se terminant au point final $(i, j)$ revient à sélectionner le premier segment de l'isoline débutant en $(i, j)$. -Cette technique réduit considérablement la quantité d'évaluations à effectuer à $32.q$ à seulement {\bf 160} évaluations par pixel, soit un total de $\mathbf{1,7.1O^8}$ dans le cas de l'exemple cité en introduction du chapitre. +Cette technique réduit considérablement la quantité d'évaluations à effectuer, la faisant passer de $32^q$ à seulement {\bf 160} évaluations par pixel (pour $q=5$), soit un total de $\mathbf{1,7.10^8}$ pour une image de 124$\times$1024 avec $a=5$. -Ce nouveau modèle, nommé PI-PD, permet donc de séparer complétement la phase de sélection par maximum de vraisemblance qui détermine le meilleur segment initial en chaque pixel, des phases d'allongement successif soumis au seul test GLRT. -Pour implémenter efficacement cet algorithme sur GPU, il faut donc répartir les calculs en deux kernels principaux : +Ce nouveau modèle, nommé PI-PD (\textit{Poly Isolines with Precomputed Directions}), permet donc de séparer complètement la phase de sélection du segment initial par maximum de vraisemblance des phases d'allongements successifs soumis au seul test GLRT. +Pour implémenter efficacement cet algorithme sur GPU, il faut alors répartir les calculs en deux kernels principaux : \begin{enumerate} -\item \texttt{kernel\_precomp()}, qui réalise la sélection du premier segment en chaque pixel $(i,j)$. La direction $d_1(i,j)$ ainsi déterminée est mémorisée dans une matrice $I_{\Theta}$. Pour effectuer les calculs relatifs au GLRT, il faut aussi connaitre, pour chaque segment $s_1$, la valeur des sommes partielles +\item \texttt{kernel\_precomp()} réalise la sélection du premier segment en chaque pixel $(i,j)$. La direction $d_1(i,j)$ ainsi déterminée est mémorisée dans une matrice $I_{\Theta}$. Pour effectuer les calculs relatifs au GLRT, il faut aussi connaître, pour chaque segment $s_1$, la valeur des sommes partielles \begin{eqnarray} C_x\left(Z\left(S_1\right)\right)= \sum_{(i,j)\in s_1} v(i,j) \label{cx} @@ -239,21 +243,20 @@ Pour implémenter efficacement cet algorithme sur GPU, il faut donc répartir le \label{cx2} \end{eqnarray} -Elles sont calculées et mémorisées dans une seconde matrice notée $I_{\Sigma}$. Remarquons que les traitements réalisés par ce kernel correspondent exactement au modèle d'\textit{isoline} à un seul segment présenté au début. Les détails de son implémentation sont donnés à l'algorithme \ref{algo-lniv-precomp}, les initialisations étant données par l'Algorithme \ref{algo-lniv-init}. -\item \texttt{kernel\_PIPD()} qui évalue les allongements successifs, qui ne nécessitent plus de sélection par maximum de vraisemblance, mais uniquememt la ealidation par GLRT. Les données nécessaires à l'évaluation du critère GLRT sont regroupées, en plus de l'image d'entrée, dans les matrices $P_d$, $I_{\Theta}$ et $I_{\Sigma}$ et ne sont donc plus à calculer à ce stade et permettent d'envisager des performances en hausse par rapport à la solution PI-LD. L'algorithme \ref{algo-lniv-pipd} fournit les détails de l'implémentation. +Elles sont calculées et mémorisées dans une seconde matrice notée $I_{\Sigma}$. Remarquons que les traitements réalisés par ce kernel correspondent exactement au modèle d'isoline à un seul segment présenté au début. Les détails de son implémentation sont donnés dans l'algorithme \ref{algo-lniv-precomp}, les initialisations étant données dans l'algorithme \ref{algo-lniv-init}. +\item \texttt{kernel\_PIPD()} évalue les allongements successifs, qui ne nécessitent plus de sélection par maximum de vraisemblance, mais uniquement la validation par GLRT. Les données nécessaires à l'évaluation du test GLRT sont regroupées, outre dans l'image d'entrée, dans les matrices $P_d$, $I_{\Theta}$ et $I_{\Sigma}$ et ne sont donc plus à calculer à ce stade. Cela permet d'envisager des performances en hausse par rapport à la solution PI-LD. L'algorithme \ref{algo-lniv-pipd} fournit les détails de l'implémentation de ce kernel. \end{enumerate} -Les schémas de la figure \ref{fig-lniv-pipd} illustrent les étapes décrites ci-dessus de l'allongement d'une \textit{isoline} par la méthode PI-PD. \begin{figure}[h] +Les schémas de la figure \ref{fig-lniv-pipd} illustrent les étapes décrites ci-dessus de l'allongement d'une isoline par la méthode PI-PD. \begin{figure}[h] \centering \subfigure[Isoline avec 2 segments $s_1$ et $s_2$ déjà validés.]{\includegraphics[width=2.3cm]{Chapters/chapter4/img/PI-PD_detail_sub1.jpg}}\quad \subfigure[La direction de $s_3$ est l'élément $(i_2,j_2)$ de $I_{\Theta}$.]{\includegraphics[width=5cm]{Chapters/chapter4/img/PI-PD_detail_sub2.jpg}}\\ \subfigure[Le motif de $s_3$ est lu dans $p_5$ et appliqué en $(i_2,j_2)$. $C_x$ et $C_{x^2}$ sont données par $I_{\Sigma}(i_2,j_2)$ et le test GLRT est effectué.]{\includegraphics[width=4cm]{Chapters/chapter4/img/PI-PD_detail_sub3.jpg}}\quad \subfigure[Si l'allongement est validé, $s_3$ est définitivement intégré.]{\includegraphics[width=2.7cm]{Chapters/chapter4/img/PI-PD_detail_sub4.jpg}} -\caption{Exemple d'application du procédé d'allongement à une \textit{isoline} comprenant initialement 2 segments. la longueur des segments est $d=5$. Le procédé se répète jusqu'à ce que le test GLRT échoue.} +\caption{Exemple d'application du procédé d'allongement à une isoline comprenant initialement 2 segments. la longueur des segments est $a=5$. Le procédé se répète jusqu'à ce que le test GLRT échoue.} \label{fig-lniv-pipd} \end{figure} - \begin{algorithm}[htb] \caption{Initialisations du modèle PI-PD, en mémoire du GPU.} \label{algo-lniv-init} @@ -301,37 +304,37 @@ $T2_{max} \leftarrow$ seuil GLRT pour la détection de bords\; $(i_1, j_1) \leftarrow (i, j)$ \tcc*[r]{premier segment} $(C_x^1, C_{x2}^1) \leftarrow I_{\Sigma}(i_1,j_1)$ \tcc*[r]{lecture depuis $I_{\Sigma}$} $d_1 \leftarrow I_{\Theta}(i,j)$ \tcc*[r]{lecture depuis $I_{\Theta}$} - $l_1 \leftarrow l$ \tcc*[r]{longueur de l'isoline} - $\sigma_1 \leftarrow (C_{x2}^1/l_1 - C_x^1)/l_1$\; - $(i_2, j_2) \leftarrow end~of~first~segment$\; + $l \leftarrow n$ \tcc*[r]{longueur de l'isoline} + $\sigma_1 \leftarrow (C_{x2}^1/l - C_x^1)/l$\; + $(i_2, j_2) \leftarrow fin~du~premier~segment$\; $(C_{x}^2, C_{x2}^2) \leftarrow I_{\Sigma}(i_2,j_2) $ \tcc*[r]{2$^{nd}$ segment} $d_2 \leftarrow I_{\Theta}(i_2,j_2)$\; - $\sigma_2 \leftarrow (C_{x2}^2/l - C_x^2)/l$ \; + $\sigma_2 \leftarrow (C_{x2}^2/n - C_x^2)/n$ \; % - \While{$GLRT(\sigma_1, \sigma_2, l_1, l) < T_{max}$}{ - $l_1 \leftarrow l_1 + l$ \tcc*[r]{allongement} + \While{$GLRT(\sigma_1, \sigma_2, l, n) < T_{max}$}{ + $l \leftarrow l + n$ \tcc*[r]{allongement} $(C_x^1, C_{x2}^1) \leftarrow (C_x^1, C_{x2}^1)+(C_x^2, C_{x2}^2)$\; - $\sigma_1 \leftarrow (C_{x2}^1/l_1 - C_x^1)/l_1$ \tcc*[r]{mise à jour} + $\sigma_1 \leftarrow (C_{x2}^1/l - C_x^1)/l$ \tcc*[r]{mise à jour} $(i_1,j_1) \leftarrow (i_2, j_2)$ \tcc*[r]{décalage} $d_1 \leftarrow d_2$\; - $(i_2, j_2) \leftarrow end~of~next~segment$\; + $(i_2, j_2) \leftarrow fin~du~segment~suivant$\; \tcc*[f]{segment suivant} $(C_{x}^2, C_{x2}^2) \leftarrow I_{\Sigma}(i_2,j_2) $\; $d_2 \leftarrow I_{\Theta}(i_2,j_2)$\; - $\sigma_2 \leftarrow (C_{s2}^2/l - C_s^2)/l$ \; + $\sigma_2 \leftarrow (C_{s2}^2/n - C_s^2)/n$ \; } } - $\widehat{I}(i, j) \leftarrow C_x^1/l_1$ \tcc*[r]{niveau de gris en sortie} + $\widehat{I}(i, j) \leftarrow C_x^1/l$ \tcc*[r]{niveau de gris en sortie} \end{algorithm} -Le processus d'allongement du modèle PI-PD est également soumis aux restrictions sur les oscillations et retours en arrière des segments, déjà enoncées pour le modèle PI-LD. Par ailleurs, nous lui avons ajouté la possibilité de gérer des segments plus épais, composés de 2 ou 3 parallèles aux motifs décrit par la matrice $P_d$. Cela a pour effet d'intégrer plus de pixels dans les calculs statistiques et d'augmenter en conséquence les gains sur le PSNR, en particulier pour traiter des images de grandes dimensions qui ne contiendraient pas de \textit{trop petits} détails que l'épaisseur des \textit{isolines} risquerait de flouter. -Cette possibilité rend cette solution encore plus versatile que notre référence BM3D dont les temps de calcul s'avèrent prohibitifs sur des grandes images, avec par exemple plus de 5 minutes pour $4096^2$ pixels (Xeon quad core E31245\@3.3GHz, 8Go RAM). +Le processus d'allongement du modèle PI-PD est également soumis aux restrictions sur les retours en arrière des segments, déjà énoncées pour le modèle PI-LD. Par ailleurs, nous lui avons ajouté la possibilité de gérer des segments plus épais, composés de 2 ou 3 segments parallèles aux motifs décrits par la matrice $P_d$. Pour l'épaisseur 3, on utilise chaque segment motif ainsi que les deux segments parallèles obtenus par translation et situés de part et d'autre du segment motif. Pour l'épaisseur 2, on ne conserve qu'un seul de ces deux segments ajoutés. Cet épaississement a pour effet d'intégrer plus de pixels dans les calculs statistiques et d'augmenter en conséquence les gains sur le PSNR, en particulier pour traiter des images de grandes dimensions qui ne contiendraient pas de \textit{trop petits} détails que l'épaisseur des isolines risquerait de flouter. +Cette possibilité rend notre solution encore plus versatile que la référence BM3D dont les temps de calcul s'avèrent prohibitifs sur des images de grandes dimensions, avec par exemple plus de 5 minutes pour 4096$\times$4096 pixels (Xeon quad core E31245\@3.3GHz, 8Go RAM). -Il demeure toutefois que l'\textit{isoline} construite n'est pas nécessairement la plus vraisemblable pour tous les pixels de l'image, les optimisations étant faites sous l'hypothèse de robustesse énoncée au paragraphe \ref{subsection-pipd-intro}. +Toutefois, il demeure que l'isoline construite n'est pas nécessairement la plus vraisemblable pour tous les pixels de l'image, les optimisations étant faites sous l'hypothèse de robustesse énoncée au paragraphe \ref{subsection-pipd-intro}. \subsection{Modèle PI-PD hybride} -Le manque de robustesse de la sélection des segments dans certaines zones provient du petit nombre de pixels impliqués dans les calculs statistiques. Ces régions sont celles où la pente de la surface définie par les niveaux de gris des pixels, pris comme élévations, est faible vis à vis du bruit qui perturbe l'image. Par souci de concision, nous nommons ces régions des LSR (pour Low Slope Region). +Le manque de robustesse de la sélection des segments dans certaines zones provient du petit nombre de pixels impliqués dans les calculs statistiques. Ces régions sont celles où la pente de la surface définie par les niveaux de gris des pixels, pris comme élévations, est faible vis à vis du bruit qui perturbe l'image. Par souci de concision, nous nommons ces régions LSR (\textit{Low Slope Regions}, régions à faible pente). Pour illustrer ce comportement du modèle PI-PD, on peut sélectionner une région de petite taille (11$\times$11 pixels) au sein d'une image de test. La zone d'étude, repérée à la figure \ref{fig-lniv-lsr1}, est choisie pour ses propriétés particulières : deux \textit{plateaux} séparés par une transition nette, que l'on observe sur la représentation en trois dimensions de la figure \ref{fig-lniv-lsr-tirages-a}. \begin{figure}[h] @@ -342,7 +345,7 @@ Pour illustrer ce comportement du modèle PI-PD, on peut sélectionner une régi \label{fig-lniv-lsr1} \end{figure} -Les figures \ref{fig-lniv-lsr-tirages-b} et \ref{fig-lniv-lsr-tirages-c} montrent la même zone de l'image après qu'elle ait été corrompue par deux tirages d'un bruit gaussien de mêmes paramètres $\mu$ et $\sigma$. Les deux diagrammes \ref{fig-lniv-lsr-tirages-d} et \ref{fig-lniv-lsr-tirages-e} représentent quant à eux les directions primaires, modulo $\pi$, des \textit{isolines} débutant en chaque pixel de la fenêtre. On observe que la détermination de la direction est robuste dans la bande de transition entre les LSR, alors que pour ces dernières, on constate une grande variabilité. +Les figures \ref{fig-lniv-lsr-tirages-b} et \ref{fig-lniv-lsr-tirages-c} montrent la même zone de l'image après qu'elle ait été corrompue par deux tirages d'un bruit gaussien de mêmes paramètres $\mu$ et $\sigma$. Les deux diagrammes \ref{fig-lniv-lsr-tirages-d} et \ref{fig-lniv-lsr-tirages-e} représentent quant à eux les directions primaires des isolines, modulo $\pi$, débutant en chaque pixel de la fenêtre. On observe que la détermination de la direction est robuste dans la bande de transition entre les LSR, alors que pour les LSR elles-mêmes, on constate une grande variabilité. \begin{figure}[h] \centering @@ -355,48 +358,47 @@ Les figures \ref{fig-lniv-lsr-tirages-b} et \ref{fig-lniv-lsr-tirages-c} montren \label{fig-lniv-lsr-tirages} \end{figure} -Dans les zones LSR, l'application du modèle PI-PD ne revêt donc que peu de sens, mais la quête de performance nous interdit d'y appliquer par exemple le modèle PI-LD décrit précédemment. Le meilleur estimateur dans une zone LSR étant la valeur moyenne, nous proposons donc à la place : +Ainsi, dans les LSR, l'application du modèle PI-PD n'a que peu de sens, mais la quête de performance nous interdit d'y appliquer, par exemple, le modèle PI-LD décrit précédemment. Considérant que le bruit est additif gaussien et que nous faisons une approximation \textit{plane} de la surface définie par les niveaux de gris des pixels à l'intérieur de la fenêtre $\omega$, le meilleur estimateur dans une zone LSR est la valeur moyenne. Nous proposons donc : \begin{enumerate} -\item d'identifier les zones à faible pente en concevant un kernel détecteur (\texttt{kernel\_LSR\_detector()}). -\item d'appliquer un simple filtre moyenneur dans les zones désignée LSR par le détecteur et le PI-PD partout ailleurs. -\item de n'appliquer le moyenneur que sur les pixels appartenant à la partie de zone LSR, lorsque la fenêtre du moyenneur se trouve à cheval sur deux zones de types différents. +\item d'identifier les zones à faible pente en concevant un kernel détecteur de bords (\texttt{kernel\_edge\_detector()}). +\item d'appliquer un simple filtre moyenneur dans les zones désignées LSR par le détecteur et le PI-PD partout ailleurs. +\item de n'appliquer le moyenneur que sur les pixels appartenant à la zone LSR lorsque la fenêtre du détecteur se trouve à cheval sur deux zones de types différents. \end{enumerate} -\subsubsection{Le détecteur de zone à faible pente} +\subsubsection{Le détecteur de bords} -Le principe retenu pour réaliser le détecteur de LSR est proche de celui mis en oeuvre pour valider les allongements des isolines : il s'agit de séparer la fenêtre d'observation autour du pixel considéré, en deux régions et d'effectuer un test GLRT pour déterminer s'il est vraisemblable ou non que ces deux régions forment un seul et même plan. Pour garantir la prise en compte d'éventuelles transitions dans toutes les directions, il faut effectuer le test avec des séparations de fenêtre dont les directions couvrent toute la plage angulaire, de $0$ à $\pi$. +Le principe retenu pour réaliser le détecteur de bords est proche de celui mis en oeuvre pour valider les allongements des isolines : il s'agit de séparer la fenêtre d'observation autour du pixel considéré en deux régions, puis d'effectuer un test GLRT pour déterminer s'il est vraisemblable ou non que ces deux régions forment un seul et même plan ou bien représentent deux zones homogènes séparées par un seuil (ou un coin). Pour garantir la prise en compte d'éventuelles transitions dans toutes les directions, il faut effectuer le test avec des séparations de fenêtre dont les directions couvrent toute la plage angulaire, de $0$ à $\pi$. -L'utilisation d'un test GLRT semblable à celui de l'équation \eqref{GLRT} sous-entend que les ensembles considérés n'ont aucun pixel en commun. Par ailleurs, afin d'optimiser les performances globales du détecteur, nous avons utilisé les motifs de la matrice $P_d$, n'ayant pas d'intersection entre eux et de direction $\Theta_{4i} = 4i\frac{\pi}{4}$. -La ligne de séparation entre les deux régions de la fenêtre est donc composée par les motifs de directions $\Theta_{4i}$ et $\Theta_{4(i+4)}$. Ces deux régions sont respectivement nommées arbitrairement $T$ et $B$, $T$ étant représentée comme la région \textit{haute} et $B$ comme la région \textit{basse} sur le schéma explicatif de la figure \ref{fig-lniv-detecteur} où $\Theta_{4i}=\frac{\pi}{4}$ et où les pixels affectés d'une élévation nulle sont les pixels de la fenêtre non impliqués dans le calcul du critère GLRT. En outre, les pixels de la limite sont supposés appartenir à la région $T$, ce qui implique qu'elle comprend au total les pixels correspondant à cinq motifs plus le pixel central, tandis que $B$ n'en comprend que l'équivalent de 3 motifs. +L'utilisation d'un test GLRT semblable à celui de l'équation \eqref{GLRT} sous-entend que les ensembles considérés n'ont aucun pixel en commun. Afin d'éviter de devoir déterminer de nouveaux ensembles de pixels pertinents, nous avons utilisé les motifs de la matrice $P_d$, n'ayant pas d'intersection entre eux et de directions $\Theta_{4i} = 4i\frac{\pi}{4}$. Ces motifs remplissent les critères pour établir l'expression d'un test GLRT. +La ligne de séparation entre les deux régions de la fenêtre est donc composée par les motifs de directions $\Theta_{4i}$ et $\Theta_{4(i+4)}$. Ces deux régions sont respectivement nommées arbitrairement $T$ et $B$, $T$ étant représentée comme la région \textit{haute} et $B$ comme la région \textit{basse} sur le schéma explicatif de la figure \ref{fig-lniv-detecteur} où $\Theta_{4i}=\frac{\pi}{4}$ et où les pixels affectés d'une élévation nulle sont les pixels non impliqués dans le calcul du test GLRT. En outre, les pixels de la limite sont supposés appartenir à la région $T$, ce qui implique qu'elle comprend au total les pixels correspondant à cinq motifs plus le pixel central, tandis que $B$ n'en comprend que l'équivalent de 3 motifs. \begin{figure}[ht] \centering - \includegraphics[width =5cm]{Chapters/chapter4/img/pattern_detecteur.jpg} - \caption{Motif de détection des zones à faible pente, pour le cas $\Theta=\Theta_4=45^{\circ}$. L'élévation des pixels permet juste de les distinguer selon 3 classes : l'élévation 1 est associée aux pixels de la région $H$, l'élévation 0.5 est associée à ceux de la région $L$ et l'élévation 0 désigne les pixels níntervnant pas dans la détection.} + \includegraphics[width =7cm]{Chapters/chapter4/img/pattern_detecteur-f.png} + \caption{Motif de détection des zones à faible pente, pour le cas $\Theta=\Theta_4=45^{\circ}$. L'élévation des pixels permet juste de les distinguer selon 3 classes : l'élévation 1 est associée aux pixels de la région $T$, l'élévation 0.5 est associée à ceux de la région $B$ et l'élévation 0 désigne les pixels n'intervenant pas dans la détection.} \label{fig-lniv-detecteur} \end{figure} -Les équations \eqref{LLNP}, \eqref{LLNP2} et \eqref{GLRT} nous permettent d'obtenir l'expression suivante pour le critère GLRT $T2$ +Les équations \eqref{LLNP}, \eqref{LLNP2} et \eqref{GLRT} nous permettent d'obtenir l'expression suivante pour le test GLRT $T2$ \begin{eqnarray} -T2 = T2_{max}- (8l+1)\left[log\left(\widehat{\sigma_3}^2\right) - log\left(\widehat{\sigma_4}^2\right) \right] +T2 = T2_{max}- (8a+1)\left[log\left(\widehat{\sigma_3}^2\right) - log\left(\widehat{\sigma_4}^2\right) \right] \label{GLRT2} \end{eqnarray} -avec $\widehat{\sigma_3}$ l'estimation de l'écart type de la situation où les deux demi-régions en formeraient une seule et $\widehat{\sigma_4}$ celle où une transition serait détectée entre les deux. Leurs expressions sont donc : +où $\widehat{\sigma_3}$ est l'estimation de l'écart type dans le cas où les deux demi régions en formeraient une seule et $\widehat{\sigma_4}$, l'estimation de l'écart type pour le cas où une transition serait détectée entre les deux. Leurs expressions sont donc : $$ \begin{array}{l} -\widehat{\sigma_3}^2 = \displaystyle\frac{1}{8l+1}\sum_{(i,j)\in T\cup B}\left( v(i,j) - \widehat{\mu_{T\cup B}} \right)^2 \\ +\widehat{\sigma_3}^2 = \displaystyle\frac{1}{8a+1}\sum_{(i,j)\in T\cup B}\left( v(i,j) - \widehat{\mu_{T\cup B}} \right)^2 \\ et\\ -\widehat{\sigma_4}^2 = \displaystyle\frac{1}{5l+1}\sum_{(i,j)\in T}\left(v(i,j) - \widehat{\mu_{T}} \right)^2 + \frac{1}{3l}\sum_{(i,j)\in B}\left(v(i,j)- \widehat{\mu_{B}} \right)^2 +\widehat{\sigma_4}^2 = \displaystyle\frac{1}{8a+1}\left(\sum_{(i,j)\in T}\left(v(i,j) - \widehat{\mu_{T}} \right)^2 + \sum_{(i,j)\in B}\left(v(i,j)- \widehat{\mu_{B}} \right)^2\right) \end{array} $$ -Le seuil de décision est noté $T2_{max}$ et d'après l'expression du critère \eqref{GLRT2}, une valeur négative du critère signifie la détection d'une transition. Ainsi, lorsque les valeurs du critère $T2$ sont connues pour toutes les 8 directions $\Theta_{4i} (i\in [0..7])$, la valeur du niveau de gris de sortie pour le pixel central est déterminée selon la stratégie suivante : +Le seuil de détection est noté $T2_{max}$ et d'après l'expression du test \eqref{GLRT2}, une valeur négative signifie la détection d'un bord. Ainsi, lorsque les valeurs de $T2$ sont connues pour toutes les 8 directions $\Theta_{4i} (i\in [0..7])$, la valeur du niveau de gris de sortie pour le pixel central est déterminée selon la stratégie suivante : \begin{itemize} -\item si plus d'une valeur du critère est négative, alors on applique la valeur issue du modèle PI-PD. -\item si une seule valeur du critère est négative, le pixel central est vraisemblablement situé sur une transition nette et on applique la valeur moyenne de la région $T$ à laquelle il appartient. Cela permet de garantir des transitions visuellement plus douces entre les zones où le PI-PD s'est appliqué et les zones moyennées. -\item si aucune valeur du critère n'est négative, alors la région autour du pixel central est vraisemblablement une LSR. En conséquence, on applique -la valeur moyenne de la zone. +\item si plus d'un bord a été détecté, alors on applique la valeur issue du modèle PI-PD. +\item si un seul bord a été détecté, le pixel central est vraisemblablement situé sur une transition nette et on applique la valeur moyenne des motifs de la région $T$ à laquelle il appartient. Cela permet de garantir des transitions visuellement plus douces entre les zones où le PI-PD est appliqué et les zones moyennées. +\item si aucun bord n'a été détecté, alors la région autour du pixel central est vraisemblablement une LSR. En conséquence, on applique la valeur moyenne de la zone. \end{itemize} La figure \ref{fig-lniv-classification} présente le résultat de la classification des pixels d'une image bruitée, pour $T2_{max}=2$. On y remarque en particulier que les pixels noirs, pour lesquels s'appliquera le PI-PD, sont situés sur des transitions bien définies. @@ -405,47 +407,48 @@ La figure \ref{fig-lniv-classification} présente le résultat de la classificat \centering \subfigure[Image bruitée]{\includegraphics{Chapters/chapter4/img/airplane_noisy_small.jpg}}\qquad \subfigure[Classification des pixels. ]{\includegraphics{Chapters/chapter4/img/img_bords_T2_small.jpg}} -\caption{Classification des pixels d'une image bruitée, pour une valeur de seuil $T2=2$ du détecteur. b) Les pixels en noir sont ceux à qui le PI-PD sera appliqué. Les pixels en blancs se verront appliquer une moyenne sur tout ou partie du voisinage.} +\caption{Classification des pixels d'une image bruitée, pour une valeur de seuil $T2=2$ du détecteur. (b) Les pixels en noir sont ceux à qui le PI-PD sera appliqué. Les pixels gris se verront appliquer une moyenne sur tout ou partie du voisinage.} \label{fig-lniv-classification} \end{figure} +Les détails d'implémentation du détecteur sont donnés par l'algorithme \ref{algo-lniv-detecteur}. Pour en optimiser les performances, les sommes individuelles $sum_{\Theta}$ sont pré-calculées aux lignes 7 à 10 pour les 8 motifs concernés. L'évaluation des 8 configurations angulaires est effectuée ensuite de la ligne 11 à la ligne 25. -Les détails d'implémentation du détecteur sont donnés par l'algorithme \ref{algo-lniv-detecteur}. Pour en optimiser les performances, les sommes partielles $sPat$ et $sPat2$ évaluées successivement pour les régions $T$ et $B$ aux lignes 11 et 12, puis 16 et 17 sont pré-calculées et mémorisées dans 8 registres. +Le choix que nous avons fait pour ce détecteur de bord sont stratégiques et ont été guidés par la recherche de performance. Il est cependant clair qu'un filtrage plus fort aurait été obtenu dans les LSR par un moyenneur intégrant l'ensemble des pixels de la fenêtre. Cette piste reste à approfondir à la lumière des résultats obtenu dans l'implémentation de l'opération de convolution. \begin{algorithm}[ht] \caption{Détecteur de zones à faible pente (LSR) \texttt{kernel\_LSR\_detector()}} \label{algo-lniv-detecteur} -\ForEach(\tcc*[f]{\textbf{in parallel}}){pixel $(i,j)$}{ +\ForEach(\tcc*[f]{\textbf{en parallèle}}){pixel $(i,j)$}{ $\Theta \leftarrow 0$\tcc*[r]{Indice de la direction} $edgeCount \leftarrow 0$\; $sumEdge \leftarrow 0$\; - $nH \leftarrow 5l+1$\; - $nL \leftarrow 3l$\; + $nT \leftarrow 5l+1$\; + $nB \leftarrow 3l$\; + \While{($\Theta < 32$) }{ + $sum_{\Theta} \leftarrow \left(\displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I_{n tex}(i+y,j+x), \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I^2_{n tex}(i+y,j+x)\right)$ \; + $\Theta \leftarrow \Theta + 4$\; + } \While{($\Theta < 32$) }{ - $sumH \leftarrow (I_{ntex}(i,j), I_{ntex}^2(i,j))$\; - $sumL \leftarrow (0, 0)$\; + $sumT \leftarrow (I_{ntex}(i,j), I_{ntex}^2(i,j))$\; + $sumB \leftarrow (0,0)$ \; \For{($\alpha=\Theta$ to $\alpha=\Theta+16$ par pas de $4$)}{ - $sPat \leftarrow \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I_{n tex}(i+y,j+x)$\; - $sPat2 \leftarrow \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I^2_{n tex}(i+y,j+x)$\; - $sumH \leftarrow sumH + (sPat, sPat2)$\; + $sumT \leftarrow sumT + sum_{\Theta}$\; } \For{($\alpha=\Theta+20$ to $\alpha=\Theta+28$ par pas de $4$)}{ - $sPat \leftarrow \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I_{n tex}(i+y,j+x)$\; - $sPat2 \leftarrow \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I^2_{n tex}(i+y,j+x)$\; - $sumL \leftarrow sumL + (sPat, sPat2)$\; + $sumB \leftarrow sumB + sum_{\Theta}$ \; } - \If{($GLRT(sumH, nH, sumL, nL) > T2_{max}$)}{ + \If{($GLRT(sumT, nT, sumB, nB) > T2_{max}$)}{ $edgeCount \leftarrow edgeCount + 1$\; - $sumEdge \leftarrow sumH.x$\; + $sumEdge \leftarrow sumT.x$\; } $\Theta \leftarrow \Theta + 4$\; - } + } \tcc{niveau de gris de l'isoline} \If{($edgeCount == 0$)}{ - $\widehat{I}(i,j) \leftarrow \dfrac{(sumH.x + sumL.x)}{nH+nL}$ \tcc*[r]{LSR} + $\widehat{I}(i,j) \leftarrow \dfrac{(sumT.x + sumB.x)}{nT+nB}$ \tcc*[r]{LSR} } \If{($edgeCount == 1$)}{ - $\widehat{I}(i,j) \leftarrow \dfrac{(sumEdge)}{nH}$ + $\widehat{I}(i,j) \leftarrow \dfrac{(sumEdge)}{nT}$ } \If{($edgeCount > 1$)}{ $\widehat{I}(i,j) \leftarrow \widehat{I_{PIPD}}(i,j)$\tcc*[r]{PI-PD} @@ -456,24 +459,23 @@ Les détails d'implémentation du détecteur sont donnés par l'algorithme \ref{ \section{Résultats} -L'implémentation du PI-PD hybride à été appliquée aux 13 images de la base de test, dans leurs versions les plus bruitées, perturbées par un bruit gaussien de moyenne nulle et d'écart type 25. -Pour ce type d'images (taille, détails), les paramètres qui ses sont avérés optimaux sont $l=5$ pour la longueur des segments avec un maximum de 5 segments. -En ce qui concerne les seuils GLRT, nous avons testé l'ensemble des combinaisons de valeurs $T_{max}$ et $T2_{max}$ variant de 1 à 10 par pas de 0.5. -La combinaison $T_{max}=1$ et $T2_{max}=2$ s'est révelée la plus appropriée, en ce sens qu'elle représente l'optimum pour 11 des 13 images, sauf \textit{peppers} et \textit{zelda}, pour lesquelles une combinaison $T_{max}=2$ et $T2_{max}=2$ permet d'améliorer l'indice de similarité MSSIM respectivement de 0.03 et 0.02. +L'implémentation du PI-PD hybride a été appliquée aux 13 images de la base de test, dans leurs versions les plus bruitées, perturbées par un bruit gaussien de moyenne nulle et d'écart type 25. +Pour ce type d'images (taille, détails), les paramètres qui se sont avérés optimaux sont $a=5$ pour la longueur des segments avec un maximum de $q=5$ segments. +En ce qui concerne les seuils GLRT, nous avons testé l'ensemble des combinaisons de valeurs $T_{max}$ et $T2_{max}$ variant de 1 à 10 par pas de 0,5. +La combinaison $T_{max}=1$ et $T2_{max}=2$ s'est révélée la plus appropriée, car elle représente le meilleur compromis PSNR/MSSIM pour 11 des 13 images, sauf \textit{peppers} et \textit{zelda}, pour lesquelles une combinaison $T_{max}=2$ et $T2_{max}=2$ permet d'améliorer l'indice de similarité MSSIM respectivement de 0,03 et 0,02. -Les images ainsi filtrées ont donc été caractérisées en termes de PSNR et de MSSIM et les résultats, regroupés dans la table \ref{tab-lniv-results}, sont comparés à ceux de la référence BM3D, ainsi qu'à ceux d'un simple filtre moyenneur GPU 5$\times$5, choisi comme référence en terme de rapidité et dont la taille de fenêtre permet des gain théoriques en PSNR du même ordre de grandeur que le PI-PD. +Les images filtrées ont été caractérisées en termes de PSNR et de MSSIM et les résultats, regroupés dans la table \ref{tab-lniv-results}, sont comparés à ceux de la référence BM3D, ainsi qu'à ceux d'un simple filtre moyenneur GPU 5$\times$5, choisi comme référence en terme de rapidité et dont la taille de fenêtre permet des gains théoriques en PSNR du même ordre de grandeur que le PI-PD. -Les mesures de qualité montrent que le PI-PD hybride améliore le PSNR de 1.5~dB et le MSSIM de 7,3\% en moyenne par rapport au moyenneur, au prix d'un t -emps de calcul multiplié par 128. Le BM3D fait encore progresser la qualité de 2.4~dB et 4,6\% en moyenne par rapport au PI-PD hybride, mais en mettant 475 fois plus de temps. -Le principal défaut du filtre proposé est la génération d'artefacts de type marches d'escalier (staircase effect), inhérente à tous les filtres de voisinage. Cependant, nous avons implémenté sur GPU la solution proposée par Buades dans \cite{BuadesCM06} et ainsi attenué nettement cet effet indésirable pour un coût de 0.2~ms. La valeur du PSNR de chaque image débruitée a ainsi été encore améliorée de 1~dB. +Les mesures de qualité montrent que le PI-PD hybride améliore en moyenne le PSNR de 1,5~dB et le MSSIM de 7,3\% par rapport au moyenneur, au prix d'un temps de calcul multiplié par 100, soit environ 7,3~ms, là où l'algorithme PI-LD prenait 35~ms. Le BM3D fait encore progresser la qualité de 2,4~dB et 4,6\% en moyenne par rapport au PI-PD hybride, mais en mettant 590 fois plus de temps que ce dernier, soit environ 4,3~s. +Le principal défaut du filtre proposé est la génération d'artefacts de type marches d'escalier (staircase effect), inhérente à tous les filtres de voisinage. Cependant, nous avons implémenté sur GPU la solution proposée par Buades dans \cite{BuadesCM06} et ainsi attenué nettement cet effet indésirable pour un coût de 0,2~ms. La valeur du PSNR de chaque image débruitée a ainsi été encore améliorée de 1~dB. La figure \ref{fig-lniv-exempleresultat} permet de constater le rendu visuel des traitements comparés, sur l'image entière ainsi que sur une zone grossie de l'image \textit{airplane}. \begin{figure}[ht] \centering \subfigure[Image \textit{airplane} bruitée.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_25_noisy.png}}\quad \subfigure[Image \textit{airplane} filtrée par moyenneur 5$\times$5.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_mean5.png}}\quad - \subfigure[Image \textit{airplane} filtrée par PI-PD hybride avec $l=5$, $n=25$, $T_{max}=2$ et $T2_{max}=2$.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_25_noisy_6_r50_T10_P2.png}}\quad - \subfigure[Image \textit{airplane} filtrée par PI-PD hybride.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_bm3d.png}}\\ + \subfigure[Image \textit{airplane} filtrée par PI-PD hybride avec $a=5$, $q=5$, $T_{max}=2$ et $T2_{max}=2$.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_25_noisy_6_r50_T10_P2.png}}\quad + \subfigure[Image \textit{airplane} filtrée par BM3D.]{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_bm3d.png}}\\ \subfigure{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_25_noisy_zoom.jpg}}\quad \subfigure{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_25_mean5_zoom.jpg}}\quad \subfigure{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_zoom_hybrid_6_r50_T10_P2.jpg}}\quad @@ -491,82 +493,83 @@ La figure \ref{fig-lniv-exempleresultat} permet de constater le rendu visuel des &(ms)&(ms) \\ \midrule Moyenneur & 0.07& 0.15 \\ -PI-PD hybride & 9.00& 0.15 \\ -BM3D & 4300& 0.00 \\ +PI-PD hybride & 7.30& 0.15 \\ +BM3D & 4300& \ldots \\ \bottomrule \end{tabular} \caption{Temps de calcul et de transfert des implémentations comparées. } \label{tab-lniv-chronos} \end{table} -Les temps de calcul des différentes implémentations testées ne dépendent que très peu du contenu de l'image, voire pas du tout pour le moyenneur. Ils sont présentés à la table \ref{tab-lniv-chronos}. Pour les implémentations GPU, il faut ajouter, dans le cas de traitements uniques (pas des séquences d'images), les temps de transfert des images vers la mémoire texture du GPU puis vers une zone de mémoire non paginée de l'hôte CPU, qui représentent un total de 0.15~ms pour les images de test. Notons que l'emploi de mémoire non-paginée pour la mémorisation des données côté CPU permet d'économiser 0.09~ms par image 8~bits. +Les temps de calcul des différentes implémentations testées dépendent très peu du contenu de l'image, voire pas du tout pour le moyenneur. Ils sont présentés à la table \ref{tab-lniv-chronos}. Pour les implémentations GPU, il faut ajouter, dans le cas de traitements uniques (hors séquences d'images), les temps de transfert des images vers la mémoire texture du GPU puis vers une zone de mémoire non paginée de l'hôte CPU, qui représentent un total de 0,15~ms pour les images de test, soit moins de 2\% du temps total du PI-PD hybride. Notons que l'emploi de mémoire pré-allouée (ne générant pas de défaut de page) pour la mémorisation des données côté CPU permet d'économiser 0,09~ms par image 8~bits, soit environ 1\% du temps total du PI-PD. Notons enfin que le traitement de séquences en haute définition (1920$\times$1080 pixels) au taux 20 images par seconde est rendu possible. \begin{table}[H] \scriptsize \centering -\begin{tabular}{crrrr} +\begin{tabular}{crrrrr} \toprule -\bf Image&\bf Noisy &\bf Moyenneur &\bf PI-PD&\bf BM3D \\ - & &\bf $5\times 5$ &\bf hybride & \\ - & \tiny{PSNR (dB)}& \tiny{gain (dB)}/noisy & \tiny{gain (dB)}/noisy& \tiny{gain (dB)}/noisy\\ - & \tiny{MSSIM} & \tiny{MSSIM} & \tiny{MSSIM}& \tiny{MSSIM}\\ +\bf Image&\bf Bruitée &\bf Moyenneur &\bf PI-LD &\bf PI-PD&\bf BM3D \\ + & &\bf $5\times 5$ & &\bf hybride & \\ + & \tiny{PSNR (dB)}& \tiny{gain (dB)}/noisy & \tiny{gain (dB)}/noisy & \tiny{gain (dB)}/noisy& \tiny{gain (dB)}/noisy\\ + & \tiny{MSSIM} & \tiny{MSSIM} & \tiny{MSSIM}&\tiny{MSSIM}& \tiny{MSSIM}\\ \midrule -airplane & 19.49 & 6.90 & 8.97 & 11.39 \\ - & 0.58 & 0.84 & 0.88 & 0.93 \\ +airplane & 19.49 & 6.90 & 8.94 &8.97 & 11.39 \\ + & 0.58 & 0.84 & 0.78 &0.88 & 0.93 \\ \midrule -barbara & 20.04 & 2.72 & 4.22 & 10.56 \\ - & 0.70 & 0.76 & 0.83 & 0.94 \\ +barbara & 20.04 & 2.72 & 4.84 &4.22 & 10.56 \\ + & 0.70 & 0.76 & 0.79 & 0.83 & 0.94 \\ \midrule -boat & 20.33 & 5.25 & 7.21 & 9.69 \\ - & 0.66 & 0.81 & 0.87 & 0.91 \\ +boat & 20.33 & 5.25 & 6.86 &7.21 & 9.69 \\ + & 0.66 & 0.81 & 0.81 &0.87 & 0.91 \\ \midrule -couple & 20.28 & 4.97 & 7.05 & 9.49 \\ - & 0.69 & 0.79 & 0.87 & 0.91 \\ +couple & 20.28 & 4.97 & 6.77 &7.05 & 9.49 \\ + & 0.69 & 0.79 & 0.82 &0.87 & 0.91 \\ \midrule -elaine & 19.85 & 8.86 & 9.09 & 10.75 \\ - & 0.59 & 0.86 & 0.87 & 0.91 \\ +elaine & 19.85 & 8.86 & 8.16 &9.09 & 10.75 \\ + & 0.59 & 0.86 & 0.79 &0.87 & 0.91 \\ \midrule -fingerprint &20.34 & 2.99 & 5.73 & 7.59 \\ - & 0.93 & 0.87 & 0.95 & 0.96 \\ +fingerprint &20.34 & 2.99 & 6.00 &5.73 & 7.59 \\ + & 0.93 & 0.87 & 0.95 &0.95 & 0.96 \\ \midrule -goldhill & 19.59 & 6.88 & 7.84 & 9.63 \\ - & 0.67 & 0.82 & 0.87 & 0.88 \\ +goldhill & 19.59 & 6.88 & 8.02 &7.84 & 9.63 \\ + & 0.67 & 0.82 & 0.81 &0.87 & 0.88 \\ \midrule -lena & 19.92 & 8.07 & 9.22 & 11.88 \\ - & 0.60 & 0.84 & 0.88 & 0.93 \\ +lena & 19.92 & 8.07 & 8.37 &9.22 & 11.88 \\ + & 0.60 & 0.84 & 0.78 &0.88 & 0.93 \\ \midrule -man & 20.38 & 4.36 & 6.36 & 7.76 \\ - & 0.71 & 0.80 & 0.86 & 0.87 \\ +man & 20.38 & 4.36 & 6.49 &6.36 & 7.76 \\ + & 0.71 & 0.80 & 0.83 &0.86 & 0.87 \\ \midrule -mandrill & 19.34 & 1.00 & 3.04 & 5.41 \\ - & 0.77 & 0.69 & 0.83 & 0.88 \\ +mandrill & 19.34 & 1.00 & 4.20 &3.04 & 5.41 \\ + & 0.77 & 0.69 & 0.83 &0.83 & 0.88 \\ \midrule -peppers & 19.53 & 7.77 & 9.15 & 11.34 \\ - & 0.61 & 0.86 & 0.87 & 0.92 \\ +peppers & 19.53 & 7.77 & 8.66 &9.15 & 11.34 \\ + & 0.61 & 0.86 & 0.79 &0.87 & 0.92 \\ \midrule -stream & 20.35 & 2.88 & 5.00 & 5.99 \\ - & 0.80 & 0.78 & 0.87 & 0.88 \\ +stream & 20.35 & 2.88 & 4.97 &5.00 & 5.99 \\ + & 0.80 & 0.78 & 0.87 &0.87 & 0.88 \\ \midrule -zelda & 17.71 & 10.42 & 10.00 & 12.78 \\ - & 0.58 & 0.87 & 0.88 & 0.93 \\ +zelda & 17.71 & 10.42 & 11.13 &10.00 & 12.78 \\ + & 0.58 & 0.87 & 0.79 &0.88 & 0.93 \\ \bottomrule \end{tabular} -\caption{Comparaison image par image de la qualité de débruitage du filtre PI-PD hybride proposé par rapport à BM3D pris comme référence de qualité et à un moyenneur GPU 5$\times$5 pris comme référence de rapidité. Les paramètres du PI-PD sont $l=5$, $n=25$, $T_{max}=1$ et $T2_{max}=2$. La colonne 'noisy' donne les mesures relatives à l'image d'entrée corrompue par un bruit gaussien de moyenne nulle et d'écart type $\sigma=25$.} +\caption{Comparaison image par image de la qualité de débruitage des filtres PI-LD et PI-PD hybride proposé par rapport à BM3D pris comme référence de qualité et à un moyenneur GPU 5$\times$5 pris comme référence de rapidité. Les paramètres du PI-LD/PI-PD sont $n=5$, $l=25$, $T_{max}=1$ et $T2_{max}=2$. La colonne 'Bruitée' donne les mesures relatives à l'image d'entrée corrompue par un bruit gaussien de moyenne nulle et d'écart type $\sigma=25$. PI-LD s'exécute en 35~ms, PI-PD en 7,3~ms et BM3D en 4,3~s.} \label{tab-lniv-results} \end{table} \section{Extension aux images couleurs} \subsection{Expression du critère} -Considérons une image couleurs à 3 canaux RVB (Rouge, Vert et Bleu). La valeur $v_k$ observée au pixel $k$ est alors un vecteur à trois éléments. -Nous faisons ici l'hypothèse de canaux décorrelés, conduisant à une matrice de covariance diagonale de la forme $R=\sigma^2\mathbb{1}_3$ où $\sigma^2$ est la puissance du bruit gaussien perturbant les trois canaux, chaque canal pouvant être corrompu par un tirage de bruit particulier. -La probabilité de $v_k$ est alors +Considérons une image couleur à 3 canaux RVB (Rouge, Vert et Bleu). La valeur $v_k$ observée au pixel $k$ est alors un vecteur à trois éléments. +À notre connaissance, la définition des lignes de niveaux n'a pas été formellement étendue aux espaces colorimétriques à plusieurs canaux. +Nous faisons ici l'hypothèse (forte) de canaux décorrelés, conduisant à une matrice de covariance diagonale de la forme $R=\sigma^2\mathbb{1}_3$ où $\sigma^2$ est la puissance du bruit gaussien IID perturbant les trois canaux. +La vraisemblance de $v_k$ est alors $$P\left(v_k|R\right) = \left(\frac{1}{2\pi^{3/2}\sqrt{|R|}}\mathrm{e}^{-\frac{1}{2}\left(v_k-\mu\right)^TR^{-1}\left(v_k-\mu\right)}\right)$$ -Pour exprimer le critère GLRT de validation des allongements, nous procédons comme précédemment, c'est à dire en distinguant les deux hypothèses : +Pour exprimer le test GLRT de validation des allongements, nous procédons comme précédemment, c'est-à-dire en distinguant les deux hypothèses : \begin{enumerate} -\item le segment candidat $S^p$ prolonge effectivement l'\textit{isoline} $S^n$ : ils partagent donc la même valeur moyenne $\mu$ et la log-vraisemblance s'écrit +\item le segment candidat $S^p$ prolonge effectivement l'isoline $S^n$ : ils partagent donc la même valeur moyenne $\mu$ et la log-vraisemblance s'écrit \begin{align} \sum_{(i,j)\in S^p\cup S^n}-\frac{3}{2}log(2\pi)-\frac{1}{2}log(|R|)-\frac{1}{2}\left(v_{(i,j)}-\mu\right)^TR^{-1}\left(v_{(i,j)}-\mu\right)\nonumber \end{align} @@ -583,7 +586,7 @@ et \begin{align} \widehat{\mu} = \frac{1}{(n+p)}\sum_{(i,j)\in S^p\cup S^n}v_{(i,j)} \end{align} -\item le segment candidat ne prolonge pas l'\textit{isoline} $S^n$ : on distingue alors leur deux valeurs moyennes $\mu_p$ et $\mu_n$ et la log-vraisemblance s'écrit +\item le segment candidat ne prolonge pas l'isoline $S^n$ : on distingue alors leur deux valeurs moyennes $\mu_p$ et $\mu_n$ et la log-vraisemblance s'écrit \begin{align} \sum_{(i,j)\in S^n}-\frac{3}{2}log(2\pi)-\frac{1}{2}log(|R|)-\frac{1}{2}\left(v_{(i,j)}-\mu_n\right)^TR^{-1}\left(v_{(i,j)}-\mu_n\right) \nonumber\\ + \sum_{(i,j)\in S^p}-\frac{3}{2}log(2\pi)-\frac{1}{2}log(|R|)-\frac{1}{2}\left(v_{(i,j)}-\mu_p\right)^TR^{-1}\left(v_{(i,j)}-\mu_p\right) @@ -606,15 +609,18 @@ et \end{enumerate} -Le critère GLRT s'obtient par la soustraction des deux expressions de \eqref{eqlv1rgb} et \eqref{eqlv0rgb} : -$$T_{rvb} = 3(n+p)\left(-log\left(\widehat{{\sigma_1}^{2}}\right)+log\left(\widehat{{\sigma_0}^{2}}\right)\right) $$ -On notera $T_{rvb-max}$ la valeur de seuil au delà de laquelle on ne validera pas l'allongement de l'\textit{isoline}. +Le test GLRT s'obtient par la soustraction des deux expressions de \eqref{eqlv1rgb} et \eqref{eqlv0rgb} : +\begin{eqnarray} +T_{rvb} = 3(n+p)\left(-log\left(\widehat{{\sigma_1}^{2}}\right)+log\left(\widehat{{\sigma_0}^{2}}\right)\right) +\end{eqnarray} + +On notera $T_{rvb-max}$ la valeur de seuil au delà de laquelle on ne validera pas l'allongement de l'isoline. -\subsection{Résultats} +\subsection{Résultats - analyse} Nous avons retenu la base d'images de test tid2008 \cite{tid2008a} pour évaluer la qualité du traitement PI-PD sur les images couleurs. Cet ensemble d'images a été utilisé avec nombre d'algorithmes de débruitage et les résultats de mesure sont disponibles. -Chacune des 25 images de référence (non bruitées) a subit 4 niveaux de distorsion, pour 17 types de bruit différents. Pour nos expérimentations, nous avons selectionné les 25 images corrompues par un bruit gaussien RVB (type 2 dans tid2008) d'écart type $\sigma = 25$ (niveau 4 dans tid2008), où chaque canal RVB est perturbé par un tirage de bruit gaussien scalaire. La figure \ref{fig-lniv-tid2008ref} présente les vignettes des 25 images de référence, soit 24 images \textit{naturelles} et une image de synthèse. +Chacune des 25 images de référence (non bruitées) a subi 4 niveaux de distorsion, pour 17 types de bruit différents. Pour nos expérimentations, nous avons selectionné les 25 images corrompues par un bruit gaussien RVB (type 2 dans tid2008) d'écart type $\sigma = 25$ (niveau 4 dans tid2008), où chaque canal RVB est perturbé par un tirage de bruit gaussien IID. La figure \ref{fig-lniv-tid2008ref} présente les vignettes des 25 images de référence, soit 24 images \textit{naturelles} et une image de synthèse. \begin{figure}[ht] \centering @@ -652,11 +658,14 @@ Chacune des 25 images de référence (non bruitées) a subit 4 niveaux de disto \end{figure} Notre référence est ici encore l'implémentation BM3D dans sa variante couleurs (CBM3D) et nous avons choisi d'exprimer la qualité de débruitage au travers la valeur du PSNR-HVS-M (voir \cite{psnrhvsm}) qui est une extension du simple PSNR prenant en compte des caractéristiques structurelles de l'image. Les expérimentations décrites dans \cite{tid2008a} montrent en outre, que pour les perturbations de la catégorie \textit{noise} à laquelle appartient le type 2 qui nous intéresse, le PSNR-HVS-M présente les meilleures corrélations avec la perception humaine de la qualité, que ce soit au sens de Spearman ou de Kendall. -Comme pour les images en niveaux de gris, notre implémentation RVB intègre la réduction de l'effet \textit{marches d'escalier}, que nous avons adpaté à la couleur en choisissant la norme 2 comme mesure de distance dans l'espace RVB. Nous avons aussi expérimenté une variante employant la norme 1, avec des résultats moins satisfaisants. Cette étape améliore le rendu visuel mais représente cette fois une proportion plus importante du temps de calcul, en raison du calcul de la norme, plus coûteux. Sur les images de 512$\times$512, cela représente environ 1~ms. +Comme pour les images en niveaux de gris, notre implémentation RVB intègre la réduction de l'effet \textit{marches d'escalier}, que nous avons adapté à la couleur en choisissant la norme 2 comme mesure de distance dans l'espace RVB. Nous avons aussi expérimenté une variante employant la norme 1, avec des résultats moins satisfaisants. Cette étape améliore le rendu visuel mais représente cette fois une proportion plus importante du temps de calcul, en raison du calcul de la norme, plus coûteux. Sur les images de 512$\times$512, cela représente environ 1~ms, soit environ 25\% du temps de calcul. + +Le PI-PD en couleur s'exécute quant à lui à la même vitesse qu'en niveaux de gris, soit environ 4,0~ms ; c'est aussi le cas de CBM3D avec une moyenne de 4,3 secondes. Sur les 25 images de test, le gain moyen apporté par PI-PD s'élève à 2,84~dB (PSNR-HVS-M) contre 7,09~dB pour CBM3D, ce qui constitue indéniablement un échelon supérieur en terme de qualité, au prix d'un temps de calcul multiplié par 1000. -Le PI-PD en couleurs s'exécute quant à lui à la même vitesse qu'en niveaux de gris, soit environ 4,0~ms ; c'est aussi le cas de CBM3D avec une moyenne de 4,3~s. Sur les 25 images de test, le gain moyen apporté par PI-PD s'elève à 2,84~dB (PSNR-HVS-M) contre 7,09~dB pour CBM3D, ce qui constitue indéniablement un échelon supérieur en terme de qualité, au prix d'un temps de calcul multiplié par 1000. +L'ensemble des résultats de mesure est consigné dans le tableau \ref{tab-lniv-rvb} et deux exemples de résultats sont reproduits en figure \ref{fig-lnivrgb-ex} pour une des images naturelles ainsi que pour l'image de synthèse. Les valeurs des paramètres sont identiques pour toutes les images et ont été déterminées empiriquement par analyse systématique des résultats produits par les combinaisons permises dans les intervalles de 3 à 7 pour la taille $n$ des segments, de 25 à 70 pour la longueur maximale $l$ des isolines et de 1 à 10 pour le seuil GLRT $T_{rvb-max}$. Cette analyse extensive a mis en évidence la combinaison $n=4$, $l=48$ et $T_{rvb-max}=5$ comme permettant au PI-PD d'apporter les meilleurs résultats d'ensemble. Certaines des images, comme l'image de synthèse n°25, bénéficieraient d'un ajustement des paramètres, mais conscients de la contrainte que cela représente, nous avons choisi de faire prévaloir un réglage unique. -L'ensemble des résultats de mesure est consigné dans le tableau \ref{tab-lniv-rvb} et deux exemples de résultats sont reproduits en figure \ref{fig-lnivrgb-xe} pour une des images naturelles ainsi que pour l'image de synthèse. Les valeurs des paramètres sont identiques pour toutes les images et ont été déterminées empiriquement par analyse systématique des résultats produits par les combinaisons permises dans les intervalles de 3 à 7 pour la taille $l$ des segments, de 25 à 70 pour la longueur maximale $n$ des \textit{isolines} et de 1 à 10 pour le seuil GLRT $T_{rvb-max}$. Cette analyse extensive a mis en évidence la combinaison $l=4$, $n=50$ et $T_{rvb-max}=5$ comme permettant au PI-PD d'apporter les meilleurs résultats d'ensemble. Certaines des images, comme l'image de synthèse n°25, bénéficieraient d'un ajustement des paramètres, mais conscients de la contrainte que représente l'ajustement des paramètres, nous avons choisi de faire prévaloir un réglage unique. +Si le modèle que nous avons employé pour étendre le principe du PI-PD aux images en couleur nous a permis de vérifier que les performances pouvaient être conservées, il repose sur des hypothèses fortes et sous-entend (pour simplifier) qu'un seuil dans l'espace RVB est associé à un seuil sur chaque composante. Pour améliorer la qualité, une autre solution consisterait à détecter un seuil dans l'espace RVB lorsqu'un seuil est détecté dans au moins un des canaux. La transposition vers un autre espace colorimétrique employant une base de composantes décorrélées pourrait également améliorer les résultats. Ces pistes n'ont pas encore été explorées, mais représenteront nécessairement un surcoût calculatoire. + \label{fig-lnivrgb-ex} \begin{table}[H] \scriptsize @@ -693,7 +702,7 @@ L'ensemble des résultats de mesure est consigné dans le tableau \ref{tab-lniv- 25& 24.46& 24.62& 31.09\\ \bottomrule \end{tabular} -\caption{Comparaison image par image de la qualité de débruitage du filtre PI-PD RVB proposé par rapport à BM3D pris comme référence de qualité. Les paramètres du PI-PD sont $l=4$, $n=50$, $T_{rvb-max}=5$. La colonne 'noisy' donne les mesures relatives à l'image d'entrée corrompue par tirage de bruit gaussien sur chaque canal ( moyenne nulle, écart type $\sigma=25$).} +\caption{Comparaison image par image de la qualité de débruitage du filtre PI-PD RVB proposé par rapport à BM3D pris comme référence de qualité. Les paramètres du PI-PD sont $n=4$, $l=48$, $T_{rvb-max}=5$. La colonne 'noisy' donne les mesures relatives à l'image d'entrée corrompue par tirage de bruit gaussien sur chaque canal ( moyenne nulle, écart type $\sigma=25$).} \label{tab-lniv-rvb} \end{table} @@ -713,7 +722,9 @@ L'ensemble des résultats de mesure est consigné dans le tableau \ref{tab-lniv- \section{Conclusion} -L'algorithme PI-PD hybride permet de débruiter 19 images en haute définition à la seconde tout en réduisant de manière importante le niveau de bruit gaussien. -La démarche adoptée pour la conception a été de se baser sur des opérations élémentaires dont nous connaissions ou avions démontré l'efficacité sur GPU. Nous jugeons cet aspect essentiel pour la conception d'algorithmes GPU performants et robustes tant le débogage peu s'avérer délicat sur ces plateformes. Par ailleurs, il nous semble qu'il faille éviter de devoir systématiquement comparer les implémentations CPU et GPU pour en déduire un facteur d'accélération comme on le rencontre trop souvent. La plupart des algorithmes qui s'avèrent rapides sur GPU ne le sont vraisemblablement pas sur CPU et il est donc tout à fait illusoire de penser qu'il en existe une implémentation optimisée. Comparer alors une implémentation GPU performante avec son pendant CPU naïf ne présente aucun intérêt. La réciproque étant généralement vraie, nous avons choisi, en particulier en ce qui concerne le filtrage dont il est question ici, de chercher à assembler des blocs fonctionnels simples mais robustes performants avec l'objectif opérationnel de réduire la puissance de bruit. +L'algorithme PI-PD hybride permet de débruiter 19 images en haute définition à la seconde tout en réduisant de manière importante le niveau de bruit gaussien. +La démarche adoptée pour sa conception a été de se baser sur des opérations élémentaires dont nous connaissions ou avions démontré l'efficacité sur GPU. Nous jugeons ce principe essentiel pour la conception d'algorithmes GPU performants et robustes tant le débogage peut s'avérer délicat sur ces plateformes. Par ailleurs, il nous semble peu pertinent de systématiquement comparer les implémentations CPU et GPU pour en déduire un facteur d'accélération comme on le rencontre trop souvent. La plupart des algorithmes qui s'avèrent rapides sur GPU ne le sont vraisemblablement pas sur CPU et il est donc tout à fait illusoire de penser qu'il en existe une implémentation optimisée. Comparer alors une implémentation GPU performante avec son pendant CPU naïf ne présente aucun intérêt. La réciproque étant généralement vraie, nous avons choisi, en particulier en ce qui concerne le filtrage dont il est question ici, de chercher à assembler des blocs fonctionnels simples mais robustes et performants avec l'objectif opérationnel de réduire la puissance de bruit. + +L'algorithme et les résultats que nous avons détaillés dans ce chapitre ont été publiés dans le \textit{Journal of real-time image processing} dans un article intitulé \textit{Fast GPU-based denoising filter using isoline levels} \cite{perrotlniv}. -L'algorithme et les résultats que nous avons détaillés dans ce chapitre ont été publiés dans le \textit{Journal of real-time image processing} dans un article intitulé \textit{Fast GPU-based denoising filter using isoline levels} \cite{}. +% LocalWords: pénalisante estimateurs colorimétriques