/************************************************************************** * * * Copyright (c) 1996 Michael Richmond and Richard Treffers * * * * This software may be copied and distributed for * * educational, research and not for profit services * * provided that this copyright and statement are * * included in all such copies. * * * **************************************************************************/ /* Fourier transform image FFT2D fname [DARK=] [SCALE=] [REAL] [IMAG] [INVERSE] [BOX=] [VERBOSE] 1/30/90 MAXSHORT/2 changed to MAXINT/2 - MWR 9/12/91 changed for RT -RRT 2/21/92 change MAXLONG to LONG_MAX to agree with limits.h MAXINT to INT_MAX 9/23/92 usage, exit(0),dead wood -rrt 10/19/92 HUGEFLOAT now defined here as float 3/19/93 added ability to input imaginary transform 2/25/94 nr.h added after math.h 4/22/94 ournr.h used instead of nr.h */ #include "pcvista.h" #include "fits.h" #include #include #include #include "ournr.h" #include #include #include /* * I don't know why we need to include this under SunOs 5.5 */ #ifndef LINUX extern double hypot(double x, double y); #endif #ifdef PROTO void main(int, char **); static void get_fft_data(char *real,char *imag,int box); static int clip(double); static double find_scale(void); static int power_of_two(int); static int de_quad(int,int); #else static void get_fft_data(); static int clip(); static double find_scale(); static int power_of_two(); static int de_quad(); #endif #define OUTNAME "fft" #define NDIM 2 /* we are doing a two-dimensional transform */ #define AMP 0 /* output mode */ #define REAL 1 #define IMAG 2 #define PHASE 3 #define POWER 4 float *data; /* floating point data array */ int nn[NDIM]; /* array of size */ int nr, nc; /* power of 2 closest size of xform */ static mode = AMP; static int verbose = 0; short *short_data; /* integer buffer for real data */ short *short_idata; /* integer buffer for imaginary data */ #define USAGE \ "usage: fft2d file [scale=] [real] [imag] [inverse] [box=] [verbose] \n [imag=] [power] [outfile=] [help]" void main(argc, argv) int argc; char *argv[]; { int i, j, boxno = -1, got_scale = 0; int index; FITS_HANDLE fout; double scale = 1.0, val; char foutname[NBUF], *gotit; static char real_file[NBUF]; static char imag_file[NBUF]; int direction = -1; if (argc == 1) { error(-1, USAGE); } if (strcmp(argv[1], "help") == 0) { printf("%s\n", USAGE); exit(0); } strcpy(foutname,OUTNAME); for(i = 2; i < argc; i++) { if (keyword("help", argv[i])) { printf("%s\n", USAGE); exit(0); } if (keyword("inverse", argv[i])) { direction = 1; continue; } if ((gotit = find("box", argv[i])) != NULL) { boxno = (int) (evaluate(gotit) + SMALL); continue; } if ((gotit = find("scale", argv[i])) != NULL) { got_scale = 1; scale = evaluate(gotit); continue; } if ((gotit = find("imag", argv[i])) != NULL) { strcpy(imag_file,gotit); continue; } if ((gotit = find("outfile", argv[i])) != NULL) { strcpy(foutname,gotit); continue; } if (keyword("real", argv[i])) { mode = REAL; continue; } if (keyword("power", argv[i])) { mode = POWER; continue; } if (keyword("imag", argv[i])) { mode = IMAG; continue; } if (keyword("verbose", argv[i])) { verbose = 1; continue; } fprintf(stderr,"unknown keyword '%s'\n",argv[i]); exit(1); } strcpy(real_file, argv[1]); get_fft_data(real_file, imag_file,boxno); fout = fits_open(foutname, "w", &nr, &nc); if (verbose) printf("transforming\n"); fourn(data - 1,nn - 1, NDIM, direction); /* this routine [1..2*n] */ if (!got_scale) scale = find_scale(); if ((scale != 1.0) && (verbose)) printf("scale factor=%.4f\n", scale); /* placate compiler */ val = 0.0; for (i = 0; i < nr; i++) { for (j = 0; j < nc; j++) { index = de_quad(i, j); switch (mode) { case AMP: val=scale*hypot((double)data[index],(double)data[index+1]); break; case POWER: val=scale* (data[index]*data[index]+data[index+1]*data[index+1]); break; case REAL: val = scale*data[index]; break; case IMAG: val = scale*data[index + 1]; break; default: error(1, "fft2d: impossible mode"); } short_data[j] = clip(val); } fits_put_data(fout, i, 0, short_data, nc); } fits_close(fout); exit(0); } /* this routine allocates the storage for 1. the floating data array 2. the integer buffer array need to read FITS data it then fills the floating buffer array in sequentially in the manner described on pg 468 of N.R.C. */ static void get_fft_data(real_file, imag_file, boxno) char *real_file,*imag_file; int boxno; { int i, j, k; int nrow, ncol; /* native size of data array */ int imag_row,imag_col; int row, tr, tc; int sr = 0, sc = 0; FITS_HANDLE fh,fimag=0; int nbytes; fh = fits_open(real_file, "r", &nrow, &ncol); /* get the data file */ if (boxno != -1) { /* find the box parameters */ tr = nrow; tc = ncol; /* temporary storage */ if(getbox(boxno, &sr, &sc, &nrow, &ncol)) error(-1, "box not found"); check_box(boxno, tr, tc); } if(imag_file[0]){ /* an imaginary file is specified */ fimag = fits_open(imag_file, "r", &imag_row, &imag_col); if(imag_row !=nrow) error(1,"imaginary file doesn't have the same number of rows"); if(imag_col !=ncol) error(1,"imaginary file doesn't have the same number of cols"); } nc = power_of_two(ncol); /* find the nearest power of two in size */ nr = power_of_two(nrow); if ((verbose) && ((nc != ncol) || (nr != nrow))) printf("padding with zeros\n"); nn[0] = nc; /* fill in the array size */ nn[1] = nr; nbytes = 2*nr*nc*sizeof(float); /* create data base array */ data = (float *) malloc(nbytes); if (data == NULL){ fprintf(stderr, "can't allocate %d bytes", nbytes); exit(1); } short_data = (short *) malloc(nc*sizeof(short)); short_idata = (short *) calloc(nc,sizeof(short)); if (short_data == NULL) error(-1, "error allocating real array"); if (short_idata == NULL) error(-1, "error allocating imaginary array"); for (row = sr, j = 0, k = 0; k < nr; k++, row++) { fits_get_data(fh, row, sc, short_data, ncol); if(fimag) fits_get_data(fimag, row, sc, short_idata, ncol); for (i = 0; i < ncol; i++) { /* spread data out */ data[j++] = (float)short_data[i] ; /* real */ data[j++] = (float) short_idata[i]; /* imaginary */ } for( ; i < nc; i++) { /* pad with zeroes */ data[j++] = (float) 0.0; data[j++] = (float) 0.0; } } fits_close(fh); } static int power_of_two(in) int in; { int i; for (i = 1; i < in; i *= 2) /* find next power of two */ ; return(i); } static int clip(x) /* converts doubles to int rationally */ double x; { if ((int) x > MAX_DATA_VAL) return(MAX_DATA_VAL); if ((int) x < MIN_DATA_VAL) return(MIN_DATA_VAL); return((int) x); } /* this routine determines the data array index so that the center frequency is in the center of the picture THIS ROUTINE IS NOT NECESSARILY CORRECT OR EFFICIENT */ static int de_quad(row, col) int row, col; { row = nr/2 - row; col = nc/2 - col; if (col < 0) col += nc; if (row < 0) row += nr; return(2*((row*nc) + col)); } /* determines the scale factor */ static double find_scale() { int i, max_ind; double max = (double) (1-LONG_MAX), min = (double) LONG_MAX; double val, bigmax, scale; /* placate compiler */ val = 0.0; max_ind = 2*nr*nc; for (i = 0; i < max_ind; i++) { switch (mode) { case AMP: val = sqrt((double)(data[i]*data[i] + data[i + 1]*data[i + 1])); break; case REAL: val = (double) data[i]; break; case IMAG: val = (double) data[i + 1]; break; default: error(1, "find_scale has impossible mode"); } if (val > max) max = val; if (val < min) min = val; } if (max > fabs(min)) bigmax = max; else bigmax = -1.0*min; if (bigmax > (double) INT_MAX) scale = (double)INT_MAX/bigmax; else scale = 1.0; return(scale); }