Mercurial > pylearn

annotate pylearn/preprocessing/fft_whiten.py @ 1484:83d3c9ee6d65

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
* changed MNIST dataset to use config.get_filepath_in_roots mechanism
author gdesjardins
date Tue, 05 Jul 2011 11:01:51 -0400
parents 28b2f17991aa
children
rev   line source
953
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
1 """
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
2 Whitening algorithm used in Olshausen & Field, 1997.
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
3
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
4 If you want to use some other images, there are a number of
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
5 preprocessing steps you need to consider beforehand. First, you
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
6 should make sure all images have approximately the same overall
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
7 contrast. One way of doing this is to normalize each image so that
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
8 the variance of the pixels is the same (i.e., 1). Then you will need
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
9 to prewhiten the images. For a full explanation of whitening see
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
10
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
11 Olshausen BA, Field DJ (1997) Sparse Coding with an Overcomplete
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
12 Basis Set: A Strategy Employed by V1? Vision Research, 37: 3311-3325.
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
13
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
14 Basically, to whiten an image of size NxN, you multiply by the filter
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
15 f*exp(-(f/f_0)^4) in the frequency domain, where f_0=0.4*N (f_0 is the
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
16 cutoff frequency of a lowpass filter that is combined with the
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
17 whitening filter). Once you have preprocessed a number of images this
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
18 way, all the same size, then you should combine them into one big N^2
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
19 x M array, where M is the number of images. Then rescale this array
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
20 so that the average image variance is 0.1 (the parameters in sparsenet
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
21 are set assuming that the data variance is 0.1). Name this array
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
22 IMAGES, save it to a file for future use, and you should be off and
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
23 running. The following Matlab commands should do this:
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
24
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
25
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
26 N=image_size;
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
27 M=num_images;
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
28
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
29 [fx fy]=meshgrid(-N/2:N/2-1,-N/2:N/2-1);
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
30 rho=sqrt(fx.*fx+fy.*fy);
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
31 f_0=0.4*N;
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
32 filt=rho.*exp(-(rho/f_0).^4);
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
33
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
34 for i=1:M
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
35 image=get_image; % you will need to provide get_image
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
36 If=fft2(image);
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
37 imagew=real(ifft2(If.*fftshift(filt)));
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
38 IMAGES(:,i)=reshape(imagew,N^2,1);
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
39 end
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
40
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
41 IMAGES=sqrt(0.1)*IMAGES/sqrt(mean(var(IMAGES)));
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
42
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
43 save MY_IMAGES IMAGES
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
44
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
45 """
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
46
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
47 import numpy as np
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
48 def whiten(X, f0=None, n=None):
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
49 """
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
50 :param X: a 3D tensor n_images x rows x cols
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
51 :param f0: filter parameter (see docs)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
52 :param n: filter parameter (see docs)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
53
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
54 :returns: 3D tensor n_images x rows x cols of filtered images.
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
55 """
1416
James Bergstra <bergstrj@iro.umontreal.ca>
parents: 953
diff changeset
56 # May be mixed up with the size2 and size1s because matlab does things
James Bergstra <bergstrj@iro.umontreal.ca>
parents: 953
diff changeset
57 # differnetly
953
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
58 R, C = X.shape[-2:]
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
59 if R %2:
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
60 raise NotImplementedError()
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
61 if C %2:
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
62 raise NotImplementedError()
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
63
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
64 if f0 is None:
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
65 f0 = .4 * min(R,C)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
66 if n is None:
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
67 n = 4
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
68
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
69
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
70 fx,fy = np.mgrid[-R/2:R/2, -C/2:C/2]
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
71 rho=np.sqrt(fx**2 + fy**2)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
72 filt=rho * np.exp( - (rho/f0)**4)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
73 If = np.fft.fft2(X)
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
74 imagew=np.real(np.fft.ifft2(If * np.fft.fftshift(filt)))
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
75 assert imagew.shape == X.shape
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
76 return imagew
2eb98a740823 Added 'preprocessing' module
James Bergstra <bergstrj@iro.umontreal.ca>
parents:
diff changeset
77