fht2d.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1989 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This low-level function computes a discrete one-dimensional Hartley      */
00009 /*   transform.  The Fourier transform may also be easily computed from       */
00010 /*   the results.                                                             */
00011 /*                                                                            */
00012 /*----------------------------------------------------------------------------*/
00013 /*-Interface Information------------------------------------------------------*/
00014 static void FHT1D (
00015 short  option,  /*  I   Switch, indicating the direction of the transform:    */
00016                 /*      FORWARD   -   forward Hartley transform is computed.  */
00017                 /*      INVERSE   -   inverse Hartley transform is computed.  */
00018 short  size,    /*  I   Number of data points.                                */
00019 short  p,       /*  I   Logarithm (base 2) of the number of data points.      */
00020 short  *a,      /*  I   Permutation index array[p/2+p%2].                     */
00021 short  *pwr,    /*  I   Array containing powers of two.                       */
00022 double *sine,   /*  I   Array containing sine function values.                */
00023 double *tang,   /*  I   Array containing tangent function values.             */
00024 double *array   /* I/O  Pointer to the data array[size].                      */
00025 /*----------------------------------------------------------------------------*/
00026 ) { register short o, n, m, l, k, j, i;
00027     register double t, u, v, w, x, y;
00028 
00029     if (TIMES) TIMING (T_FHT1D);
00030 
00031     a[0] = 0;
00032     a[1] = 1;
00033     for (i=2; i<=p/2+p%2; i++) {
00034         for (j=0; j<pwr[i-1]; j++) {
00035             a[j+pwr[i-1]] = (a[j] *= 2) + 1;
00036         }
00037     }
00038     for (i=1; i<pwr[p/2]; i++) {
00039         k = i;
00040         l = m = pwr[p/2]*a[i];
00041         t        = array[k];
00042         array[k] = array[m];
00043         array[m] = t;
00044         for (j=1; j<=a[i]-1; j++) {
00045             k += pwr[p/2];
00046             m  = l + a[j];
00047             t        = array[k];
00048             array[k] = array[m];
00049             array[m] = t;
00050         }
00051     }
00052     for (i=0; i<size-1; i+=2) {
00053         t = array[i] + array[i+1];
00054         u = array[i] - array[i+1];
00055         array[i  ] = t;
00056         array[i+1] = u;
00057     }
00058     for (i=0; i<size-3; i+=4) {
00059         t = array[i  ] + array[i+2];
00060         u = array[i+1] + array[i+3];
00061         v = array[i  ] - array[i+2];
00062         w = array[i+1] - array[i+3];
00063         array[i  ] = t;
00064         array[i+1] = u;
00065         array[i+2] = v;
00066         array[i+3] = w;
00067     }
00068     for (i=2; i<p; i++) {
00069         for (j=0; j<size; j+=pwr[i+1]) {
00070             t = array[j] + array[j+pwr[i]];
00071             u = array[j] - array[j+pwr[i]];
00072             array[j       ] = t;
00073             array[j+pwr[i]] = u;
00074             l = j + 1;
00075             m = j + pwr[i] - 1;
00076             n = l + pwr[i];
00077             o = m + pwr[i];
00078             for (k=pwr[p-i-1]; k<=size/4; k+=pwr[p-i-1]) {
00079                 t = array[n] + tang[k]*array[o];
00080                 x = array[o] - sine[k]*t;
00081                 y = t        + tang[k]*x;
00082                 t = array[l] + y;
00083                 u = array[m] - x;
00084                 v = array[l] - y;
00085                 w = array[m] + x;
00086                 array[l++] = t;
00087                 array[m--] = u;
00088                 array[n++] = v;
00089                 array[o--] = w;
00090             }
00091         }
00092     }
00093 
00094     if (option == FORWARD) for (i=0; i<size; array[i++]/=(double)size);
00095 
00096     the_end:
00097     if (TIMES) TIMING(T_EXIT);
00098 }
00099 
00100 /*-Copyright Information------------------------------------------------------*/
00101 /* Copyright (c) 1989 by the University of Arizona Digital Image Analysis Lab */
00102 /*----------------------------------------------------------------------------*/
00103 /*-General Information--------------------------------------------------------*/
00104 /*                                                                            */
00105 /*   This function computes the Hartley transform of an image.                */
00106 /*                                                                            */
00107 /*----------------------------------------------------------------------------*/
00108 /*-Background Information-----------------------------------------------------*/
00109 /*                                                                            */
00110 /*   Bracewell, R.N.:                                                         */
00111 /*   "The Hartley Transform."                                                 */
00112 /*   Oxford University Press, New York, 1986                                  */
00113 /*                                                                            */
00114 /*----------------------------------------------------------------------------*/
00115 /*-Interface Information------------------------------------------------------*/
00116 void FHT2D (
00117 IMAGE *in,      /*  I   Pointer to the input image.                           */
00118 short option,   /*  I   Switch, indicating the direction of the transform:    */
00119                 /*      FORWARD   -   forward Hartley transform is computed.  */
00120                 /*      INVERSE   -   inverse Hartley transform is computed.  */
00121 IMAGE **out     /*  O   Address of a pointer to the output image.             */
00122 /*----------------------------------------------------------------------------*/
00123 ) { register short i, j, k, p, *a=0, *pwr=0;
00124     char   msg[SLEN];
00125     double *sine=0, *tang=0, *array=0, pinc, psum;
00126     PIXEL  t, u, v, w;
00127     IMAGE  *img;
00128 
00129     if (TIMES) TIMING(T_FHT2D); 
00130     if (NAMES) {
00131         MESSAGE('I',"");
00132         MESSAGE('I',"FHT2D");
00133         MESSAGE('I',"");
00134         sprintf(msg," Input image:   %s",in->text);
00135         MESSAGE('I',msg);
00136         if (option == FORWARD) {
00137             MESSAGE('I'," Forward FHT");
00138         } else if (option == INVERSE) {
00139             MESSAGE('I'," Inverse FHT");
00140         }
00141         MESSAGE('I'," ...............");
00142     }
00143 
00144     /* check input */
00145     if (!CHECKIMG(in)) {
00146         MESSAGE('E'," Can't identify image.");
00147         goto the_end;
00148     } else if (in->nlin != 1L<<rnd(log10((double)in->nlin)/log10(2.0))) {
00149         MESSAGE('E'," Number of lines must be a power of two.");
00150         goto the_end;
00151     } else if (in->npix != 1L<<rnd(log10((double)in->npix)/log10(2.0))) {
00152         MESSAGE('E'," Number of pixels/line must be a power of two.");
00153         goto the_end;
00154     } else if (option != FORWARD  &&  option != INVERSE) {
00155         MESSAGE('E'," Transformation option must be either FORWARD or INVERSE.");
00156         goto the_end;
00157     } else if (option == FORWARD  &&  in->nlin != in->npix   ||   option == INVERSE  &&  in->nlin != in->npix/2) {
00158         MESSAGE('E'," Number of lines must be equal to number of pixels/line.");
00159         goto the_end;
00160     }
00161 
00162     /* create image of appropriate size */
00163     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,(short)((option == FORWARD  ?  2.0 : 0.5)*in->npix),out);
00164     if (!*out) goto the_end;
00165     img = option == FORWARD  ?  in : *out;
00166 
00167     /* special case */
00168     if (img->nlin == 1) {
00169         (*out)->data[0][0][0] = in->data[0][0][0];
00170         if (option == FORWARD) {
00171             (*out)->data[0][0][1] = 0.0;
00172         }
00173         goto the_end;
00174     }
00175 
00176     /* allocate permutation index array */
00177     i = max(img->nlin,img->npix);
00178     p = rnd(log10((double)i)/log10(2.0));
00179     if (!(a=(short *)malloc((1<<p/2+p%2)*sizeof(short)))) {
00180         MESSAGE('E'," Memory request failed.");
00181         goto the_end;
00182     }
00183 
00184     /* allocate and initialize table with powers of two */
00185     if (!(pwr=(short *)malloc((p+1)*sizeof(short)))) {
00186         MESSAGE('E'," Memory request failed.");
00187         goto the_end;
00188     }
00189     for (j=0; j<=p; pwr[j]=1<<j,j++);
00190 
00191     /* allocate and initialize trigonometric tables */
00192     if (!(sine=(double *)malloc((i/4+1)*sizeof(double)))) {
00193         MESSAGE('E'," Memory request failed.");
00194         goto the_end;
00195     }
00196     if (!(tang=(double *)malloc((i/4+1)*sizeof(double)))) {
00197         MESSAGE('E'," Memory request failed.");
00198         goto the_end;
00199     }
00200     for (j=0; j<=i/4; j++) {
00201         sine[j] = sin((2.0*PI*(double)j)/(double)i);
00202         tang[j] = tan((    PI*(double)j)/(double)i);
00203     }
00204             
00205     /* allocate array to hold one row/column of pixels */
00206     if (!(array=(double *)malloc(i*sizeof(double)))) {
00207         MESSAGE('E'," Memory request failed.");
00208         goto the_end;
00209     }
00210 
00211     /* initialize progress indicator */
00212     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00213     pinc = 1.0/(double)img->nbnd/(double)(img->nlin+img->npix);
00214 
00215     /* Hartley transform bands */
00216     for (i=0; i<(*out)->nbnd; i++) {
00217 
00218         if (option == FORWARD) {
00219             /* fold */
00220             for (j=0; j<(*out)->nlin/2; j++) {
00221                 for (k=0; k<(*out)->npix/4; k++) {
00222                     (*out)->data[i][j][k] = in->data[i][in->nlin/2+j][in->npix/2+k];
00223                     (*out)->data[i][(*out)->nlin/2+j][k] = in->data[i][j][in->npix/2+k];
00224                     (*out)->data[i][j][(*out)->npix/4+k] = in->data[i][in->nlin/2+j][k];
00225                     (*out)->data[i][(*out)->nlin/2+j][(*out)->npix/4+k] = in->data[i][j][k];
00226                 }
00227             }
00228         } else {
00229             /* compute Hartley from Fourier coefficients and fold */
00230             for (j=0; j<(*out)->nlin/2; j++) {
00231                 for (k=0; k<(*out)->npix/2; k++) {
00232                     (*out)->data[i][j][k] = in->data[i][in->nlin/2+j][in->npix/2+2*k] - in->data[i][in->nlin/2+j][in->npix/2+2*k+1];
00233                     (*out)->data[i][(*out)->nlin/2+j][k] = in->data[i][j][in->npix/2+2*k] - in->data[i][j][in->npix/2+2*k+1];
00234                     (*out)->data[i][j][(*out)->npix/2+k] = in->data[i][in->nlin/2+j][2*k] - in->data[i][in->nlin/2+j][2*k+1];
00235                     (*out)->data[i][(*out)->nlin/2+j][(*out)->npix/2+k] = in->data[i][j][2*k] - in->data[i][j][2*k+1];
00236                 }
00237             }
00238         }
00239 
00240         /* Hartley transform rows */
00241         for (j=0; j<img->nlin; j++) {
00242             for (k=0; k<img->npix; k++) {
00243                 array[k] = (double)(*out)->data[i][j][k];
00244             }
00245             FHT1D(option,img->npix,p,a,pwr,sine,tang,array);
00246             for (k=0; k<img->npix; k++) {
00247                 (*out)->data[i][j][k] = (PIXEL)array[k];
00248             }
00249             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00250         }
00251 
00252         /* Hartley transform columns */
00253         for (k=img->npix-1; k>=0; k--) {
00254             for (j=0; j<img->nlin; j++) {
00255                 array[j] = (double)(*out)->data[i][j][k];
00256             }
00257             FHT1D(option,img->nlin,p,a,pwr,sine,tang,array);
00258             if (option == FORWARD) {
00259                 for (j=0; j<img->nlin; j++) { 
00260                     (*out)->data[i][j][2*k] = (PIXEL)array[j];
00261                 }
00262             } else {
00263                 for (j=0; j<img->nlin; j++) { 
00264                     (*out)->data[i][j][k] = (PIXEL)array[j];
00265                 }
00266             }
00267             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00268         }
00269 
00270         if (option == FORWARD) {
00271             /* compute Fourier from Hartley coefficients */
00272             (*out)->data[i][0][1] = (*out)->data[i][(*out)->nlin/2][1] = (*out)->data[i][0][(*out)->npix/2+1] = (*out)->data[i][(*out)->nlin/2][(*out)->npix/2+1] = 0.0;
00273             for (k=2; k<(*out)->npix/2; k+=2) {
00274                 (*out)->data[i][0][(*out)->npix-k+1] = -((*out)->data[i][0][k+1] = ((*out)->data[i][0][(*out)->npix-k]-(*out)->data[i][0][k]) / 2.0);
00275                 (*out)->data[i][0][(*out)->npix-k  ] =  ((*out)->data[i][0][k  ] = ((*out)->data[i][0][(*out)->npix-k]+(*out)->data[i][0][k]) / 2.0);
00276                 (*out)->data[i][(*out)->nlin/2][(*out)->npix-k+1] = -((*out)->data[i][(*out)->nlin/2][k+1] = ((*out)->data[i][(*out)->nlin/2][(*out)->npix-k]-(*out)->data[i][(*out)->nlin/2][k]) / 2.0);
00277                 (*out)->data[i][(*out)->nlin/2][(*out)->npix-k  ] =  ((*out)->data[i][(*out)->nlin/2][k  ] = ((*out)->data[i][(*out)->nlin/2][(*out)->npix-k]+(*out)->data[i][(*out)->nlin/2][k]) / 2.0);
00278             }
00279             for (j=1; j<(*out)->nlin/2; j++) {
00280                 (*out)->data[i][(*out)->nlin-j][1] = -((*out)->data[i][j][1] = ((*out)->data[i][(*out)->nlin-j][0]-(*out)->data[i][j][0]) / 2.0);
00281                 (*out)->data[i][(*out)->nlin-j][0] =  ((*out)->data[i][j][0] = ((*out)->data[i][(*out)->nlin-j][0]+(*out)->data[i][j][0]) / 2.0);
00282                 (*out)->data[i][(*out)->nlin-j][(*out)->npix/2+1] = -((*out)->data[i][j][(*out)->npix/2+1] = ((*out)->data[i][(*out)->nlin-j][(*out)->npix/2]-(*out)->data[i][j][(*out)->npix/2]) / 2.0);
00283                 (*out)->data[i][(*out)->nlin-j][(*out)->npix/2  ] =  ((*out)->data[i][j][(*out)->npix/2  ] = ((*out)->data[i][(*out)->nlin-j][(*out)->npix/2]+(*out)->data[i][j][(*out)->npix/2]) / 2.0);
00284                 for (k=2; k<(*out)->npix/2; k+=2) {
00285                     t = (*out)->data[i][(*out)->nlin-j][(*out)->npix-k];
00286                     u = (*out)->data[i][j][k];
00287                     v = (*out)->data[i][j][(*out)->npix-k];
00288                     w = (*out)->data[i][(*out)->nlin-j][k];
00289                     (*out)->data[i][(*out)->nlin-j][(*out)->npix-k+1] = -((*out)->data[i][j][k+1] = (t-u) / 2.0);
00290                     (*out)->data[i][j][(*out)->npix-k  ] =  ((*out)->data[i][(*out)->nlin-j][k  ] = (t+u) / 2.0);
00291                     (*out)->data[i][j][(*out)->npix-k+1] = -((*out)->data[i][(*out)->nlin-j][k+1] = (v-w) / 2.0);
00292                     (*out)->data[i][(*out)->nlin-j][(*out)->npix-k  ] =  ((*out)->data[i][j][k  ] = (v+w) / 2.0);
00293                 }
00294             }
00295         } else {
00296             /* adjust Hartley coefficients */
00297             for (j=1; j<(*out)->nlin/2; j++) {
00298                 for (k=1; k<(*out)->npix/2; k++) {
00299                     t = ((*out)->data[i][j][k] - (*out)->data[i][(*out)->nlin-j][k] - (*out)->data[i][j][(*out)->npix-k] + (*out)->data[i][(*out)->nlin-j][(*out)->npix-k]) / 2.0;
00300                     (*out)->data[i][j][k]                           -= t;
00301                     (*out)->data[i][(*out)->nlin-j][k]              += t;
00302                     (*out)->data[i][j][(*out)->npix-k]              += t;
00303                     (*out)->data[i][(*out)->nlin-j][(*out)->npix-k] -= t;
00304                 }
00305             }
00306         }
00307  
00308         /* unfold */
00309         for (j=0; j<(*out)->nlin/2; j++) {
00310             for (k=0; k<(*out)->npix/2; k++) {
00311                 t = (*out)->data[i][j][k];
00312                 (*out)->data[i][j][k] = (*out)->data[i][(*out)->nlin/2+j][(*out)->npix/2+k];
00313                 (*out)->data[i][(*out)->nlin/2+j][(*out)->npix/2+k] = t;
00314                 t = (*out)->data[i][(*out)->nlin/2+j][k];
00315                 (*out)->data[i][(*out)->nlin/2+j][k] = (*out)->data[i][j][(*out)->npix/2+k];
00316                 (*out)->data[i][j][(*out)->npix/2+k] = t;
00317             }
00318         }
00319     }
00320 
00321     the_end:
00322     if (a)     free(a);
00323     if (pwr)   free(pwr);
00324     if (sine)  free(sine);
00325     if (tang)  free(tang);
00326     if (array) free(array);
00327     if (TIMES) TIMING(T_EXIT);
00328 }