fft2d.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This low-level function computes a discrete one-dimensional Fourier      */
00009 /*   transform.                                                               */
00010 /*                                                                            */
00011 /*----------------------------------------------------------------------------*/
00012 /*-Interface Information------------------------------------------------------*/
00013 void FFT1D (
00014 short  option,  /*  I   Switch, indicating the direction of the transform:    */
00015                 /*      FORWARD   -   forward Fourier transform is computed.  */
00016                 /*      INVERSE   -   inverse Fourier transform is computed.  */
00017 short  size,    /*  I   Number of complex data points.                        */
00018 double *sine,   /*  I   Array containing sine function values.                */
00019 double *cosn,   /*  I   Array containing cosine function values.              */
00020 double *array   /* I/O  Pointer to the complex data array[2*size].            */
00021 /*----------------------------------------------------------------------------*/
00022 ) { register short n, m, l, k, j, i;
00023     register double temp, e_real, e_imag, o_real, o_imag;
00024 
00025     if (TIMES) TIMING(T_FFT1D);
00026 
00027     /* scramble the array */
00028     for (i=j=0; j<size; j+=1,i+=k) {
00029         if (i > j) {
00030             temp = array[2*i];
00031             array[2*i] = array[2*j];
00032             array[2*j] = temp;
00033             temp = array[2*i+1];
00034             array[2*i+1] = array[2*j+1];
00035             array[2*j+1] = temp;
00036         }
00037         for (k=size/2; 1<=k && k<=i; i-=k,k/=2);
00038     }
00039 
00040     /* compute Fourier coefficients */
00041     for (i=2; i<=size; i*=2) {
00042         for (j=0; j<size/i; j++) {
00043             for (k=0; k<i/2; k++) {
00044                 l = 2*(j*i+k);
00045                 m = 2*(j*i+i/2+k);
00046                 n = (short)((long)k*(long)size/(long)i);
00047                 e_real = array[l  ];
00048                 e_imag = array[l+1];
00049                 o_real = array[m]*cosn[n] - array[m+1]*sine[n];
00050                 o_imag = array[m]*sine[n] + array[m+1]*cosn[n];
00051                 array[l  ] = (e_real+o_real)/2.0;
00052                 array[l+1] = (e_imag+o_imag)/2.0;
00053                 array[m  ] = (e_real-o_real)/2.0;
00054                 array[m+1] = (e_imag-o_imag)/2.0;
00055             }
00056         }
00057     }
00058 
00059     if (option == INVERSE) for (i=0; i<2*size; array[i++]*=(double)size);
00060 
00061     the_end:
00062     if (TIMES) TIMING(T_EXIT);
00063 }
00064 
00065 /*-Copyright Information------------------------------------------------------*/
00066 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00067 /*----------------------------------------------------------------------------*/
00068 /*-General Information--------------------------------------------------------*/
00069 /*                                                                            */
00070 /*   This function computes the Fourier transform of an image.                */
00071 /*                                                                            */
00072 /*----------------------------------------------------------------------------*/
00073 /*-Background Information-----------------------------------------------------*/
00074 /*                                                                            */
00075 /*   Bergland, G.D.:                                                          */
00076 /*   "A guided tour of the fast Fourier transform."                           */
00077 /*   IEEE Spectrum, Vol. 6, No. 7, July 1969                                  */
00078 /*                                                                            */
00079 /*----------------------------------------------------------------------------*/
00080 /*-Interface Information------------------------------------------------------*/
00081 void FFT2D (
00082 IMAGE *in,      /*  I   Pointer to the input image.                           */
00083 short option,   /*  I   Switch, indicating the direction of the transform:    */
00084                 /*      FORWARD   -   forward Fourier transform is computed.  */
00085                 /*      INVERSE   -   inverse Fourier transform is computed.  */
00086 IMAGE **out     /*  O   Address of a pointer to the output image.             */
00087 /*----------------------------------------------------------------------------*/
00088 ) { register short i, j, k;
00089     char   msg[SLEN];
00090     double *sine=0, *cosn=0, *array=0, pinc, psum;
00091     PIXEL  temp;
00092     IMAGE  *img;
00093 
00094     if (TIMES) TIMING(T_FFT2D);
00095     if (NAMES) {
00096         MESSAGE('I',"");
00097         MESSAGE('I',"FFT2D");
00098         MESSAGE('I',"");
00099         sprintf(msg," Input image:   %s",in->text);
00100         MESSAGE('I',msg);
00101         if (option == FORWARD) {
00102             MESSAGE('I'," Forward FFT");
00103         } else if (option == INVERSE) {
00104             MESSAGE('I'," Inverse FFT");
00105         }
00106         MESSAGE('I'," ...............");
00107     }
00108 
00109     /* check input */
00110     if (!CHECKIMG(in)) {
00111         MESSAGE('E'," Can't identify image.");
00112         goto the_end;
00113     } else if (in->nlin != 1L<<rnd(log10((double)in->nlin)/log10(2.0))) {
00114         MESSAGE('E'," Number of lines must be a power of two.");
00115         goto the_end;
00116     } else if (in->npix != 1L<<rnd(log10((double)in->npix)/log10(2.0))) {
00117         MESSAGE('E'," Number of pixels/line must be a power of two.");
00118         goto the_end;
00119     } else if (option != FORWARD  &&  option != INVERSE) {
00120         MESSAGE('E'," Transformation option must be either FORWARD or INVERSE.");
00121         goto the_end;
00122     }
00123 
00124     /* create image of appropriate size */
00125     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,(short)((option == FORWARD  ?  2.0 : 0.5)*in->npix),out);
00126     if (!*out) goto the_end;
00127     img = option == FORWARD  ?  *out : in;
00128 
00129     /* allocate and initialize trigonometric tables */
00130     i = max(img->nlin,img->npix/2);
00131     if (!(sine=(double *)malloc(i*sizeof(double)))) {
00132         MESSAGE('E'," Memory request failed.");
00133         goto the_end;
00134     }
00135     if (!(cosn=(double *)malloc(i*sizeof(double)))) {
00136         MESSAGE('E'," Memory request failed.");
00137         goto the_end;
00138     }
00139     for (j=0; j<i; j++) {
00140         sine[j] = (double)option*sin((2.0*PI*(double)j)/(double)i);
00141         cosn[j] =                cos((2.0*PI*(double)j)/(double)i);
00142     }
00143 
00144     /* allocate array to hold one row/column of pixels */
00145     if (!(array=(double *)malloc(2*i*sizeof(double)))) {
00146         MESSAGE('E'," Memory request failed.");
00147         goto the_end;
00148     }
00149 
00150     /* initialize progress indicator */
00151     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00152     pinc = 1.0/(double)img->nbnd/(double)(img->nlin+img->npix/2);
00153 
00154     /* Fourier transform bands */
00155     for (i=0; i<img->nbnd; i++) {
00156 
00157         /* convert from real to "complex" */
00158         if (option == FORWARD) {
00159             for (j=0; j<(*out)->nlin; j+=1) {
00160                 for (k=0; k<(*out)->npix; k+=2) {
00161                     (*out)->data[i][j][k  ] = in->data[i][j][k/2];
00162                     (*out)->data[i][j][k+1] = (PIXEL)0;
00163                 }
00164             }
00165         }
00166 
00167         /* fold */
00168         for (j=0; j<img->nlin/2; j++) {
00169             for (k=0; k<img->npix/2; k++) {
00170                 temp = img->data[i][j][k];
00171                 img->data[i][j][k] = img->data[i][img->nlin/2+j][img->npix/2+k];
00172                 img->data[i][img->nlin/2+j][img->npix/2+k] = temp;
00173                 temp = img->data[i][img->nlin/2+j][k];
00174                 img->data[i][img->nlin/2+j][k] = img->data[i][j][img->npix/2+k];
00175                 img->data[i][j][img->npix/2+k] = temp;
00176             }
00177         }
00178 
00179         /* Fourier transform rows */
00180         for (j=0; j<img->nlin; j++) {
00181             for (k=0; k<img->npix; k++) {
00182                 array[k] = (double)img->data[i][j][k];
00183             }
00184             FFT1D(option,img->npix/2,sine,cosn,array);
00185             for (k=0; k<img->npix; k++) {
00186                 img->data[i][j][k] = (PIXEL)array[k];
00187             }
00188             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00189         }
00190 
00191         /* Fourier transform columns */
00192         for (k=0; k<img->npix; k+=2) {
00193             for (j=0; j<img->nlin; j+=1) {
00194                 array[2*j  ] = (double)img->data[i][j][k  ];
00195                 array[2*j+1] = (double)img->data[i][j][k+1];
00196             }
00197             FFT1D(option,img->nlin,sine,cosn,array);
00198             for (j=0; j<img->nlin; j+=1) {
00199                 img->data[i][j][k  ] = (PIXEL)array[2*j  ];
00200                 img->data[i][j][k+1] = (PIXEL)array[2*j+1];
00201             }
00202             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00203         }
00204 
00205         /* unfold */
00206         for (j=0; j<img->nlin/2; j++) {
00207             for (k=0; k<img->npix/2; k++) {
00208                 temp = img->data[i][j][k];
00209                 img->data[i][j][k] = img->data[i][img->nlin/2+j][img->npix/2+k];
00210                 img->data[i][img->nlin/2+j][img->npix/2+k] = temp;
00211                 temp = img->data[i][img->nlin/2+j][k];
00212                 img->data[i][img->nlin/2+j][k] = img->data[i][j][img->npix/2+k];
00213                 img->data[i][j][img->npix/2+k] = temp;
00214             }
00215         }
00216 
00217         /* convert from "complex" to real */
00218         if (option == INVERSE) {
00219             for (j=0; j<in->nlin; j+=1) {
00220                 for (k=0; k<in->npix; k+=2) {
00221                     (*out)->data[i][j][k/2] = in->data[i][j][k];
00222                 }
00223             }
00224         }
00225     }
00226 
00227     the_end:
00228     if (sine)  free(sine);
00229     if (cosn)  free(cosn);
00230     if (array) free(array);
00231     if (TIMES) TIMING(T_EXIT);
00232 }