1CC1000 - Systèmes d'information et programmation - TD : Approfondissement de Numpy - Traitement photographique

Table des matières

L'objectif de ce TD est : (1) d'apprendre à utiliser numpy pour réaliser des calculs efficaces et (2) de mieux comprendre matplotlib.

Quand numpy est-il efficace ?

Liste de listes

Le code suivant :

crée une liste de listes A générée aléatoirement, contenant un grand nombre de 0 et 1 ;
crée une deuxième liste de listes, ayant la même taille que A, que nous appelons subtilement B ;
copie les valeurs de A vers B.

Nous pouvons mesurer le temps d'exécution du calcul des valeurs de la matrice B. Nous devrions obtenir un temps de l'ordre de quelques centièmes de secondes. Testez-le.

from random import randint
from timeit import timeit

A = [[randint(0,1) for _ in range(300)] for _ in range(600)]
B = [[0 for _ in range(300)] for _ in range(600)]

def copy(A, B):
    for i in range(0, len(A)):
        for j in range(0, len(A[0])):
            B[i][j] = A[i][j]

timeit("copy(A, B)", number=10, globals=globals())

Tableau Numpy

Le code suivant est utilisé pour effectuer une opération équivalente en utilisant numpy. Comparez son temps d'exécution avec celui obtenu précédemment. Cela vous surprend-il ?

import numpy as np

A = np.random.randint(0, 2, (300, 600))
B = np.zeros(A.shape)

def copy(A, B):
    for i in range(0, len(A)):
        for j in range(0, len(A[0])):
            B[i,j] = A[i,j]

timeit("copy(A, B)", number=10, globals=globals())

La bonne méthode

Il semblerait qu'en utilisant numpy nous obtenions de moins bonnes performances que si nous utilisions des listes. Le problème est que nous n'utilisons pas numpy correctement. En effet, cette bibliothèque est optimisée pour la "vectorisation". Il s'agit d'une technique consistant à appliquer des opérations à un ensemble de données à la fois.

Nous devons changer notre approche pour tirer pleinement parti de numpy. Par exemple, la manière correcte d'effectuer une opération de copie est la suivante. Vous devriez convenir qu'elle est nettement plus efficace.

A = np.random.randint(0, 2, (300, 600))
B = np.zeros(A.shape)

def copy(A, B):
    A[:,:] = B[:,:]

timeit("copy(A, B)", number=10, globals=globals())

Comment utiliser Numpy

Tranches de listes

Vous savez déjà comment découper des listes. Voici un exemple, on notera que :

lorsque le découpage est utilisé à gauche d'un signe égal, il est utilisé pour attribuer des valeurs dans une liste ;
lorsque le découpage est utilisé à droite d'un signe égal, il est utilisé pour dupliquer une partie d'une liste.

l = [0, 1, 2, 3, 4, 5]
print("Id of l  : ", id(l))
l[:3] = [2, 1, 0]
print("l :  ", l)

l2 = l[:4]
print("Id of l2 : ", id(l2))  #Two python objects with the same id are physically the same
l2[:] = [0, 0, 0, 0]          #Warning do not confuse with "l2 = [0, 0, 0]" which is used to change the reference l2
print("Id of l2 : ", id(l2))

print("l :  ", l)
print("l2 : ", l2)

Vues

L'utilisation de la même notation sur un tableau numpy ne crée pas de tranche. Au lieu de cela, cela crée une vue, qui est une autre façon de visualiser les données, mais ce n'est pas une copie des données.

Testez le code suivant.

L = np.array([0, 1, 2, 3, 4, 5])
print("Id of L  : ", id(L))
L[:3] = [2, 1, 0]
print("L :  ", L)

L2 = L[:4]
print("Id of L2 : ", id(L2))  #Two python objects with the same id are physically the same
L2[:] = [0, 0, 0, 0]          #Warning do not confuse with "L2 = [0, 0, 0]" which is used  to change reference L2
print("Id of L2 : ", id(L2))

print("L :  ", L)
print("L2 : ", L2)

Les tableaux et les vues peuvent souvent être utilisés indifféremment. Lorsqu'on manipule une vue, l'attribut base pointe vers le tableau correspondant.

print("Base of L :  ", L.base)
print("Base of L2 : ", L2.base)

print("Id of L : ", id(L))
print("Id of the base of L2 : ", id(L2.base))

Une vue est un lien vers les données et une organisation des données. Avec reshape, nous pouvons créer de nouvelles vues vers les données qui peuvent être organisées différemment.

L3 = np.reshape(L2, (2, 2))
L3[1, 0] = 9
print("L3 :\n", L3)
print("L :\n", L)

Indexation

Il est possible de référencer des indices comme dans les tranches (start, end, step) ; cependant, lorsqu'on manipule des tableaux numpy à plusieurs dimensions, il n'est pas nécessaire d'enchaîner les crochets pour accéder à des listes de listes (de listes...), on peut se contenter d'utiliser des virgules.

Enfin, on peut assigner les valeurs contenues dans une vue en utilisant une autre vue de même taille et de même dimension, ou un scalaire.

Pour plus d'informations : indexing.

A = np.ones((3, 4))
print("Array A :\n", A)

A[::2, :] = 2
print("Array A: modified with a scalar\n", A)

B = np.zeros((3, 4))
A[:,:] = B[:,:]
print("Array A: modified with a view\n", A)

Traitement photographique

Attention: Nous voulons répondre aux questions suivantes en utilisant la vectorisation. Il ne faut jamais utiliser explicitement des boucles. Nous fournissons des "fonctions utiles" pour chaque question ; elles sont nécessaires, mais pas toujours suffisantes.

Nous pouvons lire des images PNG avec matplotlib et renvoyer un tableau tridimensionnel contenant l'intensité de chaque couleur de chaque pixel. L'intensité est une valeur comprise entre 0 et 255.

La dimension 1 est l'abscisse d'un pixel ;
La dimension 2 est l'ordonnée d'un pixel ;
la dimension 3 est l'indice de couleur, au format RVB.

Nous avons choisi cette image de raton laveur car elle permet de visualiser correctement le traitement que nous allons décrire.
Téléchargez l'image et enregistrez-la dans votre répertoire de travail. Après avoir testé vos fonctions sur cette image, vous pouvez très bien utiliser d'autres images PNG.

Selon votre ordinateur (et les bibliothèques installées), imread() peut avoir un comportement légèrement différent. En particulier :

il peut retourner un tableau de valeurs entières (entre 0 et 255), ou un tableau de nombres à virgule flottante (entre 0. et 1.) ;
la dernière dimension du tableau peut être de dimension 3 ou 4 ; (dans ce dernier cas, chaque pixel est encodé au format RGBA : rouge, vert, bleu et transparence).

Le reste des questions suppose que vous disposez d'un tableau de nombres entiers et que la dernière dimension est 3. Si ce n'est pas le cas pour vous, veuillez décommenter les deux lignes correspondantes (regardez les commentaires) dans le code suivant :

import matplotlib.pyplot as plt
import numpy as np


img = plt.imread("raccoon.png")
#Convert to int
#img = img * 255  
#img = img.astype(int)

#Remove the transparency component
#img = img[:,:,0:3]
plt.imshow(img, norm=None)
plt.show()

Niveaux de gris

Conversion

Ecrivez une fonction grayscale(img) qui prend en paramètre une image couleur sous forme de tableau et renvoie un nouveau tableau de seulement 2 dimensions, contenant le niveau de gris de chaque pixel.

Nos yeux ne percevant pas l'intensité lumineuse de chaque couleur de la même manière, il est nécessaire d'appliquer les coefficients suivants pour la conversion : rouge 30%, vert 59% et bleu 11%.

Voici quelques fonctions / opérations utiles :

A * x, l'opération multiplie chaque valeur de A par le scalaire x ;
A.astype(int), crée une copie du tableau A contenant uniquement des valeurs arrondies ;
plt.imshow(img, cmap="gray") affiche une image en niveaux de gris.

ANSWER ELEMENTS

def grayscale(img):
    G = np.zeros(img.shape[:2])
    G = img[:,:,0] * 0.3 + img[:,:,1] * 0.59 + img[:,:,2] * 0.11
    return G.astype(int)

gray = grayscale(img)
plt.imshow(gray, cmap="gray", norm=None)
plt.show()

Histogramme

Ecrivez une fonction histogram(img) qui prend en paramètre une image en niveaux de gris et renvoie un histogramme des valeurs des intensités présentes dans l'image. La valeur de retour est donc un tableau numpy de 256 valeurs entre 0 et 1, représentant la proportion de pixels ayant une certaine valeur.

Voici quelques fonctions utiles :

A / x : l'opération consiste à diviser chaque valeur de A par le scalaire x ;
np.zeros(n) : crée un tableau à dimension unique de taille n ;
np.add.at(A, indices, x) : pour chaque valeur i dans le tableau indices, incrémente de x la i-th valeur de A ; C'est équivalent à quelque chose comme
```
for i in indices :
    A[i] += x
```
mais plus rapide.
np.sum(A) additionne la liste des éléments de A.

ANSWER ELEMENTS

def histogram(img):
    bins = np.zeros(256)
    np.add.at(bins, img, 1)
    s = np.sum(bins)
    return bins / s

histogram(gray)[:10] # print the first 10 values (avoid flooding the output).

Tracer l'histogramme

Ecrivez une fonction plot_intensity(img) qui affiche une image et l'histogramme des intensités côte à côte.

Voici quelques fonctions utiles :

np.linspace(0, 255, 256) crée un nouveau tableau avec 256 valeurs réparties uniformément entre 0 et 255 (inclus) ;
plt.fill_between(x, y) remplit la zone située entre deux courbes horizontales.
plt.subplots(figsize = (15, 5)) crée un graphique multiple d'une taille donnée ;
plt.subplot(nbr, nbc, num) indique l'emplacement de la prochaine courbe à tracer. Exemple de subplot : cliquez ici.

ANSWER ELEMENTS

def plot_intensity(img):
    bins = histogram(img)
    x = np.linspace(0, 255, 256)
    plt.subplots(figsize=(15, 5))
    plt.subplot(1, 2, 1)
    plt.title("Image")
    plt.imshow(img, cmap="gray", norm=None)
    plt.xticks([]), plt.yticks([])
    plt.subplot(1, 2, 2)
    #plt.yscale("log")
    #plt.ylim(10**-7, 1)
    plt.title("Intensity")
    plt.fill_between(x, bins)
    plt.show()

plot_intensity(gray)

Contraste

Une photo est considérée comme correctement exposée si les intensités sont bien réparties sur l'ensemble du spectre, en particulier, le noir doit être profond et le blanc très lumineux. De nombreux outils de retouche d'images permettent de choisir automatiquement le contraste et la luminosité pour obtenir une image équilibrée.

Une approche naïve pour réaliser cette opération consiste à transformer les 2% des intensités les plus faibles en noir absolu (valeur de 0) et les 2% des intensités lumineuses les plus fortes en un blanc saturé (valeur de 255). Les valeurs intermédiaires sont réparties linéairement.

La première étape consiste donc à trouver les valeurs {$p_2$} et {$p_{98}$} des 2ème et 98ème centiles. En d'autres termes, {$p_2$} est l'intensité telle que 2% des intensités totales sont inférieures. Les valeurs {$p_2$} et {$p_{98}$} sont comprises entre 0 et 255.

Une fois ces valeurs connues, il est possible de résoudre le système d'équations suivant, où {$C$} et {$I$} désignent les corrections de contraste et d'intensité : {$$C*p_2 + I = 0$$} {$$C * p_{98} + I = 255$$}

Une fois ce système résolu, il est possible de corriger chaque valeur {$i$} d'intensité de l'image en niveaux de gris par la nouvelle intensité : {$C*i + I$}.

Centiles

Ecrivez une fonction centile(img) qui, étant donné une image en niveaux de gris, calcule et renvoie les 2ème et 98ème centiles.

np.percentile(A, x) calcule le xième centile de A.

Vous pouvez utiliser la fonction ci-dessus ou, exceptionnellement, une boucle for sur l'histogramme. Le résultat pour le raton laveur devrait être (11, 213) (ou des valeurs similaires).

ANSWER ELEMENTS

def centiles(img):
    h = histogram(img)
    p2 = None
    p92 = None

    s = 0
    for i in range(256):
        s += h[i]
        if s < 0.02:
            p2 = i
        if s < 0.98:
            p98 = i

    return p2, p98

centiles(gray)

Résoudre le système

Numpy fournit un module pour résoudre les systèmes d'équations linéaires. Les équations doivent être exprimées sous forme matricielle.

{$$ \left( \begin{matrix} p_2 & 1 \\ p_{98} & 1 \end{matrix} \right) * \left( \begin{matrix} C \\ I \end{matrix} \right) = \left( \begin{matrix} 0 \\ 255 \end{matrix} \right) $$}

Écrire une fonction find_corrections(img) qui renvoie les corrections de contraste et d'intensité à effectuer. Voici les instructions utiles :

from numpy import linalg
X = linalg.solve(A, B) avec AX = B

ANSWER ELEMENTS

from scipy import linalg

def find_corrections(img):
    p2, p98 = centiles(img)
    A = np.array([[p2, 1], [p98, 1]])
    B = np.array([0, 255])
    X = linalg.solve(A, B) # AX = B
    return X  # contrast, intensity

find_corrections(gray)

Contraste automatique

Ecrivez une fonction auto_contrast(img) qui prend une image en niveaux de gris comme paramètre et renvoie une nouvelle image avec l'application de l'algorithme de contraste automatique.

Voici quelques fonctions utiles :

np.zeros_like(A) crée un nouveau tableau de la même forme et de la même taille que A et contenant des zéros ;
np.clip(A, min, max) crée un nouveau tableau dans lequel les valeurs de A inférieures à min sont remplacées par min et celles supérieures à max sont remplacées par max.

ANSWER ELEMENTS

def auto_contraste(img):
    C, I = find_corrections(img)
    corrected = np.zeros_like(img)
    corrected[...] = img[...] * C + I
    return np.clip(corrected[...], 0, 255)

corrected = auto_contraste(gray)
plot_intensity(gray)
plot_intensity(corrected)

Bravo

Vous avez terminé ce TD. Vous savez maintenant comment fonctionnent les vues numpy et comment utiliser numpy pour éviter les boucles explicites. Si vous avez encore du temps, vous pouvez répondre aux questions suivantes.

Questions optionnelles - Image couleur

Histogramme des couleurs

Ecrivez une fonction plot_color_intensities(img) qui affiche une image et, à côté, les 3 histogrammes d'intensité (R, G et B) superposés. Nous pouvons utiliser la fonction histogram(img) que nous avons créée précédemment pour récupérer l'histogramme d'intensité d'une image en niveaux de gris et l'appliquer à une vue de chaque couleur.

Fonction utile :

plt.fill_between(x, y, color="#FF000055") pour afficher une courbe pleine rouge remplie (les codes couleurs notation hexadécimale).

ANSWER ELEMENTS

def plot_color_intensities(img):  
    bins_r = histogram(img[:,:,0])
    bins_g = histogram(img[:,:,1])
    bins_b = histogram(img[:,:,2])
    x = np.linspace(0, 255, 256)

    plt.subplots(figsize=(15, 5))
    plt.subplot(1, 2, 1)
    plt.title("Image")
    plt.imshow(img, norm=None)
    plt.xticks([]), plt.yticks([])

    plt.subplot(1, 2, 2)
    plt.title("Intensities")
    plt.fill_between(x, bins_r, color="#FF000055")
    plt.fill_between(x, bins_g, color="#00FF0055")
    plt.fill_between(x, bins_b, color="#0000FF55")
    plt.show()

plot_color_intensities(img)

Conversion du modèle de couleur

Notre image est actuellement codée avec le modèle de couleur RGB. Dans ce modèle, l'intensité est "distribuée" sur chacune des trois couleurs. On pourrait imaginer un algorithme d'auto-contraste sur les images couleur en travaillant sur chaque intensité indépendamment. Ce faisant, on prend le risque de modifier la teinte d'une image.

Pour l'auto-contraste, il est préférable de changer de modèle de couleur, par exemple le modèle YCbCr. Dans ce modèle, chaque pixel est représenté par trois composantes. La première Y est l'intensité (0 pour le noir, 255 pour le blanc). Les deuxième et troisième composantes représentent la variation de la teinte.

On peut passer du modèle RVB au modèle YCbCr :

{$$ \left( \begin{matrix} Y \\ C_b \\ C_r \end{matrix} \right) = \left( \begin{matrix} 0 \\ 128 \\ 128 \end{matrix} \right) + \left( \begin{matrix} 0.299 & 0.587 & 0.114 \\ −0.168736 & −0.331264 & 0.5 \\ 0.5 & −0.418688 & −0.081312 \end{matrix} \right) * \left( \begin{matrix} R \\ G \\ B \end{matrix} \right) $$}

et inversement

{$$ \left( \begin{matrix} R \\ G \\ B \end{matrix} \right) = \left( \begin{matrix} 1 & 0 & 1.402 \\ 1 & −0.344136 & −0.714136 \\ 1 & 1.772 & 0 \end{matrix} \right) * \left( \begin{matrix} Y \\ C_b-128 \\ C_r-128 \end{matrix} \right) $$}

Ecrivez deux fonctions to_YCbCr(img) et to_RGB(img) pour passer d'un modèle à l'autre. Ecrivez aussi une fonction test_conversion(img) qui vérifie que tout est en ordre (si c'est le cas, la différence entre chaque pixel d'une image img et to_RGB(to_YCbCr(img)) doit être faible (moins de 1 si on oublie que les valeurs sont entières, quelques unités dans le cas contraire).

Conseils : La matrice de conversion d'un pixel est de dimension (3, 3). Cependant, une image est de dimension (rows, cols, 3). La documentation numpy indique à propos de la fonction de produit np.dot(a, b) : "un produit de somme sur le dernier axe de a et l'avant-dernier de b". La solution la plus simple consiste donc à inverser l'ordre des opérandes du produit et donc à transposer la matrice de conversion.

Voici quelques fonctions utiles :

np.dot(A, B), produit ;
A.T est la transposée de A ;
A.astype(int) construit un tableau avec les valeurs arrondies de A ;
np.clip(A, 0, 255) construit un tableau avec les valeurs de l'intervalle [0,255] ;
np.absolute(A) construit un tableau avec les valeurs absolues de A.

ANSWER ELEMENTS

def to_YCbCr(img):
    v = np.array([0, 128, 128])
    rgb2ycbcr = np.array([[0.299, 0.587, 0.114], [-0.168736, -0.331264, 0.5], [0.5, -0.418688, -0.081312]])
    tmp = np.dot(img, rgb2ycbcr.T)+v
    tmp = tmp.astype(int)
    return np.clip(tmp, 0, 255)

def to_RGB(img):
    v = np.array([0, 128, 128])
    ycbcr2rgb = np.array([[1, 0, 1.402], [1, -0.344136, -0.714136], [1, 1.772, 0]])
    tmp = img-v
    tmp = np.dot(tmp, ycbcr2rgb.T)
    tmp = tmp.astype(int)
    return np.clip(tmp, 0, 255)


def test_convertion(img):
    ycc = to_YCbCr(img)
    rgb = to_RGB(ycc)
    diff = np.absolute(img - rgb)
    print("Min, Max and Mean, of differences between original image and convert-reconvert image : {}, {}, {}".format(
                np.min(diff), np.max(diff), np.mean(diff)))  

test_convertion(img)

Couleur auto constrast

Ecrire une fonction auto_contrast_YCbCr(img) qui convertit img en mode YCbCr, applique l'auto-contraste à la composante Y, puis renvoie l'image convertie en mode RGB.

ANSWER ELEMENTS

def auto_contrast_YCbCr(img):
    ycbcr = to_YCbCr(img)
    ycbcr[:,:,0] = auto_contraste(ycbcr[:,:,0])
    return to_RGB(ycbcr)

plot_color_intensities(img) 
plot_color_intensities(auto_contrast_YCbCr(img))