8.2 Information a priori et Régularisation.

Nous allons considérer les images comme un ensemble de régions homogènes séparées par des bords francs. C'est pourquoi, nous envisageons la régularisation dans une approche markovienne. En effet, les champs markoviens sont bien adaptés à cette hypothèse. Nous rappelons brièvement certaines notions relatives aux champs markoviens dont nous aurons besoin.. Cette description, bien connue actuellement, s'appuie sur les travaux fondateurs de S. et D. Geman [38].

8.2.1 Voisinage.

Dans l'approche statistique, l'image $\mathbf{f}$ correspond à un champ aléatoire, de sorte qu'on peut associer la probabilité $I\! \! P(\mathbf{f})$ d'obtenir cette image. Chaque élément $f_{n_{s}}$, ou voxel, est localisé en un site $s$ (sa valeur est prise dans l'espace des états $\Lambda _{f}$). L'ensemble des sites $S$ représente une grille qui correspond, dans notre cas, au centre de chaque voxel. Le voisinage $\mathcal{V}$ associe à tout site $s$ de la grille $S$, une série de sites $t$ de cette même grille (ses voisins). Le système de voisinage $\mathcal{V}=\{\mathcal{V}_{s}\subset S,s\in S\}$ doit vérifier les deux conditions suivantes:

$$\triangleright \forall s\in S,\; s\notin \mathcal{V}_{s}$$

$$\triangleright \forall s,t\in S,\; s\in \mathcal{V}_{t}\Leftrightarrow t\in \mathcal{V}_{s}$$

Le voxel n'est pas voisin avec lui même. Et si le voxel $f_{n_{s}}$ du site $s$ est voisin du voxel $f_{n_{t}}$ au site $t$, alors $f_{n_{t}}$ est aussi voisin de $f_{n_{s}}$(réciprocité du voisinage). D'autre part, pour un site $s$ , il peut être utile de vouloir énoncer tous ses voisins. Cela se fait par l'utilisation d'une clique. Une clique $c$, c'est d'abord le site $s$ considéré puis le sous ensemble des sites de $S$ qui lui sont voisins 2 à 2. Toutes les cliques que l'on peut construire sur la grille $S$ forment l'ensemble des cliques $\mathcal{C}$:

$$c\in \mathcal{C}\Leftrightarrow \left\{ \begin{array}{l} s\\ ou\\ \forall \{s,t\}\subset c,\; s\in \mathcal{V}_{s} \end{array}\right.$$

On dit que $\mathbf{f}$ est un champ de Markov par rapport au système de voisinage $\mathcal{V}$ si l'image vérifie, pour tout $s\in S$:

• Une propriété de positivité:
  $$I\! \! P(f_{n_{s}})>0$$

• Une propriété de localité:

  $$I\! \! P(f_{n_{s}}\vert f_{n_{t}},s\neq t)=I\! \! P(f_{n_{s}}\vert f_{n_{t}},t\in \mathcal{V}_{s})$$ (8.4)

  
  

La deuxième propriété revient à dire que la connaissance du voisinage $\mathcal{V}_{s}$ d'un site $s$ est suffisante pour calculer la probabilité marginale en ce site.

Derrière ces définitions mathématiques se retrouvent des notions facilement palpables. En effet, il ne faut pas oublier que notre vecteur $\mathbf{f}$ représente une distribution de radioactivité. Il y a de grandes chances pour que cette distribution soit relativement lisse. Deux zones proches dans l'espace auront une activité ``voisine''. De surcroît, un forte augmentation locale se traduira également par une augmentation dans une région périphérique. Autrement dit, les voxels sont spatialement corrélés. Cette notion de proximité, de corrélation, se traduit aisément à l'aide d'un voisinage basé sur la distance entre deux voxels. Deux voxels sont voisins si la distance qui les sépare est inférieure à une valeur seuil. Cette distance fixe l'ordre du voisinage. La notion de champ de Markov nous dit que les répercussions d'un site sont confinées à son voisinage. La probabilité en un site, n'est pas fonction de tous les autres voxels mais seulement des sites voisins. Dans le cas où on considère le voisinage associé à la distance entre voxels, le fait qu'un voxel ait une valeur $f_{n_{s}}$ ne dépend en termes statistiques que de l'état des voxels $f_{n_{t}}$ qui sont proches du voxel $n_{s}$ ( $t\in \mathcal{V}_{s}$). Le voisinage basé sur la distance est certainement le plus courant et conduit à un voisinage dont la forme est la même quel que soit le site ou voxel considéré (aux bords du volume près). On en donne une représentation Fig.9.2. Les ronds correspondent à des voisins du site $s$. En fonction de la distance à ce site, une couleur indique l'ordre du voisinage considéré. Un voisinage basé sur 4-connexité correspond à un voisinage d'ordre 1, la 8-connexité à un voisinage d'ordre 2.

Figure 9.2: Représentation du système de voisinage basé sur la distance. Chaque disque plein représente un voisin du site $s$ central. Les disques noirs correspondent à un voisinage d'ordre 1 (4-connexité), les bleus à un voisinage d'ordre 2 (8-connexité), les rouges à un voisinage d'ordre 3, les verts à l'ordre 4. L'ordre du voisinage est basé sur la distance de voisin au site.

8.2.1.0.0.1 Remarque:

Notons, qu'un tel voisinage exclut la corrélation entre des zones éloignées dans l'objet. Par exemple, si nous savons a priori que deux régions cérébrales distantes sont corrélées de manière fonctionnelle, le voisinage construit sur la distance n'en tiendra pas compte^8.1 . Pour intégrer cette information au niveau de notre reconstruction (ce qui est possible), il faudrait définir un voisinage utilisant des informations fonctionnelles qui ne sont pas uniquement basées sur la proximité. Il faudrait définir un voisinage dédié.

La propriété de localité Eq.9.4 nous donne des indications sur les propriétés statistiques de manière ponctuelle. Or, nous cherchons à caractériser la probabilité globale $I\! \! P(\mathbf{f})$. Le lien entre ces deux grandeurs nous est fourni par le théorème de Hammersley-Clifford.

8.2.2 Potentiel d'interaction.

Ce théorème nous dit :

Théorème I : Hammersley-Clifford

Si $\mathbf{f}$ est markovien sur la grille $S$ avec le système de voisinage choisi $\mathcal{V}$, alors la probabilité d'observer $\mathbf{f}$ suit une distribution de Gibbs relativement à $\mathcal{V}$. Cette distribution est de la forme:

$$I\! \! P(\mathbf{f})=\frac{1}{Z}exp\left[ -\sum _{c\in \mathcal{C}}U_{c}(\mathbf{f})\right]$$

où $Z$ représente une fonction de partition $Z=\sum _{\mathbf{f}\in \Lambda ^{N}_{f}}exp\left[ -\sum _{c\in \mathcal{C}}U_{c}(\mathbf{f})\right]$ qui nous sert à normaliser la distribution pour qu'elle définisse bien une probabilité. Son intégration sur toutes les configurations possibles est 1.

Cette expression correspond à une distribution connue en physique statistique. Chaque valeur $f_{n}$ correspond à un état énergétique. Chaque image correspond donc à une configuration énergétique précise où chaque voxel est dans un état énergétique particulier. La somme $U=\sum U_{c}$ traduit l'énergie associée à cette configuration. La localité nous dit que cette fonction d'énergie globale est en fait la somme des énergies locales pour chacun des sites (le potentiel d'interaction $U_{c}$ est défini pour une clique qui correspond à l'énumération des voisins en un site). Nous voyons également que l'image la plus probable est celle qui rendra $U$ minimale. Toujours par analogie avec la physique, cela veut dire que la configuration la plus stable est celle d'énergie minimale.

Le potentiel d'interaction $U_{c}$ est donc une fonction construite sur une clique $c$. Plaçons nous en un site $s$, la clique $c$ construite sur ce site nous donne des voisins (des sites) qui interagissent avec $s$. Si nous prenons deux images différentes mais dont les valeurs en chacun des sites définis par cette clique $c$ sont égales, alors les potentiels d'interactions calculés sur chacune de ces deux images seront égaux:

$$_{1}\{f_{n_{s}}\}_{s\in c}=_{2}\{f_{n_{s}}\}_{s\in c}\Rightarrow U_{c}(_{1}\mathbf{f})=U_{c}(_{1}\mathbf{f})$$

On dispose ainsi d'un nouveau degré de liberté. Une fois que le système de voisinage est fixé, il est encore nécessaire de choisir le mode d'interactions entre les voisins. Le choix du potentiel d'interaction est fonction du but recherché. Fréquemment, en segmentation, on travaille avec un modèle d'Ising ou un modèle de Potts à p niveaux. Dans notre cas, nous avons dans l'esprit de régulariser (de lisser), on va donc introduire une fonction de potentiel construite sur un opérateur différentiel. Si on prend deux sites $s$ et $t$ d'une clique $c$, le potentiel d'interaction sera lié à l'opérateur différentiel d'ordre $m$ : $(D^{m}f)_{n_{s},n_{t}}$. En fonction de l'ordre, cet opérateur agit soit sur l'image elle même (ordre 0), soit sur l'approximation discrète de la dérivée première (ordre 1) ou seconde (ordre 2). L'énergie est alors [39]:

$$U_{c}(\mathbf{f})=\sum _{s,t\in c}w_{st}\varphi (\{D^{m}f\}_{n_{s},n_{t}})$$

où $w_{st}$ est un facteur de pondération associé à un couple $s,t$ particulier et $\varphi$ est une fonction de potentiel unique quel que soit le site visité et la clique considérée.

Tout ceci est bien général. Notre but n'est pas de faire une description aussi exhaustive que possible, mais juste de sensibiliser en montrant la richesse et la souplesse du modèle. Nous allons donc maintenant restreindre et illustrer le modèle retenu.

8.2.3 Cas Pratique.

On ne considère que l'opérateur différentiel d'ordre 1 ($D^{1}f=Df$) et le système de voisinage $\mathcal{V}$ basé sur la distance. Dans une image $\mathbf{f}$, envisageons un voxel (site $s)$ particulier $f_{n_{s}}=f_{i,j,k}$. Il nous faut, dans un premier temps, définir ses voisins. Pour l'illustration, nous choisissons une 6-connexité. Nous avons, en plus du site $s$, 6 voisins qui sont $\{n_{a},n_{p},n_{g},n_{d},n_{s},n_{i}\}$. A ces sites voisins correspondent respectivement les valeurs dans le volume $f_{i+1,j,k}$, $f_{i-1,j,k}$, $f_{i,j+1,k}$, $f_{i,j-1,k}$, $f_{i,j,k+1}$, $f_{i,j,k-1}$.

Dans un deuxième temps, il faut construire pour ce site $s$, le potentiel d'interaction associé à chaque clique $c$. De manière discrète, l'opérateur différentiel d'ordre 1 se traduit juste par une différence entre les valeurs de deux voisins ( $\{Df\}_{n_{s},n_{t}}=f_{n_{s}}-f_{n_{t}}$). Nous obtenons donc directement le potentiel d'interaction pour la clique construite au voxel $f_{n_{s}}$:

$$U_{c}(\mathbf{f}) = \varphi (f_{i+1,j,k}-f_{i,j,k})+\varphi (f_{i,j,k}-f_{i-1,j,k})$$

$$+ \varphi (f_{i,j+1,k}-f_{i,j,k})+\varphi (f_{i,j,k}-f_{i,j-1,k})$$

$$+ \varphi (f_{i,j,k+1}-f_{i,j,k})+\varphi (f_{i,j,k}-f_{i,j,k-1})$$

Evidemment pour avoir l'énergie totale $U$, il est nécessaire de collecter toutes les contributions locales (ensemble des cliques $C$) pour tous les voxels (sites):

$$U(\mathbf{f})=\sum _{c\in \mathcal{C}}U_{c}(\mathbf{f}) = \sum _{i,j,k}\varphi (f_{i+1,j,k}-f_{i,j,k})+\sum _{i,j,k}\varphi (f_{i,j,k}-f_{i-1,j,k})$$

$$+ \sum _{i,j,k}\varphi (f_{i,j+1,k}-f_{i,j,k})+\sum _{i,j,k}\varphi (f_{i,j,k}-f_{i,j-1,k})$$

$$+ \sum _{i,j,k}\varphi (f_{i,j,k+1}-f_{i,j,k})+\sum _{i,j,k}\varphi (f_{i,j,k}-f_{i,j,k-1})$$

Il est clair que $-\log P(\mathbf{f})=U(\mathbf{f})$, l'énergie totale correspond donc à notre fonctionnelle $J_{2}(\mathbf{f})$. Dans cette optique, la régularisation consiste à chercher l'image présentant le minimum d'énergie. Insistons encore sur le fait que l'énergie envisagée repose déjà sur deux hypothèses :

• Une qui concerne le système de voisinage (basé sur la distance entre voxel)
• Une autre sur le type d'interactions envisagées (basé sur la gradient discret de l'image).

Il reste encore à définir une fonction de potentiel $\varphi$. Le choix de cette fonction est lié à ces hypothèses :

• Nous voulons lisser notre image dans les zones où elle est homogène tout en préservant les discontinuités.
• Nous voulons utiliser un algorithme déterministe pour minimiser notre fonctionnelle.

8.2.4 Régularisation.

8.2.4.1 Hypothèses.

Pour expliciter la régularisation et donc la prise en compte des discontinuités, nous prenons l'approche de Charbonnier [17]. Dans cette approche, notre but est de trouver les propriétés que doit suivre la fonction de potentiel $\varphi$ pour préserver les discontinuités dans notre image en partant du fait que nous devons rendre le critère $J_{2}$ minimal. Nous replacerons ensuite les conditions obtenues dans le contexte des processus de lignes.

Avant d'aller plus loin, explicitons de nouvelles hypothèses:

• Par la régularisation, nous souhaitons effectuer un lissage qui vise à éliminer les effets du bruit sur la solution estimée. Ce bruit se traduit dans une zone homogène par de petites variations indésirables autour de la moyenne. Notre but est donc le lissage de ces irrégularités. C'est pourquoi, dans notre fonctionnelle, nous travaillons sur le gradient de l'image (opérateur différentiel de premier ordre) qui estime ces irrégularités. En revanche, nous voulons préserver les contours des objets présents dans le volume qui sont des éléments importants. On fait donc l'hypothèse que les forts gradients dans l'image correspondent à des contours et non à du bruit. Le contour se traduit donc par une discontinuité dans le volume. Ce n'est évidemment pas toujours vrai en pratique. De forts gradients peuvent être liés à du bruit, et inversement, des gradients faibles peuvent également correspondre à des contours inclus dans l'image. Surtout en présence d'effet de volume partiel dans le volume!
• On suppose également que l'on peut traiter de manière identique les composantes positives et négatives de notre gradient. Les fonctions de potentiel utilisées sont donc paires et travaillent sur la valeur absolue du gradient.
• Enfin, on suppose que la fonction de pondération que nous utilisons peut être la même pour les composantes en $x$, $y$, et $z$. C'est pourquoi, quand nous parlerons de gradient nous désignerons indifféremment l'image de gradient suivant $x$: $\{D_{x}f\}_{i,j,k}=f_{i,j,k}-f_{i-1,j,k}$ l'image de gradient suivant $y$: $\{D_{y}f\}_{i,j,k}=f_{i,j,k}-f_{i,j-1,k}$ et l'image de gradient suivant $z$: $\{D_{z}f\}_{i,j,k}=f_{i,j,k}-f_{i,j,k-1}$.

8.2.4.2 Lissage, préservation des discontinuités et Laplacien pondéré.

Nous avons vu que pour minimiser $J_{2}$, nous devions minimiser la somme des énergies $U_{c}(\mathbf{f})$ en chaque site $n_{s}=\{i,j,k\}$. Plaçons nous en un site particulier et minimisons localement. Supposons que nous sommes en 2D, nous avons donc une image $\mathbf{f}=\{f\}_{n_{s}}=\{f\}_{i,j}$. Si on suppose un voisinage basé sur une 4-connexité, nous avons pour chaque site, 4 voisins $\mathcal{V}_{s}=\{n_{a},n_{p},n_{g},n_{d}\}$ qui sont impliqués pour le système de voisinage considéré (Fig.9.2). La minimisation locale du critère s'écrit donc:

$$\left( \frac{\delta J_{2}(\mathbf{f})}{\delta f_{n_{s}}}\right) = \sum _{n_{t}\in \mathcal{V}_{s}}\varphi '(f_{n_{t}}-f_{n_{s}})$$

$$= \sum _{n_{t}\in \mathcal{V}_{s}}2(f_{n_{t}}-f_{n_{s}})\times \und... ...rac{\varphi '(f_{n_{t}}-f_{n_{s}})}{2(f_{n_{t}}-f_{n_{s}})}}_{\lambda _{n_{t}}}$$

$$= \sum _{n_{t}\in \mathcal{V}_{s}}2(f_{n_{t}}-f_{n_{s}})\times \lambda _{n_{t}}$$

$$= 2\times \left( f_{n_{s}}\sum _{n_{t}\in \mathcal{V}_{s}}\lambda _{n_{t}}-\sum _{n_{t}\in \mathcal{V}_{s}}\lambda _{n_{t}}f_{n_{t}}\right)$$

$$= -2\Delta _{pond}(f_{n_{s}}).\mathbf{f}$$ (8.5)

La minimisation du critère localement revient (Eq.9.5.) à filtrer l'image $\mathbf{f}$ par un Laplacien pondéré $-\Delta _{pond}(f_{n_{s}})$ négatif. Les coefficients de ce Laplacien pondéré nous sont fournis par les valeurs $\lambda _{n_{t}}$ calculées pour les voisins considérés.^8.2 Restons en 2D pour simplifier et explicitons le masque bidimensionnel correspondant à ce Laplacien pour le site $s$ (4-connexité):

$$\Delta _{pond}(f_{n_{s}})\leftrightarrow \left[ \begin{array}{ccc... ...a _{n_{t}}} & \lambda _{n_{d}}\\ 0 & \lambda _{n_{a}} & 0 \end{array}\right]$$

Supposons maintenant que tous les coefficients sont égaux à 1 ( $\frac{\varphi '(x)}{2x}=1$), alors ce Laplacien pondéré devient:

$$\Delta _{pond}(f_{n_{s}})\leftrightarrow \left[ \begin{array}{ccc} 0 & 1 & 0\\ 1 & -4 & 1\\ 0 & 1 & 0 \end{array}\right]$$

qui correspond à l'expression connue du Laplacien discret permettant de lisser de manière homogène une image bidimensionnelle. Un lissage homogène est effectué suivant les 4 directions considérées par le voisinage. Supposons maintenant que l'un de ces coefficients ( $\lambda _{n_{g}}$ par exemple) soit nul ( $\frac{\varphi '(x)}{2x}=0$). Alors notre opérateur se résume à:

$$\Delta _{pond}(f_{n_{s}})\leftrightarrow \left[ \begin{array}{ccc} 0 & 1 & 0\\ 0 & -3 & 1\\ 0 & 1 & 0 \end{array}\right]$$

Nous avons dans cet opérateur de lissage cassé l'interaction, et donc le lissage, qui existait entre le site $s$ et son voisin situé à gauche $n_{g}$.

En résumé, on constate que pour lisser dans des zones homogènes, il nous faut des coefficients proches de 1. Pour préserver une discontinuité, et donc casser l'opération de lissage, il faut que le coefficient soit nul pour cette direction .

8.2.4.3 Conditions pour un lissage préservant les discontinuités.

Cette approche, nous permet aisément de comprendre les conditions imposées, non sur la fonction de potentiel $\varphi (x)$, mais sur la fonction de pondération $\frac{\varphi '(x)}{2x}$. Nous venons de voir les conditions aux limites. en rajoutant comme contrainte pour la fonction de pondération le fait qu'elle doit varier de manière continue entre ses limites, nous aboutissons aux conditions suivantes:

Pour lisser tout en préservant les discontinuités, la fonction de potentiel $\varphi$ mais surtout la fonction de pondération utilisée doit vérifier les conditions suivantes [17]:

$$i) \lim _{u\rightarrow 0^{+}}\frac{\varphi '(u)}{2u}=1$$

$$ii) \lim _{u\rightarrow \infty }\frac{\varphi '(u)}{2u}=0$$ (8.6)

$$iii) \frac{\varphi '(u)}{2u}\; est\; continue\; et\; strictement\; d\acute{e}croissante\; sur\; [0,\infty [$$

où $u$ représente la valeur du gradient dans une direction particulière. Nous avons considéré qu'une discontinuité correspondait à une forte valeur du gradient. Cette notion est relative et dépend bien évidemment des valeurs des voxels dans l'image. Nous sommes donc amené à considérer une constante de normalisation, ou hauteur des discontinuités $\delta$, pour ramener les valeurs du gradient dans une gamme de valeurs analogues quel que soit le niveau du signal dans l'image. $u=\frac{\{Df\}_{n_{s},n_{t}}}{\delta }$ représente donc la valeur du gradient normalisée par cette hauteur des discontinuités. Par la suite $u$ représentera une valeur de gradient suivant une direction particulière normalisée par la hauteur des discontinuités.

8.2.4.4 Choix d'une fonction de potentiel.

Si on consulte la littérature pour retrouver les fonctions de potentiels qui furent utilisées pour lisser des volumes, on constate que les scientifiques ne manquent pas d'imagination. Les approches diffèrent et conduisent parfois à des conditions contradictoires sur les fonctions de potentiel ([54] versus [39]). Bien souvent, la recherche est empirique [53,42]. Dans le tableau Tab.9.1,nous donnons quelques unes des fonctions les plus utilisées dans la littérature en précisant celles susceptibles de préserver les discontinuités en accord avec le critère que nous avons retenu.

Table: Résumé des fonctions de potentiel.

Nom de la fonction.	$\varphi (u)$	$\frac{\varphi '(u)}{2u}$	Eq.9.6 vérifiées	Convexité
Tikhonov[3]	$u^{2}$	1	non	oui
Variation totale[1]	$u$	$\frac{1}{2u}$	non	oui
Quadratique tronquée[12]	$\min (u^{2},1)$	$\left\{ \begin{array}{c} 1\; u<1\\ 0\; u\geq 1 \end{array}\right.$	non	non
Huber[82]	$\min (u^{2},2u-1)$	$\left\{ \begin{array}{c} 1\; u<1\\ \frac{1}{u}\; u\geq 1 \end{array}\right.$	non	non
Bouman & Sauer[14]	$\vert u\vert^{\alpha }$ $1\leq \alpha <2$	$\frac{\alpha }{2}\vert u\vert^{\alpha -2}$	non	oui
Perona & Malik[77]	$-e^{-u^{2}}-1$	$e^{-u^{2}}$	oui	non
Geman & Mc Glure[37]	$\frac{u^{2}}{1+u^{2}}$	$\frac{1}{(1+u^{2})^{2}}$	oui	non
Hebert & Leahy[42]	$\log (1+u^{2})$	$\frac{1}{1+u^{2}}$	oui	non
Green[41]	$2\log [\cosh (u)]$	$\frac{\tanh (u)}{u}$	oui	oui
Hyper surfaces[17]	$2\sqrt{1+u^{2}}-1$	$\frac{1}{2\sqrt{1+u^{2}}}$	oui	oui

La liste n'est pas exhaustive mais traduit bien la difficulté du choix. Ce choix est lié aux hypothèses que nous faisons, et donc à l'information a priori considérée. Il va sans dire que toutes ces fonctions ne seront pas forcément adaptées à la reconstruction d'image TEP. La régularisation de Tikhonov par exemple correspond à un lissage homogène par un Laplacien classique et ne préserve pas les discontinuités. Dans notre cas, notre démarche est inverse, nous ne cherchons pas à déduire des lois partant de fonctions optimales, nous recherchons parmi toutes les fonctions mentionnées, celles qui vérifient les conditions Eq.9.6. On élimine donc la fonction de Tikhonov [3], celle de variation totale [1], celle de Bouman&Sauer [14], la Quadratique tronquée [12], la fonction de Huber [82]). Sur la figure Fig.9.3

Figure 9.3: Fonctions de pondération vérifiant Eq.9.6. Sont représentées les fonctions de pondération de Perona & Malik (PM) [77] de Geman & McGlure (GM) [37] de Hebert & Leahy (HL) [42] de Green (G) [41] des Hyper-Surfaces (HS) [17].

sont représentées les fonctions de pondération ( $\frac{\varphi '(u)}{2u}$) vérifiant les conditions de préservation des discontinuités. Il faut noter que si les fonctions de pondération se ressemblent, elles ne décroissent pas à la même vitesse. La fonction de pondération associée à Perona & Malik (PM) [77] décroît en $\frac{1}{e^{u^{2}}}$, celle de Geman & McGlure (GM) [37] en $\frac{1}{u^{4}}$, Celle de Hebert & Leahy (HL) [42] en $\frac{1}{u^{2}}$alors que celles de Green (G) [41] ou des Hyper-Surfaces (HS) [17] suivent une loi en $\frac{1}{u}$. Or la fonction de pondération établit un passage continu du lissage ( $\frac{\varphi '(u)}{2u}=1$) à l'absence de lissage ( $\frac{\varphi '(u)}{2u}=0$). La transition est d'autant plus rapide (contour plus net) que la fonction décroît rapidement. En prenant celle dont la décroissance est la plus rapide, nous avons des contours bien net! Malheureusement, les fonctions présentant une décroissance rapide ne sont pas convexes. Cela veut dire que le critère de MAP obtenu en utilisant ces fonctions non convexes présente des minima locaux. La présence de minima locaux est particulièrement gênante car elle compromet fortement l'utilisation d'algorithme de minimisation déterministes. L'algorithme restant coincé dans l'un de ces minima locaux. Bouman et Sauer [14] n'hésitent pas à dire que le problème reste mal posé , puisqu'il n'y a pas unicité de la solution^8.3. Il faut donc choisir un compromis entre la netteté de l'image (fonction qui décroît rapidement) et la présence d'un minimum global (fonction qui décroît plus lentement). On voit que l'utilisation d'un algorithme déterministe contraint l'utilisation d'une fonction de pondération soit de Green, soit des hypersurfaces. Sur la figure Fig9.3, on constate que ces deux fonctions sont très proches, on a donc tout intérêt à prendre la fonction des hypersurfaces vu que son utilisation est plus aisée.

8.2.4.5 Régularisation semi-quadratique.

Après avoir choisi le voisinage, travaillant avec l'opérateur différentiel du premier ordre avec la fonction de potentiel retenue, on cherche à minimiser le critère $J_{2}$:

$$J_{2}(\mathbf{f})=\sum _{n_{s}}\sum _{n_{t}\in \mathcal{V}_{s}}\varphi (\underbrace{f_{n_{t}}-f_{n_{s}}}_{\{Df\}_{n_{s},n_{t}}})$$

Nous aboutissons à:

$$\left( \frac{\delta J_{2}(\mathbf{f})}{\delta f}\right) =-\Delta _{pond}(\mathbf{f}).\mathbf{f}$$

Malheureusement, le Laplacien pondéré dépend de l'image! La minimisation ne conduit pas à une expression linéaire en $\mathbf{f}$. D'où l'idée de régularisation semi-quadratique introduite par Geman et Yang en 1993 [36].

On construit une nouvelle fonction d'énergie $J^{*}_{2}(\mathbf{f},\mathbf{d})$ qui, bien que définie sur un domaine plus étendu ( $\mathbf{f}$ et $\mathbf{d}$), a le même minimum global que $J_{2}(\mathbf{f})$ en $\mathbf{f}$. Elle peut être manipulée par des méthodes d'algèbre linéaire.

On introduit donc une variable auxiliaire $\mathbf{d}$ qui va rendre le critère $J_{2}$ plus facile à manipuler. On remplace le critère:

$$J_{2}(\mathbf{f})=\sum _{n_{s}}\sum _{n_{t}\in \mathcal{V}_{s}}\varphi (\{Df\}_{n_{s},n_{t}})$$

par:

$$J_{2}^{*}(\mathbf{f},\mathbf{d})=\sum _{n_{s}}\sum _{n_{t}\in \mathcal{V}_{s}}\varphi ^{*}(\{Df\}_{n_{s},n_{t}},d_{n_{s},n_{t}})$$

Pour l'instant, les choses sont devenues plus compliquées qu'elles ne l'étaient auparavant. Toutefois, comme l'on fait avant nous Geman en 1992 [39], Aubert en 1994 [6] et Charbonnier en 1994 [17], envisageons le théorème de Geman & Reynolds étendu. Ce théorème s'appuie sur les travaux de Geman & Reynolds et Geman & Yang. Charbonnier en propose une extension qui inclut les fonctions de potentiels convexes.

Théorème II : Geman & Reynolds étendu

Soit une fonction $\varphi$ définie de $[0,+\infty [\rightarrow [0,+\infty [$ qui vérifie:

$$i) \lim _{u\rightarrow 0^{+}}\frac{\varphi '(u)}{2u}=M$$

$$ii) \lim _{u\rightarrow \infty }\frac{\varphi '(u)}{2u}=L$$

$$iii) \varphi (\sqrt{u})\; strictement\; concave\; sur\; [0,+\infty [$$

Alors:

1. Il existe une fonction $\psi :[L,M]\rightarrow [\alpha ,\beta ]$ strictement convexe et décroissante telle que:
   $$\varphi (u)=\inf _{L\leq d\leq M}\left( du^{2}+\psi (d)\right)$$
   où
   

   $$\alpha = \lim _{u\rightarrow +\infty }\left( \varphi (u)-u^{2}\frac{\varphi '(u)}{2u}\right)$$

   $$\beta = \lim _{u\rightarrow 0^{+}}\varphi (u)$$

   
   

2. Pour tout $u\geq 0$ fixé, la valeur $d_{opt}$ pour laquelle le minimum est atteint:
   $$\varphi (u)=\left( d_{opt}u^{2}+\psi (d_{opt})\right)$$
   est unique et donnée par:
   $$d_{opt}=\frac{\varphi '(u)}{2u}$$

De par ce théorème, si nous posons $\varphi ^{*}(u,d)=\left( du^{2}+\psi (d)\right)$, nous arrivons à montrer que le critère étendu $J^{*}_{2}$ et le critère classique $J_{2}$ sont égaux pour une valeur de $\mathbf{d}$ précise.

$$J_{2}(\mathbf{f})=J_{2}(\mathbf{f},\mathbf{d})\vert _{\mathbf{d}_{opt}}$$

Arrêtons nous sur les conditions de ce théorème. Si $L=0$ et $M=1$ alors ces conditions sont équivalentes à celles posées sur la fonction de potentiel Eq.9.6. En effet, il y a équivalence entre la concavité d'une fonction et la stricte décroissance de sa dérivée. Or, la dérivée de $\varphi (\sqrt{u})$ est justement la fonction de pondération.

En choisissant une fonction de potentiel vérifiant les conditions Eq.9.6, nous pouvons construire une fonction étendue $\varphi ^{*}(u,d)$. Cette fonction nous permet de construire un nouveau critère $J^{*}_{2}$ étendu présentant le même minimum en $\mathbf{f}$ que $J_{2}$. Qui plus est :

• Lorsque $\mathbf{d}$ est fixé, le critère $J^{*}_{2}$ est quadratique en $\mathbf{f}$. Cela signifie que la minimisation de ce critaire conduit à des équations linéaires en $\mathbf{f}$!
• Lorsque $\mathbf{f}$ est fixé, nous avons une expression exacte de $\mathbf{d}$.

De ces propriétés découle directement un algorithme que nous décrirons ultérieurement (Par.9.4.1).

8.2.4.6 Variable auxiliaire et processus de ligne.

Dans un article fondateur, les frères Geman proposaient déjà en 1984 [38] d'introduire une variable auxiliaire $\mathbf{d}$. De booléenne en 1984, elle devient continue et trouve sa justification en 1992 (Théorème de Geman&Reynolds [39]). Cette variable, appelée processus de ligne, sert à marquer les discontinuités. Ce processus de ligne est localisé sur une grille $S^{*}$ duale à la grille de l'image $S$. $S^{*}$ correspond à une grille de sites interstitiels ou sites contours On trouve une illustration de ces grilles Fig.9.4.

Figure 9.4: Image $S$ et grille duale $S^{*}$.

On constate en 2D que la grille $S^{*}$ comporte 2 fois plus de sites que la grille $S$. Il faut marquer les discontinuités dans les deux directions (une suivant $x$ et une suivant $y$). Le système de voisinage considéré correspond alors à l'union de ces deux grilles. Avec ce type de modèle, on considère non plus l'estimation du champ $\mathbf{f}$, mais celui de l'estimation simultanée de $\mathbf{f}$ et $\mathbf{d}$. On suppose que le couple ( $\mathbf{f},\mathbf{d}$) est markovien sur le système de voisinage composé de l'union $S\bigcup S^{*}$. On peut donc construire une énergie d'interaction basée sur les deux variables: $J^{*}_{2}(\mathbf{f},\mathbf{d})$ qui utilise conjointement les cliques associées à $\mathbf{f}$ et celles associées à $\mathbf{d}$. Cette énergie d'interaction peut s'écrire:

$$\displaystyle J^{*}_{2}(\mathbf{f},\mathbf{d})= J^{*}_{inter}(\mathbf{f}\vert\mathbf{d}) +J^{*}_{disc}(\mathbf{f})$$

$$= \sum _{n_{s}}\sum _{n_{t}\in \mathcal{V}}d_{n_{s},n_{t}}\{Df\}^{2}_{n_{s},n_{t}} +\sum _{n_{s}}\sum _{n_{t}\in \mathcal{V}}\psi (d_{n_{s},n_{t}})$$

Dans cette expression le premier terme $J^{*}_{inter}$ traduit l'interaction entre l'image et le processus de ligne. Alors que le deuxième terme $J^{*}_{disc}$ traduit l'énergie qu'il faut apporter pour créer la discontinuité telle qu'elle est. C'est dans ce terme qu'est introduit l'a priori sur le processus de ligne. On montre que, sous certaines conditions (théorème de Geman&Reynolds [39]), ces modèles définis avec un processus de ligne sont équivalents à un modèle sans processus de ligne dont la fonction de potentiel est définie par:

$$\varphi (u)=\inf _{d}\left( du^{2}+\psi (d)\right)$$

Autrement dit, certaines fonctions de potentiels $\varphi$ définissent un processus de ligne implicite. .

8.2.4.7 Illustration.

La notion de processus de ligne est bien utile car elle permet de mieux sentir la signification de la variable auxiliaire. En effet, pour le système de voisinage basé sur la distance, pour un potentiel d'interaction basé sur l'opérateur différentiel de premier ordre (gradient discret) et pour une fonction de potentiel répondant au théorème de Geman & Reynolds étendu, cette variable auxiliaire $\mathbf{d}$ peut être vue comme correspondant aux images transformées par la fonction de pondération des images de gradient . En 2D, pour une 4-connexité, nous obtenons 2 images représentant le processus de ligne: une, $\mathbf{d}_{x}$, suivant la direction associé à $x$ et l'autre $\mathbf{d}_{y}$ suivant l'axe des $y$ (Fig.9.5) obtenu à partir de $\mathbf{f}$ par:

$$(\mathbf{d}_{x},\mathbf{d}_{y})=\left( \frac{\varphi '\left( \fra... ...}}}{\delta }\right) }{2\frac{\{D_{y}f\}_{n_{s}}}{\delta }}\right) _{n_{s}\in S}$$ (8.7)

Evidemment, toujours en 2D, si nous envisageons un 8-connexité, nous obtiendrons 4 images: $\mathbf{d}_{x}$, $\mathbf{d}_{y}$ mais aussi des images représentant les directions diagonales $\mathbf{d}_{xy}$ et $\mathbf{d}_{yx}$. En 3D, il est nécessaire de rajouter $\mathbf{d}_{z}$, l'image transformée suivant la direction $z$. Par la suite, lorsque nous aurons besoin d'expliciter les variables auxiliaires, nous omettrons la normalisation par la hauteur des discontinuités afin d'alléger les notations.

Figure: Représentation des discontinuités.

[Fantôme. ]

[Estimation de $d_y$.]

[Estimation de $d_x$]