lecture ch 1 à 3 27nov

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

\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. 

\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.]{
$
P_5 =
\begin{bmatrix}
(0,1)&(0,2)&(0,3)&(0,4)&(0,5)\\
(0,1)&(0,2)&(-1,3)&(-1,4)&(-1,5)\\
(0,1)&(-1,2)&(-1,3)&(-2,4)&(-2,5)\\
(-1,1)&(-1,2)&(-2,3)&(-3,4)&(-3,5)\\
(-1,1)&(-2,2)&(-3,3)&(-4,4)&(-5,5)\\
(-1,1)&(-2,1)&(-3,2)&(-4,3)&(-5,3)\\
(-1,0)&(-2,1)&(-3,1)&(-4,2)&(-5,2)\\
(-1,0)&(-2,0)&(-3,1)&(-4,1)&(-5,1)\\
\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$. }
\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}$).

\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$.}
\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.

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]
\]    

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]}
\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 

\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}
\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]}
\label{GL}
\end{eqnarray}
Que l'on peut développer, en remplaçant la probabilité 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}
\begin{eqnarray}
\displaystyle -\frac{n}{2}log\left(2\pi\right) - \frac{n}{2}log\left(\sigma^2\right) - \frac{n}{2}
\label{LL1}
\end{eqnarray}
où le vecteur des paramètres $(\mu_{S^n}, \sigma)$ est lui même obtenu par estimation au sens du maximum de vraisemblance.
$$
\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 \\
\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 \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.

\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.}
\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. 
À 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 :
\begin{enumerate}
\item Le segment $S^p$ {\bf est} une extension valide pour $S^n$. Ils forment donc une \textit{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)
\label{LLNP}
\end{eqnarray}
où \hspace{3cm}$\widehat{\sigma_1^2} = \displaystyle\frac{1}{n+p} \sum_{(i,j)\in S^nS^p} \left(v(i,j) - \widehat{\mu_{S^nS^p}}\right)^2 $.

\item Le segment $S^p$ {\bf n'est pas} une extension valide pour $S^n$. Les deux parties ont des valeurs moyennes distinctes ($\mu_{S^n}$, $\mu_{S^p}$) et la log-vraisemblance est alors la somme des log-vraisemblances des deux portions.
\begin{eqnarray}
\displaystyle -\frac{(n+p)}{2}\left(log\left(2\pi\right)+1\right) - \frac{n}{2}log\left(\widehat{\sigma_2}^2\right) - \frac{p}{2}log\left(\widehat{\sigma_2}^2\right)
\label{LLNP2}
\end{eqnarray}
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)
\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$.
 
\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.

\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$.

\begin{figure}[h] 
 \centering
\subfigure[Isoline comprenant deux segments $s_1$ et $s_2$.]{\label{pild:debut} \includegraphics{Chapters/chapter4/img/PI-LD_detail_3.jpg}}\quad
\subfigure[Premier segment évalué, associé au motif $p_{5,0}$.]{\label{pild:sub1} \includegraphics{Chapters/chapter4/img/PI-LD_detail_sub1.jpg}}\quad
\subfigure[Deuxième segment évalué, associé au motif $p_{5,1}$.]{\label{pild:sub2} \includegraphics{Chapters/chapter4/img/PI-LD_detail_sub2.jpg}}\\
\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.}
\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}.

\begin{figure}[ht]
  \centering
  \subfigure[airplane]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/airplane.png}}
  \subfigure[boat]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/boat.png}}\quad
  \subfigure[barbara]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/barbara.png}}\quad
  \subfigure[couple]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/couple.png}}\quad
  \subfigure[elaine]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/elaine.png}}\\
  \subfigure[fingerprint]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/fingerprint.png}}\quad
  \subfigure[goldhill]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/goldhill.png}}\quad
  \subfigure[lena]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/lena.png}}\quad
  \subfigure[man]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/man.png}}\quad
  \subfigure[mandrill]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/mandrill.png}}\\
  \subfigure[peppers]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/peppers.png}}\quad
  \subfigure[stream]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/stream.png}}\quad
  \subfigure[zelda]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/lansel/zelda.png}}\quad
\label{fig-lniv-imgslansel}
  \caption{Images non bruitées de la base d'images en niveaux de gris de S. Lansel.}
\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.

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}.

\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.}
\label{fig-lniv-histo-singe}
\end{figure}

 \begin{figure}[h] 
 \centering
\subfigure[Airplane]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_airplane_25_l3.png}}
\subfigure[Barbara]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_barbara_25_l3.png}}
\subfigure[Boat]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_boat_25_l3.png}}
\subfigure[Couple]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_couple_25_l3.png}}
\subfigure[Elaine]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_elaine_25_l3.png}}
\subfigure[Finger]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_finger_25_l3.png}}\\
\subfigure[Goldhill]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_goldhill_25_l3.png}}
\subfigure[Lena]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_lena_25_l3.png}}
\subfigure[Man]{\includegraphics[width=2cm]{Chapters/chapter4/img/hist_seg1_man_25_l3.png}}
\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.}
\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.

\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)$.

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.

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 :
\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
  \begin{eqnarray}
    C_x\left(Z\left(S_1\right)\right)= \sum_{(i,j)\in s_1} v(i,j)
    \label{cx}
  \end{eqnarray}
 et 
 \begin{eqnarray}
   C_{x^2}\left(Z\left(S_1\right)\right)= \sum_{(i,j)\in s_1} v(i,j)^2
   \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.
\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] 
 \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.}
\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}
$l \leftarrow$ taille segments\; 
$D \leftarrow$ nombre de motifs/directions\;
$I_n \leftarrow$ image d'entrée bruitée\;
$I_{n tex} \leftarrow I_n $\tcc*[r]{copie en texture}
$P_l \leftarrow$ kernel\_genPaths \tcc*[r]{génération de la matrice $P_l$}
$P_{l tex} \leftarrow P_l $\tcc*[r]{copie en texture}
$T_{max} \leftarrow$ seuil GLRT pour les allongements\;
$T2_{max} \leftarrow$ seuil GLRT pour la détection de bords\;
\end{algorithm}

\begin{algorithm}
  \SetKwComment{Videcomment}{}{}
\caption{\texttt{kernel\_precomp()} : génération des matrices $I_{\Theta}$ et $I_{\Sigma}$.}   
\label{algo-lniv-precomp}
\ForEach(\tcc*[f]{\textbf{en parallèle}}){pixel $(i,j)$}{
  $C_{x-best} \leftarrow  \displaystyle\sum_{(y,x)\in p_{l,0}(i,j)} I_{n tex}(i+y,j+x)$ \;
  $C_{x2-best} \leftarrow  \displaystyle\sum_{(y,x)\in p_{l,0}(i,j)} I_{n tex}^2(i+y,j+x)$ \;
  $\sigma_{best} \leftarrow$ écart type sur $p_{l,0}(i,j)$ \;
  %\Videcomment{}
  \tcc{pour chaque motif de segment}
  \ForEach{$d \in [1..D-1]$}{
    $C_x \leftarrow \displaystyle\sum_{(y,x)\in p_{l,d}(i,j)} I_{n tex}(i+y,j+x)$\;
    $C_{x2} \leftarrow \displaystyle\sum_{(y,x)\in p_{l,d}(i,j)} I_{n tex}^2(i+y,j+x)$\;
    $\sigma \leftarrow$ écart type sur $p_{l,d}(i,j)$\;
    \If(\tcc*[f]{sélection par MV}){$\sigma_d < \sigma_{best}$}{
      $C_{x-best} \leftarrow C_x$ \;
      $C_{x2-best} \leftarrow C_{x2}$ \; 
      $\Theta_{best} \leftarrow d$ \;
    }
  }
  $I_{\Sigma}(i,j) \leftarrow \left[ C_{x-best}, C_{x2-best}\right]$ \tcc*[r]{mémorisation}
  $I_{\Theta}(i,j) \leftarrow \Theta_{best}$ \tcc*[r]{dans $I_{\Theta}$ et $I_{\Sigma}$}
} 
\end{algorithm}

\begin{algorithm}[ht]
%\SetNlSty{textbf}{}{:}
\caption{\texttt{kernel\_PIPD()} : gestion du processus d'allongement.}   
\label{algo-lniv-pipd}
\ForEach(\tcc*[f]{\textbf{en parallèle}}){pixel $(i,j)$}{
  $(C_x^1, C_{x2}^1) \leftarrow z(i,j)$ \tcc*[r]{pixel de départ} 
  $(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$\; 
  $(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$ \;
  %
  \While{$GLRT(\sigma_1, \sigma_2, l_1, l) < T_{max}$}{
    $l_1 \leftarrow l_1 + l$ \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}
    $(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$\;
    \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$ \;
    }
  }
  $\widehat{I}(i, j) \leftarrow C_x^1/l_1$ \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).

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}. 

\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). 

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] 
 \centering
\subfigure[Image de référence non bruitée.]{\includegraphics[height=4cm]{Chapters/chapter4/img/airplane.png}}\quad
\subfigure[La région de 11$\times$11 pixels étudiée.]{\includegraphics[height=4cm]{Chapters/chapter4/img/zoom_windows_A.png}}
\caption{Situation de la région servant à illustrer le comportant du modèle PI-PD dans les zones à faible pente (LSR).}
\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é.

\begin{figure}[h] 
 \centering
\subfigure[Image de référence.]{\label{fig-lniv-lsr-tirages-a} \includegraphics{Chapters/chapter4/img/zoom_edge_ref2.jpg}}\quad
\subfigure[Image corrompue par le tirage de bruit $n^{\circ}1$]{\label{fig-lniv-lsr-tirages-b} \includegraphics{Chapters/chapter4/img/zoom_edge_bruit.jpg}}\quad
\subfigure[Image corrompue par le tirage de bruit $n^{\circ}2$]{\label{fig-lniv-lsr-tirages-c} \includegraphics{Chapters/chapter4/img/zoom_edge2_bruit.jpg}}\\
\subfigure[Directions déterminées par le PI-PD pour le tirage $n^{\circ}1$]{\label{fig-lniv-lsr-tirages-d} \includegraphics{Chapters/chapter4/img/zoom_edge1_2D_superpose3.jpg}}\quad
\subfigure[Directions déterminée par le PI-PD pour le tirage $n^{\circ}2$]{\label{fig-lniv-lsr-tirages-e} \includegraphics{Chapters/chapter4/img/zoom_edge2_2D_superpose2.jpg}}
\caption{Comportement du modèle PI-PD dans les zones de faible et à forte pente. On constate un manque de robustesse dans les zones à faible pente : les directions ne sont pas reproduites d'un tirage à l'autre, contrairement à celles de la zone de transition.}
\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 :
\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.
\end{enumerate}

\subsubsection{Le détecteur de zone à faible pente}

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$.

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. 

\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.}
  \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]
\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 :
$$
\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 \\
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                       
\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 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.
\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.  

\begin{figure}[h] 
 \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.}
\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 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. 

\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)$}{
  $\Theta \leftarrow 0$\tcc*[r]{Indice de la direction}
  $edgeCount \leftarrow 0$\;
  $sumEdge \leftarrow 0$\;
  $nH \leftarrow 5l+1$\;
  $nL \leftarrow 3l$\;
  \While{($\Theta < 32$) }{
    $sumH \leftarrow (I_{ntex}(i,j), I_{ntex}^2(i,j))$\;
    $sumL \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)$\;
    }
    \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)$\;
    }
    \If{($GLRT(sumH, nH, sumL, nL) > T2_{max}$)}{
      $edgeCount \leftarrow edgeCount + 1$\;
      $sumEdge \leftarrow sumH.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}
    }
    \If{($edgeCount == 1$)}{
      $\widehat{I}(i,j) \leftarrow \dfrac{(sumEdge)}{nH}$
    }
    \If{($edgeCount > 1$)}{
      $\widehat{I}(i,j) \leftarrow \widehat{I_{PIPD}}(i,j)$\tcc*[r]{PI-PD}
    }
  }
\end{algorithm}


\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.

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 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.
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{\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
  \subfigure{\includegraphics[width=3cm]{Chapters/chapter4/img/resultat/airplane_bm3d_zoom.jpg}}
  \caption{Comparaison des rendus des traitements comparés. Rangée du haut : les images complètes. Rangée du bas : Zooms sur une zone de l'ímage au dessus.}
  \label{fig-lniv-exempleresultat}
\end{figure}

\begin{table}[H]
\scriptsize
\centering
\begin{tabular}{crrrr}
\toprule
&\bf Temps de calcul &\bf Temps de transfert\\
&(ms)&(ms) \\
\midrule
Moyenneur         &  0.07&  0.15 \\
PI-PD hybride     &  9.00&  0.15 \\
BM3D              &  4300&  0.00 \\
\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.


\begin{table}[H]
\scriptsize
\centering
\begin{tabular}{crrrr}
\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}\\
\midrule
airplane &   19.49   &   6.90     &  8.97 & 11.39  \\
         &    0.58   &   0.84     &  0.88 &  0.93  \\
\midrule
barbara  &   20.04   &   2.72     & 4.22  & 10.56  \\
         &    0.70   &   0.76     &  0.83 &  0.94  \\
\midrule
boat     &   20.33   &   5.25     & 7.21  & 9.69   \\
         &    0.66   &   0.81     & 0.87  &  0.91  \\
\midrule
couple   &   20.28   &   4.97     & 7.05  & 9.49   \\
         &    0.69   &   0.79     & 0.87  &  0.91  \\
\midrule
elaine   &   19.85   &   8.86     & 9.09  & 10.75  \\
         &    0.59   &   0.86     & 0.87  &  0.91  \\
\midrule
fingerprint &20.34   &   2.99     & 5.73  &  7.59  \\
         &     0.93  &   0.87     & 0.95  &  0.96  \\
\midrule
goldhill &    19.59  &   6.88     & 7.84  &  9.63  \\
         &     0.67  &   0.82     & 0.87  &   0.88 \\
\midrule
lena     &    19.92  &   8.07     & 9.22  &  11.88 \\
         &     0.60  &   0.84     & 0.88  &   0.93 \\
\midrule
man      &    20.38  &   4.36     & 6.36  &  7.76  \\
         &     0.71  &   0.80     & 0.86  &   0.87 \\
\midrule
mandrill &    19.34  &   1.00     & 3.04  &  5.41  \\
         &     0.77  &   0.69     & 0.83  &   0.88 \\
\midrule
peppers  &    19.53  &   7.77     & 9.15  &  11.34 \\
         &     0.61  &   0.86     & 0.87  &   0.92 \\
\midrule
stream   &    20.35  &   2.88     & 5.00  &  5.99  \\
         &     0.80  &   0.78     & 0.87  &   0.88 \\
\midrule
zelda    &    17.71  &   10.42    & 10.00 &  12.78 \\  
         &     0.58  &   0.87     & 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$.}
\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
$$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 :
\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
  \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}
Soit
  \begin{align}
    \label{eqlv0rgb}
    -\frac{3}{2}(n+p)log\left(2\pi\right)-\frac{3}{2}(n+p)log\left(\widehat{{\sigma_0}^{2}}\right) -\frac{3}{2}(n+p)
  \end{align}
où 
\begin{align}
  \widehat{{\sigma_0}^{2}} = \frac{1}{3(n+p)}\sum_{(i,j)\in S^p\cup S^n}\left(v_{(i,j)}-\widehat{\mu}\right)^T\left(v_{(i,j)}-\widehat{\mu}\right)
\end{align}
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
  \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)
  \end{align}
Soit
  \begin{align}
    \label{eqlv1rgb}
    -\frac{3}{2}(n+p)log\left(2\pi\right)-\frac{3}{2}(n+p)log\left(\widehat{{\sigma_1}^{2}}\right) -\frac{3}{2}(n+p)
  \end{align}
où 
\begin{align}
  \widehat{{\sigma_1}^{2}} = \frac{1}{3(n+p)}\left(\sum_{(i,j)\in S^n}\left(v_{(i,j)}-\widehat{\mu_n}\right)^T\left(v_{(i,j)}-\widehat{\mu_n}\right)
+ \sum_{(i,j)\in S^p}\left(v_{(i,j)}-\widehat{\mu_p}\right)^T\left(v_{(i,j)}-\widehat{\mu_p}\right)\right)
\end{align}
et 
\begin{align}
  \widehat{\mu_n} = \frac{1}{n}\sum_{(i,j)\in S^n}v_{(i,j)}\\
\widehat{\mu_p} = \frac{1}{p}\sum_{(i,j)\in S^p}v_{(i,j)}
\end{align}

\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}.

\subsection{Résultats}

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. 

\begin{figure}[ht]
  \centering
  \subfigure[I01]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/0.png}}\quad
  \subfigure[I02]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/1.png}}\quad
  \subfigure[I03]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/2.png}}\quad
  \subfigure[I04]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/3.png}}\quad
  \subfigure[I05]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/4.png}}\\
 
  \subfigure[I06]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/5.png}}\quad
  \subfigure[I07]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/6.png}}\quad
  \subfigure[I08]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/7.png}}\quad
  \subfigure[I09]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/8.png}}\quad
  \subfigure[I10]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/9.png}}\\

  \subfigure[I11]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/10.png}}\quad
  \subfigure[I12]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/11.png}}\quad
  \subfigure[I13]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/12.png}}\quad
  \subfigure[I14]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/13.png}}\quad
  \subfigure[I15]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/14.png}}\\

  \subfigure[I16]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/15.png}}\quad
  \subfigure[I17]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/16.png}}\quad
  \subfigure[I18]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/17.png}}\quad
  \subfigure[I19]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/18.png}}\quad
  \subfigure[I20]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/19.png}}\\
  
  \subfigure[I21]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/20.png}}\quad
  \subfigure[I22]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/21.png}}\quad
  \subfigure[I23]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/22.png}}\quad
  \subfigure[I24]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/23.png}}\quad
  \subfigure[I25]{\includegraphics[width=2.5cm]{Chapters/chapter4/img/tid2008-petit/24.png}}
  \label{fig-lniv-tid2008ref}
\caption{Images non bruitées de la base tid2008.}
\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.  

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-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.  
\label{fig-lnivrgb-ex}
\begin{table}[H]
\scriptsize
\centering
\begin{tabular}{crrrr}
\toprule
\bf Image&\bf Noisy  &\bf PI-PD&\bf BM3D \\
         & \tiny{PSNR-HVS-M (dB)}& \tiny{gain (dB HVS-M)}/noisy& \tiny{gain (dB HVS-M)}/noisy\\
\midrule
1 &  23.91&   26.03&   31.02\\
2 &  23.39&   26.31&   28.96\\
3 &  23.31&   28.17&   33.04\\
4 &  23.21&   27.59&   31.30\\
5 &  24.52&   26.26&   30.63\\
6 &  23.91&   25.61&   29.45\\
7 &  23.68&   27.21&   32.66\\
8 &  24.51&   25.81&   30.75\\
9 &  23.70&   27.62&   33.26\\
10&   23.49&   27.16&   32.34\\
11&   23.95&   26.21&   30.63\\
12&   23.33&   27.26&   31.67\\
13&   24.17&   24.91&   29.45\\
14&   24.03&   25.83&   30.15\\
15&   23.55&   27.20&   30.63\\
16&   23.27&   27.11&   31.73\\
17&   23.74&   27.47&   32.18\\
18&   24.04&   25.67&   28.94\\
19&   24.20&   27.06&   31.45\\
20&   23.27&   26.58&   26.38\\
21&   23.75&   26.70&   31.41\\
22&   23.55&   26.22&   29.24\\
23&   23.48&   27.98&   32.11\\
24&   23.58&   26.40&   30.88\\
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$).}
\label{tab-lniv-rvb}
\end{table}
  
\begin{figure}[h]
  \centering
  \subfigure[Image I09 bruitée, PSNR-HVS-M=23,70~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I09_02_4.png}}\quad
  \subfigure[Image I09 filtrée par PI-PD, PSNR-HVS-M=27,62~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I09_02_4-4-50-5_out.png}}\quad
  \subfigure[Image I09 filtrée par CBM3D, PSNR-HVS-M=33,26~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I09_bm3d.png}}\\
  \subfigure[Image I25 bruitée, PSNR-HVS-M=24,46~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I25_02_4.png}}\quad
  \subfigure[Image I25 filtrée par PI-PD, PSNR-HVS-M=24,62~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I25_02_4-4-50-5_out.png}}\quad
  \subfigure[Image I25 filtrée par CBM3D, PSNR-HVS-M=31,09~dB.]{\includegraphics[width=4.5cm]{Chapters/chapter4/img/I25_bm3d.png}}\\
  \caption{Exemples de résultat de traitement par PI-PD RVB et par CBM3D pour deux images de la base tid2008 (une image naturelle et l'image de synthèse). Il peut être nécessaire de zoomer sur le document numérique pour visualiser les détails.}
  \label{fig-lnivgrb-ex}
\end{figure}

%\section{Extension au bruit multiplicatif de type Gamma}

\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 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{}.