---

periodogram.c

Go to the documentation of this file.

#include        "sadie.h"
#include        "proto.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function computes the periodogram of an image, using the            */
/*   Fast Fourier transform.                                                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void PERIODOGRAM (
IMAGE *in,      /*  I   Pointer to the input image.                           */
IMAGE **out     /*  O   Address of a pointer to the output image.             */
                /*      The input image may be overwritten (out = &in).       */
/*----------------------------------------------------------------------------*/
) { register short i, j, k;
    IMAGE *window=NULL;
    IMAGE *periodogram=NULL;
    IMAGE *section=NULL;
    IMAGE *scratch=NULL;
    int pheight, pwidth;
    PIXEL hwin, vwin;
    int vcount, hcount, voffset, hoffset;
    char msg[SLEN];

    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"PERIODOGRAM");
        MESSAGE('I',"");
        sprintf(msg," Input image:   %s",in->text);
        MESSAGE('I',msg);
        MESSAGE('I'," ...............");
    }


NAMES = 0;


    /* check input */
    if (!CHECKIMG(in)) {
        MESSAGE('E'," Can't identify image.");
        goto the_end;
    } else if (in->nlin < 4) {
        MESSAGE('E'," Number of lines must be at least four.");
        goto the_end;
    } else if (in->npix < 4) {
        MESSAGE('E'," Number of pixels/line must be at least four.");
        goto the_end;
    }
    
    pwidth = (int) pow(2.0,floor(log10((double)in->npix)/log10(2.0)) - 1.0);
    pheight = (int) pow(2.0,floor(log10((double)in->nlin)/log10(2.0)) - 1.0);
    
    if (in->nlin > (pheight*2)) {
        sprintf(msg, "Only the first %d lines will be used to compute the periodogram.", pheight*2);
        MESSAGE('I', msg);
    } else if (in->npix > (pwidth*2)) {
        sprintf(msg, "Only the first %d pixels in each line will be used to compute the periodogram.", pwidth*2);
        MESSAGE('I', msg);
    }
    


    /* Compute a 2-D window function assuming that lines and columns are independent */
    /* The 1-D equivalent of the window function is:

                      /             \  2
                     |       L-1     | 
                     |   j - ---     |
                     |        2      |
          W(j) = 1 - |  ----------   |      , j = 0,1,...,L-1
                     |     L+1       |
                     |     ---       |
                     |      2        |
                      \             /
     */
     
    /* create window image of appropriate size */
    if (!CHECKIMG(window)) GETMEM(1,pheight,pwidth,&window);
    if (!window) goto the_end;

    for (j=0; j<pheight; j++) {
        vwin = (PIXEL) (1.0 - ((((PIXEL)j - ((PIXEL)pheight-1.0)/2.0)/(((PIXEL)pheight+1.0)/2.0)) * (((PIXEL)j - ((PIXEL)pheight-1.0)/2.0)/(((PIXEL)pheight+1.0)/2.0))));
        for (k=0; k<pwidth; k++) {
            hwin = (PIXEL) (1.0 - ((((PIXEL)k - ((PIXEL)pwidth-1.0)/2.0)/(((PIXEL)pwidth+1.0)/2.0)) * (((PIXEL)k - ((PIXEL)pwidth-1.0)/2.0)/(((PIXEL)pwidth+1.0)/2.0))));
            window->data[0][j][k] = hwin * vwin;
        }
    }     
     
     
     
     
    /* create image of appropriate size */
    if (!CHECKIMG(*out)) GETMEM(in->nbnd,pheight,pwidth,out);
    if (!*out) goto the_end;

    /* initialize the output */
    for (i=0; i<(*out)->nbnd; i++) {
        for (j=0; j<(*out)->nlin; j++) {
            for (k=0; k<(*out)->npix; k++) {
                (*out)->data[i][j][k] = (PIXEL)0.0;
            }
        }
    }




    /* create temporary image of appropriate size to hold windowed subsections of input */
    if (!CHECKIMG(section)) GETMEM(in->nbnd,pheight,pwidth,&section);
    if (!section) goto the_end;
    
    /* create temporary image of appropriate size to hold individual periodograms */
    if (!CHECKIMG(scratch)) GETMEM(in->nbnd,pheight,pwidth,&scratch);
    if (!scratch) goto the_end;
    
    
    voffset=(int)rnd((double)pheight/2.0);
    hoffset=(int)rnd((double)pwidth/2.0);
    for (vcount=0; vcount<3; vcount++) {
        for (hcount=0; hcount<3; hcount++) {
            for (i=0; i<in->nbnd; i++) {
                for (j=0; j<pheight; j++) {
                    for (k=0; k<pwidth; k++) {
                        section->data[i][j][k] = in->data[i][j+(voffset*vcount)][k+(hoffset*hcount)] * window->data[0][j][k];;
                    }
                }
            }

            PSPECT(section,&scratch);

            for (i=0; i<(*out)->nbnd; i++) {
                for (j=0; j<(*out)->nlin; j++) {
                    for (k=0; k<(*out)->npix; k++) {
                        (*out)->data[i][j][k] += scratch->data[i][j][k];
                    }
                }
            }

        }
    }
    
    
    for (i=0; i<(*out)->nbnd; i++) {
        for (j=0; j<(*out)->nlin; j++) {
            for (k=0; k<(*out)->npix; k++) {
                (*out)->data[i][j][k] /= (PIXEL)9.0;
            }
        }
    }
    
    
    
    //RADIALAVG(*out);
    
    


    the_end:
    if (window) RELMEM(window);
    if (periodogram) RELMEM(periodogram);
    if (scratch) RELMEM(scratch);
    /* if (*out) RELMEM(*out); */
    if (TIMES)   TIMING(T_EXIT);
    
NAMES = 1;
}

---