---

roistats.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1993 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function computes the minimum, maximum, mean, standard deviation,   */
/*   and covariance matrix from a set of graylevel values.                    */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
short
ROISTATS (char *name,           /*  I                                                         */
          short nbnd,           /*  I                                                         */
          long nval,            /*  I                                                         */
          PIXEL * gval,         /*  I                                                         */
          short inc,            /*  I   Pixel subsample increment.                            */
          short nlev,           /*  I   > 0:   Number of bins in the histogram.               */
          /*      = 0:   Histogram is suppressed.                       */
          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.                      */
          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.           */
          long *zero,           /*  O   Pointer to an array[nbnd], containing the             */
          /*      number of pixels with a zero graylevel for each band. */
          double *cov,          /*  O   Pointer to an array[nbnd*nbnd], containing the        */
          /*      graylevel covariance matrix.                          */
          double *cor           /*  O   Pointer to an array[nbnd*nbnd], containing the        */
          /*      graylevel correlation matrix.                         */
/*----------------------------------------------------------------------------*/
  )
{
  register long i, j, k;
  register PIXEL p, q;
  char msg[SLEN];
  long maxval, npts, *hist = 0;
  double factor, *sum = 0, *sumsq = 0;
  PIXEL gmin_all, gmax_all;
  PIXEL **arrays = NULL;
  char minlbl[20], maxlbl[20];
  int option = NORM;

  if (TIMES)
    TIMING (T_ROISTATS);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "ROISTATS");
      MESSAGE ('I', "");
      sprintf (msg, " Input region:               %s", name);
      MESSAGE ('I', msg);
      sprintf (msg, " Pixel subsample increment:  %d", inc);
      MESSAGE ('I', msg);
      if (nlev > 0)
        {
          sprintf (msg, " Histogram resolution:       %d", nlev);
          MESSAGE ('I', msg);
        }
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (nbnd <= 0)
    {
      MESSAGE ('E', " Number of bands must be greater than zero.");
      goto the_end;
    }
  else if (nval <= 0)
    {
      MESSAGE ('E',
               " Number of values must be greater than or equal to zero.");
      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;
    }

  /* allocate space for summation arrays */
  if (!(sum = (double *) malloc (nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(sumsq = (double *) malloc (nbnd * nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }

  /* initialize summation arrays */
  for (i = 0; i < nbnd; i++)
    {
      gmin[i] = gval[i * nval];
      gmax[i] = gval[i * nval];
      zero[i] = 0L;
      sum[i] = 0.0;
      for (j = 0; j < nbnd; j++)
        {
          sumsq[i * nbnd + j] = 0.0;
        }
    }

  /* compute the sums */
  for (i = 0; i < nbnd; i += 1)
    {
      for (j = 0; j < nval; j += inc)
        {
          p = gval[i * nval + j];
          gmin[i] = min (gmin[i], p);
          gmax[i] = max (gmax[i], p);
          zero[i] += (long) (p == (PIXEL) 0);
          sum[i] += (double) p;
          sumsq[i * nbnd + i] += (double) (p * p);
        }
    }
  for (i = 0; i < nbnd - 1; i += 1)
    {
      for (j = i + 1; j < nbnd; j += 1)
        {
          for (k = 0; k < nval; k += inc)
            {
              p = gval[i * nval + k];
              q = gval[j * nval + k];
              sumsq[i * nbnd + j] = sumsq[j * nbnd + i] += (double) (p * q);
            }
        }
    }

  /* calculate the statistics */
  for (i = 0; i < nbnd; i++)
    {
      mean[i] = sum[i] / (double) (nval / inc + (nval % inc ? 1 : 0));
      dev[i] =
        sqrt (sumsq[i * nbnd + i] /
              (double) (nval / inc + (nval % inc ? 1 : 0)) -
              mean[i] * mean[i]);
    }

  /* output the statistics */
  MESSAGE ('I', "");
  MESSAGE ('I',
           " Band   Minimum     Maximum       Mean     Deviation       Zero/Non-Zero");
  for (i = 0; i < nbnd; i++)
    {
      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]);
      MESSAGE ('I', msg);
    }

  if (nlev > 0)
    {
      /* Allocate memory for plotting */
      arrays = (PIXEL **) malloc (nbnd * sizeof (PIXEL *));
      if (!arrays)
        {
          MESSAGE ('E', " Error allocating memory for plotting!");
          goto the_end;
        }
      for (i = 0; i < 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;
            }
        }

      /* allocate space for histogram array */
      if (!(hist = (long *) malloc (nbnd * nlev * sizeof (long))))
        {
          MESSAGE ('E', " Memory request failed.");
          goto the_end;
        }

      /* initialize */
      for (i = 0; i < nbnd * nlev; i++)
        {
          hist[i] = 0L;
        }
      p = gmin[0];
      q = gmax[0];
      for (i = 1; i < nbnd; i++)
        {
          p = min (p, gmin[i]);
          q = max (q, gmax[i]);
        }
      factor = (double) nlev / (double) (q - p + 1);

      /* compute histogram(s) */
      for (i = 0; i < nbnd; i++)
        {
          for (j = 0; j < nval; j += inc)
            {
              if (gval[i * nval + j] <= p)
                {
                  hist[i * nlev + 0] += 1L;
                }
              else if (p < gval[i * nval + j] && gval[i * nval + j] < q)
                {
                  hist[i * nlev +
                       (short) ((double) (gval[i * nval + j] - p) *
                                factor)] += 1L;
                }
              else
                {
                  hist[i * nlev + nlev - 1] += 1L;
                }
            }
        }

      /* find largest value */
      for (maxval = hist[0], j = 0; j < nbnd * nlev; j++)
        {
          maxval = max (maxval, hist[j]);
        }

      gmin_all = gmin[0];
      gmax_all = gmax[0];
      for (i = 0; i < nbnd; i++)
        {
          gmin_all = min (gmin_all, gmin[i]);
          gmax_all = max (gmax_all, gmax[i]);
        }

      /* output histogram(s) */
      MESSAGE ('I', "");
      sprintf (msg, " Histogram has %d levels between %.4e and %.4e", nlev, p,
               q);
      MESSAGE ('I', msg);

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

      sprintf (minlbl, "%.2f", gmin_all);
      sprintf (maxlbl, "%.2f", gmax_all);

      PLOT (arrays, nbnd, nlev, minlbl, maxlbl, option);

    }

  if (cov != NULL)
    {

      /* calculate covariance matrix */
      for (i = 0; i < nbnd; i++)
        {
          for (j = 0; j < nbnd; j++)
            {
              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));
            }
        }

      /* output covariance matrix */
      MESSAGE ('I', "");
      MESSAGE ('I', " Covariance Matrix");
      for (i = 0; i < nbnd; i++)
        {
          for (j = k = 0, msg[k++] = ' '; j < nbnd; j += 1, k = strlen (msg))
            {
              sprintf (msg + k, "%12.4e", cov[i * nbnd + j]);
            }
          MESSAGE ('I', msg);
        }
    }

  if (cov != NULL && cor != NULL)
    {

      /* calculate correlation matrix */
      for (i = 0; i < nbnd; i++)
        {
          for (j = 0; j < nbnd; j++)
            {
              cor[i * nbnd + j] =
                cov[i * nbnd +
                    j] / sqrt (cov[i * nbnd + i] * cov[j * nbnd + j]);
            }
        }

      /* output correlation matrix */
      MESSAGE ('I', "");
      MESSAGE ('I', " Correlation Matrix");
      for (i = 0; i < nbnd; i++)
        {
          for (j = k = 0, msg[k++] = ' '; j < nbnd; j += 1, k = strlen (msg))
            {
              sprintf (msg + k, "%12.4e", cor[i * nbnd + j]);
            }
          MESSAGE ('I', msg);
        }
    }

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

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

---