fft2d.cc

  1 const char *help = "\
  2 progname: fft2D.cc\n\
  3 code2html: This program reads a pgm image and computes its FFT 2D.\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 "DiskXFile.h"
 10 #include "CmdLine.h"
 11 
 12 using namespace Torch;
 13 
 14 int main(int argc, char **argv)
 15 {
 16 	char *image_filename;
 17 	bool verbose;
 18 	bool inverse;
 19 	real highpass;
 20 	real lowpass;
 21   	
 22 	CmdLine cmd;
 23 	cmd.setBOption("write log", false);
 24   	cmd.info(help);
 25   	cmd.addText("\nArguments:");
 26   	cmd.addSCmdArg("image filename", &image_filename, "image filename");
 27   	cmd.addText("\nOptions:");
 28   	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 29   	cmd.addBCmdOption("-inverse", &inverse, false, "inverse");
 30   	cmd.addRCmdOption("-hpass", &highpass, 100.0, "high pass");
 31   	cmd.addRCmdOption("-lpass", &lowpass, -100.0, "low pass");
 32 	cmd.read(argc, argv);
 33 
 34 
 35 	// load input image
 36 	Image *image_in = new ImageGray();
 37 	image_in->setBOption("verbose", verbose);
 38 	image_in->load(image_filename);
 39 
 40 	if(verbose)
 41 	{
 42 		print("Image info:\n");
 43 		print("   width = %d\n", image_in->width);
 44 		print("   height = %d\n", image_in->height);
 45 		print("   format = %s (%d)\n", image_in->coding, image_in->n_planes);
 46 	}
 47 
 48 	// check if the image if power of 2
 49 	int new_width = (int) pow(2.0, ceil(log2((double)image_in->width)));
 50 	int new_height = (int) pow(2.0, ceil(log2((double)image_in->height)));
 51 
 52 	print("width = %d -> new width = %d\n", image_in->width, new_width);
 53 	print("height = %d -> new height = %d\n", image_in->height, new_height);
 54 
 55 	// re-copy the input image into an image of size power of 2
 56 	ImageGray *newimage = new ImageGray(new_width, new_height);
 57 	pasteGray(image_in->data, 0, 0, image_in->width, image_in->height, newimage->data, new_width, new_height);
 58 	newimage->updatePixmapFromData();
 59 	newimage->save("new.pgm");
 60 
 61 	Image *image = newimage;
 62 
 63 	//
 64 	// FFT2D
 65 	ipCore *fft2d = new ipFFT2D(image->width, image->height);
 66 	fft2d->setBOption("verbose", verbose);
 67 
 68 	print("Computing FFT 2D ...\n");
 69 	fft2d->process(image);
 70 
 71 	// print coefficients
 72 	real *fft_r = fft2d->seq_out->frames[0];
 73 	real *fft_i = fft2d->seq_out->frames[1];
 74 
 75 	real energy = 0;
 76 	int n_hpass = 0;
 77 	int n_lpass = 0;
 78 	if(verbose) print("FFT 2D:\n");	
 79 	for(int i = 0; i < fft2d->seq_out->frame_size; i++)
 80 	{
 81 		if(verbose) print("[%d] = %g (%g)\n", i, fft_r[i], fft_i[i]);
 82 
 83 		energy += fft_r[i] * fft_r[i];
 84 
 85 		if(i > 0)
 86 		{
 87 			if(fft_r[i] > highpass) n_hpass++;
 88 			if(fft_r[i] < lowpass) n_lpass++;
 89 		}
 90 	}
 91 
 92 	energy = sqrt(energy);
 93 	real log_energy = energy < 1.0 ? 0.0 : (real) log(energy); 
 94 
 95 	print("DC frequency = %g + %gi\n", fft_r[0], fft_i[0]);
 96 	print("Nyquist frequency = %g + %gi\n", fft_r[fft2d->seq_out->frame_size/2], fft_i[fft2d->seq_out->frame_size/2]);
 97 	print("log Energy = %g\n", log_energy);
 98 	print("number of high passed = %d\n", n_hpass);
 99 	print("number of low passed = %d\n", n_lpass);
100 
101 
102 	// 
103 	// IFFT2D
104 	ipCore *ifft2d = NULL;
105 
106 	if(inverse)
107 	{
108 		ifft2d = new ipFFT2D(image->width, image->height, true);
109 		ifft2d->setBOption("verbose", verbose);
110 	
111 		print("Computing IFFT 2D ...\n");
112 		ifft2d->process(fft2d->seq_out);
113 
114 		// compute reconstruction error
115 		real rmse = 0.0;
116 		
117 		if(verbose) print("IFFT 2D:\n");	
118 		for(int i = 0; i < ifft2d->seq_out->frame_size ; i++)
119 		{
120 			if(verbose) print("[%d] = %g (%g)\n", i, ifft2d->seq_out->frames[0][i], ifft2d->seq_out->frames[1][i]);
121 		   	real z = image->frames[0][i] - ifft2d->seq_out->frames[0][i];
122 			rmse += z * z;
123 		}
124 		print("RMSE = %g\n", rmse);
125 
126 
127 		// save the reconstructed image
128 		ImageGray *image_out = new ImageGray(image_in->width, image_in->height);
129 		cropGray(ifft2d->seq_out->frames[0], image->width, image->height, 0, 0, image_in->width, image_in->height, image_out->data);
130 		image_out->updatePixmapFromData();		
131 		image_out->save("inverse.pgm");
132 
133 		delete ifft2d;
134 		delete image_out;
135 	}
136 	
137 	delete fft2d;
138 	delete image_in;
139 	delete newimage;
140 
141 	return(0);
142 }