 /*
  * Khoros: $Id: lvfft.c,v 1.4 1992/01/17 00:43:44 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvfft.c,v 1.4 1992/01/17 00:43:44 dkhoros Exp $";
#endif

 /*
  * $Log: lvfft.c,v $
 * Revision 1.4  1992/01/17  00:43:44  dkhoros
 * HellPatch4
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvfft.c
 >>>>
 >>>>      Program Name: vfft
 >>>>
 >>>> Date Last Updated: Wed Jan 15 10:00:35 1992 
 >>>>
 >>>>          Routines: lvfft - the library call for vfft
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvfft - library call for vfft
*
* Purpose:
*    
*    Computes the two dimensional Fast Fourier Transform (FFT)
*    
*    

* Input:
*    
*    image1         pointer to xvimage structure of type FLOAT or COM-
*                   PLEX.
*    
*    image2         pointer  to  xvimage  structure  of   type   FLOAT
*                   (optional).
*    
*    dir            integer  specifying  the  direction  of  the   FFT
*                   transform. 0 is forward, 1 is inverse.
*    
*    

* Output:
*    
*    c_flg          is an integer requesting  that  a  complex  output
*                   image  be  produiced.   Non-zero means produce the
*                   image.
*    
*    o_image1       is a pointer to the image structure  pointer  that
*                   will  receive  the  complex output image, if it is
*                   desired (see above).
*    
*    r_flg          is an integer requesting that a real output  image
*                   be produced.  Non-zero means produce the image.
*    
*    o_image2       is a pointer to the image structure  pointer  that
*                   will  receive  the  real  output  image,  if it is
*                   desired (see above).
*    
*    i_flg          is an integer requesting that a  imaginary  output
*                   image  be  produced.   Non-zero  means produce the
*                   image.
*    
*    o_image3       is a pointer to the image structure  pointer  that
*                   will  receive the imaginary output image, if it is
*                   desired (see above).
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Scott Wilson
*    
*    16-Jan-91 C.Gage -  Modified code to accept  either  a  FLOAT  or
*    COMPLEX  image  and  output a COMPLEX image, or the REAL, or IMA-
*    GINARY parts.  Routine was also modified to  perform  FORWARD  or
*    INVERSE  FFT.  Added  call  to powtwo() to check for proper image
*    size.
*    
*    8-Mar-91 Scott Wilson - Pixels sizes were not  being  transferred
*    from input to output.
*    
*    25-Aug-91 Scott Wilson - Added information about scaling and fre-
*    quency shifting to the man pages.
*    
*    12-Nov-91 Scott Wilson - Corrected memory leaks in  lvfft.  Added
*    comment about freeing the input images.
*    
*    15-Jan-92 Scott Wilson - Fixed routine so it does  not  free  the
*    input images!
*    
*    

****************************************************************/


/* -library_def */
int
lvfft(image1, image2, c_flg, o_image1, r_flg, o_image2, i_flg, o_image3, dir)
struct xvimage *image1, *image2, **o_image1, **o_image2, **o_image3;
int c_flg, r_flg, i_flg;
int dir;
/* -library_def_end */

/* -library_code */
{
    struct xvimage *c_image, *r_image, *i_image, *createimage();

    int i,j,k,nr,nc,ic,ir;      /* Indexes and sizes of rows and columns */
    int type;                   /* FFT type argument for FORTRAN call */
    float *pr, *pi;             /* Pointers to complex pairs array */

    float *fptr,                /* pointer to image1 data */
          *imptr,               /* pointer to image2 (imaginary) data */
          *xr,*xi,*zr,*zi,*tr,*ti;

    int m,m1;                    /* Stuff for Nyquist modulation */

    float         *cptr,        /* pointer to complex image data */
                  *rptr,        /* pointer to real image data */
                  *iptr;        /* pointer to imaginary image data */

    char  *program = "lvfft";

    nr = image1->number_of_rows;            /* Number of rows */
    nc = image1->number_of_cols;            /* Number of columns */

    /* Check to make sure input image is of correct type */
    if ( (!propertype(program, image1, VFF_TYP_COMPLEX, FALSE)) && (!propertype(program, image1, VFF_TYP_FLOAT, FALSE)) )
    {
      (void) fprintf(stderr, "lvfft: Input image must be of data type COMPLEX or FLOAT.\n");
      return(0);
    }

    /* Check for proper map enable */
    if (!proper_map_enable(program,image1,VFF_MAP_OPTIONAL,FALSE))
    {
      (void) fprintf(stderr, "lvfft: Can only work on map scheme NONE\n");
      return(0);
    }

    /* Check to see if optional second input image matches first input image */
    if (image2 != NULL)
    {
      if ( (!match_num_images(program, image1, image2, FALSE)) && (!match_num_bands(program, image1, image2, FALSE)) && (!match_map_enable(program, image1, image2, FALSE)) )
      {
        (void) fprintf(stderr,"lvfft: Both input images must be of same dimension.\n");
        return(0);
      }
      if ( (!matchsize(program, image1, image2, FALSE)) && (!matchtype(program, image1, image2, FALSE)) )
      {
        (void) fprintf(stderr,"lvfft: Both input images must be of same size and type.\n");
        return(0);
      }
    }

    /* Make sure size is legal for an FFT */
    if ( (powtwo(nc) == 0) || (powtwo(nr) == 0) )
    {
      (void)fprintf(stderr,"lvfft: Input image is %d rows by %d cols.\n", nr, nc);
      (void)fprintf(stderr,"Image size MUST be a power of 2!\n ");
      return(0);
    }

    /* Get space for the intermediate complex arrays */
    xr = (float *)malloc(nc*nr*sizeof(float));
    xi = (float *)malloc(nc*nr*sizeof(float));
    zr = (float *)malloc((nc+2)*sizeof(float));
    zi = (float *)malloc(nc*sizeof(float));

    if (xr == NULL || xi == NULL || zr == NULL || zi == NULL)
      {
        fprintf(stderr,"lvfft: Could not allocate enough memory!\n");
        return(0);
      }

    /****   Create the COMPLEX output image ****/
  if (c_flg)
  {
    c_image = createimage(nr,                  /* number of rows in image */
                          nc,                  /* number of cols in image */
                          VFF_TYP_COMPLEX,     /* data storage type */
                          1,                   /* number of images */
                          1,                   /* number of data bands */
                          "created by vfft", /* comment */
                          0,                   /* map row size */
                          0,                   /* map col size */
                          VFF_MS_NONE,         /* map scheme */
                          VFF_MAPTYP_NONE,
                          VFF_LOC_IMPLICIT,
                          0);
    c_image->pixsizx = image1->pixsizx;
    c_image->pixsizy = image1->pixsizy;

    if (c_image == NULL)
      {
        (void) fprintf(stderr,"lvfft: Unable to allocate enough memory for output image!\n");
        return(0);
      }
    cptr = (float *) c_image->imagedata;
  }

    /**** Create the FLOAT output image corresponding to REAL part ****/
  if (r_flg)
  {
    r_image = createimage(nr,                  /* number of rows in image */
                          nc,                  /* number of cols in image */
                          VFF_TYP_FLOAT,       /* data storage type */
                          1,                   /* number of images */
                          1,                   /* number of data bands */
                          "created by vfft", /* comment */
                          0,                   /* map row size */
                          0,                   /* map col size */
                          VFF_MS_NONE,         /* map scheme */
                          VFF_MAPTYP_NONE,
                          VFF_LOC_IMPLICIT,
                          0);
    r_image->pixsizx = image1->pixsizx;
    r_image->pixsizy = image1->pixsizy;

    if (r_image == NULL)
      {
        (void) fprintf(stderr,"lvfft: Unable to allocate enough memory for output image!\n");
        return(0);
      }
    rptr = (float *) r_image->imagedata;
  }

    /**** Create the FLOAT output image corresponding to IMAG part ****/
  if (i_flg)
  {
    i_image = createimage(nr,                  /* number of rows in image */
                          nc,                  /* number of cols in image */
                          VFF_TYP_FLOAT,       /* data storage type */
                          1,                   /* number of images */
                          1,                   /* number of data bands */
                          "created by vfft", /* comment */
                          0,                   /* map row size */
                          0,                   /* map col size */
                          VFF_MS_NONE,         /* map scheme */
                          VFF_MAPTYP_NONE,
                          VFF_LOC_IMPLICIT,
                          0);
    i_image->pixsizx = image1->pixsizx;
    i_image->pixsizy = image1->pixsizy;

    if (i_image == NULL)
      {
        (void) fprintf(stderr,"lvfft: Unable to allocate enough memory for output image!\n");
        return(0);
      }
    iptr = (float *) i_image->imagedata;
  }

             /* assign pointer to input image data */
  if (image2 == NULL)          /* optional second input image not present */
  {
    fptr = (float *)(image1->imagedata);
  }
  else if (image2 != NULL)
  {
    fptr = (float *)(image1->imagedata);
    imptr = (float *)(image2->imagedata);
  }

 if (!dir)   /*<<<<<<<<    Forward FFT, dir =  0 (FALSE)    >>>>>>>>>>>*/
 {
    type = 0;                    /* Set FFT type argument to FORWARD */
    m1 = 1;                      /* Initialize Nyquist modulator */
    for(i=0; i<nr; i++)          /* Do Nyquist modulate and FFT on each row */
      {
        m = m1; m1 = -m1;
        if ( (image1->data_storage_type == VFF_TYP_FLOAT) && (image2 == NULL) )
        {
          for (j=0; j<nc; j++)
          {
            zr[j] = m*(*fptr++);
            zi[j] = 0.0;         /* set imaginary part to 0 if a FLOAT image */
            m = -m;
          }
        }
        else if ( (image1->data_storage_type == VFF_TYP_FLOAT) && (image2 != NULL) )
        {
          for (j=0; j<nc; j++)
          {
             zr[j] = m*(*fptr++);   /* get real part from first input image */
             zi[j] = m*(*imptr++);  /* get imaginary part from second input image */
             m = -m;
          }
        }
        else
        {            
          for (j=0; j<nc; j++)
          {
             zr[j] = m*(*fptr++);   /* get real and imag parts from complex image */
             zi[j] = m*(*fptr++);
             m = -m;
          }
        }
         /* Call Fortran routine to do full complex FFT */
         fft842_(&type,&nc,zr,zi);

         pr = xr+i*nc;                 /* Store into (complex) rows) */
         pi = xi+i*nc;
        for(j=0; j<nc; j++)
        {
           *pr++ = zr[j];   /* Store the real part */
           *pi++ = zi[j];   /* Store the imag part */
        }
      }

    free(zr);               /* Give back temp row space */
    free(zi);               /* Give back temp row space */

    /* Get space for temporary complex row */
    tr = (float *)malloc(nr*sizeof(float));
    ti = (float *)malloc(nr*sizeof(float));
    if (tr == NULL || ti == NULL)
      {
        fprintf(stderr,"lvfft: Couldn't allocate enough memory!\n");
        return(0);
      }

    for(i=0; i<nc; i++)           /* Do FFT for the columns */
      {
        /* Load a column into temporary row */
        pr = xr+i;
        pi = xi+i;
        for (j=0; j<nr; j++)
          {
            tr[j] = *pr; pr += nc;
            ti[j] = *pi; pi += nc;
          } 

        /* Call Fortran routine to do FFT for the columns */
        fft842_(&type,&nr,tr,ti);

        /* Put back into the same column */
        pr = xr+i;
        pi = xi+i;
        for (j=0; j<nr; j++)
          {
            *pr = tr[j]; pr += nc;
            *pi = ti[j]; pi += nc;
          }
      }

    /* Load the data in corresponding output images */
    pr = xr;
    pi = xi;

    if (c_flg)           /* output COMPLEX image */
    {
      for (i=0; i<nr; i++)
      {
        for (j=0; j<nc; j++)
        {
          *cptr++ = *pr++;
          *cptr++ = *pi++;
        }
      }
    }

    if (r_flg)      /* output FLOAT image corresponding to REAL part */
    {
      for (i=0; i<nr; i++)
        for (j=0; j<nc; j++)
          *rptr++ = *pr++;
    }

    if (i_flg)      /* output FLOAT image corresponding to IMAG part */
    {
      for (i=0; i<nr; i++)
        for (j=0; j<nc; j++)
          *iptr++ = *pi++;
    }
    free(xr);               /* free real and imag data arrays */
    free(xi);
    free(tr);
    free(ti);
 }
 else if (dir)   /*<<<<<<<<<<<<    Inverse FFT, dir = TRUE    >>>>>>>>>>*/
 {
    type = 1;                  /* Set FFT type argument to INVERSE */
    /* Load input data from corresponding input images */
    pr = xr;               /* set pointer to real part of image */
    pi = xi;               /* set pointer to imag part of image */
    for (i=0; i<nr; i++)
    {
      if ( (image1->data_storage_type == VFF_TYP_FLOAT) && (image2 == NULL) )
      {
        for (j=0; j<nc; j++)
        {
          *pr++ = *fptr++;
          *pi++ = 0.0;        /* set imag part to zero */
        }
      }
      else if ( (image1->data_storage_type == VFF_TYP_FLOAT) && (image2 != NULL) )
      {
        for (j=0; j<nc; j++)
        {
          *pr++ = *fptr++;
          *pi++ = *imptr++;        /* get imag part from second input image */
        }
      }
      else            /* input image is COMPLEX */
      {
        for (j=0; j<nc; j++)
        {
          *pr++ = *fptr++;
          *pi++ = *fptr++;
        }
      }
    }

    /* Get space for the temporaries */
    tr = (float *)malloc(nr*sizeof(float));
    ti = (float *)malloc(nr*sizeof(float));

    if (tr == NULL || ti == NULL)
      {
        fprintf(stderr,"lvfft: Couldn't allocate enough memory!\n");
        return(0);
      }

    for(i=0; i<nc; i++)    /* Do IFFT for the columns */
    {
        /* Load a column into temporary row */
        pr = xr+i;
        pi = xi+i;
        for (j=0; j<nr; j++)
          {
            tr[j] = *pr; pr += nc;
            ti[j] = *pi; pi += nc;
          } 
        fft842_(&type,&nr,tr,ti);

        /* Put back into the same column */
        pr = xr+i;
        pi = xi+i;
        for(j=0; j<nr; j++)
          {
            *pr = tr[j]; pr += nc;
            *pi = ti[j]; pi += nc;
          }
    }

    free(tr);               /* free temporary arrays */
    free(ti);

    m1 = 1;
    for(i=0; i<nr; i++)          /* Do FFT for each row in image */
    {
        pr = xr+i*nc;             /* Store into (complex) rows) */
        pi = xi+i*nc;
        for (j=0; j<nc; j++)
        {
           zr[j] = *pr++;     /* Get the real part */
           zi[j] = *pi++;   /* Get the imag part */
        }
        fft842_(&type,&nc,zr,zi); /* Do complex IFFT */

        /* Store resulting image */
        m = m1; m1 = -m1;

        if (c_flg)              /* output COMPLEX image */
        {
          for (j=0; j<nc; j++)
          {
            *cptr++ = m*zr[j];
            *cptr++ = m*zi[j];
            m = -m;
          }
        }

        if (r_flg)              /* output FLOAT image for REAL part */
        {
          for (j=0; j<nc; j++)
          {
            *rptr++ = m*zr[j];
            m = -m;
          }
        }

        if (i_flg)              /* output FLOAT image for IMAG part */
        {
          for (j=0; j<nc; j++)
          {
            *iptr++ = m*zi[j];
            m = -m;
          }
        }
    }
    free(xr);               /* free real and imag data arrays */
    free(xi);
    free(zr);
    free(zi);
 }

    /* Assign data to new image pointers */
   if (c_flg) *o_image1 = c_image;

   if (r_flg) *o_image2 = r_image;

   if (i_flg) *o_image3 = i_image;

    return(1);

}  /*** End of main ***/
/* -library_code_end */
