 /*
  * Khoros: $Id: ldfft1d.c,v 1.2 1992/03/20 23:25:31 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: ldfft1d.c,v 1.2 1992/03/20 23:25:31 dkhoros Exp $";
#endif

 /*
  * $Log: ldfft1d.c,v $
 * Revision 1.2  1992/03/20  23:25:31  dkhoros
 * VirtualPatch5
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1992, 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 1992 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: ldfft1d.c
 >>>>
 >>>>      Program Name: dfft1d
 >>>>
 >>>> Date Last Updated: Thu Mar  5 10:24:01 1992 
 >>>>
 >>>>          Routines: ldfft1d - the library call for dfft1d
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define NO_SCALE         0
#define SCALE_BY_N       1
#define SCALE_BY_SQRT_N  2
#define CENTER           1
#define NO_CENTER        0
#define DFT_FORWARD      0
#define DFT_INVERSE      1
/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldfft1d - library call for dfft1d
*
* Purpose:
*    
*    Forward 1D Discrete Fourier Transform.
*    
*    

* Input:
*    
*    image          pointer to VIFF structure containing image data to
*                   be processed.
*    
*    center_opt     centering choice: 0 for centered, 1 for  not  cen-
*                   tered.
*    
*    scale_opt      scaling choice: 0 for no scaling, 1 for  multiply-
*                   ing  by 1/n, 2 for multiplying by 1/sqrt(n), where
*                   n is the number of data points.
*    
*    procdir        process direction:  0  indicated  vector  oriented
*                   processing, 1 indicates band oriented processing.
*    
*    

* Output:
*    
*    image          pointer to VIFF structure  containing  image  data
*                   after processing.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Ramiro Jordan and Donna Koechner
*    
*    Jeremy Worley 29 May 1990 15:50 MDT
*              fixed type cast problem in 1/sqrt() scaling option
*    
*    Jeremy Worley 07 Jun 1990 13:27 MDT
*              fixed call to lvrtoc() at the top of the library code
*    
*    Jeremy Worley 03 Feb 1992 10:45 MST
*              altered code to call fortran multi-radix fft in libvipl
*              for efficiency
*    
*    Jeremy Worley 19 Feb 1992 09:25 MST
*              Added    explicit    declarations     of     ldshift(),
*              dunload_vector(), and fft842_().
*    
*    Jeremy Worley 05 Mar 1992 10:21 MST
*              Changed an fprintf statement  so  that  the  number  of
*              arguments matched the format field in ldfft1d().
*    
*    

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


/* -library_def */
int
ldfft1d ( image, center_opt, scale_opt, process_dir )
struct xvimage  *image;
int             center_opt,
                scale_opt,
                process_dir;
/* -library_def_end */

/* -library_code */
{
  float         **data,                 /* array returned by dload_vector */
                scale_factor,
                *dr,*di;

  int           dimension,              /* number of points per data sequence */
                num_vects,              /* number of data sequences */
                i, j, k,                  /* loop variables */
                lvrtoc(),
                points,
                powtwo(),
                direction = DFT_FORWARD;

  char          **dload_vector(),
                *program = "ldfft1d";
  int           dunload_vector(),ldshift();
  void          fft842_();

 /*
  * The lfft1d routine works only on complex data, so if input data
  * is real, first convert the data to complex, where the imaginary part
  * is zero.
  */

  if(image->data_storage_type != VFF_TYP_COMPLEX){
    if(!(lvrtoc(image,NULL,1))){
       fprintf(stderr,"%s: lvrtoc failed!\n",program);
       return(0);
    }
  }

/*
 * Get data into 2-dimensional array vector format.
 */

  if((data = (float **)dload_vector(image, &num_vects, &dimension, 
                process_dir)) == NULL ) {
     fprintf(stderr,"%s: dload_vector failed\n",program);
     return(0);
  }

/*
 * if the data is complex, we will need to reorganize it.
 */

  if((dr = (float *)malloc(dimension*sizeof(float)))==NULL){
     fprintf(stderr,"%s: unable to allocate temporary data.\n",program);
     return(0);
  }

  if((di = (float *)malloc(dimension*sizeof(float)))==NULL){
     fprintf(stderr,"%s: unable to allocate temporary data.\n",program);
     return(0);
  }

  if(!dimension){
     fprintf(stderr,"ldfft1d: Number of points in sequence is zero!!\n");
     return(0);
  }

/* 
 * Verify that length of data sequence is a power of two.
 */

  if(!powtwo(dimension)){
     fprintf (stderr,"%s: Length of sequence must be ",program);
     fprintf (stderr,"a power of two value !!!\n");
     fprintf (stderr,"Sequence length is %d points\n",dimension);
     return(0);
  }

  points = dimension * 2;

/*
 * Perform the forward transform.  
 */

  for(i=0; i<num_vects; i++){
      for(j=0,k=0;j<points;j+=2,k++){
          dr[k] = data[i][j];
          di[k] = data[i][j+1];
      }
    
      (void) fft842_ ( &direction, &dimension,dr,di);

      for(j=0,k=0;j<points;j+=2,k++){
          data[i][j] = dr[k];
          data[i][j+1] = di[k];
      }
  }

/*
 * Scale output data.
 */

  switch(scale_opt)
  {
    case NO_SCALE:
         break;

    case SCALE_BY_N:
         scale_factor = dimension;
         for (i=0; i<num_vects; i++)
         {
           for (j=0; j<points; j++)
           {
             data[i][j] /= scale_factor;
           }
         }
         break;

    case SCALE_BY_SQRT_N:
         scale_factor = sqrt((float)dimension);
         for (i=0; i<num_vects; i++)
         {
           for (j=0; j<points; j++)
           {
             data[i][j] /= scale_factor;
           }
         }
         break;

    default:
         (void) fprintf(stderr,"%s: Unknown scaling option --> ",program);
         (void) fprintf(stderr,"output not scaled\n");
         break;
  }

 /*
  * Load data back into VIF format.
  */
  if (!dunload_vector((char**)data, image, (unsigned long)VFF_TYP_COMPLEX, 
        num_vects, dimension, process_dir) )
  {
    (void)fprintf(stderr,"%s: dunload_vector failed!!\n",program);
    return(0);
  }
 
  for (i=0; i<num_vects; i++)
  {
    free((float *)data[i]);
  }
  free ((float**)data);

 /*
  * Shift data if necessary.
  */
  if ( center_opt == CENTER )
  {
        if(!ldshift(image,(struct xvimage *)NULL,0,(int)(dimension/2),
                1,0.0,0.0,process_dir)){        
           fprintf(stderr,"%s:  shift routine failed.\n",program);
           return(0);
        }
  }  /* end of shift loop */

  return(1);

} /* end of ldfft1d */
/* -library_code_end */
