TP2 - Recherche par similarité par transformée de Gabor.

Objectif : comprendre comment fonctionnent les opérateurs de convolution avec application à la transformée de Gabor.

Exemple dérivées partielles

Utiliser jupyter.

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import numpy as np
import math
from scipy import signal
import time

Filtre de Gabor

Cette fonction retourne un tableau de nombres complexes correspondant à un filtre de Gabor défini selon plusieurs paramètres :

def gabor_kernel(lambd = 16.0, theta = 0.0, n = 0, sl = 0.7, st = 1.4, nl = 4.0):
	if n <= 0: n = 1+2*int(nl*lambd)
	gl = -0.5/(sl*sl)
	gt = -0.5/(st*st)
	c = math.cos(theta)/lambd
	s = math.sin(theta)/lambd
	x0 = 0.5*(n-1)*(c+s)
	y0 = 0.5*(n-1)*(c-s)
	sc = 1.0/(2*math.pi*sl*st*lambd*lambd)
	gk = np.empty( (n,n), dtype='complex64' )
	for y in range (0,n):
		for x in range (0,n):
			xr = c*x+s*y-x0  # centering, rotation and scaling
			yr = c*y-s*x-y0  # centering, rotation and scaling
			a = 2.0*math.pi*xr  # wave phase
			gk[y,x] = sc*math.exp(gl*xr*xr+gt*yr*yr)*complex(math.cos(a),math.sin(a))
	return gk

Paramètres :

lambd et theta : longueur d'onde et orientation de l'onde plane du filtre
n : si n est non nul, il détermine la taille du filtre (n x n) (sélection manuelle de la taille)
sl : facteur d'échelle de l'écart-type dans la direction de l'onde
st : facteur d'échelle de l'écart-type dans la direction perpendiculaire
nl : si n est nul, détermine la taille du filtre comme (2*int(nl*lambd)+1) (ajustement automatique de la taille)

Note : le code est optimisé pour minimiser le nombre d'opérations dans la boucle principale.
Note 2 : “lambd” parce que “lambda” est un mot-clé en python.
Note 3 : les valeurs par défaut vont bien pour toutes les valeurs autres que lambd et theta qui sont les deux paramètres utiles. On fera cependant varier les autres pour voir les différents effets si on les change.

Variante : onde seule (sans gaussienne)

def wave_kernel(lambd = 16.0, theta = 0.0, n = 0, sl = 0.7, st = 1.4, nl = 4.0):
	if n <= 0: n = 1+2*int(nl*lambd)
	gl = -0.5/(sl*sl)
	gt = -0.5/(st*st)
	c = math.cos(theta)/lambd
	s = math.sin(theta)/lambd
	x0 = 0.5*(n-1)*(c+s)
	y0 = 0.5*(n-1)*(c-s)
	sc = 1.0/(2*math.pi*sl*st*lambd*lambd)
	gk = np.empty( (n,n), dtype='complex64' )
	for y in range (0,n):
		for x in range (0,n):
			xr = c*x+s*y-x0  # centering, rotation and scaling
			yr = c*y-s*x-y0  # centering, rotation and scaling
			a = 2.0*math.pi*xr  # wave phase
			gk[y,x] = sc*complex(math.cos(a),math.sin(a))
	return gk

Variante : gaussienne seule (sans onde)

def gaussian_kernel(lambd = 16.0, theta = 0.0, n = 0, sl = 0.7, st = 1.4, nl = 4.0):
	if n <= 0: n = 1+2*int(nl*lambd)
	gl = -0.5/(sl*sl)
	gt = -0.5/(st*st)
	c = math.cos(theta)/lambd
	s = math.sin(theta)/lambd
	x0 = 0.5*(n-1)*(c+s)
	y0 = 0.5*(n-1)*(c-s)
	sc = 1.0/(2*math.pi*sl*st*lambd*lambd)
	gk = np.empty( (n,n), dtype='complex64' )
	for y in range (0,n):
		for x in range (0,n):
			xr = c*x+s*y-x0  # centering, rotation and scaling
			yr = c*y-s*x-y0  # centering, rotation and scaling
			a = 2.0*math.pi*xr  # wave phase
			gk[y,x] = sc*math.exp(gl*xr*xr+gt*yr*yr)
	return gk

Affichage d'un noyau complexe

def kernel_plot(k):
	kr = (k.view(np.float32).reshape(k.shape + (2,)))[:,:,0]  # extract real (cos) part
	ki = (k.view(np.float32).reshape(k.shape + (2,)))[:,:,1]  # extract imaginary (sin) part
	mpimg.imsave('kr.jpg',kr,cmap="gray")
	mpimg.imsave('ki.jpg',ki,cmap="gray")
	fig, (re, im) = plt.subplots(1, 2)  # real and imaginary parts
	re.imshow(kr, cmap='gray')
	re.set_title('Real part')
	re.set_axis_off()
	im.imshow(ki, cmap='gray')
	im.set_title('Imaginary part')
	im.set_axis_off()

1. Visualiser plusieurs variantes de filtre et comprendre l'effet des paramètres.

Onde plane avec longueur d'onde de 24 pixels et orientation de pi/3 (60 degrés, l'origine est en haut à gauche) :

gk = wave_kernel(lambd=24, theta = math.pi/3)
kernel_plot(gk)

Gaussienne elliptique (par défaut avec sl = 0.7 et st = 1.4)

gk = gaussian_kernel(lambd=24, theta = math.pi/3)
kernel_plot(gk)

Filtre de Gabor (produit des deux)

gk = gabor_kernel(lambd=24, theta = math.pi/3)
kernel_plot(gk)
gk.shape

Longueur d'onde de 12 pixels

gk = gabor_kernel(lambd=12, theta = math.pi/3)
kernel_plot(gk)
gk.shape

Les mêmes avec une taille de fenêtre fixe et à la même échelle :

gk = gabor_kernel(lambd=24, theta = math.pi/3, n = 256)
kernel_plot(gk)
gk.shape

gk = gabor_kernel(lambd=12, theta = math.pi/3, n = 256)
kernel_plot(gk)
gk.shape

Filtre circulaire : st = sl = 1

gk = gabor_kernel(lambd=24, theta = math.pi/3, sl = 1.0, st = 1.0)
kernel_plot(gk)
gk.shape

Filtre avec “effet de bord” (dans le filtre) :

gk = gabor_kernel(lambd=24, theta = math.pi/3, nl = 2.5)
kernel_plot(gk)
gk.shape

Faire varier l'orientation, l'échelle et l'ellipticité.

2. Application d'un filtre à une image

Lecture d'une image :

img = mpimg.imread('houses.jpg')
plt.imshow(img, cmap="gray")

Utilisation de la fonction convolve2d de scipy.signal et affichage du module de l'image transformée :

gk = gabor_kernel(lambd=4, theta = math.pi/2)
gab = signal.convolve2d(img, gk, boundary='symm', mode='valid')
plt.imshow(np.absolute(gab), cmap="gray")

Voir l'effet des changements d'orientation et d'échelle sur les images transformée.

3. Calcul de l'énergie moyenne à travers un filtre d'orientation et d'échelle donnés

Utiliser les fonctions numpy qui vont bien :

np.absolute() pour calculer le module d'un nombre complexe (élément par élément dans un tableau)
np.square() pour calculer le carré d'un nombre (élément par élément dans un tableau)
np.average() pour calculer la moyenne des valeurs dans un tableau

4. Calcul d'une transformée de Gabor

On considère 8 orientations (multiples de pi/8, partant de 0) et 4 longeurs d'ondes selon une échelle logarithmique partant de 3 et avec un facteur 2 (i.e.: 3, 6, 12 et 24 pixels). Il faut produire un tableau gt de 4 lignes par 8 colonnes dont les valeurs sont les énergies moyennes pour les orientations et les échelles correspondantes.

Fonction pour afficher le résultat :

plt.imshow(gt,cmap='gray',vmin=0,vmax=np.amax(gt))

5. Réhaussement (relatif) des courtes longeurs d'onde

On observe qu'il y a beaucoup d'énergie dans les grandes longueurs d'onde et peu dans les courtes. Cela est dû à la distribution de fréquence dans les images “naturelles”. On le compense en réhaussant les valeurs pour les petites longueurs d'onde (ici, on affaiblira plutôt les grandes, ce qui revient au même) en divisant les valeurs par un facteur proportionnel à la longueur d'onde (c'est à dire en divisant par 2 puissance n pour la nième longueur d'onde).

6. Calcul des transformées avec réhaussement pour les quatre image test et calcul de similarités par distance euclidienne.

Puis, pour chaque image, tri des autres par similarité.

Résultat attendu après réhaussement pour les quatre images :

7. Accélérer le calcul.

Un filtre de grande longueur d'onde appliqué à une grande image a approximativement le même effet qu'un filtre de taille réduite appliqué à une image réduite dans le même rapport. Remplacer le calcul sur les filtres de taille variable par un des calculs avec un seul filtre de taille fixe (longueur d'onde de trois pixels) sur des images dont la taille est divisée par 1, 2, 4 et 8. Effectuer le calcul sur deux images selon les deux méthodes, visualiser les transformées de Gabor et calculer les distances selon les deux méthodes. Comparer aussi la rapidité.

Vous pouvez utiliser la fonction resize() de PILLOW avec interpolation bilinéaire. Il faut convertir en image PILLOW.

from PIL import Image
imgp = Image.fromarray(img)
imgp = imgp.resize((img.shape[0]//2, img.shape[1]//2), resample=Image.BILINEAR)

8. Accélérer encore le calcul.

En théorie, il faut prendre une fenêtre large autour du filtre pour éviter les effets de fenêtrage (qui introduisent des fréquences artificielles et parasites). En principe, il faudrait prendre nl au mons égal à 6 dans la fonction gabor_kernel(). Nous avons pris une valeur par défaut égale à quatre. On peut encore réduire un peu cette valeur sans effet notable. Visualiser les filtres avec nl = 3 et nl = 2 et relancer les calculs avec ces valeurs. Comparer la rapidité et visualiser la différence sur les transformées de Gabor et sur les distances entre elles.

9. Implémentation manuelle d'un opérateur de pooling.

Plutôt qu'utiliser la fonction resize() de PILLOW, vous pouvez, à titre d'entraînement, implémenter un opérateur de “average pooling” 2×2 pour effectuer la réduction de taille de l'image.

10. Notes

Même si cette partie n'est pas notée, vous pouvez vous entraîner à noter les résultats obtenus et à présenter les résultats des différentes expériences dans un document (temps de calcul, effet sur les distances, images de transformées …). Il est plus simple de le faire au fur et à mesure.