---

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);
}

---