---

covar.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 band-to-band covariance and correlation       */
/*   matrices of an image.                                                    */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
COVAR (IMAGE * in,              /*  I   Pointer to the image.                                 */
       short inc,               /*  I   Pixels/line and line increment.                       */
       double *cov,             /*  O   Pointer to the covariance matrix.                     */
       double *cor              /*  O   Pointer to the correlation matrix.                    */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, l;
  char msg[SLEN];
  double *sum = 0, *sumsq = 0, count, pinc, psum;

  if (TIMES)
    TIMING (T_COVAR);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "COVAR");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:                          %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Line and pixel subsample increment:   %d", inc);
      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', " Subsample increment must be greater than zero.");
      goto the_end;
    }

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

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc =
    2.0 / (double) (in->nbnd +
                    1) / (double) in->nbnd / (double) (in->nlin / inc +
                                                       (in->nlin %
                                                        inc ? 1 : 0));

  /* initialize */
  for (i = 0; i < in->nbnd; i++)
    {
      sum[i] = 0.0;
      for (j = 0; j < in->nbnd; j++)
        {
          sumsq[i * in->nbnd + j] = 0.0;
        }
    }

  /* compute the sums */
  for (i = 0; i < in->nbnd; i++)
    {
      for (j = 0; j < in->nlin; j += inc)
        {
          for (k = 0; k < in->npix; k += inc)
            {
              sum[i] += (double) in->data[i][j][k];
              sumsq[i * in->nbnd + i] +=
                (double) in->data[i][j][k] * in->data[i][j][k];
            }
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }
    }
  for (i = 0; i < in->nbnd - 1; i++)
    {
      for (l = i + 1; l < in->nbnd; l++)
        {
          for (j = 0; j < in->nlin; j += inc)
            {
              for (k = 0; k < in->npix; k += inc)
                {
                  sumsq[i * in->nbnd + l] = sumsq[l * in->nbnd + i] +=
                    (double) (in->data[i][j][k] * in->data[l][j][k]);
                }
              if (LINES && !PROGRESS (psum += pinc))
                goto the_end;
            }
        }
    }

  /* calculate the covariance matrix */
  count =
    (double) (in->nlin / inc +
              (in->nlin % inc ? 1 : 0)) * (double) (in->npix / inc +
                                                    (in->npix % inc ? 1 : 0));
  for (i = 0; i < in->nbnd; i++)
    {
      for (j = 0; j < in->nbnd; j++)
        {
          cov[i * in->nbnd + j] =
            (sumsq[i * in->nbnd + j] - sum[i] * sum[j] / count) / (count -
                                                                   1.0);
        }
    }

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

  /* output both matrices */
  MESSAGE ('I', "");
  MESSAGE ('I', " Covariance Matrix");
  for (i = 0; i < in->nbnd; i++)
    {
      for (j = k = 0, msg[k++] = ' '; j < in->nbnd; j += 1, k = strlen (msg))
        {
          sprintf (msg + k, "%12.4e", cov[i * in->nbnd + j]);
        }
      MESSAGE ('I', msg);
    }
  MESSAGE ('I', "");
  MESSAGE ('I', " Correlation Matrix");
  for (i = 0; i < in->nbnd; i++)
    {
      for (j = k = 0, msg[k++] = ' '; j < in->nbnd; j += 1, k = strlen (msg))
        {
          sprintf (msg + k, "%12.4e", cor[i * in->nbnd + j]);
        }
      MESSAGE ('I', msg);
    }

the_end:
  if (sum)
    free (sum);
  if (sumsq)
    free (sumsq);
  if (TIMES)
    TIMING (T_EXIT);
}

---