gabor2d.cc

  1 const char *help = "\
  2 progname: gabor2D.cc\n\
  3 code2html: This program reads a pgm image and applies a Gabor filter.\n\
  4 version: Torch3 vision2.0, 2005\n\
  5 (c) Sebastien Marcel (marcel@idiap.ch)\n";
  6 
  7 #include "ImageGray.h"
  8 #include "ipFFT2D.h"
  9 #include "ipGabor2D.h"
 10 #include "matrix_addons.h"
 11 #include "DiskXFile.h"
 12 #include "CmdLine.h"
 13 
 14 using namespace Torch;
 15 
 16 int main(int argc, char **argv)
 17 {
 18 	char *image_filename;
 19 	bool verbose;
 20 	int scale;
 21 	int orientation;
 22 	real Ul;
 23 	real Uh;
 24 	int side;
 25 
 26   	
 27 	CmdLine cmd;
 28 	cmd.setBOption("write log", false);
 29   	cmd.info(help);
 30   	cmd.addText("\nArguments:");
 31   	cmd.addSCmdArg("image filename", &image_filename, "image filename");
 32   	cmd.addText("\nOptions:");
 33   	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 34   	cmd.addICmdOption("-scale", &scale, 5, "scale");
 35   	cmd.addICmdOption("-orientation", &orientation, 4, "orientation");
 36   	cmd.addRCmdOption("-ul", &Ul, 0.1, "lowest spatial frequency");
 37   	cmd.addRCmdOption("-uh", &Uh, 0.9, "highest spatial frequency");
 38   	cmd.addICmdOption("-side", &side, 12, "size of the filter");
 39 	cmd.read(argc, argv);
 40 
 41 	Allocator *allocator = new Allocator;
 42 
 43 	// load input image
 44 	Image *image_in = new(allocator) ImageGray();
 45 	image_in->setBOption("verbose", verbose);
 46 	image_in->load(image_filename);
 47 
 48 	if(verbose)
 49 	{
 50 		print("Image info:\n");
 51 		print("   width = %d\n", image_in->width);
 52 		print("   height = %d\n", image_in->height);
 53 		print("   format = %s (%d)\n", image_in->coding, image_in->n_planes);
 54 	}
 55 
 56 	Image *image = image_in;
 57 
 58 	//
 59 	// Gabor2D
 60 	
 61 	int size = 2*side+1;
 62 
 63 	print("Gabor filters:\n");
 64 	print(" scale = %d\n", scale);
 65 	print(" orientation = %d\n", orientation);
 66 	print(" highest spatial frequency = %g\n", Uh);
 67 	print(" lowest spatial frequency = %g\n", Ul);
 68 	print(" filter dimension = %d -> %dx%d\n", side, size, size);
 69 	
 70 	//
 71 	// Creates a Gabor filter bank and save them as pgm
 72 	ipGabor2D **gaborFilterBank = NULL;
 73 	gaborFilterBank = (ipGabor2D **)allocator->alloc(scale*orientation*sizeof(ipGabor2D *));
 74 
 75 	// allocate output images of Gabor wavelets
 76 	ImageGray *image_r_out = new(allocator) ImageGray(size*orientation, size*scale);
 77 	ImageGray *image_i_out = new(allocator) ImageGray(size*orientation, size*scale);
 78 	Mat *mgabor_r = new(allocator) Mat(size*scale, size*orientation);
 79 	Mat *mgabor_i = new(allocator) Mat(size*scale, size*orientation);
 80 	Mat *tmp_r = new(allocator) Mat(size, size);
 81 	Mat *tmp_i = new(allocator) Mat(size, size);
 82 
 83 	mgabor_r->zero();
 84 	mgabor_i->zero();
 85 
 86 	int f = 0;
 87 
 88 	for(int s=0;s<scale;s++)
 89 	   for(int n=0;n<orientation;n++)
 90 	   {
 91 		if(verbose) print("Computing Gabor 2D (s=%d o=%d) ...\n", s+1, n+1);
 92 
 93 		gaborFilterBank[f] = new(allocator) ipGabor2D(size, size, s+1, n+1, Ul, Uh, scale, orientation);
 94 		gaborFilterBank[f]->setBOption("verbose", verbose);
 95 
 96 		if(verbose)
 97 		{
 98 			// print coefficients
 99 			real *gabor_r = gaborFilterBank[f]->seq_out->frames[0];
100 			real *gabor_i = gaborFilterBank[f]->seq_out->frames[1];
101 
102 			print("Gabor 2D:\n");	
103 			for(int i = 0; i < gaborFilterBank[f]->seq_out->frame_size; i++)
104 				print("[%d] = %g (%g)\n", i, gabor_r[i], gabor_i[i]);
105 		}
106 
107 		mxCopy(tmp_r, gaborFilterBank[f]->Gr, 0, 0, 0, 0, 2*side, 2*side);
108 		mxCopy(tmp_i, gaborFilterBank[f]->Gi, 0, 0, 0, 0, 2*side, 2*side);
109 
110 		mxScale(tmp_r);
111 		mxScale(tmp_i);
112 
113 		mxCopy(mgabor_r, tmp_r, s*size, n*size, 0, 0, 2*side, 2*side);
114 		mxCopy(mgabor_i, tmp_i, s*size, n*size, 0, 0, 2*side, 2*side);
115 		
116 		f++;
117 	   }
118 	
119 	for(int i = 0 ; i < image_r_out->height ; i++)
120 			for(int j = 0 ; j < image_r_out->width ; j++)
121 			{
122 				image_r_out->data[i*image_r_out->width + j] = mgabor_r->ptr[i][j];
123 				image_i_out->data[i*image_i_out->width + j] = mgabor_i->ptr[i][j];
124 			}
125 	
126 	image_r_out->updatePixmapFromData();		
127 	image_r_out->save("gabor_r.pgm");
128 
129 	image_i_out->updatePixmapFromData();		
130 	image_i_out->save("gabor_i.pgm");
131 
132 	allocator->free(mgabor_r);
133 	allocator->free(mgabor_i);
134 	allocator->free(image_r_out);
135 	allocator->free(image_i_out);
136 	
137 
138 	//
139 	// Filter the input image using the filter bank in the FT domain
140 	
141 	// check if the image if power of 2
142 	int new_width = (int) pow(2.0, ceil(log2((double)(image_in->width+2.0*side))));
143 	int new_height = (int) pow(2.0, ceil(log2((double)(image_in->height+2.0*side))));
144 
145 	print("Updating image size for FFT\n");
146 	print(" width = %d -> new width = %d\n", image_in->width, new_width);
147 	print(" height = %d -> new height = %d\n", image_in->height, new_height);
148 
149 	// allocate the new image
150 	ImageGray *newimage = new(allocator) ImageGray(new_width, new_height);
151 
152 	// re-copy the input image into an image of size power of 2
153 	pasteGray(image_in->data, side, side, image_in->width, image_in->height, newimage->data, new_width, new_height);	
154 	newimage->updatePixmapFromData();
155 
156 	if(verbose) newimage->save("new.pgm");
157 
158 	// the image is the new image
159 	image = newimage;
160 	
161 	// create a FFT ipCore
162 	ipFFT2D *fft2d = new(allocator) ipFFT2D(image->width, image->height);
163 	fft2d->setBOption("verbose", verbose);
164 
165 	// create a FFT ipCore in IFFT mode
166 	ipFFT2D *ifft2d = new(allocator) ipFFT2D(image->width, image->height, true);
167 	ifft2d->setBOption("verbose", verbose);
168 
169 	// allocate snapshot images
170 	image_r_out = new(allocator) ImageGray(image_in->width*orientation, image_in->height*scale);
171 	image_i_out = new(allocator) ImageGray(image_in->width*orientation, image_in->height*scale);
172 
173 	// allocate the FFT of the input image
174 	Sequence *seqfft = new(allocator) Sequence(2, image->width * image->height);
175 
176 	// compute the FFT of the input image
177 	fft2d->process(image);
178 
179 	// copy the FFT of the input image into a sequence
180 	seqfft->copy(fft2d->seq_out);
181 	
182 	// loop on Gabor filters of the filter bank
183 	for(int f=0 ; f < scale*orientation ; f++)
184 	{
185 	   	int s = gaborFilterBank[f]->s-1;
186 	   	int n = gaborFilterBank[f]->n-1;
187 
188 		if(verbose) print("Applying Gabor 2D (s=%d o=%d) ...\n", s+1, n+1);
189 
190 		// filter the image in the FT domain 
191 		gaborFilterBank[f]->setFFT(fft2d, ifft2d, image_in->width, image_in->height);
192 
193 		// filter the image
194 		gaborFilterBank[f]->process(seqfft);
195 
196 		// scale the outputs (filtered images) both real and imaginary parts
197 		scaleGray(image_in->width, image_in->height, gaborFilterBank[f]->seq_out->frames[0]);
198 		scaleGray(image_in->width, image_in->height, gaborFilterBank[f]->seq_out->frames[1]);
199 	
200 		// save individual filtered images if needed
201 		if(verbose)
202 		{
203 			char str[250];
204 			ImageGray *image_out = new ImageGray(image_in->width, image_in->height);
205 			
206 			image_out->copyFrom(image_in->width, image_in->height, gaborFilterBank[f]->seq_out->frames[0], "gray");
207 			sprintf(str, "filter-R-s%d-o%d.pgm", s, n);
208 			image_out->save(str);
209 			
210 			image_out->copyFrom(image_in->width, image_in->height, gaborFilterBank[f]->seq_out->frames[1], "gray");
211 			sprintf(str, "filter-I-s%d-o%d.pgm", s, n);
212 			image_out->save(str);
213 	
214 			delete image_out;
215 		}
216 
217 		// paste the scaled output images into the snapshot images both real and imaginary parts
218 		pasteGray(gaborFilterBank[f]->seq_out->frames[0], n*image_in->width, s*image_in->height, image_in->width, image_in->height, image_r_out->data, image_r_out->width, image_r_out->height);
219 		pasteGray(gaborFilterBank[f]->seq_out->frames[1], n*image_in->width, s*image_in->height, image_in->width, image_in->height, image_i_out->data, image_i_out->width, image_i_out->height);
220 	}
221 		
222 	image_r_out->updatePixmapFromData();		
223 	image_r_out->save("filter_r.pgm");
224 
225 	image_i_out->updatePixmapFromData();		
226 	image_i_out->save("filter_i.pgm");
227 
228 	delete allocator;
229 
230 	return(0);
231 }