roistats.c Source File

00001 #include        "sadie.h"
00002 
00003 /*-Copyright Information------------------------------------------------------*/
00004 /* Copyright (c) 1993 by the University of Arizona Digital Image Analysis Lab */
00005 /*----------------------------------------------------------------------------*/
00006 /*-General Information--------------------------------------------------------*/
00007 /*                                                                            */
00008 /*   This function computes the minimum, maximum, mean, standard deviation,   */
00009 /*   and covariance matrix from a set of graylevel values.                    */
00010 /*                                                                            */
00011 /*----------------------------------------------------------------------------*/
00012 /*-Interface Information------------------------------------------------------*/
00013 short ROISTATS (
00014 char   *name,   /*  I                                                         */
00015 short  nbnd,    /*  I                                                         */
00016 long   nval,    /*  I                                                         */
00017 PIXEL  *gval,   /*  I                                                         */
00018 short  inc,     /*  I   Pixel subsample increment.                            */
00019 short  nlev,    /*  I   > 0:   Number of bins in the histogram.               */
00020                 /*      = 0:   Histogram is suppressed.                       */
00021 PIXEL  *gmin,   /*  O   Pointer to an array[nbnd], containing the             */
00022                 /*      minimum graylevel for each band.                      */
00023 PIXEL  *gmax,   /*  O   Pointer to an array[nbnd], containing the             */
00024                 /*      maximum graylevel for each band.                      */
00025 double *mean,   /*  O   Pointer to an array[nbnd], containing the             */
00026                 /*      mean graylevel for each band.                         */
00027 double *dev,    /*  O   Pointer to an array[nbnd], containing the             */
00028                 /*      graylevel standard deviation for each band.           */
00029 long   *zero,   /*  O   Pointer to an array[nbnd], containing the             */
00030                 /*      number of pixels with a zero graylevel for each band. */
00031 double *cov,    /*  O   Pointer to an array[nbnd*nbnd], containing the        */
00032                 /*      graylevel covariance matrix.                          */
00033 double *cor     /*  O   Pointer to an array[nbnd*nbnd], containing the        */
00034                 /*      graylevel correlation matrix.                         */
00035 /*----------------------------------------------------------------------------*/
00036 ) { register long i, j, k;
00037         register PIXEL p, q;
00038     char   msg[SLEN];
00039     long   maxval, npts, *hist=0;
00040     double factor, *sum=0, *sumsq=0;
00041     PIXEL gmin_all, gmax_all;
00042     PIXEL **arrays=NULL;
00043     char minlbl[20],maxlbl[20];
00044     int option=NORM;
00045 
00046 
00047     if (TIMES) TIMING(T_ROISTATS);
00048     if (NAMES) {
00049         MESSAGE('I',"");
00050         MESSAGE('I',"ROISTATS");
00051         MESSAGE('I',"");
00052         sprintf(msg," Input region:               %s",name);
00053         MESSAGE('I',msg);
00054         sprintf(msg," Pixel subsample increment:  %d",inc);
00055         MESSAGE('I',msg);
00056         if (nlev > 0) {
00057                 sprintf(msg," Histogram resolution:       %d",nlev);
00058                 MESSAGE('I',msg);
00059             }
00060         MESSAGE('I'," ...............");
00061     }
00062 
00063     /* check input */
00064     if (nbnd <= 0) {
00065         MESSAGE('E'," Number of bands must be greater than zero.");
00066         goto the_end;
00067     } else if (nval <= 0) {
00068         MESSAGE('E'," Number of values must be greater than or equal to zero.");
00069         goto the_end;
00070     } else if (inc <= 0) {
00071         MESSAGE('E'," Increment must be greater than zero.");
00072         goto the_end;
00073     } else if (nlev < 0) {
00074         MESSAGE('E'," Number of histogram bins must be greater than or equal to zero.");
00075         goto the_end;
00076     }
00077 
00078     /* allocate space for summation arrays */
00079     if (!(sum=(double *)malloc(nbnd*sizeof(double)))) {
00080         MESSAGE('E'," Memory request failed.");
00081         goto the_end;
00082     }
00083     if (!(sumsq=(double *)malloc(nbnd*nbnd*sizeof(double)))) {
00084         MESSAGE('E'," Memory request failed.");
00085         goto the_end;
00086     }
00087 
00088     /* initialize summation arrays */
00089     for (i=0; i<nbnd; i++) {
00090         gmin[i] = gval[i*nval];
00091         gmax[i] = gval[i*nval];
00092         zero[i] = 0L;
00093         sum[i]  = 0.0;
00094         for (j=0; j<nbnd; j++) {
00095             sumsq[i*nbnd+j] = 0.0;
00096         }
00097     }
00098     
00099     /* compute the sums */
00100     for (i=0; i<nbnd; i+=1) {
00101         for (j=0; j<nval; j+=inc) {
00102             p = gval[i*nval+j];
00103             gmin[i]          = min(gmin[i],p);
00104             gmax[i]          = max(gmax[i],p);
00105             zero[i]         += (long)(p == (PIXEL)0);
00106             sum[i]          += (double)p;
00107             sumsq[i*nbnd+i] += (double)(p*p);
00108         }
00109     }
00110     for (i=0; i<nbnd-1; i+=1) {
00111         for (j=i+1; j<nbnd; j+=1) {
00112             for (k=0; k<nval; k+=inc) {
00113                 p = gval[i*nval+k];
00114                 q = gval[j*nval+k];
00115                 sumsq[i*nbnd+j] = sumsq[j*nbnd+i] += (double)(p*q);
00116             }
00117         }
00118     }
00119 
00120     /* calculate the statistics */
00121     for (i=0; i<nbnd; i++) {
00122         mean[i] = sum[i]/(double)(nval/inc+(nval%inc ? 1:0));
00123         dev[i]  = sqrt(sumsq[i*nbnd+i]/(double)(nval/inc+(nval%inc ? 1:0)) - mean[i]*mean[i]);
00124     }
00125 
00126     /* output the statistics */
00127     MESSAGE('I',"");
00128     MESSAGE('I'," Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero");
00129     for (i=0; i<nbnd; i++) {
00130         sprintf(msg,"%5ld%12.4e%12.4e%12.4e%12.4e%10ld/%ld",i+1,(double)gmin[i],(double)gmax[i],mean[i],dev[i],zero[i],nval-zero[i]);
00131         MESSAGE('I',msg);
00132     }
00133 
00134     if (nlev > 0) {
00135         /* Allocate memory for plotting */
00136         arrays = (PIXEL**)malloc(nbnd * sizeof(PIXEL*));
00137         if (!arrays) {
00138             MESSAGE('E'," Error allocating memory for plotting!");
00139             goto the_end;
00140         }
00141         for (i=0; i<nbnd; i++) {
00142             arrays[i] = NULL;
00143             arrays[i] = (PIXEL*) malloc(nlev * sizeof(PIXEL));
00144             if (!arrays[i]) {
00145                 MESSAGE('E'," Error allocating memory for plotting!");
00146                 goto the_end;
00147             }
00148         }
00149 
00150         /* allocate space for histogram array */
00151         if(!(hist=(long *)malloc(nbnd*nlev*sizeof(long)))) {
00152             MESSAGE('E'," Memory request failed.");
00153             goto the_end;
00154         }
00155 
00156         /* initialize */
00157         for (i=0; i<nbnd*nlev; i++) {
00158             hist[i] = 0L;
00159         }
00160         p = gmin[0];
00161         q = gmax[0];
00162         for (i=1; i<nbnd; i++) {
00163                 p = min(p,gmin[i]);
00164                 q = max(q,gmax[i]);
00165         }
00166         factor = (double)nlev/(double)(q-p+1);
00167 
00168         /* compute histogram(s) */
00169         for (i=0; i<nbnd; i++) {
00170             for (j=0; j<nval; j+=inc) {
00171                 if (gval[i*nval+j] <= p) {
00172                     hist[i*nlev+0] += 1L;
00173                 } else if (p < gval[i*nval+j]  &&  gval[i*nval+j] < q) {
00174                     hist[i*nlev+(short)((double)(gval[i*nval+j]-p)*factor)] += 1L;
00175                 } else {
00176                     hist[i*nlev+nlev-1] += 1L;
00177                 }
00178             }
00179         }
00180 
00181         /* find largest value */
00182         for (maxval=hist[0],j=0; j<nbnd*nlev; j++) {
00183             maxval = max(maxval,hist[j]);
00184         }
00185         
00186         gmin_all = gmin[0];
00187         gmax_all = gmax[0];
00188         for (i=0; i<nbnd; i++) {
00189           gmin_all = min(gmin_all, gmin[i]);
00190           gmax_all = max(gmax_all, gmax[i]);
00191         }
00192 
00193         /* output histogram(s) */
00194         MESSAGE('I',"");
00195         sprintf(msg," Histogram has %d levels between %.4e and %.4e",nlev,p,q);
00196         MESSAGE('I',msg);
00197 
00198 
00199         /* Create data arrays for plotting */
00200         for (i=0; i<nbnd; i++) {
00201             for (j=0; j < nlev; j++) {
00202                 arrays[i][j] = (PIXEL) hist[i*nlev+j];
00203             }
00204         }
00205 
00206         sprintf(minlbl,"%.2f",gmin_all);
00207         sprintf(maxlbl,"%.2f",gmax_all);
00208 
00209         PLOT(arrays,nbnd,nlev,minlbl,maxlbl,option);
00210 
00211     }
00212         
00213     if (cov != NULL) {
00214 
00215         /* calculate covariance matrix */
00216         for (i=0; i<nbnd; i++) {
00217             for (j=0; j<nbnd; j++) {
00218                 cov[i*nbnd+j] = (sumsq[i*nbnd+j]-sum[i]*sum[j]/(double)(nval/inc+(nval%inc ? 1:0)))/((double)(nval/inc+(nval%inc ? 1:0)-1));
00219             }
00220         }
00221 
00222         /* output covariance matrix */
00223         MESSAGE('I',"");
00224         MESSAGE('I'," Covariance Matrix");
00225         for (i=0; i<nbnd; i++) {
00226             for (j=k=0,msg[k++]=' '; j<nbnd; j+=1,k=strlen(msg)) {
00227                 sprintf(msg+k,"%12.4e",cov[i*nbnd+j]);
00228             }
00229             MESSAGE('I',msg);
00230         }
00231     }
00232 
00233     if (cov != NULL  &&  cor != NULL) {
00234 
00235         /* calculate correlation matrix */
00236         for (i=0; i<nbnd; i++) {
00237             for (j=0; j<nbnd; j++) {
00238                 cor[i*nbnd+j] = cov[i*nbnd+j]/sqrt(cov[i*nbnd+i]*cov[j*nbnd+j]);
00239             }
00240         }
00241 
00242         /* output correlation matrix */
00243         MESSAGE('I',"");
00244         MESSAGE('I'," Correlation Matrix");
00245         for (i=0; i<nbnd; i++) {
00246             for (j=k=0,msg[k++]=' '; j<nbnd; j+=1,k=strlen(msg)) {
00247                 sprintf(msg+k,"%12.4e",cor[i*nbnd+j]);
00248             }
00249             MESSAGE('I',msg);
00250         }
00251     }
00252 
00253  the_end:
00254     /* Free Array */
00255     if (arrays) {
00256         for (i=0; i<nbnd; i++) {
00257             if (arrays[i]) free(arrays[i]);
00258             arrays[i] = NULL;
00259         }
00260         free(arrays);
00261         arrays=NULL;
00262     }
00263 
00264     if (hist)  free(hist);
00265     if (sum)   free(sum);
00266     if (sumsq) free(sumsq);
00267     if (TIMES) TIMING(T_EXIT);
00268 }