---

stats.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function computes the mean, standard deviation, minimum, maximum,   */
/*   and range of graylevels in an image, as well as the number of pixels     */
/*   with a zero graylevel, and outputs a histogram.                          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
short
STATS (IMAGE * in,              /*  I   Pointer to the input image.                           */
       short inc,               /*  I   Pixels/line and line increment.                       */
       short nlev,              /*  I   > 0:   Number of bins in the histogram.               */
       /*      = 0:   Histogram is suppressed.                       */
       short option,            /*  I   Switch, indicating how to use input pixel range:      */
       /*      INCLUDE - Use only this input range for computing     */
       /*                statistics.                                 */
       /*      EXCLUDE - Exclude this pixel range.                   */
       /*      IGNORE  - Do not use the pixel range.                 */
       PIXEL rangemin,          /*  I   Minimum pixel value in range of interest.             */
       PIXEL rangemax,          /*  I   Maximum pixel value in range of interest.             */
       double *mean,            /*  O   Pointer to an array[nbnd], containing the             */
       /*      mean graylevel for each band.                         */
       double *dev,             /*  O   Pointer to an array[nbnd], containing the             */
       /*      graylevel standard deviation for each band.           */
       PIXEL * gmin,            /*  O   Pointer to an array[nbnd], containing the             */
       /*      minimum graylevel for each band.                      */
       PIXEL * gmax,            /*  O   Pointer to an array[nbnd], containing the             */
       /*      maximum graylevel for each band.                      */
       PIXEL * grng,            /*  O   Pointer to an array[nbnd], containing the             */
       /*      graylevel range for each band.                        */
       long *zero               /*  O   Pointer to an array[nbnd], containing the             */
       /*      number of pixels with a zero graylevel for each band. */
/*----------------------------------------------------------------------------*/
  )
{
  register int i, j, k;
  char msg[SLEN];
  long *hist = 0, maxval, npts;
  double sum, sumsq;
  PIXEL **arrays = NULL;
  char minlbl[20], maxlbl[20];
  int plotoption = NORM;
  long *include, *exclude;

  if (TIMES)
    TIMING (T_STATS);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "STATS");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:                         %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Line and pixel subsample increment:  %d", inc);
      MESSAGE ('I', msg);
      if (nlev > 0)
        {
          sprintf (msg, " Histogram resolution:                %d", nlev);
          MESSAGE ('I', msg);
        }
      if (option == INCLUDE)
        {
          sprintf (msg,
                   " Include only pixels in the graylevel range: %.2f to %.2f",
                   rangemin, rangemax);
          MESSAGE ('I', msg);
        }
      else if (option == EXCLUDE)
        {
          sprintf (msg,
                   " Exclude pixels in the graylevel range: %.2f to %.2f",
                   rangemin, rangemax);
          MESSAGE ('I', msg);
        }
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (inc <= 0)
    {
      MESSAGE ('E', " Increment must be greater than zero.");
      goto the_end;
    }
  else if (nlev < 0)
    {
      MESSAGE ('E',
               " Number of histogram bins must be greater than or equal to zero.");
      goto the_end;
    }

  /* compute and output statistics for each band */
  MESSAGE ('I', "");
  if (option != IGNORE)
    {
      MESSAGE ('I',
               " Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero   Included/Excluded");
    }
  else
    {
      MESSAGE ('I',
               " Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero");
    }
  MESSAGE ('I', "");

  RANGE (in);
  include = (long *) malloc ((in->nbnd) * sizeof (long));
  exclude = (long *) malloc ((in->nbnd) * sizeof (long));
  for (i = 0; i < in->nbnd; i++)
    {
      gmin[i] = in->gmax;
      gmax[i] = in->gmin;
      zero[i] = 0L;

      if (option != IGNORE)
        {
          include[i] = exclude[i] = 0L;
        }

      for (sum = sumsq = 0.0, j = 0; j < in->nlin; j += inc)
        {
          for (k = 0; k < in->npix; k += inc)
            {
              switch (option)
                {

                case INCLUDE:
                  if (in->data[i][j][k] >= rangemin
                      && in->data[i][j][k] <= rangemax)
                    {
                      sum += (double) in->data[i][j][k];
                      sumsq +=
                        (double) (in->data[i][j][k] * in->data[i][j][k]);
                      gmin[i] = min (gmin[i], in->data[i][j][k]);
                      gmax[i] = max (gmax[i], in->data[i][j][k]);
                      if (in->data[i][j][k] == (PIXEL) 0)
                        {
                          zero[i] += 1L;
                        }
                      include[i] += 1L;
                    }
                  else
                    {
                      exclude[i] += 1L;
                    }
                  break;

                case EXCLUDE:
                  if (in->data[i][j][k] < rangemin
                      || in->data[i][j][k] > rangemax)
                    {
                      sum += (double) in->data[i][j][k];
                      sumsq +=
                        (double) (in->data[i][j][k] * in->data[i][j][k]);
                      gmin[i] = min (gmin[i], in->data[i][j][k]);
                      gmax[i] = max (gmax[i], in->data[i][j][k]);
                      if (in->data[i][j][k] == (PIXEL) 0)
                        {
                          zero[i] += 1L;
                        }
                      include[i] += 1L;
                    }
                  else
                    {
                      exclude[i] += 1L;
                    }
                  break;

                default:
                  sum += (double) in->data[i][j][k];
                  sumsq += (double) (in->data[i][j][k] * in->data[i][j][k]);
                  gmin[i] = min (gmin[i], in->data[i][j][k]);
                  gmax[i] = max (gmax[i], in->data[i][j][k]);
                  if (in->data[i][j][k] == (PIXEL) 0)
                    {
                      zero[i] += 1L;
                    }
                  break;

                }

            }
        }
      if (option != IGNORE)
        {
          if (include[i] > 0L)
            {
              mean[i] = sum / ((double) include[i]);
              dev[i] =
                sqrt (sumsq / ((double) include[i]) - mean[i] * mean[i]);
              grng[i] = gmax[i] - gmin[i];
            }
          else
            {
              mean[i] = 0.0;
              dev[i] = 0.0;
              gmin[i] = gmax[i] = 0.0;
              grng[i] = 0.0;
            }
          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]);
        }
      else
        {
          mean[i] =
            sum / ((double) (in->nlin / inc + (in->nlin % inc ? 1 : 0)) *
                   (double) (in->npix / inc + (in->npix % inc ? 1 : 0)));
          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]);
          grng[i] = gmax[i] - gmin[i];
          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]);
        }

      MESSAGE ('I', msg);
    }

  /* allocate space and compute histograms */
  if (nlev > 0)
    {

      /* Allocate memory for plotting */
      arrays = (PIXEL **) malloc (in->nbnd * sizeof (PIXEL *));
      if (!arrays)
        {
          MESSAGE ('E', " Error allocating memory for plotting!");
          goto the_end;
        }
      for (i = 0; i < in->nbnd; i++)
        {
          arrays[i] = NULL;
          arrays[i] = (PIXEL *) malloc (nlev * sizeof (PIXEL));
          if (!arrays[i])
            {
              MESSAGE ('E', " Error allocating memory for plotting!");
              goto the_end;
            }
        }

      RANGE (in);
      if (!(hist = (long *) malloc (in->nbnd * nlev * sizeof (long))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }
      MESSAGE ('I', "");

      if (option == INCLUDE)
        {
          sprintf (msg, " Histogram has %d levels between %.4e and %.4e",
                   nlev, (double) rangemin, (double) rangemax);
          MESSAGE ('I', msg);
          STRICTHISTOGRM (in, inc, option, rangemin, rangemax, nlev, hist);
        }
      else if (option == EXCLUDE)
        {
          sprintf (msg, " Histogram has %d levels between %.4e and %.4e",
                   nlev, (double) in->gmin, (double) in->gmax);
          MESSAGE ('I', msg);
          STRICTHISTOGRM (in, inc, option, rangemin, rangemax, nlev, hist);
        }
      else
        {
          sprintf (msg, " Histogram has %d levels between %.4e and %.4e",
                   nlev, (double) in->gmin, (double) in->gmax);
          MESSAGE ('I', msg);
          HISTOGRM (in, inc, in->gmin, in->gmax, nlev, hist);
        }

      /* Create data arrays for plotting */
      for (i = 0; i < in->nbnd; i++)
        {
          for (j = 0; j < nlev; j++)
            {
              arrays[i][j] = (PIXEL) hist[i * nlev + j];
            }
        }

      if (option == INCLUDE)
        {
          sprintf (minlbl, "%.2f", rangemin);
          sprintf (maxlbl, "%.2f", rangemax);
        }
      else
        {
          sprintf (minlbl, "%.2f", in->gmin);
          sprintf (maxlbl, "%.2f", in->gmax);
        }

      PLOT (arrays, in->nbnd, nlev, minlbl, maxlbl, plotoption);
    }

the_end:

  /* Free Array */
  if (arrays)
    {
      for (i = 0; i < in->nbnd; i++)
        {
          if (arrays[i])
            free (arrays[i]);
          arrays[i] = NULL;
        }
      free (arrays);
      arrays = NULL;
    }

  if (hist)
    free (hist);
  if (TIMES)
    TIMING (T_EXIT);
}

---