\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 un terme d'attache aux données lié 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 pré établi.
+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 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 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 demeuraient 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ésolution des images a traiter à, elle aussi, crû dans des proportions comparables, laissant les termes du compromis qualité/performances inchangés.
+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}
\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 sélectionné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.
+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. Le choix d'intégrer ou non d'autres segments à l'\textit{isoline} et la sélection des segments à intégrer sont effectués par évaluation d'un critère de vraisemblance généralisée dont l'obtention est détaillée dans la suite.
\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 premiers segments candidats de longueur 5 pixels. 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.]{
-$
+\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)\\
\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}
-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}$).
+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 côté $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 $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é
\[
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}
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}
-On peut développer cette expression 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 développant, 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}
$$
\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.
$$
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 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 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.
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 critère GLRT (\textit{Generalized Likelihood Ratio Test})
\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é.\\
+où $T_{max}$ est un seuil arbitrairement fixé de sorte à produire des résultats visuels et chiffrés satisfaisant.
Un allongement de $S^n$ par $S^p$ est validé si $T(S^n, S^p, T_{max}) > 0$.
\section{Modélisation des isolines pour l'implémentation parallèle sur GPU}
-Les 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.
+Les isolines sont construites segment après segment. Cela permet 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.
Il nous est également apparu pertinent de limiter le nombre de segments candidats, ce qui 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'isoline la meilleure ait 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 nous conduisons en parallèle les allongements des 32 isolines, chacune ayant l'un des motifs permis comme premier segment (voir figure \ref{fig-lniv-p5q1}).
-\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 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 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 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 être composée de segments qui se croisent.
+\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 nous conduisons en parallèle les allongements des 32 isolines, chacune ayant l'un des motifs permis comme premier segment (voir figure \ref{fig-lniv-p5q1}).
+\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 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.
+\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}
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 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.
+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 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 $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
\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 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'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 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'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{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. La présente variante y est repérée PI-LD. 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
\end{figure}
-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 divergentes à l'exécution du kernel principal qui sont sérialisées 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 isolines construites nous montre la proportion relativement faible d'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 optimales dont le premier segment s'écarte notablement de celui sélectionné 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 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}. 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 à 2 (en 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'isoline optimale et celle de l'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'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}
\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 l'impossibilité d'identifier une direction privilégiée dans ces régions.
\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. En effet, dans ce cas, étendre une 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'isoline débutant en $(i_t, j_t)$.
+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 de $32^q$ à seulement {\bf 160} évaluations par pixel, soit un total de $\mathbf{1,7.10^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 du segment initial par maximum de vraisemblance des phases d'allongements successifs soumis au seul test GLRT.
+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()} 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
\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'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()} é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 critère 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, ce qui 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.
+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 critère 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 isoline par la méthode PI-PD. \begin{figure}[h]
\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 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}
$(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$\;
+ $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 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à é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 2, on utilise chaque segment motif et le segment parallèle situé immédiatement avant (au sens trigonométrique), pour l'épaisseur 3, on ajoute le segment parallèle situé immédiatement après le motif. 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 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^2$ pixels (Xeon quad core E31245\@3.3GHz, 8Go RAM).
+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).
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 LSR (Low Slope Regions, régions à faible pente).
+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]
\label{fig-lniv-lsr-tirages}
\end{figure}
-Ainsi, 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. Le meilleur estimateur dans une zone LSR étant la valeur moyenne, nous proposons donc, à la place :
\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ées LSR par le détecteur et le PI-PD partout ailleurs.
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, 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. 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 directions $\Theta_{4i} = 4i\frac{\pi}{4}$.
+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 critère 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 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.
\begin{figure}[ht]
\centering
\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 $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.}
+ \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$
\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}
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 :
\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 est appliqué et les zones moyennées.
+\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 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 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.
\end{itemize}
\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 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.
+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.
\begin{algorithm}[ht]
\caption{Détecteur de zones à faible pente (LSR) \texttt{kernel\_LSR\_detector()}}
$\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$) }{
- $sumH \leftarrow (I_{ntex}(i,j), I_{ntex}^2(i,j))$\;
- $sumL \leftarrow (0, 0)$\;
+ $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$) }{
+ $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}
\section{Résultats}
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 $l=5$ pour la longueur des segments avec un maximum de 5 segments.
+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, 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.
-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 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 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 que ce dernier.
+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 128, soit environ 9~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 475 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}.
\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
\midrule
Moyenneur & 0.07& 0.15 \\
PI-PD hybride & 9.00& 0.15 \\
-BM3D & 4300& 0.00 \\
+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 (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. 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.
\begin{table}[H]
\scriptsize
\centering
-\begin{tabular}{crrrr}
+\begin{tabular}{crrrrr}
\toprule
-\bf Image&\bf Bruitée &\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 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-PD sont $l=5$, $n=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$.}
+\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 9~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.
+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.
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
$$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 critère 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'isoline $S^n$ : ils partagent donc la même valeur moyenne $\mu$ et la log-vraisemblance s'écrit
\begin{align}
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 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.
\begin{figure}[ht]
\centering
\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 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.
+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.
-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 $l$ des segments, de 25 à 70 pour la longueur maximale $n$ des 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.
+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.
\label{fig-lnivrgb-ex}
\begin{table}[H]
\scriptsize
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}
\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 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 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 et performants avec l'objectif opérationnel de réduire la puissance de bruit.
+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 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}.
