méthodes mathématiques pour le traitement d’image

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. 

1. 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). 

1. 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 ∈ Ω 

1. 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 

1. 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) ≤ ε .