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