/* computes 2d fft's and normalizes windows for mitosm.cpp */

#include <math.h>
#include <alloc.h>
#include <windows.h>
#include <owl.h>

void fourn(double far *data, int *nn, int ndim, int isign);

void fftnorm(unsigned char far **image1, unsigned char far **image2, int width, int height)
{
	int i, j, k, *dims;
  double mag1, mag2, mag_scale;
	double far *im1, *im2;
	HGLOBAL hglb1, hglb2;
  double huge *pt1, huge *pt2;

	hglb1 = GlobalAlloc(GMEM_FIXED, (long) 2*width*height*sizeof(double));
	im1 = (double far *) GlobalLock(hglb1);
	hglb2 = GlobalAlloc(GMEM_FIXED, (long) 2*width*height*sizeof(double));
	im2 = (double far *) GlobalLock(hglb2);

	dims = (int *) malloc (2 * sizeof(int));
	dims[0] = width;
	dims[1] = height;

	pt1 = (double huge *) im1;
	pt2 = (double huge *) im2;
	for (i=0 ; i<height ; i++)
		for (j=0 ; j<width ; j++)
		{
			*pt1++ = (double) image1[i][j];
			*pt1++ = 0.0;
			*pt2++ = (double) image2[i][j];
			*pt2++ = 0.0;
		}

	/* 2D fourier transforms */
	fourn(im1-1, dims-1, 2, 1);
	fourn(im2-1, dims-1, 2, 1);

  /* adjust magnitude of image2 */
	pt1 = im1;
  pt2 = im2;
	for (i=0 ; i<width*height ; i++)
	{
		mag1 = (*pt1) * (*pt1);
		pt1++;
		mag1 += (*pt1) * (*pt1);
		pt1++;
		mag2 = (*pt2) * (*pt2);
		pt2++;
		mag2 += (*pt2) * (*pt2);
    if (mag2 == 0.0) mag2 = 1.0;
		pt2--;
		mag_scale = sqrt(mag1/mag2);
		*pt2++ *= mag_scale;
		*pt2++ *= mag_scale;
	}

	/* inverse transform image2 */
	fourn(im2-1, dims-1, 2, -1);

	pt2 = im2;
	for (i=0 ; i<height ; i++)
		for (j=0 ; j<width ; j++)
		{
			*pt2 /= ((double) width*height);
			if (*pt2 > 255.0) *pt2 = 255.0;
      if (*pt2 < 0.0) *pt2 = 0.0;
			image2[i][j] = (unsigned char) *pt2++;
			*pt2++;
		}

	free(dims);
	GlobalUnlock(hglb1);
	GlobalFree(hglb1);
	GlobalUnlock(hglb2);
	GlobalFree(hglb2);
}