stats.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 computes the mean, standard deviation, minimum, maximum,   */
00009 /*   and range of graylevels in an image, as well as the number of pixels     */
00010 /*   with a zero graylevel, and outputs a histogram.                          */
00011 /*                                                                            */
00012 /*----------------------------------------------------------------------------*/
00013 /*-Interface Information------------------------------------------------------*/
00014 short STATS (
00015 IMAGE  *in,     /*  I   Pointer to the input image.                           */
00016 short  inc,     /*  I   Pixels/line and line increment.                       */
00017 short  nlev,    /*  I   > 0:   Number of bins in the histogram.               */
00018                 /*      = 0:   Histogram is suppressed.                       */
00019 short  option,  /*  I   Switch, indicating how to use input pixel range:      */
00020                 /*      INCLUDE - Use only this input range for computing     */
00021                 /*                statistics.                                 */
00022                 /*      EXCLUDE - Exclude this pixel range.                   */
00023                 /*      IGNORE  - Do not use the pixel range.                 */
00024 PIXEL  rangemin,/*  I   Minimum pixel value in range of interest.             */
00025 PIXEL  rangemax,/*  I   Maximum pixel value in range of interest.             */
00026 double *mean,   /*  O   Pointer to an array[nbnd], containing the             */
00027                 /*      mean graylevel for each band.                         */
00028 double *dev,    /*  O   Pointer to an array[nbnd], containing the             */
00029                 /*      graylevel standard deviation for each band.           */
00030 PIXEL  *gmin,   /*  O   Pointer to an array[nbnd], containing the             */
00031                 /*      minimum graylevel for each band.                      */
00032 PIXEL  *gmax,   /*  O   Pointer to an array[nbnd], containing the             */
00033                 /*      maximum graylevel for each band.                      */
00034 PIXEL  *grng,   /*  O   Pointer to an array[nbnd], containing the             */
00035                 /*      graylevel range for each band.                        */
00036 long   *zero    /*  O   Pointer to an array[nbnd], containing the             */
00037                 /*      number of pixels with a zero graylevel for each band. */
00038 /*----------------------------------------------------------------------------*/
00039 ) {
00040     register int i, j, k;
00041     char   msg[SLEN];
00042     long   *hist=0, maxval, npts;
00043     double sum, sumsq;
00044     PIXEL **arrays=NULL;
00045     char minlbl[20],maxlbl[20];
00046     int plotoption=NORM;
00047     long *include, *exclude;
00048 
00049 
00050     if (TIMES) TIMING(T_STATS);
00051     if (NAMES) {
00052         MESSAGE('I',"");
00053         MESSAGE('I',"STATS");
00054         MESSAGE('I',"");
00055         sprintf(msg," Input image:                         %s",in->text);
00056         MESSAGE('I',msg);
00057         sprintf(msg," Line and pixel subsample increment:  %d",inc);
00058         MESSAGE('I',msg);
00059         if (nlev > 0) {
00060             sprintf(msg," Histogram resolution:                %d",nlev);
00061             MESSAGE('I',msg);
00062         }
00063         if (option == INCLUDE) {
00064             sprintf(msg," Include only pixels in the graylevel range: %.2f to %.2f",rangemin,rangemax);
00065             MESSAGE('I',msg);
00066         } else if (option == EXCLUDE) {
00067             sprintf(msg," Exclude pixels in the graylevel range: %.2f to %.2f",rangemin,rangemax);
00068             MESSAGE('I',msg);
00069         }
00070         MESSAGE('I'," ...............");
00071     }
00072 
00073     /* check input */
00074     if (!CHECKIMG(in)) {
00075         MESSAGE('E'," Can't identify image.");
00076         goto the_end;
00077     } else if (inc <= 0) {
00078         MESSAGE('E'," Increment must be greater than zero.");
00079         goto the_end;
00080     } else if (nlev < 0) {
00081         MESSAGE('E'," Number of histogram bins must be greater than or equal to zero.");
00082         goto the_end;
00083     }
00084     
00085     /* compute and output statistics for each band */
00086     MESSAGE('I',"");
00087     if (option != IGNORE) {
00088         MESSAGE('I'," Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero   Included/Excluded");
00089     } else {
00090         MESSAGE('I'," Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero");
00091     }
00092     MESSAGE('I',"");
00093 
00094     RANGE(in);
00095     include = (long *)malloc((in->nbnd)*sizeof(long));
00096     exclude = (long *)malloc((in->nbnd)*sizeof(long));
00097     for (i=0; i<in->nbnd; i++) {
00098         gmin[i] = in->gmax;
00099         gmax[i] = in->gmin;
00100         zero[i] = 0L;
00101 
00102         if (option != IGNORE) {
00103             include[i] = exclude[i] = 0L;
00104         }
00105         
00106         for (sum=sumsq=0.0,j=0; j<in->nlin; j+=inc) {
00107             for (k=0; k<in->npix; k+=inc) {
00108                 switch (option) {
00109 
00110                 case INCLUDE:
00111                     if (in->data[i][j][k] >= rangemin && in->data[i][j][k] <= rangemax) {
00112                         sum    += (double)in->data[i][j][k];
00113                         sumsq  += (double)(in->data[i][j][k]*in->data[i][j][k]);
00114                         gmin[i] = min(gmin[i],in->data[i][j][k]);
00115                         gmax[i] = max(gmax[i],in->data[i][j][k]);
00116                         if (in->data[i][j][k] == (PIXEL)0) {
00117                             zero[i] += 1L;
00118                         }
00119                         include[i] += 1L;
00120                     } else {
00121                         exclude[i] += 1L;
00122                     }
00123                     break;
00124 
00125                 case EXCLUDE:
00126                     if (in->data[i][j][k] < rangemin || in->data[i][j][k] > rangemax) {
00127                         sum    += (double)in->data[i][j][k];
00128                         sumsq  += (double)(in->data[i][j][k]*in->data[i][j][k]);
00129                         gmin[i] = min(gmin[i],in->data[i][j][k]);
00130                         gmax[i] = max(gmax[i],in->data[i][j][k]);
00131                         if (in->data[i][j][k] == (PIXEL)0) {
00132                             zero[i] += 1L;
00133                         }
00134                         include[i] += 1L;
00135                     } else {
00136                         exclude[i] += 1L;
00137                     }
00138                     break;
00139                     
00140                 default:
00141                     sum    += (double)in->data[i][j][k];
00142                     sumsq  += (double)(in->data[i][j][k]*in->data[i][j][k]);
00143                     gmin[i] = min(gmin[i],in->data[i][j][k]);
00144                     gmax[i] = max(gmax[i],in->data[i][j][k]);
00145                     if (in->data[i][j][k] == (PIXEL)0) {
00146                         zero[i] += 1L;
00147                     }
00148                     break;
00149                         
00150                 }
00151 
00152             }
00153         }
00154         if (option != IGNORE) {
00155             if (include[i] > 0L) {
00156                 mean[i] = sum/((double)include[i]);
00157                 dev[i]  = sqrt(sumsq/((double)include[i])-mean[i]*mean[i]);
00158                 grng[i] = gmax[i]-gmin[i];
00159             } else {
00160                 mean[i] = 0.0;
00161                 dev[i]  = 0.0;
00162                 gmin[i] = gmax[i] = 0.0;
00163                 grng[i] = 0.0;
00164             }
00165             sprintf(msg,"%5d%12.4e%12.4e%12.4e%12.4e%10ld/%ld%15ld/%ld",i+1,(double)gmin[i],(double)gmax[i],mean[i],dev[i],zero[i],include[i]-zero[i],include[i],exclude[i]);
00166         } else {
00167             mean[i] = sum/((double)(in->nlin/inc+(in->nlin%inc ? 1:0))*(double)(in->npix/inc+(in->npix%inc ? 1:0)));
00168             dev[i]  = sqrt(sumsq/((double)(in->nlin/inc+(in->nlin%inc ? 1:0))*(double)(in->npix/inc+(in->npix%inc ? 1:0)))-mean[i]*mean[i]);
00169             grng[i] = gmax[i]-gmin[i];
00170             sprintf(msg,"%5d%12.4e%12.4e%12.4e%12.4e%10ld/%ld",i+1,(double)gmin[i],(double)gmax[i],mean[i],dev[i],zero[i],(long)in->nlin*(long)in->npix-zero[i]);
00171         }
00172         
00173         MESSAGE('I',msg);
00174     }
00175 
00176     /* allocate space and compute histograms */
00177     if (nlev > 0) {
00178 
00179         /* Allocate memory for plotting */
00180         arrays = (PIXEL**)malloc(in->nbnd * sizeof(PIXEL*));
00181         if (!arrays) {
00182             MESSAGE('E'," Error allocating memory for plotting!");
00183             goto the_end;
00184         }
00185         for (i=0; i<in->nbnd; i++) {
00186             arrays[i] = NULL;
00187             arrays[i] = (PIXEL*) malloc(nlev * sizeof(PIXEL));
00188             if (!arrays[i]) {
00189                 MESSAGE('E'," Error allocating memory for plotting!");
00190                 goto the_end;
00191             }
00192         }
00193 
00194         RANGE(in);
00195         if(!(hist=(long *)malloc(in->nbnd*nlev*sizeof(long)))) {
00196             MESSAGE('E'," Memory request failed.");
00197             goto the_end;
00198         }
00199         MESSAGE('I',"");
00200 
00201         if (option == INCLUDE) {
00202             sprintf(msg," Histogram has %d levels between %.4e and %.4e",nlev,(double)rangemin,(double)rangemax);
00203             MESSAGE('I',msg);
00204             STRICTHISTOGRM(in,inc,option,rangemin,rangemax,nlev,hist);
00205         } else if (option == EXCLUDE) {
00206             sprintf(msg," Histogram has %d levels between %.4e and %.4e",nlev,(double)in->gmin,(double)in->gmax);
00207             MESSAGE('I',msg);
00208             STRICTHISTOGRM(in,inc,option,rangemin,rangemax,nlev,hist);
00209         } else {
00210             sprintf(msg," Histogram has %d levels between %.4e and %.4e",nlev,(double)in->gmin,(double)in->gmax);
00211             MESSAGE('I',msg);
00212             HISTOGRM(in,inc,in->gmin,in->gmax,nlev,hist);
00213         }
00214         
00215         /* Create data arrays for plotting */
00216         for (i=0; i<in->nbnd; i++) {
00217             for (j=0; j < nlev; j++) {
00218                 arrays[i][j] = (PIXEL) hist[i*nlev+j];
00219             }
00220         }
00221 
00222 
00223         if (option == INCLUDE) {
00224             sprintf(minlbl,"%.2f",rangemin);
00225             sprintf(maxlbl,"%.2f",rangemax);
00226         } else {
00227             sprintf(minlbl,"%.2f",in->gmin);
00228             sprintf(maxlbl,"%.2f",in->gmax);
00229         }
00230         
00231         PLOT(arrays,in->nbnd,nlev,minlbl,maxlbl,plotoption);
00232     }
00233 
00234     
00235  the_end:
00236 
00237     /* Free Array */
00238     if (arrays) {
00239         for (i=0; i<in->nbnd; i++) {
00240             if (arrays[i]) free(arrays[i]);
00241             arrays[i] = NULL;
00242         }
00243         free(arrays);
00244         arrays=NULL;
00245     }
00246 
00247     if (hist)  free(hist);
00248     if (TIMES) TIMING(T_EXIT);
00249 }
00250