//
//  main.m
//  ImageApp
//
//  Created by seppke on 20.11.08.
//  Copyright __MyCompanyName__ 2008. All rights reserved.
//

#import <Cocoa/Cocoa.h>
#import <QuartzCore/QuartzCore.h>
#import <Accelerate/Accelerate.h>

/**
 *   Fast Fourier Hilfsmethoden
 */
 
// Hintransformation
int fft_forward(DSPComplex * img, int width, int height){
	return fft(img, width, height, FFT_FORWARD);
}

// Rücktransformation
int fft_backward(DSPComplex * img, int width, int height){
	return fft(img, width, height, FFT_INVERSE);
}

// Allgemeine FFT
int fft(DSPComplex * data, int width, int height, FFTDirection dir){
	
	/* Total area or volume covered by the dimensions */
	int dims = width*height;
	
	/* Get log (base 2) of each dimension */
	int logX = log2((double)width);
	int logY = log2((double)height);

	
	/* Setup and allocate some memory for performing the FFTs */
	DSPSplitComplex input;
	
	input.realp = (float*)malloc(dims*sizeof(float));
	input.imagp = (float*)malloc(dims*sizeof(float));
	
	/* Setup the weights (twiddle factors) for performing the FFT */
	FFTSetup fft_weights;
	fft_weights = create_fftsetup(MIN(logX,logY),kFFTRadix2);
	
	/* Split the complex (interleaved) data into two arrays */
	ctoz(data,2,&input,1,dims);
	
	/* Perform the FFT */
	fft2d_zip(fft_weights,&input,1,0,logX,logY,dir);
	/* Merge the separate data (real and imaginary) parts into a single array */
	ztoc(&input,1,data,2,dims);
	
	//Normalize energy 
	if(dir == FFT_INVERSE){
		int i;
		for(i=0;i<dims;i++){
			data[i].real=data[i].real/dims;
			data[i].imag=data[i].imag/dims;
		}
	}
	/* Free up the allocations */
	free(input.realp);
	free(input.imagp);
	
	destroy_fftsetup(fft_weights);
	
	return 0;
}

// Bewegt den Ursprung des Frequenzraumes in die Bildmitte
void moveDCToCenter(float* image, int width, int height)
{
	
    int w1 = width/2;
    int h1 = height/2;
    int w2 = (width+1)/2;
    int h2 = (height+1)/2;

	float temp ;
	
	int x,y;
	
    // 2. Quadrant  zum 4. und umgekehrt
	for(y=0; y<h2; ++y){
		for(x=0; x<w2; ++x){
			temp= image[(x+w1)+(y+h1)*width];
			image[(x+w1)+(y+h1)*width] = image[x+y*width];
			image[x+y*width] = temp;
		}
	}	
	
    // 1. Quadrant  zum 3. und umgekehrt
	for(y=0; y<h2; ++y){
		for(x=w2; x<width; ++x){
			temp= image[(x-w2)+(y+h1)*width];
			image[(x-w2)+(y+h1)*width] = image[x+y*width];
			image[x+y*width] = temp;
		}
	}
}

// Berechnet das Amplitudenspektrum eines Bildes
int image2spectralAmplitude(float * img, int width, int height){
	int i;
	int dims = width*height;
	
	DSPComplex *data = (DSPComplex*)malloc(dims*sizeof(DSPComplex));
	
	for(i=0; i<dims;++i){
		data[i].real = img[i];
	}
	
	int res = fft_forward(data, width, height);
	
	if(res==0){
		for(i=0; i<dims;++i){
			img[i]= sqrt(data[i].real*data[i].real + data[i].imag*data[i].imag) ;
		}
	}
	return res;
}

// Berechnet das Amplitudenspektrum eines vImage_Buffers
int vImage_Buffer2spectralAmplitude(vImage_Buffer * buf){
	return image2spectralAmplitude( buf->data, buf->width, buf->height);
}



/**
 *   Debug Hilfen: Geben das Bild (formatiert) auf der Konsole aus
 */
 
//Stellt Amplitudenspektrum eines Komplexen Bildes dar
void displayComplexImageMagnitude(DSPComplex * data, int width, int height){
	int i,j;
	for(j=0;j<height;j++){
		for(i=0;i<width;i++){
			printf("%2.02f \t",sqrt(data[i+j*width].real*data[i+j*width].real + data[i+j*width].imag*data[i+j*width].imag) );
		}
		printf("\n");
	}
}

//Stellt einfach ein Bild dar
void displayImage(float * data, int width, int height){
	int i,j;
	for(j=0;j<height;j++){
		for(i=0;i<width;i++){
			printf("%05.02f \t",data[i+j*width]);
		}
		printf("\n");
	}
}

//Stellt ein vImage_Buffer dar
void displayvImage_Buffer(vImage_Buffer * buf){
	displayImage( buf->data, buf->width, buf->height);
}


//Main-Methode
int main(int argc, char *argv[])
{	

	//load an image
	NSImage *loaded_image = [[NSImage alloc] initWithContentsOfFile:[NSString stringWithCString:"/Users/seppke/Desktop/rect.jpg"]];
	
	//determine image size
	NSSize originalSize = [loaded_image size];
	int w = originalSize.width;
	int h = originalSize.height;
	
	//determine padding width for fft padding
	int padded_w = 1<< (int)(log2((double)w)+0.5);
	int padded_h = 1<< (int)(log2((double)h)+0.5);
	padded_w = MAX(padded_w, padded_h);
	padded_h = MAX(padded_w, padded_h);
	
	//get Bitmap-Representation for loaded image
	printf("Getting reps of non- and padded image\n");
	NSData *tiffRep = [loaded_image TIFFRepresentation];
	NSBitmapImageRep *loaded_rep = [[NSBitmapImageRep alloc] initWithData:tiffRep];
	
	//get Bitmap-Representation for padded image
	NSBitmapImageRep* new_rep = [[NSBitmapImageRep alloc]	initWithBitmapDataPlanes:NULL
															pixelsWide:padded_w
															pixelsHigh:padded_h
															bitsPerSample:[loaded_rep bitsPerSample]
															samplesPerPixel:[loaded_rep samplesPerPixel]
															hasAlpha:[loaded_rep hasAlpha]
															isPlanar:[loaded_rep isPlanar]
															colorSpaceName:[loaded_rep colorSpaceName]
															bitmapFormat:[loaded_rep bitmapFormat]
															bytesPerRow:0/*[loaded_rep bytesPerRow]*/
															bitsPerPixel:0/*[loaded_rep bitsPerPixel]*/];
	//create pointer for padded image (float for each pixel) 
	float* padded_img = (float*)malloc(padded_w*padded_h*sizeof(float));
	
	//create pointer for padded image of type byte, with meaning AARRGGBB or else 
	//warning: takes [loaded_rep bitsPerSample]*[loaded_rep samplesPerPixel] bits per pixel!
    unsigned char *loaded_data;
    unsigned char *new_data;

    if ([loaded_rep bitsPerSample] != 8) {
        // usually each color component is stored
        // using 8 bits, but recently 16 bits is occasionally
        // also used ... you might want to handle that here
        printf("bitpersample!=8.... was %d\n",[loaded_rep bitsPerSample]);
		return 1;
    }
    
	//Assign to pointers
    loaded_data = [loaded_rep bitmapData];
    new_data = [new_rep bitmapData];
	
    w = [loaded_rep pixelsWide];
    h = [loaded_rep pixelsHigh];
	
    int x, y;
	float *pf = padded_img;
	//run trough loaded_data and fill new_data and padded_image on the fly
    for (y=0; y<h; y++) {
        unsigned char *loaded_p = loaded_data + [loaded_rep bytesPerRow]*y;
        unsigned char *new_p = new_data + [new_rep bytesPerRow]*y;
        for (x=0; x<w; x++) {
            // do your thing
            // p[0] is red, p[1] is green, p[2] is blue, p[3] can be alpha if spp==4
			new_p[0] = loaded_p[0];
			new_p[1] = loaded_p[1];
			new_p[2] = loaded_p[2];
			
			*pf = (float)loaded_p[0];
			//printf("setting I(%4d, %4d) =\t %05f\n",x,y,*pf);
			
			pf++;
			
            loaded_p += [loaded_rep samplesPerPixel];
			new_p += [new_rep samplesPerPixel];
        }
    }

	//Calculate the spectral amplitude
	image2spectralAmplitude(padded_img,padded_w, padded_h);
	
	//move dc to center
	moveDCToCenter(padded_img, padded_w, padded_h);
	
	//write result back to new_images data
	pf = padded_img;
	for (y=0; y<padded_h; y++) {
        unsigned char *p = new_data + [new_rep bytesPerRow]*y;
        for (x=0; x<padded_w; x++) {
            // do your thing
            // p[0] is red, p[1] is green, p[2] is blue, p[3] can be alpha if spp==4
			
			p[0] = (unsigned char) MIN(sqrt(*pf),255);
			p[1] = (unsigned char) MIN(sqrt(*pf),255);
			p[2] = (unsigned char) MIN(sqrt(*pf),255);
			//printf("setting I(%4d, %4d) =\t %05f\n",x,y,*pf);
			pf++;
			
            p += [new_rep samplesPerPixel];
        }
    }
	//write to file
	NSData *bitmapData = [new_rep representationUsingType:NSJPEGFileType properties:NULL];
	[bitmapData writeToFile:[NSString stringWithCString:"/Users/seppke/Desktop/rect2.jpg"] atomically:YES];

	//clean up
	free(padded_img);
	return NSApplicationMain(argc,  (const char **) argv);
}