// MorphFFT.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

double absolute(double a,double b){return sqrt(a*a+b*b);}

int main(int argc, char* argv[])
{
	char            *image = argv[1];
	IplImage        *img1 = 0;
    IplImage        *img2 = 0;
	uchar			*img1_data;
	uchar			*img2_data;
	
    fftw_complex    *data_in;
    fftw_complex    *fft;        
    fftw_complex    *ifft;    
    fftw_plan       plan_f;
    fftw_plan       plan_b;
    
	double          *buffer;
    int             width, height, step;
    int             i, j, k;
    
    /* load original image */
    img1 = cvLoadImage( image, CV_LOAD_IMAGE_GRAYSCALE );
    
	/* always check */
    if ( img1 == 0 ) {
        fprintf( stderr, "Cannot load file %s!\n", argv[1] );
        return 1;
    }
	
	/* create new image for IFFT result */
	img2 = cvCreateImage( cvSize( img1->width, img1->height ), IPL_DEPTH_8U, 1 );
	
	/* get image properties */
    width  	  = img1->width;
    height 	  = img1->height;
    step	  = img1->widthStep;
	img1_data = ( uchar* ) img1->imageData;
	img2_data = ( uchar* ) img2->imageData;    
    
    /*initialize arrays for fftw operations */
    data_in = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
    fft     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
    ifft    = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * width * height );
	buffer  = new double[width * height];

	/* display images */
    cvNamedWindow( "Original", CV_WINDOW_AUTOSIZE );
    cvNamedWindow( "FFT", CV_WINDOW_AUTOSIZE );
	cvShowImage( "Original", img1 );

	/* create plans */
    plan_f = fftw_plan_dft_1d( width * height, data_in, fft,  FFTW_FORWARD,  FFTW_ESTIMATE );
    plan_b = fftw_plan_dft_1d( width * height, fft,     ifft, FFTW_BACKWARD, FFTW_ESTIMATE );
    
    /* load img1's data to fftw input */
    for( i = 0, k = 0 ; i < height ; i++ ) {
        for( j = 0 ; j < width ; j++ ) {
            data_in[k][0] = ( double )img1_data[i * step + j];
            data_in[k][1] = 0.0;
            k++;
        }
    }
    
    /* perform FFT */
    fftw_execute( plan_f );
    
	double min=DBL_MAX,max=-DBL_MAX,t=0;
	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++ ) {
			t = log(absolute(fft[k][0],fft[k][1]));
			buffer[i*width+j]=t;
			if(t<min)min=t;
			if(t>max)max=t;
			k++;
		}
	}
	
	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++ ) {
			t =((((buffer[i*width+j] - min) * (255 - 0)) / (max - min)) + 0);
			img2_data[i * step + j] = ( uchar )t;
		}
	}

	cvShowImage( "FFT", img2 );
	cvWaitKey( 0 );

	// Frequency domain smoothing
	//cvSmooth(img2,img2 );

	cvDilate(img2,img2);
	cvShowImage( "FFT", img2 );
	cvWaitKey( 0 );

	for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++ ) {
			t = ( double )img2_data[i * step + j];
			t=((((t - 0) * (max - min)) / (255 - 0)) + min);
			//t = ((((t - min2) * (max - min)) / (max2 - min2)) + min);
			t= exp(t)/absolute(fft[k][0],fft[k][1]);
			fft[k][0]*=t;fft[k][1]*=t;
			k++;
		}
	}


    /* perform IFFT */
    fftw_execute( plan_b );
    
    /* normalize IFFT result */
	//double min2=DBL_MAX,max2=-DBL_MAX;
    for( i = 0 ; i < ( width * height ) ; i++ ) {
        ifft[i][0] /= ( double )( width * height );
		//t=ifft[i][0];
		//if(t<min2)min2=t;
		//if(t>max2)max2=t;
    }
    
    /* copy IFFT result to img2's data */
    for( i = 0, k = 0 ; i < height ; i++ ) {
		for( j = 0 ; j < width ; j++ ) {
			t=( uchar )ifft[k++][0];
			//t =((((t - min2) * (255 - 0)) / (max2 - min2)) + 0);
			img2_data[i * step + j] =(uchar)t;
		}
	}
    
	//cvNormalize(img2,img2,max2,min2);

    cvShowImage( "FFT", img2 );
	cvWaitKey( 0 );

	/* free memory */
	cvSaveImage("bilateral1.jpg",img2);
	cvDestroyWindow( "Original" );
    cvDestroyWindow( "FFT" );
	cvReleaseImage( &img1 );
	cvReleaseImage( &img2 );
	
    fftw_destroy_plan( plan_f );
    fftw_destroy_plan( plan_b );
    fftw_free( data_in );
    fftw_free( fft );
    fftw_free( ifft );
	return 0;
}