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