---

fconvl.c

Go to the documentation of this file.

#include        "sadie.h"

/*-Copyright Information------------------------------------------------------*/
/* Copyright (c) 1990 by the University of Arizona Digital Image Analysis Lab */
/*----------------------------------------------------------------------------*/
/*-General Information--------------------------------------------------------*/
/*                                                                            */
/*   This function convolves an image with a filter function, using the       */
/*   Fast Fourier transform.                                                  */
/*                                                                            */
/*----------------------------------------------------------------------------*/
/*-Interface Information------------------------------------------------------*/
void
FCONVL (IMAGE * in,             /*  I   Pointer to the input image.                           */
        IMAGE * filter,         /*  I   Pointer to the filter image.                          */
        short option,           /*  I   Switch, indicating the type of filter function:       */
        /*      PSF  -  spatial kernel.                               */
        /*      OTF  -  frequency transfer function.                  */
        IMAGE ** out            /*  O   Address of a pointer to the output image.             */
/*----------------------------------------------------------------------------*/
  )
{
  register short i, j, k;
  char msg[SLEN];
  short jbgn, kbgn;
  double pinc, psum;
  IMAGE *scratch1 = 0, *scratch2 = 0, *scratch3 = 0, *scratch4 = 0, *flt;

  if (TIMES)
    TIMING (T_FCONVL);
  if (NAMES)
    {
      MESSAGE ('I', "");
      MESSAGE ('I', "FCONVL");
      MESSAGE ('I', "");
      sprintf (msg, " Input image:    %s", in->text);
      MESSAGE ('I', msg);
      sprintf (msg, " Filter image:   %s", filter->text);
      MESSAGE ('I', msg);
      if (option == PSF)
        {
          MESSAGE ('I', " Filter image is spatial kernel");
        }
      else if (option == OTF)
        {
          MESSAGE ('I', " Filter image is frequency transfer function");
        }
      MESSAGE ('I', " ...............");
    }

  /* check input */
  if (!CHECKIMG (in))
    {
      MESSAGE ('E', " Can't identify image.");
      goto the_end;
    }
  else if (in->nlin != 1L << rnd (log10 ((double) in->nlin) / log10 (2.0)))
    {
      MESSAGE ('E', " Number of lines must be a power of two.");
      goto the_end;
    }
  else if (in->npix != 1L << rnd (log10 ((double) in->npix) / log10 (2.0)))
    {
      MESSAGE ('E', " Number of pixels/line must be a power of two.");
      goto the_end;
    }
  else if (!CHECKIMG (filter))
    {
      MESSAGE ('E', " Can't identify filter image.");
      goto the_end;
    }
  else if (in->nbnd != filter->nbnd)
    {
      MESSAGE ('E', " Number of bands must be identical in both images.");
      goto the_end;
    }
  if (option == PSF)
    {
      if (in->nlin < filter->nlin)
        {
          MESSAGE ('E',
                   " Image size in lines must be equal to or greater than that of filter image.");
          goto the_end;
        }
      else if (in->npix < filter->npix)
        {
          MESSAGE ('E',
                   " Image size in pixels/line must be equal to or greater than that of filter image.");
          goto the_end;
        }
    }
  else                          /* (option == OTF) */
    {
      if (in->nlin != filter->nlin)
        {
          MESSAGE ('E',
                   " Image size in lines in the frequency domain must be equal to that of filter image.");
          goto the_end;
        }
      else if (in->npix != filter->npix / 2)
        {
          MESSAGE ('E',
                   " Image size in pixels/line in the frequency domain must be equal to that of filter image.");
          goto the_end;
        }
    }

  /* Fourier transform input image */
  FFT2D (in, FORWARD, &scratch1);
  if (ABORT || !scratch1)
    goto the_end;
  sprintf (scratch1->text, "Scratch");

  if (option == PSF)
    {
      /* Zero-pad (if necessary) and Fourier transform point-spread function */
      if (in->nlin > filter->nlin || in->npix > filter->npix)
        {
          jbgn = (in->nlin - filter->nlin) / 2;
          kbgn = (in->npix - filter->npix) / 2;
          MOSAIC (&filter, 1, in->nlin, in->npix, &jbgn, &kbgn, (PIXEL) 0,
                  &scratch2);
          if (ABORT || !scratch2)
            goto the_end;
          sprintf (scratch2->text, "Scratch");
          FFT2D (scratch2, FORWARD, &scratch3);
          if (ABORT || !scratch3)
            goto the_end;
          sprintf (scratch3->text, "Scratch");
          flt = scratch3;
        }
      else
        {
          FFT2D (filter, FORWARD, &scratch2);
          if (ABORT || !scratch2)
            goto the_end;
          sprintf (scratch2->text, "Scratch");
          flt = scratch2;
        }
    }
  else                          /* (option == OTF) */
    {
      flt = filter;
    }

  /* create image of appropriate size */
  GETMEM (scratch1->nbnd, scratch1->nlin, scratch1->npix, &scratch4);
  if (!scratch4)
    goto the_end;
  sprintf (scratch4->text, "Scratch");

  /* initialize progress indicator */
  if (LINES && !PROGRESS (psum = 0.0))
    goto the_end;
  pinc = 1.0 / (double) scratch1->nbnd / (double) scratch1->nlin;

  /* multiply two complex images */
  for (i = 0; i < scratch1->nbnd; i += 1)
    {
      for (j = 0; j < scratch1->nlin; j += 1)
        {
          for (k = 0; k < scratch1->npix; k += 2)
            {
              scratch4->data[i][j][k] =
                (PIXEL) ((double) scratch1->data[i][j][k] *
                         (double) flt->data[i][j][k] -
                         (double) scratch1->data[i][j][k +
                                                       1] *
                         (double) flt->data[i][j][k + 1]);
              scratch4->data[i][j][k + 1] =
                (PIXEL) ((double) scratch1->data[i][j][k + 1] *
                         (double) flt->data[i][j][k] +
                         (double) scratch1->data[i][j][k] *
                         (double) flt->data[i][j][k + 1]);
            }
          if (LINES && !PROGRESS (psum += pinc))
            goto the_end;
        }
    }

  /* Fourier transform filtered image */
  FFT2D (scratch4, INVERSE, out);

the_end:
  if (scratch1)
    RELMEM (scratch1);
  if (scratch2)
    RELMEM (scratch2);
  if (scratch3)
    RELMEM (scratch3);
  if (scratch4)
    RELMEM (scratch4);
  if (TIMES)
    TIMING (T_EXIT);
}

---