Quelques m ́ethodes math ́ematiques pour le traitement d’image
2 janvier 2009
Ce cours est une introduction `a la th ́eorie math ́ematique de traitement de l’image . Il est donc incomplet car les m ́ethodes dans ce domaine sont nombreuses et vari ́ees. On se polarisera sur les m ́ethodes variationnelles.
- Bergounioux Master 2 – 2008-2009
2
Chapitre 1
Introduction
1.1 Qu’est-ce qu’une image num ́erique ?
Une image num ́erique est compos ́ee d’unit ́es ́el ́ementaires (appel ́ees pixels) qui repr ́esentent chacun une portion de l’image. Une image est d ́efinie par :
– le nombre de pixels qui la compose en largeur et en hauteur (qui peut varier presque `a
l’infini), – l’ ́etendue des teintes de gris ou des couleurs que peut prendre chaque pixel (on parle de
dynamique de l’image).
- Les images binaires (noir ou blanc)
Exemple, images les plus simples, un pixel peut prendre uniquement les valeurs noir ou blanc. C’est typiquement le type d’image que l’on utilise pour scanner du texte quand celui ci est compos ́e d’une seule couleur. 2. Les images en teintes de gris
En g ́en ́eral, les images en niveaux de gris renferment 256 teintes de gris. Image `a 256 couleurs, simplement chacune de ces 256 couleurs est d ́efinie dans la gamme des gris. Par convention la valeur z ́ero repr ́esente le noir (intensit ́e lumineuse nulle) et la valeur 255 le blanc (intensit ́e lumineuse maximale). 3. Les images couleurs
S’il existe plusieurs modes de repr ́esentation de la couleur, le plus utilis ́e pour le manie- ment des images num ́eriques est l’espace couleur Rouge, Vert, Bleu (R,V,B). Cet espace couleur est bas ́e sur la synth`ese additive des couleurs, c’est `a dire que le m ́elange des trois composantes (R, V, B) donne une couleur.
3
4 CHAPITRE 1. INTRODUCTION
1.1.1 Quelques d ́efinitions
Pixels et niveaux de gris
Echantillonnage et quantification
1.1. QU’EST-CE QU’UNE IMAGE NUM ́ERIQUE ? 5
L’ ́echantillonnage est le proc ́ed ́e de discr ́etisation spatiale d’une image consistant `a asso- cier `a chaque zone rectangulaire R(x, y) d’une image continue une unique valeur I(x, y). On parle de sous- ́echantillonnage lorsque l’image est d ́ej`a discr ́etis ́ee et qu’on diminue le nombre d’ ́echantillons.
Une image num ́erique est une image ́echantillonn ́ee et quantifi ́ee. La quantificationd ́esigne la limitation du nombre de valeurs diff ́erentes que peut prendre I(x, y).
L’ ́echantillonnage est une ́etape fondamentale qui doit tenir compte du contenu informa- tionnel pertinent de l’image `a analyser. Sur l’exemple ci- dessous, en 1d, le signal ́echantillonn ́e « ressemble » `a une sinusoıde de fr ́equence 8 fois plus faible :
6 CHAPITRE 1. INTRODUCTION
Ce ph ́enom`ene appel ́e aliasing est encore pire en 2D, car il affecte la fr ́equence et la direction des structures p ́eriodiques. Imaginons par exemple qu’on souhaite ́echantillonner l’image cor- respondant aux bandes noires ci-dessous :
Avec un ́echantillonnage adapt ́e, l’image num ́erique fait apparaıtre des structures conformes `a l’information pr ́esente dans l’image :
Mais en consid ́erant seulement 1 ́echantillon sur 2, une structure diff ́erente apparaıt, dont l’ana- lyse (ici des bandes verticales, plus ́epaisses) ne sera pas conforme `a la r ́ealit ́e de l’objet :
Image originale Image sous- ́echantillonn ́ee
1.1. QU’EST-CE QU’UNE IMAGE NUMERIQUE ́? 7
La quantification peut ́egalement faire apparaıtre des distortions dans les images. Comme pour l’ ́echantillonnage, il existe des r`egles pour d ́eterminer la bonne quantification (le bon nombre de bits) pour coder les images num ́eriques. L’une d ́epend du capteur, et de sa capacit ́e effective `a observer des signaux de valeurs diff ́erentes : le rapport signal sur bruit. Le rapport signal sur bruit est d ́efini `a partir du rapport entre l’amplitude des niveaux de gris mesurables par le capteur nmax − nmin et le niveau du bruit, en gros l’ ́ecart-type σn de la perturbation al ́eatoire qui affecte les niveaux de gris. En prenant le logarithme, on a le nombre de bits utile au capteur pour coder les images.
Outre les capacit ́es du capteur, le nombre de bits r ́eellement n ́ecessaires pour coder une image varie d’une image `a l’autre, en fonction de leur contenu informationnel. Ce nombre d ́epend de l’entropie, d ́efinie `a partir de la distribution des niveaux de gris de l’image (statis- tique).
E = − ∑i≤N
pi log2(pi) ,
o`u N est le nombre de niveaux de gris pr ́esents, pi est la proportion (0 < pi < 1) de points de l’image ayant pour niveau de gris i. Cette grandeur repr ́esente le nombre moyen de bits par pixel n ́ecessaires pour coder toute l’information pr ́esente. Elle est utilis ́ee dans les techniques de compression sans perte pour adapter le volume de donn ́ee des images `a leur contenu infor- mationnel.
8 CHAPITRE 1. INTRODUCTION
Profil – Histogramme
Il est possible de tracer un trait sur le flanc du z`ebre de l’image en niveaux de gris ci-dessous et obtenir le profil correspondant, c’est `a dire le niveau de gris de chaque point ou pixel travers ́e par la ligne :
L’histogramme de l’image en niveau de gris ci-dessous permet d’ ́etablir une stricte corr ́elation entre les donn ́ees num ́eriques codant la nuance de gris et la position des pixels de l’image.
1.2 Qu’est-ce que le traitement d’image ?
Les pages qui suivent sont extraites des cours en ligne [BMT, M].
1.2. QU’EST-CE QUE LE TRAITEMENT D’IMAGE ? 9
1.2.1 Quelques aspects du Traitement d’Image
– Filtrage / d ́econvolution (ou filtrage inverse) – Compression – Segmentation – Restauration / reconnaissance – Reconstruction tomographique
10 CHAPITRE 1. INTRODUCTION
1.2. QU’EST-CE QUE LE TRAITEMENT D’IMAGE ? 11
Dans ce cours nous ́evoquerons successivement trois aspects fondamentaux : – Filtrage
L’outil math ́ematique essentiel pour le filtrage est la transformation de Fourier (ou toute autre transformation du mˆeme type comme la transformation en ondelettes) – Segmentation
La segmentation fait intervenir des notions d’optimisation, des outils g ́eom ́etriques et des ́equations aux d ́eriv ́ees partielles – Restauration
Les mod`eles variationnels en restauration utilisent de l’optimisation et de l’analyse fonc- tionnelle fine (Banach, th ́eorie de la mesure).
12 CHAPITRE 1. INTRODUCTION
1.2.2 Applications
– Robotique – Industrie
– Assemblage, reconnaissance de pi`eces – Contrˆole de qualit ́e – V ́ehicule autonome – etc … – T ́el ́ed ́etection
– M ́et ́eo – Cartographie – Analyse des ressources terrestres – Astronomie – Restauration – etc… – Applications militaires – Guidage de missile – Reconnaissance (a ́erienne, sous-marine, etc …) – etc … – Imagerie m ́edicale
– Tomographie – Aide au diagnostic – Comptage (nombre de cellules) – Suivi de formes anatomiques
1.2. QU’EST-CE QUE LE TRAITEMENT D’IMAGE ? 13
– Restauration – etc… – S ́ecurit ́e
– Reconnaissance (d’empreintes, visages, signatures) – D ́etection de mouvement – etc…
14 CHAPITRE 1. INTRODUCTION
Chapitre 2
Traitement ponctuel des images num ́eriques
On s’int ́eresse d’abord aux traitements ponctuels qui consistent `a faire subir `a chaque pixel une correction ne d ́ependant que de sa valeur. On trouve dans cette cat ́egorie, les fonctions de recadrage ou d’ ́egalisation de dynamique, de binarisation …
Sauf mention particuli`ere, nous supposerons dans ce qui suit des images comportant N2 pixels cod ́es sur 256 niveaux de gris diff ́erents.
2.1 Correction ponctuelle d’une image -Recadrage de dynamique
Il s’agit d’une transformation du type f = t(f) qui permet de modifier la dynamique des niveaux de gris dans le but d’am ́eliorer l’aspect visuel de l’image. `A un niveau de gris f de l’image originale correspond le niveau t(f) dans l’image transform ́ee. On fait subir `a chaque pixel un traitement ne d ́ependant que de sa valeur. La transformation t(f) peut ˆetre r ́ealis ́ee en temps r ́eel sur l’image en cours d’acquisition `a l’aide d’une table de transcodage dans la- quelle les valeurs de la transformation sont m ́emoris ́ees. Un adressage de cette m ́emoire par une donn ́ee f fournit directement la valeur t(f).
2.1.1 Transformation de recadrage
On suppose une image de d ́epart pr ́esentant un histogramme concentr ́e dans l’intervalle [a, b]. Les valeurs a, bcorrespondent aux niveaux de gris extrˆemes pr ́esents dans cette image. Le recadrage de dynamique consiste `a ́etendre la dynamique de l’image transform ́ee `a l’ ́etendue totale [0,255]. La transformation de recadrage est donc une application affine qui s’ ́ecrit :
15
16 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
Original Recadrage : a = 30, b = 200
Variantes pour le rehaussement des contrastes
Les types de correction donn ́es ci-dessous permettent d’accentuer le contraste dans une plage pr ́ecise de niveau.
Dilatation de la dynamique des zones sombres Dilatation de la dynamique des zones claires Fonction de rehaussement de contraste.
2.1. CORRECTION PONCTUELLE D’UNE IMAGE -RECADRAGE DE DYNAMIQUE 17
t(f) =
abf pour 0 ≤ f ≤ a (255 − b)f + 255(b 255 − a − a)
pour a ≤ f ≤ 255
Dilatation de la dynamique des zones sombres
Dilatation de la dynamique des zones claires
Original Histogramme
18 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUMERIQUES ́2.1.2 Egalisation ́de l’histogramme
L’histogramme d’une image est rarement plat ce qui traduit une entropie non maximale. La transformation d’ ́egalisation est construite de telle fac ̧on que l’histogramme de l’image trans- form ́ee soit le plus plat possible. Cette technique am ́eliore le contraste et permet d’augmenter artificiellement la clart ́e d’une image grˆace `a une meilleure r ́epartition des intensit ́es relatives.
Fonction d’aplatissement continue
Consid ́erons l’histogramme continu h(f) donn ́e ci-dessous. En notant f = t(f), l’histo- gramme ́egalis ́e h(f ) doit s’approcher de la forme id ́eale d ́ecrite ci-dessous.
Deux surfaces ́el ́ementaires en correspondance dans les histogrammes initiaux et ́egalis ́es, pr ́esentent le mˆeme nombre de points ce qui permet d’ ́ecrire :
f = t(f) = 256N2
∫ f0 h(s)ds .
Histogramme d’origine Histogramme plat id ́eal Fonction id ́eale d’ ́egalisation d’un histogramme
Fonction d’aplatissement discr`ete
En remplac ̧ant l’int ́egration continue par une somme, on obtient la transformation d’ ́egalisation discr`ete suivante :
f = t(f) = 256N2
∑fi=0
h(i) .
h(i) .
2.1. CORRECTION PONCTUELLE D’UNE IMAGE -RECADRAGE DE DYNAMIQUE 19
Image ́egalis ́ee Histogramme ́egalis ́e
2.1.3 Binarisation
Le but de la binarisation d’une image est d’affecter un niveau uniforme au pixels pertinents et d’ ́eliminer les autres. Seuillage
Le seuillage consiste `a affecter le niveau 255 aux pixels dont la valeur est sup ́erieure `a un seuil S et 0 le niveau aux autres. Le graphe de la transformation correspondante est le suivant
Original Histogramme
Fonction « seuillage »
20 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
Seuillage `a 70
Extraction d’une fenêtre d’intensité
Avec la transformation d ́ecrite ci-dessous, la nouvelle image ne visualise que les pixels dont le e niveau d’intensit ́e appartient `a l’intervalle [a, b]. Sous r ́eserve d’une connaissance a priori de la distribution des niveaux de gris des objets de l’image originale, cette technique permet une segmentation d’objets particuliers de l’image.
Original Histogramme
Fonction « fenˆetre d’intensit ́e »
2.2. ANAYSE DE LA NETTET ́E D’UNE IMAGE NUM ́ERIQUE 21
Seuillage avec fenˆetre d’intensit ́e entre 30 et 100
2.2 Anayse de la nettet ́e d’une image num ́erique
La focalisation automatique, ou mise au point, des syst`emes de prise de vue traditionnels (appareils photo,, graphiques cam ́eras analogiques … ) exploite g ́en ́eralement un t ́el ́em`etre qui mesure avec pr ́ecision la distance appareil-plan objet. Les syst`emes de vision num ́eriques ac- tuels (appareils et cam ́eras num ́eriques … ) r ́ealisent automatiquement leur mise au point `a par- tir d’une analyse de l’image restitu ́ee. Pour cela un crit`ere de nettet ́e de l’image est d ́etermin ́e `a partir des valeurs f(i, .j) := fij de l’intensit ́e (luminance) des pixels.
Un grand nombre de crit`eres a ́et ́e propos ́e. Parmi les plus utilis ́es figurent les crit`eres bas ́es sur l’analyse de la distribution des niveaux d’intensit ́e pixel et ceux mesurant le contenu spec- tral de l’image.
La premi`ere cat ́egorie repose sur le fait que la d ́efocalisation engendre l’uniformisation des niveaux de gris, ou ce qui revient au mˆeme, que l’image nette pr ́esente l’histogramme le plus large.La seconde cat ́egorie repose sur le principe similaire qu’une image focalis ́ee contient plus de fr ́equences spatiales ́elev ́ees qu’une image floue.
Dans toutes ces m ́ethodes, on recherche la distance de mise au point qui assure la majora- tion du crit`ere
2.2.1 Exemples de crit`eres de nettet ́e
Crit`eres de nettet ́e bas ́es sur l’analyse de l’histogramme de l’image
On note hk la distribution statistique (histogramme) des niveaux k d’intensit ́e des pixels.
22 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUMERIQUES
́Crit`ere math ́ematique
M ́ethode de F = ∑k>T
k hk
Mendelsohn T est un param`etre choisi et Mayall proche du niveau de gris moyen de l’image
M ́ethode de F = ∑k>T(k − T)hk et T =
∑∑i
∑j fij∆ij i ∑j ∆ij Mason et Green ∆ij = 2[(fi,j+1 − fi,j−1)2 + (fi+1,j − fi−1,j)2] + (fi−1,j−1 − fi+1,j+1)2 + (fi−1,j+1 − fi+1,j−1)2
Variance des F =
∑[∑i i
∑j fij]2
niveaux de gris
N4
Crit`eres de nettet ́e bas ́es sur l’analyse spectrale de l’image
Crit`ere math ́ematique
Norme l1 du gradient1 F = ∑i,j
∑j fij 2N2 −
(|fi,j − fi,j−1| + |fi,j+1 − fi,j|)
Norme l2 du gradient F = ∑i,j
(fi+1,j − fi,j)2 + (fi,j+1 − fi,j)2
Remarques g ́en ́erales
– Les tests des crit`eres de nettet ́e r ́ealis ́es sur un grand nombre d’images diff ́erentes montrent
que leurs performances d ́ependent du type et du contenu de l’image analys ́ee ; – En r`egle g ́en ́erale, l’extremum des crit`eres est peu prononc ́e lorsque l’image contient peu
de d ́etails
2.3. CORR ELATION ́D’IMAGES NUMERIQUES ́23
– L’algorithme de Brenner est souvent utilis ́e car il pr ́esente g ́en ́eralement une bonne sen-
sibilit ́e et la charge de calculs qu’il n ́ecessite est raisonnable.
2.3 Corr ́elation d’images num ́eriques
Nous consid ́erons deux fonctions bidimensionnelles discr`etes f(i, j) etg(i, j) avec 1 ≤ i, j ≤ N.
Ces fonctions sont repr ́esentatives de deux images num ́eriques monochromes comportant N2 pixels. Leurs versions centr ́ees et r ́eduites sont d ́efinies respectivement par
fc(i, j) = f(i, j) σf − μf
et gc(i, j) = g(i, j) − μg
σg
o`u μ et σ sont respectivement la moyenne et l’ ́ecart type de chaque fonction.
La fonction d’intercorr ́elation entre fc et gc est d ́efinie par la relation :
φfg(k, l) = N12
∑i,j
fc(i, j)gc(i − k, j − l) . (2.3.1)
Il est souvent plus pratique de travailler sur des images dont les valeurs sont centr ́ees et r ́eduites. Ceci permet d’une part d’ ́eliminer de φfg(k, l) le carr ́e des moyennes des images, d’autre part de normaliser les fonctions de corr ́elation.
Pour ́eviter les effets de bord qui introduisent un biais dans le calcul de φfg(k, l) il convient : – soit de s’assurer que la double sommation est r ́ealis ́ee sur des valeurs dont les coor- donn ́ees ne sortentjainais de l’image. Cela revient `a d ́efinir un cadre d’int ́egration dont le format, inf ́erieur `a celui de l’image, est choisi en fonction du domaine de calcul de φfg(k, l); – soit de remplacer dans l’expression (2.3.1) le diviseur N2 par la valeur variable M qui
d ́epend de k et l selon la relation : M = (N − |k|)(N − |l|).
2.3.1 Application `a la reconnaissance d’empreintes digitales
Les fonctions d’intercorr ́elation bidimensionnelles peuvent ˆetre utilis ́ees en reconnaissance d’images. La comparaison de l’empreinte digitale d’un suspect avec celles contenues dans un fichier en est un exemple. les deux images repr ́esent ́ees ci-dessous repr ́esentent les empreintes digitales de deux individus que nous appellerons F et G.
24 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
Les caract ́eristiques g ́en ́erales de ces deux images sont donn ́ees dans le tableau suivant :
Individu F Individu G Format 256 x 256 256 x 256 Codage en niveaux de gris 0 `a 255 0 `a 255 Moyenne μ 132,8 127 ́Ecart type σ 33,8 36,4
Les valeurs des deux images sont centr ́ees et r ́eduites. Pour ́eviter les effets de bords, l’in- tercorr ́elation est r ́ealis ́ee entre deux zones carr ́ees de 151 x 151 pixels (voir figure ci-dessous)
La fonction est estim ́ee dans la gamme −20 ≤ k, l ≤ 20 grˆace `a la formule suivante
φfg(k, l) = 1511
2
∑ 200i,j=50fc(i, j)gc(i − k, j − l) . (2.3.2)
∑ 200i,j=50fc(i, j)gc(i − k, j − l) . (2.3.2)
La repr ́esentation 3D `a gauche dans la figure ci-dessous est celle de la fonction d’autocorr ́elation φfg(k, l) effectu ́ee sur l’empreinte F. Son maximum est ́egal `a 1 au centre du graphique qui cor- respond ici `a la position k = l = 0. La repr ́esentation de droite est la fonction d’intercorr ́elation entre les empreintes F et G.
2.3. CORR ́ELATION D’IMAGES NUM ́ERIQUES 25
Ses valeurs ne d ́epassent pas 0.15 ce qui prouve que les deux empreintes ne sont pas iden- tiques.
Dans la r ́ealit ́e, les empreintes ne sont pas toutes obtenues dans les mˆemes conditions de positionnement. Il est n ́ecessaire d’ajouter un param`etre de rotation d’image `a la fonction d’in- tercorr ́elation ce qui peut alourdir consid ́erablement les calculs.
2.3.2 Application `a l’analyse de la texture d’une image
Par d ́efinition la texture d’une image est la structure spatiale sur laquelle sont organis ́es les pixels. Les relations structurelles peuvent ˆetre :
– d ́eterministes : c’est le cas de la r ́ep ́etition quasi p ́eriodique d’un motif de base (brique,
carrelage, mailles de tissu … ) ; – al ́eatoires : il n’y a pas de motif de base (sable, mur cr ́epi…) La fonction d’autocorr ́elation bidimensionnelle est bien adapt ́ee pour mettre en ́evidence les propri ́et ́es de la texture d’une image. La figure suivante repr ́esente deux exemples de texture.
Nous donnons ci-dessous la fonction d’auto-corr ́elation des textures « sable » et « tissu ».
26 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
2.4 Transformation de Hough
[,a transformation de Hough est utilis ́ee pour d ́etecter de mani`ere syst ́ematique la pr ́esence de relations structurelles sp ́ecifiques entre des pixels dans une image Par exemple une image repr ́esentant un site urbain est compos ́ee de nombreuses lignes droites (immeubles, fenˆetres) en revanche, une vue de campagne en est quasiment d ́epourvue .
Hough propose une m ́ethode de d ́etection bas ́ee sur une transformation d’image permet- tant la reconnaissance de structures simples (droite, cercle, … ) liant des pixels entre eux. Pour limiter la charge de calcul, l’image originale est pr ́ealablement limit ́ee aux contours des objets puis binaris ́ee (2 niveaux possibles pour coder l’intensit ́e pixel).
2.4.1 Principe de la m ́ethode pour la recherche de ligne droite
Supposons que l’on suspecte la pr ́esence d’une droite ∆ reliant un certain nombre de pixesl Pi. Soit le pixel P1 de coordonn ́ees (x1,y1). Une infinit ́e de droites d’ ́equation : y1 = ax1 + b peuvent passer par P1. Cependant, dans le plan des param`etres ab l’ ́equation qui s’ ́ecrit b = −ax1 + y1 devient une droite unique D1 (voir figures ci-dessous).
Un second pixel P2 = (x2,y2) permet de d ́efinir une seconde droite D2 du type b = −ax2+y2 dans le plan ab.L’intersection de D2 avec D1 fournit le couple a ,b ) qui sont les param`etres
2.4. TRANSFORMATION DE HOUGH 27
de la droite recherch ́ee ∆ dans le plan image. Ainsi, tous les pixels Pi qui sont align ́es sur ∆ poss`ede une droite Di dnas le plan ab qui coupe les autres au point particulier (a ,b ).
2.4.2 Probl`eme pour la recherche de ligne verticale Repr ́esentation normale d’une
droite
Dans le plan image, une droite verticale poss`ede des param`etres a etb infinis, ce qui ne permet pas d’exploiter la m ́ethode pr ́ec ́edente. Pour contourner ce probl`eme, on utilise la repr ́esentation normale des droites. Cette repr ́esentation, de param`etres ρetθ, ob ́eit `a l’ ́equation
xcosθ + y sinθ = ρ . (2.4.3)
Cette ́equation repr ́esente le produit scalaire (projection) entre les vecteurs V =
( cosθ sinθ
)
et
−−→OM =
( xy
). Ainsi l’ensemble des points M d’une droite se projettent sur un mˆeme vecteur
particulier de coordonn ́ees polaires (ρ,θ) . Les cas de droites horizontales et verticales sont illustr ́es ci-dessous :
2.4.3 Transformation de Hough pour la d ́etection de droites dans une image binaire
Domaine de variation des param`etres θ et ρ
Il est `a noter que si (ρ,θ) sont les param`etres d’une droite (−ρ,θ + π) le sont ́egalement. Par cons ́equent l’intervalle [0,π[correspond au domaine de variation complet du param`etre θ. Les coordonn ́eescart ́esiennes x ety des pixels d’une image num ́erique sont g ́en ́eralement positives, l’origine ́etant plac ́ee `a un sommet de l’image. En consid ́erant une image carr ́ee comportant N × N pixels, les valeurs de ρ calcul ́ees par l’ ́equation (2.4.3) peuvent ˆetre major ́ees par √2N.
Quadrillage du plan θρ
On suspecte dans l’image la pr ́esence d’une ou plusieurs structures caract ́eris ́ees par l’ali- gnement d’un certain nombre de pixels. Soit les deux intervalles [θmin,θmax] dans lesquels on cherche `a calculer les param`etres de la repr ́esentation normale d’une ou plusieurs droites.
28 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
Le plan θρ limit ́e aux intervalles de recherche, est subdivis ́e en cellules par quantification des param`etres θ et ρ avec les pas respectifs ∆θ et ∆ρ (indexation par les indices p et q). La figure ci-dessous donne un exemple de quadrillage du plan θρ. Un tableau A(p, q) dont les valeurs sont initialement mises `a z ́ero, est associ ́e `a ce quadrillage
Proc ́edure de transformation
– Pour chaque pixel non nul Pi(xi,yi) de l’image, on balaye l’axe des θ de θmin `a θmax
suivant la trame du tableau et pour chaque θp on r ́esout l’ ́equation (2.4.3) ; – le r ́esultat ρ obtenu est arrondi `a la valeur ρq du tableau la plus proche ; – si `a la valeur θp correspond la solution ρq la valeur A(p, q) est incr ́ement ́ee d’une unit ́e. `A la fin de cette proc ́edure, A(p, q) = M signifie que M points de l’image sont align ́es sur la droite de param`etres approximatifs (θp,ρq). Apr`es transformation compl`ete, le tableau peut ˆetre repr ́esent ́e sous la forme d’une image en niveaux de gris (exemple plus loin) ou celle d’un graphique 3D. La d ́ecision sur la d ́etection de droites peut ˆetre prise apr`es recherche des coor- donn ́ees des valeurs significatives du tableau.
La charge de calcul n ́ecessaire pour r ́ealiser la transformation de Hough est importante. Elle d ́epend du nombre de param`etres recherch ́es : – recherche de droites : 2 param`etres ; – recherche de cercles : 3 param`etres.
Exemple pour la recherche de droite : pour N pixels non nuls de l’image binaire et K sub- divisions de l’axe θ.ilya NK d ́eterminations de l’ ́equation (2.4.3). Une r ́eduction du temps de calcul peut ˆetre obtenue par l’utilisation de tables pr ́eenregistr ́ees des conversions sinθp et cosθp
Exemple
Nous consid ́erons dans cet exemple une image binaire comportant 6 x 6 pixels dont les valeurs sont donn ́ees dans le tableau suivant
2.4. TRANSFORMATION DE HOUGH 29
Deux droites ∆1 et ∆2 comportant chacune 6 pixels align ́es apparaissent clans cette image. La proc ́edure de calcul d ́ecrite au paragraphe pr ́ec ́edent est utilis ́ee avec les param`etres sui- vants : θp varie de 0 `a 180o par pas de 1o. Plus et ρmax = g ́en ́eralement √(N1 si on consid`ere − 1)2 + (N2 − 1)2. une Si on image de tailleN1×N2, ρmin = −√(N1 − 1)2 + (N2 − 1)2
discr ́etise l’espace de Hough avec une r ́esolution de hθ (par exemple 1o) pour θ et une r ́esolution de hρ pour ρ le tableau de r ́ef ́erence est donn ́e par
θp = −90 + phθ avec p ∈ {0,Nθ} et − 90 + Nθhθ = 90 ,
ρq = ρmin + qhρ avec q ∈ {0,Nρ} et ρmin + Nρhρ = ρmax.
Le programme it ́eratif suivant permet le calcul de la transformation de Hough :
Mi,j sont les valeurs binaires de l’image originale : Mi,j ∈ {0,1} Pour p de 1 `a Nθ
Pour q de 1 `a Nρ
Ap,q = 0 Pour i de 1 `a N1
Pour j de 1 `a N2
Pour p de 1 `a Nθ q = Arrondi
{ 1hρ
[icos( p180π)
+ j sin( p180π)]
− ρmin} Ap,q = Ap,q + 1 si Mi,j = 0
30 CHAPITRE 2. TRAITEMENT PONCTUEL DES IMAGES NUM ́ERIQUES
Le mˆeme type de programme appliqu ́e `a l’image binaire (64 × 64) donne le r ́esultat ci-dessous :
Chapitre 3
Filtrage
3.1 Filtrage lin ́eaire des signaux 1D
L’exemple le plus courant de signal 1D est le signal sonore.
D ́efinition 3.1.1 (Filtre) Un filtre est un syst`eme lin ́eaire continu et invariant. En d’autres termes, on se donne deux espaces normés X et Y ainsi qu’un opérateur linéaire continu, invariant par translation A : X → Y.
Supposons qu’ on puisse d ́efinir la transform ́ee de Fourier des signaux d’entrée f et de sortie g (par exemple si ces fonctions sont dans L1(R) ou dans L2(R) et que le signal de sortie le filtrage de f par un filtre lin ́eaire est tel queˆg = H f ˆo`u ˆg d ́esigne la transform ́ee de Fourier de g. La fonction H est appel ́ee fonction de transfert du filtre. Si la fonction de transfert H est dans L2(R) ∩ L∞(R), elle admet une transform ́ee de Fourier inverse h = F−1H, born ́ee et `a d ́ecroissance rapide, continue (sauf peut-ˆetre `a l’origine) et la relation
ˆg = h ˆf ˆentraıne
g = h ∗ f .
La r ́eponse `a l’entr ́ee f est donc un produit de convolution de l’entr ́ee avec une fonction fixe h que l’on appelle r ́eponse impulsionnelle.
D ́efinition 3.1.2 On appelle r ́eponse indicielle d’un filtre sa réponse `a l’échelon unité (fonction de Heaviside) :
h1(t) =
∫ t−∞ h(s)ds .
31
32 CHAPITRE 3. FILTRAGE
Notons que la plupart des filtres courants sont des filtres de convolution. Il est courant de d ́efinir un filtre par la fac ̧on dont il modifie les fr ́equences du signal d’entr ́ee, c’est-`a-dire par sa fonction de transfert H(λ) puisque les fr ́equences de l’entr ́ee f et de la sortie g sont li ́ees par
ˆg(λ) = H(λ) ˆf(λ) .
Un filtre id ́eal est un filtre qui ́elimine totalement les bandes de fr ́equence ind ́esirables sans transition et sans d ́ephasage dans les bandes conserv ́ees.
Selon la bande rejet ́ee, on rencontre 4 grandes cat ́egories de filtres.
La r ́egion o`u les fr ́equences sont coup ́ees est en gris ci-dessus. Par exemple le passe-bas id ́eal est le filtre qui ne modifie par les fr ́equences λ telles que λ ≤ λc (fr ́equence de coupure) et supprime les autres. D’o `u
H(λ) =
{ 1 si |λ| ≤ λc
0 sinon
On reconnaıt h ∈ L2(R) telle que h ˆ= H. C’est
h(t) = sin 2πλπt ct
.
Si sait on que se limite g est continue, `a des signaux born ́ee d’entr ́ee et nulle d’ ́energie `a l’infini. finie, h ∗ f f, est h dans et H sont L2(R) dans comme L2(R) f.
ˆet g = h∗f. On Toutefois un tel filtre n’est pas r ́ealisable. On se contente de filtres passe-bas approch ́es (voir ci-dessous un dispositif physique de r ́ealisation d’un tel filtre) qui att ́enuent les hautes fr ́equences au lieu de les supprimer (voir exercices) et qui entraınent un d ́ephasage :
3.2. FILTRAGE 2D : CONVOLUTION /M ́EDIAN 33
3.2 Filtrage 2D : convolution /m ́edian
Le c ́el`ebre format bitmap, qui tire son nom de l’anglais « bitmap » pour « carte de bits » montre qu’une image est avant tout un domaine spatial sur lequel on peut se promener avec la souris de l’ordinateur : les distances en pixels dans l’image I sont d`es lors li ́ees aux distances r ́eelles en m`etres dans la sc`ene r ́eelle S.
La fr ́equence spatiale est un concept d ́elicat qui d ́ecoule du fait que les images appar- tiennent au domaine spatial. Pour commencer on peut rappeler que la fr ́equence est une gran- deur qui caract ́erise le nombre de ph ́enom`enes qui se d ́eroulent au cours d’un temps donn ́e : en voiture le long d’une route vous voyez 2 bandes blanches PAR seconde : c’est une fr ́equence temporelle. Il est ensuite facile de comprendre que ce concept de fr ́equence «temporelle» peut aussi se traduire en disant qu’il y a 200 bandes blanches PAR kilom`etre : c’est une fr ́equence spatiale.
Dans une image, les d ́etails se r ́ep`etent fr ́equemment sur un petit nombre de pixels, on dit qu’ils ont une fr ́equence ́elev ́ee : c’est le cas pour les bords et les contours dans une image.
Au contraire, les fr ́equences basses correspondent `a des variations qui se r ́ep`etent peu car, dilu ́ees sur de grandes parties de l’image, par exemple des variations de fond de ciel.
Nous verrons dans la suite que la plupart des filtres agissent s ́electivement sur ces fr ́equences pour les s ́electionner, en vue de les amplifier ou de les r ́eduire.
34 CHAPITRE 3. FILTRAGE
3.2.1 Filtrage spatial
Filtres de convolution
Le filtrage spatial est essentiellement une op ́eration de convolution (2D). Si f est l’image `a filtrer (ou `a rehausser) et g le filtre spatial (ou PSF ou masque)ona:
f(x, y) ∗ g(x, y) = F−1
F(f(x, y)) · F(g(x, } {{ y))} }
G(u,v)
.
G est la fonction de transfert du filtre. On peut distinguer trois types de filtrage :
Une image num ́erique ́etant essentiellement discr`ete (pixels et niveaux de gris) nous allons pr ́esenter les filtres dans le cas discret. Dans tout ce qui suit x et y sont des entiers (coordonn ́ees des pixels) et f est `a valeurs enti`eres (dans {0,··· ,255}).
3.2. FILTRAGE 2D : CONVOLUTION /M ́EDIAN 35
On ne fait pas en g ́en ́eral une convolution globale mais une transformation bas ́ee sur le voisinage d’un point (x, y) :
Le noyau de convolution du filtre κ est `a support compact inclus dans [x1,x2] × [y1,y2] :
g(x, y)=(f ∗ κ)(x, y) =
∑x2∑y2i=x1 j=y1
f(x − i, y − j)κ(i, j) .
G ́en ́eralement le filtre est de dimension d impaire et est sym ́etrique. Dans ce cas
[x1,x2]=[y1,y2]=[−d/2, d/2] ,
(f ∗ κ)(x, y) =
(d−1)/2 ∑
i=−(d−1)/2
(d−1)/2 ∑ j=−(d−1)/2f(x + i, y + j)κ(i, j) . (3.2.1)
Filtre(i,j)
w1 w2 w3 ← y − 1 w4 w5 w6 ← y w7 w8 w9 ← y + 1
↑ ↑ ↑ x − 1 x x + 1
Ici d = 3. On ne filtre pas les bords pour ́eviter des distorsions ; donc κ(0,0) = w5. Sur cet exemple on a pr ́ecis ́ement
g(x, y) = w1f(x − 1,y − 1) + w2f(x, y − 1) + w3f(x + 1,y − 1) +w4f(x − 1,y) + w5f(x, y) + w6f(x + 1,y) +w7f(x − 1,y + 1) + w8f(x, y + 1) + w9f(x + 1,y + 1) .
Afin 1 :∑de conserver la moyenne de l’image f, la somme des ́el ́ements du filtre est normalis ́ee `a
i
wi = 1.
36 CHAPITRE 3. FILTRAGE
Un filtre 2D est dit s ́eparable s’il est possible de d ́ecomposer le noyau de convolution h2D en deux filtres 1D appliqu ́es successivement en horizontal puis en vertical (ou inversement) :
h2D = hV1D ⊗ hH1D , o `u le symbole ⊗ d ́esigne le produit de convolution. On peut alors traiter s ́epar ́ement les lignes et les colonnes de l’image. Pour qu’un filtre 2D soit s ́eparable il faut et il suffit que les coefficients de ses lignes et de ses colonnes soient proportionnels.
Exemple 3.2.1 (Filtres s ́eparables )
a b c ⊗
αβγ
aα bα cα aβ bβ cβ aγ bγ cγ
Exemple 3.2.2 (Filtre de moyenne passe -bas )
19 ·
=
1 1 1 1 1 1 1 1 1 1 1 1 1 1
125 ·
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Filtre 3 x 3 Filtre 5 x 5
Ce sont des filtres séparables.
3.2. FILTRAGE 2D : CONVOLUTION /M ́EDIAN 37
Exemple 3.2.3 (Filtre gaussien)
1 2 1 2 4 2 1 2 1
116 ·
Id ́ealement, on devrait pr ́evoir un filtre de taille (6σ+1)×(6σ+1). En g ́en ́eral un filtre gaussien avec σ < 1 est utilis ́e pour r ́eduire le bruit, et si σ > 1 c’est dans le but de fabriquer une image qu’on va utiliser pour faire un « masque flou » personnalis ́e. Il faut noter que plus σ est grand, plus le flou appliqu ́e `a l’image sera marqu ́e.
Exemple 3.2.4 (Filtre binˆomial)
Les coefficients de ce filtre sont obtenus par le binˆome de Newton. Un filtre 1D binˆomial d’ordre 4
est un filtre par v v :
séparable donné par le vecteur v = 161(1 4 6 4 1). Un filtre 2D binˆomial d’ordre 4 est donné
1256 ·
1 4 6 4 1 4 16 24 16 4 6 24 36 24 6 4 16 24 16 4 1 4 6 4 1
Filtres m ́edians
Ce ne sont pas des filtres de convolution, ni des filtres lin ́eaires.
g(x, y) = m ́edian{f(n, m) | (n, m) ∈ S(x, y) } ,
o`u S(x, y) est un voisinage de (x, y).
38 CHAPITRE 3. FILTRAGE
Exemple :
30 10 20 10 250 25 20 25 30
bruit ↓ →
10 10 20 20 25 25 30 30 250
↑ m ́ediane
On remplace la valeur du pixel par la valeur m ́ediane ou la valeur moyenne. Ce filtre est utile pour contrer l’effet « Poivre et Sel » (P& S) c’est-`a-dire des faux « 0» et « 255 » dans l’image.
Image bruit ́ee « Poivre et Sel »
3.2. FILTRAGE 2D : CONVOLUTION /M ́EDIAN 39
Filtre de moyenne : rayon 3 Filtre de moyenne : rayon 5
Filtre m ́edian : rayon 3 Filtre m ́edian : rayon 5
Filtre de moyenne : rayon 7
40 CHAPITRE 3. FILTRAGE
Filtres passe-haut
L’image obtenue par un filtre passe-haut correspond en g ́en ́eral `a ce qui « reste » apr`es un filtrage passe-bas.
Filtre m ́edian : rayon 7
Si le bruit P& S est sup ́erieur `a la moiti ́e de la dimension du filtre, le filtrage est inefficace.
3.2. FILTRAGE 2D : CONVOLUTION /MEDIAN ́41
3.2.2 Filtrage fr ́equentiel
Transformation de Fourier 2D
Avant d’envisager le traitement num ́erique des signaux il est n ́ecessaire de donner l’in- terpr ́etation des signaux bidimensionnels dans le domaine des fr ́equences.
La repr ́esentation fr ́equentielle des signaux 2D est l’extension directe de celle des signaux monodimensionnels. La transform ́ee F(u, v) de Fourier d’un signal 2D f(x, y) est
F(u, v) =
∫∫R2 f(x, y)e−2iπ(xu+yv) dx dy . (3.2.2)
La reconstitution du signal spatial se fait par transformation inverse :
f(x, y) =
∫∫R2 F(u, v)e2iπ(xu+yv) du dv . (3.2.3)
La transform ́ee de Fourier est une fonction complexe, qui a pour chaque composante un module et une phase.
Elle poss`ede les mˆemes propri ́et ́es que la transformation de Fourier 1D (lin ́earit ́e, d ́ecalage, d ́erivation, convolution, extension `a L2(R2).
On a vu qu’il est ́equivalent de multiplier la transform ́ee de Fourier par une fonction de transfert et d’agir sur la fonction dans la domaine spatial par convolution. Il est parfois plus facile de mettre en œuvre le filtrage en agissant dans le domaine fr ́equentiel (apr`es transform ́ee de Fourier ).
Filtre passe-bas
Nous en avons d ́ej`a parl ́e en 1D du filtre passe-bas id ́eal. Ici on d ́efinit une fr ́equence de coupure δc au dessus de laquelle les fr ́equences sont annul ́ees (filtre id ́eal) . La fonction de transfert est alors
H(λ,μ) =
{ 1 si √λ2 + μ2 ≤ δc
0 sinon
Le cr ́eneau centr ́e H admet une transform ́ee de Fourier inverse qui est le sinus cardnal qui pr ́esente d’autant plus d’ondulations que la fr ́equence de coupure est petite. Cela entraıne un flou qui sera d’autant plus r ́eduit que δc est grand.
42 CHAPITRE 3. FILTRAGE
Application d’un cr ́eneau « id ́eal » (δc ≃ 15% de la taille de l’image) : on voit clairement les ondulations
On a vu que le filtre passe-bas id ́eal n’est pas r ́ealisable et on fait donc une approximation de la fonction H pr ́ec ́edente qui aura pour effet, non pas de couper les hautes fr ́equences mais de les att ́enuer fortement. Le filtre suivant est le filtre passe-bas de Butterworth La fonction de
3.2. FILTRAGE 2D : CONVOLUTION /MEDIAN ́43
transfert est alors
H(λ,μ) = 1
1 +
(√λ2+μ2
δc
)2n
o`u δc est encore la fr ́equence de coupure.
En image un filtre passe-bas att ́enue les hautes fr ́equences : le r ́esultat obtenu apr`es un tel filtrage est un adoucissement des d ́etails et une r ́eduction du bruit granuleux.
Filtres passe-haut
Le filtre passe-haut id ́eal est obtenu de mani`ere sym ́etrique au passe- bas par
H(λ,μ) =
{ 1 si √λ2 + μ2 ≥ δc
0 sinon
Le filtre passe-haut de Butterworth est donn ́e par
H(λ,μ) = 1
1 +
(
√λδc 2+μ2)2n
44 CHAPITRE 3. FILTRAGE
Un filtre passe haut favorise les hautes fr ́equences spatiales, comme les d ́etails, et de ce fait, il am ́eliore le contraste. Toutefois, il produit des effets secondaires
– Augmentation du bruit : dans les images avec un rapport Signal/ Bruit faible, le filtre
augmente le bruit granuleux dans l’image. – Effet de bord : il est possible que sur les bords de l’image apparaisse un cadre. Mais cet
effet est souvent n ́egligeable et peut s’ ́eliminer en tronquant les bords de l’image.
Filtrage passe-bas et passe-haut avec un filtre de Butterworth (n = 4 et δc ≃ 0.15∗ taille de l’image)
3.3. FILTRAGE DIFFERENTIEL ́45
Filtres passe-bande
Ils permettent de ne garder que les fr ́equences comprises dans un certain intervalle :
H(λ,μ) =
{ 1 0 si sinon δc − 2 ε≤ √λ2 + μ2 ≤ δc + 2 εε est la largeur de bande et δc la fr ́equence de coupure.
3.3 Filtrage diff ́erentiel
Dans les mod`eles diff ́erentiels, on consid`ere l’image comme une fonction continue f : I × I → [0,255] dont on ́etudie le comportement local `a l’aide de ses d ́eriv ́ees. Une telle ́etude n’a de sens que si la fonction f est assez r ́eguli`ere. Ce n’est pas toujours le cas ! ! une image noir et blanc sera discontinue (en fait continue par morceaux) les zones de discontinuit ́e ́etant par essence les contours.
Au premier ordre on peut calculer en chaque point M(x, y), le gradient de l’image :
∇f(x, y)=(∂f∂x, ∂f∂y) .
Grˆace au plongement dans l’espace continu, un grand nombre d’op ́erations d’analyse peuvent s’exprimer en termes d’ ́equations aux d ́eriv ́ees partielles. Ceci permet de donner un fondement math ́ematique satisfaisant aux traitements et aussi de fournir des m ́ethodes pour les calculer, par des sch ́emas num ́eriques de r ́esolution.
Les filtres diff ́erentiels permettent de mettre en ́evidence certaines variations spatiales de l’image . Ils sont utilis ́es comme traitements de base dans de nombreuses op ́erations comme le r ́ehaussement de contraste ou la d ́etection de contours.
3.3.1 Calcul par convolution
En pratique, il faut approcher les gradients pour travailler avec des gradients discrets. Les approximations les plus simples des d ́eriv ́ees directionnelles se font par diff ́erences finies cal- cul ́ees par convolution avec des noyaux tr`es simples : par exemple, l’approximation de ∂f∂x se
46 CHAPITRE 3. FILTRAGE
fait par convolution avec [0 − 1 1]. En effet, dans ce cas, la formule g ́en ́erale de convolution discr`ete (3.2.1) donne (avec dx = 3 et dy = 0 )
g(x, y) =
∑1i=−1∑j=0
f(x + i, y + j)κ(i, j) .
κ(i, j)
0 −1 1 ← y ↑ ↑ ↑ x − 1 x x + 1
et
g(x, y) = −f(x, y) + f(x + 1,y) ≃ ∂f∂x(x, y)
De mˆeme l’approximation de ∂f∂y se fait par convolution avec
0−11
:
κ(i, j)
0 ← y − 1 −1 ← y
1 ← y + 1 ↑x
et
g(x, y) = −f(x, y) + f(x, y + 1) ≃ ∂f∂y (x, y).
On utilise plus souvent [−1 0 1] et
−101
qui produisent des fronti`eres plus ́epaisses mais qui
sont bien centr ́ees. Ces op ́erations sont tr`es sensibles au bruit et on les combine g ́en ́eralement avec un filtre lisseur dans la direction orthogonale `a celle de d ́erivation, par exemple par le noyau suivant (ou sa transpos ́ee) : [1 2 1]. Le calcul des d ́eriv ́ees directionnelles en x et y revient finalement `a la convolution avec les noyaux suivants :
∂f∂x(x, y) ≃ (f ∗ hx)(x, y) et ∂f∂y (x, y) ≃ (f ∗ hy)(x, y)
avec
hx =
−1 0 1 −2 0 2 −1 0 1
et hy =
−1 −1 −1
0 0 0 1 2 1
Ce sont les masques de Sobel.
3.3. FILTRAGE DIFF ́ERENTIEL 47
De la mˆeme fac ̧on , l’approximation par diff ́erences finies la plus simple de la d ́eriv ́ee se-
conde est la convolution par le noyau [1 − 2 1] pour l’approximation de ∂∂x2f
2 et
1−21
pour
l’approximation de ∂∂y2f
2 . Le laplacien
∆f = ∂∂x2f
2 + ∂∂y2f 2
peut donc ˆetre approch ́e par l’un des deux op ́erateurs lin ́eaires suivant
0 1 0 1 −4 1 0 1 0
ou
1 1 1 1 −8 1 1 1 1
48 CHAPITRE 3. FILTRAGE
Gradient en x Gradient en y D ́etection de contours par une convolution de type Sobel
Laplacien 4-connexe Laplacien 8-connexe
Norme du gradient
3.3. FILTRAGE DIFFERENTIEL ́49
TYPES DE MASQUE GRADIENTS PARTIELS AMPLITUDE DIRECTION
Masques de Roberts -1 0
G1, G2 Substitution du pixel A = √G21 + G22 θ = π4 0 1
0 -1 1 0 sup ́erieur gauche + arctan(GG21)
Masques de Sobel
1 0 -1 2 0 -2 1 0 -1
√G2x + G2y θ = arctan(GGxy)
Masques de Prewitt
1 0 -1 1 0 -1 1 0 -1
1 2 1 0 0 0 -1 -2 -1
Gx, Gy A =
1 1 1 √0 0 0 -1 -1 -1
G2x + G2y θ = arctan(GGxy)
Masques de Kirsh
Direction 5 5 5 -3 0 -3 -3 -3 -3
Gx, Gy A =
Gi pour maximum des |Gi] correspondant
+ les 7 autres masques obtenus par i de1`a8 au Gi s ́electionn ́e permutation circulaire des coefficients
Masques de Robinson
1 1 1 1 -2 1 -1 -1 -1
Gi pour maximum des |Gi] Idem
+ les 7 autres masques obtenus par i de1`a8 permutation circulaire des coefficients
Laplacien discret Laplacien de Robinson
1 1 1 1 -8 1 1 1 1
1 -2 1 -2 4 -2 1 -2 1
50 CHAPITRE 3. FILTRAGE
3.3.2 Equation de la chaleur
Consid ́erons un filtrage gaussien dans le cadre continu. On sait que si l’image de d ́epart est une fonction uo d ́efinie sur R2 (mais L∞ `a support compact), l’image filtr ́ee est la convol ́ee de uo avec un noyau gaussien
Gσ(x) = Gσ(x1,x2) = 2πσ1
2 exp(−x21 2σ+ 2
x22
)
= 2πσ1
2 exp(− 2σx 2 2
)
.
On pose u(t, x)=(h(t,·) ∗ uo)(x) o`u h(t, x) = G√2t(x) = 4πt 1exp(− x 2
4t
)
. Comme h(t,·) ∈
C∞(R2) a ses d ́eriv ́ees born ́ees et uo ∈ L1(R2), la convol ́ee u(t,·) est aussi C∞ et on peut calculer ∆u :
∀t > 0,∀x ∈ R2 ∆u(t, x) = ∂∂x2u
21(t, x) + ∂∂x2u
22(t, x)=(∆h(t,·) ∗ uo)(x) .
Un rapide calcul montre que
∆h(t, x) =
(− 4πt1
2 + 16πtx 2
3)exp(− x 2
4t
)
=
(−1t + x 2 4t2
)h(t, x) ,
et on obtient
∆u(t, x) =
(−1t + x 2 4t2
)u(t, x) .
D’autre part, pour t > 0 on peut d ́eriver directement u par rapport `a t :
∂u∂t (t, x) =
∫∫R
∂h∂t (t, y)uo(x − y)dy
et on obtient finalement
∂u∂t (t, x) − ∆u(t, x)=0 sur ]0,t[×R2 .
D’autre part, avec
u(t, x) =
∫∫1R
4πt exp(− y 4t
2
)uo(x − y)dy
on obtient
t→0lim u(t, x) =< δx,uo >= uo(x) ,
car la famille de Gaussiennes converge au sens des distributions vers la mesure de Dirac.
3.3. FILTRAGE DIFF ́ERENTIEL 51
La fonction « filtr ́ee » u v ́erifie l’ ́equation aux d ́eriv ́ees partielles suivante ( ́equation de la chaleur) : { u(0,x) ∂u∂t (t, x) = − u∆u(t, o(x) x)=0 dans ∀x ∈ Ω ]0,T[×Ω
(3.3.4)
o`u Ω est le « cadre » de l’image, i.e. l’ouvert de R2 o `u la fonction u est d ́efinie . On peut alors imposer soit des conditions aux limites au bord de Ω (niveau de gris fix ́e) soit des conditions aux limites p ́eriodiques en p ́eriodisant la fonction u si le cadre est rectangle par exemple.
On peut alors utiliser un sch ́ema aux diff ́erences finies pour calculer la solution de l’EDP. Suivant le temps d’ ́evolution, on obtient une version plus ou moins liss ́ee de l’image de d ́epart.
Mod`ele de Peronna-Malik
Sur une image liss ́ee, on peut plus facilement essayer de d ́etecter les bords (ou contours). On peut utiliser le d ́etecteur de Hildrett-Marr : on cherche les z ́eros du Laplacien d’une image u. Si en un point x, ∆u change de signe et si ∇u est non nul, on consid`ere alors que l’image u poss`ede un bord en x.
Pour am ́eliorer les r ́esultats obtenus par l’EDP de la chaleur, Peronna et Malik ont propos ́e de modifier l’ ́equation en y int ́egrant le processus de d ́etection des bords :
u(0,x) ∂u∂u∂n ∂t (t, = x) = div (c(|∇u|)∇u)(t, 0 sur ]0,T[×∂Ω
= uo(x) ∀x ∈ Ω
- x) dans ]0,T[×Ω
(3.3.5)
o`u c est une fonction d ́ecroissante de R+ dans R+.
52 CHAPITRE 3. FILTRAGE
Si c = 1, on retrouve l’ ́equation de la chaleur. On impose souvent lim t→+∞c(t)=0 et c(0) = 1. Ainsi, dans les r ́egions de faible gradient, l’ ́equation agit essentiellement comme l’EDP de la chaleur, et dans les r ́egions de fort gradient, la r ́egularisation est stopp ́ee ce qui permet de pr ́eserver les bords.
Filtrage (EDP Chaleur) Filtrage de Peronna-Malik
3.3.3 Mise en œuvre num ́erique
La mise en œuvre num ́erique se fait avec une discr ́etisation en diff ́erences finies. La condi- tion de Neumann est assur ́ee grˆace `a une r ́eflexion de l’image par rapport `a ses bords. En image on consid`ere souvent que la taille est donn ́ee par le nombre de pixels de sorte que le pas de discr ́etisation est h = 1. On peut discr ́etiser le gradient de diff ́erentes mani`eres (centr ́ee, `a droite, `a gauche )
δxui,j = ui+1,j − ui−1,j
2 , δyui,j = ui,j+1 − ui,j−1
2 , (3.3.6) δ+x ui,j = ui+1,j − ui,j, δ+y ui,j = ui,j+1 − ui,j,
Image originale Bords donn ́es par Hildrett-Marr
(apr`es filtrage)
3.3. FILTRAGE DIFFERENTIEL ́53
δx −ui,j = ui,j − ui−1,j, δy −ui,j = ui,j − ui,j−1. La norme du gradient peut se calculer par un sch ́ema ENO (Essentially Non Oscillatory).Deux approximations possibles de |∇u| sont
∇+ui,j =
√max(δx −ui,j,0)2 + min(δx +ui,j,0)2 + max(δy −ui,j,0)2 + min(δy +ui,j0)2 , (3.3.7)
ou
∇−ui,j =
√max(δx +ui,j,0)2 + min(δx −ui,j,0)2 + max(δy +ui,j,0)2 + min(δy −ui,j0)2 . (3.3.8)
Si l’op ́erateur gradient est discr ́etis ́e par diff ́erences finies `a droite, alors une discr ́etisation pos- sible de la divergence est donn ́ee par
(div p)i,j =
p1i,j − p1i−1,j si 1 <i<N
p1i,j si i = 1 −p1i−1,j si i = N
p2i,j − p2i,j−1 si 1 <j<N +
p2i,j si j = 1
(3.3.9)
−p2i,j−1 si j = N
avec p = (p1,p2) et N est la taille de l’image (carr ́ee).
54 CHAPITRE 3. FILTRAGE
Chapitre 4
Quelques mod`eles de restauration d’image
Etant donn ́ee une image originale , on suppose qu’elle a ́et ́e d ́egrad ́ee par un bruit additif v, et ́eventuellement par un op ́erateur R . Un tel op ́erateur est souvent mod ́elis ́e par un produit de convolution , et n’est pas n ́ecessairement inversible (et mˆeme lorsqu’il est inversible, son inverse est souvent num ́eriquement difficile `a calculer). A partir de l’image observ ́ee f = Ru+v (qui est donc une version d ́egrad ́ee de l’image originale u), on cherche `a reconstruire u. Si on suppose que le bruit additif v est gaussien , la m ́ethode du Maximum de vraisemblance nous conduit `a chercher u comme solution du probl`eme de minimisation
infu f − Ru 22 , o`u · 2 d ́esigne la norme dans L2 . Il s’agit d’un probl`eme inverse mal pos ́e. Pour le r ́esoudre num ́eriquement, on est amen ́e `a introduire un terme de r ́egularisation, et `a consid ́erer le probl`eme
infu } f − {{ Ru 22 }
ajustement aux donn ́ees+ L(u) }{{}
R ́egularisation
.
Dans toute la suite, nous ne consid`ererons que le cas o `u est R est l’op ́erateur identit ́e (Ru = u). Commenc ̧ons par un proc ́ed ́e de r ́egularisation classique : celui de Tychonov.
4.1 R ́egularisation de Tychonov
C’est un proc ́ed ́e de r ́egularisation tr`es classique mais trop sommaire dans le cadre du trai- tement d’image. Nous le pr ́esentons sur un exemple.
Soit V = Ho1(Ω) et H = L2(Ω) : on consid`ere le probl`eme de minimisation originel (ajuste- ment aux donn ́ees) :
(P) min u∈V u − ud 2H , o`u ud est l’image observ ́ee (donn ́ees) et le probl`eme r ́egularis ́e suivant : pour tout α > 0
(Pα) min u∈V u − ud 2H + α ∇u 2H .
55
56 CHAPITRE 4. QUELQUES MODELES `DE RESTAURATION D’IMAGE
Non seulement on veut ajuster u `a la donn ́ee ud, mais on impose ́egalement que le gradient soit «assez petit» (cela d ́epend du param`etre α). Une image dont le gradient est petit est «liss ́ee», estomp ́ee. Les bords sont ́erod ́es et la restauration donnera une image flout ́ee. Il est facile de coercive sur V :
voir sur l’exemple suivant que la fonctionnnelle J(u) = u − ud 2H n’est pas
Ω =]0,1[, un(x) = xn, ud = 0 .
On voit que
un 2 = √12n , u n 2 = √2n n − 1 . On a donc
n→+∞ lim un V = +∞ et n→+∞lim J(un)=0 . Il n’est donc mˆeme pas clair (a priori) que le probl`eme (P) ait une solution.
Proposition 4.1.1 Supposons que (P) admet au moins une solution ̃u. Le probl`eme (Pα) admet une solution unique uα. De la famille (uα) on peut extraire une sous-suite qui converge (faiblement ) dans V vers une solution u∗ de (P) lorsque α → 0.
Démonstration – Le probl`eme (Pα) admet une solution unique uα car la fonctionnelle
Jα = u − ud 2H + α ∇u 2H ,
est coercive et strictement convexe (c’est en gros la norme de V `a une partie affine pr`es) . Mon- trons maintenant que la famille (uα) est born ́ee dans V . On pourra ainsi en extraire une sous- suite faiblement convergente dans V vers u∗ ∈ V .
∀u ∈ V Jα(uα) ≤ Jα(u) .
En particulier pour ̃u :
J( ̃u) } ≤ {{ J(uα) } ̃u est solution de (P)
≤ J} α(uα) = J(uα) + {{ α ∇uα 2H ≤ Jα( ̃u) }
uα est solution de (Pα)
= J( ̃u) + α ∇ ̃u 2H . (4.1.1)
Par cons ́equent, pour α ≤ αo, Jα(uα) est born ́e ind ́ependamment de α. Ceci entraıne que la famille (uα)α≤αo est born ́ee dans H. D’autre part, avec (4.1.1) on a aussi
α ∇uα 2H ≤ J( ̃u) + α ∇ ̃u 2H − J(uα) ≤ J( ̃u) + α ∇ ̃u 2H − J( ̃u) = α ∇ ̃u 2H ;
par cons ́equent la famille (uα)α≤αo qui converge (faiblement ) dans V est born ́ee dans V . On peut donc en extraire une sous-suite
vers u∗ . D’autre part l’ ́equation (4.1.1) montre que
α→0lim Jα(uα) = J( ̃u) = inf(P) .
Par semi-continuit ́e inf ́erieure de J il vient
J(u∗) ≤ lim α→0 inf J(uα) = lim α→0 inf Jα(uα) ≤ inf(P) .
4.2. L’ESPACE DES FONCTIONS A `VARIATION BORNEE ́BV (Ω) 57
Par cons ́equent u∗ est une solution de (P). 2 Cherchons maintenant le moyen de calculer uα. Comme la fonctionnelle est strictement convexe, une condition n ́ecessaire et suffisante d’optimalit ́e est
J α(uα)=0 .
Un calcul assez standard montre que
∀u ∈ V J α(uα) · u =
∫Ω(uα − ud)(x)u(x)dx +
∫Ω ∇uα(x)∇u(x)dx
=
∫Ω(uα − ud − ∆uα)(x)u(x)dx .
Par cons ́equent l’ ́equation d’Euler qui fournit la solution uα est la suivante :
uα − ud − ∆uα = 0 , uα ∈ Ho1(Ω) .
Comme pour la m ́ethode des contours actifs d ́ecrite dans le chapitre suivant, on peut se conten- ter d’approcher la solution uα en ́ecrivant la formulation dynamique :
∂u∂t − ∆u + u = ud .
Le terme de r ́egularisation classique L(u) = ∇u adapt ́e au probl`eme de restauration d’images : l’image 22( r ́egularisation restaur ́ee u de Tychonov) n’est pas est alors beaucoup trop liss ́ee (en consid ́erer particulier, les bords sont ́erod ́es). Une la variation totale , c’est `a dire `a prendre approche L(u) = ∫ beaucoup |Du|.
plus efficace consiste `a
Cela conduit `a une minimisation de fonctionnelle dans un espace de Banach particulier, mais bien adapt ́e au probl`eme : l’espace des fonctions `a variation born ́ee que nous commenc ̧ons par pr ́esenter.
4.2 L’espace des fonctions `a variation born ́ee BV (Ω)
4.2.1 G ́en ́eralit ́es
Dans ce qui suit Ω un ouvert born ́e de R2 de fronti`ere Lipschitz fonctions C1 `a support compact dans Ω `a valeurs dans R2.
et Cc1(Ω,R2) est l’espace des
D ́efinition 4.2.1 Une fonction f de L1(Ω) (`a valeurs dans R) est `a variation born ́ee dans Ω si J(f) < +∞ o`u
J(f) = sup{∫Ω f(x)div φ(x)dx | φ ∈ Cc1(Ω,R2) , φ ∞ ≤ 1
}. (4.2.2)
On note
BV (Ω) = {f ∈ L1(Ω) | J(f) < +∞}
l’espace de telles fonctions.
58 CHAPITRE 4. QUELQUES MODELES `DE RESTAURATION D’IMAGE
Remarque 4.2.1 On rappelle que si φ = (φ1,φ2) ∈ Cc1(Ω,R2) alors
div φ(x) = ∂φ∂x11
(x) + ∂φ∂x22
(x) .
Donc, par intégration par parties
∫Ω f(x)div φ(x)dx =
(f(x) ∂φ∂x1 1
(x) + ∂φ∂x2 2 ∫)
(x)dx = −Ω
∫Ω
( ∂x∂f
1(x)φ1(x) + ∂x∂f
2)
(x)φ2(x)dx = −∫Ω ∇f(x) · φ(x) dx , o`u · désigne le produit scalaire de R2.
D ́efinition 4.2.2 (P ́erim`etre) Un ensemble E mesurable (pour la mesure de Lebesgue) dans R2 est de périm`etre (ou de longueur) fini si sa fonction indicatrice χE est dans BV (Ω).
Commenc ̧ons par donner une propri ́et ́e structurelle des fonctions BV.
Th ́eor`eme 4.2.1 Soit f ∈ BV (Ω). Alors il existe une mesure de Radon positive μ sur Ω et une fonction μ-mesurable σ : Ω → R telle que (i) |σ(x)| = 1, μ p.p. , et
(ii)
∫Ω f(x)div φ(x)dx = −∫Ω φ(x)σ(x)dμ pour toute fonction φ ∈ Cc1(Ω,R2) La relation (ii) est une formule d’int ́egration par parties « faible ». Ce th ́eor`eme indique que la d ́eriv ́ee faible (au sens des distributions ) d’une fonction BV est une mesure de Radon.
On rappelle qu’une mesure de Radon est une mesure finie sur tout compact et que grˆace au th ́eor`eme de Riesz, toute forme L lin ́eaire continue sur Cco(Ω) (fonctions continues `a support compact) est de la forme
L(f) =
∫Ω f(x) dμ , o`u μ est une (unique) mesure de Radon associ ́ee `a L. Plus pr ́ecis ́ement
Th ́eor`eme 4.2.2 ([Rudin] p. 126, [EG] p. 49 ) A toute forme linéaire bornée Φ sur Cco(Ω,R2), c’est- `a-dire
∀K compact de Ω sup{Φ(φ) | φ ∈ Cco(Ω,R2) , φ ∞ ≤ 1 , suppφ ⊂ K } < +∞ , il correspond une unique mesure de Radon positive μ sur Ω et une fonction μ-mesurable σ (fonction « signe ») telle que (i) |σ(x)| = 1, μ p.p. , et
(ii) Φ(φ) =
∫Ω φ(x)σ(x)dμ pour toute fonction φ ∈ Cco(Ω,R2) . (iii) De plus μ est la mesure de variation et vérifie
μ(Ω) = sup{ Φ(φ) | φ ∈ Cco(Ω,R2) , φ ∞ ≤ 1 , suppφ ⊂ Ω } . (4.2.3)
4.2. L’ESPACE DES FONCTIONS A `VARIATION BORNEE ́BV (Ω) 59
Démonstration du Théor`eme 4.2.1 – Soit f un ́el ́ement de BV (Ω). On consid`ere la forme lin ́eaire L d ́efinie sur Cc1(Ω,R2) par
L(φ) =
∫Ω f(x)div φ(x)dx . Comme f ∈ BV (Ω),
sup{L(φ) | φ ∈ Cc1(Ω,R2) , φ ∞ ≤ 1 } := CL < +∞
pour toute fonction φ ∈ Cc1(Ω,R2). Par cons ́equent
∀φ ∈ Cc1(Ω,R2) L(φ) ≤ CL φ ∞ . (4.2.4)
Soit K un compact de Ω. Pour toute fonction φ trouver (par densit ́e) une suite de fonctions φk ∈ ∈ CCcc1o(Ω,R(Ω,R22) ) qui `a support converge compact uniform ́ement dans K, on vers peut φ. Posons alors L(φ) ̄= k→+∞lim L(φk) . Grˆace ainsi ́etendre `a (4.2.4) L cette par densit ́e limite existe en une et forme elle est lin ́eaire ind ́ependante L ̄sur Ccode (Ω,Rla suite 2) telle (φque
- k) choisie. On peut donc
sup{ L(φ) ̄| φ ∈ Cco(Ω,R2) , φ ∞ ≤ 1 , suppφ ⊂ K } < +∞ .
On conclut alors grˆace au th ́eor`eme de Riesz (que nous avons rappel ́e). 2 D’apr`es la propri ́et ́e (4.2.3), J(u) = μ(Ω) ≥ 0 : c’est la variation totale de f.
Proposition 4.2.1 L’applicationBV (Ω) → R+
u ↦→ u BV (Ω) = u L1 + J(u) .
est une norme.
Démonstration – La d ́emonstration est facile et laiss ́ee en exercice. 2 On munit d ́esormais l’espace BV (Ω) de cette norme.
Exemple 4.2.1 Supposons que
f ∈ W 1,1(Ω) = { f ∈ L1(Ω) | Df ∈ L1(Ω) } ,
o`u Df est la dérivée de f au sens des distributions. Soit φ ∈ Cc1(Ω,R2) telle que φ ∞ ≤ 1. Alors
∫Ω f div φdx = −∫Ω Df · φdx ≤ φ ∞
∫Ω |Df|dx ≤ Df L1 < +∞ .
Donc f ∈ BV (Ω). De plus
J(f) = sup{−∫Ω Df · φdx | φ ∞ ≤ 1 } = Df L1 ,
60 CHAPITRE 4. QUELQUES MODELES `DE RESTAURATION D’IMAGE
et
σ =
Df |Df| si Df = 0 , 0 sinon.
Donc
W 1,1(Ω) ⊂ BV (Ω) .
En particulier, comme Ω est borné∀1 ≤ p ≤ +∞ W 1,p(Ω) ⊂ BV (Ω) .
Toute fonction d’un espace de Sobolev est `a variation bornée.
Remarque 4.2.2 D’apr`es le théor`eme de Radon-Nikodym de décomposition des mesures, pour toute fonction u ∈ BV (Ω), nous avons la décomposition suivante de Du (dérivée au sens des distributions) :
Du = ∇u dx + Dsu ,
o`u ∇u dx est la partie absolument continue de Du par rapport `a la mesure de Lebesgue et Dsu est la partie singuli`ere.
4.2.2 Approximation et compacit ́e
Th ́eor`eme 4.2.3 (Semi-continuit ́e inf ́erieure de la variation totale) L’application u ↦→ J(u) de BV (Ω) dans R+ est semi-continue inférieurement (sci) pour la topologie séquentielle de L1(Ω). Plus précisément, si (uk) est une suite de fonctions de BV (Ω) qui converge vers u fortement dans L1(Ω) alors
J(u) ≤ lim k→+∞ inf J(uk) .
Démonstration – Soit φ ∈ Cc1(Ω,R2) telle que φ ∞ ≤ 1. Alors
∫Ω u(x)div φ(x)dx = k→+∞lim ∫Ω uk(x)div φ(x)dx .
Donc, ∀ε > 0,∃k(φ,ε) tel que
∀k ≥ k(φ,ε)
∫Ω u(x)div φ(x)dx − ε ≤
∫Ω uk(x)div φ(x)dx ≤
∫Ω u(x)div φ(x)dx + ε .
Comme ∫Ω uk(x)div φ(x)dx ≤ J(uk) il vient
∀k ≥ k(φ,ε)
∫Ω u(x)div φ(x)dx − ε ≤ J(uk) ,
4.2. L’ESPACE DES FONCTIONS `A VARIATION BORN ́EE BV (Ω) 61
et donc ∫Ω u(x)div φ(x)dx ≤ lim inf k→+∞ J(uk) . Comme c’est le cas pour tout φ, on obtient
J(u) ≤ lim inf k→+∞ J(uk) .
2 Nous admettrons le r ́esultat suivant
Th ́eor`eme 4.2.4 (Approximation r ́eguli`ere) Pour toute fonction u ∈ BV (Ω), il existe une suite de fonctions (uk)k∈N de BV (Ω) ∩ C∞(Ω) telle que (i) uk → u dans L1(Ω) et (ii) J(uk) → J(u) (dans R).
La d ́emonstration est technique et utilise un proc ́ed ́e classique de r ́egularisation par convolu- tion. On peut se r ́ef ́erer `a [EG] p.172. Remarquons que le r ́esultat ci-dessus n’est pas un r ́esultat de densit ́e de BV (Ω) ∩ C∞(Ω) dans BV (Ω) car on n’a pas J(uk − u) → 0.
Th ́eor`eme 4.2.5 L’espace BV (Ω) muni de la norme
u BV (Ω) = u L1 + J(u) ,
est un espace de Banach.
Démonstration – Soit (un)n∈N une suite de Cauchy dans BV (Ω). C’est donc une suite de Cauchy dans L1(Ω) : elle converge donc vers u ∈ L1(Ω). D’autre part, elle est born ́ee dans BV (Ω) (toute suite de Cauchy est born ́ee), donc∃M > 0,∀n J(un) ≤ M .
D’apr`es le th ́eor`eme 4.2.3,
J(u) ≤ lim inf n→∞ J(un) ≤ M < +∞ .
Par cons ́equent u ∈ BV (Ω). Soit ε > 0 et N tel que
∀n, k ≥ N un − uk BV (Ω) ≤ ε .
Donc
∀n, k ≥ N J(un − uk) ≤ ε et avec la semi-continuit ́e inf ́erieure de J en fixant n on obtient
∀n ≥ N J(un − u) ≤ lim inf k→∞ J(un − uk) ≤ ε .
Ceci prouve que J(un − u) → 0. 2 Terminons par un r ́esultat important de compacit ́e que nous admettrons.
62 CHAPITRE 4. QUELQUES MOD`ELES DE RESTAURATION D’IMAGE
Th ́eor`eme 4.2.6 (Compacit ́e) L’espaceBV (Ω)s’injecte dansL1(Ω)de mani`ere compacte. Plus précisément : si (un)n∈N une suite bornée de BV (Ω)sup
n∈N un BV (Ω) < +∞ ,
alors, il existe une sous-suite (unk)k∈N et une fonction u ∈ BV (Ω) telle que unk converge fortement vers u dans L1(Ω).
Démonstration – [EG] p.176. 2 Pour plus de d ́etail sur les fonctions BV, on pourra se reporter `a [EG]. Dans ce qui suit, nous allons pr ́esenter plusieurs mod`eles qui se diff ́erencient par la nature du terme de r ́egularisation. En pratique, on p ́enalise par la norme de tout ou partie de la fonction u de l’image et le cadre fonctionnel est tr`es important.
4.3 Le mod`ele continu de Rudin-Osher-Fatemi
4.3.1 Pr ́esentation
Il faut trouver une alternative `a la r ́egularisation de Tychonov qui est trop violente. La premi`ere id ́ee consiste `a remplacer le terme de r ́egularisation ∇u 2H qui est en r ́ealit ́e une p ́enalisation par une norme V , par une norme moins contraignante et r ́egularisante. Rudin- Osher et Fatemi ont propos ́e un mod`ele o `u l’image est d ́ecompos ́ee en deux parties : ud = u+v o`u v est le bruit et u la partie «r ́eguli`ere». On va donc a priori chercher la solution du probl`eme sous la forme u + v avec u ∈ BV (Ω) et v ∈ L2(Ω) et ne faire porter la r ́egularisation que sur la partie « bruit ». Cela conduit `a :
(PROF) min{ J(u) + 12ε v 22 | u ∈ BV (Ω), v ∈ L2(Ω) , u + v = ud } .
Ici J(u) d ́esigne la variation totale de u et ε > 0.
Th ́eor`eme 4.3.1 Le probl`eme (PROF) admet une solution unique.
Démonstration – Soit un ∈ BV (Ω), vn ∈ L2(Ω) une suite minimisante. Comme vn est born ́ee dans L2(Ω) on peut en extraire une sous-suite (not ́ee de la mˆeme fac ̧on) faiblement convergente vers v∗ dans L2(Ω). Comme la norme de L2(Ω) est convexe, sci, il vient
v∗ 22 ≤ lim inf n→+∞ vn 22 .
De la mˆeme fac ̧on un = ud − vn est born ́ee dans L2(Ω) et donc dans L1(Ω) puisque Ω est born ́e. Comme J(un) est born ́e , il s’ensuit que un est born ́ee dans BV (Ω). Grˆace au r ́esultat de compacit ́e du th ́eor`eme 4.2.6 , cela entraıne que un converge (`a une sous-suite pr`es) fortement dans L1(Ω) vers u∗ ∈ BV (Ω). D’autre part J est sci (th ́eor`eme (4.2.3) , donc
J(u∗) ≤ lim inf n→+∞ J(un) ,
4.4. MODELE `DISCRET 63
et finalement
J(u∗) + 2ε 1v∗ 22 ≤ lim n→+∞ inf J(un) + 2ε 1vn 22 = inf(PROF) . Comme un + vn = ud pour tout n, on a u∗ + v∗ = ud. Par cons ́equent u∗ est une solution du probl`eme (PROF). La fonctionnelle est strictement convexe par rapport au couple (u, v) et la contrainte est affine. On a donc unicit ́e. 2 Nous aurons besoin d’ ́etablir des conditions d’optimalit ́e pour la ou les solutions optimales des mod`eles propos ́es. Toutefois les fonctionnelles consid ́er ́ees (en particulier J) ne sont en g ́en ́eral pas Gˆateaux-diff ́erentiables et nous devons utiliser des notions d’analyse non lisse (voir chapitre 3).
4.3.2 Condition d’optimalit ́e du premier ordre
Le probl`eme (PROF) peut s’ ́ecrire de la mani`ere ( ́equivalente) suivante
u∈BV min (Ω)F(u) := J(u) + 2ε 1u − ud 22. (4.3.5)
La fonctionnelle F est convexe et
̄u solution de (PROF) ⇐⇒ 0 ∈ ∂F( ̄u) .
On aimerait utiliser la relation
0 ∈ ∂F( ̄u) ⇐⇒ ̄u ∈ ∂F∗(0) ,
mais elle n’est vraie que si l’espace est r ́eflexif, ce qui n’est pas le cas de BV (Ω). Toutefois ce sera vrai apr`es discr ́etisation (ce que nous ferons dans la section suivante). On peut utiliser le th ́eor`eme J diff ́erentiable est finie A.3.4 sur sur BV pour L(Ω) 2(Ω). calculer `a valeurs Le calcul ∂F(u). dans du L’application R sous-diff ́erentiel ∪ {+∞} e. u D’autre ↦→ pouvant u − part ud se u 22 faire est ↦→ continue u via − la ud transformation 22 sur est LGˆateaux 2(Ω) est
de Fenchel-Legendre dans le cas o`u l’espace est r ́eflexif (voir Annexe). Nous allons donc faire la discr ́etisation du probl`eme.
4.4 Mod`ele discret
On va maintenant consid ́erer des images discr`etes (ce qui est le cas en pratique). Une image discr`ete est une matrice N × N que nous identifierons `a une vecteur de taille N2 (par exemple en la rangeant ligne par ligne). On note X l’espace euclidien RN×N et Y = X ×X. On munit X du produit scalaire usuel
(u, v)X = ∑
1≤i,j≤N
uijvij ,
et de la norme associ ́ee : · X.
64 CHAPITRE 4. QUELQUES MODELES `DE RESTAURATION D’IMAGE
Nous allons donner une formulation discr`ete de ce qui a ́et ́e fait auparavant et en particulier d ́efinir une variation totale discr`ete que nous noterons de la mˆeme fac ̧on. Pour cela nous intro- duisons une version discr`ete de l’op ́erateur gradient. Si u ∈ X, le gradient ∇u est un vecteur de Y donn ́e par
(∇u)i,j = ((∇u)1i,j,(∇u)2i,j) , avec
(∇u)1i,j =
{ ui+1,j − ui,j si i<N
0 si i = N (4.4.6)
(∇u)2i,j =
{ ui,j+1 − ui,j si j<N
0 si j = N (4.4.7) La variation totale discr`ete est alors donn ́ee par J(u) = ∑
1≤i,j≤N
|(∇u)i,j| , (4.4.8)
o`u |(∇u)i,j| =
√|(∇u)1i,j|2 + |(∇u)2i,j|2. On introduit ́egalement une version discr`ete de l’op ́erateur de divergence. On le d ́efinit par analogie avec le cadre continu en posant
div = −∇∗ ,
o`u ∇∗ est l’op ́erateur adjoint de ∇, c’est-`a-dire
∀p ∈ Y,∀u ∈ X (−div p, u)X = (p,∇u)Y = (p1,∇1u)X + (p2,∇2u)X .
On peut alors v ́erifier que
(div p)i,j =
p1i,j − p1i−1,j si 1 <i<N
p1i,j si i = 1 −p1i−1,j si i = N
p2i,j − p2i,j−1 si 1 <j<N +
p2i,j si j = 1
(4.4.9)
−p2i,j−1 si j = N
Nous utiliserons aussi une version discr`ete du laplacien d ́efinie par
∆u = div ∇u .
4.4.1 Le probl`eme discret
On va remplacer le probl`eme (PROF) (dans la version (4.3.5)) par le probl`eme obtenu apr`es discr ́etisation suivant
min u∈X J(u) + 2ε 1u − ud 2X. (4.4.10) Il est facile de voir que ce probl`eme a une solution unique que nous allons caract ́eriser. On rappelle que |gi,j| =
√(gi,j1)2 + (gi,j2)2 et que la version discr`ete de la variation totale est donn ́ee (de mani`ere analogue au cas continu) par
J(u) = sup
ξ∈K 〈u,ξ〉 .
4.4. MOD`ELE DISCRET 65
o`u
K := {ξ = div (g) | g ∈ Y, |gi,j| ≤ 1 ,∀i, j } , est la version « discr`ete »de{ξ = div φ | φ ∈ C1c(Ω,R2) , φ ∞ ≤ 1 } . Donnons tout d’abord une caract ́erisation de la conjugu ́ee de Fenchel de J.
Th ́eor`eme 4.4.1 La transformée de Fenchel J∗ de la fonctionnelle J « variation totale » approchée définie sur X, est l’indicatrice de l’ensemble ̄K, o`u
K = {ξ = div (g) | g ∈ Y, |gi,j| ≤ 1 ,∀i, j } ,
Démonstration – Comme J est positivement homog`ene la conjugu ́ee J∗ de J est l’indicatrice d’un ensemble convexe ferm ́e ̃K (proposition A.3.3). Montrons d’abord que K ⊂ ̃K : soit u ∈ K. Par d ́efinition de J
J(u) = sup
ξ∈K 〈ξ,u〉 ,
o`u 〈·,·〉 d ́esigne le produit de dualit ́e entre BV (Ω) et son dual. Par cons ́equent 〈ξ,u〉−J(u) ≤ 0 pour tous ξ ∈ K et u ∈ BV (Ω). On d ́eduit donc que pour tout u∗ ∈ K
J∗(u∗) = sup
u∈K 〈u∗,u〉 − J(u) ≤ 0 . Comme J∗ ne prend qu’une seule valeur finie on a J∗(u∗)=0, et donc u∗ ∈ ̃K. Par cons ́equent K ⊂ ̃K et comme ̃K est ferm ́e : ̄K ⊂ ̃K .
En particulierJ(u) = sup
ξ∈K 〈u,ξ〉 ≤ sup
ξ∈ ̄K 〈u,ξ〉 ≤ sup
ξ∈ ̃K 〈u,ξ〉 = sup
ξ∈ ̃K 〈u,ξ〉 − J∗(ξ) = J∗∗(u) . Comme J∗∗ = J, il vient
sup ξ∈K 〈u,ξ〉 ≤ sup
ξ∈ ̃K 〈u,ξ〉 ≤ sup
ξ∈K 〈u,ξ〉 ,
et donc
sup ξ∈K 〈u,ξ〉 = sup
ξ∈ ̄K 〈u,ξ〉 = sup
ξ∈ ̃K 〈u,ξ〉 . (4.4.11) Supposons maintenant qu’il existe u∗ ∈ ̃K tel que u∗ /∈ ̄K. On peut alors s ́eparer strictement u∗ et le convexe ferm ́e ̄K. Il existe α ∈ R et uo tels que
〈uo,u∗〉 > α ≥ sup
v∈ ̄K 〈uo,v〉 . D’apr`es (4.4.11) il vient
sup ξ∈ ̃K 〈uo,ξ〉≥〈uo,u∗〉 > α ≥ sup
v∈ ̄K 〈uo,v〉 = sup
v∈ ̃K 〈uo,v〉. On a donc une contradiction : ̃K = ̄K . 2
66 CHAPITRE 4. QUELQUES MODELES `DE RESTAURATION D’IMAGE
4.4.2 Algorithme de projection de Chambolle
Le r ́esultat suivant donne la caract ́erisation attendue de la solution [Cham] :
Th ́eor`eme 4.4.2 La solution de (4.4.10 ) est donnée par
u = ud − PεK(ud) , (4.4.12)
o`u PK est le projecteur orthogonal sur K.
Démonstration – D’apr`es la d ́efinition du sous-diff ́erentiel, u est une solution de (4.4.10 ) est ́equivalent
0 ∈ ∂
(J(u) + 2ε 1u − ud 2X)
= u − ud
ε + ∂J(u) . Comme J est convexe, sci, propre on peut appliquer le corollaire A.3.3. Donc
ud ε − u
∈ ∂J(u) ⇐⇒ u ∈ ∂J∗(ud ε − u
) ⇐⇒ 0 ∈ −u + ∂J∗(ud ε − u
) .
Ceci est aussi ́equivalent `a
0 ∈ ud ε − u
− uε d+ 1ε∂J∗(ud ε − u
) .
On en d ́eduit que w = ud ε − u
est un minimiseur de
w − uε d2 2
+ 1εJ∗(w) .
Comme sur K. Comme J∗ est l’indicatrice PK(uε d) = 1εde PεK(uK, cela d), on implique peut conclure.
que ud ε − u
est la projection orthogonale de uε
2 dTout revient maintenant `a calculer
PεK(ud) = argmin { ε div (p) − ud 2X | |pi,j ≤ 1, i,j = 1,··· ,N } . On peut le r ́esoudre par une m ́ethode de point fixe : po = 0.
pn+1 i,j = pni,j + ρ (∇[ div pn − ud/ε])i,j
1 + ρ
∣∣∣(∇[ div pn − ud/ε])i,j∣∣∣
. (4.4.13)
Th ́eor`eme 4.4.3 Si le param`etre ρ dans (4.4.13) vérifie ρ ≤ 1/8, alors
ε div pn → PεK(ud) .
La solution du probl`eme est donc donn ́ee par
u = ud − ε div p∞ o`u p∞ = n→+∞lim pn .
Chapitre 5
M ́ethode des contours actifs
5.1 Rappels sur la g ́eom ́etrie des courbes planes
Dans ce qui suit on consid`erera des courbes planes param ́etr ́ees connexes, compactes.
5.1.1 Abscisse curviligne – longueur
On consid`ere une courbe param ́etr ́ee (Γ,Φ) : Φ est une application d’un intervalle I de R2. C’est une param ́etrisation de Φ :
Γ = {M(x, y) ∈ R2 | (x(t),y(t)) = Φ(t), t ∈ I }.
On suppose que I = [a, b] est un intervalle compact de R. Soit Σ l’ensemble des subdivisions σ de [a, b] :
σ = {to = a, t1,··· ,tn−1,tn = b } .
On pose
lσ(Γ) =
n−1∑k=0d(Φ(tn+1),Φ(tn)) ,
o`u d est la distance euclidienne dans R2.
D ́efinition 5.1.1 On dit que la courbe Γ est rectifiable si
sup σ∈Σlσ(Γ) < +∞ .
Ce nombre est la longueur de la courbe Γ : l(Γ). Elle est indépendante de la paramétrisation (réguli`ere) choisie.
Remarquons que la longueur v ́erifie la relation de Chasles. On suppose maintenant que la courbe (c’est-`a-dire Φ) est de classe C1.
67
68 CHAPITRE 5. METHODE ́DES CONTOURS ACTIFS
D ́efinition 5.1.2 Soit Γ une courbe rectifiable, connexe, compacte de classe C1. Soit to ∈ I. On pose
S(t) =
{ l(Γt,to) si t ≥ to, −l(Γto,t) si t ≤ to,
o`u Γa,b désigne la courbe obtenue lorsque t ∈ [a, b].
Th ́eor`eme 5.1.1 La fonction réelle S est dérivable sur I et
S (t) = |Φ (t)| = |dMdt (t)| ,
o`u |·| désigne la norme euclidienne de R2 et M(t) = M(Φ(t)) = M(x(t),y(t)).
Corollaire 5.1.1 Sous les hypoth`eses précédentes la longueur de la courbe s’exprime par la formule :
l(Γ) =
∫ ba |dMdt (t)|dt . (5.1.1)
La fonction S est continue dans l’intervalle I. Son image J = S(I) ⊂ R est donc un intervalle et la variable s = S(t) d ́ecrivant J est appel ́ee abscisse curviligne de la courbe Γ. Le point Mo = M(to) est l’origine de l’abscisse curviligne. On rappelle qu’un point r ́egulier de Γ est un point o `u Φ (t) = 0. Si l’ensembles des valeurs du param`etre t telles que M(t) soit un point r ́egulier de Γ est dense dans I, la fonction S est un hom ́eomorphisme de I dans J = S(I). On peut donc red ́efinir un param ́etrage Ψ = Φ◦S−1 de Γ appel ́e param ́etrage de Γ par une abscisse curviligne.
Supposons maintenant que la courbe Γ soit param ́etr ́ee par une abscisse curviligne : Φ : [a, b] → R2. La longueur de la courbe est une quantit ́e intrins`eque mais on peut l’ ́ecrire en fonction de la param ́etrisation : l(Γ) = l(Φ). On alors le r ́esultat suivant qui nous servira par la suite
Th ́eor`eme 5.1.2 Soit X un ensemble de fonctions C1 de R sur R2 représentant des courbes fermées (valeurs de la fonction et de sa dérivée égales aux bornes) . Alors la fonctionnelle
l : X → R+
Φ ↦→ l(Φ)
est Gˆateaux différentiable en toute fonction Φ dont le gradient ∇Φ est non nul et
∀h ∈ X dΦdl(Φ) · h =
∫ ba div( ∇Φ
|∇Φ|)(s)h(s)ds ,
o`u |·| désigne la norme euclidienne de R2.
5.1. RAPPELS SUR LA GEOM ́ETRIE ́DES COURBES PLANES 69
Démonstration – L’application en tout point x = 0 et
N : R2 → R+ d ́efinie par N(x) = |x| = √x21 + x22 est d ́erivable
∇N(x) = |x| x. Si Ψ est une fonction C1non nulle de [a, b| dans R2, le th ́eor`eme des fonctions compos ́ees donne
d|Ψ| dΨ (Ψ) = |Ψ| Ψ.
Soit h ∈ X : calculons la Gˆateaux-d ́eriv ́ee de l :
dΦdl(Φ) · h =
∫ ba
∇Φ |∇Φ|(s) · ∇h(s)ds .
Une int ́egration par parties coupl ́ee `a la condition aux limites donne (la courbe est ferm ́ee)
dΦdl(Φ) · h =
∫ ba div( |∇Φ|∇Φ
(s))h(s)ds .
2
5.1.2 Etude ́g ́eom ́etrique locale d’une courbe param ́etr ́ee
On rappelle que si M(t) est un point r ́egulier d’une courbe C1 un vecteur tangent `a la courbe en M est d ́efini par
T = dMdt (t) = Φ (t) . On appelle vecteur normal `a la courbe un vecteur N orthogonal `a T.
Supposons maintenant que la courbe est C2 et que tous ses points sont r ́eguliers. Le pa- ram ́etrage par l’abscisse curviligne est alors aussi C2. On suppose donc que la courbe Γ est param ́etr ́ee par une abscisse curviligne. Alors le vecteur tangent
T(s) = dMds (s) est d ́erivable et le vecteur dTds (s) est ind ́ependant du choix de l’abscisse curviligne s. D ́efinition 5.1.3 La courbure de Γ au point M(so) est le nombre réel
ρ(so) = |dTds (so)| ≥ 0 .
Un point M de Γ est bi-r ́egulier si et seulement si sa courbure est non nulle. Pr ́ecisons enfin tout cela dans un rep`ere du plan. Soit Γ une courbe C2, connnexe, dont tous les points sont r ́eguliers. On la param`etre par une abscisse curviligne s. Etant ́donn ́e un rep`ere ( ı, j) on pose :
T(s) = cos(φ(s)) ı + sin(φ(s)) j. La d ́eriv ́ee dφds(s) est ind ́ependante du choix du rep`ere et de φ.
70 CHAPITRE 5. M ́ETHODE DES CONTOURS ACTIFS
D ́efinition 5.1.4 La courbure alg ́ebrique de Γ au point M(s) est le nombre réel
ρa(s) = dφds(s) .
On peut alors montrer que
ρ(s) = |ρa(s)| et dTds (s) = ρ(s)N(s) .
5.2 M ́ethodes des contours actifs
5.2.1 Introduction
La segmentation permet d’isoler certaines parties de l’image qui pr ́esentent une forte corr ́elation avec les objets contenus dans cette image, g ́en ́eralement dans l’optique d’un post-traitement. Les domaines d’application sont nombreux : m ́edecine, g ́eophysique, g ́eologie, etc … Dans le domaine m ́edical, la segmentation d’images est extrˆemement compliqu ́ee. En effet, pour chaque organe (cerveau, cœur, etc …), l’approche est diff ́erente : l’outil de segmentation doit donc pouvoir s’adapter `a un organe particulier, suivant une modalit ́e d’acquisition particuli`ere (scanners, radiographie, Imagerie par R ́esonance Magn ́etique, …) et pour une s ́equence de donn ́ees particuli`ere. L’objectif est la quantification de l’information, par exemple, la volum ́etrie : volume d’une tumeur dans le cerveau, ́etude de la cavit ́e ventriculaire cardiaque, etc … C’est `a ce niveau que la segmentation de l’image est utilis ́ee. En g ́eophysique, la segmentation peut permettre d’isoler des objets du sous-sol (failles, horizons …) `a partir de donn ́ees sismiques dans le but, par exemple, de mod ́eliser ou d’exploiter un gisement.
En imagerie math ́ematique, deux types de segmentation sont exploit ́es :
- la segmentation par r ́egions qui permet de caract ́eriser les r ́egions d’une image pr ́esentant une structure homog`ene
- la segmentation par contours qui permet de d ́elimiter les diff ́erentes r ́egions par leurs fronti`eres. C’est `a ce dernier type de contourage que nous allons nous int ́eresser
Les points de contours sont les points de l’image pour lesquels la norme du gradient, dans la direction de ce gradient, est maximale. Un seuillage est r ́ealis ́e pour ne conserver que les points de variation de niveau de gris significative. Une question l ́egitime qui ́emerge alors, est comment d ́efinir le seuil, autrement dit, pour quel crit`ere de variation de niveau de gris un point sera qualifi ́e de point de contour ?
Un seuil choisi trop faible accordera l’ ́etiquette « point de contour » `a un nombre trop ́elev ́e de points tandis qu’un seuil trop ́elev ́e ne permettra d’extraire que les points de fort contraste et les contours d ́etect ́es ne seront plus connexes. La repr ́esentation math ́ematique des fronti`eres d’un objet ne sera plus r ́ealis ́ee dans ce cas.
Pour pallier cette difficult ́e, une r ́egularit ́e sur la mod ́elisation des contours doit ˆetre intro- duite : les contours seront assimil ́es `a des courbes poss ́edant des propri ́et ́es de r ́egularit ́e et satisfaisant le crit`ere de d ́etection ́enonc ́e pr ́ec ́edemment.
Dans ce chapitre on s’int ́eresse `a la m ́ethode des contours actifs (encore appel ́es « snakes») introduits par Kass, Witkin et Terzopoulos [KWT], m ́ethode qui int`egre cette notion de r ́egularit ́e
5.2. METHODES ́DES CONTOURS ACTIFS 71
des points de contour en introduisant une fonctionnelle interpr ́et ́ee en terme d’ ́energie pour les propri ́et ́es m ́ecaniques qu’elle revˆet. En effet, cette m ́ethode permet de faire ́evoluer en temps et en espace la repr ́esentation du mod`ele vers la solution du probl`eme de minimisation introduit dans la mod ́elisation. Ces m ́ethodes de contours actifs font appel `a la notion de corps ́elastique subissant des contraintes ext ́erieures. La forme prise par l’ ́elastique est li ́ee `a une minimisation d’ ́energie compos ́ee de deux termes :
- un terme d’ ́energie interne et
- un terme d’ ́energie externe.
Le minimum local obtenu par minimisation de cette fonctionnelle non-convexe est li ́e `a la condition initiale qui d ́efinit un voisinage de recherche du minimum.
De mani`ere g ́en ́erale, les difficult ́es ma jeures rencontr ́ees dans le processus de segmenta- tion par contours actifs sont les suivantes :
– le mod`ele est non-intrins`eque du fait de la param ́etrisation et n’est donc pas li ́e `a la
g ́eom ́etrie de l’objet `a segmenter. – forte d ́ependance du mod`ele `a la condition initiale, qui n’autorise pas `a choisir une condi-
tion initiale ́eloign ́ee de la solution. – connaissance de la topologie de/ou/des objets `a segmenter n ́ecessaire, ce qui implique (lorsqu’il y a plusieurs objets `a segmenter) l’utilisation de proc ́edures particuli`eres du fait de la param ́etrisation. – complexit ́e des images / ambiguıt ́e des donn ́ees correspondantes, en particulier lorsque les donn ́ees des images manquent et/ou que deux r ́egions de texture similaire sont adja- centes. – bruit sur les donn ́ees.
5.2.2 Mod ́elisation du probl`eme
On peut d ́efinir un contour actif ou « snake »comme une courbe ferm ́ee qui minimise son ́energie, influenc ́ee par une contrainte interne et guid ́ee par une force d’image qui pousse la courbe vers les contours pr ́esents dans l’image. L’espace des formes est l’ensemble Φ des courbes param ́etr ́ees (par l’abscisse curviligne) suivant :
Φ =
{
v | v : [0,1] s → R2
↦→ (x(s),y(s)) ,v(0) = v(1)
}
.
L’ ́energie du « snake » peut s’ ́ecrire sous la forme :
Es(v) = Ei(v) + Ee(v) .
Elle est constitu ́ee d’un terme de r ́egularisation interne (Ei ) et d’un terme de potentiel d’at- traction ou ajustement aux donn ́ees (Ee ). L’ ́energie interne peut se d ́ecomposer comme suit :
Ei(v) = 12
∫ 10
[α(s)|v (s)|2 + β(s)|v”(s)|2] ds (5.2.2)
72 CHAPITRE 5. METHODE ́DES CONTOURS ACTIFS
o`u |·|d ́esigne la norme euclidienne dans R2. α ∈ L∞(0,1)est le coefficient d’ ́elasticit ́e (r ́esistance `a l’allongement) et β ∈ L∞(0,1) le coefficient de rigidit ́e. Le premier terme de Ei(v) :
∫ 10 α(s)|v (s)|2 ds
p ́enalise la longueur du «snake» ( augmenter α(s) tend `a ́eliminer les boucles en r ́eduisant la longueur du contour). Le second terme ∫ 10 β(s)|v”(s)|2 ds
p ́enalise la courbure ( augmenter β(s) tend `a rendre le « snake » moins flexible). On fait donc apparaıtre dans cette expression des propri ́et ́es m ́ecaniques du comportement d’un ́elastique (d ́eriv ́ee du premier ordre) et d’une poutre (d ́eriv ́ee du second ordre). Les mod`eles d ́eformables se comportent comme des corps ́elastiques qui r ́epondent naturellement aux forces et aux contraintes qui leur sont appliqu ́ees. Notons que le choix β = 0 autorise n ́eanmoins les discontinuit ́es du second ordre. En ce qui concerne l’expression de l’ ́energie externe, plusieurs expressions li ́ees `a une fonction potentielle sont `a notre disposition. Les contours que l’on souhaite d ́eterminer sont :
– soit assimil ́es aux points de fort gradient de I , image donn ́ee : l’expression de Ee fait
alors apparaıtre la norme du gradient :
Ee(v) =
∫ 10 P(v(s))ds := −λ∫ 10 |∇I(v(s))|2 ds , (5.2.3)
avec λ > 0. – soit assimil ́es aux points de d ́eriv ́ee seconde nulle : si Gσ d ́esigne un filtre gaussien, on
peut choisir
Ee(v) =
∫ 10 P(v(s))ds := λ∫ 10 |Gσ ∗ ∆I(v(s))|n ds , (5.2.4)
La mod ́elisation du probl`eme que l’on propose revient donc `a trouver une fonction ̄v ∈ Φ ∩ X
min{Es(v) | v ∈ Φ ∩ X }. (5.2.5)
Ici X est un espace fonctionnel adapt ́e `a la formulation du probl`eme de minimisation. En effet, nous avons d ́ecrit les choses formellement jusqu’`a pr ́esent. Il est clair que pour d ́efinir les diff ́erentes ́energies, la fonction v doit ˆetre d ́erivable au sens des distributions et de d ́eriv ́ees int ́egrables. En pratique
X = H1(0,1) ou H2(0,1) .
Il s’agit donc d’un probl`eme de minimisation de fonctionnelle (minimisation de l’ ́energie) que l’on peut r ́esoudre au cas par cas en utilisant les th ́eor`emes et les m ́ethodes classiques d’opti- misation pr ́esent ́ees dans le chapitre ??, section A.1.3 par exemple.
5.2. METHODES ́DES CONTOURS ACTIFS 73
5.2.3 Conditions d’optimalit ́e
Supposons avoir d ́emontr ́e que le probl`eme (5.2.5) admet au moins une solution ̄v. Nous allons utiliser le th ́eor`eme d’Euler-Lagrange (condition d’optimalit ́e du premier ordre : cha- pitre ??, Th ́eor`eme A.1.5 ) pour d ́eterminer l’ ́equation aux d ́eriv ́ees partielles que satisfait ̄v. Supposons (pour simplifier) que X = H2(0,1) de sorte que v ∈ L2(0,1) et v” ∈ L2(0,1). On supposera ́egalement que la fonction P : X → L1(0,1) est d ́erivable (au sens des distributions ) et que sa d ́eriv ́ee est dans L1(0,1).(Cet ensemble se note W 1,1(0,1)). La fonctionnelle coˆut consid ́er ́ee est
Es(v) = 12
∫ 10
∫ 10 P(v(s))ds , sur l’espace X. Elle n’est en g ́en ́eral pas convexe de sorte que si l’on peut ( ́eventuellement) montrer l’existence d’une solution ̄u , celle ci n’est en g ́en ́eral pas unique. D’apr`es le Th ́eor`eme A.1.5, si ̄u est solution du probl`eme alors
∀v ∈ X dEdv s
( ̄u) · v = 0 .
Lemme 5.2.1 La fonctionnelle Es est Gˆateaux-différentiable et sa Gˆateaux-dérivée en ̄u, dans la direc- tion v est
dEdv s
( ̄u) · v =
[α(s)|v (s)|2 + β(s)|v”(s)|2] ds +
∫ 10
[α(s) ̄u (s)v (s) + β(s) ̄u”(s)v”(s) + ∇vP( ̄u(s))v(s)] ds .
Démonstration – La d ́emonstration, facile, est laiss ́ee en exercice. 2 Soit alors v ∈ D(0,1) : on obtient
∫ 10
[ dds
(αd ̄uds)(s)v(s) + dsd2 2
(βdds2 ̄u 2
)]
(s)v(s) + ∇vP( ̄u)(s)v(s)ds = 0 ,
c’est-`a-dire dds
(αd ̄uds)
+ dsd2 2
(βdds2 ̄u 2
)
+ ∇vP( ̄u)=0 ,
au sens des distributions. Il s’agit ensuite de donner des conditions aux limites : ici elles pro- viennent de la param ́etrisation par abscisse curviligne :
v(0) = v(1) = dvds(0) = dvds(1) . En d ́efinitive on est ramen ́e `a l’ ́etude et `a la r ́esolution (num ́erique) d’une ́equation aux d ́eriv ́ees partielles :
ds d(αd ̄uds)
+ dsd2 2
(βdds2 ̄u 2
) v(0) = v(1) = v ∈ Φ ∩ H2(0,1) + ∇vP( ̄u)=0 dans D (0,1)
dvds.
(0) = dvds(1)
(5.2.6)
74 CHAPITRE 5. M ́ETHODE DES CONTOURS ACTIFS
5.2.4 Un mod`ele dynamique
L’ ́equation ́ecrite ci-dessus exprime un ́etat statique ou ́etat d’ ́equilibre quand on est au mi- nimum d’ ́energie. Mais bien souvent, alors mˆeme qu’on connaıt l’existence d’un infimum on n’est pas capable de montrer l’existence d’un minimum et il faut approcher cet ́etat d’ ́equilibre plutˆot que le chercher exactement (il n’est parfois pas atteignable). Le principe des mod`eles d ́eformables et de consid ́erer que le contour que l’on cherche est l’ ́etat d’ ́equilibre d’un contour qui va ́evoluer avec le temps. Le probl`eme stationnaire est transform ́e en un probl`eme dyna- mique : on cherche la solution u(x, t) d’un probl`eme d’ ́evolution. La courbe cherch ́ee ̄u est alors donn ́ee par
̄u(x) = lim t→+∞u(t, x) . Une mani`ere simple d’imposer un mouvement au contour est d’imposer sa vitesse d’ ́evolution ∂u∂t (x, t) en posant ∂u∂t (x, t) = −∇uEs(u(x, t)) . (5.2.7) En effet, la famille u(t,·) de courbes (index ́ee par t) que l’on cherche doit ˆetre choisie de fac ̧on `a faire d ́ecroıtre l’ ́energie Es qui se rapprochera ainsi de son infimum. Si on fait (formellement) un d ́eveloppement de Es au premier ordre, on a pour δt > 0
E(u(t + δt,·)) − E(u(t,·)) ≃ 〈∇uE(u(t,·)),u(t + δt,·) − u(t,·)〉 . (5.2.8)
Le choix « u(t + δt,·) − u(t,·) = −δt∇uE(u(t,·) »montre qu’on fait bien d ́ecroıtre l’ ́energie. En faisant tendre δt vers 0 on obtient (5.2.7). Ceci conduit `a une ́equation aux d ́eriv ́ees partielles d’ ́evolution (en g ́en ́eral parabolique `a laquelle il convient d’ajouter une condition initiale
u(t = 0,x) = uo(x) (on se donne le contour de d ́epart)
et des conditions aux limites issues de l’analyse du probl`eme stationnaire.
On calcule ensuite cette solution et on s’int ́eresse `a son comportement symptotique (en temps long). En pratique, num ́eriquement il suffit de s’arrˆeter lorsque t est assez grand (c’est-`a-dire quand deux valeurs cons ́ecutives sont assez voisines ).
5.2.5 Un exemple
Nous allons illustrer ce qui pr ́ec`ede par un exemple. Supposons que la fonction α soit constante, que β soit nulle (on ne contraint pas la rigidit ́e du contour actif ). L’espace fonc- tionnel des courbes est alors X = H2(0,1) ∩ H1o(0,1) (en prenant en compte la condition aux limites v(0) = v(1) = 0. On supposera que l’image I est `a d ́eriv ́ee dans L2 et que
v ↦→ ∇I(v)
est continue de X dans L2(0,1). Choisissons dans un premier temps
P(v) = −12λ|∇I(v)|2 .
5.2. METHODES ́DES CONTOURS ACTIFS 75
La fonctionnelle d’ ́energie est
Es(v) = 12
∫ 10
[α|v (s)|2 − λ|∇I(v(s))|2] ds , v ∈ H2(0,1) ∩ Ho1(0,1).
On remarque imm ́ediatement que la fonctionnelle Es peut ne pas ˆetre minor ́ee ! ! On peut alors prendre pour P(v) une fonction dite de détection de contours moins « sommaire » que celle que nous venons de choisir, par exempleP(v) = 1 + |∇I(v)|1
2 .
Cette fonction a le m ́erite d’ˆetre positive et la minimiser revient `a maximiser son d ́enominateur donc le gradient de l’image. La fonctionnelle d’ ́energie est `a pr ́esent
Es(v) = 12
[α|v (s)|2 + 1 + |∇I(v)|λ
2]
ds , v ∈ H2(0,1) ∩ Ho1(0,1).
Cette fonctionnelle est minor ́ee, donc l’infimum existe. Toutefois il est encore extrˆemement d ́elicat de montrer qu’il est atteint, les propri ́et ́es de semi-continuit ́e et surtout de coercivit ́e de Es n’ ́etant pas ́evidentes. Dans ce cas, on adopte la d ́emarche « dynamique » comme dans la section pr ́ec ́edente. Calculons 〈∇uE(u(t,·)),u(t + δt,·) − u(t,·)〉pourδ > 0. Un calcul analogue `a la section pr ́ec ́edente donne∫ 1〈∇uE(u(t,·)),u(t + δt,·) − u(t,·)〉 = 0
∫ 10
[αu (t, s)[u (t + δt, s) − u (t, s)] + ∇uP(u)(t, s)[u(t + δt, s) − u(t, s)]] ds .
Une int ́egration par parties coupl ́ee aux conditions aux limites de la forme (courbe ferm ́ee) :
v(0) = v(1) = dvds(0) = dvds(1) donne
〈∇uE(u(t,·)),u(t + δt,·) − u(t,·)〉 =
∫ 10
[−αu (t, s) + ∇uP(u)(t, s)] [u(t + δt, s) − u(t, s)]ds .
Si on choisit
u(t + δt, s) − u(t, s) = −δt[−αu (t, s) + ∇uP(u)(t, s)] p.p. l’ ́energie va d ́ecroıtre (pour δt assez petit) grˆace `a l’approximation (5.2.8). On obtient
u(t + δt, s) δt − u(t, s)
= αu (t, s) − ∇uP(u)(t, s) p.p.
et par passage `a la limite lorsque δt → 0
∂u∂t − α∂∂s2u 2 ∀t > 0 u(t,0) + ∇vP(u)=0 dans D (0,1)
= u(t,1) = duds(t,0) = duds(t,1) , ∀s ∈ [0,1] u(0,s) = uo(s) donn ́ee.
(5.2.9)
C’est une EDP parabolique non lin ́eaire dont l’ ́etude est en g ́en ́eral classique.
76 CHAPITRE 5. M ́ETHODE DES CONTOURS ACTIFS
5.3 La m ́ethode des lignes de niveau (« Level set »)
Le principe des contours actifs est de faire ́evoluer une courbe. On a vu dans la section pr ́ec ́edente une formulation dynamique qui fait intervenir ∂u∂t c’est-`a-dire la vitesse d’ ́evolution du contour. On va donc s’int ́eresser `a la fac ̧on de faire ́evoluer la courbe et plus g ́en ́eralement `a la notion de propagation de fronts.
5.3.1 Propagation de fronts
On se donne une courbe plane (on peut aussi se placer en 3D avec une surface) ferm ́ee que l’on supposera r ́eguli`ere (on pr ́ecisera cela plus tard). Cette courbe partage le plan en deux r ́egions, l’int ́erieur et l’ext ́erieur. On oriente la courbe de fac ̧on `a d ́efinir une normale ext ́erieure n. On se donne la vitesse de propagation de la courbe : cette vitesse est port ́ee par la normale :
v = F n ,
car on supposera qu’il n’y pas de d ́eplacement tangentiel (pas de rotation ou de glissement par exemple). La fonction F d ́epend de plusieurs facteurs :
- locaux d ́etermin ́es par l’information g ́eom ́etrique locale (comme la courbure et la nor-
male)
- globaux : ce sont ceux qui d ́ependent de la forme (globale) ou de la position du front. La
vitesse peut par exemple inclure une int ́egrale sur tout le domaine.
Le probl`eme le plus important est la mod ́elisation du front de d ́eplacement de la courbe, c’est- `a-dire la formulation de l’expression de la vitesse F. Nous allons supposer ici que la vitesse est connue. On donnera ensuite quelques exemples.
On se donne donc une courbe r ́eguli`ere (par exemple C2 dont tous les points sont r ́eguliers) Γo de R2. La famille obtenue par d ́eplacement de long de la normale `a la vitesse F est not ́ee Γt. Si φ = (x, y) est la position d’un point de la courbe et n la normale ext ́erieure (unitaire ) `a la courbe, on aura F = n · φ. On consid`ere une param ́etrisation de la courbe Γt par une abscisse curviligne :
φ(s, t)=(x(s, t),y(s, t)) o`u s ∈ [0,S] et φ(0,t) = φ(S, t) .
On param`etre la courbe de fac ̧on que l’int ́erieur soit `a gauche de la direction des s croissants. n(s, t) est une param ́etrisation de la normale ext ́erieure et κ(t, s) une param ́etrisation de la courbure.
5.3. LA M ́ETHODE DES LIGNES DE NIVEAU (« LEVEL SET ») 77
On rappelle que dans ce cas la courbure est donn ́ee par
κ(t, s) =
[yssxs − xssys (x2s + ys2)3/2
](t, s)
o`u xs d ́esigne la d ́eriv ́ee par rapport `a s et xss d ́esigne la d ́eriv ́ee seconde par rapport `a s.
On se restreint dans ce qui suit `a des vitesses ne d ́ependant que de la courbure locale
F = F(κ) .
On obtient alors les ́equations donnant le mouvement de Γt :
(yssxs−xssys (x2s+ys2)3/2
)( ys
(x2s+ys2)1/2) xt = F
yt = F
(yssxs−xssys
(5.3.10)
(x2s+ys2)3/2
)( (x2s+yxs s2)1/2)
.
78 CHAPITRE 5. METHODE ́DES CONTOURS ACTIFS
La variation totale de la courbe Γt (ou de la fonction qui la param`etre φ) est aussi sa lon- gueur. Elle est donn ́ee par
l(t) =
∫ S0 |κ(s, t)|(x2s + ys2)1/2 ds .
La proposition suivante donne une id ́ee de l’ ́evolution du front lorsque la courbure s’annule (courbe non convexe).
Proposition 5.3.1 On consid`ere un front évoluant `a la vitesse F(κ) le long du champ de vecteurs normaux. Supposons que la courbe initiale Γo est simple, réguli`ere et non-convexe, de sorte que κ(s,0) change de signe. Supposons que F est deux fois différentiable et que κ l’est aussi pour 0 ≤ s ≤ S et 0 ≤ t ≤ T. Alors, pour tout 0 ≤ t ≤ T
- Si F (κ) ≤ 0 (resp. ≥ 0) d`es que κ = 0 alors
dldt ≤ 0 (resp. ≥ 0 )
- Si F (κ) < 0 (resp. > 0) et κ = 0 d`es que κ = 0 alors
dldt < 0 (resp. > 0 )
5.3.2 La m ́ethode « Level set »
La m ́ethode des lignes de niveaux ( «Level set») permet de s’affranchir de la param ́etrisation des courbes en terme d’abscisse curviligne. Le prix `a payer est une augmentation de la dimen- sion de l’espace dans lequel on travaille.
L’id ́ee consiste `a consid ́erer une courbe plane comme une ligne de niveau d’une surface 3D d’ ́equation z − Φ(x, y)=0. On choisit le ligne de niveau 0, c’est-`a-dire l’intersection de la surface 3D, avec le plan z = 0.
5.3. LA M ́ETHODE DES LIGNES DE NIVEAU (« LEVEL SET ») 79
Plus pr ́ecis ́ement on va chercher une fonction
Ψ : R3 → R
v ́erifiant
Γt = { X(t)=(x(t),y(t)) ∈ R2 | Ψ(t, X(t)) = 0 } . (5.3.11) Pour cela, on va partir d’une courbe initiale Γo d’ ́equation ψo, on pose
Ψ(0,X(0)) = ψo(X(0)) ,
et on va faire ́evoluer la surface en prenant en compte la vitesse d’ ́evolution du front donn ́ee par X (t) = dXdt . L’ ́equation (5.3.11) se traduit par
∀t ≥ 0 Ψ(t, X(t)) = 0 .
Comme nous cherchons une surface r ́eguli`ere (au moins C1) nous pouvons d ́eriver et la formule de composition donne∀t ≥ 0 ∂Ψ∂t (t, X(t)) + ∇Ψ(t, X(t)) · X (t)=0 ,
o`u ∇Ψ d ́esigne le gradient par rapport `a X. Or nous connaissons la vitesse de propagation du front :
X (t) · n = F .
80 CHAPITRE 5. METHODE ́DES CONTOURS ACTIFS
Le vecteur normal est donn ́e par
n = ∇Ψ
|∇Ψ|(X(t)) .
On obtient donc une ́equation d’ ́evolution pour Ψ :
∂Ψ∂t (t, X) + F(t, X)|∇Ψ(t, X)| = 0 , Ψ(0,X) donn ́ee
(5.3.12)
Une fois l’ ́equation r ́esolue ( ! !) on obtient la courbe Γt en r ́esolvant
Ψ(t, X)=0 .
Pour certaines valeurs de F, l’ ́equation (5.3.12) est une ́equation d’Hamilton-Jacobi :
∂Ψ∂t + H(∇Ψ)=0 ,
o`u H est le hamiltonien.
5.3.3 Application `a l’image
Reprenons l’exemple de la section 5.2.5. L’ ́equation (5.2.9) donne la vitesse d’ ́evolution du front :
F = ∂u∂t = α∂∂s2u
2 − ∇vP(u) .
de ∂∂s2u 2 Ψ est : P(Ψ). ̃`a une constante L’ ́equation pr`es issue la de courbure la m ́ethode κ au point des consid ́er ́e et P(u) peut se r ́e ́ecrire en fonction
lignes de niveau est donc la suivante ∂Ψ∂t + kdiv( ∇Ψ
|∇Ψ|)|∇Ψ| − P(Ψ)|∇Ψ| ̃= 0 .
5.4 Le mod`ele des ballons ( « Balloons »)
Le mod`ele des ballons ( « Balloons ») a ́et ́e introduit par Laurent Cohen [Coh]. En effet un des principaux probl`emes des « snakes » provient de la condition initiale. En effet, si la condition initiale n’est pas assez proche de la solution, le contour n’ ́evolue pas suffisamment et tend `a se localiser sur un minimum local non significatif. L’int ́erˆet du mod`ele des ballons r ́eside dans la r ́esolution de ce probl`eme. On ajoute une force suppl ́ementaire que l’on peut appeler « force de gonflage ». La courbe est assimil ́ee `a un ballon que l’on gonfle. Deux possibilit ́es sont alors envisageables :
- Soit la nuance (dans le cas des intensit ́es) est assez forte et la courbe s’arrˆete.
- Soit la nuance est trop faible et la courbe la surmonte pour aller chercher plus loin.
télécharger gratuitement méthodes mathématiques pour le traitement d’image
