Main Page   Data Structures   File List   Data Fields   Globals  

filter.c

Go to the documentation of this 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 spatially filters an image, using the Fast Fourier         */
00009 /*   transform and a user-supplied frequency filter function. The filter      */
00010 /*   must be real (amplitude only) and bilaterally symmetric.                 */
00011 /*                                                                            */
00012 /*----------------------------------------------------------------------------*/
00013 /*-Background Information-----------------------------------------------------*/
00014 /*                                                                            */
00015 /*   General form of "flt":                                                   */
00016 /*                                                                            */
00017 /*      double flt (double fx, double fy, double xradius, double yradius)     */
00018 /*      {                                                                     */
00019 /*          return (filter amplitude at fx and fy);                           */
00020 /*      }                                                                     */
00021 /*                                                                            */
00022 /*   where fx and fy are spatial frequencies in the range [0.0,0.5]           */
00023 /*   cycles/pixel.  "Flt" is called by this function for all values of        */
00024 /*   fx and fy and should return a filter function which is bilaterally       */
00025 /*   symmetric about the fx and fy axes. "xradius" and "yradius" are the      */
00026 /*   radii of the filter (in cycles/pixel) along the fx and fy frequency      */
00027 /*   axes.  "Flt" may be made radially symmetric by setting "xradius"         */
00028 /*   equal to "yradius".                                                      */
00029 /*                                                                            */
00030 /*----------------------------------------------------------------------------*/
00031 /*-Background Information-----------------------------------------------------*/
00032 /*                                                                            */
00033 /*   Schowengerdt, R.A.:                                                      */
00034 /*   "Techniques for Image Processing and Classification in Remote Sensing."  */
00035 /*   Academic Press, New York, 1983, pp. 75-79                                */
00036 /*                                                                            */
00037 /*----------------------------------------------------------------------------*/
00038 /*-Interface Information------------------------------------------------------*/
00039 void
00040 FILTER (IMAGE * in,             /*  I   Pointer to the input image.                           */
00041         double (*flt) (double, double, double, double),
00042         /*  I   Pointer to a filter function.                         */
00043         double xradius,         /*  I   Radius of filter along fx in cycles/pixel.            */
00044         double yradius,         /*  I   Radius of filter along fy in cycles/line.             */
00045         double option,          /*  I   Switch, indicating the type of filter:                */
00046         /*      = 0.0   -    low-pass   (e.g. filter = flt).          */
00047         /*      = 1.0   -    high-pass  (e.g. filter = 1.0-flt).      */
00048         /*      = 2.0   -    high-boost (e.g. filter = 2.0-flt).      */
00049         /*      else    -                     filter = option-flt.    */
00050         IMAGE ** out            /*  O   Address of a pointer to the output image.             */
00051         /*      The input image may be overwritten (out = &in).       */
00052 /*----------------------------------------------------------------------------*/
00053   )
00054 {
00055   register short i, j, k;
00056   char msg[SLEN];
00057   double fx, fy, ampl, pinc, psum;
00058   IMAGE *scratch = 0;
00059 
00060   if (TIMES)
00061     TIMING (T_FILTER);
00062   if (NAMES)
00063     {
00064       MESSAGE ('I', "");
00065       MESSAGE ('I', "FILTER");
00066       MESSAGE ('I', "");
00067       sprintf (msg, " Input image:                   %s", in->text);
00068       MESSAGE ('I', msg);
00069       if (flt == FLTBOX)
00070         {
00071           MESSAGE ('I', " Box filter");
00072         }
00073       else if (flt == FLTCONE)
00074         {
00075           MESSAGE ('I', " Cone filter");
00076         }
00077       else if (flt == FLTCYL)
00078         {
00079           MESSAGE ('I', " Cylinder filter");
00080         }
00081       else if (flt == FLTEXP)
00082         {
00083           MESSAGE ('I', " Exponential filter");
00084         }
00085       else if (flt == FLTGAUSS)
00086         {
00087           MESSAGE ('I', " Gaussian filter");
00088         }
00089       sprintf (msg, " Filter radius: cycles/line:  %8.4f", yradius);
00090       MESSAGE ('I', msg);
00091       sprintf (msg, "                cycles/pixel: %8.4f", xradius);
00092       MESSAGE ('I', msg);
00093       if (option == 0.0)
00094         {
00095           MESSAGE ('I', " Low-pass filter");
00096         }
00097       else if (option == 1.0)
00098         {
00099           MESSAGE ('I', " High-pass filter");
00100         }
00101       else
00102         {
00103           sprintf (msg, " Parametric high-boost filter, c = %5.2f", option);
00104           MESSAGE ('I', msg);
00105         }
00106       MESSAGE ('I', " ...............");
00107     }
00108 
00109   /* check input */
00110   if (!CHECKIMG (in))
00111     {
00112       MESSAGE ('E', " Can't identify image.");
00113       goto the_end;
00114     }
00115   else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
00116     {
00117       MESSAGE ('E', " Number of lines must be a power of two.");
00118       goto the_end;
00119     }
00120   else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
00121     {
00122       MESSAGE ('E', " Number of pixels/line must be a power of two.");
00123       goto the_end;
00124     }
00125   else if (!flt)
00126     {
00127       MESSAGE ('E', " Can't identify filter function.");
00128       goto the_end;
00129     }
00130 
00131   /* Fourier transform image */
00132   FFT2D (in, FORWARD, &scratch);
00133   if (ABORT || !scratch)
00134     goto the_end;
00135   sprintf (scratch->text, "Scratch");
00136 
00137   /* initialize progress indicator */
00138   if (LINES && !PROGRESS (psum = 0.0))
00139     goto the_end;
00140   pinc = 1.0 / (double) scratch->nbnd / (double) (scratch->nlin / 2 - 1);
00141 
00142   /* filter image */
00143   for (i = 0; i < scratch->nbnd; i++)
00144     {
00145       for (j = 1; j < scratch->nlin / 2; j++)
00146         {
00147           fy = 0.5 - (double) j / (double) scratch->nlin;
00148           for (k = 2; k < scratch->npix / 2 - 1; k += 2)
00149             {
00150               fx = 0.5 - (double) k / (double) (scratch->npix);
00151               ampl = (*flt) (fx, fy, xradius, yradius);
00152               if (option != 0.0)
00153                 ampl = option - ampl;
00154 
00155               /* real components */
00156               scratch->data[i][j][k] *= ampl;
00157               scratch->data[i][j][scratch->npix - k] *= ampl;
00158 
00159               /*  conjugate symmetry */
00160               scratch->data[i][scratch->nlin - j][k] =
00161                 scratch->data[i][j][scratch->npix - k];
00162               scratch->data[i][scratch->nlin - j][scratch->npix - k] =
00163                 scratch->data[i][j][k];
00164 
00165               /* imaginary components */
00166               scratch->data[i][j][k + 1] *= ampl;
00167               scratch->data[i][j][scratch->npix - k + 1] *= ampl;
00168 
00169               /* conjugate symmetry */
00170               scratch->data[i][scratch->nlin - j][k + 1] =
00171                 -scratch->data[i][j][scratch->npix - k + 1];
00172               scratch->data[i][scratch->nlin - j][scratch->npix - k + 1] =
00173                 -scratch->data[i][j][k + 1];
00174             }
00175           if (LINES && !PROGRESS (psum += pinc))
00176             goto the_end;
00177         }
00178     }
00179 
00180   /* zero frequencies */
00181   for (i = 0; i < scratch->nbnd; i++)
00182     {
00183       for (j = 1; j < scratch->nlin; j++)
00184         {
00185           fx = 0.0;
00186           fy = 0.5 - (double) j / (double) scratch->nlin;
00187           ampl = (*flt) (fx, fy, xradius, yradius);
00188           if (option != 0.0)
00189             ampl = option - ampl;
00190           scratch->data[i][j][scratch->npix / 2] *= ampl;
00191           scratch->data[i][j][scratch->npix / 2 + 1] *= ampl;
00192         }
00193     }
00194   for (i = 0; i < scratch->nbnd; i++)
00195     {
00196       for (k = 2; k < scratch->npix; k += 2)
00197         {
00198           fx = 0.5 - (double) k / (double) scratch->npix;
00199           fy = 0.0;
00200           ampl = (*flt) (fx, fy, xradius, yradius);
00201           if (option != 0.0)
00202             ampl = option - ampl;
00203           scratch->data[i][scratch->nlin / 2][k] *= ampl;
00204           scratch->data[i][scratch->nlin / 2][k + 1] *= ampl;
00205         }
00206     }
00207 
00208   /* folding frequencies */
00209   for (i = 0; i < scratch->nbnd; i++)
00210     {
00211       for (j = 0; j < scratch->nlin; j++)
00212         {
00213           fx = 0.5;
00214           fy = 0.5 - (double) j / (double) scratch->nlin;
00215           ampl = (*flt) (fx, fy, xradius, yradius);
00216           if (option != 0.0)
00217             ampl = option - ampl;
00218           scratch->data[i][j][0] *= ampl;
00219           scratch->data[i][j][1] *= ampl;
00220         }
00221     }
00222   for (i = 0; i < scratch->nbnd; i++)
00223     {
00224       for (k = 2; k < scratch->npix; k += 2)
00225         {
00226           fx = 0.5 - (double) k / (double) scratch->npix;
00227           fy = 0.5;
00228           ampl = (*flt) (fx, fy, xradius, yradius);
00229           if (option != 0.0)
00230             ampl = option - ampl;
00231           scratch->data[i][0][k] *= ampl;
00232           scratch->data[i][0][k + 1] *= ampl;
00233         }
00234     }
00235 
00236   /* Fourier transform image */
00237   FFT2D (scratch, INVERSE, out);
00238 
00239 the_end:
00240   if (scratch)
00241     RELMEM (scratch);
00242   if (TIMES)
00243     TIMING (T_EXIT);
00244 }
00245 
00246 /*-Copyright Information------------------------------------------------------*/
00247 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00248 /*----------------------------------------------------------------------------*/
00249 /*-General Information--------------------------------------------------------*/
00250 /*                                                                            */
00251 /*   This library function computes a box spatial frequency filter            */
00252 /*                                                                            */
00253 /*----------------------------------------------------------------------------*/
00254 /*-Interface Information------------------------------------------------------*/
00255 double
00256 FLTBOX (double fx,              /*  I   Spatial frequency.                                    */
00257         double fy,              /*  I   Spatial frequency.                                    */
00258         double xradius,         /*  I   Radius of filter along fx.                            */
00259         double yradius          /*  I   Radius of filter along fy.                            */
00260 /*----------------------------------------------------------------------------*/
00261   )
00262 {
00263   double ampl;
00264 
00265   if (TIMES)
00266     TIMING (T_FLTBOX);
00267 
00268   ampl = fabs (fx / xradius) <= 1.0 && fabs (fy / yradius) <= 1.0 ? 1.0 : 0.0;
00269 
00270 the_end:
00271   if (TIMES)
00272     TIMING (T_EXIT);
00273   return (ampl);
00274 }
00275 
00276 /*-Copyright Information------------------------------------------------------*/
00277 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00278 /*----------------------------------------------------------------------------*/
00279 /*-General Information--------------------------------------------------------*/
00280 /*                                                                            */
00281 /*   This library function computes a cone spatial frequency filter.          */
00282 /*                                                                            */
00283 /*----------------------------------------------------------------------------*/
00284 /*-Interface Information------------------------------------------------------*/
00285 double
00286 FLTCONE (double fx,             /*  I   Spatial frequency.                                    */
00287          double fy,             /*  I   Spatial frequency.                                    */
00288          double xradius,        /*  I   Radius of filter along fx.                            */
00289          double yradius         /*  I   Radius of filter along fy.                            */
00290 /*----------------------------------------------------------------------------*/
00291   )
00292 {
00293   double ampl;
00294 
00295   if (TIMES)
00296     TIMING (T_FLTCONE);
00297 
00298   ampl =
00299     max (0.0,
00300          1.0 - sqrt ((fx * fx) / (xradius * xradius) +
00301                      (fy * fy) / (yradius * yradius)));
00302 
00303 the_end:
00304   if (TIMES)
00305     TIMING (T_EXIT);
00306   return (ampl);
00307 }
00308 
00309 /*-Copyright Information------------------------------------------------------*/
00310 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00311 /*----------------------------------------------------------------------------*/
00312 /*-General Information--------------------------------------------------------*/
00313 /*                                                                            */
00314 /*   This library function computes a cylinder spatial frequency filter       */
00315 /*                                                                            */
00316 /*----------------------------------------------------------------------------*/
00317 /*-Interface Information------------------------------------------------------*/
00318 double
00319 FLTCYL (double fx,              /*  I   Spatial frequency.                                    */
00320         double fy,              /*  I   Spatial frequency.                                    */
00321         double xradius,         /*  I   Radius of filter along fx.                            */
00322         double yradius          /*  I   Radius of filter along fy.                            */
00323 /*----------------------------------------------------------------------------*/
00324   )
00325 {
00326   double ampl;
00327 
00328   if (TIMES)
00329     TIMING (T_FLTCYL);
00330 
00331   ampl =
00332     sqrt ((fx * fx) / (xradius * xradius) +
00333           (fy * fy) / (yradius * yradius)) <= 1.0 ? 1.0 : 0.0;
00334 
00335 the_end:
00336   if (TIMES)
00337     TIMING (T_EXIT);
00338   return (ampl);
00339 }
00340 
00341 /*-Copyright Information------------------------------------------------------*/
00342 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00343 /*----------------------------------------------------------------------------*/
00344 /*-General Information--------------------------------------------------------*/
00345 /*                                                                            */
00346 /*   This library function computes a negative exponential spatial frequency  */
00347 /*   filter.                                                                  */
00348 /*                                                                            */
00349 /*----------------------------------------------------------------------------*/
00350 /*-Interface Information------------------------------------------------------*/
00351 double
00352 FLTEXP (double fx,              /*  I   Spatial frequency.                                    */
00353         double fy,              /*  I   Spatial frequency.                                    */
00354         double xradius,         /*  I   Radius of filter along fx.                            */
00355         double yradius          /*  I   Radius of filter along fy.                            */
00356 /*----------------------------------------------------------------------------*/
00357   )
00358 {
00359   double ampl;
00360 
00361   if (TIMES)
00362     TIMING (T_FLTEXP);
00363 
00364   ampl =
00365     exp (-sqrt
00366          ((fx * fx) / (xradius * xradius) + (fy * fy) / (yradius * yradius)));
00367 
00368 the_end:
00369   if (TIMES)
00370     TIMING (T_EXIT);
00371   return (ampl);
00372 }
00373 
00374 /*-Copyright Information------------------------------------------------------*/
00375 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00376 /*----------------------------------------------------------------------------*/
00377 /*-General Information--------------------------------------------------------*/
00378 /*                                                                            */
00379 /*   This library function computes a gaussian spatial frequency filter.      */
00380 /*                                                                            */
00381 /*----------------------------------------------------------------------------*/
00382 /*-Interface Information------------------------------------------------------*/
00383 double
00384 FLTGAUSS (double fx,            /*  I   Spatial frequency.                                    */
00385           double fy,            /*  I   Spatial frequency.                                    */
00386           double xradius,       /*  I   Radius of filter along fx.                            */
00387           double yradius        /*  I   Radius of filter along fy.                            */
00388 /*----------------------------------------------------------------------------*/
00389   )
00390 {
00391   double ampl;
00392 
00393   if (TIMES)
00394     TIMING (T_FLTGAUSS);
00395 
00396   ampl =
00397     exp (-
00398          ((fx * fx) / (xradius * xradius) + (fy * fy) / (yradius * yradius)));
00399 
00400 the_end:
00401   if (TIMES)
00402     TIMING (T_EXIT);
00403   return (ampl);
00404 }

Generated on Sun May 18 15:36:09 2003 for tclSadie by doxygen1.3