---

sconvl.c

Go to the documentation of this file.

#include        "sadie.h"
#include        "math.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1988 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function computes a two-dimensional spatial convolution between     */
/*   a user-supplied spatial kernel and a single-band image.                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
SCONVL (IMAGE * in,             /*  I   Pointer to the input image.                           */
        PIXEL * psf,            /*  I   Pointer to the convolution mask matrix[jsize][ksize]. */
        short jsize,            /*  I   Number of lines in the convolution mask.              */
        short ksize,            /*  I   Number of pixels/line in the convolution mask.        */
        IMAGE ** out            /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
  )
{
  register short j, k, l, m, n, jmarg = jsize / 2, kmarg = ksize / 2, box;
  char msg[SLEN];
  double sum, c1, cn, pinc, psum;

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

  progmeter = ProgMeter_Create ("Convolution");
  progress = 0.0;

  if (TIMES)
    TIMING (T_SCONVL);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "SCONVL");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:           %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Kernel size: lines:    %d", jsize);
      MESSAGE ('I', msg);
      sprintf (msg, "              pixels:   %d", ksize);
      MESSAGE ('I', msg);
      if (jsize <= SLEN / 16 && ksize <= SLEN / 16)
        {
          MESSAGE ('I', " Spatial kernel:");
          for (n = j = 0; j < jsize; j++)
            {
              for (m = k = 0, msg[m++] = ' '; k < ksize;
                   k++, m = strlen (msg))
                {
                  sprintf (msg + m, "%12.4e", (double) psf[n++]);
                }
              MESSAGE ('I', msg);
            }
        }
      else
        {
          MESSAGE ('I', " Kernel too large to list.");
        }
      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;
    }
  else if (jsize < 0 || jsize % 2 == 0)
    {
      MESSAGE ('E',
               " Number of lines in the convolution mask must be an odd number greater than zero.");
      goto the_end;
    }
  else if (jsize > in->nlin)
    {
      MESSAGE ('E',
               " Image size must be equal to or greater than convolution mask size in the line dimension.");
      goto the_end;
    }
  else if (ksize < 0 || ksize % 2 == 0)
    {
      MESSAGE ('E',
               " Number of pixels/line in the convolution mask must be an odd number greater than zero.");
      goto the_end;
    }
  else if (ksize > in->npix)
    {
      MESSAGE ('E',
               " Image size must be equal to or greater than convolution mask size in the pixel dimension.");
      goto the_end;
    }

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

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc = 1.0 / (double) (in->nlin - 2 * jmarg);
  progtotal = in->nlin - 2 * jmarg;
  progcount = 0;
  progincr = rint ((double) progtotal / 20.0);

  /* see if it's a box-filter PSF */
  for (box = TRUE, k = 1; box && k < jsize * ksize; k++)
    {
      box = psf[k] == psf[0];
    }

  if (box)
    {
      for (j = jmarg; j < in->nlin - jmarg; j++)
        {

          /* initialize box */
          for (c1 = sum = 0.0, l = j - jmarg; l <= j + jmarg; l++)
            {
              c1 += (double) in->data[0][l][0];
              for (k = 1; k < ksize; k++)
                {
                  sum += (double) in->data[0][l][k];
                }
            }
          (*out)->data[0][j][kmarg] = (PIXEL) ((sum += c1) * (double) psf[0]);

          /* box-filter algorithm */
          for (k = kmarg + 1; k < in->npix - kmarg; k++)
            {
              for (sum -= c1, c1 = cn = 0.0, l = j - jmarg; l <= j + jmarg;
                   l++)
                {
                  c1 += (double) in->data[0][l][k - kmarg];
                  cn += (double) in->data[0][l][k + kmarg];
                }
              (*out)->data[0][j][k] = (PIXEL) ((sum += cn) * (double) psf[0]);
            }
          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;
                }
            }
        }

    }
  else                          /* compute normal convolution */
    {
      for (j = jmarg; j < in->nlin - jmarg; j++)
        {
          for (k = kmarg; k < in->npix - kmarg; k++)
            {
              for (sum = 0.0, n = jsize * ksize, l = 0; l < jsize; l++)
                {
                  for (m = 0; m < ksize; m++)
                    {
                      sum +=
                        (double) in->data[0][j - jmarg + l][k - kmarg +
                                                            m] *
                        (double) psf[--n];
                    }
                }
              (*out)->data[0][j][k] = (PIXEL) sum;
            }
          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 < jmarg; j++)
    {
      for (k = 0; k < in->npix - kmarg; k++)
        {
          (*out)->data[0][j][k] = (*out)->data[0][jmarg][k];
          (*out)->data[0][in->nlin - j - 1][k] =
            (*out)->data[0][in->nlin - jmarg - 1][k];
        }
    }

  /* fill left and right margins */
  for (j = 0; j < in->nlin; j++)
    {
      for (k = 0; k < kmarg; k++)
        {
          (*out)->data[0][j][k] = (*out)->data[0][j][kmarg];
          (*out)->data[0][j][in->npix - k - 1] =
            (*out)->data[0][j][in->npix - kmarg - 1];
        }
    }

the_end:
  if (TIMES)
    TIMING (T_EXIT);

  if (progmeter)
    {
      ProgMeter_Destroy (progmeter);
    }

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

---