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