\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.
+Bertaux \textit{et al.} ont plus récemment proposé un algorithme de réduction du \textit{speckle} dans les images éclairées en lumière cohérente en introduisant, pour les pixels de l'image observée, une contrainte d'appartenance aux lignes de niveaux du modèle d'image non bruitée \cite{bertaux2004speckle}. L'image observée étant perturbée, les lignes de niveaux ne sont pas accessibles et il s'agit donc d'en estimer, localement par morceaux, la valeur et la forme, en se basant sur un modèle de forme pré établi.
Pour un pixel dont on cherche à estimer la valeur du niveau de gris, la contrainte d'appartenance à une ligne de niveau demeure locale, avec cependant un voisinage de forme et de taille (en nombre de pixels) variables en fonction des propriétés de l'image bruitée dans la zone concernée.
Ce voisinage, dont la forme, l'étendue et le niveau de gris sont déterminés par maximum de vraisemblance, appelé une \textit{isoline}, est une estimation locale de la ligne de niveau à laquelle appartient le pixel concerné.
Cette technique a montré qu'elle permettait de réduire très significativement le niveau de bruit tout en préservant les contours des objets. Elle s'est en revanche averée gourmande en ressources, ce qui a initialement conduit les auteurs à réduire la résolution de calcul des \textit{isolines} par application d'un maillage sur l'image bruitée.
\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.
+Pour chaque pixel de l'image d'entrée, le premier segment est choisi comme celui présentant la meilleure vraisemblance parmi les 32 possibles. Un test statistique GLRT (\textit{Generalized Likelihood Ratio Test}) détermine la pertinence d'intégrer ou non un segment supplémentaire à l'\textit{isoline} et le segment effectivement sélectionné pour l'allongement est le meilleur au sens du maximum de vraisemblance, parmi tous les segements satisfaisant au test GLRT.
\begin{figure}[h]
\subfigure[Les 8 premières lignes de la table $P_5$. Les éléments sont les positions relatives des pixels de chaque motif par rapport au pixel central.]{$
\caption{\label{fig-lniv-p5q1}Détail des motifs et de leur représentation interne, pour la taille $a=5$. }
\end{figure}
-\subsubsection{Isolines à un seul segment}
+\subsubsection{Détermination du premier segment}
-Pour chacun des pixels $(i,j)$ de l'image corrompue, on calcule la vraisemblance associée à chaque segment candidat de la table $P_{n-1}$ dans la région carrée $\omega$ centrée en $(i,j)$ et de côté $2n-1$. La région $\omega$ est l'union des deux sous régions $S^n$ et $\overline{S^n}$ telles que $S^n$ décrit le segment candidat à évaluer comme un ensemble de $n$ pixels de coordonnées $(i_q,j_q)$ où $q\in [0..n[$.
+Nous avons arbitrairement choisi un modèle simple où, au sein de la région carrée $\omega$ centrée en $(i,j)$ et de côté $2n-1$, les pixels du segment de ligne de niveau composent la région $S^n$ et possèdent la même valeur de niveu de gris. Les pixels de $\omega$ n'appartenant pas à $S^n$ composent la région $\overline{S^n}$ et aucun modèle ne leur est attaché. La région $\omega$ est ainsi l'union des des sous régions $S^n$ et $\overline{S^n}$ et le segment à évaluer est l'ensemble des $(a+1)$ pixels de coordonnées $(i_q,j_q)$ où $q\in [0..n[$.
+
+Pour chacun des pixels $(i,j)$ de l'image observée, on calcule la vraisemblance associée à chaque segment candidat de la table $P_{n-1}$ dans la région $\omega$.
La figure \ref{fig-lniv-regions} montre cette répartition pour $a=5$ et le motifs $p_{5,3}$.
\begin{figure}[h]
Soit $Z$ l'ensemble des niveaux de gris des pixels de $\omega$ et $\{\mu_{ij}\}_{\overline{S^n}}$ l'ensemble des valeurs moyennes des pixels de $\overline{S^n}$. On peut écrire la probabilité
\[
-P[Z|S^n, \mu_{S^n}, \{\mu_{ij}\}_{S^n}, \sigma]
+P[Z|S^n, \mu_{S^n}, \{\mu_{ij}\}_{\{overline{S^n}}, \sigma]
\]
qui se développe comme suit, en distinguant les contributions de $S^n$ et $\overline{S^n}$
\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}$.
+Or, d'après notre modèle, sur $\overline{S^n}$, les niveaux de gris $z(i,j)$ sont aussi les estimations $\widehat{\mu_{ij}}$ de leurs moyennes $\mu_{ij}$.
Le second terme de l'expression \eqref{LL2} devient donc
\begin{eqnarray}
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.
+Pour un bruit indépendant, la variance varie de manière inversement proportionnelle à la racine du nombre d'échantillons.
+L'objectif est alors d'améliorer la force du filtrage en intégrant plus de valeurs de niveaux de gris dans le calcul de la moyenne qui deviendra la valeur de sortie filtrée.
+Pour cela nous permettons à chaque isoline, comportant initialement un seul segment, d'être prolongée par d'autres segments, chaque allongement faisant l'objet d'une validation par un test (GLRT) issu de la théorie de la détection statistique (voir \cite{van2004detection}).
L'évaluation de l'ensemble des isolines pouvant être construite sur ce modèle présente un coût prohibitif en temps de calcul et l'idée en a donc été abandonnée.
À la place, nous effectuons une sélection à chaque étape d'allongement : on évalue l'ensemble des 32 allongements possibles et si au moins un des motifs est accepté, on retient l'\textit{isoline} ayant la meilleure vraisemblance. Ce processus est répété tant qu'au moins un motif représente un allongement valide.
\begin{figure}[h]
\center
\includegraphics[height=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.}
+\caption{\label{fig-lniv-allongement}Allongement du segment $S^n$. Deux candidats $S^{p'}$ et $S^{p''}$ sont évalués au travers du test GLRT de l'équation \eqref{GLRT} que seul $S^{p''}$ s'avère satisfaire. a) Représentation dans le plan de l'image. b) Évolution des niveaux de gris en fonction de la position des pixels dans les lignes brisées ainsi formées.}
\end{figure}
Soit $S^n$ une 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 :
+Pour formaliser, nous comparons les log-vraisemblances des deux situations suivantes :
\begin{enumerate}
\item Le segment $S^p$ {\bf est} une extension valide pour $S^n$. Ils forment donc 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
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})
+La différence entre \eqref{LLNP} et \eqref{LLNP2} nous donne l'expression du test GLRT
\begin{eqnarray}
T(S^n, S^p, T_{max}) = T_{max}- (n+p)\left[log\left(\widehat{\sigma_1}^2\right) - log\left(\widehat{\sigma_2}^2\right) \right]
\label{GLRT}
\end{eqnarray}
-où $T_{max}$ est un seuil arbitrairement fixé 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$.
+Dans cette étude, la valeur du seuil de détection $T_{max}$ est déterminée empiriquement, de sorte à produire des résultats visuels et chiffrés satisfaisants. Une détermination rigoureuse de la valeur adéquate de $T_{max}$ est envisageable, en fonction de la probabilité de fausse alarme visée, par le calcul de la PDF de $T$, mais dépasse le cadre de nos travaux. L'allongement de $S^n$ par $S^p$ est finalement validé si $T(S^n, S^p, T_{max}) > 0$.
\section{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.
+Les isolines sont construites segment après segment afin de pouvoir approcher des formes courbes. La validité d'un segment et son éventuelle sélection sont soumises au tests statistiques décrits dans le paragraphe précédent.
+La solution la plus évidente, que nous qualifierons de \textit{globale}, consiste à soumettre l'ensemble des isolines possibles aux tests statistiques pour sélectionner la meilleure. Toutefois, évaluer systématiquement les 32 motifs pour chaque extension peut rendre l'algorithme très coûteux. En effet, si $q$ est le nombre de segments maximum autorisés pour une isoline, le nombre d'évaluations à effectuer se monte à $32^q$ par pixel. Cela représente par exemple un total de $\mathbf{3,5.10^{13}}$ évaluations pour des isolines de $q=5$ segments dans une image de 1024$\times$1024 pixels.
+
+De plus, cette solution génère potentiellement des artefacts pour le cas où des isolines s'enrouleraient sur elles-mêmes au delà du simple rebouclage, ou bien feraient des aller-retours entre les deux extrémités d'un segment, ou encore présenteraient des segments qui se croisent.
+
+Pour toutes ces raisons, nous avons limité le nombre de segments candidats à chaque étape d'allongement, mais n'effectuons pas de sélection directe du premier segment. Nous conduisons en parallèle les allongements de 32 isolines, chacune ayant l'un des motifs de $P_a$ comme premier segment (voir figure \ref{fig-lniv-p5q1}).
+Cela permet d'éviter une erreur sur la direction primaire potentiellement très pénalisante pour la qualité du traitement, tout en conservant une quantité d'évaluations réalisable dans un temps assez court.
+% \begin{enumerate}
+% \item la sélection du premier segment est cruciale mais il n'est pas prouvé que la meilleure isoline soit celle qui a pour premier segment celui qui a été effectivement sélectionné en premier. Une telle erreur sur la direction primaire peut s'avérer très pénalisante pour la qualité du traitement. C'est pourquoi
+% \item
+% \item permettre à tout allongement de se faire dans chacune des 32 directions risque de générer des isolines oscillant entre les deux extrémités de l'un de ses segments, ou bien s'enroulant sur elles-même au delà du simple rebouclage.
+% \item une ligne de niveau ne peut pas se couper, donc une isoline ne peut pas être composée de segments qui se croisent.
+% \end{enumerate}
+La solution retenue consiste à limiter la déviation angulaire pouvant résulter de toute procédure d'allongement. Nous notons $\Delta d_{max}$ l'écart maximal toléré entre les indices des motifs de deux segments successifs.
Le choix d'une valeur de $\Delta d_{max}$ adaptée dépend de la taille des segments ainsi que du nombre maximal de segments que peut comporter une isoline.
L'autre conséquence de cette limitation est la diminution du nombre total d'évaluations nécessaires. Si $\Delta d_{max} = 2$, le nombre d'évaluations effectuées dans l'exemple du point 2 passe ainsi à $1024^2\times 32\times 5^{q-1} = \mathbf{2,0.10^{10}}$ soit 1500 fois moins (avec $q=5$).
\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 $a=5$ et $\Delta d_{max}=2$. Chaque segment évalué est soumis au critère GLRT. Si au moins un des segments présente un test GLRT positif, alors l'allongement est réalisé avec le segment qui forme l'isoline la plus vraisemblable.}
+\caption{Processus de sélection lors de l'allongement d'une isoline comportant initialement deux segments $s_1$ et $s_2$. Dans cet exemple $a=5$ et $\Delta d_{max}=2$. Chaque segment évalué est soumis au test GLRT. Si au moins un des segments présente un test GLRT positif, alors l'allongement est réalisé avec le segment qui forme l'isoline la plus vraisemblable.}
\label{fig-lniv-pild}
\end{figure}
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 sélectionné l'absence d'allongement, c'est-à-dire par PI-LD avec $q=1$.
-L'exemple représentatif de la figure \ref{fig-lniv-histo-singe} montre l'histogramme des différences constatée pour l'image du singe. Les autres images de l'ensemble de test fournissent des histogrammes très semblables qui sont reproduits en petit format à la figure \ref{fig-lniv-histo-autres}. On y observe que pour environ 60\% des pixels de l'image, il y a correspondance des directions et que pour 80\% des pixels, l'écart angulaire reste inférieur à 2 (en indices des motifs).
+Une analyse plus poussée des isolines construites nous montre qu'il y a une proportion relativement faible d'isolines déterminées par la méthode \textit{globale} dont le premier segment s'écarte notablement de celui sélectionné en l'absence d'allongement, c'est-à-dire par PI-LD avec $q=1$.
+L'exemple représentatif de la figure \ref{fig-lniv-histo-singe} montre l'histogramme des écarts angulaires constatés entre les deux estimateurs, pour l'image du singe en version bruitée ($\sigma=25$). Les autres images de l'ensemble de test fournissent des histogrammes très semblables qui sont reproduits en petit format à la figure \ref{fig-lniv-histo-autres}. On y observe que pour environ 60\% des pixels de l'image, il y a correspondance des directions et que pour 80\% des pixels, l'écart angulaire reste inférieur à 22,5 degrés (soit 2 indices des motifs).
\begin{figure}[h]
\centering
\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.
+On observe également que les pixels pour lesquels la sélection du premier segment n'est pas robuste sont situés dans les zones de l'image ne contenant pas de forts gradients de niveaux de gris, ce qui est cohérent avec la difficulté d'identifier une direction privilégiée dans ces régions.
\subsection{Isolines à segments 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)$.
\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.
+\item \texttt{kernel\_PIPD()} évalue les allongements successifs, qui ne nécessitent plus de sélection par maximum de vraisemblance, mais uniquement la validation par GLRT. Les données nécessaires à l'évaluation du test GLRT sont regroupées, outre dans l'image d'entrée, dans les matrices $P_d$, $I_{\Theta}$ et $I_{\Sigma}$ et ne sont donc plus à calculer à ce stade. Cela permet d'envisager des performances en hausse par rapport à la solution PI-LD. L'algorithme \ref{algo-lniv-pipd} fournit les détails de l'implémentation de ce kernel.
\end{enumerate}
Les schémas de la figure \ref{fig-lniv-pipd} illustrent les étapes décrites ci-dessus de l'allongement d'une isoline par la méthode PI-PD. \begin{figure}[h]
$\widehat{I}(i, j) \leftarrow C_x^1/l$ \tcc*[r]{niveau de gris en sortie}
\end{algorithm}
-Le processus d'allongement du modèle PI-PD est également soumis aux restrictions sur les oscillations et retours en arrière des segments, déjà énoncées pour le modèle PI-LD. Par ailleurs, nous lui avons ajouté la possibilité de gérer des segments plus épais, composés de 2 ou 3 segments parallèles aux motifs décrits par la matrice $P_d$. Pour l'épaisseur 2, on utilise chaque segment motif et le segment parallèle situé immédiatement avant (au sens trigonométrique), pour l'épaisseur 3, on ajoute le segment parallèle situé immédiatement après le motif. Cela a pour effet d'intégrer plus de pixels dans les calculs statistiques et d'augmenter en conséquence les gains sur le PSNR, en particulier pour traiter des images de grandes dimensions qui ne contiendraient pas de \textit{trop petits} détails que l'épaisseur des isolines risquerait de flouter.
+Le processus d'allongement du modèle PI-PD est également soumis aux restrictions sur les retours en arrière des segments, déjà énoncées pour le modèle PI-LD. Par ailleurs, nous lui avons ajouté la possibilité de gérer des segments plus épais, composés de 2 ou 3 segments parallèles aux motifs décrits par la matrice $P_d$. Pour l'épaisseur 3, on utilise chaque segment motif ainsi que les deux segments parallèles obtenus par translation et situés de part et d'autre du segment motif. Pour l'épaisseur 2, on ne conserve qu'un seul de ces deux segments ajoutés. Cet épaississement a pour effet d'intégrer plus de pixels dans les calculs statistiques et d'augmenter en conséquence les gains sur le PSNR, en particulier pour traiter des images de grandes dimensions qui ne contiendraient pas de \textit{trop petits} détails que l'épaisseur des isolines risquerait de flouter.
Cette possibilité rend notre solution encore plus versatile que la référence BM3D dont les temps de calcul s'avèrent prohibitifs sur des images de grandes dimensions, avec par exemple plus de 5 minutes pour 4096$\times$4096 pixels (Xeon quad core E31245\@3.3GHz, 8Go RAM).
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}.
\label{fig-lniv-lsr-tirages}
\end{figure}
-Ainsi, dans les LSR, l'application du modèle PI-PD n'a que peu de sens, mais la quête de performance nous interdit d'y appliquer par exemple le modèle PI-LD décrit précédemment. Le meilleur estimateur dans une zone LSR étant la valeur moyenne, nous proposons donc, à la place :
+Ainsi, dans les LSR, l'application du modèle PI-PD n'a que peu de sens, mais la quête de performance nous interdit d'y appliquer, par exemple, le modèle PI-LD décrit précédemment. Considérant que le bruit est additif gaussien et que nous faisons une approximation \textit{plane} de la surface définie par les niveaux de gris des pixels à l'intérieur de la fenêtre $\omega$, le meilleur estimateur dans une zone LSR est la valeur moyenne. Nous proposons donc :
\begin{enumerate}
-\item d'identifier les zones à faible pente en concevant un kernel détecteur (\texttt{kernel\_LSR\_detector()}).
+\item d'identifier les zones à faible pente en concevant un kernel détecteur de bords (\texttt{kernel\_edge\_detector()}).
\item d'appliquer un simple filtre moyenneur dans les zones désignées LSR par le détecteur et le PI-PD partout ailleurs.
\item de n'appliquer le moyenneur que sur les pixels appartenant à la zone LSR lorsque la fenêtre du détecteur se trouve à cheval sur deux zones de types différents.
\end{enumerate}
-\subsubsection{Le détecteur de zone à faible pente}
+\subsubsection{Le détecteur de bords}
-Le principe retenu pour réaliser le détecteur de LSR est proche de celui mis en oeuvre pour valider les allongements des isolines : il s'agit de séparer la fenêtre d'observation autour du pixel considéré en deux régions, 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$.
+Le principe retenu pour réaliser le détecteur de bords est proche de celui mis en oeuvre pour valider les allongements des isolines : il s'agit de séparer la fenêtre d'observation autour du pixel considéré en deux régions, puis d'effectuer un test GLRT pour déterminer s'il est vraisemblable ou non que ces deux régions forment un seul et même plan ou bien représentent deux zones homogènes séparées par un seuil (ou un coin). Pour garantir la prise en compte d'éventuelles transitions dans toutes les directions, il faut effectuer le test avec des séparations de fenêtre dont les directions couvrent toute la plage angulaire, de $0$ à $\pi$.
-L'utilisation d'un test GLRT semblable à celui de l'équation \eqref{GLRT} sous-entend que les ensembles considérés n'ont aucun pixel en commun. 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.
+L'utilisation d'un test GLRT semblable à celui de l'équation \eqref{GLRT} sous-entend que les ensembles considérés n'ont aucun pixel en commun. Afin d'éviter de devoir déterminer de nouveaux ensembles de pixels pertinents, nous avons utilisé les motifs de la matrice $P_d$, n'ayant pas d'intersection entre eux et de directions $\Theta_{4i} = 4i\frac{\pi}{4}$. Ces motifs remplissent les critères pour établir l'expression d'un test GLRT.
+La ligne de séparation entre les deux régions de la fenêtre est donc composée par les motifs de directions $\Theta_{4i}$ et $\Theta_{4(i+4)}$. Ces deux régions sont respectivement nommées arbitrairement $T$ et $B$, $T$ étant représentée comme la région \textit{haute} et $B$ comme la région \textit{basse} sur le schéma explicatif de la figure \ref{fig-lniv-detecteur} où $\Theta_{4i}=\frac{\pi}{4}$ et où les pixels affectés d'une élévation nulle sont les pixels non impliqués dans le calcul du test GLRT. En outre, les pixels de la limite sont supposés appartenir à la région $T$, ce qui implique qu'elle comprend au total les pixels correspondant à cinq motifs plus le pixel central, tandis que $B$ n'en comprend que l'équivalent de 3 motifs.
\begin{figure}[ht]
\centering
\label{fig-lniv-detecteur}
\end{figure}
-Les équations \eqref{LLNP}, \eqref{LLNP2} et \eqref{GLRT} nous permettent d'obtenir l'expression suivante pour le critère GLRT $T2$
+Les équations \eqref{LLNP}, \eqref{LLNP2} et \eqref{GLRT} nous permettent d'obtenir l'expression suivante pour le test GLRT $T2$
\begin{eqnarray}
T2 = T2_{max}- (8a+1)\left[log\left(\widehat{\sigma_3}^2\right) - log\left(\widehat{\sigma_4}^2\right) \right]
\label{GLRT2}
\end{array}
$$
-Le seuil de décision est noté $T2_{max}$ et d'après l'expression du critère \eqref{GLRT2}, une valeur négative du critère signifie la détection d'une transition. Ainsi, lorsque les valeurs du critère $T2$ sont connues pour toutes les 8 directions $\Theta_{4i} (i\in [0..7])$, la valeur du niveau de gris de sortie pour le pixel central est déterminée selon la stratégie suivante :
+Le seuil de détection est noté $T2_{max}$ et d'après l'expression du test \eqref{GLRT2}, une valeur négative signifie la détection d'un bord. Ainsi, lorsque les valeurs de $T2$ sont connues pour toutes les 8 directions $\Theta_{4i} (i\in [0..7])$, la valeur du niveau de gris de sortie pour le pixel central est déterminée selon la stratégie suivante :
\begin{itemize}
-\item si plus d'une valeur du critère est négative, alors on applique la valeur issue du modèle PI-PD.
-\item si une seule valeur du critère est négative, le pixel central est vraisemblablement situé sur une transition nette et on applique la valeur moyenne des motifs de la région $T$ à laquelle il appartient. Cela permet de garantir des transitions visuellement plus douces entre les zones où le PI-PD est appliqué et les zones moyennées.
-\item si aucune valeur du critère n'est négative, alors la région autour du pixel central est vraisemblablement une LSR. En conséquence, on applique la valeur moyenne de la zone.
+\item si plus d'un bord a été détecté, alors on applique la valeur issue du modèle PI-PD.
+\item si un seul bord a été détecté, le pixel central est vraisemblablement situé sur une transition nette et on applique la valeur moyenne des motifs de la région $T$ à laquelle il appartient. Cela permet de garantir des transitions visuellement plus douces entre les zones où le PI-PD est appliqué et les zones moyennées.
+\item si aucun bord n'a été détecté, alors la région autour du pixel central est vraisemblablement une LSR. En conséquence, on applique la valeur moyenne de la zone.
\end{itemize}
La figure \ref{fig-lniv-classification} présente le résultat de la classification des pixels d'une image bruitée, pour $T2_{max}=2$. On y remarque en particulier que les pixels noirs, pour lesquels s'appliquera le PI-PD, sont situés sur des transitions bien définies.
\centering
\subfigure[Image bruitée]{\includegraphics{Chapters/chapter4/img/airplane_noisy_small.jpg}}\qquad
\subfigure[Classification des pixels. ]{\includegraphics{Chapters/chapter4/img/img_bords_T2_small.jpg}}
-\caption{Classification des pixels d'une image bruitée, pour une valeur de seuil $T2=2$ du détecteur. (b) Les pixels en noir sont ceux à qui le PI-PD sera appliqué. Les pixels en blancs se verront appliquer une moyenne sur tout ou partie du voisinage.}
+\caption{Classification des pixels d'une image bruitée, pour une valeur de seuil $T2=2$ du détecteur. (b) Les pixels en noir sont ceux à qui le PI-PD sera appliqué. Les pixels gris se verront appliquer une moyenne sur tout ou partie du voisinage.}
\label{fig-lniv-classification}
\end{figure}
-
Les détails d'implémentation du détecteur sont donnés par l'algorithme \ref{algo-lniv-detecteur}. Pour en optimiser les performances, les sommes individuelles $sum_{\Theta}$ sont pré-calculées aux lignes 7 à 10 pour les 8 motifs concernés. L'évaluation des 8 configurations angulaires est effectuée ensuite de la ligne 11 à la ligne 25.
+Le choix que nous avons fait pour ce détecteur de bord sont stratégiques et ont été guidés par la recherche de performance. Il est cependant clair qu'un filtrage plus fort aurait été obtenu dans les LSR par un moyenneur intégrant l'ensemble des pixels de la fenêtre. Cette piste reste à approfondir à la lumière des résultats obtenu dans l'implémentation de l'opération de convolution.
+
\begin{algorithm}[ht]
\caption{Détecteur de zones à faible pente (LSR) \texttt{kernel\_LSR\_detector()}}
\label{algo-lniv-detecteur}
L'implémentation du PI-PD hybride a été appliquée aux 13 images de la base de test, dans leurs versions les plus bruitées, perturbées par un bruit gaussien de moyenne nulle et d'écart type 25.
Pour ce type d'images (taille, détails), les paramètres qui se sont avérés optimaux sont $a=5$ pour la longueur des segments avec un maximum de $q=5$ segments.
En ce qui concerne les seuils GLRT, nous avons testé l'ensemble des combinaisons de valeurs $T_{max}$ et $T2_{max}$ variant de 1 à 10 par pas de 0,5.
-La combinaison $T_{max}=1$ et $T2_{max}=2$ s'est révélée la plus appropriée, 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.
+La combinaison $T_{max}=1$ et $T2_{max}=2$ s'est révélée la plus appropriée, car elle représente le meilleur compromis PSNR/MSSIM pour 11 des 13 images, sauf \textit{peppers} et \textit{zelda}, pour lesquelles une combinaison $T_{max}=2$ et $T2_{max}=2$ permet d'améliorer l'indice de similarité MSSIM respectivement de 0,03 et 0,02.
Les images filtrées ont été caractérisées en termes de PSNR et de MSSIM et les résultats, regroupés dans la table \ref{tab-lniv-results}, sont comparés à ceux de la référence BM3D, ainsi qu'à ceux d'un simple filtre moyenneur GPU 5$\times$5, choisi comme référence en terme de rapidité et dont la taille de fenêtre permet des gains théoriques en PSNR du même ordre de grandeur que le PI-PD.
\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
+À notre connaissance, la définition des lignes de niveaux n'a pas été formellement étendue aux espaces colorimétriques à plusieurs canaux.
+Nous faisons ici l'hypothèse (forte) de canaux décorrelés, conduisant à une matrice de covariance diagonale de la forme $R=\sigma^2\mathbb{1}_3$ où $\sigma^2$ est la puissance du bruit gaussien IID perturbant les trois canaux.
+La vraisemblance de $v_k$ est alors
$$P\left(v_k|R\right) = \left(\frac{1}{2\pi^{3/2}\sqrt{|R|}}\mathrm{e}^{-\frac{1}{2}\left(v_k-\mu\right)^TR^{-1}\left(v_k-\mu\right)}\right)$$
-Pour exprimer le critère GLRT de validation des allongements, nous procédons comme précédemment, c'est-à-dire en distinguant les deux hypothèses :
+Pour exprimer le test GLRT de validation des allongements, nous procédons comme précédemment, c'est-à-dire en distinguant les deux hypothèses :
\begin{enumerate}
\item le segment candidat $S^p$ prolonge effectivement l'isoline $S^n$ : ils partagent donc la même valeur moyenne $\mu$ et la log-vraisemblance s'écrit
\begin{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) $$
+Le test GLRT s'obtient par la soustraction des deux expressions de \eqref{eqlv1rgb} et \eqref{eqlv0rgb} :
+\begin{eqnarray}
+T_{rvb} = 3(n+p)\left(-log\left(\widehat{{\sigma_1}^{2}}\right)+log\left(\widehat{{\sigma_0}^{2}}\right)\right)
+\end{eqnarray}
+
On notera $T_{rvb-max}$ la valeur de seuil au delà de laquelle on ne validera pas l'allongement de l'isoline.
-\subsection{Résultats}
+\subsection{Résultats - analyse}
Nous avons retenu la base d'images de test tid2008 \cite{tid2008a} pour évaluer la qualité du traitement PI-PD sur les images couleurs.
Cet ensemble d'images a été utilisé avec nombre d'algorithmes de débruitage et les résultats de mesure sont disponibles.
-Chacune des 25 images de référence (non bruitées) a 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.
+Chacune des 25 images de référence (non bruitées) a subi 4 niveaux de distorsion, pour 17 types de bruit différents. Pour nos expérimentations, nous avons selectionné les 25 images corrompues par un bruit gaussien RVB (type 2 dans tid2008) d'écart type $\sigma = 25$ (niveau 4 dans tid2008), où chaque canal RVB est perturbé par un tirage de bruit gaussien IID. La figure \ref{fig-lniv-tid2008ref} présente les vignettes des 25 images de référence, soit 24 images \textit{naturelles} et une image de synthèse.
\begin{figure}[ht]
\centering
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 $n$ des segments, de 25 à 70 pour la longueur maximale $l$ des isolines et de 1 à 10 pour le seuil GLRT $T_{rvb-max}$. Cette analyse extensive a mis en évidence la combinaison $n=4$, $l=48$ et $T_{rvb-max}=5$ comme permettant au PI-PD d'apporter les meilleurs résultats d'ensemble. Certaines des images, comme l'image de synthèse n°25, bénéficieraient d'un ajustement des paramètres, mais conscients de la contrainte que cela représente, nous avons choisi de faire prévaloir un réglage unique.
+L'ensemble des résultats de mesure est consigné dans le tableau \ref{tab-lniv-rvb} et deux exemples de résultats sont reproduits en figure \ref{fig-lnivrgb-ex} pour une des images naturelles ainsi que pour l'image de synthèse. Les valeurs des paramètres sont identiques pour toutes les images et ont été déterminées empiriquement par analyse systématique des résultats produits par les combinaisons permises dans les intervalles de 3 à 7 pour la taille $n$ des segments, de 25 à 70 pour la longueur maximale $l$ des isolines et de 1 à 10 pour le seuil GLRT $T_{rvb-max}$. Cette analyse extensive a mis en évidence la combinaison $n=4$, $l=48$ et $T_{rvb-max}=5$ comme permettant au PI-PD d'apporter les meilleurs résultats d'ensemble. Certaines des images, comme l'image de synthèse n°25, bénéficieraient d'un ajustement des paramètres, mais conscients de la contrainte que cela représente, nous avons choisi de faire prévaloir un réglage unique.
+
+Si le modèle que nous avons employé pour étendre le principe du PI-PD aux images en couleur nous a permis de vérifier que les performances pouvaient être conservées, il repose sur des hypothèses fortes et sous-entend (pour simplifier) qu'un seuil dans l'espace RVB est associé à un seuil sur chaque composante. Pour améliorer la qualité, une autre solution consisterait à détecter un seuil dans l'espace RVB lorsqu'un seuil est détecté dans au moins un des canaux. La transposition vers un autre espace colorimétrique employant une base de composantes décorrélées pourrait également améliorer les résultats. Ces pistes n'ont pas encore été explorées, mais représenteront nécessairement un surcoût calculatoire.
+
\label{fig-lnivrgb-ex}
\begin{table}[H]
\scriptsize
\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.
+La démarche adoptée pour sa conception a été de se baser sur des opérations élémentaires dont nous connaissions ou avions démontré l'efficacité sur GPU. Nous jugeons ce principe essentiel pour la conception d'algorithmes GPU performants et robustes tant le débogage peut s'avérer délicat sur ces plateformes. Par ailleurs, il nous semble peu pertinent de systématiquement comparer les implémentations CPU et GPU pour en déduire un facteur d'accélération comme on le rencontre trop souvent. La plupart des algorithmes qui s'avèrent rapides sur GPU ne le sont vraisemblablement pas sur CPU et il est donc tout à fait illusoire de penser qu'il en existe une implémentation optimisée. Comparer alors une implémentation GPU performante avec son pendant CPU naïf ne présente aucun intérêt. La réciproque étant généralement vraie, nous avons choisi, en particulier en ce qui concerne le filtrage dont il est question ici, de chercher à assembler des blocs fonctionnels simples mais robustes et performants avec l'objectif opérationnel de réduire la puissance de bruit.
L'algorithme et les résultats que nous avons détaillés dans ce chapitre ont été publiés dans le \textit{Journal of real-time image processing} dans un article intitulé \textit{Fast GPU-based denoising filter using isoline levels} \cite{perrotlniv}.
-% LocalWords: pénalisante
+% LocalWords: pénalisante estimateurs colorimétriques
