---

treeringmag.c

Go to the documentation of this file.

#include "sadie.h"

/*
static PIXEL firstdiff_y[3][3] = { { 0.0,  0.0,  0.0 }, 
                                   { 0.0,  1.0,  1.0 }, 
                                   { 0.0, -1.0, -1.0 } };

static PIXEL firstdiff_x[3][3] = { { 0.0,  0.0,  0.0 }, 
                                   { 0.0, -1.0,  1.0 }, 
                                   { 0.0, -1.0,  1.0 } };
*/

/* Use sobel filter */

static PIXEL firstdiff_x[3][3] = { { -1.0,  0.0,  1.0 }, 
                                   { -2.0,  0.0,  2.0 }, 
                                   { -1.0,  0.0,  1.0 } };

static PIXEL firstdiff_y[3][3] = { { 1.0,  2.0,  1.0 }, 
                                   { 0.0,  0.0,  0.0 }, 
                                   {-1.0, -2.0, -1.0 } };


/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/* This function detects tree rings.                                          */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void TREERINGMAG (
IMAGE  *in,      /*  I   Pointer to the input image.                          */
double sigma,    /*  I   Std. deviation of gaussian smoothing filter.         */
int    size,     /*  I   Size of gaussian smoothing filter.                   */
IMAGE  **out     /*  O   Address of a pointer to the output image             */
                 /*      containing an edge map.                              */
/*----------------------------------------------------------------------------*/
) { register short i, j, k, n;
    char msg[SLEN];
    IMAGE *gaussian, *smoothed;
    IMAGE *FoG_mag, *FoG_dir;
    PIXEL *psf=0, factor;
    PIXEL theta;
    int sectorval;
        

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

    progmeter = ProgMeter_Create("Compute Tree Ring Magnitudes");
    progress = 0.0;



    /* if (TIMES) TIMING(T_TREERINGMAG); */
    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"TREERINGMAG");
        MESSAGE('I',"");
        sprintf(msg," Input image:                       %s",in->text);
        MESSAGE('I',msg);
        sprintf(msg," Gaussian smoothing filter size:    %d",size);
        MESSAGE('I',msg);
        sprintf(msg,"             Standard Deviation:    %f",sigma);
        MESSAGE('I',msg);
        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;
    }


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


    if (size > 0) {
        /* create gaussian smoothing filter */
        CREATEGAUSS(size, sigma, &gaussian);
        if (CHECKIMG(gaussian)) {
            sprintf(gaussian->text, "%dx%d Gaussian -- sigma=%f", size,size,sigma);
        } else {
            MESSAGE('E'," Can't identify gaussian filter.");
            goto the_end;
        }
        if (progmeter) {
                progress = 10.0;
                cancel = ProgMeter_Update(progmeter,progress);
                if (cancel != 0) goto the_end;
        }

        /* perform gaussian smoothing on input image */
        if (!(psf = (PIXEL *) malloc(size*size*sizeof(PIXEL)))) {
            MESSAGE('E'," Memory allocation failed.");
            goto the_end;
        }
        for (n=0, j=0, factor=0.0; j<size; j++) {
            for (k=0; k<size; k++) {
                psf[n] = gaussian->data[0][j][k];
                factor += psf[n++];
            }
        }
        for (n=0, j=0; j<size; j++) {
            for (k=0; k<size; k++) {
                psf[n++] /= factor;
            }
        }
        SCONVL (in, psf, size, size, &smoothed);
    } else {
        IMGCOPY(in,&smoothed);
    }

    if (progmeter) {
            progress = 30.0;
            cancel = ProgMeter_Update(progmeter,progress);
            if (cancel != 0) goto the_end;
    }


    if (CHECKIMG(smoothed)) {
        sprintf(smoothed->text, "Smoothed image");
    } else {
        MESSAGE('E'," Can't identify smoothed image.");
        goto the_end;
    }
    
    /* create the First Order Gaussian filtered image */
    GRADIENT(smoothed,(PIXEL *)firstdiff_y,(PIXEL *)firstdiff_x,3,&FoG_mag,&FoG_dir);

    if (progmeter) {
            progress = 70.0;
            cancel = ProgMeter_Update(progmeter,progress);
            if (cancel != 0) goto the_end;
    }
    

    /* perform non-maxima supression, only detecting low to high transitions */
    if (progmeter) {
        localpercent = 0.3;
        localprogress = 0.0;
        progcount = 0;
        progtotal = in->npix * in->nlin;
        progincr = rint((double)progtotal/8.0);
    }
    for (j=0; j<in->nlin; j++) {
        for (k=0; k<in->npix; k++) {

            if (progmeter) {
                progcount++;
                if ((progcount%progincr) == 0) {
                    localprogress = ((double)progcount/(double)progtotal)*100.0;
                    cancel = ProgMeter_Update(progmeter,(localprogress*localpercent)+progress);
                    if (cancel != 0) goto the_end;
                }
            }

            if (((k>0)&&(FoG_mag->data[0][j][k-1] > FoG_mag->data[0][j][k]))
                || ((k<in->npix-1)&&(FoG_mag->data[0][j][k+1] > FoG_mag->data[0][j][k]))) {
                    /* Supress this non-maxima edge pixel */
                    (*out)->data[0][j][k] = 0.0;
            } else if ((k<in->npix-1)&&(k>0)
                &&(smoothed->data[0][j][k+1] < smoothed->data[0][j][k-1])) {
                    /* Supress this non-ring edge */
                    (*out)->data[0][j][k] = 0.0;
            } else {
                (*out)->data[0][j][k] = FoG_mag->data[0][j][k];
            }
        }
    }
    
        
    the_end:
    if (psf) free(psf);
    if (CHECKIMG(gaussian))  RELIMG(&gaussian);
    if (CHECKIMG(smoothed))  RELIMG(&smoothed);
    if (CHECKIMG(FoG_mag))  RELIMG(&FoG_mag);
    if (CHECKIMG(FoG_dir))  RELIMG(&FoG_dir);
    /* if (TIMES) TIMING(T_EXIT); */

    if (progmeter) {
        ProgMeter_Destroy(progmeter);
    }

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

---