Master MSI

Année universitaire 2005/2006

Modélisation markovienne en
imagerie

Vincent Barra

Institut Supérieur d’Informatique, de Modélisation et de leurs

Applications
Campus des Cézeaux - B.P. 1025 - 63173 AUBIERE CEDEX
Table des matières
1 L’image 2

2 Modélisation probabiliste de l’image 3

3 Champs de Markov et de Gibbs 3

4 Echantillonnage des processus de Markov 5

4.1 Première approche : échantillonneur de Gibbs . . . . . . . . . . . 5
4.2 Seconde approche : algorithme de Métropolis . . . . . . . . . . . 6

5 L’algorithme du recuit simulé 6

5.1 Distribution de Gibbs avec température . . . . . . . . . . . . . . 7
5.2 Algorithme du recuit simulé . . . . . . . . . . . . . . . . . . . . 7
5.3 Algorithme Iterative Conditional Mode . . . . . . . . . . . . . . 8

6 Exemples de champs de Markov 9

6.1 Modèle d’Ising . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
6.2 Modèle de Potts . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
6.3 Modèle gaussien . . . . . . . . . . . . . . . . . . . . . . . . . . 9

7 Quelques applications en traitement d’images 10

7.1 Spécification de U(x | y) pour le problème de restauration d’images 11
7.2 Spécification de U(x | y) pour le problème de segmentation d’images 12

8 Estimateurs dans un cadre markovien 13

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
8.2 Modélisation bayésienne et fonction de coût . . . . . . . . . . . . 13
8.3 Estimateur MAP . . . . . . . . . . . . . . . . . . . . . . . . . . 14
8.4 Estimateur MPM . . . . . . . . . . . . . . . . . . . . . . . . . . 14
8.5 Estimateur TPM . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
8.6 Comparaison . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

9 Estimation de paramètres 16
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
9.2 Données complètes . . . . . . . . . . . . . . . . . . . . . . . . . 17
9.2.1 Méthode des codages . . . . . . . . . . . . . . . . . . . . 17
9.2.2 Pseudo vraisemblance . . . . . . . . . . . . . . . . . . . 18
9.2.3 Algorithme du gradient stochastique . . . . . . . . . . . . 19
9.3 Données incomplètes . . . . . . . . . . . . . . . . . . . . . . . . 20
9.3.1 Position du problème . . . . . . . . . . . . . . . . . . . . 20
 9.3.2 Estimation des paramètres d’attache aux données . . . . . 20
9.3.3 Algorithme EM . . . . . . . . . . . . . . . . . . . . . . . 21

10 Processus de bords 22
10.1 Processus implicites et explicites . . . . . . . . . . . . . . . . . . 23
10.2 Algorithmes de minimisation associés . . . . . . . . . . . . . . . 24

11 Quelques illustrations 24
11.1 Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
11.1.1 Imagerie radar . . . . . . . . . . . . . . . . . . . . . . . 24
11.1.2 Imagerie sismique . . . . . . . . . . . . . . . . . . . . . 25
11.1.3 Imagerie de test . . . . . . . . . . . . . . . . . . . . . . . 26
11.2 Restauration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
11.2.1 Imagerie de test . . . . . . . . . . . . . . . . . . . . . . . 27
11.2.2 Restauration avec prise en compte de processus de bords . 28

L’information véhiculée par une image est portée par bien d’autres données que
les seuls niveaux de gris. En effet, une image représente une suite d’objets dans
un espace 3D et des relations spatiales évidentes existent entre chaque site de
résolution de l’image (pixel). Le niveau de gris d’un pixel n’est donc souvent pas
significatif en lui-même, mais dans ses relations avec ses voisins. Cette notion
de relation de voisinage va nous permettre d’introduire un contexte markovien
pour l’analyse d’images, dans des domaines qui vont de la segmentation d’objets
à la restauration d’images bruitées. Dans la suite, un processus markovien sera
indifféremment dénoté comme tel, ou comme un champ de Markov, ou par son
abréviation anglaise MRF (Markov Random Field).

1 L’image
L’image est formée d’un ensemble fini S de sites s correspondant aux pixels, aux
groupes de pixels, aux régions dans l’image ou à tout autre attribut haut niveau.
S est donc essentiellement un réseau discret, fini de Z2 ou Z3 (même Z4 si l’on
considère des volumes évoluant dans le temps), muni d’une topologie V définie
par :

s∈/ Vs
Vs = t tels que :
t ∈ Vs ⇔ s ∈ Vt
A partir de cette topologie, un système de cliques peut être défini : une clique est
soit un singleton de S, soit un ensemble de sites tous voisins au sens de V. L’en-
semble des cliques relatif à V est noté C, et l’ensemble des cliques de cardinal k,

2
Ck .

Enfin, à chaque site s est associé un descripteur xs , représentant l’état du site :

si dans la plupart des cas il s’agit du niveau de gris, un pixel peut néanmoins être
caractérisé par un vecteur quelconque, prenant ses valeurs dans un espace E. Les
interactions locales entre descripteurs de sites voisins s’expriment alors comme
un potentiel de clique. Pour c ∈ C, un potentiel Uc est associé à c, dont la valeur
dépend des descripteurs des pixels constituant la clique. De même, le potentiel de
l’image est défini comme la somme des potentiels de toutes les cliques :
X
U= Uc
c∈C

et le potentiel local en un site est calculé comme la somme des potentiels de

toutes les cliques auxquelles il appartient :
X
Us = Uc
c∈C,s∈c

2 Modélisation probabiliste de l’image

A chaque site s on associe une variable aléatoire Xs prenant ses valeurs dans E.
Le descripteur xs est une réalisation de Xs . On définit alors un processus aléatoire
X = (X1 , · · · , XCard(S) ) prenant ses valeurs dans Ω = E Card(S) . L’image est
alors une réalisation de X, et la probabilité P (X = x) quantifie la vraisemblance
de l’image, les probabilités conditionnelles locales d’une valeur en un site permet-
tant quant à elles de mesurer le lien statistique entre un descripteur et le reste de
l’image. L’hypothèse markovienne va permettre d’évaluer ces quantités.

3 Champs de Markov et de Gibbs

Sous l’hypothèse markovienne, on a

P (Xs = xs | xr , r 6= s) = P (Xs = xs | xr , r ∈ Vs )

Autement dit, la probabilité d’observer un descripteur en s ne dépend pas de toutes

les valeurs des descripteurs de l’image, mais seulement de ceux situés dans le voi-
sinage de s, ce qui intuitivement est souvent satisfaisant lorsqu’on parle de for-
mation d’une image numérique. Cette hypothèse va permettre, via le théorème
d’Hammersley-Clifford, d’accéder aux expressions de ces probabilités condition-
nelles locales. Il faut pour cela introduire la notion de Champ de Gibbs.

3
Définition 1. La mesure de Gibbs d’hamiltonien U : Ω → IR est la probabilité
P définie sur Ω par :
e−U (x)
P (X = x) = X
e−U (y)
y∈Ω
avec X
U (x) = Uc (x)
c∈C
Le dénominateur est appelé fonction de partition de Gibbs, et permet de nor-
maliser la probabilité. Il est cependant quasi impossible à calculer, même pour des
configurations simples : Card(Ω) = 2262144 par exemple pour une image simple
512 × 512 à deux niveaux de gris !
Définition 2. Le champ de Gibbs de potentiel associé au système de voisinage V
est le processus aléatoire X dont la probabilité est une mesure de Gibbs associée
à V
L’énergie totale d’un champ de Gibbs se décompose donc sous forme d’une
somme d’énergies locales.
Le lien entre champ de Gibbs et processus Markovien est effectué par le théorème
fondamental suivant :
Théorème 1. : théorème d’Hammersley-Clifford
Soient S un ensemble fini ou dénombrable, V une topologie sur S et E un espace
d’états discret. Soit X un processus aléatoire à valeurs dans E Card(S) alors :
X est un champ de Markov relativement à V et P (X = x) > 0 pour tout x ∈ Ω si
et seulement si X est un champ de Gibbs de potentiel associé à V
Ainsi, la probabilité conditionnelle P (xs | Xr = xr , r 6= s) s’écrit :
P (X = x) e−U (xs ,xr ,r6=s)
P (Xs = xs | Xr = xr , r =
6 s) = =X
P (Xr = xr , r 6= s) e−U (e,xr ,r6=s)
e∈E

En définissant l’énergie locale du site s par :

X
Us (Xs = xs | Xr = xr , r ∈ Vs ) = Uc (xs , Xr = xr , r ∈ Vs )
c∈C/s∈C

on remarque que Us ne fait intervenir que les voisins de s, et l’énergie globale

U (x) s’écrit alors :
X X
U (x) = Uc (x) + Uc (x)
c∈C/s∈C
/ c∈C/s∈C
X
= Uc (x) + Us (Xs = xs | Xr = xr , r ∈ Vs )
c∈C/s∈C
/

4
 et en passant aux champs de Gibbs :
e−Us (Xs =xs |Xr =xr ,r∈Vs )
P (Xs = xs | Xr = xr , r 6= s) = X
e−Us (Xs =e|Xr =xr ,r∈Vs )
e∈E

L’expression obtenue ne fait donc intervenir que les potentiels de cliques dans
lesquelles est impliqué s (où au passage on retrouve l’hypothèse markovienne),
et le calcul du dénominateur de la fraction devient maintenant possible, contraire-
ment à l’expression générale de la fonction de partition de Gibbs.

4 Echantillonnage des processus de Markov

Ainsi, il suffit de définir séquentiellement une topologie, un système de cliques,
des potentiels de clique Uc et une fonction d’énergie U pour préciser le processus
de Markov attaché à l’image. L’ensemble de ces données permet, comme précisé
précedemment, d’accéder à la probabilité d’une configuration (i.e. pour le traite-
ment d’images, le fait d’observer l’image que l’on voit), et aux probabilités condi-
tionnelles locales.
Se pose donc maintenant le problème de la réalisation du tirage d’une configura-
tion donnée suivant la loi de Gibbs caractéristique (via le théorème d’Hammersley-
Clifford) du champ de Markov ainsi défini.

4.1 Première approche : échantillonneur de Gibbs

L’échantillonneur de Gibbs est un algorithme de relaxation probabiliste, qui
se propose de réaliser ce tirage par une méthode itérative, convergeant vers une
configuration suivant la loi de Gibbs.
A l’itération n, la construction de la réalisation courante se fait à l’aide des étapes
suivantes :

(i) Choix d’un site s (par tirage selon une loi uniforme sur E)
(ii) Calcul de la probabilité conditionnelle

e−U (xs ,xr ,r∈Vs )

P (Xs = xs | Xr = xr , r ∈ Vs ) = X
e−U (e,xr ,r∈Vs )
e∈E

suivant la configuration des sites r ∈ Vs à l’itération précédente

(iii) mise à jour du site s par tirage aléatoire selon P (Xs = xs | Xr = xr , r 6= s)

Algorithme 1: Echantillonneur de Gibbs

5
 Dans le cas où chaque site s est visité une infinité de fois, on peut montrer que les
réalisations générées sont des réalisations de la loi de Gibbs, indépendamment de
la configuration initiale utilisée pour l’algorithme. En pratique, on fera converger
l’algorithme en itérant le procédé ”un nombre suffisant de fois”, ou en détectant
le faible nombre de changement dans les mises à jour de sites d’une itération à
l’autre.

4.2 Seconde approche : algorithme de Métropolis

L’algorithme de Métropolis, issu de la physique statistique, est également un
algorithme de relaxation probabiliste. Le principe est donc itératif, et repose sur
les étapes suivantes :

(i) Choix d’un site s (par tirage selon une loi uniforme sur E)
(ii) Tirage aléatoire d’un descripteur δ ∈ E selon une loi uniforme
(iii) Calcul de la variation d’énergie induite par le changement d’état du site s (de
son état à l’itération précédente à l’état δ) :

∆U = Us (Xs = δ | Xr = xn−1
r , r ∈ Vs ) − Us (Xs = xsn−1 | Xr = xrn−1 , r ∈ Vs )
Si (∆U < 0) Alors
Le changement d’état est accepté pour le site s
Sinon
Le changement est accepté avec une probabilité e−∆U
Fin Si
Algorithme 2: Algorithme de Metropolis

Les critères de convergence sont similaires à ceux de l’échantillonneur de Gibbs,

mais cette dernière intervient cependant plus tardivement car les transitions ne
sont ici pas toujours acceptées. En revanche, l’itération est plus rapide que celle
de l’échantillonneur de Gibbs, car la variation d’énergie ne se calcule qu’entre
deux configurations données.

5 L’algorithme du recuit simulé

Dans les approches précédentes, chaque itération donnait accès à une réalisation,
qui convergeait vers une configuration suivant la loi de Gibbs. Dans ce paragraphe,
nous nous intéressons à la recherche des configurations les plus probables qui
correspondent à des états d’énergie minimale.

6
5.1 Distribution de Gibbs avec température
Définition 3. Une probabilité de Gibbs avec paramètre de température T > 0 est
une loi de probabilité de la forme :
U (x)
e− T
PT (X = x) = X U (y)
e− T
y∈Ω

Lorsque la température tend vers l’infini, il est facile de voir que PT converge
vers la loi uniforme sur Ω, et donc toutes les configurations sont équiprobables.

Lorsque T tend vers 0, on note Ω∗ = {x1 · · · xk } l’ensemble des configurations

atteignant l’énergie minimale U ∗ . Alors :
U (x)
e− T
PT (X = x) = X U (y)
e− T
y∈Ω
U (x)−U ∗
e− T
= X U (y)−U ∗
e− T
y∈Ω
U (x)−U ∗
e− T
= X U (y)−U ∗ X
e− T + 1
/ ∗
y ∈Ω y∈Ω∗

et deux cas se présentent alors :

– s’il existe i ∈ {1 · · · k} tel que x = xi , PT (x) = k1 ( il y a une infinité de
termes qui tendent vers 0 au dénominateur)
– sinon, U (x) − U ∗ > 0, l’exponentielle tend vers 0 lorsque T tend vers 0 et PT
aussi.
Lorsque la température est nulle, PT est donc uniformément distribuée sur les
configurations d’énergie minimale, et donc sur les configurations les plus pro-
bables.

5.2 Algorithme du recuit simulé

L’algorithme du recuit simulé exploite ce résultat. C’est une méthode itérative
qui procède, après avoir choisi une température initiale T 0 et une configuration
initiale x0 , selon les étapes suivantes :

7
 U (x)
(i)Simulation d’une configuration xn pour la loi de Gibbs d’énergie T
à partir
de la configuration (n − 1).
(i)Mise à jour de la température, (e.g. théorème d’Hajek :
c
Tn >
log(2 + n)

, c constante dépendant de la variation énergétique globale maximale sur Ω.

Algorithme 3: Algorithme du recuit simulé

La convergence est réalisée lorsque le nombre de changements entre deux confi-

gurations est faible. Contrairement aux algorithmes de l’échantillonneur de Gibbs
et de Métropolis qui échantillonnent selon la loi de Gibbs et qui sont en mesure de
donner toutes les configurations possibles, les images obtenues par recuit simulé
sont uniques et doivent en théorie correspondre aux minima globaux de l’énergie.

5.3 Algorithme Iterative Conditional Mode

Si l’algorithme du recuit simulé permet d’atteindre un optimum global, il peut
être coûteux en temps de calcul car il nécessite la génération d’un grand nombre
de configurations au fur et à mesure des itérations. Des algorithmes sous-optimaux
ont donc été proposés pour accélérer la convergence, et c’est le cas par exemple
de l’algorithme des modes conditionnels itérés.
Le principe est toujours itératif, mais cette fois-ci la modification du descripteur
xs du s se fait de manière déterministe.
On recherche ici une approximation du maximum a posteriori x b, résultat de la
convergence des itérations. Pour une itération k donnée, le déroulement de la mise
à jour de l’ensemble d’une configuration se fait, pour tout site s, en mettant à jour
xs par Arg max P (Xs = e | x br (k), r ∈ Vs ). L’algorithme converge lorsque le
e∈E
nombre de changements d’états des sites est faible (par exemple en pourcentage
de Card(S)).
Si on montre que l’énergie globale de x
b diminue, on ne peut affirmer que la conver-
gence aboutit à un optimum global, l’algorithme dépendant fortement de son ini-
tialisation.

8
6 Exemples de champs de Markov
6.1 Modèle d’Ising
L’espace d’état est un espace binaire {-1,1}, et la topologie est indifféremment
la 4- ou la 8-connexité dans une image 2D. Les potentiels de cliques d’ordre 1 (un
seul site) sont de la forme −Bxs , et les potentiels d’ordre 2 sont définis par :

−β si xs = xt
Uc={s,t} (xs , xt ) = −βxs xt =
β sinon
, où β est la constante de couplage entre s et t. L’énergie globale d’une configura-
tion est donc X X
U (x) = − βxs xt − Bxs
c={s,t} s∈S
.
β régularise le modèle : lorsque β > 0, les configurations d’énergie minimale
sont celles pour lesquelles les sites ont le même état, et si β < 0 les configurations
privilégiées sont celles où les sites sont d’états opposés.

6.2 Modèle de Potts

L’espace d’état est un espace de n valeurs {1,... n}, et la topologie est indifféremment
la 4- ou la 8-connexité dans une image 2D. Seuls les potentiels d’ordre 2 sont
définis par exemple par :

−β si xs = xt
Uc={s,t} (xs , xt ) = −βxs xt =
β sinon
Lorsque β > 0, les configurations les plus probables correspondent à des sites voi-
sins ayant des descripteurs égaux, ce qui dans le cas du traitement d’image donne
des réalisations avec de grandes zones homogènes (au sens de la similarité entre
descripteurs). La taille de ces régions est gouvernée par la valeur de la constante
de couplage.
Il est possible de raffiner ce modèle en considérant des valeurs de β différentes
suivant les directions explorées dans la clique (vertical/horizontal en 4-connexité
par exemple), et donc de privilégier la détection de zones homogènes direction-
nelles.

6.3 Modèle gaussien

Là encore, l’espace d’état est un espace de n valeurs {1,... n}, et la topologie est
indifféremment la 4- ou la 8-connexité dans une image 2D. Les potentiels d’ordre

9
1 sont de la forme α(xs − µs ), où µs est une donnée extérieure qui, dans le cas
du traitement d’images, peut être une moyenne de niveaux de gris attendue en
s lorsque les descripteurs sont les niveaux de gris. Les potentiels d’ordre 2 sont
quant à eux donnés par β(xs − xt ), et régularisent les configurations, favorisant
les faibles différences de descripteurs entre s et t lorsque β > 0. Le rapport αβ
pondère les influences respectives des termes de régularisation et d’attache aux
données, les valeurs absolues de ces paramètres décrivant le caractère équiréparti
ou localisé de la distribution.

7 Quelques applications en traitement d’images

Dans la suite, y est l’image observée, considérée comme la réalisation d’un
processus aléatoire Y . Le problème posé ici est de trouver une réalisation x,
modélisée par un processus markovien X, correspondant à :
– une image correspondant à y, dans laquelle le bruit a été corrigé : c’est le
problème de restauration d’images
– une image d’étiquettes, chaque étiquette correspondant à un objet présent dans
y : c’est le problème de segmentation d’images
Etant donné que l’on cherche à remonter à une image idéale de la scène (étiquettée
ou débruitée) connaissant seulement une observation y de cette scène, le problème
se résout à l’aide de chaı̂nes de Markov cachées.
Dans un cadre bayésien, on peut rechercher la configuration x b maximisant la pro-
babilité P (X = x | Y = y), qui s’écrit :
P (Y = y | X = x)P (X = x)
P (X = x | Y = y) =
P (Y = y)
– Le premier terme du numérateur décrit la probabilité d’observer y, sachant
que l’image idéale est x : il modélise donc l’acquisition de l’image. Sous
l’hypothèse (pas toujours vérifiée) d’indépendance des pixels,
Y
P (Y = y | X = x) = P (Ys = ys | Xs = xs )
s∈S

– le second terme du numérateur décrit la probabilité d’existence de l’image

idéale x, qui dans le cadre qui nous intéresse, répond à l’hypothèse marko-
vienne selon un système de voisinage V et un modèle φ dépendant de l’appli-
cation (restauration ou segmentation) :
e−U (x)
P (X = x) = X
e−U (z)
z∈Ω

10
 – le dénominateur est constant, et en particulier indépendant de x
Ainsi,

P (X = x | Y = y) = KP (Y = y | X = x)P (X = x)
= Keln(P (Y =y|X=x))−U (x)
= Ke−U (x|y)

où X X
U(x | y) = − ln (P (Ys = ys | Xs = xs )) + Uc (x)
s∈S c∈C

La distribution P (X = x | Y = y) est donc une distribution de Gibbs, et par

le théorème d’Hammersley-Clifford, le champ X conditionnellement à y est un
champ de Markov. Il est donc possible de simuler ce processus, par exemple à
l’aide du recuit simulé qui permet d’aboutir à x
b, qui maximise la probabilité a
posteriori (ou de manière équivalente qui minimise U(x | y).

7.1 Spécification de U(x | y) pour le problème de restauration

d’images
En supposant que la chaı̂ne d’acquisition de l’image introduit un bruit blanc
gaussien de variance σ 2 , la probabilité d’observer en chaque site s le descripteur
ys (supposé ici être le niveau de gris), sachant que le descripteur idéal est xs est :
1 (xs −ys )2
P (Ys = ys | Xs = xs ) = √ e− 2σ 2
2πσ 2
En faisant l’hypothèse que X est un processus markovien, la probabilité P (X =
x) peut être contrainte par des hypothèses de régularité, et avec un modèle de
potentiel d’ordre 2 :
X
−β φ(xs , xt )
s,t∈C2
e
P (X = x) = X
−β φ(zs , zt )
z∈Ω,s,t∈C2
e
U(x | y) vaut alors :
X (xs − ys )2 X
U(x | y) = +β φ(xs , xt )
s∈S
2σ 2 s,t∈C2

Le champ X conditionné par y est un champ de Gibbs, pour le même système

de voisinage que X. La constante β pondère le terme d’attache aux données qui

11
impose des descripteurs xs proches de ys par rapport au terme regularisation (uti-
lisation des cliques d’ordre 2) qui impose lui une image constituée de larges zones
homogènes.

7.2 Spécification de U(x | y) pour le problème de segmentation

d’images
Dans ce cas, X est défini sur l’espace E des étiquettes, et Y sur l’espace de ses
descripteurs. Le terme P (Y = y | X = x) traduit donc la probabilité d’observa-
tion de l’image y connaissant l’appartenance de chaque pixel à un objet présent
dans la scène. En supposant l’indépendance de chaque site, et en supposant que
ys ne dépend que de xs , on peut écrire :
Y
P (Y = y | X = x) = P (Ys = ys | Xs = xs )
s∈S

Les valeurs de probabilité conditionnelles sont par exemple données par la fréquen-
ce d’observation des niveaux de gris pour une classe donnée, et si chaque classe k
suit une distribution gaussienne de moyenne mk et d’écart-type σk :
(y −m ) 2
1 − s 2k
P (Ys = ys | xs = k) = p 2
e 2σ
k
2πσk

Si on fait encore une fois l’hypothèse markovienne sur X et si on se restreint aux

cliques d’ordre 2 X
−β φ(xs , xt )
s,t∈C2
e
P (X = x) = X
−β φ(zs , zt )
z∈Ω,s,t∈C2
e
et donc :
X (xs − mx )2 X
s
U(x | y) = + log (2π) σx s + β φ(xs , xt )
s∈S
2σx2s s,t∈C 2

Le terme d’ordre 1 exprime l’attache aux données, et le terme d’ordre 2 régularise

l’image. Souvent φ est le modèle de Potts, pour favoriser les grandes régions ho-
mogènes dans l’image.

12
8 Estimateurs dans un cadre markovien
8.1 Introduction
Nous avons vu précédemment comment il était possible d’utiliser le formalisme
markovien à des fins de restauration et de segmentation. On se situe alors dans
le cadre de données incomplètes (on parle aussi de champs de Markov cachés)
car la réalisation dont on dispose est une réalisation bruitée (ou plus généralement
vue au travers du système d’acquisition) du champ de Markov originel. En notant
Y le champ dont on observe une réalisation, et X le champ initial, l’objectif est
alors d’obtenir la meilleure réalisation xb de X connaissant l’observation y, autre-
ment dit, reconstruire x de manière optimale vis-à-vis d’un certain critère. Dans
le précédent paragraphe, nous nous étions intéressés à la réalisation maximisant la
probabilité a posteriori P (X = x|Y = y) et nous avions vu que le recuit simulé
permetttait d’accéder à cette réalisation. D’autres choix sont possibles, auxquels
correspondent d’autres méthodes de résolution.

8.2 Modélisation bayésienne et fonction de coût

En appliquant la règle de Bayes,

P (Y = y|X = x)P (X = x)
P (X = x|Y = y) =
P (Y = y)

Sous certaines hypothèses, la distribution a posteriori est une distribution de Gibbs

et le champ X conditionnellement à la donnée y est markovien. Cette propriété
n’est pas nécessaire pour les notions suivantes, mais elle sera primordiale pour les
algorithmes de résolution.
Le problème est alors de déterminer une estimation x b = Φ(y) optimisant un cer-
tain critère. L’estimation bayésienne construit une fonction L, estimant le coût de
remplacer x par Φ(y), qui vérifie :
– (∀x, x0 )L(x, x0 ) ≥ 0
– L(x, x0 ) = 0 ⇒ x = x0
L’estimateur optimal, i.e. la fonction Φ optimale est alors la fonction minimisant
X
IE(L(X, Φ(y))|Y = y) = L(x, Φ(y))P (x|y)

Suivant les fonctions de coût envisagées, on obtient différents estimateurs et différentes

méthodes de résolution associées.

13
8.3 Estimateur MAP
Si L(x, x0 ) = 1{x6=x0 } , L pénalise toute différence entre deux configurations, et
ce, quel que soit le nombre de sites en lesquels elles diffèrent. Ainsi
X
IE(L(X, Φ(y))|Y = y) = L(x, Φ(y))P (x|y)
= 1 − P (X = Φ(y)|y)
et la fonction Φ minimisant l’espérance de la fonction de coût est celle qui maxi-
mise la probabilité a posteriori. Il faut donc trouver x b, fonction de y, maximisant
P (X|y). On parle de l’estimateur MAP (maximum a posteriori) ou de maximum
de vraisemblance a posteriori. Les solutions algorithmiques associées à cet esti-
mateur sont le recuit simulé et l’ICM, utilisés avec la distribution a posteriori avec
paramètre de température.

8.4 Estimateur MPM

X X
Si L(x, x0 ) = L(xs , x0s ) = 1{xs 6=x0s } , L pénalise une configuration pro-
s∈S s∈S
portionnellement au nombre de différences entre deux configurations. Elle paraı̂t
donc plus naturelle que la fonction précédente. On montre alors que :
X
IE(L(X, Φ(y))|Y = y) = L(x, Φ(y))P (x|y)
X
= IE(L(Xs , Φ(y)s |y)
s∈S

On passe donc de la probabilité conditionnelle globale d’une configuration à la

probabilité conditionnelle en un site. Il s’agit d’une somme de termes positifs, et
par conséquent la fonction Φ optimale minimise en chaque site l’espérance condi-
tionnelle du coût local IE(L(Xs , Φ(y)s |y)). Ce résultat est valable pour toutes les
fonctions de coût définies par une somme de coûts en chaque site.
Comme précédemment, on obtient des estimateurs du maximum a posteriori lo-
caux, à calculer en chaque pixel contrairement à la recherche précédente qui était
globale. L’estimateur est appelé maximum de vraisemblance a posteriori local ou
maximum posterior marginal (pour maximum a posteriori de la marginale) abrégé
en MPM.
D’un point de vue algorithmique, la taille de l’espace des configurations ne per-
met pas un calcul direct des quantités P (xs |y). Aussi réalise-t-on en pratique des
approximations de type Monte-Carlo. En effet, supposons que l’on soit capable de
tirer des réalisations de X selon sa loi conditionnelle à y, et notons les x1 · · · xn . Il
est alors possible de calculer une approximation de l’estimateur MPM. Le tirage
des réalisations ne pose quant à lui pas de problème particulier (échantillonneur
de Gibbs et algorithme de Métropolis).

14
8.5 Estimateur TPM
X
Si L(x, x0 ) = (xs − x0s )2 , L est l’erreur quadratique et pénalise la somme
s∈S
des carrés des différences entre les deux configurations. Elle peut donc être plus
adaptée dans certains cas que les précédentes, puisqu’elle tient compte non seule-
ment du nombre de différences comme le MPM, mais aussi de leurs valeurs.
On montre alors que le minimum de L est atteint pour Φ telle que Φ(y)s =
IE(Xs |y), Cet estimateur consistant alors à prendre en chaque site la moyenne
conditionnelle locale donnée par la loi a posteriori, d’où le nom de TPM (Thre-
sholded Posteriori Mean).
D’un point de vue algorithmique, la démarche est similaire à celle effectuée dans
le paragraphe précédent pour le MPM. On approche l’espérance conditionnelle en
chaque site par la moyenne empirique en ce site de N échantillons tirés selon la
loi a posteriori.

8.6 Comparaison
Ces trois estimateurs sont ici comparés dans le cadre de la restauration. Dans le
cas du MAP, les résultats sont obtenus par ICM et par recuit simulé. L’image à
restaurer est une image bruitée par un bruit blanc gaussien. L’énergie a posteriori
utilisée s’écrit :
X (xs − ys )2 X
U(x | y) = +β φ(xs − xt )
s∈S
σ2 s,t∈C2

avec
u2
Φ(u) =
1 + u2
L’initialisation est donnée par l’image à restaurer. Pour tous les estimateurs, on
réalise 600 itérations. Dans le cas du recuit simulé, la température initiale est de
6.
Visuellement, le meilleur résultat est obtenu par le MAP du recuit simulé. Comme
l’image originale est une assez bonne initialisation, il n’y a pas de grandes différences
entre les algorithmes de recuit simulé et d’ICM pour l’estimateur du MAP. Ce
n’est pas vrai dès que l’initialisation s’écarte du résultat à obtenir et les différences
avec le recuit simulé peuvent être très importantes. Par ailleurs, on constate que
l’estimateur TPM, qui est par définition plus local, donne un résultat plus bruité et
moins régularisé. Cette analyse visuelle est confirmée par l’étude statistique qui
peut être effectuée sur des zones homogènes de l’image.

15
9 Estimation de paramètres
9.1 Introduction
Le problème de l’estimation de paramètres (encore appelés hyperparamètres
dans la littérature), revient très fréquemment en traitement d’image par champs
de Markov, et par exemple :
– on se donne une réalisation d’un champ de Markov associé au modèle d’Ising,
mais on ne connaı̂t pas ses paramètres.
– On veut généraliser ceci à une image de texture donnée dont on connaı̂t le
modèle sous-jacent (par exemple un modèle gaussien en 4-connexité) mais
pas les paramètres, qui sont du type : moyenne locale, variance locale, poids
de la régularisation locale. Quels sont-ils ? Leur connaissance pourrait bien
en effet servir à la classification d’images composées de zones texturées en se

16
 basant sur l’estimation locale de tels paramètres. On classifierait alors selon
les valeurs de ces attributs locaux.
– On veut segmenter une image, et pour cela apprendre les paramètres de chaque
classe, ainsi que le coefficient optimal du modèle de régularisation adapté à
cette tâche. On sait en effet que le résultat d’une segmentation par estima-
teur MAP par exemple dépend fondamentalement du poids respectif de la
régularisation par rapport à celui de l’attache aux données. Il faut donc là
aussi estimer ce poids d’une façon optimale dans un sens à définir.
L’ensemble forme un problème réputé difficile. Les solutions les plus utilisées se
décomposent en deux classes fondamentales :
– le cas des données dites complètes, correspondant aux deux premiers problèmes
cités plus haut : un échantillon d’une distribution de Gibbs est connu. Il s’agit
de remonter aux paramètres de cette distribution.
– le cas des données dites incomplètes. Là non seulement le résultat de traite-
ment est inconnu, mais les paramètres sont également à estimer.

9.2 Données complètes

Soit x une configuration observée relativement à une distribution de Gibbs donnée
Pθ , dont l’énergie associée peut s’écrire sous la forme d’une fonction linéaire d’un
paramètre θ , par exemple U (x) = θΦ(x), Φ étant un potentiel donné. Un principe
naturel en vue de la recherche de θ est d’écrire la vraisemblance de la donnée x :
exp−θΦ(x)
L(θ) = Pθ (x) =
Zθ

et de chercher par exemple la valeur de l’hyperparamètre θb maximisant cette vrai-

semblance. Le problème essentiel est que l’on ne sait en général pas calculer exac-
tement la fonction de partition Zθ . Même pour des modèles aussi simples et fon-
damentaux que ceux d’Ising et de Potts, le résultat (analytique) est obtenu après
des calculs excessivement compliqués. Dans les autres cas, on est amené :
– soit à effectuer des approximations de la fonction de partition globale au
moyen des fonctions de partition conditionnelles locales
– soit à employer des algorithmes itératifs (gradient stochastique) à partir de la
vraisemblance exacte, mais dont il s’agit alors de prouver la convergence ainsi
que le type d’optimum trouvé

9.2.1 Méthode des codages

Le principe de la méthode des codages est le suivant. Une fois défini un système
de voisinage pour un champ de Markov, on est capable de définir un certain
nombre de sous réseaux, chacun formé de sites/pixels indépendants les uns des

17
autres : chacun de ces sous réseaux est appelé un codage. Par exemple avec un
voisinage en 4 connexité il existe deux codages, et 4 codages différents dans le
cas de la 8-connexité.

L’estimation est alors posée dans le cadre de chaque codage. Pour un codage
donné, les différents sites/pixels le constituant sont indépendants les uns des autres
puisqu’ils ne sont pas des voisins pour le champ de Markov de départ. La proba-
bilité globale d’un codage se trouvera donc être le produit des probabilités indi-
viduelles de chacun des sites du codage. Or du fait de la structure de Markov du
champ de départ cette probabilité individuelle se trouve être la probabilité condi-
tionnelle locale du site/pixel dans le champ de Markov. Pour un codage Codn
donné, on a alors :
Y
Pθ ({Xs = xs }s∈Codn |{Xr = xr }s∈Cod/ n) = Pθ (Xs = xs |Vs )
s∈Codn

Etant donnée la structure des probabilités locales, la fonction de vraisemblance de-

vient calculable, puique les fonctions de partition conditionnelles locales le sont.
Dans le cas où la dépendance des énergies locales est linéaire vis à vis des pa-
ramètres, la log-vraisemblance associée,
X X
log (Pθ ({Xs = xs }s∈Codn |{Xr = xr }s∈Cod
/ n )) = −θ Uc (x) − log(Zs )
c∈C s∈Codn

est une fonction concave du (des) paramètre(s), comme somme de fonctions concaves.
Elle se prête donc bien à la recherche d’un optimum par une méthode classique
de type gradient. On peut également montrer qu’il s’agit dans ce cas d’un simple
problème de moindres carrés.

9.2.2 Pseudo vraisemblance

Expérimentalement, la méthode des codages n’est pas fiable. La méthode du
maximum de vraisemblance vrai paraı̂t quant à elle incalculable. Des algorithmes
ont cependant été étudiés pour tenter de résoudre ce problème. En fait on utilise
souvent une méthode intermédiaire qui a de bonnes propriétés : la méthode du
pseudo-maximum de vraisemblance.

18
Du maximum de vraisemblance vrai, cette méthode conserve l’idée de travailler
sur l’ensemble de l’image et non séparément sur des réseaux indépendants. De
la méthode de codage, elle garde l’idée de manipuler une fonction de vraisem-
blance produit des probabilités locales de chacun des sites/pixels. Cette fonction
est appelée pseudo-maximum de vraisemblance, et elle s’écrit
Y
P Lθ (X = x) = P (Xs = xs |Vs )
s∈S

ou, encore, en considérant le logarithme de cette fonction, et l’expression de la

probabilité locale de chaque site/pixel :
X X
logP Lθ (X = x) = −θ Uc (x) − log(Zs )
c∈C s∈S

log(Zs ) devient calculable, puisque relié à la fonction de normalisation de la pro-

babilité conditionnelle locale. Donc un raisonnement identique à celui du para-
graphe précédent conduit au fait que la log-pseudo-vraisemblance est une fonc-
tion concave des paramètres lorsque l’énergie en dépend de façon linéaire. Les
algorithmes usuels de type gradient ou gradient conjugué s’appliquent donc aussi
ici naturellement à la recherche de l’optimum unique. Il s’agit alors de qualifier la
valeur des paramètres obtenus par cette méthode par rapport à la valeur vraie, et
des auteurs ont montré que la méthode de la pseudovraisemblance est consistante
et convergente.

9.2.3 Algorithme du gradient stochastique

−θΦ(x)
La vraisemblance exacte du paramètre θ s’écrit L(θ) = Pθ (x) = e Zθ . La va-
leut de l’hyperparamètre θb satisfaisant au principe du maximum de vraisemblance
doit donc vérifier :
∂logPθ (x) ∂logZθ
= −Φ(x) − =0
∂θ ∂θ

On peut montrer que θb est unique et satisfait à l’équation stochastique

IEθb(Φ) = Φ(x)

La résolution est itérative, fondée sur des méthodes de type Newton-Raphson.

L’idée est de remplacer les grandeurs statistiques mises en jeu par leurs valeurs
empiriques approchées. Ainsi pour l’espérance du potentiel de régularisation Φ,
on prend sa moyenne empirique au cours d’une seule itération (valeur effective
obtenue !) d’un échantillonneur de Gibbs ou de Metropolis mené avec la valeur

19
courante du paramètre. Quant à la variance de ce potentiel, on l’estime encore plus
crûment par une grandeur positive fixée V . On peut montrer que le prix à payer
pour cette approximation est l’introduction d’un terme correctif supplémentaire
en 1/n + 1 dans le schéma iteratif.
On peut alors montrer que cet algorithme stochastique converge presque sûrement,
en termes de probabilité, vers la valeur optimale θb lorsque V est suffisamment
grand

9.3 Données incomplètes

9.3.1 Position du problème
Dans ce cas on connaı̂t une observation y , échantillon de la variable aléatoire Y ,
appelée incomplète ou dégradée, car reliée à une scène originale x, non-dégradée,
dont le champ aléatoire correspondant sera noté X. La relation entre y et x s’ef-
fectue via une loi de probabilité conditionnelle représentant l’attache aux données
explicitée la dépendance par rapport à un paramètre λ > 0 :

e−Uλ (y|x)
Pλ (Y = y|X = x) =
Zλ
On suppose également disposer d’une connaissance a priori sur la scène à re-
trouver x, que ce soit en segmentation ou restauration, via la distribution de Gibbs :

e−θΦ(x)
Pθ (x) =
Zθ
et on connaı̂t la forme de la distribution de Gibbs a posteriori

e−Uλ (y|x)−θΦ(x)
Pθ,λ (X = x|Y = y) =
Zθ,λ

L’estimation du meilleur paramètre de régularisation θ peut être réalisée, lorsque

la loi d’observation (attache aux données) est complètement connue, par une
méthode type gradient stochastique généralisé. Lorsque le loi d’observation est
paramétrique, il est nécessaire d’estimer les paramètres de l’attache aux données.
Nous détaillons ce cas particulier, car important, dans le cas du problème de la
segmentation d’images.

9.3.2 Estimation des paramètres d’attache aux données

On cherche ici à estimer λ. Dans le cas de la segmentation d’images, supposant
par exemple que le niveau de gris (descripteur) de chaque région 1 ≤ i ≤ m (site)

20
de l’image suive une loi gaussienne de paramètre µi , σi2 , alors λ = (µi , σi2 )1≤i≤m .
Pour l’exemple présenté ici, nous étudierons l’estimation de λ dans le cas gaussien
par raison de commodité. La probabilité des observations conditionellement aux
labels s’écrit donc :
(y −µ ) 2
Y 1 − s 2xs
Pλ (Y = y|X = x) = √ e 2σxs
s∈S
2πσxs

La distribution a posteriori s’écrit alors :

e−Eθ,λ (x|y)
Pθ,λ (X = x|Y = y) =
Zθ,λ

où
X  (ys − µx )2 
s
Eθ,λ (x|y) = 2
+ logσxs + θΦ(x)
s∈S
2σx s

La vraisemblance des paramètres d’attache aux données s’écrit donc

Zθ,λ
L(λ) = Pθ,λ (Y = y) = √
2π |S| Zθ
Si on souhaite maximiser L(θ) par rapport à l’un des paramètres de λ, on montre
en dérivant que
X
ys Pθ,λ (Xs = i)
s∈S
(∀i ∈ {1 · · · m})µi = X
Pθ,λ (Xs = i)
s∈S

et X
(ys − µi )2 Pθ,λ (Xs = i)
s∈S
(∀i ∈ {1 · · · m})σi = X
Pθ,λ (Xs = i)
s∈S

ces solutions sont complexes, et peuvent avantageusement être remplacées par

des formes itératives simples, en utilisant un algorithme de type EM (Expectation-
Maximisation).

9.3.3 Algorithme EM
L’algorithme EM est une méthode itérative. L’idée fondamentale est de considérer
que les données observées ne correspondent qu’à une connaissance partielle des

21
données complétées, et que la maximisation de la vraisemblance associée à ces
données est simple.
Supposons connus les paramètres θn et λn à l’étape n. Supposons avoir également
accès aux statistiques liées à la loi a posteriori des données cachées courantes
(étape E - Expectation). On définit alors
Q(θ, λ, θn , λn ) = IEθn ,λn (log(Pθ,λ (X = x|Y = y))
La partie optimisation (étape M - Maximisation) consiste à rechercher
(θn+1 , λn+1 ) = arg max Q(θ, λ, θn , λn )
θ,λ

ce qui peut être réalisé de manière séparée sur θ et λ, grâce à la séparabilité de la

loi jointe :
(
θn+1 = arg max IEθn ,λn (log(Pθ (X = x))
θ
λn+1 = arg max IEθn ,λn (log(Pλ (Y = y|X = x))
λ

On montre alors que la vraisemblance des paramètres croı̂t en fonction de n.

Le principe de l’algorithme EM est donc d’alterner une phase E et une phase M,
avec à chaque itération n :
X
ys Pθn ,λn (Xs = i)
s∈S
(∀i ∈ {1 · · · m})µn+1
i = X
Pθn ,λn (Xs = i)
s∈S

et X
(ys − µn+1
i )2 Pθn ,λn (Xs = i)
s∈S
(∀i ∈ {1 · · · m})σin+1 = X
Pθn ,λn (Xs = i)
s∈S
Ce sont bien les distributions a posteriori qui sont mises en jeu dans ces estima-
tions itératives. Elles sont donc à re-estimer à chaque itération n. Une manière bien
naturelle est de remplacer ces probabilités à l’étape courante du processus EM par
la fréquence empirique d’apparition des labels lors d’une série d’échantillonnages
de la distribution de Gibbs a posteriori courante menés à l’aide d’un échantillonneur
(Gibbs, Metropolis).

10 Processus de bords
Les modèles markoviens peuvent être utilisés, comme nous l’avons indiqué, en
restauration d’images : connaissant une réalisation y, on cherche une solution x

22
minimisant U(x|y) = U1 (x, y) + βU2 (x, y). Le terme de régularisation U2 per-
met d’agir sur les dérivées de la solution recherchée. Dans le cas d’un potentiel
quadratique, cela entraı̂ne un fort lissage de la solution au détriment des disconti-
nuités naturelles présentes dans l’image (bord des objets).
Il est possible de contraindre la dérivée première de la solution pour préserver ces
discontinuités naturelles, et gérer ainsi ces processus de bords.

10.1 Processus implicites et explicites

On s’intéresse ici aux cliques d’ordre 2 et donc à la prise en compte de disconti-
nuités naturelles entre deux sites voisins.
Un processus de bord, noté B est défini sur une grille duale de celle des sites et
est décomposé en deux processus, l’un, Bh modélisant les effets de bords horizon-
taux, l’autre, Bv les effets de bords verticaux. Le terme de régularisation s’écrit
alors pour un processus booléen (b=1 si discontinuité, 0 sinon) selon le modèle
weak membrane :
X  X 
(xs − xt )2 (1 − bvst ) + γbvst )+ (xs − xt )2 (1 − bhst ) + γbhst )
 
U2 (x, b) =
(s,t)∈Cv (s,t)∈Ch

En l’absence de discontinuités, on retrouve un potentiel quadratique qui lisse la

solution. En présence de discontinuités, la pénalité infligée à la discontinuité vaut
γ (qui peut prendre des valeurs différentes suivant Bh et Bv .
En réalité, il n’est pas nécessaire d’introduire explicitement un processus de bords
et un choix de fonction judicieux peut procurer le même résultat. Si φv (sx , st , bvst ) =
(xs −xt )2 (1−bvst )+γbvst (idem pour φh (sx , st , bvst )), on cherche (x, b) qui minimise
U2 (x, b), ce qui conduit à :
X X
min U2 (x, b) = min ψ(xs − xt ) + ψ(xs − xt )
x,b x
(s,t)∈C v (s,t)∈C h

avec
ψ(z) = min(z 2 , γ)
La quadratique est remplacée dans ce modèle par une quadratique tronquée. L’uti-
lisation d’un processus de bords explicite est dans ce cas équivalente à l’utilisa-
tion de la fonction ψ (on parle de processus de bords implicite). La valeur de γ
détermine alors à partir de quelle valeur du gradient on introduit une discontinuité
dans l’image.
Le choix d’un processus de bords implicite ou explicite va donner lieu à différents
algorithmes de minimisation. En effet, le recuit simulé permettant d’accéder au
minimum global de l’énergie peut être long et coûteux. Il peut dans de nombreux
cas être remplacé par un algorithme déterministe pour accél´erer la recherche de
la solution.

23
10.2 Algorithmes de minimisation associés
Parmi les algorithmes de minimisation associés aux processus de bords, citons
l’algorithme Graduated Non Convexity (GNC) qui minimise le critère sous sa
forme implicite. Le principe de l’algorithme consiste à approcher le critère par une
fonction convexe qui permet de définir une bonne solution initiale. Puis le critère
est graduellement modifié perdant sa propriété de convexité pour se rapprocher du
critère initial. A chaque étape, une solution est trouvée de manière déterministe
en utilisant pour initialisation la solution donnée par l’étape précédente. Si des
preuves de convergence peuvent être obtenues dans certains cas particuliers, cette
démarche n’assure cependant pas de trouver le minimum global pour toutes les
fonctions d’énergie.
D’autres approches utilisent un recuit par champ moyen MFA (Mean Field An-
nealing) sur le processus de bords sous sa forme explicite. L’algorithme met en
oeuvre une descente de température au cours de laquelle à chaque étape une solu-
tion (image restaurée et processus de bords) est estimée au sens du champ moyen.
Des expressions explicites sont dérerminées pour le processus de bords, tandis que
x est estimée de manière itérative.

11 Quelques illustrations
Dans cette section, nous donnons quelques résultats en images des champs de
Markov en segmentation et restauration.

11.1 Segmentation
11.1.1 Imagerie radar
L’image suivante (a) (source : INRIA, Sophia-Antipolis) représente une zone
rurale prise par le satellite SPOT. Sa segmentation en une dizaine de régions de
nature différente est donnée en (b).

24
11.1.2 Imagerie sismique
L’image suivante (a) est une image de réflexion sismique. Le but est ici de seg-
menter les couches géologiques, les champs de Markov fournissant les séparateurs
optimaux (b) reportés sur l’image originale (c).

25
11.1.3 Imagerie de test
L’exemple quivant montre, sur une image test classique, l’influence de la méthode
de recherche des configurations les plus probables qui correspondent à des états
d’énergie minimale. A partir d’une image originale bruitée, l’ICM et l’échantillonneur
de Gibbs fournissent deux résultats de segmentation bien distincts (problème de
minimum local pour ICM)

26
11.2 Restauration
11.2.1 Imagerie de test
La première ligne de l’image suivante montre une image originale, bruitée par
un bruit gaussien. L’image restaurée est donnée sur la seconde ligne, sans prise en
compte des effets de bords .

27
11.2.2 Restauration avec prise en compte de processus de bords
L’image suivante montre la restauration effectuée sur l’image (a) avec (d) et sans
(c) prise en compte de processus de bords. L’image (b) est un lissage de (a) qui
peut être effectué pour diminuer l’influence du bruit dans l’image

28