---

gradient.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 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 TREES_GRADIENT (
IMAGE *in,      /*  I   Pointer to the input image.                           */
LOCAL_INDEX *localindex,                   
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 int 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 output image */
    for(j = 0; j < in->nlin; j++) {
        for(k = 0; k < in->npix; k++) {
            (*out1)->data[0][j][k] = (PIXEL)0;
            (*out2)->data[0][j][k] = (PIXEL)0;
        }
    }

    /* 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);

    /* Compute convolution */
    for (k = localindex->hstart+n; k <= localindex->hend-n; k++) {
        for (j = localindex->vstart[k]+n; j <= localindex->vend[k]-n; j++) {
            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[0][j-n+l][k-n+m];
                    dh += (double)horz[l*size+m]*(double)in->data[0][j-n+l][k-n+m];
                }
            }
            (*out1)->data[0][j][k] = (PIXEL)sqrt(dv*dv+dh*dh)/1443.0*255.0;
            if (dh == 0.0) {
                if (dv > 0.0) {
                    (*out2)->data[0][j][k] = (PIXEL)(PI/2.0);
                } else if (dv < 0.0) {
                    (*out2)->data[0][j][k] = (PIXEL)(-PI/2.0);
                } else /* (dv == 0.0 ) */ {
                    (*out2)->data[0][j][k] = (PIXEL)0.0;
                }
            } else {
                (*out2)->data[0][j][k] = (PIXEL)atan2(dv,dh);
            }
            (*out2)->data[0][j][k] = (unsigned char)(128 + (int)(127.0*(*out2)->data[0][j][k]/PI));
        }
        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 (*out1) RELMEM(*out1);
        *out1 = NULL;
        if (*out2) RELMEM(*out2);
        *out2 = NULL;
    }
}

---