---

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);
}

---