#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "nrutil.h"
#include "binaryMatrixFile.h"

//Note that the first component must be set to 1.

void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, int isign);

/* example1 This fragment shows how one might filter a 256 by 256 digital image. */
int main(void) {
	float ***data, **speq;

	//Here the image would be loaded into data.	
	double *bmfData = NULL;
	int nX=0, nY=0;
	if(bmfLoad("221-12-512x512.bin", &nX, &nY, &bmfData, 0)){ printf("Read failed\n"); return -1; }
	printf("Loaded %ix%i\n", nX, nY); fflush(stdout);

	data = f3tensor(1,1,1,nX,1,nY);
	speq = matrix(1,1,1,2*nX);
	printf("init done\n"); fflush(stdout);

	rlft3(data, speq, 1, nX, nY, 1);
	printf("fft done\n"); fflush(stdout);

	// Here the arrays data and speq would be multiplied by a suitable filter function (of frequency).
	rlft3(data, speq,1,nX,nY,-1);
	printf("ifft done\n"); fflush(stdout);

	//Here the filtered image would be unloaded from data.
	if(bmfWrite("filtered.bin", nX, nY, NULL, NULL, bmfData, 0)){ printf("Write failed\n"); return -1; }

	free_matrix(speq,1,1,1,2*nX);
	free_f3tensor(data,1,1,1,nX,1,nY);
printf("free done\n"); fflush(stdout);

	return 0;
}