TNImages 2012 Vol2 V1 .pdf



Nom original: TNImages_2012_Vol2_V1.pdf
Titre: Microsoft Word - TNImages_2012_Vol2_V1.doc
Auteur: Florence Rossant

Ce document au format PDF 1.4 a été généré par PScript5.dll Version 5.2 / Acrobat Distiller 8.1.0 (Windows), et a été envoyé sur fichier-pdf.fr le 07/06/2015 à 17:19, depuis l'adresse IP 88.181.x.x. La présente page de téléchargement du fichier a été vue 534 fois.
Taille du document: 3.7 Mo (84 pages).
Confidentialité: fichier public




Télécharger le fichier (PDF)










Aperçu du document


SITe – Signal Image Télécoms

Traitement et Analyse des
Images Numériques
2011-2012 – V2

Florence Rossant

Florence Rossant

08/03/2012

TNImages_2012_Vol2_V1

SITe – Signal Image Télécoms

Traitement et Analyse des
Images Numériques
2011-2012

VOLUME 2

Florence Rossant

-3-

-4-

Table des Matières

VOLUME 2 .......................................................................................................... 3
CHAPITRE IV ......................................................................................................................... 9
SEGMENTATION – EXTRACTION DE CONTOURS ........................................................................ 9
1.
Introduction à la segmentation ......................................................................................... 9
2.
Méthodes locales différentielles..................................................................................... 11
2.1.
2.2.
2.2.1.
2.2.2.
2.2.3.
2.3.
2.3.1.
2.3.2.
2.3.3.
2.4.

Principe ...................................................................................................................................................... 11
Approche gradient ...................................................................................................................................... 14
Noyaux de Prewitt ................................................................................................................................... 14
Opérateur de Sobel.................................................................................................................................. 15
Opérateurs de Kirsch et compas gradient : filtres directionnels............................................................. 16
Approche Laplacien ................................................................................................................................... 16
Opérateur Laplacien ............................................................................................................................... 17
Opérateur LOG : approche Marr-Hildreth............................................................................................. 18
Opérateur DoG ....................................................................................................................................... 20
Discussion .................................................................................................................................................. 20

3.

Approches analytiques ................................................................................................... 22

3.1. Filtre de Canny ........................................................................................................................................... 22
3.1.1. Modélisation du contour ......................................................................................................................... 22
3.1.2. Formalisation du critère 1 : bonne détection .......................................................................................... 22
3.1.3. Formalisation du critère 2 : bonne localisation...................................................................................... 23
3.1.4. Formalisation du critère 3 : unicité de la réponse .................................................................................. 24
3.1.5. Résolution................................................................................................................................................ 24
3.1.6. Mise en œuvre du filtre de Canny............................................................................................................ 26
3.2. Filtre de Deriche......................................................................................................................................... 26
3.2.1. Filtre mono-dimensionnel ....................................................................................................................... 26
3.2.2. Filtre bi-dimensionnel ............................................................................................................................. 27
3.2.3. Implémentation récursive ........................................................................................................................ 28
3.2.4. Exemples ................................................................................................................................................. 30

4.

Extraction des points de contour et post-traitements ..................................................... 31

4.1. Suppression des non-maxima ..................................................................................................................... 31
4.2. Seuillage par hystérésis .............................................................................................................................. 32
4.3. Fermeture des contours .............................................................................................................................. 34
4.3.1. Algorithme 1: recherche locale, critère de ressemblance ....................................................................... 34
4.3.2. Algorithme 2: gradient moyen maximal .................................................................................................. 34

5.

Contours actifs................................................................................................................ 35

5.1.
5.2.
5.2.1.
5.2.2.
5.2.3.
5.3.
5.4.
5.5.
5.5.1.
5.5.2.
5.6.

Modélisation............................................................................................................................................... 35
Résolution (cas d’un contour fermé, α et β constants) ............................................................................... 37
Equation d’Euler ..................................................................................................................................... 37
Modélisation discrète .............................................................................................................................. 37
Approche variationnelle pour la résolution ............................................................................................ 39
Algorithme (cas d’un contour fermé, α et β constants) .............................................................................. 39
Exemples .................................................................................................................................................... 41
Cas des contours non fermés ...................................................................................................................... 41
Contours à extrémités libres.................................................................................................................... 41
Contours à extrémités fixes ..................................................................................................................... 42
Gradient Vector Flow................................................................................................................................. 42

-5-

Traitement et Analyse des Images Numériques

6.

Transformée de Hough................................................................................................... 45

CHAPITRE V......................................................................................................................... 49
SEGMENTATION DE REGIONS ................................................................................................... 49
1.
Segmentation par seuillage............................................................................................. 50
1.1. Principe ...................................................................................................................................................... 50
1.2. Choix du seuil............................................................................................................................................. 51
1.2.1. Valeur minimale de l’histogramme entre les deux pics principaux......................................................... 51
1.2.2. Méthode de Otsu : maximisation de la variance inter-classes ................................................................ 51

2.
3.
4.

k-moyennes .................................................................................................................... 53
Croissance de régions..................................................................................................... 55
Division et rassemblement ............................................................................................. 56

4.1.
4.2.

Phase 1 : division........................................................................................................................................ 56
Phase 2 : rassemblement............................................................................................................................. 58

5.
6.

Traitement multi-résolution............................................................................................ 61
Classification statistique................................................................................................. 66

6.1. Apprentissage ............................................................................................................................................. 67
6.1.1. Estimation non paramétrique par la méthode des tM plus proches voisins ............................................. 67
6.1.2. Estimation paramétrique :....................................................................................................................... 68
6.2. Décision...................................................................................................................................................... 70
6.3. Relaxation .................................................................................................................................................. 70

7.

Etiquetage et post-traitement.......................................................................................... 71

7.1.
7.2.

Etiquetage à partir des contours ................................................................................................................. 71
Améliorations, fusion de régions................................................................................................................ 73

ANNEXES .................................................................................................................................. 75
1.
Discrétisation des opérateurs différentiels ..................................................................... 75
1.1.
1.2.
1.3.

Discrétisation des dérivées ......................................................................................................................... 75
Grille .......................................................................................................................................................... 76
Schéma numérique ..................................................................................................................................... 76

2.

Filtre de Canny ............................................................................................................... 77

2.1. Rappel : filtrage d’un bruit gaussien .......................................................................................................... 77
2.2. Calcul des critères de Canny dans le cas général ....................................................................................... 77
2.2.1. Critère 1 : bonne détection...................................................................................................................... 78
2.2.2. Critère 2 : bonne localisation.................................................................................................................. 78
2.2.3. Critère 3 : unicité de la réponse.............................................................................................................. 79
2.3. Calcul des critères de Canny dans le cas d’un échelon .............................................................................. 80

REFERENCES ............................................................................................................................ 83
1.
Livres.............................................................................................................................. 83
2.
Articles ........................................................................................................................... 83
3.
Liens internet.................................................................................................................. 83

08/03/2012

-6-

Table des Matières

08/03/2012

-7-

-8-

Chapitre IV
Segmentation – Extraction de contours

1.

Introduction à la segmentation

Le but de la segmentation d’image est de détecter et d’isoler les différents objets
composant l’image. Prenons un exemple très simple :

FIG. IV-1. Image à segmenter

Cette image est manifestement constituée de 3 composantes, les 2 pions et le fond de
l’image. Pour les isoler, on peut imaginer deux techniques : soit on considère que les objets
sont délimités par de fortes transitions de niveaux de gris, et on cherche des opérations de
filtrage de type gradient permettant de les mettre en évidence ; soit on considère que les pixels
constituant un même objet ont des propriétés de similarité, par exemple des niveaux de gris
voisins, et on définit un critère de ressemblance permettant de les réunir dans une même
composante. Dans le premier cas, il s’agit d’une segmentation par extraction de contour, dans
le second on parle d’une segmentation de régions.
La figure FIG. IV-2 illustre la détection de contours. Cet exemple ultra simple montre
cependant que la segmentation est une opération délicate. Colonne de gauche, on a appliqué
les noyaux de Prewitt pour calculer le module du gradient (a). Puis on a seuillé pour mettre en
blanc les pixels de gradient supérieur au seuil et à zéro les autres (b). On observe la difficulté
de choisir un seuil qui permette de sélectionner les pixels de contour et rien que ces pixels :

-9-

Traitement et Analyse des Images Numériques




Soit il existe des contours « parasites », comme des pixels blancs internes au pion
de gauche,
Soit les contours ne sont pas fermés, et donc les composantes de l’image ne sont
pas correctement délimitées.

(a) Détection des contours
(b) Seuillage de l’image (a) pour 2 seuils différents
FIG. IV-2. Segmentation par extraction de contour

Considérons maintenant la segmentation de région obtenue par seuillage de
l’histogramme de l’image source. L’histogramme (FIG. IV-3) montre 3 pics, correspondant
aux trois composantes de l’image. On applique sur l’image source im la transformation
suivante pour obtenir l’image seuillée ims :
⎧1 si im(i , j ) < S
ims(i , j ) = ⎨
⎩0 si im(i , j ) ≥ S

Eq. IV-1

FIG. IV-3. Histogramme de l’image source normalisée (FIG. IV-1)

(a) S=0.47

(b) S=0.50

FIG. IV-4. Histogramme de l’image source et fonction de seuillage

08/03/2012

- 10 -

Chapitre IV : Segmentation – Extraction de contours

On voit cette fois la difficulté à fixer un seuil qui permette de détecter les pions dans leur
intégralité, et rien que les pions. Comme les pics de l’histogramme ne sont pas séparés, ce
seuil n’existe pas, et on ne peut que minimiser le taux de pixels mal classés. Cet exemple est
particulièrement révélateur du problème des ombres et des reflets.
Dans la suite de ce chapitre, nous explorerons les méthodes de segmentation fondées
sur l’extraction des contours. La détection de régions sera exposée dans le Chapitre V. Dans
les deux cas, nous verrons plusieurs méthodes et leurs limitations. Il est en effet impossible de
mettre au point un unique algorithme qui fonctionnerait dans tous les cas, et ceux qui seront
proposés devront en pratique être adaptés au problème particulier traité.
Nous verrons également dans le Chapitre V les algorithmes qui permettent d’étiqueter
les régions trouvées, à partir d’une détection de contours ou de région. Etiqueter consiste à
donner une même étiquette (label) à tous les pixels appartenant à une même composante de
pixels connexes. Cela aboutit à une carte, où tous les pixels appartenant à une même région
sont affectés d’une couleur identique (étiquette identique) (Fig. IV-5). C’est à partir de cette
carte de région et de l’image d’origine que l’on pourra passer à une description de plus haut
niveau, permettant de représenter les régions par un codage plus compact, d’en extraire des
caractéristiques, etc. (Volume 3, Chapitre VIII).

Fig. IV-5 : Etiquetage des composantes connexes

Il faut également garder à l’esprit que la segmentation d’une image est une étape
primordiale et très délicate, et qu’elle conditionne le succès des opérations ultérieures mises
en œuvre dans tout système d’analyse automatique d’images.

2.

Méthodes locales différentielles

2.1.

Principe

Un contour est défini comme la frontière entre les objets d’une image. Il correspond à
une forte transition d’intensité lumineuse (niveau de gris). Mathématiquement un point de
contour va se caractériser :
• dans le domaine fréquentiel, par de hautes fréquences.
• dans le domaine spatial, par un maximum du gradient ou un passage à zéro de la

08/03/2012

- 11 -

Traitement et Analyse des Images Numériques

dérivée seconde dans la direction normale au gradient.

FIG. IV-6. Contours et dérivées

Rappelons la définition du gradient pour la fonction f(x,y) :
⎡ ∂f ( x , y ) ∂ f ( x , y ) ⎤
,
∇f ( x , y )=⎢
∂ y ⎥⎦
⎣ ∂x

Eq. IV-2

2

⎛ ∂ f ( x, y ) ⎞ ⎛ ∂ f ( x, y ) ⎞
⎟⎟ + ⎜⎜
⎟⎟
Module : G ( x , y ) = ∇f ( x , y ) = ⎜⎜
⎝ ∂x ⎠ ⎝ ∂y ⎠

Argument :

⎛ ∂f ( x , y )
⎝ ∂y

ϕ ( x , y ) = tan −1 ⎜⎜

∂f ( x , y ) ⎞

∂x ⎟⎠

2

Eq. IV-3

Eq. IV-4

Considérons donc une image f(x,y). On calcule suivant Eq. IV-2 le gradient en chaque
coordonnée (x,y). Les points de contours sont détectés en recherchant les maxima locaux de la
norme du gradient (Eq. IV-3). La direction du gradient, c'est-à-dire l’angle formé par la
normale du contour avec l’axe des x, est donnée par Eq. IV-4. Soit M un point de contour. Le
r ∇f
vecteur n =
est un vecteur unitaire normal au contour en M(x,y) :
∇f

FIG. IV-7. Contour et orientation

08/03/2012

- 12 -

Chapitre IV : Segmentation – Extraction de contours

r ∇f ⎡cos ϕ ⎤
n=
=
∇f ⎢⎣ sin ϕ ⎥⎦

Eq. IV-5

L’orientation du contour correspond à l’angle formé par la tangente au contour et l’axe des x.
Cet angle vaut :

θ=

π
2

Eq. IV-6

−ϕ

Calculons maintenant les dérivées du gradient dans la direction normale au contour,
r
c’est-à-dire dans la direction n :
r ∂f
∂f
∂f
= ∇f .n = cos ϕ + sin ϕ
∂n
∂x
∂y

Eq. IV-7

∂2 f ∂2 f
∂2 f
∂2 f
2
=
cos
ϕ
+
2
cos
ϕ
sin
ϕ
+
sin 2 ϕ
2
2
2
∂n
∂x
∂x∂y
∂y

Eq. IV-8

La recherche d’un extremum du gradient revient à rechercher le passage par zéro de la dérivée
seconde directionnelle (Eq. IV-8). On peut montrer que pour des courbures faibles, cela
revient à rechercher le passage par zéro du Laplacien défini par :
∂2 f ∂2 f
Δf ( x , y ) = ∇ f ( x , y ) = 2 + 2
∂x
∂y

Eq. IV-9

2

Donnons une explication intuitive de cela. Le calcul du Laplacien est invariant par
changement de repère euclidien :

Δf ( x , y ) = Δf (n ,t ) =

∂2 f ∂2 f
+
∂n 2 ∂t 2

Eq. IV-10

La dérivée seconde suivant t est en relation avec la courbure locale du contour. Pour
information :
∂2 f ∂2 f
∂2 f
∂2 f
2
=
sin
ϕ

2
cos
ϕ
sin
ϕ
+
cos 2 ϕ
∂t 2
∂x 2
∂x∂y
∂y 2

Pour une courbure faible du contour, on a

Δf (x , y ) ≈ 0 .

Eq. IV-11

∂2 f
∂2 f

0
.
Par
ailleurs
≈ 0 et donc
∂t 2
∂n 2

Une seconde méthode de détection de contours consiste donc à rechercher les passages
par 0 du Laplacien Δf ( x , y ) (Eq. IV-10).
Nous allons maintenant étudier comment on peut concrètement estimer le gradient et
le Laplacien d’une image discrète, par des opérateurs différentiels.

08/03/2012

- 13 -

Traitement et Analyse des Images Numériques

2.2.

Approche gradient

On considère un repère de l’image f défini par l’origine au coin en haut à droite, l’axe
x vertical orienté de haut en bas, l’axe y horizontal orienté de gauche à droite.
L’estimation du vecteur gradient est réalisée en chaque point de l'image. La dérivée est
approximée par la méthode des différences finies (Annexes, Section 1). Par exemple la
dérivée par rapport à x au point f(x,y) est approximée par :

∂f ( x , y ) 1
= ( f ( x + 1, y ) − f ( x − 1, y ))
∂x
2

Eq. IV-12

⎛ -1 ⎞
⎜ ⎟
Ceci correspond, au facteur près, à une convolution spatiale de noyau ⎜ 0 ⎟ .
⎜ + 1⎟
⎝ ⎠
De même la dérivée partielle suivant y correspond à une convolution spatiale de noyau
(− 1 0 1) .

2.2.1. Noyaux de Prewitt
En traitement d’image, on réalise cette approximation sur des voisinages locaux de
3*3 pixels. Les noyaux 3 x 3 de Prewitt sont définis par :
⎛ − 1 0 + 1⎞


S y = ⎜ − 1 0 + 1⎟
⎜ − 1 0 + 1⎟



⎛ − 1 − 1 − 1⎞


Sx = ⎜ 0
0
0 ⎟
⎜ + 1 + 1 + 1⎟


Dérivée suivant x

Dérivée suivant y

L’algorithme de détection de contour est le suivant :
1. Convoluer l’image avec le noyau Sx, et le noyau Sy. Soient imx et imy les images
résultats.
2. Calculer le module : ∀( x , y ),img( x , y ) = imx( x , y )2 + imy( x , y )2
3. Rechercher les maxima locaux.
On peut remarquer que ces noyaux réalisent un lissage dans la direction horizontale pour Sx et
dans la direction verticale pour Sy. La figure FIG. IV-8 illustre les résultats obtenus, après
calibration de chaque image :

Image source

imx

imy

FIG. IV-8. Extraction des contours par les noyaux de Prewitt

08/03/2012

- 14 -

img

Chapitre IV : Segmentation – Extraction de contours

Notons que les images imx et imy ont des valeurs positives et négatives. La figure FIG. IV-8
montre ces images après calibration : le zéro correspond au gris moyen, les pixels plus
sombres correspondent à des valeurs négatives de la dérivée, les pixels plus clairs à des
valeurs positives.
Le problème majeur est que le bruit est également caractérisé par de fortes transitions,
et sera donc également très bien détecté ! Un deuxième problème est la recherche des maxima
locaux. Le plus simple consiste à seuiller l’image du module de gradient afin de ne conserver
que les plus fortes valeurs. On voit sur l’exemple ci dessous (FIG. IV-9) que ce choix est
critique : un fort seuil ne permet pas d’obtenir des contours fermés. Un seuil trop faible donne
beaucoup de points de contours parasites. Nous verrons dans les paragraphes 4.1 et 4.2 une
méthode de seuillage plus robuste.

img(x,y)

imc obtenue par seuillage de img: imc(x,y)=1 si img(x,y)>S, 0 sinon
(pour 4 valeurs croissantes de S)
FIG. IV-9. Image de gradient seuillée pour des seuils croissants

2.2.2. Opérateur de Sobel
Cet opérateur donne des résultats plus satisfaisants que les noyaux de Prewitt dans le
cas d’une image bruitée. La méthode est identique mais les noyaux utilisés sont définis par :
⎛ − 1 − 2 − 1⎞
⎛ − 1 0 + 1⎞




Sy = ⎜− 2 0 + 2⎟
Sx = ⎜ 0
0
0 ⎟
⎜ + 1 + 2 + 1⎟
⎜ − 1 0 + 1⎟




Ceci revient à appliquer un noyau intégrateur avant ou après l’opérateur différentiel :
moyenne d’un opérateur gradient appliqué aux 4 sous-fenêtre 2*2 du voisinage local 3*3 de
chaque pixel.

En effet, soit un voisinage 3*3 de pixels de l’image source, prenant les valeurs
a, b, …, i. Appliquons le noyau de Sobel Sx. Le pixel central de valeur e prend la nouvelle
valeur g-a+2(h-b)+i-c :
a

b

d

e

g

h

⎛ − 1 − 2 − 1⎞


f *⎜ 0
0
0 ⎟ = g − a + 2(h − b ) + i − c
⎜ + 1 + 2 + 1⎟
i



c

Appliquons maintenant un noyau dérivateur 2*2 dans les quatre sous-voisinages 2*2 :
a

b

d

e

g

h

08/03/2012

c

⎛ − 1 − 1 ⎞ d − a + e − b e − b + f − c somme
f * ⎜⎜
⎟⎟ =
⎯⎯⎯→ g − a + 2(h − b ) + i − c
g

d
+
h

e
h

e
+
i

f
+
1
+
1

i ⎝
- 15 -

Traitement et Analyse des Images Numériques

La somme des 4 dérivées est égale au résultat obtenu par application directe de Sx. On
comprend donc que le noyau de Sobel peut s’interpréter comme la moyenne de quatre
dérivées locales. Cette moyenne lisse les fortes dérivées dues à un bruit non directif et très
local. L’opérateur de Sobel est donc moins sensible au bruit.
2.2.3. Opérateurs de Kirsch et compas gradient : filtres directionnels
Cette fois, on cherche à définir un ensemble de noyaux qui représentent une
modélisation discrète d’un contour idéal dans une direction spécifique. Pour un noyau de
taille 3*3, il y a 8 directions principales. On convolue l’image par ces 8 masques et on retient
la réponse maximale. Cela revient en fait à choisir le masque qui ressemble localement le plus
à l’image, la convolution revenant à un calcul de corrélation. L’algorithme de filtrage est
donc :
1. Convoluer par les huit masques Mk en chaque point (x,y) de l’image f :
S k ( x , y ) = f ( x , y )* M k , k = 1,...,8
2. Image de module : choisir la réponse maximale :
Mag ( x , y ) = Max S k ( x , y ) = S n ( x , y )
k

3. Image de direction : numéro du masque qui a donné la réponse maximale
Arg( x , y ) = n
Pour l’opérateur de Kirsch, les masques se déduisent du premier par rotations de 45°
⎛− 5 − 5 − 5⎞


0
3 ⎟
⎜ 3
⎜ 3
3
3 ⎟⎠

3
3 ⎞
⎛ 3


0
3 ⎟
⎜ 3
⎜− 5 − 5 − 5⎟



⎛3 − 5 − 5⎞


⎜3 0 − 5⎟
⎜3 3
3 ⎟⎠

3 3⎞
⎛ 3


⎜ − 5 0 3⎟
⎜ − 5 − 5 3⎟



⎛3 3 − 5⎞


⎜3 0 − 5⎟
⎜3 3 − 5⎟


⎛ − 5 3 3⎞


⎜ − 5 0 3⎟
⎜ − 5 3 3⎟



3 ⎞
⎛3 3


⎜3 0 − 5⎟
⎜3 − 5 − 5⎟


⎛ − 5 − 5 3⎞


⎜ − 5 0 3⎟
⎜ 3
3 3 ⎟⎠


Il en est de même pour l’opérateur compas gradient :
⎛− 1

⎜ 1
⎜ 1

⎛ 1

⎜ 1
⎜− 1


− 1 − 1⎞

−2 1 ⎟
1
1 ⎟⎠
1 ⎞

−2 1 ⎟
− 1 − 1 ⎟⎠
1

⎛ 1 − 1 − 1⎞


⎜ 1 − 2 − 1⎟
⎜1 1
1 ⎟⎠

1 1⎞
⎛ 1


⎜ − 1 − 2 1⎟
⎜ − 1 − 1 1⎟



⎛ 1 1 − 1⎞


⎜ 1 − 2 − 1⎟
⎜ 1 1 − 1⎟


⎛ − 1 1 1⎞


⎜ − 1 − 2 1⎟
⎜ − 1 1 1⎟



1 ⎞
⎛1 1


⎜ 1 − 2 − 1⎟
⎜ 1 − 1 − 1⎟


⎛ − 1 − 1 1⎞


⎜ − 1 − 2 1⎟
⎜ 1
1 1 ⎟⎠


Remarquons que la somme des coefficients des masques dérivateurs présentés est
toujours nulle.

2.3.

Approche Laplacien

Nous avons vu précédemment (paragraphe 2.1) que les contours correspondent aux
passages par zéro de la dérivée seconde directionnelle (normale au contour), et que ces

08/03/2012

- 16 -

Chapitre IV : Segmentation – Extraction de contours

passages par zéros coïncident avec ceux du Laplacien en cas de faible courbure. En pratique,
on préfère calculer le Laplacien, beaucoup moins complexe à mettre en œuvre que la dérivée
seconde directionnelle, malgré la petite erreur commise. Nous allons dans ce paragraphe
développer cette approche.
2.3.1. Opérateur Laplacien
Cette fois, on recherche une estimation de Δf (x , y ) (Eq. IV-10). Considérons un
voisinage local de l’image f(x,y) :
a

b

c

d

e

f

g

h

i

Calculons une approximation du Laplacien au point (x,y) de niveau de gris e :


:
∂y

:
∂x

(e-d)
(f-e)

∂2
:(f-e) − (e-d) = d + f − 2e :
∂y 2

(e-b)
(h-e)


: (h − e ) − (e − b ) = b + h − 2e
∂x 2

Δf : d + f + b + h − 4e

2

Le Laplacien de l’image est donc approximé en convoluant l’image
⎛ 0 −1 0 ⎞
⎛− 1 − 1



L1 = ⎜ − 1 4 − 1 ⎟ . En 8-connexité, on utilise le masque : L2 = ⎜ − 1 8
⎜ 0 −1 0 ⎟
⎜− 1 − 1




avec le noyau
− 1⎞

− 1⎟
− 1 ⎟⎠

Ces masques sont caractérisés par un coefficient central positif, des coefficients
périphériques négatifs et une somme des coefficients nulle. L’image obtenue peut prendre des
valeurs positives ou négatives. On doit donc la calibrer pour obtenir une visualisation
correcte. Dans ce cas, le zéro correspond à un gris moyen (FIG. IV-10).
Le problème est que l’opérateur Laplacien est très sensible au bruit, plus que la
dérivée première. C’est pourquoi on prétraite généralement l’image en la filtrant par un filtre
gaussien (paragraphe 2.3.2).
D’autre part, le Laplacien ne donne aucune indication sur la direction du contour.
C’est pourquoi on combine souvent gradient et Laplacien, le Laplacien permettant de savoir si
un pixel fait partie du côté le plus clair ou le plus sombre du contour (FIG. IV-6, FIG.
IV-10(c)). Soit f(x,y) l’image source dont on calcule le gradient et le Laplacien. L’image s(x,y)
déduite de :
⎧0 si ∇ f < T

s( x , y ) = ⎨ A si ∇ f ≥ T et ∇ ² f ≥ 0
⎪ B si ∇ f ≥ T et ∇ ² f < 0


Eq. IV-13

donne une valeur A pour les pixels de contour du côté clair, B pour les pixels du côté sombre,

08/03/2012

- 17 -

Traitement et Analyse des Images Numériques

0 pour les pixels n’appartenant pas au contour.

(a) Image source (cameraman)

(b) Laplacien (L1) après calibration : le zéro
correspond au gris moyen

(c) Zoom : on observe le changement de signe de la dérivée
FIG. IV-10. Extraction de contours par la méthode du Laplacien

Enfin rappelons que les zéros du Laplacien de l’image ne coïncident pas exactement
avec les zéros de la dérivée seconde directionnelle, pour les contours de forte courbure.
2.3.2. Opérateur LOG : approche Marr-Hildreth
Pour obtenir des résultats plus robustes au bruit, on peut préalablement filtrer l’image
par un filtre passe-bas. Le filtre le plus couramment utilisé est un filtre gaussien, dont on
ajuste l’écart-type σ en fonction des détails que l’on veut conserver (plus l’écart-type est
grand, plus l’image est lissée). Le filtrage global (gaussien+Laplacien) d’une image f s’écrit
∇² ( f * g ) , ce qui peut se transformer, grâce aux propriétés de l’opérateur de convolution en :
∇² ( f * g ) = ∇² ( f )* g = f * ∇² (g )

Eq. IV-14

Cela revient donc à filtrer en une seule fois l’image par le Laplacien d’une gaussienne, que
nous allons expliciter.

08/03/2012

- 18 -

Chapitre IV : Segmentation – Extraction de contours

1
Soit g ( x , y ) =
e
2π σ

la fonction h( x , y ) = e



x2 + y 2
2σ 2

(x
∇ g (x , y ) = K
2

2

x2 + y2
2σ 2

le filtre gaussien d’écart-type σ. Le calcul du Laplacien de

nous donne ∇ 2 h(x , y ) =

+ y 2 − 2σ 2

σ

4

)e



(x

+ y 2 − 2σ 2 )

2

σ4

h( x , y ) d’où :

x2 + y2

Eq. IV-15

2σ 2

avec K un facteur de proportionnalité. Ce filtre est appelé LoG (Laplacian of Gaussian) ou
« Chapeau mexicain » (FIG. IV-11). Remarquons que ∇ 2 g ( x , y ) = 0 pour x 2 + y 2 = 2σ 2 .
L’écart-type définit donc bien la résolution du filtre : les détails de taille inférieure à cette
quantité ne seront pas détectés.
En pratique, les coefficients du noyau LoG sont obtenus en discrétisant les variables x
et y. D’autre part, on souhaite que la somme des coefficients soit nulle, afin d’obtenir une
réponse nulle sur les régions uniformes de l’image. Supposons un filtre de taille N=2n+1, on
peut définir le masque de convolution g de l’opérateur LoG par :
1. g0 (x , y ) =

e
n



x2 + y2
2σ 2

n

∑ ∑e



i2 + j2


2

i =−n j =−n

2. g1 ( x , y ) =

(x

2

+ y 2 − 2σ 2

σ

4

⎧ x ∈ {− n ,..., n}
pour ⎨
⎩ y ∈ {− n ,..., n}

(filtre gaussien normalisé)

) g (x , y )

(calcul du Laplacien)

0

⎛ n n
2⎞
3. g log ( x , y ) = g1 ( x , y ) − ⎜⎜ ∑ ∑ g1 (i , j ) / (2 n + 1) ⎟⎟
⎝ i=−n j =−n


(somme des coefficients nulle)

La figure suivante montre la réponse continue de l’opérateur pour σ=0.5 (on a pris le
négatif pour une meilleure visualisation) et le masque de convolution discret glog pour n=2
(N=2n+1=5) :

0.0448
0.0468
0.0564
0.0468
0.0448

Chapeau mexicain

0.0468
0.3167
0.7146
0.3167
0.0468

0.0564
0.7146
-4.9048
0.7146
0.0564

0.0468
0.3167
0.7146
0.3167
0.0468

0.0448
0.0468
0.0564
0.0468
0.0448

Opérateur LoG pour N=5 et sigma = 0.5

FIG. IV-11. Opérateur LoG (Laplacian of Gaussian) ou chapeau mexicain

Une fois l’opérateur LoG appliqué sur l’image, l’étape suivante consiste à rechercher
les passages par zéro dans l’image résultat. Pour cela, on compare le signe de chaque pixel
08/03/2012

- 19 -

Traitement et Analyse des Images Numériques

avec celui de son voisinage. Notons iml l’image filtrée par le noyau LoG, S un seuil défini par
l’utilisateur. Les pixels marqués comme points de contour sont ceux qui satisfont à l’une des
conditions suivantes :
iml( i , j ) < 0 & iml( i , j + 1 ) > 0 & iml( i , j ) − iml( i , j + 1 ) > S

iml( i , j ) < 0 & iml( i , j − 1 ) > 0 & iml( i , j ) − iml( i , j − 1 ) > S
iml( i , j ) < 0 & iml( i + 1, j ) < 0 & iml( i , j ) − iml( i + 1, j ) > S
iml( i , j ) < 0 & iml( i − 1, j ) < 0 & iml( i , j ) − iml( i − 1, j ) > S
iml( i , j ) == 0 & iml( i − 1, j ) < 0 & iml( i + 1, j ) > 0 & iml( i − 1, j ) − iml( i + 1, j ) > 2 S
iml( i , j ) == 0 & iml( i − 1, j ) > 0 & iml( i + 1, j ) < 0 & iml( i − 1, j ) − iml( i + 1, j ) > 2 S
iml( i , j ) == 0 & iml( i , j − 1 ) < 0 & iml( i , j + 1 ) > 0 & iml( i , j − 1 ) − iml( i , j + 1 ) > 2 S
iml( i , j ) == 0 & iml( i , j − 1 ) > 0 & iml( i , j + 1 ) < 0 & iml( i , j − 1 ) − iml( i , j + 1 ) > 2 S

Appliquons maintenant l’opérateur LoG sur l’image « cameraman » (FIG. IV-10(a))
pour différentes valeurs de σ . La figure FIG. IV-12 illustre les résultats. A titre indicatif, on
montre également les images filtrées par un noyau gaussien. On remarque l’influence du
paramètre σ qui détermine le niveau de détail détectable.
2.3.3. Opérateur DoG
Cet opérateur est une approximation de l’opérateur LoG par la différence de deux
gaussiennes dont les écarts-types sont dans un rapport 1.6.
g DOG ( x , y ) =

1
2πσ 1

2

e



x2 + y2
2σ 1

2



1
2πσ 2

2

e



x2 + y2
2σ 2 2

avec

σ1
= 1.6
σ2

Eq. IV-16

Les opérateurs LoG et DoG ont les mêmes passages par zéros si on respecte la relation
entre les écarts-types donnée par l’équation Eq. IV-17. Une fois l’amplitude en zéro
normalisée, les profils des deux opérateurs sont alors presque confondus.

σ2 =

σ 1 ²σ 2 ² ⎛ σ 1 ² ⎞

ln⎜
σ 1 ² − σ 2 ² ⎜⎝ σ 2 ² ⎟⎠

Eq. IV-17

Dans les deux cas, les filtres sont séparables. On peut donc implémenter le filtrage par
2 convolutions 1D au lieu d’une convolution 2D, ce qui réduit le coût de calcul.

2.4.

Discussion

Nous avons donc vu différentes méthodes agissant par convolution spatiale. La mise
en œuvre est très simple et la complexité de calcul est faible. Cependant, différents problèmes
se posent :
1. le bruit est interprétable comme un contour (point de transition, haute fréquence), d’où
de fausses détections. Pour pallier ce problème, certains opérateurs ont été conçus
comme la combinaison d’un filtre passe-bas et d’un filtre passe-haut, le filtre passe-

08/03/2012

- 20 -

Chapitre IV : Segmentation – Extraction de contours

bas servant à réduire le bruit. On peut aussi prétraiter l’image par un filtre de lissage
non linéaire (filtre médian, filtre de diffusion non linéaire, …) avant d’appliquer
l’opérateur de détection de contour. Enfin, il peut être avantageux de combiner
plusieurs méthodes d’extraction de contour pour obtenir un résultat final moins
sensible au bruit.
2. les contours ne sont pas fermés (points manquants), ni uniques (épaisseur > 1 pixel).
Pour améliorer les résultats, on post-traite l’image de gradient en supprimant les pixels
qui ne sont pas des maxima locaux. On utilise aussi le seuillage par hystérésis qui donne de
meilleurs résultats que le seuillage simple. Enfin, des procédures de prolongement de contours
peuvent être appliqué pour tenter de les fermer. Ces méthodes seront exposées dans la section
4 de ce chapitre.

Σ=0.5 ;N = 5, S = 0.3

σ=2, N = 13 ; S = 0.01

σ=5, N = 31, S = 0.0005

FIG. IV-12. Application de l’opérateur LOG pour différentes valeurs de σ. Première ligne : image filtrée
par un filtre gaussien, deuxième ligne : image filtrée par l’opérateur LoG, 3ème ligne : passages par 0.

08/03/2012

- 21 -

Traitement et Analyse des Images Numériques

3.

Approches analytiques

3.1.

Filtre de Canny

Canny a proposé en 1986 une formalisation du problème de la détection d’un contour
[11]. Il se place dans le cas 1D, et modélise le contour par un signal A(x) (par exemple un
échelon) bruité par un bruit blanc gaussien. La démarche consiste à déterminer la réponse
impulsionnelle h(x) d’un filtre satisfaisant trois critères :
• Bonne détection
• Bonne localisation
• Réponse unique
Ces trois critères ont été mathématiquement formulés comme suit. Les détails des
calculs sont donnés dans l’annexe, section 2.
3.1.1. Modélisation du contour
Le contour est modélisé par un signal A(x) monodimensionnel, somme d’un échelon
U(x) d’amplitude U0 et d’un bruit gaussien N(x) de moyenne nulle et de densité spectrale de
2

puissance constante N0 :

FIG. IV-13. Modélisation du contour (1D) par un échelon bruité par un bruit gaussien

A(x ) = U 0U ( x ) + N ( x ) , avec

⎧1 si x ≥ 0
U (x ) = ⎨
⎩0 sinon

Eq. IV-18

On recherche la réponse impulsionnelle h(x) d’un filtre qui permet de détecter ce
contour de manière optimale au sens des trois critères de Canny. La réponse de ce filtre au
signal d’entrée A(x) s’écrit :
+∞

C ( x ) = ( A * h )( x ) = ∫ A(u )h( x − u )du

Eq. IV-19

−∞

3.1.2. Formalisation du critère 1 : bonne détection
On aura une bonne détection du contour en maximisant le rapport signal à bruit, c’est-

08/03/2012

- 22 -

Chapitre IV : Segmentation – Extraction de contours

à-dire le rapport entre la réponse due au signal non bruité (l’échelon) et la réponse due au
bruit seul (racine carrée de la puissance du bruit en sortie du filtre).
+∞

SNR =

U 0 ∫ h( x − u )du
0



N0

Eq. IV-20

+∞

−∞

h 2 (u )du

Cette réponse doit être maximale en 0, la coordonnée du contour. Le premier critère de Canny
consiste donc à maximiser le rapport
+∞

U 0 ∫ h(− u )du
U0
0
Σ=
+∞
N0
N0 ∫ h 2 (u )du

Eq. IV-21

−∞

3.1.3. Formalisation du critère 2 : bonne localisation
Dans le cas idéal, c’est-à-dire sans bruit, le maximum de la réponse du filtre doit se
trouver en x=0. A cause du bruit, ce maximum peut se déplacer en x0 ≠ 0 . Le filtre aura de
bonnes propriétés de localisation si x0 est proche de 0. Pour évaluer le critère de localisation,
2
on est donc amené à mesurer la variance de x0, soit E x0 .

{ }

{ }
2

E x0 ≈

N0

2



+∞

−∞
2
0

h' 2 (u )du

Eq. IV-22

U h' 2 (0 )

Plus la variance est faible, meilleure est la localisation. On prendra donc comme
critère de bonne localisation la racine carrée de l’inverse de l’expression Eq. IV-22.

U 0 h' (0 )
U0
Λ=
+∞
N0
N0 ∫ h' 2 (u )du

Eq. IV-23

−∞

Le produit ΣΛ (Eq. IV-21, Eq. IV-23) combine une bonne détection et une bonne
localisation :
+∞

ΣΛ =

∫ h(− u )du
∫ h (u )du ∫

h' (0 )

0

+∞

−∞

2

+∞

−∞

Eq. IV-24

h' 2 (u )du

On remarque que ce produit ne dépend pas de l’amplitude de l’échelon ni de la densité
spectrale de puissance du bruit. Il ne dépend pas non plus du facteur d’échelle. En effet, si on
Λ
⎛x⎞
pose hκ (x ) = h⎜ ⎟ , on obtient Σ κ = κ Σ et Λκ =
. En d’autres termes, si la réponse
κ
⎝κ ⎠
impulsionnelle du filtre s’étend, le rapport signal à bruit s’améliore au détriment de la
localisation, et inversement. Il n’est pas possible d’optimiser simultanément les deux critères.

08/03/2012

- 23 -

Traitement et Analyse des Images Numériques

3.1.4. Formalisation du critère 3 : unicité de la réponse
En présence de bruit, la sortie du filtre peut présenter plusieurs maxima locaux. Si on
prend la dérivée de la sortie, on peut donc avoir plusieurs passages par zéros.
Canny a choisi de déterminer un filtre h(x) à réponse impulsionnelle finie, de taille
2M. La réponse du filtre est donc concentrée sur une fenêtre de largeur 2M. Il impose ensuite,
comme dernière contrainte, que la distance xmax entre deux maxima locaux consécutifs de la
réponse due au bruit soit égale à une fraction de M.
1

xmax = 2 xmoy

⎛ + M h' (u )2 du ⎞ 2

⎜∫
= 2π ⎜ −+MM
⎟ = km .M
2
⎜ ∫ h' ' (u ) du ⎟

⎝ −M

Eq. IV-25

3.1.5. Résolution
Canny propose donc de minimiser le produit ΛΣ (Eq. IV-24) sous la contrainte
xmax = km .M (Eq. IV-25). Nous ne donnerons pas ici tous les développements mathématiques
(très complexes [11]…). Nous indiquons seulement les résultats principaux.
• La minimisation débouche sur une équation différentielle dont la solution est de la
forme :

h( x ) = a1eαx sin ωx + a2 eαx cos ωx + a3e −αx sin ωx + a4 e −αx cos ωx






Eq. IV-26

On recherche un filtre de réponse impulsionnelle finie, définie sur l’intervalle
[-M,M]. Canny impose 4 conditions aux limites pour h :
h(0)=0, h(M)=0, h’(0)=S, h’(M)=0)
Canny peut ainsi trouver une expression analytique définissant les paramètres a1 à
a4. De plus, h étant impaire, la solution peut être étendue aux x négatif avec
h(-x)=- h(x).
Par optimisation numérique, Canny montre ensuite que le filtre h le plus
performant correspond à l’indice de performance ΛΣ = 1.12 , obtenu pour une
combinaison particulière des constantes définissant a1 à a4.
Le filtrage étant complexe à mettre en œuvre, Canny propose de l’approximer par
l’opérateur dérivée première de la Gaussienne, pour lequel ΛΣ = 0.92 et k=0.51
(dégradation des performances de 20%).

h( x ) ∝ − xe

⎛ x² ⎞
⎜−

⎝ 2σ ² ⎠

Eq. IV-27

On retrouve donc l’opérateur de Marr-Hildreth (§2.3.2), mais appliqué au niveau de la
dérivée première…
La figure FIG. IV-14 représente la réponse de ce filtre à un échelon pour différentes
valeurs de σ. Le paramètre σ est un paramètre d’échelle qui indique la distance à partir de

08/03/2012

- 24 -

Chapitre IV : Segmentation – Extraction de contours

1.5

1

0.5

(a)

0

-0.5

-1

-1.5

0

20

40

60

80

100

120

140

160

180

200

s igma = 0.50
3

2

1

(b)

0

-1

-2

-3

0

20

40

60

80

100

120

140

160

180

200

140

160

180

200

140

160

180

200

140

160

180

200

s igma = 1.00
2.5
2
1.5
1
0.5

(c)

0
-0.5
-1
-1.5
-2
-2.5

0

20

40

60

80

100

120

s igma = 2.00
2.5
2
1.5
1
0.5

(d)

0
-0.5
-1
-1.5
-2
-2.5

0

20

40

60

80

100

120

s igma = 5.00
2.5
2
1.5
1

(e)

0.5
0
-0.5
-1
-1.5
-2

0

20

40

60

80

100

120

FIG. IV-14. Filtre de Canny appliqué à une marche d’escalier pour plusieurs écarts-types σ.
En bleu le signal d’entrée, en vert le signal lissé, en rouge la dérivée du signal lissé. (a) : signal non bruité, σ =
1, (b)(c)(d)(e) Signal bruité par un bruit gaussien (écart-type 0.4) pour (b) σ = 0.5, (c) σ = 1,(d) σ = 2, (e) σ = 5.

08/03/2012

- 25 -

Traitement et Analyse des Images Numériques

laquelle deux contours proches peuvent être séparés. On voit qu’il y a une sorte de compromis
à faire : une grande valeur de σ lisse bien le bruit et permet une bonne détection de contour.
En revanche deux contours proches, d’une distance inférieure à σ ne pourront être séparés.
3.1.6. Mise en œuvre du filtre de Canny
Le filtre précédent est un filtre 1D, qu’il faudrait théoriquement appliquer dans la
direction orthogonale au contour de l’image. La direction du contour n’étant pas connue a
priori, il faudrait donc en principe appliquer le filtre dans toutes les directions.
En pratique, on préfère réaliser l’approximation suivante :
• Lisser l’image dans toutes les directions par un filtre Gaussien G(x,y) (isotrope)
(Eq. IV-28)
• Calculer le gradient suivant les deux axes, horizontal et vertical, par les techniques
classiques déjà vues au paragraphe 2.2. Cette opération est effectuée en convoluant
l’image lissée par les dérivées suivant x et y du noyau gaussien G(x,y) (Eq. IV-29).
• Calculer le module (Eq. IV-3) et l’argument (Eq. IV-4) du gradient pour détecter
les points de contour et estimer la direction du contour.
G(x , y ) ∝ e



x2 + y2

Eq. IV-28

2σ 2

∂(G ( x , y )* I ( x , y )) ∂ (G (x , y ))
=
* I (x , y )
∂x
∂x

Eq. IV-29

Le filtre gaussien 2D est séparable : on peut donc appliquer successivement 2 filtres
1D, le premier sur les lignes de l’image I, le second sur les colonnes de l’image en sortie du
premier filtre. Cette implémentation réduit le coût de calcul.
Naturellement les contours obtenus seront plus ou moins épais. Pour améliorer les
résultats, on applique une procédure d’amincissement (suppression de pixels qui ne sont pas
des maxima locaux dans la direction du gradient, voir le paragraphe 4.1 de ce chapitre).
La dernière étape consiste à détecter les contours en appliquant un seuillage par
hystérésis (voir paragraphe 4.1 de ce chapitre).

3.2.

Filtre de Deriche

Deriche a repris les critères de Canny, mais il a cherché un filtre à réponse
impulsionnelle infinie [12]. Il a ensuite proposé une implémentation de ce filtre pour le
traitement des images.
3.2.1. Filtre mono-dimensionnel
La résolution des critères de Canny pour un filtre à réponse impulsionnelle infinie
aboutit à la fonction :

08/03/2012

- 26 -

Chapitre IV : Segmentation – Extraction de contours

h( x ) =

S

ω

e

−α x

sin(ωx )

Eq. IV-30

Deriche a ensuite évalué ce filtre en calculant les critères de performance définis par
Canny, pour différents paramétrages :

Λ = 2α

Σ=


2
α + ω2

k=

α 2 + ω2
5α 2 + ω 4 + 6α 2ω 2

Eq. IV-31

Il démontre ainsi la supériorité de son filtre par rapport au filtre de Canny : indice de
performance ΛΣ supérieur pour une même valeur de k. En particulier, les meilleures
performances sont obtenues en faisant tendre ω vers 0. Ce cas limite correspond à :
h( x ) = cxe

−α x

Eq. IV-32

avec c une constante de normalisation : c = −α 2 . Dans ce cas l’indice de performance est égal
à ΛΣ = 2 pour k=0.44.
3.2.2. Filtre bi-dimensionnel
L’idée est de mettre en œuvre les techniques usuelles de détection de contour en
utilisant le filtre optimal de Deriche : lissage pour enlever du bruit, calcul du gradient pour
mettre en évidence les contours.
Pour le lissage, Deriche propose de prendre l’intégrale du filtre optimal :

f ( x ) = b(α x + 1) e

−α x

Eq. IV-33

b étant une constante calculée pour donner une réponse 1 à un signal constant valant 1 :
b =α / 4.
Le filtre séparable de lissage s’écrit :
f ( x , y ) = b² (α x + 1) e

−α x

(α y + 1) e α


y

Eq. IV-34

L’image A lissée s’écrit
B(x , y ) = A * f (x , y )

Eq. IV-35

On calcule ensuite les dérivées suivant x et y de l’image lissée B ( x , y ) :

∂B
(x , y ) = ∂( A * f ) (x , y ) = A * ∂f (x , y )
∂x
∂x
∂x

Eq. IV-36

∂f
(x , y ) = ηxe−α x (α y + 1)e−α y
∂x

Eq. IV-37

Or :

avec η une constante de normalisation pour fournir une réponse maximale égale à 1 à un

08/03/2012

- 27 -

Traitement et Analyse des Images Numériques

échelon d’amplitude 1. η = −α 3 / 4 . De Eq. IV-32, Eq. IV-33 et Eq. IV-37 on déduit :

∂f
( x , y ) = h( x ) f ( y )
∂x

Eq. IV-38

En remplaçant dans Eq. IV-36 :

∂B
(x , y ) = A * h(x )* f ( y ) = A * (h(x ). f ( y ))
∂x

Eq. IV-39

Ce résultat est très important. Il montre que la dérivée suivant x de l’image lissée peut
être obtenue en lissant l’image suivant y puis en dérivant suivant x. Et réciproquement pour
obtenir la dérivée suivant y de l’image lissée. La suite de la démarche pour l’extraction de
contour est de calculer classiquement la norme et l’argument du gradient.
3.2.3. Implémentation récursive

On peut montrer, grâce à la transformée en Z, que le filtre intégrateur (Eq. IV-33) et le
filtre dérivateur (Eq. IV-32) peuvent s’implémenter sous la forme d’un filtre récursif d’ordre
2, et ceci quel que soit la valeur du paramètre α .
Le filtre dérivateur selon x (Eq. IV-39), discrétisé a pour expression :

SS x (m ,n ) = kme

−α m

k (α n + 1)e

−α n

Eq. IV-40

De même pour le filtre dérivateur selon y :

SS y (m ,n ) = k (α m + 1)e

−α m

kne

−α n

Eq. IV-41

Chacun de ces deux filtres peut être implémenté de façon récursive, en deux phases.
La première traite les lignes de l’image et la seconde des colonnes, ce qu’il est possible de
faire puisque les filtres sont séparables. Pour chacune des phases, on filtre l’image par deux
filtres récursifs d’ordre 2, le premier suivant les indices croissants, le second suivant les
indices décroissants, puis on réunit les deux résultats y1 et y2. Nous présentons ci-dessous
l’algorithme :
Soit H et W le nombre de lignes et le nombre de colonnes de l’image A.
Phase 1 : filtrage de chaque ligne de l’image source A…
- Pour m variant de 1 à H
% …traiter chaque ligne
- Pour n variant de 1 à W :
% … par index de colonne croissant
y1 (m ,n ) = a1 A(m ,n ) + a2 A(m ,n − 1) + b1 y1 (m ,n − 1) + b2 y1 (m , n − 2 )

- Pour n variant de W à 1 :
% … par index de colonne décroissant
y2 (m ,n ) = a3 A(m ,n + 1) + a4 A(m ,n + 2 ) + b1 y2 (m ,n + 1) + b2 y2 (m ,n + 2 )
- Pour n variant de 1 à W :
% réunion des résultats
R (m ,n ) = c1 ( y1 (m , n ) + y2 (m ,n ))

08/03/2012

- 28 -

Chapitre IV : Segmentation – Extraction de contours

Phase 2 : filtrage de chaque colonne de l’image résultat R de la phase 1…
- Pour n variant de 1 à W
% … par index de ligne croissant
- Pour m variant de 1 à H :
y1 (m ,n ) = a5 R(m ,n ) + a6 R(m − 1,n ) + b1 y1 (m − 1,n ) + b2 y1 (m − 2 ,n )

- Pour m variant de H à 1 :
% … par index de ligne décroissant
y2 (m ,n ) = a7 R(m + 1,n ) + a8 A(m + 2 ,n ) + b1 y2 (m + 1,n ) + b2 y2 (m + 2 ,n )
- Pour n variant de 1 à H :
F (m , n ) = c2 ( y1 (m , n ) + y2 (m ,n ))

% réunion des résultats

Les valeurs non définies (index négatifs ou dépassant la taille des images), on prendra 0.
L’image F est le résultat du filtrage.
Le tableau Tab. IV.1 indique les coefficients utilisés pour chacun des filtres : lissage
dans les deux directions, dérivée selon x, dérivée selon y.
(a) Lissage

(1 − e )

−α 2

K

1 + 2αe −α − e − 2α
a1
a2

k

a5

(1 − e )

−α 2

1 + 2αe −α − e − 2α

k

0

ke (α − 1)
ke −α (α + 1)
− ke−2α

k

0

−α

ke (α − 1)
ke −α (α + 1)
− ke−2α
−α

(1 − e )

1 + 2αe −α − e − 2α

−α

a6
a7
a8

(c) Dérivée selon y
−α 2

ke (α − 1)
ke −α (α + 1)
− ke−2α
−α

a3
a4

(b) Dérivée selon x

1
-1
0
k

0

ke (α − 1)
ke −α (α + 1)
− ke−2α
2 e −α
− e −2α

−α

1
-1

b1
b2
c1

2e
− e −2α

2 e −α
− e −2α

1

1

c2

1

− (1 − e −α )

− (1 − e −α )

2

2

1

Tab. IV.1 : coefficients des filtres récursifs de Deriche (1) lissage (2)(3) dérivation

L’implémentation récursive privilégie la précision numérique. D’autre part, le nombre
de calculs effectués par pixel est indépendant de la taille du filtre (spécifiée par α ).
-

La méthode générale de détection des contours est la suivante :
Déterminer la dérivée directionnelle selon x : lissage suivant la direction y puis
dérivation suivant x.
Déterminer la dérivée directionnelle selon y : lissage suivant la direction x puis
dérivation suivant y
Calculer la norme du gradient en fonction des dérivées directionnelles (Eq. IV-3).
Extraire les maxima locaux dans la direction du gradient (§4.1).
Seuillage par hystérésis (§4.2)

08/03/2012

- 29 -

Traitement et Analyse des Images Numériques

3.2.4. Exemples

La figure Fig. IV-15 montre des exemples de lissage (coefficients donnés dans la
colonne (1) du tableau Tab. IV.1), et la figure Fig. IV-16 des exemples d’extraction de
contours en appliquant les filtres récursifs de Deriche (colonnes (2) et (3) du tableau Tab.
IV.1) et en calculant le module de la dérivée suivant Eq. IV-3.

Image originale

α=2

α=1

α=0.75

α=0.5

α=0.25

Fig. IV-15 : Filtre de lissage de Deriche, implémentation récursive pour plusieurs valeurs de α

Module, α=3

Module, α=2

Module, α=0.5

Fig. IV-16: Filtre dérivateur de Deriche, implémentation récursive pour plusieurs valeurs de α

08/03/2012

- 30 -

Chapitre IV : Segmentation – Extraction de contours

On voit l’effet du paramètre α qui agit dans le sens inverse du paramètre σ du filtre de
Canny : α = π / σ . Il agit donc comme un paramètre d’échelle : plus α est faible, plus
l’image est lissée (Fig. IV-15) et plus les détails sont ignorés (Fig. IV-16). Dans les images de
contours (Fig. IV-16), plus α est faible, plus les contours sont épais. Ils sont donc bien
détectables mais aussi moins bien localisés. Si on revient aux équations Eq. IV-31, diminuer α
fait augmenter Σ (meilleure détection) mais diminuer Λ (moins bonne localisation).

4.

Extraction des points de contour et post-traitements

Nous avons décrit des méthodes fondées sur l’application de filtres, permettant de
mettre en évidence les contours : maxima de la dérivée première ou passage par zéro de la
dérivée seconde.
L’étape suivante consiste à classer les pixels en « contour » ou « non contour », le
résultat de cette opération étant une image binaire dans laquelle les contours sont à 1.
Nous avons déjà évoqué la méthode la plus intuitive, le seuillage. Cette méthode
simple ne permet clairement pas d’opérer une classification correcte, car on se heurte à des
problèmes de fausses détections, de contours mal localisés car épais, de points de contours
non détectés, de contours non fermés.
Dans la suite nous allons présenter trois techniques qui permettent d’obtenir des
résultats plus fiables et plus précis : la suppression des non maxima, le seuillage par
hystérésis, la fermeture des contours.

4.1.

Suppression des non-maxima

Cet algorithme détecte les crêtes de l’image du module du gradient en mettant à zéro
les points qui ne sont pas des maxima locaux dans la direction perpendiculaire à la crête, cette
dernière étant indiquée par l’image de l’argument du gradient. Concrètement il s’agit donc de
comparer le module du gradient au point P(x,y) avec le module du gradient de ses deux
voisins P1 et P2, ces deux points étant situés dans la direction du gradient, à une distance
unitaire du point P (FIG. IV-17(a)).
L’image étant discrète et le maillage généralement un carré, on définit 4 directions
principales dk : − 45° ± 22.5°,0° ± 22.5°,45° ± 22.5°,90° ± 22.5° (FIG. IV-17(b)).

(a)

(b)

FIG. IV-17. (a) comparaison du module du gradient au voisinage de P dans la direction du gradient ; (b)
quantification de la direction du gradient sur 4 directions principales

08/03/2012

- 31 -

Traitement et Analyse des Images Numériques

Soient M ( x , y ) le module du gradient au point P(x,y) et α (x , y ) l’argument du gradient. En
chaque pixel P, on applique l’algorithme suivant :
• Déterminer la direction dk la plus proche de α (x , y )



Si le pixel (x,y) a au moins un voisin dans la direction dk de module M supérieur à
M ( x , y ) , alors M ' ( x , y ) = 0 , sinon M ' (x , y ) = M ( x , y )
L’image M ' ( x , y ) représente donc le module du gradient après mise à zéro des pixels qui ne
sont pas des maxima locaux dans la direction du gradient.
La figure FIG. IV-18 illustre le procédé sur un exemple très simple. La figure FIG.
IV-20 montre un autre résultat sur une image plus complexe.

Image source

Module du gradient

M (x , y )

Argument du gradient
quantifié sur 4 directions
principales

Module après
suppression des nonmaxima

FIG. IV-18. Suppression des pixels qui ne sont pas des maxima locaux dans la direction du gradient

Pour être plus précis, on peut éviter de quantifier les directions. Dans ce cas, il faut
calculer les modules du gradient des deux voisins de P par interpolation linéaire.

4.2.

Seuillage par hystérésis

Nous avons vu au paragraphe 2.2, FIG. IV-9 que les points de contour peuvent être
extraits par seuillage de l’image de gradient. Les résultats sont très dépendants du seuil choisi,
et un seul seuil ne suffit généralement pas à discriminer les pixels du contour et d’autres
pixels de gradient assez élevé correspondant à des transitions de niveau de gris non
recherchées (dues à du bruit par exemple).
Afin d’obtenir de meilleurs résultats, on peut appliquer un seuillage à hystérésis. Le
principe est d’appliquer deux seuils. Le premier seuil, bas, est choisi pour que tous les points
de contour soient détectés. Le second seuil, haut, est choisi pour que seuls des points de
contour soient détectés. Ces derniers sont étiquetés « point de contour ». Tous les pixels dont
le module du gradient est compris entre ces deux seuils, et qui possèdent dans leur voisinage
un point de contour dans la direction du gradient, sont également étiquetés comme points de
contours. L’algorithme est itéré tant que de nouveaux pixels peuvent être étiquetés « point de
contour ».
La figure FIG. IV-19 montre un exemple. Cette approche donne de meilleurs résultats
qu’un seuillage simple. L’automatisation des seuils reste cependant très délicate. Cette
méthode peut fonctionner dans le cas d’images assez simples. Pour des images naturelles

08/03/2012

- 32 -

Chapitre IV : Segmentation – Extraction de contours

Image source

Module du gradient

Direction du gradient

Module gradient
seuillé S1 = 0.4

Module gradient
seuillé S2 = 0.95

Seuillage par hystérésis
Différence angulaire
max du gradient: π/20

Seuillage par hystérésis
Différence angulaire
max du gradient: π/4

Module gradient
seuillé S1 = 0.55

Module gradient
seuillé S2 = 0.99

Seuillage par hystérésis
Différence angulaire
max du gradient: π/20

Seuillage par hystérésis
Différence angulaire
max du gradient: π/4

FIG. IV-19. Seuillage d’une image de gradient par hystérésis (Les images sont calibrées entre 0 et 1)

Module du gradient

Module du gradient après
suppression des non-maxima
locaux

Contours après seuillage par
hystérésis

FIG. IV-20. Extraction du contour, par calcul du gradient local et post-traitements

08/03/2012

- 33 -

Traitement et Analyse des Images Numériques

complexes, on obtient rarement un partitionnement de l’image en objets cohérents. La figure
FIG. IV-20 montre un exemple qui combine algorithme de suppression des non-maxima du
module du gradient et le seuillage par hystérésis.

4.3.

Fermeture des contours

Le but est de créer des régions connexes à partir des contours extraits d'une image de
niveaux de gris. Nous présentons ici 2 algorithmes réalisant cette opération.
4.3.1. Algorithme 1: recherche locale, critère de ressemblance

Soient :

∇f (x,y) : image du module du gradient

α (x,y)
c (x,y)

: image de la direction du gradient
: image binaire des contours c (x,y)
= 0 sinon.

= 1 si c∈C

L’algorithme suivant permet d’inclure au contour C de nouveaux points de contour,
sur critère local de ressemblance au voisinage des extrémités de contour, et ainsi de fermer
certains d’entre eux.
1. Pour chaque extrémité de contour : c ∈ C (i.e. c(x,y) = 1 )
2. Trouver dans un voisinage 3 x 3 de c les points de coordonnées (x’,y’) qui ne sont pas
des points de contours : c’(x’,y’)=0
3. Retenir, s’il existe, le point c’ de coordonnées (x’,y’) qui satisfait au mieux les
conditions suivantes : |∇ f (x,y) -∇ f (x',y')|≤ S (1)
|α (x,y) -α (x',y')|< A
(2)
4. Inclure c’ à C : c(x’,y’) = 1
5. Retourner en 1.
4.3.2. Algorithme 2: gradient moyen maximal

L’algorithme précédent est très local, et on lui préfèrera généralement celui-ci. Le
principe est d’élaborer un chemin de crête à partir de chaque extrémité de contour, afin de
réaliser une analyse arborescente sur une profondeur P (FIG. IV-21). On calcule sur chaque
chemin une fonction de coût, par exemple la moyenne des normes des gradients, et on retient
le chemin qui maximise le coût. Le contour est prolongé suivant le chemin trouvé. On itère à
partir des nouvelles extrémités tant qu’on ne rencontre pas de points de contours ou que l’on
n’atteint pas le critère d’arrêt (coût inférieur à un seuil par exemple).
1. Pour chaque extrémité Ei du contour C de l'image binaire
2. Création des chemins Ei (n) sur une profondeur P
3. Recherche du chemin dont le gradient moyen est maximal et supérieur à Gseuil :
Gmoy > Gseuil. Soit Ej la nouvelle extrémité.
• Si Gmoy ≤ Gseuil , aucune prolongation, retourner en 1.
• Si Gmoy > Gseuil et Ej ∉ C, Ej est une nouvelle extrémité. Retourner en 2.
• Si Gmoy > Gseuil et Ej ∈ C, Ej le coutour est fermé. Retourner en 1.

08/03/2012

- 34 -

Chapitre IV : Segmentation – Extraction de contours

FIG. IV-21. Analyse arborescente pour la recherche du gradient moyen maximal (P=2)

Contours ajoutés par l’algorithme de
fermeture (rouge)
Contours initiaux (blanc)

FIG. IV-22. Fermeture de contours

5.

Contours actifs

Les méthodes proposées initialement pour la détection des contours reposent
essentiellement sur une analyse très locale de l’image. On recherche les pixels de fort
gradient, puis on applique des post-traitements permettant d’extraire les points de contour
« sûrs », d’amincir les morceaux de contours obtenus et de combler les contours manquants
par des algorithmes de fermeture. Sans ajout d’informations a priori, il est assez difficile
d’obtenir des résultats fiables, puisqu’il est à peu près impossible de faire la distinction entre
les pixels de fort gradient correspondant aux contours recherchés, et les pixels de fort gradient
correspondant à du bruit ou autre transition de niveau de gris « parasite ».
A la fin des années 1980, d’autres méthodes ont été proposées, permettant d’analyser
l’image de manière plus globale. Il s’agit des contours actifs, introduits par Kass, Witkin et
Terzopoulos (1988) [14]. [14]. Une courbe, généralement fermée, est initialisée au voisinage
du contour recherché. L’initialisation peut être faite manuellement (extraction de contour
semi-automatique) ou automatiquement, à partir d’une analyse préalable, d’informations a
priori, etc. On applique ensuite à cette courbe des déformations, qui vont l’amener à évoluer
progressivement vers le contour recherché. Ce processus itératif correspond à la minimisation
d’une fonctionnelle d’énergie, que nous allons expliciter dans le paragraphe suivant.

5.1.

Modélisation
Considérons une image I et un contour C de cette image. Ce contour est modélisé par

08/03/2012

- 35 -

Traitement et Analyse des Images Numériques

une fonction de l’abscisse curviligne s:
C = {v(s ) = ( x(s ), y (s )), s ∈ [0 ,1]}

Eq. IV-42

On associe à la courbe une énergie, qui s’écrit, dans le cas général, comme la somme de 3
termes :
E (C ) = ∫ Eint erne (v(s )) + Eexterne (v(s )) + Eimage (v(s ))ds
1

Eq. IV-43

0

L’énergie interne est liée à la courbe uniquement. Kass, Witkin et Terzopoulos [14]
ont choisi la formulation suivante :

∂ 2 v(s )
∂v(s )
+ β (s )
= α (s )
∂s 2
∂s
2

Eint erne

2

Eq. IV-44

L’énergie interne est un terme de régularisation de la courbe. Dans l’équation Eq.
IV-44, le premier terme prend en compte les variations de longueur de la courbe. C’est donc
un terme de tension (résistance à la rupture) contrôlé le paramètre α(s) d’élasticité. Lorsque
α(s)=0, la courbe peut présenter des discontinuités. Le second terme exprime les variations de
courbure. C’est donc un terme de flexion. Le paramètre, β(s), contrôle la capacité de la courbe
s’incurver. Par exemple, si β(s)=0, la courbe peut être discontinue au second ordre et
développer un coin. En général, on prend α et β indépendants de s, mais il peut être intéressant
de les faire varier localement.
L’énergie image est directement calculée à partir de l’image I. Pour que la courbe
converge vers les zones de fort gradient, on peut choisir une énergie liée au gradient de
l’image (Eq. IV-45), ou au gradient de l’image filtrée par un filtre gaussien (Eq. IV-46) :

Eimage = − ∇I (v(s ))

2

Eimage = − ∇(I * gσ )(v(s ))

Eq. IV-45

2

Eq. IV-46

Néanmoins, si on recherche des crêtes dans l’image, l’énergie image peut être tout
simplement l’image d’intensité elle-même, avec un signe qui dépend du type de crête
recherché, clair ou sombre.

Eimage = ± I (v(s ))

Eq. IV-47

En ce qui concerne l’énergie externe, beaucoup de possibilités peuvent être
envisagées, suivant le problème. L’objectif est d’introduire des informations a priori sur les
contours que l’on recherche. Cette énergie est donc aussi appelée énergie de contexte. Si
aucune énergie externe n’entre en jeu, alors Eexterne = 0 .
On recherche la courbe C qui minimise l’énergie du contour (Eq. IV-43). Cette
équation peut être vue comme la somme de deux forces qui s’annulent à l’équilibre : une

08/03/2012

- 36 -

Chapitre IV : Segmentation – Extraction de contours

force externe, correspondant à l’énergie externe et l’énergie de contexte, qui vont déformer le
contour dans la direction voulue (vers les contours cherchés), et une force interne,
correspondant à l’énergie interne, qui a tendance à s’opposer à ce mouvement pour préserver
la régularité de la courbe.

5.2.

Résolution (cas d’un contour fermé, α et β constants)

On recherche la courbe C qui minimise l’énergie du contour (Eq. IV-43). On peut réécrire cette énergie sous une forme plus générale :
∂v(s )
∂ 2 v (s )
E (C ) = ∫ α (s )
+ β (s )
+ P(v(s ))ds
0
∂s
∂s 2
1

2

2

Eq. IV-48

P (v ) représente la somme des énergies image et externe. Par exemple, P(v ) = − ∇I (v ) si on
2

considère une énergie image déduite du gradient image et une énergie externe nulle.
5.2.1. Equation d’Euler

Trouver un minimum de E(C) (Eq. IV-48) revient à résoudre l’équation d’Euler :
− (α (s )v' (s ))' +(β (s )v' ' (s ))' ' +∇P(v ) = 0

Eq. IV-49

Ce qui s’écrit aussi :
− (α (s )v' (s ))' +(β (s )v' ' (s ))' ' = F (v ) avec F (v ) = −∇P(v )

Eq. IV-50

F (v ) peut s’interpréter comme un champ de forces déduites de l’image (et éventuellement de

connaissance a priori), qui poussent la courbe vers le contour recherché. On peut aussi
réécrire l’équation Eq. IV-49 sous la forme [15] :
⎧ F (v ) = −∇P(v )
Fint + Fext = 0 avec ⎨ ext
⎩ Fint (v ) = (α (s )v' (s ))' + (β (s )v' ' (s ))' '

Eq. IV-51

Les forces externes poussent la courbe vers le contour recherché, tandis que les forces internes
s’opposent à la déformation de la courbe. La solution correspond à un état d’équilibre entre
les deux.
L’équation Eq. IV-50 peut avoir plusieurs solutions puisque l’énergie n’est
généralement pas convexe et peut donc présenter plusieurs minima locaux. Nous reviendrons
sur ce point plus tard.
5.2.2. Modélisation discrète

Pour résoudre l’équation d’Euler, on commence par discrétiser la courbe. Considérons
T
donc les N points V = [v0 v1 ... vN −1 ] de la courbe discrétisée avec un pas h. On
approxime les dérivées par la méthode des différences finies (Annexes, paragraphe 1).

08/03/2012

- 37 -

Traitement et Analyse des Images Numériques

FIG. IV-23. Discrétisation d’un contour

Tout d’abord, supposons que les paramètres α et β sont invariants le long de la courbe.
L’équation Eq. IV-50 s’écrit alors :
− αv' ' (s ) + β v' ' ' ' ' (s ) = F (v )

Eq. IV-52

On discrétise :
v ' ' −2vi' ' +vi −1' '
⎛ v ' −v ' ⎞
− α ⎜ i +1 i ⎟ + β i +1
= fi
h2
⎝ h ⎠
α ⎛ v − v v − v ⎞ β ⎛ v − 2vi + 1 + xi
v − 2vi + vi −1 xi − 2vi −1 + vi − 2 ⎞
− 2 i +1
+
− ⎜ i +1 i − i i −1 ⎟ + 2 ⎜ i + 2
⎟ = fi
2
h2
h2
h⎝ h
h ⎠ h ⎝
h

α ⎛ v − v v − v ⎞ β ⎛ v − 2vi + 1 + vi
v − 2vi + vi −1 vi − 2vi −1 + vi − 2 ⎞
− ⎜ i + 1 i − i i −1 ⎟ + 2 ⎜ i + 2
− 2 i +1
+
⎟ = fi
2
h⎝ h
h ⎠ h ⎝
h
h2
h2


βvi − 2 + (− 4 β − h 2α )vi −1 + (6 β + 2 h 2α )vi + (− 4 β − h 2α )vi +1 + βvi + 2 = h4 fi

Eq. IV-53

Considérons le vecteur colonne formé des points de la courbe : V = [v0

... vN −1 ] et le

[

]

T

T

vecteur F = h 4 f0 ... h 4 f N −1
représentant les forces appliquées en chaque point.
Supposons le contour fermé ( vN −1 = v0 ). L’équation Eq. IV-53 peut se mettre sous une forme
matricielle (on pose h=1) :
AV = F

Eq. IV-54

avec :
⎡ 2α + 6 β
⎢− α − 4 β

⎢ β

0

A = ⎢ ...

0


0

⎢ β
⎢− α − 4 β


− α − 4β
2α + 6 β
− α − 4β

β
...
0
0
0

β

β

0

− α − 4β
2α + 6 β
− α − 4β
...

β

0
0
0
0

− α − 4β
2α + 6 β
...
0
0
0

...
...
...
...
...

0
0
0
0
...

... − α − 4 β
... 2α + 6 β
... − α − 4 β
...
β

β
0
0
0
...

β
− α − 4β
2α + 6 β
− α − 4β

− α − 4β ⎤

β


0

0

... ⎥

0


β

− α − 4β ⎥
2α + 6 β ⎥⎦

Notons (xi , yi ) les coordonnées cartésiennes du point vi (FIG. IV-23). L’équation Eq. IV-54
s’applique sur chacune des coordonnées x et y, ce qui signifie que l’on doit résoudre le
système :

08/03/2012

- 38 -

Chapitre IV : Segmentation – Extraction de contours

⎧ AX = Fx

⎩ AY = Fy

∂P

T
⎪⎪ Fx = − ∂x = [Fx (0 ) ... Fx ( N − 1)]

∂P
T
= Fy (0 ) ... Fy ( N − 1)
⎪ Fy = −
⎪⎩
∂y

avec

[

⎧ X = [x0

⎩Y = [ y0

]

Evq. IV-55

x1 ... xN −1 ]
T
y1 ... y N −1 ]

T

5.2.3. Approche variationnelle pour la résolution

On résout l’équation par une approche variationnelle. On introduit la variable temps t
qui modélise l’évolution de la courbe. L’équation d’Euler (Eq. IV-49, Eq. IV-50) devient :

γ

∂v
− (α (s )v' (s ))' +(β (s )v' ' (s ))' ' +∇P(v ) = 0
∂t

γ

∂v
− (α (s )v' (s ))' +(β (s )v' ' (s ))' ' = F
∂t

Eq. IV-56

Notons Vt le vecteur décrivant la courbe discrète à l’instant t. En reprenant les résultats établis
précédemment (Eq. IV-50) et en discrétisant le terme d’évolution introduit en Eq. IV-56, il
vient :

γ (V t − V t −1 ) + AV t = F t −1

( A + γI )V t = F t −1 + γV t −1

Eq. IV-57

Dans cette équation, F t −1 désigne les forces appliquées à la courbe définie à l’instant t-1,
c’est-à-dire à V t −1 . Le schéma numérique est implicite et stable. On peut montrer que la
matrice ( A + γI ) est inversible. On a donc abouti ainsi à une formulation explicite qui permet
de déduire l’état V t de la courbe à l’instant t en fonction de l’état V t −1 de la courbe à l’instant
t-1 :

(

V t = ( A + γI ) F t −1 + γV t −1
−1

5.3.

)

Eq. IV-58

Algorithme (cas d’un contour fermé, α et β constants)

Afin de fixer les idées, considérons le cas très courant où on recherche un contour
fermé, avec :
• Les paramètres α et β définissent énergie interne (Eq. IV-44) et sont supposés
constants,
• l’énergie image Eimage est définie comme indiquée ci-dessus (Eq. IV-45, Eq. IV-46,
Eq. IV-47),
• Il n’y a pas d’énergie de contexte : Eexterne=0
Dans ce cas :

08/03/2012

- 39 -

Traitement et Analyse des Images Numériques

⎧⎡ x0 (t ) ⎤
⎛ ⎡ f x ( x0 (t − 1), y0 (t − 1)) ⎤
⎡ x0 (t − 1) ⎤ ⎞
⎜⎢
⎪⎢



⎥⎟
...
...
...

⎪⎢




⎥⎟
−1 ⎜
⎪⎢ xi (t ) ⎥ = ( A + γI ) ⎢ f x ( xi (t − 1), yi (t − 1)) ⎥ + γ ⎢ xi (t − 1) ⎥ ⎟
⎜⎢
⎪⎢


⎥⎟

⎜⎢
...
...
⎪⎢ ... ⎥


⎥⎟
⎜⎢
⎪⎢ x (t )⎥

⎢⎣ xN −1 (t − 1)⎥⎦ ⎟

⎪⎣ N −1 ⎦
⎝ ⎣ f x (xN −1 (t − 1), y N −1 (t − 1))⎦

⎛ ⎡ f x ( x0 (t − 1), y0 (t − 1)) ⎤
⎡ y0 (t − 1) ⎤ ⎞
⎪⎡ y0 (t ) ⎤
⎜⎢


⎥⎟
⎪⎢ ... ⎥
...
...





⎥⎟
⎪⎢
−1 ⎜

⎪⎢ yi (t ) ⎥ = ( A + γI ) ⎜ ⎢ f x ( xi (t − 1), yi (t − 1)) ⎥ + γ ⎢ yi (t − 1) ⎥ ⎟




⎪⎢ ... ⎥
⎜⎢
...
...



⎥⎟
⎪⎢



⎢⎣ y N −1 (t − 1)⎥⎦ ⎟
⎪⎩⎢⎣ y N −1 (t )⎥⎦
⎝ ⎣ f x (xN −1 (t − 1), y N −1 (t − 1))⎦


avec

∂Eimage
(xi (t − 1), yi (t − 1))
∂x
• A est la matrice NxN définie par Eq. IV-54 qui ne dépend que des paramètres α et
β.
• γ est un paramètre de viscosité qui contrôle la vitesse de déplacement de la courbe
d’une itération à l’autre.
Il faut noter que si α β et γ sont indépendants du temps, alors la matrice ( A + γI ) n’est

f x ( xi (t − 1), yi (t − 1)) = −



inversée qu’une seule fois (du moment bien sûr que l’on ne change pas la discrétisation de la
courbe puisque la taille de la matrice A est liée au nombre N de points).
Sous ces conditions, indiquons maintenant les étapes de l’algorithme appliqué à une
image I (i , j ) de taille HxW. Dans la suite, (i , j ) désignent des coordonnées entières alors que
(x , y ) désignent des coordonnées flottantes.
3. Définir les paramètres α β et γ.
2
4. Calculer l’énergie image. Par exemple P(i , j ) = − ∇I (i , j ) .
5. Calculer champ de forces, i.e. le gradient de l’énergie de l’image : dérivée suivant x et
dérivée suivant y. Cela conduit à deux images, de taille HxW, que l’on notera f x (i , j )
et f y (i , j ) , définissant un champ de vecteurs. Dans l’exemple précédent,

[ f (i , j ), f (i , j ))] = −∇(− ∇I (i , j )
x

y

2

) = ∇(∇I (i , j ) ).
2

6. Initialiser la courbe : par des méthodes de traitement d’image ou via une interface
utilisateur. Cette courbe doit être définie par N points régulièrement espacés.
Soulignons que les coordonnées (xi , yi ) des points de la courbe ne sont pas

nécessairement entières. Dans la suite, on notera x = [xi ]i =0 ,..., N −1 le vecteur colonne des
T

abscisses et y = [ yi ]i =0 ,..., N −1 le vecteur colonne des ordonnées.
T

7. Calculer A (taille NxN) et l’inverse de la matrice ( A + γI ) , notée invAI.
8. Définir le nombre d’itérations nb_iter effectuées avant de redistribuer les points qui

définissent la courbe discrète.

08/03/2012

- 40 -

Chapitre IV : Segmentation – Extraction de contours

9. Tant que « le contour évolue ! » (i.e. on n’a pas atteint un minimum d’énergie)
o Pour iter = 0 → nb_iter, appliquer la déformation du snake :
ƒ Calculer les forces f x ( xi , yi ) et f y ( xi , yi ) en chaque point de la courbe

(xi , yi ) .

Ces valeurs sont déduites des images f x (i , j ) et f y (i , j ) par

interpolation, puisque les coordonnées (xi , yi ) ne sont pas entières.
Notons
ƒ

f x = [ f x ( xi , yi )]i=0 ,...,N −1
T

[

]

f y = f y ( xi , yi )

et

T
i =0 ,...,N −1

les vecteurs

stockant ces valeurs.
Calculer la nouvelle position de chaque point de la courbe :
⎧x = invAI * (γx + f x )

⎩y = invAI * (γy + f y )

o Interpoler pour redéfinir une courbe discrète avec des points régulièrement
espacés.
o Mettre à jour A et l’inverse de la matrice ( A + γI ) , notée invAI.

Un programme Matlab de cet algorithme peut-être trouvé sous [15]

5.4.

Exemples

Prenons une forme géométrique simple, en blanc sur fond noir. Le contour est
initialisé à un cercle. La figure montre le champ de forces qui converge vers le contour, et qui
va itérativement déformer le contour initial pour le faire coller au contour recherché.

Image source I

Gradient ∇I (x , y )

(f

2

x

, f y ) = ∇ ∇I ( x , y )

2

(à droite, zoom sur le champ de force autour du contour supérieur)

Déformation du contour (α = β =1, 175 itérations )
FIG. IV-24. Exemple de segmentation par contour actifs

5.5.

Cas des contours non fermés

5.5.1. Contours à extrémités libres

On se place dans le cas d’un contour ouvert à extrémités libres, avec les conditions aux
limites :

08/03/2012

- 41 -

Traitement et Analyse des Images Numériques

⎧⎪v' ' (0 ) = v( 3 ) (0 )

⎪⎩v' ' (1) = v( 3 ) (1)

Eq. IV-59

La formalisation du problème est similaire. La résolution montre que la matrice A s’écrit :
⎡ β
⎢− α − 2 β

⎢ β

0
A=⎢
⎢ ...

0


0

0
⎣⎢

− 2β

β

0

0

0

... 0

0

0

0

2α + 5 β
− α − 4β

− α − 4β
2α + 6 β

β
− α − 4β

0
0

0
0

... 0
... 0

0
0

0
0

0
0

β






0
−⎥

...

β

− α − 2β ⎥

β
⎦⎥
0

0
0

β

− α − 4β

2α + 6 β

− α − 4β

... 0

0

0

0

...

...

...

...

... ... ...

...

...

...

0
0

0
0

0
0

0
0

0
0

... β
... 0

− α − 4β

β

2α + 5 β
− α − 4β

− α − 2β
2α + 5 β

0

0

0

0

0

... 0

0

β

− 2β

Eq.
IV-60

5.5.2. Contours à extrémités fixes

La matrice A prend une nouvelle forme :
⎡ 2α + 6 β
⎢− α − 4 β


β

0
A=⎢

...

0


0

0
⎢⎣

5.6.

− α − 4β

β

0

2α + 6 β
− α − 4β

− α − 4β
2α + 6 β

β

0

− α − 4β

β

0

β

− α − 4β

2α + 6 β

− α − 4β

...
0

...
0

...
0

...
0

0
0

0
0

0
0

0
0

0 ⎤
0 ... 0 0 0 0 0 ⎥⎥
0 ... 0 0 0 0 0 ⎥

β ... 0 0 0 0 0 ⎥

... ... ... ... ... ... ... ⎥

0 ...


0 ... 0

0 ... 0 0
⎥⎦
0

... 0

0

0

0

Eq. IV-61

Gradient Vector Flow

La méthode proposée par Kass et al [14] a constitué les fondements de la segmentation
par contour actif. Avec cette approche, on peut obtenir directement des contours continus et
même fermés. Néanmoins tous les problèmes ne sont pas résolus. Les principaux sont :
• La nécessité de définir un contour initial proche du contour recherché, parce que la
fonctionnelle d’énergie n’est pas convexe et que les forces images dérivées du
gradient agissent localement (FIG. IV-25).
• La difficulté du contour à converger dans les concavités des contours.

Image source I(x,y)

Carte de contour f(x,y)

Champ de vecteurs F = ∇f

FIG. IV-25. Champ de force pour un snake « traditionnel »

08/03/2012

- 42 -

Chapitre IV : Segmentation – Extraction de contours

Ces deux problèmes ont été adressés par Xu et al. [15][16]. Les auteurs proposent une
nouvelle définition des forces images, que nous allons introduire ci-dessous.
Dans la formulation classique, on résout l’équation (Eq. IV-50) :
− (α (s )v' (s ))' +(β (s )v' ' (s ))' ' = F (v ) avec F (v ) = −∇P(v )

avec F (v ) une force définie comme le gradient d’une « carte de contour », typiquement le
2
2
module au carré du gradient de l’image f ( x , y ) = ∇I (x , y ) . Ainsi F = ∇f = ∇ ∇I ( x , y ) .
Les vecteurs du champ de forces obtenu pointent vers les contours. En revanche, ils sont nuls
sur les zones assez uniformes, ce qui explique la nécessité d’initialiser le contour actif aussi
près que possible du contour que l’on recherche. La figure FIG. IV-25 illustre cela.

(

)

Xu et al. proposent de calculer le champ de forces v( x , y ) = [u ( x , y ), w( x , y )] qui
minimise une fonction d’énergie définie par :

ε = ∫∫ μ (u x2 + u y2 + wx2 + wy2 ) + ∇f v − ∇f dxdy
2

2

Eq. IV-62

Dans cette équation, ux (resp. uy) désigne la dérivée partielle de la fonction u par rapport à x
(resp. y).
Le champ de vecteurs solution est appelé Gradient Vector Flow (GVF). En quoi cette
nouvelle formulation va résoudre les problèmes mentionnés en introduction, nous le verrons
en pratique à la fin de ce paragraphe… Dans l’immédiat, nous nous concentrons sur la
résolution du problème, afin de trouver un algorithme qui permette de calculer
numériquement le GVF.
On peut montrer que minimiser l’énergie Eq. IV-62 revient à résoudre les équations
d’Euler :

(
μ∇ w − (w − f )( f

)
)= 0

μ∇ 2u − (u − f x ) f x 2 + f y 2 = 0
2

2

y

x

+ fy

Eq. IV-63

2

∇ 2 désignant le Laplacien. Ces équations sont résolues en introduisant la variable t et en
traitant les fonctions u et w comme des fonctions du temps :

(
w ( x , y ,t ) = ∇ w( x , y ,t ) − (w( x , y ,t ) − f ( x , y ))( f ( x , y )

)
+ f (x , y ) )

ut ( x , y ,t ) = μ∇ 2u ( x , y ,t ) − (u ( x , y ,t ) − f x ( x , y )) f x ( x , y ) + f y ( x , y )
2

2

t

2

y

x

2

2

Eq. IV-64

y

Ces nouvelles équations sont appelées équations de diffusion. On peut les simplifier en
utilisant les notations suivantes :
ut ( x , y ,t ) = μ∇ 2u ( x , y ,t ) − b( x , y )u ( x , y ,t ) + c1 (x , y )

wt ( x , y ,t ) = μ∇ 2 w(x , y ,t ) − b(x , y )w(x , y ,t ) + c 2 ( x , y )

08/03/2012

- 43 -

Eq. IV-65

Traitement et Analyse des Images Numériques

⎧b( x , y ) = f x ( x , y )2 + f y ( x , y )2
⎪⎪
avec ⎨c1 ( x , y ) = b( x , y ) f x ( x , y )
⎪ 2
⎪⎩c ( x , y ) = b( x , y ) f y ( x , y )

On discrétise Eq. IV-65 par la méthode des différences finies. Notons i et j les
coordonnées discrètes correspondant à x et y, n les itérations modélisant l’évolution
temporelle, Δx, Δy, Δt les pas de discrétisation. Les dérivées sont approximées par :

1 n +1
(
ui , j − uin, j )
Δt
1
wt = (win, +j 1 − win, j )
Δt
1
(u + u + u + u − 4ui , j )
∇ 2u =
ΔxΔy i + 1, j i , j + 1 i −1, j i , j −1
1
∇2w =
(w + wi , j +1 + wi −1, j + wi , j −1 − 4 wi , j )
ΔxΔy i + 1, j
ut =

Eq. IV-66

En remplaçant ces approximations dans Eq. IV-65, on obtient itérativement la solution
de Eq. IV-64 par :

uin, +j 1 = (1 − bi , j Δt )uin, j + r (uin+ 1, j + uin, j + 1 + uin−1, j + uin, j −1 − 4uin, j ) + ci1, j Δt

win, +j 1 = (1 − bi , j Δt )win, j + r (win+ 1, j + win, j + 1 + win−1, j + win, j −1 − 4 win, j ) + ci2, j Δt

avec :

Δt
r=μ
ΔxΔy

On peut montrer que le schéma numérique converge pour Δt ≤

Eq. IV-67

ΔxΔy
.


L’équation Eq. IV-67 permet une implémentation algorithmique directe pour le calcul
du GVF. Les valeurs initiales ui0, j et wi0, j sont données par le gradient de la carte de contours
∇f . L’algorithme prend en entrée ∇f , μ et le nombre d’itérations pour la diffusion.

Voyons maintenant concrètement l’effet de cet algorithme de diffusion appliqué au
gradient de la carte de contour ∇f (même exemple qu’en figure FIG. IV-25), en fonction du
nombre d’itérations. On a pris μ=0.2.

0 itération
5 itérations
10 itérations
50 itérations
FIG. IV-26. GVF obtenu en fonction du nombre d’itérations, pour l’image de la figure FIG. IV-25

08/03/2012

- 44 -

Chapitre IV : Segmentation – Extraction de contours

On observe bien le processus de diffusion qui est d’autant plus important que le nombre
d’itérations est élevé, augmentant la zone de « capture » de l’algorithme de contour actif.
Ainsi il n’est plus nécessaire d’avoir une initialisation fine du contour. Il est aussi intéressant
de remarquer comment les vecteurs pointent vers les contours à l’intérieur de la concavité
(FIG. IV-27). La diffusion permet d’obtenir des champs de forces qui sont une sorte
d’interpolation des différentes forces dues aux différents contours en présence. Dans la figure
FIG. IV-27, la direction des vecteurs montre bien que le contour peut converger à l’intérieur
de la concavité, ce qui n’était pas le cas avant la diffusion.

FIG. IV-27. Zoom sur le GVF obtenu avec 50 itérations (FIG. IV-26), dans la concavité

Les figures FIG. IV-28 et FIG. IV-29 comparent la segmentation avec et sans
diffusion. L’algorithme qui permet de faire évoluer le contour actif est le même dans les deux
cas : il s’agit de celui décrit au paragraphe 5.3. On constate que le champ de forces GVF est
beaucoup plus cohérent et homogène que le champ de forces calculé par le simple gradient,
qui lui ne peut agir que très localement (FIG. IV-28). Le contour converge correctement avec
le GVF, alors qu’il ne converge pas avec le champ de forces non diffusé (FIG. IV-29).
Les GVF snakes apportent donc une solution séduisante au problème d’extraction de
contours. Il faut néanmoins garder en mémoire que cette méthode ne permet d’extraire qu’un
seul contour à la fois. Elle est assez coûteuse en calculs. Enfin, les paramètres sont nombreux
et le réglage optimal de ces paramètres est complexe.

6.

Transformée de Hough

La transformée de Hough permet de détecter la présence de formes géométriques
modélisables par une équation. Par exemple,
- droite :
- cercle :

y = ax + b

(x − a ) + ( y − b )
2

2

ou ρ = x cos(θ ) + y sin(θ ) ,
= c2 .

Prenons le cas de la recherche de droites (FIG. IV-30). Considérons un point Pi de
l’image, de coordonnées (xi,yi). Une infinité de droites, définies par des couples (a,b), passent
par ce point. Considérons maintenant le plan (a,b) définissant l’espace des paramètres a et b.
Toutes les droites passant par Pi satisfont à l’équation b = − xi a + yi , et sont donc également
représentables par une droite dans l’espace des paramètres. Si on considère un deuxième point

08/03/2012

- 45 -

Traitement et Analyse des Images Numériques

Image source

Image de gradient

(

Champ des forces ∇ ∇I (i , j )

2

) (cas 1)

Initialisation du contour actif

Gradient Vector Flow, 75 itérations (cas 2)

FIG. IV-28. Segmentation par contour actif, comparaison des résultats obtenus sans et avec diffusion

Evolution du contour actif « classique » (cas 1) : itérations 10,20,30,40,50

Evolution du GVF snake (cas 2) : itérations 10,20,30,40,50
FIG. IV-29. Segmentation par contour actif, comparaison des résultats obtenus sans et avec diffusion

08/03/2012

- 46 -

Chapitre IV : Segmentation – Extraction de contours

Pj(xj,yj), alors de nouveau on a une infinité de couple (a’,b’) modélisant les droites qui passent
par ce point, représentables dans l’espace des paramètre par la droite b = − x j a + y j . En
revanche, un seul couple (aD,bD) correspond à la droite passant par Pi(xi,yi) et Pj(xj,yj). Ce
couple est à l’intersection des droites d’équations b = − xi a + yi et b = − x j a + y j . On peut
poursuivre ce raisonnement en prenant un troisième point sur la droite (Pi Pj ) . Dans l’espace
des paramètres, ce point définit une nouvelle droite qui passe elle aussi par (aD,bD).
Ceci introduit naturellement la notion de vote. Soit V(a,b) un accumulateur, initialisé à
0, chaque élément de cet accumulateur représentant le vote pour une unique droite. L’image
source est par exemple une carte de contours, où les points de contour sont à 1, les autres à 0.
Pour chaque pixel de contour, on considère tous les couples (a,b) modélisant les droites
passant par ce point et on incrémente l’accumulateur pour chacun de ces couples. Les maxima
de l’accumulateur correspondent aux paramètres des droites présentes dans l’image. En
pratique, il faut bien sûr discrétiser l’espace des paramètres pour définir l’accumulateur.

FIG. IV-30 : principe de la transformée de Hough

L’algorithme de détection des droites, appliqué à une carte de contour, est donc le
suivant :
- Créer un accumulateur V(i,j)
- Pour chaque point non nul ( xk , yk ) du plan image
Pour chaque valeur de aq permise (valeur quantifiée)
Calculer b = − xk aq + yk
Arrondir b à la valeur quantifiée bq la plus proche
Transformer les valeurs quantifiées en indices (i,j) dans l’accumulateur
Accumuler : V(i,j)= V(i,j)+1
- Rechercher les maxima de l’accumulateur, en déduire les équations de droites
correspondantes.
En pratique, la représentation des droites par l’équation y = ax + b n’est pas très bien
adaptée à la quantification, car le paramètre a peut prendre toutes les valeurs. On préfère donc
la forme polaire (FIG. IV-31) ρ = x cos(θ ) + y sin(θ ) . Les deux paramètres sont dans ce cas
bornés :

[

]

⎧ ρ ∈ 0 , 2 N avec N la dimension maximale de l' image


⎡ π ⎤
⎪θ ∈ ⎢− 2 ,π ⎥




08/03/2012

- 47 -

Eq. IV-68

Traitement et Analyse des Images Numériques

FIG. IV-31. Codage d’une droite en polaire

La détection des maxima dans l’accumulateur s’effectue généralement après lissage,
afin de rassembler les votes similaires. On applique ensuite un seuil afin de ne retenir que les
maxima significatifs.
La figure FIG. IV-32 montre un exemple de transformée de Hough appliquée à une
image simple de contours. Le principe est généralisable à toute forme géométrique qui peut
être décrite par une équation paramétrique. La dimension de l’accumulateur est alors égale au
nombre de paramètres à définir.

FIG. IV-32. Détection de droites par transformée de Hough. (a) image de contour, (c) accumulateur de
Hough (représentation normalisée) : on remarque la présence de 4 maxima (en blanc) correspondant aux 4
droites dans l’image.

08/03/2012

- 48 -

Chapitre V
Segmentation de régions

Nous abordons maintenant la méthode duale de l’extraction de contours : la
segmentation de régions. Cette fois, il s’agit de réaliser une partition de l’image en régions
distinctes en se basant sur des propriétés de similarité des pixels appartenant à une même
région. Nous traitons dans ce cours six techniques importantes de segmentation : le seuillage,
les k-moyennes, la croissance de région, l’approche division-rassemblement, l’approche
multi-résolution, les méthodes statistiques de classification des pixels en régions. Les deux
dernières sont plutôt appliquées sur les images texturées.
Le seuillage consiste à comparer le niveau de gris de chaque pixel de l’image à des
seuils, en vue d’attribuer le point considéré à une certaine classe (région). Cette technique
traite des cas simples, telle que la classification de pixels en « fond » ou « objet ».
L’algorithme des k-moyennes est un algorithme itératif qui permet de classer les pixels
en k-classes, en comparant leur valeur à la moyenne des régions. Les moyennes des régions et
le classement des pixels évoluent itérativement, d’une initialisation quelconque jusqu’à
stabilisation. Cet algorithme est très populaire, à cause de sa simplicité et des très bons
résultats qu’il permet d’obtenir sur des images pas trop complexes.
La technique de croissance de régions regroupe récursivement des points adjacents qui
possèdent un ou plusieurs attributs homogènes, par exemple un niveau de gris similaire. Cette
méthode permet d’obtenir un ensemble de pixels connexes à partir d’une unique initialisation.
L’approche division/rassemblement s’apparente à la précédente, car elle regroupe
également des pixels voisins similaires, mais avec une méthode différente.
La segmentation multi-résolution s’appuie sur une analyse des images à des niveaux
de résolution différents. Les niveaux de résolution dégradés sont, par exemple, obtenus par
des moyennages successifs de pixels. La classification se fait au niveau de résolution le plus
bas, et est affinée aux niveaux de résolution supérieurs. Cette technique donne de bons
résultats sur les images texturées.
Les méthodes de classification statistiques sont basées sur des fonctions de densité de
- 49 -



Documents similaires


tnimages 2012 vol2 v1
rapport tp1 m2 analyse image sauvageon vellutini
a robust region based multiscale image fusion scheme for mis
foreground background segmentation using temporal and spatial
lidar topographie
tutobrick


Sur le même sujet..