Ce notebook a pour but de présenter une utilisation possible de l'algorithme OMP. D'un point de vue théorique, l'algorithme OMP vise à trouver la meilleure approximation sparse d'un vecteur $x$ à partir d'un dictionnaire $D$ donné. Mathématiquement, on cherche à résoudre en $\gamma$ le problème d'optimisation suivant: $$\min_\gamma \left||x - D\gamma \right||_{2}^{2}$$ sous la contrainte que tous les éléments de $\gamma$ sont soit nuls soit égaux à 1 et que: $$\sum_{i = 1}^{n} \gamma _i \leq K$$ avec $K$ le nombre maximal de features du dictionnaire $D$ présentes dans le vecteur $x$.
Dans notre cas, le vecteur $x$ sera représenté par une image pouvant contenir plusieurs formes géométriques. L'ensemble des formes géométriques se stocké dans le dictionnaire et nous verront comment les retrouver avec l'implémentation de l'algorithme OMP que propose scikit-learn.
Commençons, à l'aide de la bibliothèqe scikit-image, par générer quelques formes géométriques...
In [1]:
import numpy as np
from skimage.draw import circle_perimeter, ellipse_perimeter, polygon, line
In [2]:
cercle = np.zeros((10, 10), np.uint8)
cercle[circle_perimeter(4, 4, 3)] = 1
ellipse = np.zeros((10, 10), np.uint8)
ellipse[ellipse_perimeter(4, 4, 3, 5)] = 1
square = np.zeros((10, 10), np.uint8)
square[polygon(np.array([1, 4, 4, 1]), np.array([1, 1, 4, 4]))] = 1
dline = np.zeros((10, 10), np.uint8)
dline[line(1, 1, 8, 8)] = 1
shapes = [cercle, ellipse, square, dline]
Une petite représentation graphique de l'une des formes ci-dessus...
In [3]:
import matplotlib.pyplot as plt
%matplotlib inline
In [4]:
plt.imshow(cercle, cmap='gray')
plt.title("Un cercle obtenu avec scikit-image")
Out[4]:
Il ne reste plus qu'à constituer le dictionnaire:
In [5]:
D = np.zeros((100, 4))
for i in range(len(shapes)):
D[:,i] = shapes[i].ravel()
D[:,i] = D[:,i]/np.linalg.norm(D[:,i])
On définit une fonction permettant de générer une combinaison aléatoire des signaux images regroupés dans notre dictionnaire.
In [6]:
import random
In [7]:
def generate_img(i):
"""generates combination of i shapes from D dictionnary"""
output = np.zeros(100)
for _ in range(i):
output += D[:,random.randint(0, len(shapes)-1)]
return output
Quelques démonstrations...
In [8]:
plt.imshow(generate_img(1).reshape((10, 10)), cmap='gray')
#dline
Out[8]:
In [10]:
plt.imshow(generate_img(2).reshape((10, 10)), cmap='gray')
#ellispe + carré
Out[10]:
In [11]:
from sklearn.linear_model import orthogonal_mp
A partir du moment où l'on a identifié les signaux utilisés dans l'image à analyser, on peut en déduire les formes géométriques présentes dans l'image. C'est ce que fait la fonction ci-dessous.
In [12]:
shapes_names = np.array(['circle', 'ellipse', 'square', 'line'])
def recognize(img, n_coefs):
gamma = orthogonal_mp(D, img, n_nonzero_coefs = n_coefs)
plt.imshow(img.reshape((10, 10)), cmap = 'gray')
return shapes_names[gamma.astype('bool')]
Quelques exemples pour finir... Remarque: parfois la fonction génératrice d'images peut choisir deux fois le même signaux, on peut donc avoir un nombre de formes inférieur à $n$.
In [13]:
n = 1
recognize(generate_img(n), n)
Out[13]:
In [14]:
n = 2
recognize(generate_img(n), n)
Out[14]:
In [15]:
n = 3
recognize(generate_img(n), n)
Out[15]:
Grâce à l'algorithme OMP, on est bien en mesure de retrouver les signaux utilisés pour construire une image. Tout le challenge consiste donc pour nous, à partir d'un certain nombre d'images à construire un dictionnaire de signaux aussi pertinent que possible. C'est là qu'intervient l'algorithme K-SVD !
In [16]:
import imp
ksvd = imp.load_source('ksvd', '../Source/ksvd.py')
On construit une grosse matrice avec pleins d'images comme ci-dessus mais toujours avec deux symboles pour corser un peu les choses !
In [20]:
X = []
for _ in range(500):
X.append(generate_img(2))
X = np.array(X).T
On construit à partir de cette matrcie un dictionnaire grâce à notre implémentation de la classe KSVD
In [24]:
model = ksvd.KSVD((100, 4), K = 2, n_iter=400)
model.fit(X)
Et on représente les atoms du dictionnaire obtenu!
In [25]:
plt.figure(figsize=(12, 12))
for i in range(4):
plt.subplot(2, 2, i+1)
plt.imshow(model.D[:,i].reshape((10,10)), cmap = 'gray')
plt.title(u"Feature n°{}".format(i+1))
Et quelques exemples de la façon dont on parvient à reconstituer des images à l'aide de leur décomposition sparse, issue du dictionnaire préceddement obtenu.
In [23]:
gamma = model.sparse_rep(X)
plt.figure(figsize=(10, 25))
for i in range(5):
plt.subplot(5, 2, i*2 + 1)
plt.imshow(X[:,i].reshape((10,10)), cmap = 'gray')
plt.title("Original image")
plt.subplot(5, 2, i*2 + 2)
plt.imshow(model.D.dot(gamma[:,i]).reshape((10,10)), cmap = 'gray')
plt.title("Sparse representation")
Trop beau pour être vrai? L'algorithme retrouve parfaitement les images avec lesquelles a été constituée la matrice donc nous avons triché? Un peu... Essayons par souci d'honnêteté intellectuelle de refaire la même opération mais en mettant 2 et oarfois 3 des images de base dans les images de la matrice:
In [26]:
X = []
for _ in range(250):
X.append(generate_img(2))
for _ in range(250):
X.append(generate_img(3))
X = np.array(X).T
In [27]:
model = ksvd.KSVD((100, 4), K = 3, n_iter=500)
model.fit(X)
In [28]:
plt.figure(figsize=(12, 12))
for i in range(4):
plt.subplot(2, 2, i+1)
plt.imshow(model.D[:,i].reshape((10,10)), cmap = 'gray')
plt.title(u"Feature n°{}".format(i+1))
Ce n'est pas si mal mais c'est un peu moins bien (on peut dire que le problème était sans doute plus difficile). Voyons tout de même comment sont reconstituées les images:
In [29]:
gamma = model.sparse_rep(X)
plt.figure(figsize=(10, 25))
for i in range(5):
plt.subplot(5, 2, i*2 + 1)
plt.imshow(X[:,i].reshape((10,10)), cmap = 'gray')
plt.title("Original image")
plt.subplot(5, 2, i*2 + 2)
plt.imshow(model.D.dot(gamma[:,i]).reshape((10,10)), cmap = 'gray')
plt.title("Sparse representation")
Très correct non ?
Tant que j'y pense, choses qu'il pourrait être intéressant de faire avec les tests sur les images générées: