Université de Toulouse III, FSI – Département de Mathématiques

Licence 2, module “Calcul scientifique”, 2019–2020

Projet: “Traitement d’images, détection de

contours”
Février 2020

Détecter les contours d’une image constitue une étape préliminaire à de nombreuses applica-
tions pratiques en traitement d’images : analyse d’images médicales, reconnaissance de formes,
restauration... Le but de ce projet est d’aborder quelques méthodes simples de détection de
contours.

1 Fichiers fournis et installation

Avant tout, il faut installer la bibliothèque image de Python grâce à cette commande :
pip install image -- user

L’ensemble des fichiers de travail peut être téléchargé à cette adresse :

[Link]

Vous trouverez dans cette archive la bibliothèque imtools (voir par exemple le listing 1) qui
simplifie les tâches courantes de gestion des images, ainsi que quelques images pour vos expé-
riences numériques. Pour les utiliser, une fois image installée comme expliqué plus haut, il suffit
de recopier tous les fichiers dans votre répertoire de travail.

2 Images digitales

Figure 1 – Échantillonnage d’une image sur une grille de taille 32 × 32.

1
 Une manière de définir une image digitale u de taille M × N est de la voir comme une
matrice  
um,n 0≤m<M ∈ MM,N (R). (1)
0≤n<N

On notera qu’à la différence des conventions habituelles en mathématiques, on a numéroté

les lignes et les colonnes à partir de zéro, afin de respecter l’indexation des tableaux numpy de
Python. Pour simplifier, on ne considérera que des images en niveaux de gris. Chaque case de la
matrice est un pixel (abréviation de picture element). Habituellement une image digitale prend
des valeurs comprises entre zéro (pixel noir) et 255 (pixel blanc), même si d’autres conventions
sont possibles.
Une image digitale correspond à l’échantillonnage d’une image analogique f , vue comme une
application f : R2 → R représentant le flux lumineux théorique reçu au niveau d’un plan, en
fonction de la position sur ce plan. Un récepteur placé derrière ce plan ne perçoit qu’une portion
bornée du flux, sur le pavé [a, b] × [c, d]. Sur ce pavé on construit une grille de points (xm , ym )
uniformément répartis (pour simplifier, on va supposer que la grille d’échantillonnage utilisée
a un pas égal à 1), et on considère que l’image échantillonnée est obtenue tout simplement en
posant (voir figure 1) :
um,n = f (xm , yn ). (2)
Dans les calculs qu’on sera amené à traiter, on fera parfois référence à des pixels d’indices
situés à l’extérieur du pavé {0, . . . , M − 1} × {0, . . . , N − 1}. Dans ce cas, plusieurs conventions
sont possibles (voir figure 2) :

• considérer que l’image est constante (par exemple égale à zéro) à l’extérieur de la zone
d’échantillonnage ;

• prolonger l’image en attribuant la même valeur que le pixel du bord situé à distance
minimale du point considéré ;

• prolonger l’image par parité et périodicité ;

• prolonger l’image par périodicité, en posant um,n = um0 ,n0 où m0 est le reste de la division
euclidienne de m par M et n0 celui de la division euclidienne de n par N .

On préférera l’utilisation de cette dernière convention, car elle est facile à implémenter en
Python : il suffit par exemple d’écrire u[m % M, u % N] pour effectuer cette périodisation.

Figure 2 – Prolongements possibles d’une image hors de son support (en rouge).

On fournit une petite bibliothèque [Link] qui introduit quelques fonctions permettant
de manipuler des images en Python :

• open_image(filename) permet d’ouvrir un fichier image en niveaux de gris (habituel-

lement au format PNG) et d’obtenir les données sous forme d’un tableau numpy à deux
dimensions ;

2
 • save_image(u, filename) permet d’enregistrer un tableau numpy dans un fichier PNG.
Les données sont converties en nombres entiers compris entre 0 et 255 : il peut donc être
nécessaire de les ajuster avant d’enregistrer (si par exemple on manipule des tableaux
avec des valeurs négatives) ;

• display_image(u) affiche une image à l’écran. On notera que par défaut, Python va
automatiquement interpréter les données de manière à ce que la valeur la plus basse du
tableau corresponde à du noir et la valeur la plus haute à du blanc. Ce n’est pas toujours
souhaitable. On pourra, si nécessaire, utiliser par exemple

display_image(u, vmin = 0, vmax = 255)

pour ajuster l’échelle (ici : le noir à zéro et le blanc à 255). Il peut être également nécessaire
d’ajuster la taille d’affichage, car par défaut les notebook affichent les figures avec une
petite taille, ce qui empêche de percevoir les détails. On pourra utiliser par exemple :

[Link](figsize = (10, 10))

avant l’utilisation de display_image pour augmenter la taille d’affichage.

• display_images(l) affiche une liste d’images en une seule fois.

3 Produit de convolution
3.1 Filtrage des images
Si u et h sont deux images de même taille M × N , on définit le produit de convolution u ∗ h
par X
(u ∗ h)m,n = um−k,n−l hk,l . (3)
0≤k<M
0≤l<N

On peut vérifier que le produit de convolution ainsi défini est bilinéaire, commutatif, associatif,
et que l’élément neutre est la masse de Dirac δ0,0 définie par
 
1 0 ... 0
0 0 ... 0


δ0,0 =
 .. .. .
..  (4)
. . .
0 0 ... 0

De manière similaire, on notera δm,n la masse de Dirac en (m, n), qui vaut zéro partout sauf à
l’indice (m, n) où elle vaut 1.
Mathématiquement, u et h jouent des rôles symétriques dans la définition (3). Toutefois
dans le contexte du traitement d’image, on considère souvent que l’un des opérandes (noté en
général u) est une image tandis que l’autre (noté en général h) est un filtre. Avec ces conventions,
le calcul de u ∗ h est appelé un filtrage de l’image.

3
3.2 Calcul en pratique lorsque le support du filtre est petit
Lorsque h a peu de coefficients non nuls, la façon la plus efficace de calculer le produit de
convolution u ∗ h en Python repose sur l’observation suivante. Si on note τk,l l’opérateur linéaire
de décalage (avec la convention de périodisation des indices comme expliqué précédemment)

M
M,N (R) → MM,N (R)
τk,l :  (5)
(um,n ) 7→ (um−k,n−l ),
alors on peut écrire : X
u∗h= hk,l τk,l (u). (6)
(k,l)∈Spt h

Ici, Spt h désigne le support de h, c’est à dire l’ensemble des indices (k, l) tels que hk,l 6= 0.
La bibliothèque imtools fournit une fonction pour calculer τk,l (u) : il suffit d’appeler
shift_image(u, k, l)
pour obtenir le tableau correspondant. Cette fonction est beaucoup plus rapide que le calcul
direct avec des boucles for, notoirement lentes en Python.
Pour simplifier les notations, lorsqu’un filtre a un support petit, on se contentera de le noter
sous la forme d’un tableau réduit à la taille de son support. Pour lever d’éventuelles ambiguïtés,
on encadrera le coefficient situé à l’indice (0, 0). Ainsi, le filtre δ0,0 peut se noter sous forme
abrégée :  
1 0 ... 0  
0 0 0
0 0 . . . 0
   
δ0,0 = 
 .. .. ..  = 0 1 0 = 1 . (7)
  
 
. . . 0 0 0
0 0 ... 0
Le listing 1 donne un exemple d’implémentation efficace en Python du filtre

−2

0 0
h = −3δ1,−1 + 5δ0,0 − 2δ−1,1 = 0 5 0 . (8)
 

−3 0 0

from i m t o o l s import ∗

def f i l t r a g e ( u ) :
return −3 ∗ s h i f t _ i m a g e ( u , 1 , −1) + 5 ∗ u − 2 ∗ s h i f t _ i m a g e ( u , −1, 1 )

u = open_image ( ’ l e n a . png ’ )
v = f i l t r a g e (u)
p l t . f i g u r e ( f i g s i z e = (10 , 5) )
display_images ( [ u , v ] )

Listing 1 – Implémentation efficace du filtre −3δ1,−1 + 5δ0,0 − 2δ−1,1

4 Discrétisation de quelques opérateurs différentiels

La plupart des méthodes de détection de contours sont basées sur l’étude des dérivées de
l’image. Mais il reste à définir ce qu’on entend par là, puisque les images digitales sont des
objets de nature discrète. Prenons l’exemple de la première dérivée selon la première variable.
Si on regarde l’équation (2), en supposant que f est suffisamment régulière, il s’agit donc de
∂f
∂x
. On se pose alors la question d’estimer l’échantillonnage de la dérivée ∂f
∂x
(xm , yn ) connaissant
seulement l’échantillonnage u de f .

4
4.1 Schémas de différences finies d’ordre 1
Plusieurs schémas numériques ont été proposés pour résoudre ce problème. L’une des possibilités
consiste à supposer que sur un voisinage U de (xk , yk ) suffisamment grand, l’image analogique f
est affine :
∀(x, y) ∈ U : f (x, y) = f (xm , yn ) + α(x − xm ) + β(y − yn ). (9)
Dans ce cas, en supposant que U contient les pixels voisins, on a :
∂f
(xm , yn ) = α = f (xm + 1, yn ) − f (xm , yn ) = um+1,n − um,n (10)
∂x
∂f
(xm , yn ) = β = f (xm , yn + 1) − f (xm , yn ) = um,n+1 − um,n , (11)
∂y
et dans ces conditions le calcul des dérivées de u peut s’obtenir en effectuant le produit de
convolution de u respectivement par les deux filtres :
!
  1
δ−1,0 − δ0,0 = 1 −1 δ0,−1 − δ0,0 = . (12)
−1
On parle alors de dérivée par différences finies décentrées à droite. Mais d’autres choix sont
possibles, comme par exemple :
!
  1
δ0,0 − δ1,0 = 1 −1 δ0,0 − δ0,1 = . (13)
−1
On parle cette fois de dérivées décentrées à gauche. On peut également obtenir des dérivées
centrées en considérant la moyenne des deux versions :
 
1
∂ 1 1   ∂ 1 1 
: (δ−1,0 − δ1,0 ) = 1 0 −1 : (δ0,−1 − δ0,1 ) =  0  . (14)
∂x 2 2 ∂y 2 2
−1

4.2 Schémas de différences finies d’ordre 2

Il est possible d’itérer les schémas d’ordre 1 pour calculer les dérivées d’ordre supérieur. Étant
donnée l’associativité du produit de convolution, le filtre correspondant à une double différen-
tiation est obtenu en convoluant les deux filtres de différentiation. Voici ce qu’on obtient à
l’ordre 2 :
∂2  
: (δ −1,0 − δ 0,0 ) ∗ (δ0,0 − δ 1,0 ) = 1 −2 1 (15)
∂x2  
1
∂2
: (δ 0,−1 − δ 0,0 ) ∗ (δ0,0 − δ 0,1 ) = 
−2

(16)
∂y 2
 
1

−1

2 1 0
∂ 1 1
: (δ0,−1 − δ0,1 ) ∗ (δ−1,0 − δ1,0 ) =  0 0 0 . (17)

∂x∂y 4 4
−1 0 1
Un cas particulier est celui de l’opérateur laplacien. En faisant comme précédemment une
moyenne entre différents décalages possibles, on peut donner cette version qui a un comporte-
∂2 ∂2
ment un peu plus anisotrope que la simple addition des filtres de ∂x 2 et ∂y 2 :

 
2 1 2 1 1
∂ ∂ 1
∆= + : 1 −8 1 . (18)

∂x2 ∂y 2 8
1 1 1
On en trouvera un exemple d’implémentation dans le listing 2

5
def l a p l a c i a n ( u ) :
v = 0.0 ∗ u
f o r x in range ( −1 , 2 ) :
f o r y in range ( −1 , 2 ) :
i f x == 0 and y == 0 :
r = −1
else :
r = 1 / 8
v += r ∗ s h i f t _ i m a g e ( u , x , y )
return v

Listing 2 – Calcul du laplacien

4.3 Sensibilité au bruit

En pratique, l’échantillonnage décrit dans l’égalité (2) fait souvent intervenir un bruit, qui peut
être dû à différents facteurs (imprécision des capteurs, mouvement lors de la mesure, etc...). Pour
simplifier, on supposera que le bruit est additif et décrit par une famille de variables aléatoires
bm,n indépendantes, gaussiennes, centrées de lois N (0, σ 2 ). La formule d’échantillonnage devient
alors
um,n = f (xm , yn ) + bm,n , (19)
et en pratique on considère que l’image obtenue correspond à une réalisation des variables bm,n .
On pourra simuler l’ajout d’un bruit à une image avec la fonction donnée dans le listing 3.

Figure 3 – Effet de l’ajout d’un bruit sur le calcul des dérivées de l’image

def n o i s e ( u , sigma ) :
return u + np . random . normal ( 0 , sigma , s i z e = u . shape )

u = open_image ( ’ p e p p e r s . png ’ )
v = n o i s e (u , 30)
p l t . f i g u r e ( f i g s i z e = (15 , 10) )
display_images ( [ u , v ]

Listing 3 – Exemple de simulation du bruitage d’une image

6
 Comme on peut le constater sur la figure 3, les opérateurs de dérivation sont très fortement
sensibles au bruit, ce qui va nuire à la détection des contours. Il est donc préférable d’effectuer
un traitement préalable de l’image pour réduire le bruit.

4.4 Réduction du bruit

Une méthode classique de réduction du bruit en traitement d’images est l’utilisation de l’EDP
(équation aux dérivées partielles) dite de la chaleur :

 ∂v

− ∆v = 0
∂t (20)
v(0, x) = f (x).


Ici v : R+ × R2 → R est une fonction des deux variables t ≥ 0 et x ∈ R2 . En partant de la

condition initiale v(0, ·) = f on fait évoluer v au cours du temps de la même manière que la
chaleur se répartit au sein d’un solide homogène. Il se produit alors une régularisation de la
fonction initiale f qui estompe peu à peu le bruit.
On peut montrer par exemple que sous certaines conditions, notamment de régularité de f ,
alors v(t, ·) converge uniformément vers la fonction constante égale à la valeur moyenne de f
lorsque t → ∞. Il n’est pas question d’aller jusque là dans les cas qui nous intéressent, on se
contentera d’aller jusqu’à t = T , où T est une constante à ajuster selon le degré d’atténuation
du bruit désiré (voir figure 4). Il faut trouver le bon compromis entre le floutage des contours
et la réduction du bruit.
Pour résoudre l’équation (20) avec des images digitales, on pourra utiliser le schéma numé-
rique discret suivant : 
u(0) = u
(21)
u(tk+1 ) = u(tk ) + δt∆u(tk ),

où ∆ désigne la convolution par le laplacien discrétisé donné en (18).

Figure 4 – Réduction du bruit par l’équation de la chaleur pour différentes valeurs de T

7
 f (x) f 0 (x) f 00 (x)

x
x0

x x
x0 x0

Figure 5 – Exemple d’un contour en dimension 1 localisé en x = x0 avec le comportement des

deux premières dérivées

5 Méthodes de détection des contours

Avant d’aller plus loin, il faut définir un peu plus précisément ce qu’on entend par la notion
de « contour ». Intuitivement, l’œil percevra un contour là où la luminosité de l’image varie
brusquement. En dimension 1, on peut considérer qu’un contour correspond à un extremum
local de la dérivée première, ou encore à une annulation de la dérivée seconde (voir figure 5).

5.1 Seuillage du gradient

Pa analogie avec la dimension 1, on peut considérer que les contours d’une image sont localisés
là où la norme du gradient présente un maximum local.

def t h r e s h o l d ( u , s ) :
return ( u > s ) . a s t y p e ( f l o a t )

Listing 4 – Implémentation du seuillage : la fonction renvoie un tableau contenant des zéros et

des uns, qui correspondent à la fonction caractéristique de l’ensemble considéré.

Figure 6 – Seuillage de la norme du gradient de Lena pour différentes valeurs de s

8
 Une première idée pour trouver leur localisation consiste à seuiller la norme du gradient,
c’est à dire à considérer que les pixels qui font partie des contours de l’image sont ceux pour
lesquels la norme du gradient est plus grande que s, où s > 0 est un seuil fixé arbitrairement.
Formellement, on peut écrire :

(m, n) est sur un contour de l’image ⇐⇒ k∇um,n k > s. (22)

Une fonction de seuillage peut se définir très simplement en Python (voir listing 4). On peut
voir quelques exemples de seuillage du gradient sur la figure 6.
On rappelle qu’avant d’utiliser la méthode, il peut être nécessaire de débruiter l’image
comme expliqué dans la section 4.4, car le calcul approché du gradient par différences finies est
très sensible au bruit.

5.2 Annulations du laplacien

La méthode précédente est facile à implémenter, mais présente plusieurs défauts :

• il est difficile de savoir quel seuil s choisir a priori ;

• certains contours rectilignes sont mal localisés : on trouve des lignes épaisses là où on
voudrait une ligne ponctuelle.

On peut améliorer les résultats en cherchant les annulations du laplacien de l’image.

Pour simplifier, on considère que le laplacien s’annule entre deux pixels voisins si les valeurs
qu’il prend sur chaque pixel sont de signes opposés. Toutefois, cette méthode risque de détecter
les petites oscillations locales dues au bruit, même après débruitage (qui atténue le bruit sans le
faire disparaître totalement). Pour éviter cela, il faut appliquer un double seuillage du laplacien.
On fixe un seuil  > 0 et on définit :

+1
si ∆um,n > 


vm,n = 0

si − ≤ ∆um,n ≤  (23)


−1 si ∆um,n < −.

Une fois que v a été calculée, il ne reste plus qu’à détecter les changements de signe comme sur
la figure 7. On fournit (voir listing 5) l’implémentation d’une fonction en Python qui utilise ce
procédé pour détecter les annulations d’une image.

− − − − +

− − + − +

− + + + +

− − − + +

− − − − +

Figure 7 – Détection des changements de signe d’une image (en rouge).

9
def d e t e c t _ v a n i s h ( u , e p s i l o n ) :
v1 = u < −e p s i l o n
v2 = u > e p s i l o n
w = np . z e r o s ( u . shape , bool )
f o r ( x , y ) in ( ( 0 , 1 ) , ( 1 , 0 ) , ( 1 , 1 ) ) :
w += v1 ∗ s h i f t _ i m a g e ( v2 , x , y ) + v2 ∗ s h i f t _ i m a g e ( v1 , x , y )
return w . a s t y p e ( f l o a t )

Listing 5 – Dt́ection des annulations d’une image.

5.3 Annulation de la dérivée seconde dans la direction du gradient

Il est possible d’améliorer encore la méthode précédente. On rappelle que la matrice hessienne
d’une fonction réelle f de deux variables et de classe C 2 est définie par :

∂ 2f ∂ 2f
 

2
(a, b) (a, b)
Hessf (a, b) =  ∂x ∂x∂y

. (24)
 
 ∂ 2f ∂ 2f 
(a, b) 2
(a, b)
∂x∂y ∂y

Il s’agit de la matrice d’une forme bilinéaire (qui est la différentielle seconde de f au point (a, b)
considéré). On va alors chercher les changements de signe non plus du Laplacien de l’image,
mais de la dérivée seconde dans le sens du gradient, qui est donnée (dans le cas d’une fonction)
par :
>
∇f Hessf ∇f. (25)

10