fftwandopencv - FFTW@wiki

/**
* @brief basic form of fftw with opencv
*/
 
#include <stdio.h>
 
#include <cv.h>
#include <cxcore.h>
#include <highgui.h>
 
#define _USE_MATH_DEFINES
#include "math.h"
#include "fftw3.h"
 
#pragma comment( lib, "cv.lib" )
#pragma comment( lib, "cxcore.lib" )
#pragma comment( lib, "highgui.lib" )
#pragma comment( lib, "fftw3.lib" )
 
#define SIZEX 128 
#define SIZEY 64
 
int maGetPowerSpectol2D( IplImage *src, IplImage *dst );
int maShiftOrigin2D( IplImage *src, IplImage *dst );
 
int main( void )
{
	IplImage *img = NULL; // image for data manipulation
	IplImage *out = NULL; // image for result visualization
	fftw_complex *data = NULL;
	fftw_plan plan= NULL;
	CvScalar tmp;
	int i,j;
 
	data = (fftw_complex*)fftw_malloc( sizeof(fftw_complex) * SIZEX * SIZEY );
	img = cvCreateImageHeader( cvSize( SIZEX, SIZEY ), IPL_DEPTH_64F, 2 );
	img->imageData = (char*)data;
 
	for( j=0; j<SIZEY; j++ ){
		for( i=0; i<SIZEX; i++ ){
			tmp = cvScalar( 1 + 2*sin(2*M_PI*i/SIZEX) + 4*sin(4*M_PI*j/SIZEY), 0 );
			cvSet2D( img, j, i, tmp );
		}
	}
 
	plan = fftw_plan_dft_2d( SIZEY, SIZEX, data, data, FFTW_FORWARD, FFTW_ESTIMATE );
 
	fftw_execute( plan );
 
	out = cvCreateImage( cvSize( SIZEX, SIZEY ), IPL_DEPTH_8U, 1 );
 
	maGetPowerSpectol2D( img, out );
	maShiftOrigin2D( out, out );
 
	cvNamedWindow( "result" );
	cvShowImage( "result", out );
	cvWaitKey( 0 );
	cvDestroyWindow( "result" );
 
	if( plan ) fftw_destroy_plan( plan );
 
	if( out  ) cvReleaseImage( &out );
	if( img  ) cvReleaseImageHeader( &img );
	if( data ) fftw_free( data );
 
	return true;
}
/*
* Get the powerspectol
* src must be IPL_DEPTH_64F(double), dst must be IPL_DEPTH_8U(unsigned char)
* @param src source image
* @param dst destination image
*/
int maGetPowerSpectol2D( IplImage *src, IplImage *dst )
{
	IplImage *re = NULL;
	IplImage *im = NULL;
	CvSize size;
	double max,min;
 
	size = cvGetSize( src );
	re = cvCreateImage( size, IPL_DEPTH_64F, 1 );
	im = cvCreateImage( size, IPL_DEPTH_64F, 1 );
 
	// get the log scaled powerspectol : log( 1 + sqrt(re*re+im*im) )
	cvSplit( src, re, im, NULL, NULL );
	cvPow( re, re, 2.0 );
	cvPow( im, im, 2.0 );
	cvAdd( re, im, re );
	cvPow( re, re, 0.5 );
 
	cvAddS( re, cvScalarAll(1.0), re );
	cvLog( re, re );
 
	// scailing to 0-255
	cvMinMaxLoc( re, &min, &max );
	cvScale( re, re, 1.0 / (max - min), 1.0 * (-min) / (max - min ) );
	cvConvertScale( re, dst, 255 );
 
	if( re ) cvReleaseImage( &re );
	if( im ) cvReleaseImage( &im );
 
	return true;
}
/*
* Shift axis origin ( FFT result, and so on )
* This function applies for both in-place and out-place transform.
* Input and Output takes any size, channel, depth. They have only to be the same size, channel, depth
* @param src source image
* @param dst destination image
*/
int maShiftOrigin2D( IplImage *src, IplImage *dst )
{
  int cx, cy;  
  CvMat *tmp = NULL;
  CvMat q1stub, q2stub;
  CvMat q3stub, q4stub;
  CvMat d1stub, d2stub;
  CvMat d3stub, d4stub;
  CvMat *q1, *q2, *q3, *q4;
  CvMat *d1, *d2, *d3, *d4;
 
  CvSize size = cvGetSize (src);
  CvSize dst_size = cvGetSize (dst);
 
  if (src == dst)  tmp = cvCreateMat (size.height/2, size.width/2, cvGetElemType (src) );
 
  cx = size.width / 2;
  cy = size.height / 2;
 
  q1 = cvGetSubRect (src, &q1stub, cvRect (0,  0,  cx, cy));
  q2 = cvGetSubRect (src, &q2stub, cvRect (cx, 0,  cx, cy));
  q3 = cvGetSubRect (src, &q3stub, cvRect (cx, cy, cx, cy));
  q4 = cvGetSubRect (src, &q4stub, cvRect (0,  cy, cx, cy));
 
  d1 = cvGetSubRect (dst, &d1stub, cvRect (0,  0,  cx, cy));
  d2 = cvGetSubRect (dst, &d2stub, cvRect (cx, 0,  cx, cy));
  d3 = cvGetSubRect (dst, &d3stub, cvRect (cx, cy, cx, cy));
  d4 = cvGetSubRect (dst, &d4stub, cvRect (0,  cy, cx, cy));
 
  if (src != dst) {
	  cvCopy (q3, d1, 0);
	  cvCopy (q4, d2, 0);
	  cvCopy (q1, d3, 0);
	  cvCopy (q2, d4, 0);
  }
  else {
	  cvCopy (q3, tmp, 0);
	  cvCopy (q1, q3, 0);
	  cvCopy (tmp, q1, 0);
	  cvCopy (q4, tmp, 0);
	  cvCopy (q2, q4, 0);
	  cvCopy (tmp, q2, 0);
  }
 
  if( tmp ) cvReleaseMat( &tmp );
 
  return true;
}