---

decorstr.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1992 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function performs a decorrelation stretch on a multi-band           */
/*   image. The transformed image is uncorrelated, and each band has a        */
/*   variance of the average of the original band variances. The mean         */
/*   of each band remains unchanged.                                          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
DECORSTR (IMAGE * in,           /*  I   Pointer to the input image                            */
          short inc,            /*  I   Pixels/line and line increment to be used             */
          /*      for calculating the covariance matrix.                */
          IMAGE ** out          /*  O   Address of a pointer to the output image.             */
          /*      The input image may be overwritten (out = &in).       */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k, l;
  double *mean = 0, *dev = 0, *cov = 0, *cor = 0, *eig = 0, *pct =
    0, avgdev, sum, pinc, psum;
  PIXEL *gmin = 0, *gmax = 0, *grng = 0;
  long *zero = 0;

  if (TIMES)
    TIMING (T_DECORSTR);

  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "DECORSTR");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (in->nbnd == 1)
    {
      MESSAGE ('E', " Image must be multi-band.");
      goto the_end;
    }
  else if (inc <= 0)
    {
      MESSAGE ('E', " Subsample increment must be greater than zero.");
      goto the_end;
    }

  /* create image of appropriate size */
  if (!CHECKIMG (*out))
    GETMEM (in->nbnd, in->nlin, in->npix, out);
  if (!*out)
    goto the_end;

  /* allocate space */
  if (!(mean = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(dev = (double *) malloc (in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(gmin = (PIXEL *) malloc (in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(gmax = (PIXEL *) malloc (in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(grng = (PIXEL *) malloc (in->nbnd * sizeof (PIXEL))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(zero = (long *) malloc (in->nbnd * sizeof (long))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(cov = (double *) malloc (in->nbnd * in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(cor = (double *) malloc (in->nbnd * in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(eig = (double *) malloc (in->nbnd * in->nbnd * sizeof (double))))
    {
      MESSAGE ('E', " Memory request failed.");
      goto the_end;
    }
  if (!(pct = (double *) malloc (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 = 1.0 / (double) in->nlin;

  /* calculate statistics, covariance matrix and eigenvectors of the input image */
  STATS (in, inc, (short) 0, IGNORE, 0.0, 0.0, mean, dev, gmin, gmax, grng,
         zero);
  if (ABORT)
    goto the_end;
  avgdev = 0.0;
  for (i = 0; i < in->nbnd; i++)
    {
      avgdev += dev[i];
    }
  avgdev /= in->nbnd;
  COVAR (in, inc, cov, cor);
  if (ABORT)
    goto the_end;
  EIGEN (cov, in->nbnd, eig);
  if (ABORT)
    goto the_end;

  for (j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < in->npix; k++)
        {

          /* perform principal component transform on zero mean image */
          for (i = 0; i < in->nbnd; i++)
            {
              for (pct[i] = 0.0, l = 0; l < in->nbnd; l++)
                {
                  pct[i] +=
                    (double) ((in->data[l][j][k] -
                               mean[l]) * eig[l * in->nbnd + i]);
                }
            }

          /* equalize variance of separate principal component scores and do inverse PCT.  Add mean back in */
          for (i = 0; i < in->nbnd; i++)
            {
              for (sum = 0.0, l = 0; l < in->nbnd; l++)
                {
                  sum +=
                    eig[i * in->nbnd +
                        l] * ((pct[l] * avgdev) /
                              sqrt (cov[l * in->nbnd + l]));
                }
              (*out)->data[i][j][k] = (PIXEL) (sum + mean[i]);
            }

        }
      if (LINES && !PROGRESS (psum += pinc))
        goto the_end;
    }

the_end:
  if (cov)
    free (cov);
  if (cor)
    free (cor);
  if (eig)
    free (eig);
  if (pct)
    free (pct);
  if (mean)
    free (mean);
  if (dev)
    free (dev);
  if (gmin)
    free (gmin);
  if (gmax)
    free (gmax);
  if (grng)
    free (grng);
  if (zero)
    free (zero);
  if (TIMES)
    TIMING (T_EXIT);
}

---