[Image] Gradient et Hessienne par Convolution

[pseudocode]

Bonjour,

Parmis toutes les objets mathématiques qu'il est souvent necessaire de calculer sur les pixels d'une image, il y en a deux qui reviennent assez souvent:

- Le vecteur Gradient, c'est a dire les derivées partielles du 1er ordre de l'intensité:

Code :

1
2
3
4
           dI(x,y)    dI(x,y)
G(x,y) = ( ------- , -------- )
             dx         dy

- La matrice Hessienne, c'est a dire les derivées partielles du 2nd ordre de l'intensité:

Code :

1
2
3
4
5
6
7
8
         | d²I(x,y)    d²I(x,y) |
         | --------    -------- |
         |    dx²        dxdy   |
H(x,y) = |                      |
         | d²I(x,y)    d²I(x,y) |
         | --------    -------- |
         |   dxdy         dy²   |

Remarque: la trace de cette matrice est appelé "Laplacien"


De part la nature discrete d'une image numérisée, les derivées sont approximées par des différences d'intensité au voisinage du pixel. Ces calculs vont dont etre tres sensibles au bruit. Pour limiter les effets du bruit, il est d'usage de lisser préalablement l'image (en la convoluant par une Gaussienne) avant de calculer les derivées.

Et c'est la que la magie de la convolution intervient car:

Code :

1
2
3
4
      d                        dF(t) 
     ---- [ F(t) * I(t) ] =  [ ----- ] * I(t)
      dt                        dt

Donc on peut inverser le calcul: d'abord calculer les derivées de la Gaussienne, puis convoluer ce résultat avec l'image. Cela nous donne deux avantages:

- Réduire les calculs de derivées partielles de l'intensité à des convolutions.
- Précalculer les noyaux de convolution, car ils sont indépendants de l'image.


Calculons les dérivées partielles de la Gaussienne:

Code :

1
2
G(x,y) = Exp( - (x²+y²) / 2*sigma² )

Derivées partielles du 1er ordre

Code :

1
2
3
4
5
6
7
8
dG(x,y)   
------- = - (x/sigma²) * Exp( - (x²+y²) / 2*sigma² )
  dx     

dG(x,y)   
------- = - (y/sigma²) * Exp( - (x²+y²) / 2*sigma² )
  dy

Derivées partielles du 2nd ordre

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
d²G(x,y)   
-------- = (x²-sigma²)/(sigma^4) * Exp( - (x²+y²) / 2*sigma² )
  dx²     


d²G(x,y)   
-------- = (y²-sigma²)/(sigma^4) * Exp( - (x²+y²) / 2*sigma² )
  dy²     


d²G(x,y)   
-------- = (x*y)/(sigma^4) * Exp( - (x²+y²) / 2*sigma² )
  dxdy

Quelques remarques sur les noyaux de convolutions obtenus:

1. Modifiez la valeur centrale du noyau afin que la somme des coefficients de la matrice soit égale a zero. Cela permet de s'assurer que la dérivée partielle sera nulle sur une zone d'intensité uniforme.

2. Choisissez une taille minimum de 3x3 pour le Gradient et 5x5 pour la Hessienne. En effet, plus la dérivée est d'un ordre élevé, plus il faut de points pour la calculer.

3. La valeur de sigma² permet de moduler la force du lissage (generalement une valeur entre 1 et 3 pour le Gradient, et entre 2 et 5 pour la Hessienne)


Exemple de code en Java:

Code Java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
 
public class HessianMatrix {
 
	int halfwindow = 3; // 7x7 kernel
	double sigma2 = 2.0;
 
	double[][] kernelGxx = new double[2*halfwindow+1][2*halfwindow+1];
	double[][] kernelGyy = new double[2*halfwindow+1][2*halfwindow+1];
	double[][] kernelGxy = new double[2*halfwindow+1][2*halfwindow+1];
 
	// Constructor
	public HessianMatrix() {
		for(int y=-halfwindow;y<=halfwindow;y++) {
			for(int x=-halfwindow;x<=halfwindow;x++) {
				kernelGxx[halfwindow+y][halfwindow+x] = Gxx(x, y);
				kernelGyy[halfwindow+y][halfwindow+x] = Gyy(x, y);
				kernelGxy[halfwindow+y][halfwindow+x] = Gxy(x, y);
			}
		}
	}
 
	// Kernel functions (gaussian 2nd order partial derivatives)
	private double Gxx(int x, int y) {
		double t = (x*x+y*y)/(2*sigma2);
		double d2t = (x*x-sigma2) / (sigma2*sigma2);
		double e = d2t * Math.exp( -t );
		return e;
	}
	private double Gyy(int x, int y) {
		double t = (x*x+y*y)/(2*sigma2);
		double d2t = (y*y-sigma2) / (sigma2*sigma2);
		double e = d2t * Math.exp( -t );
		return e;
	}
	private double Gxy(int x, int y) {
		double t = (x*x+y*y)/(2*sigma2);
		double d2t = (x*y) / (sigma2*sigma2);
		double e = d2t * Math.exp( -t );
		return e;
	}
 
	// return the 2x2 Hessian Matrix for pixel(x,y) 
	public double[][] getMatrix(Channel c, int x, int y) {
		double fxx=0, fyy=0, fxy=0;
		for(int dy=-halfwindow;dy<=halfwindow;dy++) {
			for(int dx=-halfwindow;dx<=halfwindow;dx++) {
				int xk = x + dx;
				int yk = y + dy;
				double vk = c.getValue(xk,yk); // <-- value of the pixel
				fxx += kernelGxx[halfwindow-dy][halfwindow-dx] * vk;
				fyy += kernelGyy[halfwindow-dy][halfwindow-dx] * vk;
				fxy += kernelGxy[halfwindow-dy][halfwindow-dx] * vk;
			}
		}
 
		double[][] hessianMatrix = new double[][] {
				new double[] { fxx, fxy },
				new double[] { fxy, fyy },
		};
 
		return hessianMatrix;
	}	
}

Code Java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
 
public class GradientVector {
 
	int halfwindow = 1; // 3x3 kernel
	double sigma2 = 1.4;
 
	double[][] kernelGx = new double[2*halfwindow+1][2*halfwindow+1];
	double[][] kernelGy = new double[2*halfwindow+1][2*halfwindow+1];
 
	// Constructor
	public GradientVector() {
		for(int y=-halfwindow;y<=halfwindow;y++) {
			for(int x=-halfwindow;x<=halfwindow;x++) {
				kernelGx[halfwindow+y][halfwindow+x] = Gx(x, y);
				kernelGy[halfwindow+y][halfwindow+x] = Gy(x, y);
			}
		}
	}
 
	// Kernel functions (gaussian 1st order partial derivatives)
	private double Gx(int x, int y) {
		double t = (x*x+y*y)/(2*sigma2);
		double d2t = -x / sigma2;
		double e = d2t * Math.exp( -t );
		return e;
	}
	private double Gy(int x, int y) {
		double t = (x*x+y*y)/(2*sigma2);
		double d2t = -y / sigma2;
		double e = d2t * Math.exp( -t );
		return e;
	}
 
	// return the Gradient Vector for pixel(x,y) 
	public double[] getVector(Channel c, int x, int y) {
		double gx=0, gy=0;
		for(int dy=-halfwindow;dy<=halfwindow;dy++) {
			for(int dx=-halfwindow;dx<=halfwindow;dx++) {
				int xk = x + dx;
				int yk = y + dy;
				double vk = c.getValue(xk,yk); // <-- value of the pixel
				gx += kernelGx[halfwindow-dy][halfwindow-dx] * vk;
				gy += kernelGy[halfwindow-dy][halfwindow-dx] * vk;
			}
		}
 
		double[] gradientVector = new double[] { gx, gy };
 
		return gradientVector;
	}	
}

[pseudocode]

A gauche: la norme du vecteur Gradient, Noyau 3x3, sigma²=1.4 (echelle Logarithmique)

A droite: la trace de la matrice Hessienne (Laplacien), Noyau 5x5, sigma²=2.0 (effet emboss)

[Methode]

Bonjour,

Par rapport à ces codes peux-tu m'expliquer d'où provient la classe channel ?
Par la même occasion, connais-tu les histogrammes pondérés et si ces codes pourraient s'y appliquer ?

Merci.

[pseudocode]

Envoyé par Methode

Par rapport à ces codes peux-tu m'expliquer d'où provient la classe channel ?

Ah. La classe Channel est une classe de ma librairie personnelle qui représente une image 8-bit, c'est à dire un tableau 2D d'entiers. Les accesseurs sont chargés de gérer les débordements (border-extender).

Une implémentation minimaliste possible est:

Code java :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
 
public class Channel {
 
	private int width,height;
	private int[][] values = null;
 
	public Channel(int width,int height,int[][] array) {
		this.width = width;
		this.height = height;
		this.values = array;
	}
 
	public int getValue(int x,int y) {
		// gestion des débordements
		if (x<0) x=0;
		if (x>=this.width) x=this.width-1;
		if (y<0) y=0;
		if (y>=this.height) y=this.height-1;
 
		return this.values[x][y];
	}
 
	public void setValue(int x,int y,int v) {
		// gestion des débordements
		if (x<0) return;
		if (x>=width) return;
		if (y<0) return;
		if (y>=height) return;
 
		this.values[x][y]=v;
	}
 
	public int getWidth() {
		return this.width;
	}
 
	public int getHeight() {
		return this.height;
	}
 
	public int[][] asArray() {
		return this.values;
	}
}

[albert1983]

Bonjour,

Si j'ai bien compris.. dans un contexte d'une image :
G(x,y) represente la valeur du pixel au point (x,y).

Le gradient de chaque pixel nous donne: grad(G(x,y)) = (dG/dx,dG/dy)

Donc, pour chaque valeur de pixel j'aurais un couple de valeur!! Ce que je ne comprends pas, comment reproduire l'imagedu "gradient de l'image" avec ces nouveaux valeurs..

Merci d'avance

[pseudocode]

Envoyé par albert1983

Donc, pour chaque valeur de pixel j'aurais un couple de valeur!! Ce que je ne comprends pas, comment reproduire l'image du "gradient de l'image" avec ces nouveaux valeurs..

Ce qu'on représente généralement ce n'est pas l'image du "gradient", mais l'image de la "norme du gradient".

Donc à partir d'un couple de valeur (a,b) on calcule la valeur n=racine(a²+b²). C'est cette valeur "n" qui est utilisée dans l'image.

[albert1983]

Superrr

Merci beaucoup...

Je passe mnt au codage..

[albert1983]

Hello..

Je me suis basé sur le code disponible dans ce post pour programmer une fonction qui évalue le gradient d'une image.

Mon but est de trouver le point de recalage avec une grande précision en utilisant les transformées de Fourier des images.

Lorsque j'applique le gradient sur les deux images, j'ai un phénomène bizarroïde:

quand j'évalue la norme du gradient comme vous avez suggéré n=racine(a²+b²). ça donne de mauvaise résultats.

lorsque je prend seulement le gradient vertical n = a; j'obtiens le résultat voulu dans un axe, et résultat entier dans l'autre.
ex: pour (dx = 13.89 et dy = 4.33) j'obtiens (dx = 14 et dy = 4.35)
pour (dx = 15.39 et dy = 27.7) j'obtiens (dx = 15.40 et dy = 28)
etc..



Quelq'un comprend ce qui se passe??

Code C :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
double ApplyGradient(double* pixels, int pos_x, int pos_y, int largeur, int hauteur) 
	{
		double getValue(int, int, double*, int, int); 
		double Gx(int, int);
		double Gy(int, int);
 
		double gx=0, gy=0, gradient=0;
		int halfwindow = 1; // 3x3 kernel
 
		double kernelGx[2*halfwindow+1][2*halfwindow+1];
		double kernelGy[2*halfwindow+1][2*halfwindow+1];
 
		for(int y=-halfwindow;y<=halfwindow;y++) 
		{
			for(int x=-halfwindow;x<=halfwindow;x++) 
			{
				kernelGx[halfwindow+y][halfwindow+x] = Gx(x, y);
				kernelGy[halfwindow+y][halfwindow+x] = Gy(x, y);
			}
		}
 
		for(int dy=-halfwindow;dy<=halfwindow;dy++) 
		{
			for(int dx=-halfwindow;dx<=halfwindow;dx++) 
			{
				int xk = pos_x + dx;
				int yk = pos_y + dy;
				double vk = getValue(xk,yk, pixels, largeur, hauteur); // <-- value of the pixel
				gx += kernelGx[halfwindow-dy][halfwindow-dx] * vk;
				gy += kernelGy[halfwindow-dy][halfwindow-dx] * vk;
			}
		}
		gradient = gx;
		return gradient;
	}

[pseudocode]

Si c'est un problème de précision des valeurs retournées, tu peux prendre un noyau plus grand (halfwindow=2, 3, ...).

Par contre cela va atténuer les hautes frequences de l'image, et donc tu risques de perdre en précision sur les petites variations.

[albert1983]

Je ne pense pas que c'est un probleme de precision.
Car la precision est tres bonne (0.01 a 0.05)
Cependant, lorsque la precision est bonne dans un axe, elle ne l'est pas dans l'autre qui affiche bizarrement un "entier" au lieu d'une valeur "double".

Voici un exemple des resultats :

dxTheo dyTheo dxExp dyExp Err dx Err dy
(0.33, 3.4) (0.00,3.45) Erreur : [0.33 0.05]
(1.08 3.29) (1) 3.31) Erreur [0.08 0.02]
(2 5.85) (2) 5.79) Erreur [0 0.06]
(18.38 3.12) (18.39 3) Erreur [0.01 0.12]
(22.25 16.37) (22 16.42) Erreur [0.25 0.05]
(13.89 4.33) (14 4.35) Erreur [0.11 0.02]
(36.2 47.31) (36 47.33) Erreur [0.2 0.02]
(15.39 27.7) (15.4 28 Erreur [0.01 0.3]
(16.27 15.28) (16.37 15) Erreur [0.1 0.28]
(0.31 4) (0.31 4) Erreur [0 0]

[pseudocode]

Envoyé par albert1983

Je ne pense pas que c'est un probleme de precision.
Car la precision est tres bonne (0.01 a 0.05)
Cependant, lorsque la precision est bonne dans un axe, elle ne l'est pas dans l'autre qui affiche bizarrement un "entier" ou une valeur "double".

Tu peux donner les valeurs de ton image (les voisinages 3x3) utilisés pour tes calculs ?

[albert1983]

Voici un exemple des valeurs de l'image :

22813.000000 22813.000000 22812.999997
22813.000000 22813.000000 22812.999997
22812.141112 22812.141112 22812.137753

[pseudocode]

Voila ce que je trouve:

Code :

1
2
3
4
5
6
7
8
9
10
11
12
Value (*) KernelGx:

22813,000000 22813,000000 22812,999997     0,367879 0,000000 -0,367879
22813,000000 22813,000000 22812,999997 (*) 0,606531 0,000000 -0,606531 = 0.0012386302769300528
22812,141112 22812,141112 22812,137753     0,367879 0,000000 -0,367879

Value (*) KernelGy:

22813,000000 22813,000000 22812,999997      0,367879441  0,606530660  0,367879441
22813,000000 22813,000000 22812,999997 (*)  0,000000000  0,000000000  0,000000000 = 1.1541109835980023
22812,141112 22812,141112 22812,137753     -0,367879441 -0,606530660 -0,367879441

Ce qui nous donne donc:

Code :

1
2
3
4
Gradient = (gx,gy) = (0.0012386302769300528,1.1541109835980023)
Norme    = racine(gx²+gy²) = 1.1541116482674072
Angle    = arctan(gy/gx) = 89.93850828819237°

[albert1983]

Moi j'obtiens :

Pour Gx:
22812.137753 22812.141112 22812.141112 -0.349673 0.000000 0.349673
22812.999997 22813.000000 22813.000000 * -0.499766 0.000000 0.499766
22812.999997 22813.000000 22813.000000 -0.349673 0.000000 0.349673

Pour Gy:
22812.137753 22812.141112 22812.141112 -0.349673 -0.499766 -0.349673
22812.999997 22813.000000 22813.000000 * 0.000000 0.000000 0.000000
22812.999997 22813.000000 22813.000000 0.349673 0.000000 0.349673

gx : 0.001177
gy : 1.031076

norme : 1.031077
angle : 89.934590

[pseudocode]

Il y a sans doute une erreur de copier/coller dans ton noyau Gy (il n'est pas symetrique). Sinon, on obtient le meme genre de résultat aux arrondis et à la valeur de sigma2 près.

Donc je ne vois pas dans quel cas tu obtiens "toujours" un nombre entier.

[albert1983]

Je pense que j'ai trouve ou est le probleme..
En regardant ma matrice de correlation je vois un pic qui apparait au point (0,0) ce qui me fausse mes calculs..
Je pense que c a cause des conditions aux limites de l'image.. vu que j'utilise un tableau d'une seule dimension tab[largeur*hauteur] au lieu d'une matrice tab[largeur][hauteur]..
je reverifie les conditions au bords..

[albert1983]

J'ai verifie le code.. et j'ai mis les conditions suivantes pour la gestion des bords :
(j'ai un tableau unidimensionnel qui represente l'image)

if (x<0) x=largeur-1;
if (x>=largeur) x=largeur-1;
if (y<0) y=hauteur-1;
if (y>=hauteur) y=0;

return pixels[y*largeur+x];


quand j'ai les trois pics correctes c merveilleux mais parfois j'ai tjrs pic au point (0,0)..
Pensez vous que les conditions sont erronees?

[pseudocode]

Envoyé par albert1983

En regardant ma matrice de correlation je vois un pic qui apparait au point (0,0) ce qui me fausse mes calculs..

Correlation ? Quel rapport avec le calcul du Gradient/Hessienne ?

Envoyé par albert1983

J'ai verifie le code.. et j'ai mis les conditions suivantes pour la gestion des bords

Code C :

1
2
3
4
5
6
7
 
if (x<0) x=largeur+x;
if (x>=largeur) x=x-largeur;
if (y<0) y=hauteur+y;
if (y>=hauteur) y=y-hauteur;
 
return pixels[y*largeur+x];

[albert1983]

[QUOTE=pseudocode;3718427]Correlation ? Quel rapport avec le calcul du Gradient/Hessienne ?

Merci pseudo code.. Je mesure de shift entre deux images en utilisant la methode de correlation de phase. Pour des meilleurs resultats, je dois travailler sur le gradient des images au lieu de travailler sur les images elles meme..

[pseudocode]

Envoyé par albert1983

Merci pseudo code.. Je mesure de shift entre deux images en utilisant la methode de correlation de phase. Pour des meilleurs resultats, je dois travailler sur le gradient des images au lieu de travailler sur les images elles meme..

Ah... Ton problème n'a donc plus rien à voir avec le calcul du gradient ?