window.c Source File

00001 #include        "sadie.h"
00002 
00003 extern short  csize;
00004 extern PIXEL  thresh;
00005 extern double weight;
00006 
00007 /*-Copyright Information------------------------------------------------------*/
00008 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00009 /*----------------------------------------------------------------------------*/
00010 /*-General Information--------------------------------------------------------*/
00011 /*                                                                            */
00012 /*   This function performs a user-specified neighborhood operation using     */
00013 /*   a moving window. Several library operations (WTF...) are available.      */
00014 /*                                                                            */
00015 /*----------------------------------------------------------------------------*/
00016 /*-Interface Information------------------------------------------------------*/
00017 void WINDOW (
00018 IMAGE *in,      /*  I   Pointer to the input image.                           */
00019 PIXEL (*wtf)(PIXEL *, short, short),
00020                 /*  I   Address of the filter function.                       */
00021 short jsize,    /*  I   Number of lines in the window.                        */
00022 short ksize,    /*  I   Number of pixels/line in the window.                  */
00023 IMAGE **out     /*  O   Address of a pointer to the output image.             */
00024 /*----------------------------------------------------------------------------*/
00025 ) { register short i, j, k, l, m, n, jmarg=(jsize-1)/2, kmarg=(ksize-1)/2;
00026     char   msg[SLEN];
00027     PIXEL  *wnd=0;
00028     double pinc, psum;
00029 
00030     if (TIMES) TIMING(T_WINDOW);
00031     if (NAMES) {
00032         MESSAGE('I',"");
00033         MESSAGE('I',"WINDOW");
00034         MESSAGE('I',"");
00035         sprintf(msg," Input image:           %s",in->text);
00036         MESSAGE('I',msg);
00037         sprintf(msg," Window size: lines:    %d",jsize);
00038         MESSAGE('I',msg);
00039         sprintf(msg,"              pixels:   %d",ksize);
00040         MESSAGE('I',msg);
00041         if (wtf == WTFCAVG) {
00042             MESSAGE('I'," Conditional average filter");
00043             sprintf(msg," Set size threshold:               %d",csize);
00044             MESSAGE('I',msg);
00045             sprintf(msg," Graylevel difference threshold: %12.4e",thresh);
00046             MESSAGE('I',msg);
00047         } else if (wtf == WTFDIVER) {
00048             MESSAGE('I'," Diversity map");
00049         } else if (wtf == WTFMAJ) {
00050             MESSAGE('I'," Majority filter");
00051             sprintf(msg," Center pixel weight: %12.4e",weight);
00052             MESSAGE('I',msg);
00053         } else if (wtf == WTFMAX) {
00054             MESSAGE('I'," Maximum filter");
00055         } else if (wtf == WTFMIN) {
00056             MESSAGE('I'," Minimum filter");
00057         } else if (wtf == WTFVAR) {
00058             MESSAGE('I'," Variance map");
00059         }
00060         MESSAGE('I'," ...............");
00061     }
00062 
00063     /* check input */
00064     if (!CHECKIMG(in)) {
00065         MESSAGE('E'," Can't identify image.");
00066         goto the_end;
00067     } else if (!wtf) {
00068         MESSAGE('E'," Can't identify window function.");
00069         goto the_end;
00070     } else if (jsize <= 0  ||  jsize%2 == 0) {
00071         MESSAGE('E'," Window line size must be an odd number greater than zero.");
00072         goto the_end;
00073     } else if (jsize > in->nlin) {
00074         MESSAGE('E'," Image line size must be equal to or greater than window line size.");
00075         goto the_end;
00076     } else if (ksize <= 0  ||  ksize%2 == 0) {
00077         MESSAGE('E'," Window pixel size must be an odd number greater than zero.");
00078         goto the_end;
00079     } else if (ksize > in->npix) {
00080         MESSAGE('E'," Image pixel size must be equal to or greater than window pixel size.");
00081         goto the_end;
00082     }
00083 
00084     /* create image of appropriate size */
00085     if (!CHECKIMG(*out)) GETMEM(in->nbnd,in->nlin,in->npix,out);
00086     if (!*out) goto the_end;
00087 
00088     /* allocate space for 1-D array */
00089     if (!(wnd=(PIXEL *)malloc(jsize*ksize*sizeof(PIXEL)))) {
00090         MESSAGE('E'," Memory request failed.");
00091         goto the_end;
00092     }
00093 
00094     /* initialize progress indicator */
00095     if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
00096     pinc = 1.0/(double)in->nbnd/(double)(in->nlin-2*jmarg);
00097 
00098     /* put pixels inside window into array and call the window function */
00099     for (i=0; i<in->nbnd; i++) {
00100         for (j=jmarg; j<in->nlin-jmarg; j++) {
00101             for (k=kmarg; k<in->npix-kmarg; k++) {
00102                 for (n=0,l=j-jmarg; l<=j+jmarg; l++) {
00103                     for (m=k-kmarg; m<=k+kmarg; m++) {
00104                         wnd[n++] = in->data[i][l][m];
00105                     }
00106                 }
00107                 (*out)->data[i][j][k] = (*wtf)(wnd,jsize,ksize);
00108             }
00109             if (LINES  && !PROGRESS(psum+=pinc)) goto the_end;
00110         }
00111 
00112         /* fill top and bottom margins */
00113         for (j=0; j<jmarg; j++) {
00114             for (k=kmarg; k<in->npix-kmarg; k++) {
00115                 (*out)->data[i][j][k]            = (*out)->data[i][jmarg][k];
00116                 (*out)->data[i][in->nlin-j-1][k] = (*out)->data[i][in->nlin-jmarg-1][k];
00117             }
00118         }
00119 
00120         /* fill left and right margins */
00121         for (j=0; j<in->nlin; j++) {
00122             for (k=0; k<kmarg; k++) {
00123                 (*out)->data[i][j][k]            = (*out)->data[i][j][kmarg];
00124                 (*out)->data[i][j][in->npix-k-1] = (*out)->data[i][j][in->npix-kmarg-1];
00125             }
00126         }
00127     }
00128 
00129     the_end:
00130     if (wnd)   free(wnd);
00131     if (TIMES) TIMING(T_EXIT);
00132 }
00133  
00134 /*-Copyright Information------------------------------------------------------*/
00135 /* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
00136 /*----------------------------------------------------------------------------*/
00137 /*-General Information--------------------------------------------------------*/
00138 /*                                                                            */
00139 /*   This library function computes a local conditional average filter.       */
00140 /*                                                                            */
00141 /*----------------------------------------------------------------------------*/
00142 /*-Background Information-----------------------------------------------------*/
00143 /*                                                                            */
00144 /*   Ryan, T.W.:                                                              */
00145 /*   "Image Segmentation Algorithms."                                         */
00146 /*   SPIE, Vol. 534, pp. 172-178, 1985                                        */
00147 /*                                                                            */
00148 /*----------------------------------------------------------------------------*/
00149 /*-Interface Information------------------------------------------------------*/
00150 extern short csize;     /*  I   Threshold for size of center pixel set.       */
00151 extern PIXEL thresh;    /*  I   Graylevel difference threshold.               */
00152 PIXEL WTFCAVG (
00153 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00154 short jsize,    /*  I   Number of lines in the window.                        */
00155 short ksize     /*  I   Number of pixels/line in the window.                  */
00156 /*----------------------------------------------------------------------------*/
00157 ) { register short i, count=0;
00158     PIXEL center=wnd[jsize/2*ksize+ksize/2], sum1=(PIXEL)0, sum2=(PIXEL)0, minval=wnd[0], maxval=wnd[0];
00159 
00160     if (TIMES) TIMING(T_WTFCAVG);
00161 
00162     /* compute conditional average */
00163     for (i=0; i<jsize*ksize; i++) {
00164         if ((PIXEL)fabs((double)(wnd[i]-center)) <= thresh) {
00165             sum1  += wnd[i];
00166             count += 1;
00167         } else if (csize != 0) {
00168             sum2  += wnd[i];
00169             minval = min(minval,wnd[i]);
00170             maxval = max(maxval,wnd[i]);
00171         }
00172     }
00173     if (csize == 0  ||  count > csize) {
00174         center = (PIXEL)((double)sum1/(double)count);
00175     } else if (csize != 0  &&  (PIXEL)fabs((double)(maxval-minval)) <= thresh) {
00176         center = (PIXEL)((double)sum2/(double)(jsize*ksize-count));
00177     }
00178 
00179     the_end:
00180     if (TIMES) TIMING(T_EXIT);
00181     return(center);
00182 }
00183  
00184 /*-Copyright Information------------------------------------------------------*/
00185 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00186 /*----------------------------------------------------------------------------*/
00187 /*-General Information--------------------------------------------------------*/
00188 /*                                                                            */
00189 /*   This library function computes the local graylevel diversity.            */
00190 /*                                                                            */
00191 /*----------------------------------------------------------------------------*/
00192 /*-Interface Information------------------------------------------------------*/
00193 PIXEL WTFDIVER (
00194 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00195 short jsize,    /*  I   Number of lines in the window.                        */
00196 short ksize     /*  I   Number of pixels/line in the window.                  */
00197 /*----------------------------------------------------------------------------*/
00198 ) { register short i, j, diversity=0;
00199 
00200     if (TIMES) TIMING(T_WTFDIVER);
00201 
00202     for (i=0; i<jsize*ksize; i++) {
00203         for (++diversity,j=0; j<i; j++) {
00204             if (wnd[i] == wnd[j]) {
00205                 diversity -= 1;
00206                 break;
00207             }
00208         }
00209     }
00210 
00211     the_end:
00212     if (TIMES) TIMING(T_EXIT);
00213     return((PIXEL)diversity);
00214 }
00215 
00216 /*-Copyright Information------------------------------------------------------*/
00217 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00218 /*----------------------------------------------------------------------------*/
00219 /*-General Information--------------------------------------------------------*/
00220 /*                                                                            */
00221 /*   This library function computes a local majority filter.                  */
00222 /*                                                                            */
00223 /*----------------------------------------------------------------------------*/
00224 /*-Background Information-----------------------------------------------------*/
00225 /*                                                                            */
00226 /*   Schowengerdt, R.A.:                                                      */
00227 /*   "Techniques for Image Processing and Classification in Remote Sensing."  */
00228 /*   Academic Press, New York, 1983, pp. 188-190                              */
00229 /*                                                                            */
00230 /*----------------------------------------------------------------------------*/
00231 /*-Interface Information------------------------------------------------------*/
00232 extern short  nlev;     /*  I   Number of graylevels in the input image.      */
00233 extern PIXEL  gmin;     /*  I   Minimum graylevel in the input image.         */
00234 extern double weight;   /*  I   Class count weight of the center pixel.       */
00235 extern double *count;   /*  O   Pointer to an array[nlev], containing         */
00236                         /*      the weighted number of pixels for each class. */
00237 PIXEL WTFMAJ (
00238 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00239 short jsize,    /*  I   Number of lines in the window.                        */
00240 short ksize     /*  I   Number of pixels/line in the window.                  */
00241 /*----------------------------------------------------------------------------*/
00242 ) { register short i, pclass, center=(short)(wnd[jsize/2*ksize+ksize/2]-gmin);
00243     double plurality;
00244 
00245     if (TIMES) TIMING(T_WTFMAJ);
00246 
00247     /* initialize counter array */
00248     for (i=0; i<nlev; count[i++]=0.0);
00249 
00250     /* count number of pixels in each class */
00251     for (i=0; i<jsize*ksize; i++) {
00252         count[(short)(wnd[i]-gmin)] += 1.0;
00253     }
00254 
00255     /* weight center pixel class */
00256     count[center] *= weight;
00257 
00258     /* find plurality class and make sure that plurality > simple majority and that plurality is unique */
00259     for (plurality=count[pclass=0],i=1; i<nlev; i++) {
00260         if (count[i] > plurality  &&  i != center) plurality = count[pclass=i];
00261     }
00262     if (plurality <= (double)(jsize*ksize)/2.0  ||  plurality <= count[center]) pclass = center;
00263 
00264     the_end:
00265     if (TIMES) TIMING(T_EXIT);
00266     return((PIXEL)pclass+gmin);
00267 }
00268 
00269 /*-Copyright Information------------------------------------------------------*/
00270 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00271 /*----------------------------------------------------------------------------*/
00272 /*-General Information--------------------------------------------------------*/
00273 /*                                                                            */
00274 /*   This library function computes a local maximum filter.                   */
00275 /*                                                                            */
00276 /*----------------------------------------------------------------------------*/
00277 /*-Interface Information------------------------------------------------------*/
00278 PIXEL WTFMAX (
00279 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00280 short jsize,    /*  I   Number of lines in the window.                        */
00281 short ksize     /*  I   Number of pixels/line in the window.                  */
00282 /*----------------------------------------------------------------------------*/
00283 ) { register short i;
00284     PIXEL maxval=wnd[0];
00285 
00286     if (TIMES) TIMING(T_WTFMAX);
00287 
00288     /* find local maximum */
00289     for (i=1; i<jsize*ksize; i++) {
00290         maxval = max(maxval,wnd[i]);
00291     }
00292 
00293     the_end:
00294     if (TIMES) TIMING(T_EXIT);
00295     return(maxval);
00296 }
00297 
00298 /*-Copyright Information------------------------------------------------------*/
00299 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00300 /*----------------------------------------------------------------------------*/
00301 /*-General Information--------------------------------------------------------*/
00302 /*                                                                            */
00303 /*   This library function computes a local minimum filter.                   */
00304 /*                                                                            */
00305 /*----------------------------------------------------------------------------*/
00306 /*-Interface Information------------------------------------------------------*/
00307 PIXEL WTFMIN (
00308 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00309 short jsize,    /*  I   Number of lines in the window.                        */
00310 short ksize     /*  I   Number of pixels/line in the window.                  */
00311 /*----------------------------------------------------------------------------*/
00312 ) { register short i;
00313     PIXEL minval=wnd[0];
00314 
00315     if (TIMES) TIMING(T_WTFMIN);
00316 
00317     /* find local minimum */
00318     for (i=1; i<jsize*ksize; i++) {
00319         minval = min(minval,wnd[i]);
00320     }
00321 
00322     the_end:
00323     if (TIMES) TIMING(T_EXIT);
00324     return(minval);
00325 }
00326 
00327 /*-Copyright Information------------------------------------------------------*/
00328 /* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
00329 /*----------------------------------------------------------------------------*/
00330 /*-General Information--------------------------------------------------------*/
00331 /*                                                                            */
00332 /*   This library function computes a local variance map.                     */
00333 /*                                                                            */
00334 /*----------------------------------------------------------------------------*/
00335 /*-Interface Information------------------------------------------------------*/
00336 PIXEL WTFVAR (
00337 PIXEL *wnd,     /*  I   Pointer to the window matrix[jsize][ksize].           */
00338 short jsize,    /*  I   Number of lines in the window.                        */
00339 short ksize     /*  I   Number of pixels/line in the window.                  */
00340 /*----------------------------------------------------------------------------*/
00341 ) { register short i;
00342     double area=jsize*ksize;
00343     PIXEL  sum=(PIXEL)0, sumsq=(PIXEL)0;
00344 
00345     if (TIMES) TIMING(T_WTFVAR);
00346 
00347     for (i=0; i<jsize*ksize; i++) {
00348         sum   += wnd[i];
00349         sumsq += wnd[i]*wnd[i];
00350     }
00351 
00352     the_end:
00353     if (TIMES) TIMING(T_EXIT);
00354     return((PIXEL)((double)sumsq/area-((double)sum/area)*((double)sum/area)));
00355 }