v1.2 19 décembre

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

\message{ !name(chapter4.tex)}
\message{ !name(chapter4.tex) !offset(-2) }
\section{Introduction}
Le concept de ligne de niveau dans les images a été introduit dès 1975 par Matheron \cite{matheron75}, puis Caselles \textit{et al.} \cite{caselles97} l'ont exploité et proposé le cadre définissant les  \textit{images naturelles} comme les scènes photographiées, en intérieur ou en extérieur, à l'aide d'un appareil standard. Ces images vérifient alors l'hypothèse de gradient à valeurs bornées et peuvent être décomposées en un ensemble de lignes de niveaux.   
Bertaux \textit{et al.} ont plus récemment proposé un algorithme de réduction du \textit{speckle} dans les images éclairées en lumière cohérente en introduisant, pour les pixels de l'image observée, une contrainte  d'appartenance aux lignes de niveaux du modèle d'image non bruitée \cite{bertaux2004speckle}. L'image observée étant perturbée, les lignes de niveaux ne sont pas accessibles et il s'agit donc d'en estimer, localement  par morceaux, la valeur et la forme, en se basant sur un modèle 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 bruitée. 
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. Cependant, la résolution des images a 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 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$. 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)\\
(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}
$
}\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 $a=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$ ($a=5$).}
\end{figure}

La densité de probabilité des valeurs des niveaux de gris des pixels de  $S^n$ est supposée gaussienne de paramètres $\mu_{S^n}$ (moyenne) et $\sigma$ (écart type) inconnus. Les moyennes des niveaux de gris de chaque pixel sur $S^n$ sont inconnues et supposées indépendantes.

Soit $Z$ l'ensemble des niveaux de gris des pixels de $\omega$ et $\{\mu_{ij}\}_{\overline{S^n}}$ l'ensemble des valeurs moyennes des pixels de $\overline{S^n}$. On peut écrire la probabilité 
\[
P[Z|S^n, \mu_{S^n}, \{\mu_{ij}\}_{S^n}, \sigma]
\]    

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 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.}
\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}
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 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}
\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 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.
À la place, nous effectuons une sélection à chaque étape d'allongement : on évalue l'ensemble des 32 allongements possibles et si au moins un des motifs est accepté, on retient l'\textit{isoline} ayant la meilleure vraisemblance. Ce processus est répété tant qu'au moins un motif représente un allongement valide.

\begin{figure}[h]
\center
\includegraphics[height=5cm]{Chapters/chapter4/img/exemple_extension_1.jpg}
\caption{\label{fig-lniv-allongement}Allongement du segment $S^n$. Deux candidats $S^{p'}$ et $S^{p''}$ sont évalués au travers du 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 isoline précédemment validée et $S^p$ un segment connecté à $S^n$ de telle sorte qu'il représente un allongement potentiel de $S^n$. Une situation de cette nature est représentée à  la figure \ref{fig-lniv-allongement} avec un premier segment valide et deux candidats $S^{p'}$ et $S^{p''}$.
À gauche de la figure est reproduite une petite zone d'image réelle, suffisamment grossie pour permettre de bien individualiser les pixels et  à laquelle ont été superposés les trois segments en question. Les deux relevés de la partie droite montrent quant à eux l'évolution des valeurs des niveaux de gris des pixels des deux isolines possibles que sont $S^nS^{p'}$ et $S^nS^{p''}$. On a également identifié les différents sous-ensembles de pixels $S^n$, $S^{p'}$ et $S^{p''}$ ainsi que les valeurs moyennes de chacun. 
À la lecture de ces deux représentations, on peut aisément imaginer que $S^{p'}$ ne soit pas retenu comme extension valide de $S^n$, au contraire de $S^{p''}$.

Pour formaliser ce que notre intuition semble nous dicter dans l'exemple précédent, 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 tous deux une isoline $S^nS^p$.\\
Dans ce cas et par hypothèse, la valeur moyenne des niveaux de gris est définie sur $S^nS^p$ et vaut $\mu_{S^nS^p}$. D'après \eqref{LL2}, la log-vraisemblance est alors donnée par
\begin{eqnarray}
\displaystyle -\frac{(n+p)}{2}\left(log\left(2\pi\right)+1\right) - \frac{(n+p)}{2}log\left(\widehat{\sigma_1}^2\right)
\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 (\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é 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. 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 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 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 (\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[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 isoline comportant initialement deux segments $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.}
\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. 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 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 qu'il y a une proportion relativement faible d'isolines optimales dont le premier segment s'écarte notablement de celui qui présente la plus grande vraisemblance en l'absence d'allongement. 
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, l'écart est nul.}
\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'isoline optimale et celle de l'isoline sélectionnée, pour les images de l'ensemble de test de S. Lansel. La répartition des erreurs est semblable dans toutes ces images, mais également dans toute image naturelle.}
\label{fig-lniv-histo-autres}
\end{figure}

On observe également que les pixels pour lesquels la sélection du premier segment n'est pas robuste sont situés dans les zones de l'image ne contenant pas de forts gradients de niveaux de gris, 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. 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, 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 (\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
  \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'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] 
 \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 isoline comprenant initialement 2 segments. la longueur des segments est $a=5$. Le procédé se répète jusqu'à ce que le test GLRT échoue.}
\label{fig-lniv-pipd}
\end{figure}

\begin{algorithm}[htb]
\caption{Initialisations du modèle PI-PD, en mémoire du GPU.}   
\label{algo-lniv-init}
$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 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$ \;
  %
  \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 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$ \;
    }
  }
  $\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à é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$\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  (\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] 
 \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 des isolines, modulo $\pi$, débutant en chaque pixel de la fenêtre. On observe que la détermination de la direction est robuste dans la bande de transition entre les LSR, alors que pour les LSR elles-mêmes, on constate une grande variabilité.

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

Ainsi, dans les 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ées LSR par le détecteur et le PI-PD partout ailleurs.
\item de n'appliquer le moyenneur que sur les pixels appartenant à la zone LSR lorsque la fenêtre du détecteur se trouve à cheval sur deux zones de types différents.
\end{enumerate}

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

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

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 \\
et\\
\widehat{\sigma_4}^2 = \displaystyle\frac{1}{8l+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 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 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()}}   
\label{algo-lniv-detecteur}
\ForEach(\tcc*[f]{\textbf{en parallèle}}){pixel $(i,j)$}{
  $\Theta \leftarrow 0$\tcc*[r]{Indice de la direction}
  $edgeCount \leftarrow 0$\;
  $sumEdge \leftarrow 0$\;
  $nT \leftarrow 5l+1$\;
  $nB \leftarrow 3l$\;
  \While{($\Theta < 32$) }{
    $sum_{\Theta} \leftarrow \left(\displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I_{n tex}(i+y,j+x), \displaystyle\sum_{(y,x)\in P_{l,\alpha}(i,j)} I^2_{n tex}(i+y,j+x)\right)$ \;
    $\Theta \leftarrow \Theta + 4$\;
  }
  \While{($\Theta < 32$) }{
    $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$)}{
      $sumT \leftarrow sumT + sum_{\Theta}$\; 
    }
    \For{($\alpha=\Theta+20$ to $\alpha=\Theta+28$ par pas de $4$)}{
      $sumB \leftarrow sumB + sum_{\Theta}$ \;
    }
    \If{($GLRT(sumT, nT, sumB, nB) > T2_{max}$)}{
      $edgeCount \leftarrow edgeCount + 1$\;
      $sumEdge \leftarrow sumT.x$\;
    }
    $\Theta \leftarrow \Theta + 4$\;
  }
    \tcc{niveau de gris de l'isoline}
    \If{($edgeCount == 0$)}{
      $\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)}{nT}$
    }
    \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 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 $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 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 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 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 . 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}.

\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 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
  \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&  \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 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 non paginée 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}
\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}\\
\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 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$.}
\label{tab-lniv-results}
\end{table}


\section{Extension aux images couleurs}
\subsection{Expression du critère}
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 :
\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}
    \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'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'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 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
  \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 adapté à la couleur en choisissant la norme 2 comme mesure de distance dans l'espace RVB. Nous avons aussi expérimenté une variante employant la norme 1, avec des résultats moins satisfaisants. Cette étape améliore le rendu visuel mais représente cette fois une proportion plus importante du temps de calcul, en raison du calcul de la norme, plus coûteux. Sur les images de 512$\times$512, cela représente environ 1~ms, soit environ 25\% du temps de calcul.  

Le PI-PD en couleur s'exécute quant à lui à la même vitesse qu'en niveaux de gris, soit environ 4.0~ms ; c'est aussi le cas de CBM3D avec une moyenne de 4.3 secondes. Sur les 25 images de test, le gain moyen apporté par PI-PD s'élève à 2.84~dB (PSNR-HVS-M) contre 7.09~dB pour CBM3D, ce qui constitue indéniablement un échelon supérieur en terme de qualité, au prix d'un temps de calcul multiplié par 1000. 

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

% LocalWords:  pénalisante

\message{ !name(chapter4.tex) !offset(-720) }