---

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

---