---

gradient.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 a two-dimensional gradient (magnitude and         */
/*   direction) of an image, using two user-supplied convolution kernels.     */
/*   The magnitude is computed as the vector magnitude of the output          */
/*   of the two kernels, and the direction is computed as the angle           */
/*   between the two orthogonal gradient vectors.                             */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Robinson, G.S.:                                                          */
/*   "Detection and Coding of Edges Using Directional Masks."                 */
/*   Optical Engineering, Vol. 16, No. 6 (Nov/Dec 1977), pp. 580-585          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void GRADIENT (
IMAGE *in,      /*  I   Pointer to the input image.                           */
PIXEL *vert,    /*  I   Pointer to a square convolution kernel[size][size] for*/
                /*      the y-derivative. It should return positive values    */
                /*      for gradients increasing to the "top" of the image.   */
PIXEL *horz,    /*  I   Pointer to a square convolution kernel[size][size] for*/
                /*      the x-derivative. It should return positive values    */
                /*      for gradients increasing to the "right" of the image. */
short size,     /*  I   Kernel size in lines and pixels/line.                 */
IMAGE **out1,   /*  O   Address of a pointer to the output image              */
                /*      containing the gradient magnitude.                    */
IMAGE **out2    /*  O   Address of a pointer to the output image              */
                /*      containing the gradient direction.                    */
/*----------------------------------------------------------------------------*/
) { register short i, j, k, l, m, n=size/2;
    char   msg[SLEN];
    double dv, dh, pinc, psum;

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

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


    if (TIMES) TIMING(T_GRADIENT);
    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"GRADIENT");
        MESSAGE('I',"");
        sprintf(msg," Input image:                    %s",in->text);
        MESSAGE('I',msg);
        sprintf(msg," Kernel size: lines and pixels:  %d",size);
        MESSAGE('I',msg);
        if (size <= SLEN/16) {
            MESSAGE('I'," Horizontal kernel:");
            for (i=0; i<size; i++) {
                for (j=k=0,msg[k++]=' '; j<size; j+=1,k=strlen(msg)) {
                    sprintf(msg+k,"%12.4e",(double)horz[i*size+j]);
                }
                MESSAGE('I',msg);
            }
            MESSAGE('I'," Vertical kernel:");
            for (i=0; i<size; i++) {
                for (j=k=0,msg[k++]=' '; j<size; j+=1,k=strlen(msg)) {
                    sprintf(msg+k,"%12.4e",(double)vert[i*size+j]);
                }
                MESSAGE('I',msg);
            }
        } else {
            MESSAGE('I'," Kernels too large to list.");
        }
        MESSAGE('I'," ...............");
    }

    /* check input */
    if (!CHECKIMG(in)) {
        MESSAGE('E'," Can't identify image.");
        goto the_end;
    } else if (size < 2  ||  size%2 == 0) {
        MESSAGE('E'," Size of convolution mask must be an odd number greater than two.");
        goto the_end;
    } else if (size > in->nlin  ||  size > in->npix) {
        MESSAGE('E'," Image size must be equal to or greater than convolution mask size.");
        goto the_end;
    }

    /* create images of appropriate size */
    if (!CHECKIMG(*out1)) GETMEM(in->nbnd,in->nlin,in->npix,out1);
    if (!*out1) goto the_end;
    if (!CHECKIMG(*out2)) GETMEM(in->nbnd,in->nlin,in->npix,out2);
    if (!*out2) goto the_end;

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

    for (i=0; i<in->nbnd; i++) {

        /* compute convolution */
        for (j=n; j<in->nlin-n; j++) {
            for (k=n; k<in->npix-n; k++) {
                dv = dh = 0.0;
                for (l=0; l<size; l++) {
                    for (m=0; m<size; m++) {
                        dv += (double)vert[l*size+m]*(double)in->data[i][j-n+l][k-n+m];
                        dh += (double)horz[l*size+m]*(double)in->data[i][j-n+l][k-n+m];
                    }
                }
                (*out1)->data[i][j][k] = (PIXEL)sqrt((dv*dv+dh*dh));
                if (dh == 0.0) {
                    if (dv > 0.0) {
                        (*out2)->data[i][j][k] = (PIXEL)(PI/2.0);
                    } else if (dv < 0.0) {
                        (*out2)->data[i][j][k] = (PIXEL)(-PI/2.0);
                    } else /* (dv == 0.0 ) */ {
                        (*out2)->data[i][j][k] = (PIXEL)0.0;
                    }
                } else {
                    (*out2)->data[i][j][k] = (PIXEL)atan2(dv,dh);
                }
            }
            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<n; j++) {
            for (k=n; k<in->npix-n; k++) {
                (*out1)->data[i][j][k]            = (*out1)->data[i][n][k];
                (*out1)->data[i][in->nlin-j-1][k] = (*out1)->data[i][in->nlin-n-1][k];
                (*out2)->data[i][j][k]            = (*out2)->data[i][n][k];
                (*out2)->data[i][in->nlin-j-1][k] = (*out2)->data[i][in->nlin-n-1][k];
            }
        }

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

    the_end:
    if (TIMES) TIMING(T_EXIT);

    if (progmeter) {
        ProgMeter_Destroy(progmeter);
    }

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

---