---

fconvl.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function convolves an image with a filter function, using the       */
/*   Fast Fourier transform.                                                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void FCONVL (
IMAGE  *in,     /*  I   Pointer to the input image.                           */
IMAGE  *filter, /*  I   Pointer to the filter image.                          */
short  option,  /*  I   Switch, indicating the type of filter function:       */
                /*      PSF  -  spatial kernel.                               */
                /*      OTF  -  frequency transfer function.                  */
IMAGE  **out    /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
) { register short  i, j, k;
    char   msg[SLEN];
    short  jbgn, kbgn;
    double pinc, psum;
    IMAGE  *scratch1=0, *scratch2=0, *scratch3=0, *scratch4=0, *flt;

    if (TIMES) TIMING(T_FCONVL);
    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"FCONVL");
        MESSAGE('I',"");
        sprintf(msg," Input image:    %s",in->text);
        MESSAGE('I',msg);
        sprintf(msg," Filter image:   %s",filter->text);
        MESSAGE('I',msg);
        if (option == PSF) {
            MESSAGE('I'," Filter image is spatial kernel");
        } else if (option == OTF) {
            MESSAGE('I'," Filter image is frequency transfer function");
        }
        MESSAGE('I'," ...............");
    }

    /* check input */
    if (!CHECKIMG(in)) {
        MESSAGE('E'," Can't identify image.");
        goto the_end;
    } else if (in->nlin != 1L<<rnd(log10((double)in->nlin)/log10(2.0))) {
        MESSAGE('E'," Number of lines must be a power of two.");
        goto the_end;
    } else if (in->npix != 1L<<rnd(log10((double)in->npix)/log10(2.0))) {
        MESSAGE('E'," Number of pixels/line must be a power of two.");
        goto the_end;
    } else if (!CHECKIMG(filter)) {
        MESSAGE('E'," Can't identify filter image.");
        goto the_end;
    } else if (in->nbnd != filter->nbnd) {
        MESSAGE('E'," Number of bands must be identical in both images.");
        goto the_end;
    }
    if (option == PSF) {
        if (in->nlin < filter->nlin) {
            MESSAGE('E'," Image size in lines must be equal to or greater than that of filter image.");
            goto the_end;
        } else if (in->npix < filter->npix) {
            MESSAGE('E'," Image size in pixels/line must be equal to or greater than that of filter image.");
            goto the_end;
        }
    } else /* (option == OTF) */ {
        if (in->nlin != filter->nlin) {
            MESSAGE('E'," Image size in lines in the frequency domain must be equal to that of filter image.");
            goto the_end;
        } else if (in->npix != filter->npix/2) {
            MESSAGE('E'," Image size in pixels/line in the frequency domain must be equal to that of filter image.");
            goto the_end;
        }
    }
    
    /* Fourier transform input image */
    FFT2D(in,FORWARD,&scratch1);
    if (ABORT  ||  !scratch1) goto the_end;
    sprintf(scratch1->text,"Scratch");

    if (option == PSF) {
        /* Zero-pad (if necessary) and Fourier transform point-spread function */
        if (in->nlin > filter->nlin  ||  in->npix > filter->npix) {
            jbgn = (in->nlin-filter->nlin)/2;
            kbgn = (in->npix-filter->npix)/2;
            MOSAIC(&filter,1,in->nlin,in->npix,&jbgn,&kbgn,(PIXEL)0,&scratch2);
            if (ABORT  ||  !scratch2) goto the_end;
            sprintf(scratch2->text,"Scratch");
            FFT2D(scratch2,FORWARD,&scratch3);
            if (ABORT  ||  !scratch3) goto the_end;
            sprintf(scratch3->text,"Scratch");
            flt = scratch3;
        } else {
            FFT2D(filter,FORWARD,&scratch2);
            if (ABORT  ||  !scratch2) goto the_end;
            sprintf(scratch2->text,"Scratch");
            flt = scratch2;
        }
    } else /* (option == OTF) */ {
        flt = filter;
    }  

    /* create image of appropriate size */
    GETMEM(scratch1->nbnd,scratch1->nlin,scratch1->npix,&scratch4);
    if (!scratch4) goto the_end;
    sprintf(scratch4->text,"Scratch");

    /* initialize progress indicator */
    if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
    pinc = 1.0/(double)scratch1->nbnd/(double)scratch1->nlin;

    /* multiply two complex images */
    for (i=0; i<scratch1->nbnd; i+=1) {
        for (j=0; j<scratch1->nlin; j+=1) {
            for (k=0; k<scratch1->npix; k+=2) {
                scratch4->data[i][j][k  ] = (PIXEL)((double)scratch1->data[i][j][k]*(double)flt->data[i][j][k] - (double)scratch1->data[i][j][k+1]*(double)flt->data[i][j][k+1]);
                scratch4->data[i][j][k+1] = (PIXEL)((double)scratch1->data[i][j][k+1]*(double)flt->data[i][j][k] + (double)scratch1->data[i][j][k]*(double)flt->data[i][j][k+1]);
            }
            if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
        }
    }
    
    /* Fourier transform filtered image */
    FFT2D(scratch4,INVERSE,out);

    the_end:
    if (scratch1) RELMEM(scratch1);
    if (scratch2) RELMEM(scratch2);
    if (scratch3) RELMEM(scratch3);
    if (scratch4) RELMEM(scratch4);
    if (TIMES)    TIMING(T_EXIT);
}

---