---

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

---