[image] Pyramide Gaussienne (scale-space)

Version imprimable

Voici une implémentation Java pour calculer une pyramide gaussienne d'images, utilisée notamment pour modéliser l'espace d'échelle (scale-space).

http://xphilipp.developpez.com/contr...scalespace.png

Une petite classe pour stocker un niveau de la pyramide

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
//Pyramid Level structure
public class PyramidLevel {
	// size of the data
	public int width,height;
	// the values
	public float[][] data;
	// the scale value (relative to original image)
	public double sigma2;
	// constructor allocating a new array
	public PyramidLevel(int width,int height,double sigma2) {
		this.width=width;
		this.height=height;
		this.sigma2=sigma2;
		this.data=new float[width][height];
	}
	// return the data at position(x,y)
	public float getData(int x, int y) {
		if (x<0) x=0; else if (x>=this.width)  x=this.width-1;
		if (y<0) y=0; else if (y>=this.height) y=this.height-1;
		return this.data[x][y];
	}
}

une autre pour stocker un noyau de convolution gaussien

Code:

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
public class GaussianKernel {
	// the kernel values (integer)
	public final float[] data;
	// the normalizer (sum of elements)
	public final float normalizer;
	// the radius of the kernel
	public final int radius;
 	// the sigma2 value
	public final double sigma2;
 
	// public constructor
	public GaussianKernel(double sigma2) {
		this.sigma2 = sigma2;
		// radius = 1*sigma, 2*sigma and 3*sigma 
		// represent respectivly 68%, 95% and 99% of the distribution 
		radius = (int)Math.round(2*Math.sqrt(sigma2));
		// compute gaussian values
		data = new float[1+2*radius];
		for(int r=-radius;r<=radius;r++)
			data[r+radius] = (float)Math.exp(-(r*r)/(2.0*sigma2));
		// compute the normalizer
		float sum=0;
		for(int i=0;i<data.length;i++)	sum+=data[i];
		normalizer=sum;
	}
}

Et enfin le code principal de construction de la pyramide

Code:

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
/**
 * Construction d'une pyramide gaussienne d'images
 *
 * @author Xavier PHILIPPEAU
 */
public class GaussianPyramid {
 
	// the levels of the pyramid 
	private PyramidLevel[][] levels = null;
 
	// number of octaves and scales
	private int octaves=0, scales=0;
 
	// public constructor
	public GaussianPyramid(int octaves, int scales) {
		this.octaves = octaves;
		this.scales = scales;
		this.levels = new PyramidLevel[this.octaves][this.scales+1];
	}
 
	// ---- private methodes ------------------------------------------------------------------
 
	// create a downsampled version of the level "src" 
	// L(x,y,s) = L(x/2,y/2,s/4)
	private PyramidLevel newDownSampled(PyramidLevel src) {
		int w = src.width/2, h = src.height/2;
		PyramidLevel dest = new PyramidLevel(w,h,src.sigma2/4);
		for(int y=0;y<h;y++)
			for(int x=0;x<w;x++)
				dest.data[x][y]=src.getData(2*x+1,2*y+1);
		return dest;
	}
 
	// create a downscaled version of the level "src" using the given kernel
	// L(x,y,s) = L(x,y,t) * G(x,y,w) 
	private PyramidLevel newDownScaled(PyramidLevel src, GaussianKernel kernel) {
		PyramidLevel dest = new PyramidLevel(src.width,src.height,src.sigma2+kernel.sigma2);
		convolve2D(src,dest,kernel);
		return dest;
	}
 
	// 2D convolution (using kernel separation) 
	private void convolve2D(PyramidLevel src, PyramidLevel dest, GaussianKernel kernel) {
		PyramidLevel buffer = new PyramidLevel(src.width,src.height,0.0);
		float normalizer=kernel.normalizer*kernel.normalizer;
		float v;
		// horizontal
		for(int y=0;y<src.height;y++) {
			for(int x=0;x<src.width;x++) {
				v=0;
				for(int k=-kernel.radius;k<=kernel.radius;k++)
					v+=kernel.data[kernel.radius+k]*src.getData(x+k,y);
				buffer.data[x][y]=v;
			}
		}
		// vertical
		for(int y=0;y<src.height;y++) {
			for(int x=0;x<src.width;x++) {
				v=0;
				for(int k=-kernel.radius;k<=kernel.radius;k++)
					v+=kernel.data[kernel.radius+k]*buffer.getData(x,y+k);
				dest.data[x][y]=v/normalizer;
			}
		}
	}
 
	// ---- public methodes -------------------------------------------------------------------
 
	// build the pyramid using the given image
	public void build(int[][] img, int width, int height) {
 
		// copy initial image into our custom structure
		this.levels[0][0] = new PyramidLevel(width,height,1.0);
		for(int y=0;y<height;y++) 
			for(int x=0;x<width;x++)
				this.levels[0][0].data[x][y]=img[x][y]; 
 
		// precompute the gaussian kernels used at each scale such that
		GaussianKernel[] kernels = new GaussianKernel[this.scales+1];
		// the scales increase geometricaly from 1.0 to 4.0 => K = 4^(1/scales)
		double K = Math.pow(4,1.0/this.scales);
		for(int s=1; s<this.scales+1; s++) {
			// Level[s] = Level[s-1]*G(w) => K^s = K^(s-1) + w
			double w = Math.pow(K, s) - Math.pow(K, s-1);
			kernels[s] = new GaussianKernel(w);
		}
 
		// for each octave
		for(int o=0; o<this.octaves; o++) {
 			// first scale : create a downsampled copy of the previous octave
			if (o>0) this.levels[o][0] = newDownSampled(this.levels[o-1][this.scales]);
			// subsequent scales : create a convolved copy of the previous scale 
			for(int s=1; s<this.scales+1; s++)
				this.levels[o][s] = newDownScaled(this.levels[o][s-1],kernels[s]);
		}
	}
 
	// return value of pixel(x,y) at scale(o,s) using bilinear interpolation
	public float getValue(int x, int y, int o, int s) {
		PyramidLevel p = this.levels[o][s];
		// coords of the pixel (x,y) in the (possibly downsampled) data 
		int xp=x>>o, yp=y>>o;
		// size of the "original" area of pixels (before downsampling)
		int size = 1<<o;
		// position of the pixel (x,y) in the "original" area
		int ypos=y-(yp<<o);
		int xpos=x-(xp<<o);
		// bilinear interpolation
		float v = 0;
		v+=(size-ypos)*( (size-xpos)*p.getData(xp,yp)   + (xpos)*p.getData(xp+1,yp)   );
		v+=(     ypos)*( (size-xpos)*p.getData(xp,yp+1) + (xpos)*p.getData(xp+1,yp+1) );
		return v/(size*size);
	}
 
	// ---- getters / setters -----------------------------------------------------------------
 
	public PyramidLevel getLevel(int o,int s) {
		return this.levels[o][s];
	}
 
	public double getSigma2(int o,int s) {
		return Math.pow(4, o)*this.levels[o][s].sigma2;
	}
 
	public int getOctaves() {
		return octaves;
	}
 
	public int getScales() {
		return scales;
	}
}

Voici un exemple d'utilistation de cette classe:

Code:

1
2
3
4
5
6
7
8
9
// l'image initiale en niveaux de gris
int[][] data = new int[width][height];
 
// construction d'une pyramide constituée de 4 octaves, avec 3 niveaux par octaves
GaussianPyramid gp = new GaussianPyramid(4, 3);
gp.build(data, width, height);
 
// récupération de la valeur d'un pixel a une position/echelle donnée
int pixel = gp.getValue(x, y, o, s);

Un exemple d'utilisation de la pyramide Gaussienne pour construire un détecteur de Harris-Laplace :

http://xphilipp.developpez.com/contr...risLaplace.png

Code:

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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
public class HarrisLaplace {
 
	private GaussianPyramid gp = null;
 
	// compute Gradient (sobel) 
	private float[] sobel(int x, int y, int o, int s, float[] sobel) {
		PyramidLevel p = this.gp.getLevel(o, s);
		float v00 = p.getData(x-1, y-1), v10 = p.getData(x, y-1), v20 = p.getData(x+1, y-1);
		float v01 = p.getData(x-1, y  ),                          v21 = p.getData(x+1, y  );
		float v02 = p.getData(x-1, y+1), v12 = p.getData(x, y+1), v22 = p.getData(x+1, y+1);
		float sx = (v20+2*v21+v22)-(v00+2*v01+v02);
		float sy = (v02+2*v12+v22)-(v00+2*v10+v20);
		if (sobel==null) sobel=new float[2];
		sobel[0]=sx/4; sobel[1]=sy/4;
		return sobel;
	}
 
	// compute Laplacian
	private float laplacian(int x,int y, int o, int s) {
		PyramidLevel p = this.gp.getLevel(o, s);
		float                            v10 = p.getData(x, y-1);
		float v01 = p.getData(x-1, y  ), v11 = p.getData(x, y  ), v21 = p.getData(x+1, y  );
		float                            v12 = p.getData(x, y+1);
		float laplacian = -v11 + 0.25f*(v10+v01+v21+v12);
		return laplacian;
	}
 
	// compute Scale-Normalized Laplacian
	private float laplacianNormalized(int x,int y, int o, int s) {
		float sigma2 = (float)this.gp.getLevel(o, s).sigma2;
		return sigma2*laplacian(x,y,o,s);
	}
 
	// compute Harris measure
	private float harris(int x,int y, int o, int s) {
		float m00=0, m01=0, m10=0, m11=0;
		float[] sobel = new float[2];
		// Gaussian kernel of the integration window
		float[][] kernel = {{1,2,1},{2,4,2},{1,2,1}};
		float normalizer=16;
		// for each pixel in the integration window
		for(int dy=-1;dy<=1;dy++) {
			for(int dx=-1;dx<=1;dx++) {
				int xk = x + dx;
				int yk = y + dy;
				// gradient value
				sobel(xk,yk,o,s,sobel);
				float gx = sobel[0];
				float gy = sobel[1];
				// gaussian weight
				float w = kernel[1+dx][1+dy];
				// second-moment matrix elements
				m00 += gx * gx * w;
				m01 += gx * gy * w;
				m10 = m01;
				m11 += gy * gy * w;
			}
		}
		// harris measure
		float harris = m00*m11 - m01*m10 - 0.06f*(m00+m11)*(m00+m11);
		harris = harris/(255*255*normalizer);
		return harris;
	}
 
	// return true if the pixel (x,y,o,s) is a scale extrema 
	private boolean isScaleExtrema(int x,int y, int o, int s) {
		float lap = laplacianNormalized(x,y,o,s);
		if (Math.abs(lap)<1) return false;
 
		float lap_prev = laplacianNormalized(x,y,o,s-1);
		float lap_next = laplacianNormalized(x,y,o,s+1);
		if (lap_prev<lap && lap>lap_next) return true;
		if (lap_prev>lap && lap<lap_next) return true;
		return false;
	}
 
	// return true if the pixel (x,y,o,s) is a spatial maxima 
	private boolean isSpatialMaxima(int x,int y, int o, int s) {
		double h = harris(x,y,o,s);
		if (Math.abs(h)<1) return false;
 
		if (h<=harris(x-1,y-1,o,s) || h<=harris(x  ,y-1,o,s) || h<=harris(x+1,y-1,o,s)) return false;
		if (h<=harris(x-1,y  ,o,s) ||                           h<=harris(x+1,y  ,o,s)) return false;
		if (h<=harris(x-1,y+1,o,s) || h<=harris(x  ,y+1,o,s) || h<=harris(x+1,y+1,o,s)) return false;
		return true;			
	}
 
	// main algorithm
	public void detect(int[][] image, int width, int height) {
		int octaves=5, scales=9;
 
		this.gp = new GaussianPyramid(octaves, scales);
		this.gp.build(image, width, height);
 
		// for each octave+scale
		for(int o=0;o<gp.getOctaves();o++) {
			for(int s=1;s<gp.getScales()-1;s++) {
				int w = gp.getLevel(o, 0).width;
				int h = gp.getLevel(o, 0).height;
 
				// explore all pixels of the level
				for (int y=1; y<h-1; y++) {
					for (int x=1; x<w-1; x++) {
 
						// keep only spatial maxima
						if (!isSpatialMaxima(x, y, o, s)) continue;
 
						// keep only scale extrema
						if (!isScaleExtrema(x, y, o, s)) continue;
 
						// localization in original image (average)
						int xc = (int)Math.round((0.5+x)*(1<<o)); 
						int yc = (int)Math.round((0.5+y)*(1<<o)); 
						int sc = (int)Math.round(gp.getSigma2(o,s)); 
 
						System.out.printf("Found points x=%d, y=%d, sigma²=%d\n",xc,yc,sc);
					}
				}
			}
		}
	}
}

01/07/2009, 15h38
godhiva

Bonjour,

J'ai une question relative aux images d'illustration du détecteur de Harris-Laplace.

Je suppose que les centre des cercles rouges correspondent aux points d'intérêts détectés et que les rayon des cercles ont rapport avec le niveau de la pyramide ou les points ont été trouvés. Mais comment les as tu déterminé ? Est-ce arbitraire ?
01/07/2009, 17h40
pseudocode

Citation:

Envoyé par godhiva

Je suppose que les centre des cercles rouges correspondent aux points d'intérêts détectés et que les rayon des cercles ont rapport avec le niveau de la pyramide ou les points ont été trouvés. Mais comment les as tu déterminé ? Est-ce arbitraire ?

J'ai utilisé la formule : rayon = 2*racine(1+s), ou "s" et l'échelle du point d'intérêt.

Cela s'explique en notant que j'ai calculé les valeurs du Laplacien et de Harris avec un petit noyau 3x3 sigma²=1 sur l'image réduite. C'est équivalent à utiliser un grand noyau sigma²=(1+s) sur l'image de départ, ce qui correspondrait a prendre un noyau gaussien de rayon 2*racine(sigma²) dans la classe GaussianKernel.