---

relief.c

Go to the documentation of this file.

#include "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1993 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function calculates a shaded relief image from digital terrain data */
/*   or any other surface data.                                               */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Background Information-----------------------------------------------------*/
/*                                                                            */
/*   Woodham, Robert J.:                                                      */
/*   "Using digital terrain data to model image formation in remote sensing"  */
/*   SPIE, Vol. 238 (1988), pp. 361-369                                       */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void RELIEF (
IMAGE  *in,     /*  I   Pointer to the input image.                           */
double phi,     /*  I   Horizontal illumination angle,                        */
                /*      measured in radians counterclockwise from east.       */
                /*      Note - for standard NCIC terrain data,                */
                /*      east appears at top of the image.  In this case,      */
                /*      phi for a due north azimuth would be PI/2.            */
double theta,   /*  I   Vertical illumination angle,                          */
                /*      measured in radians from the horizon.                 */
double dist,    /*  I   Distance between adjacent data points,                */
                /*      measured in the same units as the elevation data.     */
                /*      Note - for standard NCIC terrain data,                */
                /*      this value is 208.3 feet.                             */
IMAGE  **out    /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
) { register short i, j, k;
    char msg[SLEN];
    double ps, qs, sn, p, q;

    if (TIMES) TIMING(T_RELIEF);
    if (NAMES) {
        MESSAGE('I',"");
        MESSAGE('I',"RELIEF");
        MESSAGE('I',"");
        sprintf(msg," Input image:                       %s",in->text);
        MESSAGE('I',msg);
        sprintf(msg," Horizontal illumination angle:   %12.4e",180.0*phi/PI);
        MESSAGE('I',msg);
        sprintf(msg," Vertical illumination angle:     %12.4e",180.0*theta/PI);
        MESSAGE('I',msg);
        sprintf(msg," Distance between pixels:         %12.4e",dist);
        MESSAGE('I',msg);
        MESSAGE('I'," ...............");
    }

    /* check input */
    if (!CHECKIMG(in)) {
        MESSAGE('E'," Can't identify image.");
        goto the_end;
    } else if (phi < -2.0*PI  ||  2.0*PI < phi) {
        MESSAGE('E'," Horizontal illumination angle must be between -360 and +360 degrees.");
        goto the_end;
    } else if (theta < -0.5*PI  ||  0.5*PI < theta) {
        MESSAGE('E'," Vertical illumination angle must be between -90 and +90 degrees.");
        goto the_end;
    }

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

    /* initialize irradiance parameters */
    ps = sin(phi) / tan(theta);
    qs = cos(phi) / tan(theta);
    sn = sqrt(1.0 + ps*ps + qs*qs);

    /* build relief image */
    for (i=0; i<in->nbnd; i++) {
        for (j=1; j<in->nlin-1; j++) {
            for (k=1; k<in->npix-1; k++) {
                p = (double)(in->data[i][j][k+1] - in->data[i][j][k-1]);
                q = (double)(in->data[i][j-1][k] - in->data[i][j+1][k]);
                (*out)->data[i][j][k] = (PIXEL)((dist + ps*p + qs*q) / (sn * sqrt(dist*dist + p*p + q*q)));
            }
        }
    }

    /* zero borders */
    for (i=0; i<in->nbnd; i++) {
        for (j=0; j<in->nlin; j++) {
            (*out)->data[i][j][         0] = (PIXEL)0;
            (*out)->data[i][j][in->npix-1] = (PIXEL)0;
        }
        for (k=0; k<in->npix; k++) {
            (*out)->data[i][         0][k] = (PIXEL)0;
            (*out)->data[i][in->nlin-1][k] = (PIXEL)0;
        }
    }
        
    the_end:
    if (TIMES) TIMING(T_EXIT);
}

---