TP : Classification de tableaux

---

Dans ce TP, nous voulons classifier des tableaux de Claude Monet en utilisant l’algorithme des k-moyennes.

Le volume de données étant relativement important, je vous invite à travailler en local sur vos machines. Commençons par télécharger les données ici

import numpy as np
from matplotlib import pyplot as plt
import glob

file = '*.jpg' 
glob.glob(file)
# Using List Comprehension to read all images
#images = np.array([plt.imread(image) for image in glob.glob(file)])
images = [plt.imread(image) for image in glob.glob(file)]

images est une liste contenant tous nos tableaux de Monet.

images[0] # la 1ère matrice

array([[[  0, 110, 124],
        [ 18, 134, 147],
        [ 29, 135, 148],
        ...,
        [  9, 115, 157],
        [  8, 113, 158],
        [ 26, 131, 176]],

       [[ 25, 137, 149],
        [ 31, 139, 151],
        [ 22, 123, 133],
        ...,
        [  4, 110, 152],
        [ 12, 117, 162],
        [  8, 113, 158]],

       [[ 52, 151, 157],
        [ 51, 146, 152],
        [ 41, 130, 134],
        ...,
        [  8, 114, 156],
        [ 23, 130, 174],
        [  2, 109, 153]],

       ...,

       [[183, 164,  98],
        [191, 173, 111],
        [183, 166, 114],
        ...,
        [179, 161, 137],
        [121, 124, 115],
        [152, 166, 166]],

       [[173, 161,  75],
        [184, 171,  92],
        [150, 136,  71],
        ...,
        [183, 145, 124],
        [147, 131, 132],
        [165, 159, 173]],

       [[175, 166,  71],
        [227, 217, 128],
        [211, 201, 129],
        ...,
        [201, 153, 133],
        [215, 188, 195],
        [186, 171, 192]]], dtype=uint8)

Chacune de ces matrices est en fait une image au format .jpg contenant chaque pixel d’un tableau de Monet.
Chaque élément de la matrice représente un pixel coloré représenté par trois valeurs entières entre 0 et 255 : le niveau de rouge, vert et bleu. Par exemple, affichons l’image correspondant à la matrice d’indice 0 :

plt.imshow(images[0])
plt.show()

Redimensionner les images

Nous avons tout d’abord besoin de redimensionner nos images à la même taille. Une manière simple de s’y prendre consiste à couper nos images, de sorte que la taille finale de chacune corresponde à la hauteur minimale parmi les images du jeu de données par la largeur minimale trouvée.

Question : Définir une fonction crop_to_min standardisant la taille des images.

def crop_to_min(images):
    #recherche de minimums
    h_min = len(images[0])
    l_min = len(images[0][0])
    for i, image in enumerate(images):
        #Pas utile à la fonction mais permet de detecter quelles images sont "trop grandes"
        if h_min < len(image) or l_min < len(image[0]):
            print("Trop grande:", i)
        if h_min > len(image):
            h_min = len(image)       
        if l_min > len(image[0]):
            l_min = len(image[0])
    print("Taille min:",h_min, l_min)
    #Découpage des images
    l=[]
    for image in images:
        l.append(image[:h_min,:l_min])
    return l

digits_matrix = crop_to_min(images)
plt.imshow(images[296])
plt.figure()
plt.imshow(digits_matrix[296])
plt.show() 

Trop grande: 296
Taille min: 256 256

Pour appliquer l’algorithme des \(k\)-moyennes, il faut que les données soient des vecteurs. On transforme donc une matrice en vecteur en concaténant les lignes de la matrice :
$\( \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \longrightarrow \begin{pmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \end{pmatrix} \)$

Question : Écrire une fonction to_vector(m) qui prend en argument une matrice m et renvoie le vecteur correspondant, sous forme de liste.

def to_vector(m):
    v = []
    for i in range(0,len(m)):
        for j in range(0,len(m[0])):
            v.append(m[i][j])
    return v

print(len(to_vector(digits_matrix[0])))
assert len(to_vector(digits_matrix[0])) == len(digits_matrix[0])*len(digits_matrix[0][0]) # les images 256*256*3 sont transformées en vecteur 65536*3

65536

Question : Écrire une fonction to_matrix(v) qui prend en argument un vecteur v et renvoie la matrice correspondante, sous forme de liste de listes. On suppose que la matrice est carrée.

def to_matrix(vector):
    mat = [[[0,0,0] for i in range(0, 256)] for j in range(0, 256)]
    for i in range(0, 256):
        for j in range(0, 256):
            mat[i][j] = vector[i*256+j]
    return mat

assert (to_matrix(to_vector(digits_matrix[0])) == digits_matrix[0]).all() 
# to_matrix et to_vector sont inverses l'une de l'autre

Question : Définir une liste X contenant les vecteurs correspondant aux images dans digits_matrix. Cette variable X sera utilisée dans toute la suite.

X = []
for i in range(0, len(digits_matrix)):
    X.append(to_vector(digits_matrix[i]))

Algorithmes des \(k\)-moyennes

Question : Écrire une fonction d(u, v) qui prend en argument deux images u et v et renvoie la distance euclidienne entre ces deux vecteurs.

On risque d’obtenir un warning lors du calcul de soustraction (overflow encountered in scalar subtract), dans ce cas vérifier le type des données avec type(ma_donnée). Que remarque-t-on ? Proposer une modification de la fonction pour éviter ce problème.

def d(u, v):
    r = np.int64(0)
    g = np.int64(0)
    b = np.int64(0)
    #print(type(u[0][0]))
    for i in range(0, len(u)):
        r += (np.int64(u[i][0]) - np.int64(v[i][0]))**2
        g += (np.int64(u[i][1]) - np.int64(v[i][1]))**2
        b += (np.int64(u[i][2]) - np.int64(v[i][2]))**2
    return (r + g + b)**.5

def d(u, v):
    u = np.asarray(u, dtype=np.int64)
    v = np.asarray(v, dtype=np.int64)
    diff = u - v                       # vectorized int64 subtraction
    sq = diff * diff                   # elementwise square
    return np.sqrt(sq.sum())           # final Euclidean distance

d(X[0], X[1])

np.float64(28317.940320581227)

Comme dans le cours, nous utilisons une liste centres pour stocker les centres des classes. centres[i] est le centre de la classe \(i\).

Question : Écrire une fonction centres_aléatoires(X, k) qui prend en argument une liste de vecteurs X et un entier k et renvoie une liste de k vecteurs aléatoires choisis dans X.
On utilisera random.sample dont on pourra consulter la documentation ici. On pourra préalablement utiliser random.seed(0) pour avoir les mêmes résultats que le corrigé du TP.

import random

random.seed(0) # pour avoir toujours les mêmes résultats
def centres_aléatoires(X, k):
    return random.sample(X, k)

centres = centres_aléatoires(X, 5)
assert len(centres) == 5

Question : Écrire une fonction plus_proche(x, centres) qui prend en argument un vecteur x et une liste de vecteurs centres et renvoie l’indice du centre le plus proche de x.

def plus_proche(x, centres):
    return min(range(len(centres)), key=lambda i: d(x, centres[i]))

plus_proche(X[1], centres)

3

Contrairement à ce qui a été fait dans le cours, on utilise dans la suite une liste classes telle que classes[i] est le numéro de classe de X[i] (dans le cours, classes[i] était la liste des vecteurs de la classe \(i\)).

Question : Écrire une fonction calculer_classes(X, centres) qui renvoie la liste classes correspondant à X et à centres. Il faut donc que classes[i] soit l’indice du centre (dans centres) le plus proche de X[i].

def calculer_classes(X, centres):
    classes = []
    for x in X:
        classes.append(plus_proche(x, centres))
    return classes

classes = calculer_classes(X[:20], centres_aléatoires(X[:20], 5))
#assert len(classes) == len(X) and max(classes) < 10 and min(classes) >= 0

print(classes)

[2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 4, 1, 2, 4, 0, 2, 4, 2, 3]

Question : Écrire une fonction centre(X, classes, i) qui renvoie le centre de la classe \(i\).

def centre(X, classes, i):
    p = len(X[0])
    c = [[0]*3 for _ in range(p)]
    n = 0
    for j in range(len(X)):
        if classes[j] == i:
            n += 1
            for k in range(p):
                for l in range(3):
                    c[k][l] += np.float64(X[j][k][l])
    if n==0:
        return c
    for k in range(p):
        for l in range(3):
            c[k][l] = np.uint8(c[k][l]/n)
    return c

plt.imshow(to_matrix(centre(X[:20], classes, 2)))
plt.show() 

Question : Écrire une fonction calculer_centres(X, classes, k) qui renvoie la liste centres correspondant à X et à classes.

def calculer_centres(X, classes, k):
    centres = []
    for i in range(0, k):
        centres.append(centre(X, classes, i))
    return centres

#classes = calculer_classes(X[:20], centres_aléatoires(X, 5))
centres = calculer_centres(X[:20], classes, 5)

plt.imshow(to_matrix(centres[0]))
plt.show()
plt.imshow(to_matrix(centres[1]))
plt.show()
plt.imshow(to_matrix(centres[2]))
plt.show()
plt.imshow(to_matrix(centres[3]))
plt.show()
plt.imshow(to_matrix(centres[4]))
plt.show()

Question : Écrire une fonction kmeans(X, centres, k) qui, à partir de centres initiaux donnés en argument, applique l’algorithme des \(k\)-moyennes et renvoie un couple correspondant à la liste (centres, classes) obtenue à la fin de l’algorithme.

def kmeans(X, centres, k):
    classes = calculer_classes(X, centres)
    centres = calculer_centres(X, classes, k)
    if classes != calculer_classes(X, centres):
        return kmeans(X, centres, k)
    return centres, classes

k = 5
centres, classes = kmeans(X[:20], centres_aléatoires(X, k), k) # ceci prend quelques secondes
fig, axes = plt.subplots(5, figsize=(20, 20))
for ax, center in zip(axes.ravel(), centres):
    ax.matshow(np.array(to_matrix(center)))
    ax.axis('off')
plt.show()

Calcul de l’inertie

Dans cette partie, on réutilise les listes centres et classes obtenues précédemment par un appel à kmeans.

Question : Écrire une fonction inertie(X, centres, classes) qui renvoie l’inertie de la partition définie par classes. On rappelle que l’inertie d’une partition est la somme des distances au carré des points à leur centre de classe.

def inertie(X, centres, classes):
    s = 0
    for i in range(len(X)):
        s += d(X[i], centres[classes[i]])**2
    return s

inertie(X[:20], centres, classes)

np.float64(8000858827.0)

Question : Réécrire la fonction kmeans précédente de façon à ce qu’elle renvoie aussi la liste des inerties obtenues à chaque itération de l’algorithme. Utiliser le code ci-dessous pour afficher aussi l’évolution de l’inertie au cours des itérations. Commenter.

from copy import deepcopy
def kmeans(X, centres, k):
    inerties = []
    centres2 = deepcopy(centres)
    centres2[0][0]=257
    while not np.array_equal(centres ,centres2):
        centres2 = deepcopy(centres)
        classes = calculer_classes(X, centres2)
        centres = calculer_centres(X, classes, k)
        inerties.append(inertie(X, centres, classes))
    return centres, classes, inerties

k=4
centres, classes, inerties = kmeans(X[:20], centres_aléatoires(X[:20], k), k) # ceci prend quelques minutes
plt.plot(inerties)
plt.show()

inerties

[np.float64(8431736040.0),
 np.float64(7941687198.0),
 np.float64(7790241383.0),
 np.float64(7790241383.0)]

Choix de la valeur de k

Question : Proposer une méthode permettant de déterminer la valeur de k à choisir et la mettre en oeuvre

inerties=[]
max_k=20
c=centres_aléatoires(X[:20], max_k)
for k in range(2,max_k+1):
    centres, classes, i = kmeans(X[:20], c[:k], k) # ceci prend quelques minutes
    inerties.append(i[-1])

plt.plot(range(2,max_k+1),inerties)

[<matplotlib.lines.Line2D at 0x7a4361066240>]

Question : Quelle valeur de k à l’air pertinente pour notre echantillon de tableaux ?

Amélioration de la fonction de distance

import numpy as np
from scipy.ndimage import gaussian_filter

def d(u, v, window_sigma=1.5):
    """
    SSIM-based distance for RGB data (last dimension = 3).
    Returns a distance in [0, 1].
    """

    # Convert to float, keep shape (H,W,3) or (N,3)
    u = np.asarray(u, dtype=np.float64)
    v = np.asarray(v, dtype=np.float64)

    # dynamic range of values
    L = max(u.max() - u.min(), v.max() - v.min())
    K1, K2 = 0.01, 0.03
    C1 = (K1 * L) ** 2
    C2 = (K2 * L) ** 2

    # Compute SSIM for each channel separately
    ssim_channels = []
    for c in range(3):
        U = u[..., c]
        V = v[..., c]

        mu_u = gaussian_filter(U, window_sigma)
        mu_v = gaussian_filter(V, window_sigma)

        mu_u2 = mu_u * mu_u
        mu_v2 = mu_v * mu_v
        mu_uv = mu_u * mu_v

        sigma_u2 = gaussian_filter(U * U, window_sigma) - mu_u2
        sigma_v2 = gaussian_filter(V * V, window_sigma) - mu_v2
        sigma_uv = gaussian_filter(U * V, window_sigma) - mu_uv

        num = (2 * mu_uv + C1) * (2 * sigma_uv + C2)
        den = (mu_u2 + mu_v2 + C1) * (sigma_u2 + sigma_v2 + C2)

        ssim_value = np.mean(num / den)
        ssim_channels.append(ssim_value)

    # mean across R,G,B
    ssim_value = np.mean(ssim_channels)

    # similarity → distance
    return (1.0 - ssim_value)*10

inerties=[]
max_k=8
c=centres_aléatoires(X[:20], max_k)
for k in range(2,max_k+1):
    centres, classes, i = kmeans(X[:20], c[:k], k) # ceci prend quelques minutes
    inerties.append(i[-1])

plt.plot(range(2,max_k+1),inerties)

[<matplotlib.lines.Line2D at 0x7a435c93eb40>]

fig, axes = plt.subplots(8, figsize=(20, 20))
for ax, center in zip(axes.ravel(), centres):
    ax.matshow(np.array(to_matrix(center)))
    ax.axis('off')
plt.show()

classes

[5, 5, 5, 5, 6, 5, 5, 5, 2, 5, 5, 7, 5, 1, 3, 4, 6, 7, 0, 5]

Amélioration de la fonction crop

Précédement, nous avions utilisé une fonction qui découpait les images trop grandes en ne gardant que les pixels correspondant au rectangle en haut à gauche de l’image.

Question : Proposer une amélioration de la fonction crop qui prenne mieux en compte l’image dans sa globalité.