---

sconvl.c

Go to the documentation of this file.

#include "trees.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 TREES_SCONVL (
IMAGE *in,      /*  I   Pointer to the input image.                           */
LOCAL_INDEX *localindex,                 
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 int j, k, l, m, n, jmarg=jsize/2, kmarg=ksize/2, box;
    char   msg[SLEN];
    double sum, r1, rn, 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 output image */
    for(j = 0; j < in->nlin; j++) {
        for(k = 0; k < in->npix; k++) {
            (*out)->data[0][j][k] = (PIXEL)0;
        }
    }

    /* 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) { /* Box filter implementation is yet to be tested */
        for (k = localindex->hstart + kmarg; k <= localindex->hend - kmarg; k++) {

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

            /* box-filter algorithm */
            for (j = localindex->vstart[k] + 1; j <= localindex->vend[k]; j++) {
                for (sum -= r1, r1 = rn = 0.0, l = k - kmarg; l <= k + kmarg; l++) {
                    r1 += (double)in->data[0][j - jmarg][l];
                    rn += (double)in->data[0][j + jmarg][l];
                }
                (*out)->data[0][j][k] = (PIXEL)((sum += rn)*(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 (k = localindex->hstart; k <= localindex->hend; k++) {
            for (j = localindex->vstart[k]; j <= localindex->vend[k]; j++) {
                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;
                    }
            }
        }
    }

    
    the_end:
    if (TIMES) TIMING(T_EXIT);

    if (progmeter) {
        ProgMeter_Destroy(progmeter);
    }

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

---