filter.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 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 FILTER (
00040 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 ) { register short i, j, k;
00054     char   msg[SLEN];
00055     double fx, fy, ampl, pinc, psum;
00056     IMAGE  *scratch=0;
00057 
00058     if (TIMES) TIMING(T_FILTER);
00059     if (NAMES) {
00060         MESSAGE('I',"");
00061         MESSAGE('I',"FILTER");
00062         MESSAGE('I',"");
00063         sprintf(msg," Input image:                   %s",in->text);
00064         MESSAGE('I',msg);
00065         if (flt == FLTBOX) {
00066             MESSAGE('I'," Box filter");
00067         } else if (flt == FLTCONE) {
00068             MESSAGE('I'," Cone filter");
00069         } else if (flt == FLTCYL) {
00070             MESSAGE('I'," Cylinder filter");
00071         } else if (flt == FLTEXP) {
00072             MESSAGE('I'," Exponential filter");
00073         } else if (flt == FLTGAUSS) {
00074             MESSAGE('I'," Gaussian filter");
00075         }
00076         sprintf(msg," Filter radius: cycles/line:  %8.4f",yradius);
00077         MESSAGE('I',msg);
00078         sprintf(msg,"                cycles/pixel: %8.4f",xradius);
00079         MESSAGE('I',msg);
00080         if (option == 0.0) {
00081             MESSAGE('I'," Low-pass filter");
00082         } else if (option == 1.0) {
00083             MESSAGE('I'," High-pass filter");
00084         } else {
00085             sprintf(msg," Parametric high-boost filter, c = %5.2f",option);
00086             MESSAGE('I',msg);
00087         }
00088         MESSAGE('I'," ...............");
00089     }
00090 
00091     /* check input */
00092     if (!CHECKIMG(in)) {
00093         MESSAGE('E'," Can't identify image.");
00094         goto the_end;
00095     } else if (in->nlin != 1L<<rnd(log10((double)in->nlin)/log10(2.0))) {
00096         MESSAGE('E'," Number of lines must be a power of two.");
00097         goto the_end;
00098     } else if (in->npix != 1L<<rnd(log10((double)in->npix)/log10(2.0))) {
00099         MESSAGE('E'," Number of pixels/line must be a power of two.");
00100         goto the_end;
00101     } else if (!flt) {
00102         MESSAGE('E'," Can't identify filter function.");
00103         goto the_end;
00104     }
00105 
00106     /* Fourier transform image */
00107     FFT2D(in,FORWARD,&scratch);
00108     if (ABORT  ||  !scratch) goto the_end;
00109     sprintf(scratch->text,"Scratch");
00110 
00111     /* initialize progress indicator */
00112     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00113     pinc = 1.0/(double)scratch->nbnd/(double)(scratch->nlin/2-1);
00114 
00115     /* filter image */
00116     for (i=0; i<scratch->nbnd; i++) {
00117         for (j=1; j<scratch->nlin/2; j++) {
00118             fy = 0.5-(double)j/(double)scratch->nlin;
00119             for (k=2; k<scratch->npix/2-1; k+=2) {
00120                 fx = 0.5-(double)k/(double)(scratch->npix);
00121                 ampl = (*flt)(fx,fy,xradius,yradius);
00122                 if (option != 0.0) ampl = option - ampl;
00123 
00124                 /* real components */
00125                 scratch->data[i][j][k]               *= ampl;
00126                 scratch->data[i][j][scratch->npix-k] *= ampl;
00127 
00128                 /*  conjugate symmetry */
00129                 scratch->data[i][scratch->nlin-j][k] = scratch->data[i][j][scratch->npix-k];
00130                 scratch->data[i][scratch->nlin-j][scratch->npix-k] = scratch->data[i][j][k];
00131 
00132                 /* imaginary components */
00133                 scratch->data[i][j][k+1]               *= ampl;
00134                 scratch->data[i][j][scratch->npix-k+1] *= ampl;
00135 
00136                 /* conjugate symmetry */
00137                 scratch->data[i][scratch->nlin-j][k+1] = -scratch->data[i][j][scratch->npix-k+1];
00138                 scratch->data[i][scratch->nlin-j][scratch->npix-k+1] = -scratch->data[i][j][k+1];
00139             }
00140             if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
00141         }
00142     }
00143 
00144     /* zero frequencies */
00145     for (i=0; i<scratch->nbnd; i++) {
00146         for (j=1; j<scratch->nlin; j++) {
00147             fx = 0.0;
00148             fy = 0.5 - (double)j/(double)scratch->nlin;
00149             ampl = (*flt)(fx,fy,xradius,yradius);
00150             if (option != 0.0) ampl = option - ampl;
00151             scratch->data[i][j][scratch->npix/2]   *= ampl;
00152             scratch->data[i][j][scratch->npix/2+1] *= ampl;
00153         }
00154     }
00155     for (i=0; i<scratch->nbnd; i++) {
00156         for (k=2; k<scratch->npix; k+=2) {
00157             fx = 0.5 - (double)k/(double)scratch->npix;
00158             fy = 0.0;
00159             ampl = (*flt)(fx,fy,xradius,yradius);
00160             if (option != 0.0) ampl = option - ampl;
00161             scratch->data[i][scratch->nlin/2][k]   *= ampl;
00162             scratch->data[i][scratch->nlin/2][k+1] *= ampl;
00163         }
00164     }
00165 
00166     /* folding frequencies */
00167     for (i=0; i<scratch->nbnd; i++) {
00168         for (j=0; j<scratch->nlin; j++) {
00169             fx = 0.5;
00170             fy = 0.5-(double)j/(double)scratch->nlin;
00171             ampl = (*flt)(fx,fy,xradius,yradius);
00172             if (option != 0.0) ampl = option - ampl;
00173             scratch->data[i][j][0] *= ampl;
00174             scratch->data[i][j][1] *= ampl;
00175         }
00176     }
00177     for (i=0; i<scratch->nbnd; i++) {
00178         for (k=2; k<scratch->npix; k+=2) {
00179             fx = 0.5-(double)k/(double)scratch->npix;
00180             fy = 0.5;
00181             ampl = (*flt)(fx,fy,xradius,yradius);
00182             if (option != 0.0) ampl = option - ampl;
00183             scratch->data[i][0][k]   *=ampl;
00184             scratch->data[i][0][k+1] *= ampl;
00185         }
00186     }
00187 
00188     /* Fourier transform image */
00189     FFT2D(scratch,INVERSE,out);
00190 
00191     the_end:
00192     if (scratch) RELMEM(scratch);
00193     if (TIMES)   TIMING(T_EXIT);
00194 }
00195 
00196 /*-Copyright Information------------------------------------------------------*/
00197 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00198 /*----------------------------------------------------------------------------*/
00199 /*-General Information--------------------------------------------------------*/
00200 /*                                                                            */
00201 /*   This library function computes a box spatial frequency filter            */
00202 /*                                                                            */
00203 /*----------------------------------------------------------------------------*/
00204 /*-Interface Information------------------------------------------------------*/
00205 double FLTBOX (
00206 double fx,      /*  I   Spatial frequency.                                    */
00207 double fy,      /*  I   Spatial frequency.                                    */
00208 double xradius, /*  I   Radius of filter along fx.                            */
00209 double yradius  /*  I   Radius of filter along fy.                            */
00210 /*----------------------------------------------------------------------------*/
00211 ) { double ampl;
00212 
00213     if (TIMES) TIMING(T_FLTBOX);
00214 
00215     ampl = fabs(fx/xradius) <= 1.0  &&  fabs(fy/yradius) <= 1.0  ?  1.0 : 0.0;
00216 
00217     the_end:
00218     if (TIMES) TIMING(T_EXIT);
00219     return(ampl);
00220 }
00221 
00222 /*-Copyright Information------------------------------------------------------*/
00223 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00224 /*----------------------------------------------------------------------------*/
00225 /*-General Information--------------------------------------------------------*/
00226 /*                                                                            */
00227 /*   This library function computes a cone spatial frequency filter.          */
00228 /*                                                                            */
00229 /*----------------------------------------------------------------------------*/
00230 /*-Interface Information------------------------------------------------------*/
00231 double FLTCONE (
00232 double fx,      /*  I   Spatial frequency.                                    */
00233 double fy,      /*  I   Spatial frequency.                                    */
00234 double xradius, /*  I   Radius of filter along fx.                            */
00235 double yradius  /*  I   Radius of filter along fy.                            */
00236 /*----------------------------------------------------------------------------*/
00237 ) { double ampl;
00238 
00239     if (TIMES) TIMING(T_FLTCONE);
00240 
00241     ampl = max(0.0,1.0-sqrt((fx*fx)/(xradius*xradius)+(fy*fy)/(yradius*yradius)));
00242 
00243     the_end:
00244     if (TIMES) TIMING(T_EXIT);
00245     return(ampl);
00246 }
00247 
00248 /*-Copyright Information------------------------------------------------------*/
00249 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00250 /*----------------------------------------------------------------------------*/
00251 /*-General Information--------------------------------------------------------*/
00252 /*                                                                            */
00253 /*   This library function computes a cylinder spatial frequency filter       */
00254 /*                                                                            */
00255 /*----------------------------------------------------------------------------*/
00256 /*-Interface Information------------------------------------------------------*/
00257 double FLTCYL (
00258 double fx,      /*  I   Spatial frequency.                                    */
00259 double fy,      /*  I   Spatial frequency.                                    */
00260 double xradius, /*  I   Radius of filter along fx.                            */
00261 double yradius  /*  I   Radius of filter along fy.                            */
00262 /*----------------------------------------------------------------------------*/
00263 ) { double ampl;
00264 
00265     if (TIMES) TIMING(T_FLTCYL);
00266 
00267     ampl = sqrt((fx*fx)/(xradius*xradius)+(fy*fy)/(yradius*yradius)) <= 1.0  ?  1.0 : 0.0;
00268 
00269     the_end:
00270     if (TIMES) TIMING(T_EXIT);
00271     return(ampl);
00272 }
00273 
00274 /*-Copyright Information------------------------------------------------------*/
00275 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00276 /*----------------------------------------------------------------------------*/
00277 /*-General Information--------------------------------------------------------*/
00278 /*                                                                            */
00279 /*   This library function computes a negative exponential spatial frequency  */
00280 /*   filter.                                                                  */
00281 /*                                                                            */
00282 /*----------------------------------------------------------------------------*/
00283 /*-Interface Information------------------------------------------------------*/
00284 double FLTEXP (
00285 double fx,      /*  I   Spatial frequency.                                    */
00286 double fy,      /*  I   Spatial frequency.                                    */
00287 double xradius, /*  I   Radius of filter along fx.                            */
00288 double yradius  /*  I   Radius of filter along fy.                            */
00289 /*----------------------------------------------------------------------------*/
00290 ) { double ampl;
00291 
00292     if (TIMES) TIMING(T_FLTEXP);
00293 
00294     ampl = exp(-sqrt((fx*fx)/(xradius*xradius)+(fy*fy)/(yradius*yradius)));
00295 
00296     the_end:
00297     if (TIMES) TIMING(T_EXIT);
00298     return(ampl);
00299 }
00300 
00301 /*-Copyright Information------------------------------------------------------*/
00302 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00303 /*----------------------------------------------------------------------------*/
00304 /*-General Information--------------------------------------------------------*/
00305 /*                                                                            */
00306 /*   This library function computes a gaussian spatial frequency filter.      */
00307 /*                                                                            */
00308 /*----------------------------------------------------------------------------*/
00309 /*-Interface Information------------------------------------------------------*/
00310 double FLTGAUSS (
00311 double fx,      /*  I   Spatial frequency.                                    */
00312 double fy,      /*  I   Spatial frequency.                                    */
00313 double xradius, /*  I   Radius of filter along fx.                            */
00314 double yradius  /*  I   Radius of filter along fy.                            */
00315 /*----------------------------------------------------------------------------*/
00316 ) { double ampl;
00317 
00318     if (TIMES) TIMING(T_FLTGAUSS);
00319 
00320     ampl = exp(-((fx*fx)/(xradius*xradius)+(fy*fy)/(yradius*yradius)));
00321 
00322     the_end:
00323     if (TIMES) TIMING(T_EXIT);
00324     return(ampl);
00325 }