fft2d.cc

torch tracking code, it is a good code

const char *help = "\
progname: fft2D.cc\n\
code2html: This program reads a pgm image and computes its FFT 2D.\n\
version: Torch3 vision2.0, 2005\n\
(c) Sebastien Marcel ([e-mail])\n";

#include "ImageGray.h"
#include "ipFFT2D.h"
#include "DiskXFile.h"
#include "CmdLine.h"

using namespace Torch;

int main(int argc, char **argv)
{
	char *image_filename;
	bool verbose;
	bool inverse;
	real highpass;
	real lowpass;
  	
	CmdLine cmd;
	cmd.setBOption("write log", false);
  	cmd.info(help);
  	cmd.addText("\nArguments:");
  	cmd.addSCmdArg("image filename", &image_filename, "image filename");
  	cmd.addText("\nOptions:");
  	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
  	cmd.addBCmdOption("-inverse", &inverse, false, "inverse");
  	cmd.addRCmdOption("-hpass", &highpass, 100.0, "high pass");
  	cmd.addRCmdOption("-lpass", &lowpass, -100.0, "low pass");
	cmd.read(argc, argv);


	// load input image
	Image *image_in = new ImageGray();
	image_in->setBOption("verbose", verbose);
	image_in->load(image_filename);

	if(verbose)
	{
		print("Image info:\n");
		print("   width = %d\n", image_in->width);
		print("   height = %d\n", image_in->height);
		print("   format = %s (%d)\n", image_in->coding, image_in->n_planes);
	}

	// check if the image if power of 2
	int new_width = (int) pow(2.0, ceil(log2((double)image_in->width)));
	int new_height = (int) pow(2.0, ceil(log2((double)image_in->height)));

	print("width = %d -> new width = %d\n", image_in->width, new_width);
	print("height = %d -> new height = %d\n", image_in->height, new_height);

	// re-copy the input image into an image of size power of 2
	ImageGray *newimage = new ImageGray(new_width, new_height);
	pasteGray(image_in->data, 0, 0, image_in->width, image_in->height, newimage->data, new_width, new_height);
	newimage->updatePixmapFromData();
	newimage->save("new.pgm");

	Image *image = newimage;

	//
	// FFT2D
	ipCore *fft2d = new ipFFT2D(image->width, image->height);
	fft2d->setBOption("verbose", verbose);

	print("Computing FFT 2D ...\n");
	fft2d->process(image);

	// print coefficients
	real *fft_r = fft2d->seq_out->frames[0];
	real *fft_i = fft2d->seq_out->frames[1];

	real energy = 0;
	int n_hpass = 0;
	int n_lpass = 0;
	if(verbose) print("FFT 2D:\n");	
	for(int i = 0; i < fft2d->seq_out->frame_size; i++)
	{
		if(verbose) print("[%d] = %g (%g)\n", i, fft_r[i], fft_i[i]);

		energy += fft_r[i] * fft_r[i];

		if(i > 0)
		{
			if(fft_r[i] > highpass) n_hpass++;
			if(fft_r[i] < lowpass) n_lpass++;
		}
	}

	energy = sqrt(energy);
	real log_energy = energy < 1.0 ? 0.0 : (real) log(energy); 

	print("DC frequency = %g + %gi\n", fft_r[0], fft_i[0]);
	print("Nyquist frequency = %g + %gi\n", fft_r[fft2d->seq_out->frame_size/2], fft_i[fft2d->seq_out->frame_size/2]);
	print("log Energy = %g\n", log_energy);
	print("number of high passed = %d\n", n_hpass);
	print("number of low passed = %d\n", n_lpass);


	// 
	// IFFT2D
	ipCore *ifft2d = NULL;

	if(inverse)
	{
		ifft2d = new ipFFT2D(image->width, image->height, true);
		ifft2d->setBOption("verbose", verbose);
	
		print("Computing IFFT 2D ...\n");
		ifft2d->process(fft2d->seq_out);

		// compute reconstruction error
		real rmse = 0.0;
		
		if(verbose) print("IFFT 2D:\n");	
		for(int i = 0; i < ifft2d->seq_out->frame_size ; i++)
		{
			if(verbose) print("[%d] = %g (%g)\n", i, ifft2d->seq_out->frames[0][i], ifft2d->seq_out->frames[1][i]);
		   	real z = image->frames[0][i] - ifft2d->seq_out->frames[0][i];
			rmse += z * z;
		}
		print("RMSE = %g\n", rmse);


		// save the reconstructed image
		ImageGray *image_out = new ImageGray(image_in->width, image_in->height);
		cropGray(ifft2d->seq_out->frames[0], image->width, image->height, 0, 0, image_in->width, image_in->height, image_out->data);
		image_out->updatePixmapFromData();		
		image_out->save("inverse.pgm");

		delete ifft2d;
		delete image_out;
	}
	
	delete fft2d;
	delete image_in;
	delete newimage;

	return(0);
}