/*********************************************************************/
/* /vision/usr/barron/2dfft.c					     */
/* This program provides tools for performing 2D Fourier Transforms  */
/* on square images that are powers of 2 and doing lowpass, highpass */
/* bandpass and bandreject filtering. A histogram equalization       */
/* routine, routines to do phase and magnitude suppression, etc. are */
/* also included. Note that pixel values are not rescaled to be      */
/* between 0 and 255 before the images are saved except in equalize  */
/* where the function rescale_pixels is used to do this.	     */
/*								     */
/*       Written by John Barron, December 1989.			     */
/*								     */
/*********************************************************************/
#include <math.h>
#include <stdio.h>
#include <fcntl.h>
#include <rasterfile.h>

#define maxn 256
#define epsilon 0.001
#define image_header_size 32
#define numrows 512    /* Image size */     
#define numcols 512    /* Work size  */    
#define maxrows 512
#define maxcols 512
#define maxlutsize 256
#define TRUE 1
#define FALSE 0
#define WHITE 255
#define BLACK 0
#define NUMBER 999999.9


float pi = M_PI;
float twopi = 6.2831853;

float ffttransform = -1.0;
float fftinverse   = 1.0;

float rasterfile[maxrows][maxcols][2];
float pic1[maxrows][maxcols][2],pic2[maxrows][maxcols][2];
float pic3[maxrows][maxcols][2];
float energyfile[maxrows][maxcols][2];
float edges[maxrows][maxcols][2];

int sign();
void NLOGN();
void fft_2d(),add(),sub();
void gauss_smoothing(),transpose(),mult(),laplacian();
void equalize(),tom_high_pass();
float sinc();


/*=========================================================================*/
/*  Main Routine                                                           */
/*=========================================================================*/


main()
{
char header[32];

loadcolour("630clip09.ras",pic1,pic2,pic3);
savecolour("630clip09-1.ras",pic1,pic2,pic3);
if(TRUE) exit(1);

loadraster ("630image",rasterfile,header);
translate_to_origin( rasterfile );
fft_2d (8,rasterfile, ffttransform );
butterworth_highpass(rasterfile);
fft_2d (8,rasterfile, fftinverse );
translate_to_origin( rasterfile );
rescale_pixels(rasterfile,8); 
saveraster ("630out",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    compute_energy_spectrum(rasterfile,energyfile);
    saveraster ("630energy",energyfile,header);
if(TRUE) exit(1);



loadraster("630image",pic1,header);
translate_to_origin(pic1);
fft_2d(8,pic1,ffttransform);
copy_data(pic1,pic2);
copy_data(pic1,pic3);
butterworth_lowpass(pic1,0.05,1);
butterworth_bandpass(pic2,0.05,1,0.1);
butterworth_highpass(pic3,0.05,1);
fft_2d(8,pic1,fftinverse);
fft_2d(8,pic2,fftinverse);
fft_2d(8,pic3,fftinverse);
translate_to_origin(pic1);
translate_to_origin(pic2);
translate_to_origin(pic3);
rescale_pixels(pic1,8); 
rescale_pixels(pic2,8); 
rescale_pixels(pic3,8); 
savecolour("630colour",pic1,pic2,pic3);
if(TRUE) exit(1);

test0();
if(TRUE) exit(1);




if(TRUE) exit(1);
    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile ); 
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.0, 1.0);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile ); 
    saveraster ("630smooth",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile ); 
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.15,1);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile ); 
    saveraster ("630butterhighpass0.15",rasterfile,header);
if(TRUE) exit(1);

test0();

loadraster ("630image",rasterfile,header);
set_ideal_lowpass(rasterfile,0.05);
fft_2d(8,rasterfile,fftinverse);
rearrange(rasterfile,8);
rescale_pixels(rasterfile,8);
saveraster ("ideal",rasterfile,header);
if(TRUE) exit(1);

/*  test1();
    test2();
    test3(); 
    test4();
    if(TRUE) exit(1);  */

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.0, 1.0);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630smooth2.0",rasterfile,header);
if(TRUE)exit(1);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.0, 1.0); 
    laplacian( rasterfile );
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    compute_zero_crossings(rasterfile,edges,8);
    rescale_pixels(rasterfile,8);
    saveraster ("630laplacian2.0",rasterfile,header);
    saveraster ("630edges2.0",edges,header);
    printf("630laplacian2.0 saved\n");
if(TRUE) exit(1);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.25,1);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630butterhighpass0.25-1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.25,5);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630butterhighpass0.25-5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.40,1);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630butterhighpass0.4-1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.40,5);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630butterhighpass0.4-5",rasterfile,header);
    if(TRUE) exit(1);


    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.0, 1.0); 
    laplacian( rasterfile );
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    compute_zero_crossings(rasterfile,edges,8);
    rescale_pixels(rasterfile,8);
    saveraster ("630laplacian2.0",rasterfile,header);
    saveraster ("630edges2.0",edges,header);
    printf("630laplacian2.0 saved\n");
if(TRUE) exit(1);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.0, 1.0);
    fft_2d (8,rasterfile, fftinverse );
    translate_to_origin( rasterfile );
    saveraster ("630smooth2.0",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    translate_to_origin( rasterfile );
    fft_2d (8,rasterfile, ffttransform );
    compute_energy_spectrum(rasterfile,energyfile);
    saveraster ("630energy",energyfile,header);
if(TRUE) exit(1);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.5, 1.0);
    fft_2d (8,rasterfile, fftinverse);
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.5, -1.0);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630original",rasterfile,header);
    system("print_postscript 630original");
    if(TRUE) exit(1);

    loadraster ("630blurred2.5",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.5, -1.0);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630unblurred2.5",rasterfile,header);

    loadraster ("630smooth2.5",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    gauss_smoothing ( rasterfile, 2.5, -1.0);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630unsmooth2.5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    equalize(rasterfile,8);
    saveraster ("630equal",rasterfile,header);

    loadraster ("630blurred2.5",rasterfile,header);
    equalize(rasterfile,8);
    saveraster ("630equal2.5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_mag(rasterfile,8,0.1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630mag0.1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_mag(rasterfile,8,0.5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630mag0.5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_mag(rasterfile,8,0.7);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630mag0.7",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_mag(rasterfile,8,0.9);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630mag0.9",rasterfile,header);


    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_phase(rasterfile,8,0.1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630phase0.1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_phase(rasterfile,8,0.5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630phase0.5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_phase(rasterfile,8,0.7);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630phase0.7",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    scale_phase(rasterfile,8,0.9);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630phase0.9",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    compute_energy_spectrum(rasterfile,energyfile);
    saveraster ("630energy",energyfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    duda_hart_low_pass(rasterfile);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630lowpass",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    duda_hart_high_pass(rasterfile);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630highpass",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    duda_hart_high_emphasis(rasterfile);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630highemphasis",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_lowpass(rasterfile,0.15,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterlowpass0.15",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_lowpass(rasterfile,0.25,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterlowpass0.25",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_lowpass(rasterfile,0.4,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterlowpass0.4",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.25,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.25-1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.25,5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.25-5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.40,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.4-1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.40,5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.4-5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandreject(rasterfile,0.25,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandreject0.25",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandreject(rasterfile,0.1,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandreject0.1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandreject(rasterfile,0.4,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandreject0.4",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.25,5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.25-5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.15,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.15",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.15,3);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.15-3",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_highpass(rasterfile,0.15,5);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterhighpass0.15-5",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandpass(rasterfile,0.25,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandpass0.25",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandpass(rasterfile,0.1,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandpass0.1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    butterworth_bandpass(rasterfile,0.4,1,0.25);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630butterbandpass0.4",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,1.0,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc1",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,1.0,2);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc2",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,1.0,3);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc3",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,1.0,4);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc4",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,2.0,1);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc1amp2",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,2.0,2);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc2amp2",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,2.0,3);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc3amp2",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,2.0,4);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630sinc4amp2",rasterfile,header);

    loadraster ("630image",rasterfile,header);
    fft_2d (8,rasterfile, ffttransform );
    sinc_filter(rasterfile,4.0,4);
    fft_2d (8,rasterfile, fftinverse );
    equalize(rasterfile,8);
    saveraster ("630sinc4amp4equal",rasterfile,header);

    loadraster ("630image1",rasterfile,header);
    loadraster ("630image2",energyfile,header);
    fft_2d (8,rasterfile, ffttransform );
    fft_2d (8,energyfile, ffttransform );
    average_phase_and_mag(rasterfile,energyfile,rasterfile,8);
    fft_2d (8,rasterfile, fftinverse );
    saveraster ("630average",rasterfile,header);
}   /* main */

/*=========================================================================*/


/****************************************************************/
/* This function transposes the complex image			*/
/****************************************************************/
void transpose ( matrix )
    float matrix [maxrows][maxcols][2];

{
    int i, j;
    float temp;

    for (i=0; i<maxrows; i++)
        for (j=0; j<i; j++)  {
            temp = matrix [i][j][0];
            matrix [i][j][0] = matrix [j][i][0];
            matrix [j][i][0] = temp;

            temp = matrix [i][j][1];
            matrix [i][j][1] = matrix [j][i][1];
            matrix [j][i][1] = temp;
        }

}   /* transpose */


/*=========================================================================*/
/*=========================================================================*/


/****************************************************************/
/* This function computes the 2D Fourier Transform using the    */
/* 1D FFT allpied to rows and the columns if direct=-1.0 If     */
/* direct=1.0 then the invserse 2D Fourier Transform is         */
/* computed by applying the 1D inverse FFT to the columns and   */
/* then the rows of the complex image. direct must be float.    */
/****************************************************************/
void fft_2d(image_dim,matrix,direct)

    float matrix [maxrows][maxcols][2];
    float direct;
    int image_dim;

{
    int i,dim;
    double t1,t2;
    float direction;

    direction = direct;
    t1 = 2.0;
    t2 = image_dim*1.0;
    dim = pow(t1,t2);
    if (direction==ffttransform)  {
        printf ( "Fast fourier transform\n" );
	fflush(stdout);

        for (i=0; i<dim; i++)
	    {
            NLOGN(image_dim,&matrix[i][0][0],direction );
	    }

	printf("Processed %d rows\n",dim);
	fflush(stdout);

        transpose ( matrix );

        for (i=0; i<dim; i++)
	    {
            NLOGN(image_dim,&matrix[i][0][0], direction );
	    }
	printf("Processed %d columns\n",dim);
	fflush(stdout);
        transpose ( matrix );
        /* rearrange(matrix,image_dim); */
    }
    else  {
        printf ( "Inverse fast fourier transform\n" );
	fflush(stdout);
        /* rearrange(matrix,image_dim);*/
        transpose ( matrix );

        for (i=0; i<dim; i++)
	    {
            NLOGN(image_dim,&matrix[i][0][0], direction );
	    }
	printf("Processed %d columns\n",dim);
	fflush(stdout);

        transpose ( matrix );

        for (i=0; i<maxrows; i++)
	    {
            NLOGN(image_dim,&matrix[i][0][0],direction );
	    }
	printf("Processed %d rows\n",dim);
	fflush(stdout);
    }
    fflush(stdout);

}   /* fft_2d */


/****************************************************************/
/* This function blurrs the image by using the FT of the        */
/* Gaussian function. If direction is 1.0 positive blurring     */
/* occurs but if direction is -1.0 negative blurring, i.e.      */
/* deblurring, is attempted.					*/
/****************************************************************/
void gauss_smoothing ( image, sigma, direction )

    float image [maxrows][maxcols][2];
    float sigma;
    float direction;

{
    int u,v,uu,vv,i,j;
    float halfsigma;
    float G[maxrows][maxcols],constant;


    printf ( "Guassian smoothing direction=%f\n",direction);
    halfsigma = -0.5 * sigma * sigma * direction * 4.0 * pi * pi; 
    constant = maxrows*maxrows; 

    for (u=0; u<maxrows; u++)
    for (v=0; v<maxcols; v++)
	{
	uu = u-(maxrows/2);
	vv = v-(maxcols/2);
	G[u][v] = exp(((uu*uu+vv*vv)/constant)*halfsigma);
	}
     for (u=0; u<maxrows; u++)
     for (v=0; v<maxcols; v++)  
        {
        image[u][v][0] = image[u][v][0] * G[u][v];
        image[u][v][1] = image[u][v][1] * G[u][v];
        }

}   /* gauss_smoothing */


/****************************************************************/
/* This function reads in the image in the rasterfile specified */
/* by filename and sets the image array called rasterfile to    */
/* hold the image in its real part. header is the 32 byte 	*/
/* header information for the rasterfile.			*/
/****************************************************************/
loadraster(filename,rasterfile,header)
char filename[100];
float rasterfile[maxrows][maxcols][2];
char header[32];
{
unsigned char pic[maxrows][maxcols];
int fd,i,j;

fd = open(filename,O_RDONLY,0755);
read(fd,header,32);
read(fd,pic,sizeof(pic));
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	rasterfile[i][j][0] = pic[i][j];
	rasterfile[i][j][1] = 0;
	}
close(fd);
}


/****************************************************************/
/* This function writes the real part of an image in rasterfile */
/* into a file specified by filename with header to make it a   */
/* rasterfile.							*/
/****************************************************************/
saveraster(filename,rasterfile,header)
char filename[100];
float rasterfile[maxrows][maxcols][2];
char header[32];
{
unsigned char pic[maxrows][maxcols],pic2[maxrows][maxcols];
float mag,max_mag,min_mag,sum_mag; /* Use to compute imaginary magnitudes */
float sum_mag2,min_mag2,max_mag2,mag2;
int fd,i,j;

sum_mag = 0.0;
min_mag =  HUGE;
max_mag = -HUGE;
sum_mag2 = 0.0;
min_mag2 =  HUGE;
max_mag2 = -HUGE;
fd = creat(filename,0644);
printf("File %s created\n",filename);
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	pic[i][j] = (unsigned) rasterfile[i][j][0];
	pic2[i][j] = (unsigned) rasterfile[i][j][1];
	mag = rasterfile[i][j][1];
	mag2 = rasterfile[i][j][0];
	sum_mag += mag;
	sum_mag2 += mag2;
	if(mag < min_mag) min_mag = mag;
	if(mag > max_mag) max_mag = mag;
	if(mag2 < min_mag2) min_mag2 = mag2;
	if(mag2 > max_mag2) max_mag2 = mag2;
	}
printf("Average Real residual: %f\n",sum_mag2/(maxrows*maxcols*1.0));
printf("Minimum Real individual residual: %f\n",min_mag2);
printf("Maximum Real individual residual: %f\n",max_mag2);
printf("Average Imaginary residual: %f\n",sum_mag/(maxrows*maxcols*1.0));
printf("Minimum Imaginary individual residual: %f\n",min_mag);
printf("Maximum Imaginary individual residual: %f\n",max_mag);
write(fd,header,32);
write(fd,pic,sizeof(pic));
close(fd);
fd = creat("imaginary_image",0644);
write(fd,header,32);
write(fd,pic2,sizeof(pic2));
close(fd);
}

/****************************************************************/
/* This function reads the real part of 3 complex image as     */
/* 24-bit colour image  					*/
/* The input image has to be maxrows by maxcols                 */
/****************************************************************/
loadcolour(filename,pic1,pic2,pic3)
char filename[100];
float pic1[maxrows][maxcols][2];
float pic2[maxrows][maxcols][2];
float pic3[maxrows][maxcols][2];
{
int header[8];
unsigned char data[maxrows][maxcols][3];
unsigned char image1[maxrows][maxcols];
unsigned char image2[maxrows][maxcols];
unsigned char image3[maxrows][maxcols];
int fd,i,j;

fd = open(filename,O_RDONLY);
printf("File %s opened\n",filename);
read(fd,header,32);
read(fd,data,maxrows*maxcols*3);
/* Store blue, green and red planes as complex floats */
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	pic1[i][j][0] = (float) ((int) (data[i][j][0]));
	pic2[i][j][0] = (float) ((int) (data[i][j][1]));
	pic3[i][j][0] = (float) ((int) (data[i][j][2]));
	pic1[i][j][1] = 0.0;
	pic2[i][j][1] = 0.0;
	pic3[i][j][1] = 0.0;
	}
printf("File %s read\n",filename);
}

/****************************************************************/
/* This function writes the real part of 3 complex image as     */
/* 24-bit colour image  					*/
/****************************************************************/
savecolour(filename,pic1,pic2,pic3)
char filename[100];
float pic1[maxrows][maxcols][2];
float pic2[maxrows][maxcols][2];
float pic3[maxrows][maxcols][2];
{
int header[8];
unsigned char data[maxrows][maxcols][3];
unsigned char image1[maxrows][maxcols];
unsigned char image2[maxrows][maxcols];
unsigned char image3[maxrows][maxcols];
int fd,i,j;

fd = creat(filename,0644);
printf("File %s created\n",filename);
/* Convert float reals to unsigned char, ignore complex part */
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	image1[i][j] = (unsigned char) ((int) pic1[i][j][0]);
	image2[i][j] = (unsigned char) ((int) pic2[i][j][0]);
	image3[i][j] = (unsigned char) ((int) pic3[i][j][0]);
	}
header[0] = RAS_MAGIC;
header[1] = maxcols;
header[2] = maxrows;
header[3] = 24;
header[4] = maxrows*maxcols*3;
header[5] = RT_STANDARD;
header[6] = RMT_NONE;
header[7] = 0;
write(fd,header,32);
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	data[i][j][0] = image1[i][j];
	data[i][j][1] = image2[i][j];
	data[i][j][2] = image3[i][j];
	}
write(fd,data,maxrows*maxcols*3);
close(fd);
printf("File %s written\n",filename);
}

/****************************************************************/
/* This function computed a normalized erergy map of the        */
/* complex image stored in rasterfile and saves that image in   */
/* energyfile. Energy is computed as the log of 1+magnitude of  */
/* complex FFT results and then normalized to be between 0 and  */
/* 255.								*/
/****************************************************************/
compute_energy_spectrum(rasterfile,energyfile)
float rasterfile[maxrows][maxcols][2];
float energyfile[maxrows][maxcols][2];
{
int i,j;
float min,max,diff;
int i_min,i_max,j_min,j_max;

min = HUGE;
max = 0.0;
i_min = i_max = j_min = j_max = 0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	energyfile[i][j][0] = log(1.0+rasterfile[i][j][0]*rasterfile[i][j][0]
				     +rasterfile[i][j][1]*rasterfile[i][j][1]);
	energyfile[i][j][1] = 0.0;
	if(min > energyfile[i][j][0]){ min = energyfile[i][j][0]; i_min=i; j_min=j; }
	if(max < energyfile[i][j][0]){ max = energyfile[i][j][0]; i_max=i; j_max=j; }
	}
printf("Minimum energy is %f at %d,%d\n",min,i_min,j_min);
printf("Maximum energy is %f at %d,%d\n",max,i_max,j_max);
/* scale the energy image */
diff = max-min;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	energyfile[i][j][0] = ((energyfile[i][j][0]-min)/diff)*255.0;
	}
}

/****************************************************************/
/* This function computes the maximum real and imaginary parts  */
/* of some image.						*/
/****************************************************************/
compute_max_real_imaginary(string,rasterfile,maxreal,maximag)
char string[100];
float rasterfile[maxrows][maxcols][2];
float *maxreal,*maximag;
{
int i,j;
float max1,max2;
max1 = 0.0;
max2 = 0.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	if(fabs(rasterfile[i][j][0])>max1) max1 = fabs(rasterfile[i][j][0]);
	if(fabs(rasterfile[i][j][1])>max2) max2 = fabs(rasterfile[i][j][1]);
	}
*maxreal = max1;
*maximag = max2;
printf("\n%s\n",string);
printf("\nMaximum abs of real: %f \nMaximum abs of imaginary: %f\n",
	    *maxreal,*maximag);
fflush(stdout);
}

/****************************************************************/
/* This function takes some complex image and scales the real   */
/* part (z=0) or the imaginary part (z=1) so that all numbers   */
/* are between 0 and 255: in this way the data can be displayed */
/* as an image.							*/
/****************************************************************/
make_image(filename,rasterfile,header,z)
char filename[100],header[32];
float rasterfile[maxrows][maxcols][2];
int z;
{
int fd,i,j;
unsigned char pic[maxrows][maxcols];
float min,max,diff;

min = 99999999.9;
max = 0.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	if(fabs(rasterfile[i][j][z]) > max) max = fabs(rasterfile[i][j][z]);
	if(fabs(rasterfile[i][j][z]) < min) min = fabs(rasterfile[i][j][z]);
	}
diff = max-min;

fd = creat(filename,0644);
write(fd,header,32);
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	pic[i][j] = (fabs(rasterfile[i][j][z])-min)/diff*255;
write(fd,pic,sizeof(pic));
close(fd);
}

/****************************************************************/
/* This function performs low pass filtering as specified in    */
/* Duda and Hart, Pattern Classification and Scene Analysis,    */
/* John Wiley and Sons, 1973, pp313.				*/
/****************************************************************/
duda_hart_low_pass(rasterfile)
float rasterfile[maxrows][maxcols][2];
{
float H[maxrows][maxcols];
int i,j;
float u,v;
double sixteen;

sixteen = 16.0;

for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
H[i][j] = pow((cos(pi*u)*cos(pi*v)),sixteen);
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs high pass filtering as specified in   */
/* Duda and Hart, Pattern Classification and Scene Analysis,    */
/* John Wiley and Sons, 1973, pp314.				*/
/****************************************************************/
duda_hart_high_pass(rasterfile)
float rasterfile[maxrows][maxcols][2];
{
float H[maxrows][maxcols];
int i,j;
float u,v;
double four;

four = 4.0;

for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
H[i][j] = 1.5-pow((cos(pi*u)*cos(pi*v)),four);
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs high emphasis frequency filtering as  */
/* specified in Duda and Hart, Pattern Classification and Scene */
/* Analysis, John Wiley and Sons, 1973, pp314.			*/
/****************************************************************/
duda_hart_high_emphasis(rasterfile)
float rasterfile[maxrows][maxcols][2];
{
float H[maxrows][maxcols];
int i,j;
float u,v;
double four;

four = 4.0;

for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
H[i][j] = 2.0-pow((cos(pi*u)*cos(pi*v)),four);
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}
/****************************************************************/
/* This function performs low pass Butterworth filtering as     */
/* specified in Gonzalez and Wintz, Digital Image Processing,   */
/* 1977, pp145.							*/
/****************************************************************/
butterworth_lowpass(rasterfile,D0,n)
float rasterfile[maxrows][maxcols][2],D0;
int n;
{
float H[maxrows][maxcols];
double two_n,Duv;
float u,v;
int i,j;

two_n = 2.0*n;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
H[i][j] = 1.0/(1.0+pow((Duv/D0),two_n));
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs high pass Butterworth filtering as    */
/* specified in Gonzalez and Wintz, Digital Image Processing,   */
/* 1977, pp162.							*/
/* Mov, 1994 - highpass filtering removes all dc component      */
/* either rescale (you'll get both negative and positive        */
/* numbers) or let some very low frequency information pass     */
/****************************************************************/
butterworth_highpass(rasterfile,D0,n)
float rasterfile[maxrows][maxcols][2],D0;
int n;
{
float H[maxrows][maxcols];
double two_n,Duv;
float u,v;
int i,j;


two_n = 2.0*n;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
if(Duv != 0.0) H[i][j] = 1.0/(1.0+pow((D0/Duv),two_n));
else H[i][j] = 0.0;
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs band reject Butterworth filtering as  */
/* specified in Gonzalez and Wintz, Digital Image Processing,   */
/* 1977, pp181.							*/
/****************************************************************/
butterworth_bandreject(rasterfile,D0,n,W)
float rasterfile[maxrows][maxcols][2],D0,W;
int n;
{
float H[maxrows][maxcols],Duv;
double two_n;
float u,v;
int i,j;


two_n = 2.0*n;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
if(Duv!=D0) H[i][j] = 1.0/(1.0+pow((Duv*W)/(Duv*Duv-D0*D0),two_n));
else H[i][j] = HUGE;
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs bandpass Butterworth filtering as     */
/* specified in Gonzalez and Wintz, Digital Image Processing,   */
/* 1977, pp180.							*/
/****************************************************************/
butterworth_bandpass(rasterfile,D0,n,W)
float rasterfile[maxrows][maxcols][2],D0,W;
int n;
{
float H[maxrows][maxcols],Duv;
double two_n;
float u,v;
int i,j;


two_n = 2.0*n;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
if(Duv!=D0) H[i][j] = -(1.0/(1.0+pow((Duv*W)/(Duv*Duv-D0*D0),two_n))-1.0);
else H[i][j] = HUGE;
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/****************************************************************/
/* This function performs ideal low pas filtering		*/
/****************************************************************/
ideal_lowpass(rasterfile,D0)
float rasterfile[maxrows][maxcols][2],D0;
{
float H[maxrows][maxcols];
double two_n,Duv;
float u,v;
int i,j;


for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
if(Duv <= D0) H[i][j] = 1.0;
else H[i][j] = 0.0;
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}
/**********************************************/
/* Set the image to the ideal low pass filter */
/**********************************************/
set_ideal_lowpass(rasterfile,D0)
float rasterfile[maxrows][maxcols][2],D0;
{
float H[maxrows][maxcols];
double two_n,Duv;
float u,v;
int i,j,ct1,ct2;


ct1 = ct2 = 0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
Duv = sqrt(u*u+v*v);
if(Duv <= D0) { H[i][j] = 100.0; ct1++; }
else { H[i][j] = 0.0; ct2++; }
}
printf("Ideal filter coefficients computed\n %d non-zero coefficients\n %d zero coefficients\n",ct1,ct2);

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = H[i][j];
	rasterfile[i][j][1] = 0.0;
	}
}

/****************************************************************/
/* This function performs low pass filtering using a sinc       */
/* function  							*/
/****************************************************************/
sinc_filter(rasterfile,A,n)
float rasterfile[maxrows][maxcols][2],A;
int n;
{
float H[maxrows][maxcols];
float u,v;
float min,max,diff;
int i,j;


for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = (i*1.0/maxrows-0.5)*pi*n;
v = (j*1.0/maxcols-0.5)*pi*n;
H[i][j] = A*sinc(u)*sinc(v);
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}

/*****************************************************************/
/* Rescale image so between 0 and 255				 */
/*****************************************************************/
rescale_pixels(rasterfile,dim)
float rasterfile[maxrows][maxcols][2];
int dim;
{
float min,max,diff;
int i,j,size;
int mini,minj,maxi,maxj;

size = pow(2.0,dim*1.0);
printf("\nIn rescale_pixels: size=%d\n",size);
/* Make all pixel values are between 0 and 255. */
min = rasterfile[0][0][0];
max = rasterfile[0][0][0];
mini=maxi=minj=maxj=0;
for(i=0;i<size;i++)
for(j=0;j<size;j++)
	{
	if(min > rasterfile[i][j][0]) 
		{
		min = rasterfile[i][j][0];
		mini = i;
		maxj = j;
		}
	if(max < rasterfile[i][j][0]) 
		{
		max = rasterfile[i][j][0];
		maxi = i;
		maxj = j;
		}
	}
/* scale the energy image */
diff = max-min;
for(i=0;i<size;i++)
for(j=0;j<size;j++)
	{
	rasterfile[i][j][0] = ((rasterfile[i][j][0]-min)/diff)*255.0;
	}
printf("Minimum real value: %f at %d,%d\n",min,mini,minj);
printf("Maximum real value: %f at %d,%d\n",max,maxi,maxj);
}



/****************************************************************/
/* This function tests the 1D FFT routine.			*/
/****************************************************************/
test()
{
int n,i,j,k;
float in[256][2];
float out[256][2],out2[256][2];
float lambda[13];
float direction;
float size1,size2,x,C,sample,sign;

n = 8;
C = 127.5;
direction = -1.0;
lambda[0] = 1.1892071;
lambda[1] = 1.4142136;
lambda[2] = 1.6817928;
lambda[3] = 2.0;
lambda[4] = 4.0;
lambda[5] = 8.0;
lambda[6] = 16.0;
lambda[7] = 32.0;
lambda[8] = 64.0;
lambda[9] = 128.0;
lambda[10] = 256.0;
lambda[11] = 512.0;
lambda[12] = 1024.0;
sample = 1.0;
for(j=0;j<7;j++) sample = sample/2.0;

for(j=0;j<=10;j++)
{
sample = 2.0*sample;
for(i=0;i<256;i++) 
	{
	x = sample*i;
	sign=1.0;
	if(i%2==1) sign= -1.0;
	in[i][0] = out[i][0] = sign*C*sin(2.0*M_PI*x/2.0);
	in[i][1] = out[i][1] = 0.0;
	}

direction = 1.0;
NLOGN(n, out, direction );
for(i=0;i<256;i++)
	{
	out2[i][0] = out[i][0];
	out2[i][1] = out[i][1];
	}
direction = -1.0;
NLOGN(n, out, direction );
printf("j: %3d  --- lambda is %f --- sample is %f\n",j,2.0,sample);
for(k=0;k<256;k++)
	printf("k=%3d in: %9.4f %9.4f out: %9.4f %9.4f out2: %9.4f %9.4f\n",k,in[k][0],in[k][1],out[k][0],out[k][1],out2[k][0],out2[k][1]);

size1 = size2 = 0.0;
for(k=0;k<256;k++)
	{
	size1 += (in[k][0]-out[k][0])*(in[k][0]-out[k][0]);
	size2 += (in[k][1]-out[k][1])*(in[k][1]-out[k][1]);
	}
size1 = sqrt(size1);
size2 = sqrt(size2);
printf("j=%d size1: %f size2: %f\n",j,size1,size2);
}
printf("\nFinished aliasing test\n");
fflush(stdout);
}

/****************************************************************/
/* This function tests the 1D FFT routine.			*/
/****************************************************************/
test1()
{
int n,i,j;
float row[16][2];
float direction;
n = 4;
direction = -1.0;
row[0][0] = 62; row[1][0] = 1; row[2][0] = 35; row[3][0] = 16;
row[4][0] = 45;  row[5][0] = 62; row[6][0] = 100; row[7][0] = 92;
row[8][0] = 123; row[9][0] = 101; row[10][0] = 145; row[11][0] = 164;
row[12][0] = 150; row[13][0] = 183; row[14][0] = 200; row[15][0] = 197;
for(i=0;i<16;i++) row[i][1] = 0.0;

printf("\bIn test1:\n");
printf("Original reals and imaginaries\n");
for(j=0;j<16;j++) printf("%f %f\n",row[j][0],row[j][1]);

NLOGN(n, row, direction );

printf("\nFT reals and imaginaries\n");
for(j=0;j<16;j++) printf("%f %f\n",row[j][0],row[j][1]);

direction = direction*-1.0;
NLOGN(n, row, direction );
printf("\nFT Inverse reals and imaginaries\n");
for(j=0;j<16;j++) printf("%f %f\n",row[j][0],row[j][1]);

printf("\nFinished test1\n");
fflush(stdout);
}

/****************************************************************/
/* This function tests the 2D FFT routine.			*/
/* The same input as in Rosenfeld and Kak, Digitial Picture     */
/* Processing, Volume 1, 2nd Edition, 1982, pp22-23 and the     */
/* output of the real and imaginary components is also the      */
/* same.							*/
/****************************************************************/
test2()
{
int i,j,image_dim;
float direction;

image_dim = 4;
direction = -1.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	rasterfile[i][j][0] = rasterfile[i][j][1] = 0.0;
rasterfile[6][6][0] = rasterfile[6][7][0] = rasterfile[6][8][0] = 1.0;
rasterfile[8][6][0] = rasterfile[8][7][0] = rasterfile[8][8][0] = 1.0;
rasterfile[7][6][0] = rasterfile[7][8][0] = 1.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	if((i+j)%2!=0) rasterfile[i][j][0] = -rasterfile[i][j][0];

fft_2d ( image_dim,rasterfile, direction );
printf("\n\nBlock of real components of 2D fft\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][0]);
	}
printf("\n\n");
}
printf("\n\nBlock of imaginary components of 2D fft\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][1]);
	}
printf("\n\n");
}

direction = -direction;
fft_2d ( image_dim,rasterfile, direction );
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	if((i+j)%2!=0) 
		{
		rasterfile[i][j][0] = -rasterfile[i][j][0];
		rasterfile[i][j][1] = -rasterfile[i][j][1];
		}
printf("\n\nBlock of real components of image\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][0]);
	}
printf("\n\n");
}
printf("\n\nBlock of imaginary components of image\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][1]);
	}
printf("\n\n");
}
fflush(stdout);
}

test4()
{
int i,j,image_dim;
float direction;

image_dim = 4;
direction = -1.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	rasterfile[i][j][0] = rasterfile[i][j][1] = 0.0;

rasterfile[6][6][0] = rasterfile[6][7][0] = rasterfile[6][8][0] = 1.0;
rasterfile[8][6][0] = rasterfile[8][7][0] = rasterfile[8][8][0] = 1.0;
rasterfile[7][6][0] = rasterfile[7][8][0] = 1.0;
rasterfile[4][3][0] = 23.9;
rasterfile[5][7][0] = 12.74;

for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	if((i+j)%2!=0) rasterfile[i][j][0] = -rasterfile[i][j][0];  

fft_2d ( image_dim,rasterfile, direction );
printf("\n\nBlock of real components of 2D fft\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][0]);
	}
printf("\n");
}
printf("\n\nBlock of imaginary components of 2D fft\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][1]);
	}
printf("\n");
}
direction = -direction;
fft_2d ( image_dim,rasterfile, direction );
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	if((i+j)%2!=0) 
		{
		rasterfile[i][j][0] = -rasterfile[i][j][0];  
		rasterfile[i][j][1] = -rasterfile[i][j][1];  
		}
printf("\n\nBlock of real components of image\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][0]);
	}
printf("\n\n");
}
printf("\n\nBlock of imaginary components of image\n");
for(i=0;i<16;i++)
{
for(j=0;j<16;j++)
	{
	if(j==8) printf("\n");
	printf("%f ",rasterfile[i][j][1]);
	}
printf("\n\n");
}
fflush(stdout);
}

/****************************************************************/
/* This function test a 64*64 set of data supplied by Rick      */
/* Gurnsey at UWO.						*/
/****************************************************************/
test3()
{
int i,j,image_dim;
float direction;

printf("\nIn test3\n");
fflush(stdout);
image_dim = 6;
direction = -1.0;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	rasterfile[i][j][0] = rasterfile[i][j][1] = 0.0;
for(i=15;i<47;i++)
	rasterfile[i][31][0] = 1.0;
printf("image_dim=%d\n",image_dim);
fft_2d ( image_dim,rasterfile, direction );
printf("\n\nBlock of real components of 2D fft\n");
for(i=0;i<64;i++)
{
for(j=0;j<64;j++)
	printf("%f ",rasterfile[i][j][0]);
printf("\n");
}
printf("\n\nBlock of imaginary components of 2D fft\n");
for(i=0;i<64;i++)
{
for(j=0;j<64;j++)
	printf("%f ",rasterfile[i][j][1]);
printf("\n");
}
printf("\n\nMagnitude Spectrum\n");
for(i=0;i<64;i++)
{
for(j=0;j<64;j++)
	printf("%f ",sqrt(rasterfile[i][j][0]*rasterfile[i][j][0]+
			  rasterfile[i][j][1]*rasterfile[i][j][1]));
printf("\n");
}
fflush(stdout);
}

/****************************************************************/
/* This function rearranges the 2D FFT results by swaping the   */
/* four blocks of real and imaginary data that comprise the     */
/* corners of the output data. The result is that the upper     */
/* left hand corner has the lowest frequency value (-1/2,-1/2), */
/* the middle has frequency (0,0) and the lower right corner    */
/* has the largest frequency (1/2-1/N,1/2-1/N), where N is the  */
/* dimension of the image, say N=256. This is the form of the   */
/* FFT output in Rosenfeld and Kak, Digital Image Processing,   */
/* Volume 1, 1982, pp22-23.					*/
/****************************************************************/
rearrange(rasterfile,n)
float rasterfile[maxrows][maxcols][2];
int n;
{
int i,j,halfi,halfj,dim;
float temp;

if(FALSE) printf("\nRearranging FT data:\n");

dim = pow(2.0,(n)*1.0);
halfi = dim/2;
halfj = dim/2;
printf("\nIn rearrange: dim %d half dim %d\n",dim,halfi);
for(i=0;i<dim/2;i++)
for(j=0;j<dim/2;j++)
{
temp = rasterfile[i][j][0];
rasterfile[i][j][0] = rasterfile[i+halfi][j+halfj][0];
rasterfile[i+halfi][j+halfj][0] = temp;
temp = rasterfile[i][j][1];
rasterfile[i][j][1] = rasterfile[i+halfi][j+halfj][1];
rasterfile[i+halfi][j+halfj][1] = temp;

temp = rasterfile[i+halfi][j][0];
rasterfile[i+halfi][j][0] = rasterfile[i][j+halfj][0];
rasterfile[i][j+halfj][0] = temp;
temp = rasterfile[i+halfi][j][1];
rasterfile[i+halfi][j][1] = rasterfile[i][j+halfj][1];
rasterfile[i+halfi][j][1] = temp;
}
}


/****************************************************************/
/* This function scales the magnitude of the complex image by   */
/* some scale factor.					        */
/****************************************************************/
scale_mag(rasterfile,n,scale)
int n;
float scale;
float rasterfile[maxrows][maxcols][2];
{
int i,j,dim;

dim = pow(2.0,(n)*1.0);
for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*scale;
	rasterfile[i][j][1] = rasterfile[i][j][1]*scale;
	}
}

/****************************************************************/
/* This function scales the phase of the complex image by some  */
/* scale factor.						*/
/****************************************************************/
scale_phase(rasterfile,n,scale)
float rasterfile[maxrows][maxcols][2];
float scale;
int n;
{
int i,j,dim;
float phase[maxrows][maxcols],mag[maxrows][maxcols];

dim = pow(2.0,(n)*1.0);
for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
	{
	mag[i][j] = sqrt(rasterfile[i][j][0]*rasterfile[i][j][0]+
			 rasterfile[i][j][1]*rasterfile[i][j][1]);
	phase[i][j] = atan2(rasterfile[i][j][1],rasterfile[i][j][0])*scale;
	}
for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
	{
	rasterfile[i][j][0] = mag[i][j]*cos(phase[i][j]);
	rasterfile[i][j][1] = mag[i][j]*sin(phase[i][j]);
	}
}


/****************************************************************/
/* Add 2 complex numbers.					*/
/****************************************************************/
void add(x,y,z)
float x[2],y[2],z[2];
{
z[0] = x[0] + y[0];
z[1] = x[1] + y[1];
}   /* add */


/****************************************************************/
/* Subtract 2 complex numbers.					*/
/****************************************************************/
void sub(x,y,z)
float x[2],y[2],z[2];
{
z[0] = x[0] - y[0];
z[1] = x[1] - y[1];
}   /* sub */



/****************************************************************/
/* Multiply 2 complex numbers.					*/
/****************************************************************/
void mult(x,y,z)
float x[2],y[2],z[2];
{
z[0] = x[0] * y[0] - x[1] * y[1];
z[1] = x[0] * y[1] + x[1] * y[0];
}   /* mult */



/****************************************************************/
/* 1D FFT algorithm, converted from Fortran. The original       */
/* Fortran algorithm comes from E.A. Roberson, Multichannel     */
/* Time Series Analysis with Digital Computer Programs, Goose   */
/* Pond Press, pp63-64.						*/
/****************************************************************/
void NLOGN(N,X,SIGN)
int N;
float X[maxn][2];
float SIGN;
{
float WK[2],HOLD[2],Q[2];
int LX,FLX;
int II,I,J,K,L,M[25],JH;
int IBLOCK,NBLOCK,LBLOCK,LBHALF,ISTART;

float FK,V;

LX = pow(2.0,(double)N);

for (I=0;I<N;I++)
     M[I] = pow(2.0,(double)(N-I-1));


for (L=0;L<N;L++)  {
     NBLOCK = pow (2.0,(double)L);
     LBLOCK=LX/NBLOCK;
     LBHALF=LBLOCK/2;
     K = 0.0;

for (IBLOCK=0;IBLOCK<NBLOCK;IBLOCK++)  {
     FK = K;
     FLX = LX;

     V = SIGN*6.2831853*FK/FLX;

     WK[0] = cos(V);
     WK[1] = sin(V);

     ISTART = LBLOCK*(IBLOCK);

     for(I=0;I<LBHALF;I++)  {
           J = ISTART + I;
           JH = J + LBHALF;
           mult (X[JH],WK,Q);
           sub  (X[J],Q,X[JH]);
           add  (X[J],Q,X[J]);
           }


     for(I=1;I<N;I++)  {
           II=I;
           if (K < M[I])
           break;
           K = K - M[I];
           }
     K = K + M[II];
     }
}

K = -1.0;

for (J=0;J<LX;J++)  {

    if (K < J)  {
    }
    else
    for(I=0;I<2;I++)  {
           HOLD[I] = X[J][I];
           X[J][I]   = X[K+1][I];
           X[K+1][I] = HOLD[I];
           }

    for(I=0;I<N;I++)  {
           II = I;
           if (K+1 < M[I]) break;
           K = K - M[I];
           }

     K = K + M[II];
     }

if (SIGN < 0.0);
else
   for (I=0;I<LX;I++)  {
        X[I][0] = X[I][0]/FLX;
        X[I][1] = X[I][1]/FLX;
        }
}  /* NLOGN */

/****************************************************************/
/* This function compute the sinc function.			*/
/****************************************************************/
float sinc(x)
float x;
{
if(fabs(x) < 0.00001) return(1.0);
return(sin(x)/x);
}

/****************************************************************/
/* Given two complex images, raster1 and raster2, compute their */
/* magnitude and phase information, average that magnitude and  */
/* phase information and then convert the average magnitude and */
/* phase information into another new complex image.		*/
/****************************************************************/
average_phase_and_mag(raster1,raster2,raster3,n)
float raster1[maxrows][maxcols][2];
float raster2[maxrows][maxcols][2];
float raster3[maxrows][maxcols][2];
int n;
{
int i,j,dim;
float mag1[maxrows][maxcols],phase1[maxrows][maxcols];
float mag2[maxrows][maxcols],phase2[maxrows][maxcols];


dim = pow(2.0,n*1.0);

for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
	{
	mag1[i][j] = sqrt(raster1[i][j][0]*raster1[i][j][0]+
			 raster1[i][j][1]*raster1[i][j][1]);
	phase1[i][j] = atan2(raster1[i][j][1],raster1[i][j][0]);
	mag2[i][j] = sqrt(raster2[i][j][0]*raster2[i][j][0]+
			 raster2[i][j][1]*raster2[i][j][1]);
	phase2[i][j] = atan2(raster2[i][j][1],raster2[i][j][0]);
	}
for(i=0;i<dim;i++)
for(j=0;j<dim;j++)
	{
	raster3[i][j][0] = (mag1[i][j]+mag2[i][j])*0.5*cos((phase1[i][j]+phase2[i][j])/2.0);
	raster3[i][j][1] = (mag1[i][j]+mag2[i][j])*0.5*sin((phase1[i][j]+phase2[i][j])/2.0);
	}
}


/****************************************************************/
/* Perform histogram equalization.				*/
/****************************************************************/
void equalize(image,dim)
float image[maxrows][maxcols][2];
int dim;
/* this routine does the histogram equalization of the image */
{
   int	ctr;		 /* a counter */
   int	data[2][maxrows];/* an array to store the number of pixels and the
			   accumulated totals for each intensity */
   int	result[maxcols];/* an array with the result */
   int	max;		/* the total number of pixels */
   int	ptr;		/* a reference pointer */
   int size;
   int row,col;


   rescale_pixels(image,dim); /* Ensure all pixels between 0 and 255 */
   /* get the total number of pixels */
   size = pow(2.0,dim*1.0);
   max = size*size;
   printf("\nIn equalize: size=%d max = %d\n",size,max);
   fflush(stdout);

   /* initialize the array */
   for (ctr=0; ctr<size; ctr++) 
	data[0][ctr] = 0;

   /* calculate the total number of pixels for each intensity */
   for (row=0; row<size; row++)
	for (col=0; col<size; col++)
	{
	ptr = (int)image[row][col][0];
	data[0][ptr]++;
	}

   /* calculate the accumulated number of pixels for each intensity */
   data[1][0] = data[0][0];
   for (ctr=1; ctr<size; ctr++)
	data[1][ctr] = data[0][ctr] + data[1][ctr-1];

   /* calculate the result */
   for (ctr=0; ctr<size; ctr++)
	result[ctr] = (int)((float)(data[1][ctr] * 255.0) / (float)max);
 
   /* update the image */
   for (row=0; row<size; row++)
        for (col=0; col<size; col++)
	{
	   ptr = (int)image[row][col][0];
	   image[row][col][0] = result[ptr];
	}

} /* end of equalize */

/****************************************************************/
/* This function computes the Laplacian of a Gaussian in        */
/* frequency space by multiplying the FT of an image by         */
/* -(u^2 + v^2) (2 pi)^2					*/
/****************************************************************/
void laplacian( image )

    float image [maxrows][maxcols][2];

{
    int u,v,uu,vv,i,j,min_u,min_v,max_u,max_v;
    float L[maxrows][maxcols],constant,min,max;


    constant = 4.0*pi*pi * 1.0/(maxrows*maxrows); 
    printf("\nComputing Laplacian in frequency space constant=%f\n",constant);

    min =  HUGE;
    max = -HUGE;
    min_u=min_v=max_u=max_v=0;
    for (u=0; u<maxrows; u++)
    for (v=0; v<maxcols; v++)
	{
	uu = u-maxrows/2.0;
	vv = v-maxcols/2.0;
	L[u][v] = -(uu*uu+vv*vv)*constant;
	if(L[u][v] < min) { min = L[u][v]; min_u=u; min_v=v; }
	if(L[u][v] > max) { max = L[u][v]; max_u=u; max_v=v; }
	}
     printf("Computing Laplacian...\n");
     printf("Maximun Laplacian value %f at %d,%d\n",max,max_u,max_v);
     printf("Miniumm Laplacian value %f at %d,%d\n",min,min_u,min_v);
     for (u=0; u<maxrows; u++)
     for (v=0; v<maxcols; v++)  
        {
        image[u][v][0] = image[u][v][0] * L[u][v];
        image[u][v][1] = image[u][v][1] * L[u][v];
        }

}   /*  laplacian */

/**********************************************************/
/* Multiply each pixel (x,y) by -1^(x,y)                  */
/**********************************************************/
translate_to_origin( image )
float image[maxrows][maxcols][2];
{
int u,v,sign,ct1,ct2;
ct1 = ct2 = 0;
for (u=0; u<maxrows; u++)
for (v=0; v<maxcols; v++)  
        {
	if((u+v)%2==1) sign = -1.0; else sign = 1.0;
	if(sign==-1) ct1++;
	if(sign== 1) ct2++;
        image[u][v][0] = image[u][v][0]*sign;
        image[u][v][1] = image[u][v][1]*sign;
        }
printf("Sign changes for translation\n");
printf("%d values of function have changed sign\n",ct1);
printf("%d values of function have remained the same\n\n",ct2);
}

translate_from_origin( image )
float image[maxrows][maxcols][2];
{
int u,v,sign;
for (u=0; u<maxrows; u++)
for (v=0; v<maxcols; v++)  
        {
	if(u+v % 2 !=0) sign = 1.0; else sign = -1.0;
        image[u][v][0] = image[u][v][0]*sign;
        image[u][v][1] = image[u][v][1]*sign;
        }
}
 
/******************************************************************/
/* Compute zero crossings of laplacian 			   */
/******************************************************************/
compute_zero_crossings(rasterfile,complex_edges,dim)
float rasterfile[numrows][numcols][2],complex_edges[numrows][numcols][2];
int dim;
{
unsigned char raster_edges[numrows][numcols];
float fpic[numrows][numcols],edges[numrows][numcols];
int pic_x,pic_y,i,j,offset;

offset = 0;
pic_x = pow(2.0,dim*1.0);
pic_y = pow(2.0,dim*1.0);
printf("Computing Zero Crossings\n");
printf("pic_x=%d pic_y=%d\n",pic_x,pic_y);

for(i=0;i<pic_x;i++)
for(j=0;j<pic_y;j++)
	{
	fpic[i][j] = rasterfile[i][j][0]; /* fpic is the laplacian values */
	}
/* Find zero crossings */
zero_cross(edges,fpic,offset,pic_x,pic_y); 
prepare_edges_for_output(edges,raster_edges,pic_x,pic_y);
for(i=0;i<pic_x;i++)
for(j=0;j<pic_y;j++)
	{
	complex_edges[i][j][0] = raster_edges[i][j];
	complex_edges[i][j][1] = 0;
	}
}


/**************************************************************************/
/* Prepare the edge maps for output as raster images			  */
/**************************************************************************/
prepare_edges_for_output(edges,raster_edges,pic_x,pic_y)
int pic_x,pic_y;
float edges[numrows][numcols];
unsigned char raster_edges[numrows][numcols];
{
int i,j,k;
float max,min;
int count,no_count,pixel;

printf("\nEdge Report:\n");
/* Flip edge map: 0 goes to WHITE and 1 goes to BLACK */
count = no_count = 0;
for(i=0;i<pic_x;i++)
for(j=0;j<pic_y;j++)
	{
	raster_edges[i][j] = WHITE;
	if(edges[i][j]!=0.0) 
		{ 
		pixel = (1.0-edges[i][j])*WHITE;
		if(pixel > 255) pixel = 255;
		if(pixel < 0) pixel = 0;
		raster_edges[i][j] = pixel;
		count++; 
		}
	else no_count++;
	}
printf("%d edgels and %d non-edgels in the image\n",count,no_count);
/* Print Black Border */
for(i=0;i<pic_x;i++) 
    { raster_edges[i][0] = raster_edges[i][pic_y-1] = BLACK; }
for(j=0;j<pic_y;j++) 
    { raster_edges[0][j] = raster_edges[pic_x-1][j] = BLACK; } 
printf("Raster edge map created\n");
fflush(stdout);
}



/**************************************************************/
/* Compute zero crossing of the Laplacian.		      */
/* Older version - gives best edge maps			      */
/**************************************************************/
zero_cross(edges,laplacian,offset,pic_x,pic_y)
float edges[numrows][numcols];
float laplacian[numrows][numcols];
int offset,pic_x,pic_y;
{
int sign1,sign2,j,p1,p2,p3;

for(j=offset;j<pic_x-offset;j++) 
	{ 
        p1 = offset;
        while(p1<pic_y-offset && sign(laplacian[j][p1])==0) p1++; 
        sign1 = sign(laplacian[j][p1]);

	while(sign1 && p1<pic_y-offset) 
		{
		p2 = p1+1;
               	while(p2 < pic_y-offset &&
		      sign(laplacian[j][p2]) != -sign1) p2++;
		p3 = p2-1;
		sign1 = sign(laplacian[j][p2]);
		while(p3 > p1 && sign1 != -sign(laplacian[j][p3])) p3--; 
		
		/* Check for case when edge of image reached, 
		   but no edge present */
		if(sign(laplacian[j][p3]) != -sign(laplacian[j][p2])) break;	

		/* Set edge */
		if((p2-p3) % 2 == 0) edges[j][p3+(p2-p3)/2] = 1.0;
		else 
                if(fabs(laplacian[j][p3])<fabs(laplacian[j][p2]))
                       	edges[j][p3+(p2-p3)/2] = 1.0;
                else edges[j][p3+(p2-p3)/2+1] = 1.0;
		p1 = p2;
                sign1 = sign(laplacian[j][p1]);
	    	}
	}	

for(j=offset;j<pic_y-offset;j++)
	{
        p1 = offset;
        while(p1<pic_x-offset && sign(laplacian[p1][j])==0) p1++; 
        sign1 = sign(laplacian[p1][j]);
 
        while(sign1 && p1 < pic_x-offset) 
		{
               	p2 = p1+1;
		while(p2 < pic_x-offset && 
		      sign1 != -sign(laplacian[p2][j])) p2++; 
		p3 = p2-1;
		sign1 = sign(laplacian[p2][j]);
		while(p3>p1 && sign1!=-sign(laplacian[p3][j])) p3--; 
                 
                /* Check for case when edge of image reached, 
		   but no edge present **/ 
                if(sign(laplacian[p3][j]) != -sign(laplacian[p2][j])) break;   
 
                /* Set edge */ 
                if((p2-p3) % 2 == 0) edges[p3+(p2-p3)/2][j] = 1.0;
                else   
                if(fabs(laplacian[p3][j])<fabs(laplacian[p2][j]))
                	edges[p3+(p2-p3)/2][j] = 1.0;
                else edges[p3+(p2-p3)/2+1][j] = 1.0;
                p1 = p2; 
                sign1 = sign(laplacian[p1][j]);
            	} 
       	}  
}

/**************************************************************/
/* Determine sign of x: -1 for negative, +1 for positive and  */
/* 0 for zero.						      */
/**************************************************************/
int sign(x)
float x;
{
if(x> 0.000001) return( 1);
if(x<-0.000001) return(-1);
return(0);
}

test0()
{
int i,inc,j;
float x,start;
float data[32][2],y[32];

for(j=1;j<=20;j++)
{
start = 0.0;
inc = j;
x = start;
for(i=0;i<32;i++)
	{
	data[i][0] = y[i]=127.5*sin(2.0*M_PI/16.0 * x);
	if(i%2==1) data[i][0] *= -1;
	data[i][1] = 0.0;
	x += inc;
	}
NLOGN(5, data, ffttransform );
printf("j=%d\n",j);
printf("    x        y        Freq          Real         Imag\n");
x = start;
for(i=0;i<32;i++)
	{
	printf("%8.4f  %8.4f  %8.4f  %8.4f %8.4f\n",x,y[i],(-0.5+i*1.0/32.0)/(j*1.0),data[i][0],data[i][1]);
	x += inc;
	}
printf("\n");
}
}

/************************************************************************/
/* Copy rasterfile1 into rasterfile2					*/
/************************************************************************/
copy_data(rasterfile1,rasterfile2)
float rasterfile1[maxrows][maxcols][2];
float rasterfile2[maxrows][maxcols][2];
{
int i,j;
for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
	{
	rasterfile2[i][j][0] = rasterfile1[i][j][0];
	rasterfile2[i][j][1] = rasterfile1[i][j][1];
	}
}

/****************************************************************/
/* This function performs low pass filtering as specified in    */
/* Duda and Hart, Pattern Classification and Scene Analysis,    */
/* John Wiley and Sons, 1973, pp313.				*/
/* Written by Tom Emberson, 1994 				*/
/****************************************************************/
void tom_high_pass(rasterfile)
float rasterfile[maxrows][maxcols][2];
{
float H[maxrows][maxcols];
int i,j;
float u,v;
double sixteen;

sixteen = 16.0;

for(i=0;i<maxrows;i++)
for(j=0;j<maxcols;j++)
{
u = i*1.0/maxrows-0.5;
v = j*1.0/maxcols-0.5;
if(u*u+v*v !=0.0)
H[i][j] = 1.0 / (1.0 + pow (0.15 / sqrt(u*u + v*v), (double) 2.0));
else H[i][j] = 0.0;
}

for(i=0;i<maxrows;i++)
for(j=0;j<maxrows;j++)
	{
	rasterfile[i][j][0] = rasterfile[i][j][0]*H[i][j];
	rasterfile[i][j][1] = rasterfile[i][j][1]*H[i][j];
	}
}