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
00015 FHT1D (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   )
00027 {
00028   register short o, n, m, l, k, j, i;
00029   register double t, u, v, w, x, y;
00030 
00031   if (TIMES)
00032     TIMING (T_FHT1D);
00033 
00034   a[0] = 0;
00035   a[1] = 1;
00036   for (i = 2; i <= p / 2 + p % 2; i++)
00037     {
00038       for (j = 0; j < pwr[i - 1]; j++)
00039         {
00040           a[j + pwr[i - 1]] = (a[j] *= 2) + 1;
00041         }
00042     }
00043   for (i = 1; i < pwr[p / 2]; i++)
00044     {
00045       k = i;
00046       l = m = pwr[p / 2] * a[i];
00047       t = array[k];
00048       array[k] = array[m];
00049       array[m] = t;
00050       for (j = 1; j <= a[i] - 1; j++)
00051         {
00052           k += pwr[p / 2];
00053           m = l + a[j];
00054           t = array[k];
00055           array[k] = array[m];
00056           array[m] = t;
00057         }
00058     }
00059   for (i = 0; i < size - 1; i += 2)
00060     {
00061       t = array[i] + array[i + 1];
00062       u = array[i] - array[i + 1];
00063       array[i] = t;
00064       array[i + 1] = u;
00065     }
00066   for (i = 0; i < size - 3; i += 4)
00067     {
00068       t = array[i] + array[i + 2];
00069       u = array[i + 1] + array[i + 3];
00070       v = array[i] - array[i + 2];
00071       w = array[i + 1] - array[i + 3];
00072       array[i] = t;
00073       array[i + 1] = u;
00074       array[i + 2] = v;
00075       array[i + 3] = w;
00076     }
00077   for (i = 2; i < p; i++)
00078     {
00079       for (j = 0; j < size; j += pwr[i + 1])
00080         {
00081           t = array[j] + array[j + pwr[i]];
00082           u = array[j] - array[j + pwr[i]];
00083           array[j] = t;
00084           array[j + pwr[i]] = u;
00085           l = j + 1;
00086           m = j + pwr[i] - 1;
00087           n = l + pwr[i];
00088           o = m + pwr[i];
00089           for (k = pwr[p - i - 1]; k <= size / 4; k += pwr[p - i - 1])
00090             {
00091               t = array[n] + tang[k] * array[o];
00092               x = array[o] - sine[k] * t;
00093               y = t + tang[k] * x;
00094               t = array[l] + y;
00095               u = array[m] - x;
00096               v = array[l] - y;
00097               w = array[m] + x;
00098               array[l++] = t;
00099               array[m--] = u;
00100               array[n++] = v;
00101               array[o--] = w;
00102             }
00103         }
00104     }
00105 
00106   if (option == FORWARD)
00107     for (i = 0; i < size; array[i++] /= (double) size);
00108 
00109 the_end:
00110   if (TIMES)
00111     TIMING (T_EXIT);
00112 }
00113 
00114 /*-Copyright Information------------------------------------------------------*/
00115 /* Copyright (c) 1989 by the University of Arizona Digital Image Analysis Lab */
00116 /*----------------------------------------------------------------------------*/
00117 /*-General Information--------------------------------------------------------*/
00118 /*                                                                            */
00119 /*   This function computes the Hartley transform of an image.                */
00120 /*                                                                            */
00121 /*----------------------------------------------------------------------------*/
00122 /*-Background Information-----------------------------------------------------*/
00123 /*                                                                            */
00124 /*   Bracewell, R.N.:                                                         */
00125 /*   "The Hartley Transform."                                                 */
00126 /*   Oxford University Press, New York, 1986                                  */
00127 /*                                                                            */
00128 /*----------------------------------------------------------------------------*/
00129 /*-Interface Information------------------------------------------------------*/
00130 void
00131 FHT2D (IMAGE * in,              /*  I   Pointer to the input image.                           */
00132        short option,            /*  I   Switch, indicating the direction of the transform:    */
00133        /*      FORWARD   -   forward Hartley transform is computed.  */
00134        /*      INVERSE   -   inverse Hartley transform is computed.  */
00135        IMAGE ** out             /*  O   Address of a pointer to the output image.             */
00136 /*----------------------------------------------------------------------------*/
00137   )
00138 {
00139   register short i, j, k, p, *a = 0, *pwr = 0;
00140   char msg[SLEN];
00141   double *sine = 0, *tang = 0, *array = 0, pinc, psum;
00142   PIXEL t, u, v, w;
00143   IMAGE *img;
00144 
00145   if (TIMES)
00146     TIMING (T_FHT2D);
00147   if (NAMES)
00148     {
00149       MESSAGE ('I', "");
00150       MESSAGE ('I', "FHT2D");
00151       MESSAGE ('I', "");
00152       sprintf (msg, " Input image:   %s", in->text);
00153       MESSAGE ('I', msg);
00154       if (option == FORWARD)
00155         {
00156           MESSAGE ('I', " Forward FHT");
00157         }
00158       else if (option == INVERSE)
00159         {
00160           MESSAGE ('I', " Inverse FHT");
00161         }
00162       MESSAGE ('I', " ...............");
00163     }
00164 
00165   /* check input */
00166   if (!CHECKIMG (in))
00167     {
00168       MESSAGE ('E', " Can't identify image.");
00169       goto the_end;
00170     }
00171   else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
00172     {
00173       MESSAGE ('E', " Number of lines must be a power of two.");
00174       goto the_end;
00175     }
00176   else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
00177     {
00178       MESSAGE ('E', " Number of pixels/line must be a power of two.");
00179       goto the_end;
00180     }
00181   else if (option != FORWARD && option != INVERSE)
00182     {
00183       MESSAGE ('E',
00184                " Transformation option must be either FORWARD or INVERSE.");
00185       goto the_end;
00186     }
00187   else if (option == FORWARD && in->nlin != in->npix || option == INVERSE
00188            && in->nlin != in->npix / 2)
00189     {
00190       MESSAGE ('E',
00191                " Number of lines must be equal to number of pixels/line.");
00192       goto the_end;
00193     }
00194 
00195   /* create image of appropriate size */
00196   if (!CHECKIMG (*out))
00197     GETMEM (in->nbnd, in->nlin,
00198             (short) ((option == FORWARD ? 2.0 : 0.5) * in->npix), out);
00199   if (!*out)
00200     goto the_end;
00201   img = option == FORWARD ? in : *out;
00202 
00203   /* special case */
00204   if (img->nlin == 1)
00205     {
00206       (*out)->data[0][0][0] = in->data[0][0][0];
00207       if (option == FORWARD)
00208         {
00209           (*out)->data[0][0][1] = 0.0;
00210         }
00211       goto the_end;
00212     }
00213 
00214   /* allocate permutation index array */
00215   i = max (img->nlin, img->npix);
00216   p = rnd (log10 ((double) i) / log10 (2.0));
00217   if (!(a = (short *) malloc ((1 << p / 2 + p % 2) * sizeof (short))))
00218     {
00219       MESSAGE ('E', " Memory request failed.");
00220       goto the_end;
00221     }
00222 
00223   /* allocate and initialize table with powers of two */
00224   if (!(pwr = (short *) malloc ((p + 1) * sizeof (short))))
00225     {
00226       MESSAGE ('E', " Memory request failed.");
00227       goto the_end;
00228     }
00229   for (j = 0; j <= p; pwr[j] = 1 << j, j++);
00230 
00231   /* allocate and initialize trigonometric tables */
00232   if (!(sine = (double *) malloc ((i / 4 + 1) * sizeof (double))))
00233     {
00234       MESSAGE ('E', " Memory request failed.");
00235       goto the_end;
00236     }
00237   if (!(tang = (double *) malloc ((i / 4 + 1) * sizeof (double))))
00238     {
00239       MESSAGE ('E', " Memory request failed.");
00240       goto the_end;
00241     }
00242   for (j = 0; j <= i / 4; j++)
00243     {
00244       sine[j] = sin ((2.0 * PI * (double) j) / (double) i);
00245       tang[j] = tan ((PI * (double) j) / (double) i);
00246     }
00247 
00248   /* allocate array to hold one row/column of pixels */
00249   if (!(array = (double *) malloc (i * sizeof (double))))
00250     {
00251       MESSAGE ('E', " Memory request failed.");
00252       goto the_end;
00253     }
00254 
00255   /* initialize progress indicator */
00256   if (LINES && !PROGRESS (psum = 0.0))
00257     goto the_end;
00258   pinc = 1.0 / (double) img->nbnd / (double) (img->nlin + img->npix);
00259 
00260   /* Hartley transform bands */
00261   for (i = 0; i < (*out)->nbnd; i++)
00262     {
00263 
00264       if (option == FORWARD)
00265         {
00266           /* fold */
00267           for (j = 0; j < (*out)->nlin / 2; j++)
00268             {
00269               for (k = 0; k < (*out)->npix / 4; k++)
00270                 {
00271                   (*out)->data[i][j][k] =
00272                     in->data[i][in->nlin / 2 + j][in->npix / 2 + k];
00273                   (*out)->data[i][(*out)->nlin / 2 + j][k] =
00274                     in->data[i][j][in->npix / 2 + k];
00275                   (*out)->data[i][j][(*out)->npix / 4 + k] =
00276                     in->data[i][in->nlin / 2 + j][k];
00277                   (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 4 +
00278                                                         k] =
00279                     in->data[i][j][k];
00280                 }
00281             }
00282         }
00283       else
00284         {
00285           /* compute Hartley from Fourier coefficients and fold */
00286           for (j = 0; j < (*out)->nlin / 2; j++)
00287             {
00288               for (k = 0; k < (*out)->npix / 2; k++)
00289                 {
00290                   (*out)->data[i][j][k] =
00291                     in->data[i][in->nlin / 2 + j][in->npix / 2 + 2 * k] -
00292                     in->data[i][in->nlin / 2 + j][in->npix / 2 + 2 * k + 1];
00293                   (*out)->data[i][(*out)->nlin / 2 + j][k] =
00294                     in->data[i][j][in->npix / 2 + 2 * k] -
00295                     in->data[i][j][in->npix / 2 + 2 * k + 1];
00296                   (*out)->data[i][j][(*out)->npix / 2 + k] =
00297                     in->data[i][in->nlin / 2 + j][2 * k] -
00298                     in->data[i][in->nlin / 2 + j][2 * k + 1];
00299                   (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 +
00300                                                         k] =
00301                     in->data[i][j][2 * k] - in->data[i][j][2 * k + 1];
00302                 }
00303             }
00304         }
00305 
00306       /* Hartley transform rows */
00307       for (j = 0; j < img->nlin; j++)
00308         {
00309           for (k = 0; k < img->npix; k++)
00310             {
00311               array[k] = (double) (*out)->data[i][j][k];
00312             }
00313           FHT1D (option, img->npix, p, a, pwr, sine, tang, array);
00314           for (k = 0; k < img->npix; k++)
00315             {
00316               (*out)->data[i][j][k] = (PIXEL) array[k];
00317             }
00318           if (LINES && !PROGRESS (psum += pinc))
00319             goto the_end;
00320         }
00321 
00322       /* Hartley transform columns */
00323       for (k = img->npix - 1; k >= 0; k--)
00324         {
00325           for (j = 0; j < img->nlin; j++)
00326             {
00327               array[j] = (double) (*out)->data[i][j][k];
00328             }
00329           FHT1D (option, img->nlin, p, a, pwr, sine, tang, array);
00330           if (option == FORWARD)
00331             {
00332               for (j = 0; j < img->nlin; j++)
00333                 {
00334                   (*out)->data[i][j][2 * k] = (PIXEL) array[j];
00335                 }
00336             }
00337           else
00338             {
00339               for (j = 0; j < img->nlin; j++)
00340                 {
00341                   (*out)->data[i][j][k] = (PIXEL) array[j];
00342                 }
00343             }
00344           if (LINES && !PROGRESS (psum += pinc))
00345             goto the_end;
00346         }
00347 
00348       if (option == FORWARD)
00349         {
00350           /* compute Fourier from Hartley coefficients */
00351           (*out)->data[i][0][1] = (*out)->data[i][(*out)->nlin / 2][1] =
00352             (*out)->data[i][0][(*out)->npix / 2 + 1] =
00353             (*out)->data[i][(*out)->nlin / 2][(*out)->npix / 2 + 1] = 0.0;
00354           for (k = 2; k < (*out)->npix / 2; k += 2)
00355             {
00356               (*out)->data[i][0][(*out)->npix - k + 1] =
00357                 -((*out)->data[i][0][k + 1] =
00358                   ((*out)->data[i][0][(*out)->npix - k] -
00359                    (*out)->data[i][0][k]) / 2.0);
00360               (*out)->data[i][0][(*out)->npix - k] = ((*out)->data[i][0][k] =
00361                                                       ((*out)->
00362                                                        data[i][0][(*out)->
00363                                                                   npix - k] +
00364                                                        (*out)->
00365                                                        data[i][0][k]) / 2.0);
00366               (*out)->data[i][(*out)->nlin / 2][(*out)->npix - k + 1] =
00367                 -((*out)->data[i][(*out)->nlin / 2][k + 1] =
00368                   ((*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] -
00369                    (*out)->data[i][(*out)->nlin / 2][k]) / 2.0);
00370               (*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] =
00371                 ((*out)->data[i][(*out)->nlin / 2][k] =
00372                  ((*out)->data[i][(*out)->nlin / 2][(*out)->npix - k] +
00373                   (*out)->data[i][(*out)->nlin / 2][k]) / 2.0);
00374             }
00375           for (j = 1; j < (*out)->nlin / 2; j++)
00376             {
00377               (*out)->data[i][(*out)->nlin - j][1] = -((*out)->data[i][j][1] =
00378                                                        ((*out)->
00379                                                         data[i][(*out)->nlin -
00380                                                                 j][0] -
00381                                                         (*out)->
00382                                                         data[i][j][0]) / 2.0);
00383               (*out)->data[i][(*out)->nlin - j][0] = ((*out)->data[i][j][0] =
00384                                                       ((*out)->
00385                                                        data[i][(*out)->nlin -
00386                                                                j][0] +
00387                                                        (*out)->
00388                                                        data[i][j][0]) / 2.0);
00389               (*out)->data[i][(*out)->nlin - j][(*out)->npix / 2 + 1] =
00390                 -((*out)->data[i][j][(*out)->npix / 2 + 1] =
00391                   ((*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] -
00392                    (*out)->data[i][j][(*out)->npix / 2]) / 2.0);
00393               (*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] =
00394                 ((*out)->data[i][j][(*out)->npix / 2] =
00395                  ((*out)->data[i][(*out)->nlin - j][(*out)->npix / 2] +
00396                   (*out)->data[i][j][(*out)->npix / 2]) / 2.0);
00397               for (k = 2; k < (*out)->npix / 2; k += 2)
00398                 {
00399                   t = (*out)->data[i][(*out)->nlin - j][(*out)->npix - k];
00400                   u = (*out)->data[i][j][k];
00401                   v = (*out)->data[i][j][(*out)->npix - k];
00402                   w = (*out)->data[i][(*out)->nlin - j][k];
00403                   (*out)->data[i][(*out)->nlin - j][(*out)->npix - k + 1] =
00404                     -((*out)->data[i][j][k + 1] = (t - u) / 2.0);
00405                   (*out)->data[i][j][(*out)->npix - k] =
00406                     ((*out)->data[i][(*out)->nlin - j][k] = (t + u) / 2.0);
00407                   (*out)->data[i][j][(*out)->npix - k + 1] =
00408                     -((*out)->data[i][(*out)->nlin - j][k + 1] =
00409                       (v - w) / 2.0);
00410                   (*out)->data[i][(*out)->nlin - j][(*out)->npix - k] =
00411                     ((*out)->data[i][j][k] = (v + w) / 2.0);
00412                 }
00413             }
00414         }
00415       else
00416         {
00417           /* adjust Hartley coefficients */
00418           for (j = 1; j < (*out)->nlin / 2; j++)
00419             {
00420               for (k = 1; k < (*out)->npix / 2; k++)
00421                 {
00422                   t =
00423                     ((*out)->data[i][j][k] -
00424                      (*out)->data[i][(*out)->nlin - j][k] -
00425                      (*out)->data[i][j][(*out)->npix - k] +
00426                      (*out)->data[i][(*out)->nlin - j][(*out)->npix -
00427                                                        k]) / 2.0;
00428                   (*out)->data[i][j][k] -= t;
00429                   (*out)->data[i][(*out)->nlin - j][k] += t;
00430                   (*out)->data[i][j][(*out)->npix - k] += t;
00431                   (*out)->data[i][(*out)->nlin - j][(*out)->npix - k] -= t;
00432                 }
00433             }
00434         }
00435 
00436       /* unfold */
00437       for (j = 0; j < (*out)->nlin / 2; j++)
00438         {
00439           for (k = 0; k < (*out)->npix / 2; k++)
00440             {
00441               t = (*out)->data[i][j][k];
00442               (*out)->data[i][j][k] =
00443                 (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 + k];
00444               (*out)->data[i][(*out)->nlin / 2 + j][(*out)->npix / 2 + k] = t;
00445               t = (*out)->data[i][(*out)->nlin / 2 + j][k];
00446               (*out)->data[i][(*out)->nlin / 2 + j][k] =
00447                 (*out)->data[i][j][(*out)->npix / 2 + k];
00448               (*out)->data[i][j][(*out)->npix / 2 + k] = t;
00449             }
00450         }
00451     }
00452 
00453 the_end:
00454   if (a)
00455     free (a);
00456   if (pwr)
00457     free (pwr);
00458   if (sine)
00459     free (sine);
00460   if (tang)
00461     free (tang);
00462   if (array)
00463     free (array);
00464   if (TIMES)
00465     TIMING (T_EXIT);
00466 }