#define USE_FFTW_DLL

#include <windows.h>
#include "common.h"

#include "fftw.h"



// Function to display the Fourier Transform

Image_data * FourierTransform(Image_data *piData)
{
	int i;
	Image_data *pproData;
	fftwnd_plan plan;  
	int width, height;
	double magme; 
	float *magnitude;
	
	
	FFTW_COMPLEX *in;
	FFTW_COMPLEX *out;
	
	width = piData->n_cols;
	height = piData->n_rows;
	
	in = calloc(width*height,sizeof(FFTW_COMPLEX));
	out = calloc(width*height,sizeof(FFTW_COMPLEX));
	magnitude = calloc(width*height,sizeof(float));
	

	// ----------- Memory allocation stuff --------------//
	
	pproData = malloc(sizeof(Image_data));
	pproData->image_G_ptr = malloc(piData->n_rows * piData->n_cols);
	pproData->n_cols = piData->n_cols;
	pproData->n_rows = piData->n_rows;
	lstrcpy(pproData->filename, "Fourier Transform");
	
	// -------------------------------------------------//
	
	
	plan = fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_USE_WISDOM);
	if (plan==NULL)    
	{      
		MessageBox(0, "Error Creating Plan", "", MB_OK);
		return NULL;    
	}  
	
	
	
	// Fill in the 'in' array with intensity values
	
	for (i=0; i < width*height; i++)     
	{     
		c_re(in[i]) = (int)piData->image_G_ptr[i];
		c_im(in[i])=0;    
	}  
	
	fftwnd(plan,1,in,1,0,out,1,0);
	


	// Scale the transform so that it becomes visible

	for (i=0; i < width*height; i++)     
	{
		magme = sqrt(c_re(out[i])*c_re(out[i]) + c_im(out[i])*c_im(out[i]));
		magme = magme / 100.0;
		if (magme > 255.0)
			magme = 255.0;
		pproData->image_G_ptr[i]=(unsigned char)(magme);
 	}



	// Free up allocated memory	
	
	fftwnd_destroy_plan(plan);  
	fftw_free(in);  
	fftw_free(out);
	free(magnitude);
	
	return pproData;
	
}



// The Low Pass Filter implementation

Image_data * LowPass(Image_data *piData, int x, int y)
{
	int i;
	Image_data *pproData;
	fftwnd_plan plan;  
	int width, height;
	double magme, scale=0; 
	double threshold, d1, d2, d3, d4;
	
	
	FFTW_COMPLEX *in;
	FFTW_COMPLEX *out;
	
	width = piData->n_cols;
	height = piData->n_rows;
	
	in = calloc(width*height,sizeof(FFTW_COMPLEX));
	out = calloc(width*height,sizeof(FFTW_COMPLEX));
	
	// ----------- Memory allocation stuff --------------//
	
	pproData = malloc(sizeof(Image_data));
	pproData->image_G_ptr = malloc(piData->n_rows * piData->n_cols);
	pproData->n_cols = piData->n_cols;
	pproData->n_rows = piData->n_rows;
	lstrcpy(pproData->filename, "Low Passed Image");
	
	// -------------------------------------------------//
	
	
	plan = fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_USE_WISDOM);
	if (plan==NULL)    
	{      
		MessageBox(0, "Error Creating Plan", "", MB_OK);
		return NULL;    
	}  
	
	

	// Fill 'in' array with intensity values of the pixels
	
	for (i=0; i < width*height; i++)     
	{     
		c_re(in[i]) = (int)piData->image_G_ptr[i];
		c_im(in[i])=0;    
	}  
	
	fftwnd(plan,1,in,1,0,out,1,0);
		// Compute transform
	
	fftwnd_destroy_plan(plan);
	


	// Compute the threshold which the the shortest distance
	// to any of the 4 corners

	if ((x<width/2) && (y<height/2))		// top left quadrant
		threshold = sqrt(x*x + y*y);
	else if ((x<width/2) && (y>height/2))	// bottom left quadrant
		threshold = sqrt(x*x + (height-y)*(height-y));
    	else if ((x>width/2) && (y<height/2))	// top right quadrant	
		threshold = sqrt((width-x)*(width-x) + y*y);
	else									// bottom right quadrant
		threshold = sqrt((width-x)*(width-x) + (height-y)*(height-y));
	
	
	
	// Calculate distances from the 4 corners and set 
	// appropriate pixels to 0.0
	
	for (x = 0; x < width; x++)
		for (y = 0; y < height; y++)
		{
			d1 = sqrt(x*x + y*y);
			d2 = sqrt((width-x)*(width-x) + y*y);
			d3 = sqrt(x*x + (height-y)*(height-y));
			d4 = sqrt((width-x)*(width-x) + (height-y)*(height-y));
			
			if ((d1 > threshold) && (d2 > threshold) && (d3 > threshold) && (d4 > threshold))
			{
				c_re(out[y*width+x]) = 0.0;
				c_im(out[y*width+x]) = 0.0;
			}
		}
		
		
		plan = fftw2d_create_plan(height, width, FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_USE_WISDOM);
		if (plan == NULL)    
		{      
			MessageBox(0, "Error Creating Plan", "", MB_OK);
			return NULL;    
		} 
		
		fftwnd(plan,1,out,1,0,in,1,0);		// Inverse Transform
		
		scale = c_re(in[0]);
			// Compute appropriate scale
		for (i=0; i < width*height; i++)
		{
			if (scale < c_re(in[i]))
				scale = c_re(in[i]);
		}
		
		
		// Retrieve pixel values from inverse transform
		
		for (i=0; i < width*height; i++)
		{
			magme = sqrt(c_re(in[i])*c_re(in[i]) + c_im(in[i])*c_im(in[i]));
			pproData->image_G_ptr[i]=(unsigned char)((magme/scale)*255.0);
		}
		
		
		
		// Free up allocated memory
		
		fftwnd_destroy_plan(plan);
		fftw_free(in);  
		fftw_free(out);
		
		
		return pproData;
		
}



// The High Pass Filter implementation

Image_data * HighPass(Image_data *piData, int x, int y)
{
	int i, j;
	Image_data *pproData;
	fftwnd_plan plan;  
	int width, height;
	double magme, scale=0; 
	double threshold, d1, d2, d3, d4;
	
	
	FFTW_COMPLEX *in;
	FFTW_COMPLEX *out;
	
	width = piData->n_cols;
	height = piData->n_rows;
	
	in = calloc(width*height,sizeof(FFTW_COMPLEX));
	out = calloc(width*height,sizeof(FFTW_COMPLEX));
	
	// ----------- Memory allocation stuff --------------//
	
	pproData = malloc(sizeof(Image_data));
	pproData->image_G_ptr = malloc(piData->n_rows * piData->n_cols);
	pproData->n_cols = piData->n_cols;
	pproData->n_rows = piData->n_rows;
	lstrcpy(pproData->filename, "High Passed Image");
	
	// -------------------------------------------------//
	
	
	plan = fftw2d_create_plan(height, width, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_USE_WISDOM);    
	if (plan == NULL)    
	{      
		MessageBox(0, "Error Creating Plan", "", MB_OK);
		return NULL;    
	}  
	
	
	// Fill 'in' array with pixel intensiy values
	
	for (i=0; i < width*height; i++)     
	{     
		c_re(in[i]) = (int)piData->image_G_ptr[i];
		c_im(in[i])=0;    
	}  
	
	fftwnd(plan,1,in,1,0,out,1,0);
		// Compute transform
	
	fftwnd_destroy_plan(plan);
	


	// Compute the threshold which the the shortest distance
	// to any of the 4 corners
	
	if ((x<width/2) && (y<height/2))		// top left quadrant
		threshold = sqrt(x*x + y*y);
	else if ((x<width/2) && (y>height/2))	// bottom left quadrant
		threshold = sqrt(x*x + (height-y)*(height-y));
  	else if ((x>width/2) && (y<height/2))	// top right quadrant	
		threshold = sqrt((width-x)*(width-x) + y*y);
	else									// bottom right quadrant
		threshold = sqrt((width-x)*(width-x) + (height-y)*(height-y));
	

	
	// Check distances from the 4 corners
	// and set approproate pixels to 0.0
	
	for (i = 0; i < width; i++)
		for (j = 0; j < height; j++)
		{
			d1 = sqrt(i*i + j*j);
			d2 = sqrt((width-i)*(width-i) + j*j);
			d3 = sqrt(i*i + (height-j)*(height-j));
			d4 = sqrt((width-i)*(width-i) + (height-j)*(height-j));
			
			if ((d1 < threshold) || (d2 < threshold) || (d3 < threshold) || (d4 < threshold))
			{
				c_re(out[j*width+i]) = 0.0;
				c_im(out[j*width+i]) = 0.0;
			}
		}
		
		
		plan = fftw2d_create_plan(height, width, FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_USE_WISDOM);
		if (plan == NULL)    
		{      
			MessageBox(0, "Error Creating Plan", "", MB_OK);
			return NULL;    
		} 
		
		fftwnd(plan,1,out,1,0,in,1,0);		// Inverse Transform
		
		scale = c_re(in[0]);
			// Compute scale
		for (i=0; i < width*height; i++)
		{
			if (scale < c_re(in[i]))
				scale = c_re(in[i]);
		}
		
		
		
		// Retrieve pixel intensities from the inverse transform
		
		for (i=0; i < width*height; i++)
		{
			magme = sqrt(c_re(in[i])*c_re(in[i]) + c_im(in[i])*c_im(in[i]));
			pproData->image_G_ptr[i]=(unsigned char)((magme/scale)*255.0);
		}
		
		
		
		// Free up allocated memory
		
		fftwnd_destroy_plan(plan);
		fftw_free(in);  
		fftw_free(out);
		
		
		return pproData;
		
}



// Function to execute High Pass on certain area of image

Image_data * LocalHighPass(Image_data *piData, int startx, int starty, int endx, int endy, int x, int y)
{
	Image_data *ptempData, *psubimage, *pproData;
	int i, j, k, width, height;
	
	
	width = piData->n_cols;
	height = piData->n_rows;
	
	if (startx > endx)	// If endx < startx, swap them, otherwise for loop will fail
	{
		i = startx;
		startx = endx;
		endx = i;
	}
	
	
	if (starty > endy)	// If endy < starty, swap them, otherwise for loop will fail
	{
		i = starty;
		starty = endy;
		endy = i;
	}
	
	
	// ----------- Memory allocation stuff --------------//
	
	pproData = malloc(sizeof(Image_data));
	pproData->image_G_ptr = malloc(width * height);
	memcpy(pproData->image_G_ptr, piData->image_G_ptr, width*height);
	pproData->n_cols = piData->n_cols;
	pproData->n_rows = piData->n_rows;
	lstrcpy(pproData->filename, "High Passed (Local)");
	
	psubimage = malloc(sizeof(Image_data));
	psubimage->image_G_ptr = malloc((endx-startx) * (endy-starty));
	psubimage->n_cols = endx-startx;
	psubimage->n_rows = endy-starty;
	lstrcpy(psubimage->filename, "Subimage");
	
	//---------------------------------------------------//
	
	
	
	// Create a subimage to pass to the HighPass() function
	
	k = 0;
	for (j = starty; j < endy; j++)
		for (i = startx; i < endx; i++)
		{
			psubimage->image_G_ptr[k] = piData->image_G_ptr[j*width + i];
			++k;
		}
		
	ptempData = HighPass(psubimage, x, y);
		
	

	// Superimpose processed subimage on the original image

	k = 0;
	for (j = starty; j < endy; j++)
		for (i = startx; i < endx; i++)
		{
			pproData->image_G_ptr[j*width + i] = ptempData->image_G_ptr[k];
			++k;
		}
			
	
		
	// Free up allocated memory
		
	free(psubimage->image_G_ptr);
	free(psubimage);
			
	free(ptempData->image_G_ptr);
	free(ptempData);
			
	return pproData;
}



// Function to execute Low Pass on certain area of image

Image_data * LocalLowPass(Image_data *piData, int startx, int starty, int endx, int endy, int x, int y)
{
	Image_data *ptempData, *psubimage, *pproData;
	int i, j, k, width, height;
	
	
	width = piData->n_cols;
	height = piData->n_rows;
	
	if (startx > endx)	// If endx < startx, swap them, otherwise for loop will fail
	{
		i = startx;
		startx = endx;
		endx = i;
	}
	
	
	if (starty > endy)	// If endy < starty, swap them, otherwise for loop will fail
	{
		i = starty;
		starty = endy;
		endy = i;
	}
	
	
	// ----------- Memory allocation stuff --------------//
	
	pproData = malloc(sizeof(Image_data));
	pproData->image_G_ptr = malloc(width * height);
	memcpy(pproData->image_G_ptr, piData->image_G_ptr, width*height);
	pproData->n_cols = piData->n_cols;
	pproData->n_rows = piData->n_rows;
	lstrcpy(pproData->filename, "High Passed (Local)");
	
	psubimage = malloc(sizeof(Image_data));
	psubimage->image_G_ptr = malloc((endx-startx) * (endy-starty));
	psubimage->n_cols = endx-startx;
	psubimage->n_rows = endy-starty;
	lstrcpy(psubimage->filename, "Subimage");
	
	//----------------------------------------------------//
	
	

	// Create a subimage to pass to the LowPass() function

	k = 0;
	for (j = starty; j < endy; j++)
		for (i = startx; i < endx; i++)
		{
			psubimage->image_G_ptr[k] = piData->image_G_ptr[j*width + i];
			++k;
		}
		
	ptempData = LowPass(psubimage, x, y);
	
	

	// Superimpose processed subimage on the original image

	k = 0;
	for (j = starty; j < endy; j++)
		for (i = startx; i < endx; i++)
		{
			pproData->image_G_ptr[j*width + i] = ptempData->image_G_ptr[k];
			++k;
		}
		
		

	// Free up allocated memory

	free(psubimage->image_G_ptr);
	free(psubimage);
			
	free(ptempData->image_G_ptr);
	free(ptempData);
			
	return pproData;
}