---

sconvl.c

Go to the documentation of this file.

#include        "sadie.h"
#include        "math.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function computes a two-dimensional spatial convolution between     */
/*   a user-supplied spatial kernel and a single-band image.                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void SCONVL (
IMAGE *in,      /*  I   Pointer to the input image.                           */
PIXEL *psf,     /*  I   Pointer to the convolution mask matrix[jsize][ksize]. */
short jsize,    /*  I   Number of lines in the convolution mask.              */
short ksize,    /*  I   Number of pixels/line in the convolution mask.        */
IMAGE **out     /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
) { register short j, k, l, m, n, jmarg=jsize/2, kmarg=ksize/2, box;
    char   msg[SLEN];
    double sum, c1, cn, pinc, psum;

    int progmeter;
    double progress;
    int progcount, progtotal, progincr;
    int cancel=0;

    progmeter = ProgMeter_Create("Convolution");
    progress = 0.0;


    if (TIMES) TIMING(T_SCONVL);
    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"SCONVL");
        MESSAGE('I',"");
        sprintf(msg," Input image:           %s",in->text);
        MESSAGE('I',msg);
        sprintf(msg," Kernel size: lines:    %d",jsize);
        MESSAGE('I',msg);
        sprintf(msg,"              pixels:   %d",ksize);
        MESSAGE('I',msg);
        if (jsize <= SLEN/16  &&  ksize <= SLEN/16) {
            MESSAGE('I'," Spatial kernel:");
            for (n=j=0; j<jsize; j++) {
                for (m=k=0,msg[m++]=' '; k<ksize; k++,m=strlen(msg)) {
                    sprintf(msg+m,"%12.4e",(double)psf[n++]);
                }
                MESSAGE('I',msg);
            }
        } else {
            MESSAGE('I'," Kernel too large to list.");
        }
        MESSAGE('I'," ...............");
    }

    /* check input */
    if (!CHECKIMG(in)) {
        MESSAGE('E'," Can't identify image.");
        goto the_end;
    } else if (in->nbnd > 1) {
        MESSAGE('E'," Image must be single-band.");
        goto the_end;
    } else if (jsize < 0  ||  jsize%2 == 0) {
        MESSAGE('E'," Number of lines in the convolution mask must be an odd number greater than zero.");
        goto the_end;
    } else if (jsize > in->nlin) {
        MESSAGE('E'," Image size must be equal to or greater than convolution mask size in the line dimension.");
        goto the_end;
    } else if (ksize < 0  ||  ksize%2 == 0) {
        MESSAGE('E'," Number of pixels/line in the convolution mask must be an odd number greater than zero.");
        goto the_end;
    } else if (ksize > in->npix) {
        MESSAGE('E'," Image size must be equal to or greater than convolution mask size in the pixel dimension.");
        goto the_end;
    }

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

    /* initialize progress indicator */
    if (LINES  &&  !PROGRESS(psum=0.0)) goto the_end;
    pinc = 1.0/(double)(in->nlin-2*jmarg);
    progtotal = in->nlin-2*jmarg;
    progcount = 0;
    progincr = rint((double)progtotal/20.0);

    /* see if it's a box-filter PSF */
    for (box=TRUE,k=1; box && k<jsize*ksize; k++) {
        box = psf[k] == psf[0];
    }
    
    if (box) {
        for (j=jmarg; j<in->nlin-jmarg; j++) {

            /* initialize box */
            for (c1=sum=0.0,l=j-jmarg; l<=j+jmarg; l++) {
                c1 += (double)in->data[0][l][0];
                for (k=1; k<ksize; k++) {
                    sum += (double)in->data[0][l][k];
                }
            }
            (*out)->data[0][j][kmarg] = (PIXEL)((sum+=c1)*(double)psf[0]);

            /* box-filter algorithm */
            for (k=kmarg+1; k<in->npix-kmarg; k++) {
                for (sum-=c1,c1=cn=0.0,l=j-jmarg; l<=j+jmarg; l++) {
                    c1 += (double)in->data[0][l][k-kmarg];
                    cn += (double)in->data[0][l][k+kmarg];
                }
                (*out)->data[0][j][k] = (PIXEL)((sum+=cn)*(double)psf[0]);
            }
            if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
            if (progmeter) {
                    progcount++;
                    if ((progcount%progincr) == 0) {
                        progress = ((double)progcount/(double)progtotal)*100.0;
                        cancel = ProgMeter_Update(progmeter,progress);
                        if (cancel != 0) goto the_end;
                    }
            }
        }

    } else /* compute normal convolution */ {
        for (j=jmarg; j<in->nlin-jmarg; j++) {
            for (k=kmarg; k<in->npix-kmarg; k++) {
                for (sum=0.0,n=jsize*ksize,l=0; l<jsize; l++) {
                    for (m=0; m<ksize; m++) {
                        sum += (double)in->data[0][j-jmarg+l][k-kmarg+m]*(double)psf[--n];
                    }
                }
                (*out)->data[0][j][k] = (PIXEL)sum;
            }
            if (LINES  &&  !PROGRESS(psum+=pinc)) goto the_end;
            if (progmeter) {
                    progcount++;
                    if ((progcount%progincr) == 0) {
                        progress = ((double)progcount/(double)progtotal)*100.0;
                        cancel = ProgMeter_Update(progmeter,progress);
                        if (cancel != 0) goto the_end;
                    }
            }
        }
    }

    /* fill top and bottom margins */
    for (j=0; j<jmarg; j++) {
        for (k=0; k<in->npix-kmarg; k++) {
            (*out)->data[0][j][k]            = (*out)->data[0][jmarg][k];
            (*out)->data[0][in->nlin-j-1][k] = (*out)->data[0][in->nlin-jmarg-1][k];
        }
    }

    /* fill left and right margins */
    for (j=0; j<in->nlin; j++) {
        for (k=0; k<kmarg; k++) {
            (*out)->data[0][j][k]            = (*out)->data[0][j][kmarg];
            (*out)->data[0][j][in->npix-k-1] = (*out)->data[0][j][in->npix-kmarg-1];
        }
    }

    the_end:
    if (TIMES) TIMING(T_EXIT);

    if (progmeter) {
        ProgMeter_Destroy(progmeter);
    }

    if (cancel != 0) {
        if (*out) RELMEM(*out);
        *out = NULL;
    }
}

---