En Python les matrices sont réprésentées par des listes de listes:
M = [ [1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12] ]
print(M)
[[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
for ligne in M: # affichage ligne par ligne
print(ligne)
[1, 2, 3, 4] [5, 6, 7, 8] [9, 10, 11, 12]
len(M) # nombre de lignes
3
len(M[0]) # nombre de colones
4
M[1][2] # coefficient ligne 2 colonne 3
7
On peut créer une matrice de 0 de taille $H\times L$ avec l'instruction
M = [ [ 0 for j in range(L) ] for i in range(H) ]
M = [ [ 0 for j in range(5) ] for i in range(4) ]
for ligne in M:
print(ligne)
[0, 0, 0, 0, 0] [0, 0, 0, 0, 0] [0, 0, 0, 0, 0] [0, 0, 0, 0, 0]
Le but de ce TP est de détecter un contour dans une image en utilisant l'algorithme de Dijkstra. L'idée est que l'utilisateur doit désigner deux pixels dans une image et l'algorithme doit trouver un contour qui permet de relier ces deux points.
Pour plus de simplicité, on se limite à des images en niveaux de gris, ainsi chaque pixel de l'image est identifié avec un entier entre 0 et 255 qui représente son niveau de gris.
Pour rappel, l'algorithme de Dijkstra permet, dans un graphes à poids positifs ou nuls, de calculer les chemins de poids minimal d'un sommet donné vers tous les autres sommets du graphe.
On donne les fonctions suivantes permettant:
from PIL import Image
def img2array(fic):
"""fic est le nom d'un fichier contenant une image
la fonction renvoie une liste de listes d'entiers représentants des niveaux de gris"""
im = Image.open(fic).convert('L') # chargement et passage en niveaux de gris
tab = [ [0 for _ in range(im.height)] for _ in range(im.width)]
i = 0
j = -1
for c in im.getdata():
if i % im.width == 0: # debut d'une nouvelle ligne dans l'image
i = 0
j += 1
tab[i][j] = c
i += 1
im.close()
return tab # variable locale donc effacée après exécution
# on ouvre l'image plage.jpeg et on l'enregistre dans la matrice tableau
tableau = img2array('plage.jpeg')
len(tableau) # nombre de lignes de l'image
320
len(tableau[0]) # nombre de colonnes de l'image
240
Fonction d'affichage:
from PIL import Image
from IPython.display import display
def affiche_image(img):
"""img : image sous forme d'une liste de listes de niveaux de gris
affichage de l'image"""
im = Image.new(size=(len(img), len(img[0])), mode='L')
for i in range(len(img)):
for j in range(len(img[i])):
im.putpixel((i, j), img[i][j])
display(im)
# on affiche l'image précédente donnée par la matrice tableau
affiche_image(tableau)
Une image est un tableau bidimensionnel de pixels. Ainsi à toute image on peut associer un graphe non orienté dont les sommets sont les pixels de l'image et tel que deux sommets sont adjacents si et seulement si les pixels correspondants sont voisins dans l'image. En dehors des sommets représentant des pixels sur le bord de l'image, chaque sommet possède donc 8 voisins.
La présence d'une arête entre deux sommets signifie simplement que les deux pixels sont voisins.
Voici ce graphe pour une image de taille $4\times 6$.
Une ligne sur une image est une suite de pixels voisins. Cette notion correspond donc à celle de chemin dans le graphe associé à l'image.
Sur l'image de taille $4\times 6$ précédente voici la ligne formée des pixels $(2,1)$, $(1,2)$, $(1,3)$, $(2,4)$, $(3,3)$, $(4,4)$, $(4,5)$ et $(4,6)$.
Intuitivement, un contour dans une image représente une suite de pixels qui marquent une variation brusque de couleur entre deux régions (ou une variation brusque d'intensité lumineuse dans le cas d'une image en noir et blanc).
Pour que le chemin de poids minimal entre deux sommets représente un contour, il faut bien choisir les poids sur les arêtes du graphe. Nous allons commencer par prendre les poids les plus naïfs possibles pous simplifier l'implémentation des algorithmes.
On identifie une image avec une matrice d'entiers de taille $H\times L$ (où on numérote les lignes et les colonnes à partir de 0 comme en Python): $$M=(m_{i,j})_{\genfrac{}{}{0pt}{}{0\leq i< H}{0\leq j< L}}$$
Si P = (a, b)
est un pixel de l'image ($0\leq a < H$ et $0\leq b < H$) sa couleur est donc donnée par M[a][b]
.
Le poids d'une arête d'un pixel $P=(a,b)$ vers un pixel $Q(c,d)$ est pour le moment défini comme la valeur absolue de la différence entre les couleurs de ces deux pixels: $$\big|\;M[a][b] - M[c][d]\; \big|$$
Écrire une fonction poids(M, P, Q)
qui renvoie le poids de l'arête entre les pixels P
et Q
dans l'image M
.
# Cellule à compléter
def poids(M, P, Q):
# Tests unitaires
M = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
P = (1,2)
Q = (2,3)
print( poids(M,P,Q) == 5 )
# Tests libres
Écrire une fonction voisins
qui prend en argument une matrice représentant une image M
(pour ses dimensions) et un pixel P
(sous la forme d'un tuple d'entiers) et renvoie la liste des voisins du pixel dans l'image.
# cellule à compléter
def voisins(M, P):
# Tests unitaires
from collections import Counter
# pour tester l'égalité de deux listes sans tenir compte de l'ordre des éléments
M = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
print( Counter( voisins(M, (2,3)) ) == Counter( [(1, 2), (1, 3), (2, 2)] ) )
# cellule pour tests libres
Pour implémenter l'algorithme de Dijkstra, nous allons utiliser des files de priorité.
Une file de priorité (priority queue en anglais) est une structure de donnée linéaire où chaque élément de la file possède une priorité. Intuitivement, c’est avec une telle structure de données que vous fonctionnez lorsqu’on vous donne du travail. Supposons que vous ayez plusieurs tâches à accomplir, listées par importance décroissante : travailler le cours de maths du jour, finir le DM de physique pour la semaine prochaine et faire un gros score à flappy bird. On se représente cet ensemble de tâches comme ceci
la flèche nous indiquant que la prochaine tâche à accomplir est de relire votre cours de maths. Les nombres flottants placés en dessous de chaque tâche sont des priorités : plus ce nombre est bas, plus la tâche est prioritaire. Si vous apprenez que votre prochaine colle de français a lieu dans trois jours, il faut vous ménager du temps pour la préparer. Vous insérez donc cette nouvelle tâche dans la file avec une priorité de $3.0$.
Le module heapq
fournit tout ce qu'il faut pour cela:
# Cellule à exécuter
import heapq
#help(heapq)
filep = [] # Pour créer une file de priorité vide
n = len(filep) # Pour connaitre la longueur de la file
p = 1
x = 'Cours de Maths'
heapq.heappush( filep, (p, x)) # Pour insérer x de priorité p
heapq.heappush( filep, (7., 'DM de Physique') ) # Pour insérer x de priorité p
heapq.heappush( filep, (1000., 'Flappy Bird') ) # Pour insérer x de priorité p
print('\n') # saut de ligne
print(filep)
heapq.heappush( filep, (3., 'Colle de Français') ) # Pour insérer x de priorité p
print('\n') # saut de ligne
print(filep)
p, x = heapq.heappop( filep ) # Pour défiler un élément de priorité MINIMALE
print('\n') # saut de ligne
print('p =',p)
print('x =', x)
[(1, 'Cours de Maths'), (7.0, 'DM de Physique'), (1000.0, 'Flappy Bird')] [(1, 'Cours de Maths'), (3.0, 'Colle de Français'), (1000.0, 'Flappy Bird'), (7.0, 'DM de Physique')] p = 1 x = Cours de Maths
Pour mettre un pixel $P$ avec un poids $p$, on utilise donc un tuple $(p, P)$ car les couples sont ordonnés en fonction de leur première composante.
Dans toutes les fonctions qui suivent, les paramètres qui apparaissent sont les suivants, si M
est la matrice associée à une image:
P
, Q
: points (sous le forme de couples d'entiers)dist
: tableau des distances calculéespred
: tableau des prédécesseurs calculésfile
: file d'attenteinitialisation
qui prend en argument une matrice représentant une image et renvoie un couple de deux tableaux bidimensionnels de mêmes dimensions que l'image, le premier contenant des math.inf
dans toutes ses cases (ce sont les distances initiales), et le deuxième des None
dans toutes ses cases (ce sont les prédécesseurs initiaux).# cellule à compléter
import math
def initialisation(M):
# Tests unitaires
import random
c1, c2 = random.randint(10, 20), random.randint(10, 20)
# nombres de lignes et de colonnes aléatoires entre 10 et 20
dist, pred = initialisation([[0]*c1]*c2)
# image toute noire (couleur = 0)
print(len(dist)==c2 and len(dist[0])==c1 and len(pred)==c2 and len(pred[0])==c1)
print(dist[c1//2][c2//2] == math.inf and pred[c1//2][c2//2] == None)
# cellule pour tests libres
mise_a_jour_file(M, P, Q, dist, pred, file)
qui teste si la longueur du enregistré dans dist
pour aller à Q
est plus long que la longueur du dhemin pour aller à P
enregistré dans dist
auquel on ajoute le poids de l'arête de P
à Q
. Si c'est le cas, la fonction met à jour dist
et pred
et insère dans la file de priorité le point Q
avec la nouvelle distance qui lui est associée (on peut ignorer les éventuelles autres instances de Q
dans la file car les valeurs qui leur sont associées sont supérieures).# cellule à compléter
def mise_a_jour_file(M, P, Q, dist, pred, file):
poidsArete = ......
if dist[Q[0]][Q[1]] > dist[P[0]][P[1]] + poidsArete:
dist[Q[0]][Q[1]] = ......
pred[Q[0]][Q[1]] = ......
heapq.heappush( ...... , ...... )
Dijkstra(M, S, T)
qui implémente l'algorithme de Dijkstra. Pour accélérer: on donne un point initial S
et un point final T
. Quand le point final est extrait de la file de priorité, on s'arrête.import heapq
def Dijkstra(M, S, T):
dist, pred = initialisation(M)
dist[S[0]][S[1]] = ......
file = ......
heapq.heappush( ...... , ...... )
while len(file) > 0:
p, P = ......
if p == math.inf or P == T:
return ......
v = voisins(M, P)
for Q in v:
............
return pred
contour(M, P, Q, couleur)
qui ajoute le contour calculé entre les points P
et Q
, dans la couleur spécifiée (entier entre 0 pour noir et 255 pour blanc).def contour(M, P, Q, couleur):
pred = Dijkstra(M, P, Q)
p = Q
while p != None:
M[p[0]][p[1]] = ......
p = ......
On peut maintenant regarder le résultat du calcul du contour sur l'image suivante:
img = img2array("plage.jpeg")
affiche_image(img)
contour(img, (281,57), (269,235), 255) # un peu long, c'est normal
affiche_image(img)
Comme on pouvait s'y attendre on voit que le contour est complètement délirant ! On va donc decvoir utiliser des poids plus sophistiqués.
Pour que le chemin de poids minimal entre deux sommets représente un contour, il faut bien choisir les poids sur les arêtes du graphe. Nous allons prendre les poids suggérés par l'article Intelligent scissors for image composition de Eric N. Mortensen et William A. Barrett.
Pour calculer ces poids, il faut utiliser l'opération de convolution sur des matrices.
Si on identifie une image avec une matrice d'entiers de taille $H\times L$ (où on numérote les lignes et les colonnes à partir de 0 comme en Python): $$M=(m_{i,j})_{\genfrac{}{}{0pt}{}{0\leq i< H}{0\leq j< L}}\:,$$ on peut lui appliquer un filtre (matrice de taille $(2I+1)\times (2I+1)$ où $I$ est un entier naturel) $$F = (f_{i, j})_{\genfrac{}{}{0pt}{}{0\leq i\leq 2I}{0\leq j\leq 2I}}$$ par produit de convolution. Le produit de convolution de $M$ et de $F$ noté $M\star F$ est une matrice de taille $H\times L$ dont les coefficients sont donnés par: $$(M\star F)_{k,\ell} = \sum_{i=0}^{2I}\sum_{j=0}^{2I} m_{k+i-I, \ell+j-I}\times f_{i, j}$$ (Bien sûr on se limite aux coefficients de $M$ qui ont un sens dans la somme ci-dessus. En général $I$ est nettement plus petit que $H$ et $L$, cette limitation ne concerne donc qu'un mince bord de l'image.)
convolution
qui prend en argument une matrice et un filtre et renvoie le produit de convolution des deux. (À noter que dans la vraie vie, on utilise une transformée de Fourier rapide pour cette opération, ce qui permet de gagner deux ordres de grandeur en complexité; ici nous ferons le calcul naïvement et cela explique la lenteur sur de vraies images.)#cellule à compléter
def convolution(M, F):
H = len(M)
L = len(M[0])
I = len(F)//2
N = [ [ 0 for j in range(L) ] for i in range(H) ]
for k in range(H):
for l in range(L):
somme = 0
for i in range(2*I+1):
for j in range(2*I+1):
if ............ :
somme = ............
N[k][l] = ......
return N
#cellule pour vérifier
F = [[0, 0, 0], [0, 1, 0], [0, 0, 0]]
M = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
print(convolution(M, F) == M)
F = [[-1,0,1], [1,0,-1], [0,1,-1]]
print(convolution(M, F) == [[-3, -3, -3, 11], [-5, -1, -1, 16], [-4, 0, 0, 4]])
# cellule pour tests libres
norme(X)
qui prend en argument un tuple X=(x,y)
de réels représentant un vecteur et renvoie sa norme euclidienne.# cellule à compléter
def norme(P):
# cellule pour tester
import math
epsilon = .000001
print( math.fabs(norme((0,0))) < epsilon )
print( math.fabs(norme((4,7)) - 8.06225774829855) < epsilon )
print( math.fabs(norme((-5,8)) - 9.433981132056603) < epsilon )
# cellule pour tests libres
produit_scalaire(X1,X2)
qui prend en argument deux tuples de réels représentant des vecteurs et renvoie le produit scalaire de ces deux vecteurs.# cellule à compléter
def produit_scalaire(X1, X2):
# cellule pour tester
print(produit_scalaire((1,3),(-2,9)) == 25)
# cellule pour tests libres
Pour déterminer le poids d'une arête on va utiliser des convolutions avec des matrices particuilères définies ci-dessous.
Laplacien. Le laplacien est un opérateur qui mesure la différence entre l'intensité lumineuse d'un pixel et la valeur moyenne des intensités lumineuses des pixels autour. Il peut être approché en convoluant l'image avec le filtre
Laplacien = [
[-1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1],
[-1, -1, 24, -1, -1],
[-1, -1, -1, -1, -1],
[-1, -1, -1, -1, -1]
]
Gradient. Le gradient est une généralisation vectorielle de la dérivée. La dérivée d'une fonction définie sur ${\mathbb R}$ (c'est-à-dire une ligne) nous donne des renseignements sur les variations de cette fonction, le gradient nous donne des renseignements sur les variations d'une fonction définie sur un plan. On peut considérer qu'une image en niveaux de gris est une fonction qui à chaque point associe comme valeur l'intensité lumineuse de ce point.
Le gradient est dirigé perpendiculairement aux lignes de niveau, vers les valeurs croissantes.
La composante horizontale (suivant l'axe des abscisses) du gradient peut être approchée par une convolution avec le filtre
Gx = [
[1, 0, -1],
[2, 0, -2],
[1, 0, -1]
]
La composante verticale (suivant l'axe des ordonnées) du gradient peut être approchée par une convolution avec le filtre transposé.
Gy = [
[1, 2, 1],
[0, 0, 0],
[-1, -2, -1]
]
Dans toutes les fonctions qui suivent, les paramètres qui apparaissent sont les suivants, si M
est la matrice associée à une image:
laplM
: laplacien de M
gradMx
: composante horizontale du gradient de M
gradMy
: composante vertitale du gradient de M
gradMax
: maximum des normes des gradients de M
P
, Q
: pixels (sous la forme de tuples d'entiers)gradMx
s'obtient avec l'instruction gradMx = convolution(M, Gx)
et de même gradMy
s'obtient avec l'instruction gradMy = convolution(M, Gy)
.
grad_max
qui prend en argument la composante horizontale gradMx
et la composante verticale gradMy
du gradient d'une matrice, et renvoie gradMAx
le maximum de la norme du gradient sur tous les points de l'image correspondante.# Cellule à compléter
def grad_max(gradMx, gradMy):
# Tests unitaires
import math
epsilon = 0.000001
M = [[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]]
gradMx, gradMy = convolution(M, Gx), convolution(M,Gy)
print(math.fabs(grad_max(gradMx, gradMy) - 37.013511046643494) < epsilon)
# cellule pour tests libres
Pour calculer le poids d'une arête on donne les fonctions suivantes:
# Cellule à valider
def fZ(laplM, Q):
if laplM[Q[0]][Q[1]] < 0.00001:
return 0
else:
return 1
def fG(gradMx, gradMy, gradMax, Q):
i = Q[0]
j = Q[1]
X= ( gradMx[i][j], gradMy[i][j] )
return 1 - norme(X)/gradMax
def dX(gradMx, gradMy, P, Q, X):
epsilon = 0.000000001
gradX = (gradMx[X[0]][X[1]], gradMy[X[0]][X[1]])
norme_gradX = norme(gradX)
if norme_gradX < epsilon:
return 0
DX = (gradX[1]/norme_gradX, -gradX[0]/norme_gradX)
gradP = (gradMx[P[0]][P[1]], gradMy[P[0]][P[1]])
DP = (gradP[1], -gradP[0])
normePQ = norme((Q[0]-P[0], Q[1]-P[1]))
LPQ = ((Q[0]-P[0])/normePQ, (Q[1]-P[1])/normePQ) if produit_scalaire(DP, [Q[0]-P[0], Q[1]-P[1]])>=0\
else ((P[0]-Q[0])/normePQ, (P[1]-Q[1])/normePQ)
return produit_scalaire(DX, LPQ)
def fD(gradMx, gradMy, P, Q):
return 2/(3*math.pi)*(math.acos(dX(gradMx, gradMy, P, Q, P)) + math.pi - math.acos(dX(gradMx, gradMy, P, Q, Q)))
def poids(laplM, gradMx, gradMy, gradMax, P ,Q):
return 0.43*fZ(laplM, Q) + 0.43*fG(gradMx, gradMy, gradMax, Q) + 0.14*fD(gradMx, gradMy, P, Q)
Seule la fonction poids(laplM, gradMx, gradMy, gradMax, P ,Q)
sera utilisée.
Pour accélerer les calculs, au lieu de transmettre l'image M
en argument des fonctions, on transmet seulement les valeurs liées au laplacien et au gradient de l'image.
Dans toutes les fonctions qui suivent, les paramètres qui apparaissent sont donc les suivants, si M
est la matrice associée à une image:
laplM
: laplacien de M
gradMx
: composante horizontale du gradient de M
gradMy
: composante vertitale du gradient de M
gradMax
: maximum des gradients de M
P
, Q
: points (sous le forme de couples d'entiers)dist
: tableau des distances calculéespred
: tableau des prédécesseurs calculésfile
: file d'attentemise_à_jour_file_v2(laplM, gradMx, gradMy, gradMax, P, Q, dist, pred, file)
qui teste si le chemin donné par dist
pour aller à Q
est plus long que le chemin pour aller à P
donné par dist
auquel on ajoute le poids de l'arête de P
à Q
. Si c'est le cas, elle met à jour dist
et pred
et insère dans la file de priorité $Q$ avec la nouvelle distance qui lui est associée (on peut ignorer les éventuelles autres instances de Q
dans la file car les valeurs qui leur sont associées sont supérieures).# cellule à compléter
def mise_a_jour_file_v2(laplM, gradMx, gradMy, gradMax, P, Q, dist, pred, file):
Dijkstra_v2(M, S, T)
qui implémente l'algorithme de Dijkstra. Pour accélérer: on donne un point initial S
et un point final T
. Quand le point final est extrait de la file de priorité, on s'arrête.import heapq
def Dijkstra_v2(M, S, T):
laplM = convolution(M, Laplacien)
gradMx, gradMy = convolution(M, Gx), convolution(M, Gy)
gradMax = grad_max(gradMx, gradMy)
dist, pred = ......
dist[S[0]][S[1]] = ......
file = ......
heapq.heappush( ...... , ...... )
while len(file) > 0:
p, P = ......
if P == T:
return ......
v = ......
for Q in v:
............
return ......
contour_v2(M, P, Q, couleur)
qui ajoute le contour calculé entre les points P
et Q
, dans la couleur spécifiée (entier entre 0 pour noir et 255 pour blanc).def contour_v2(M, P, Q, couleur):
On peut maintenant regarder le résultat du calcul du contour sur l'image suivante:
img = img2array("plage.jpeg")
affiche_image(img)
contour_v2(img, (281,57), (269,235), 255) # un peu long, c'est normal
affiche_image(img)