 /*
  * Khoros: $Id: ldhart1d.c,v 1.2 1992/03/20 23:27:26 dkhoros Exp $
  */

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

 /*
  * $Log: ldhart1d.c,v $
 * Revision 1.2  1992/03/20  23:27:26  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: ldhart1d.c
 >>>>
 >>>>      Program Name: dhart1d
 >>>>
 >>>> Date Last Updated: Thu Mar  5 10:38:57 1992 
 >>>>
 >>>>          Routines: ldhart1d - the library call for dhart1d
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#define FSWAP(a,b) {float tempr;tempr=(a);(a)=(b);(b)=tempr;}
#define TWOPI     (double)6.283185307179586
#define PI        (double)3.14159265359
#define SQRT2     (double)1.414213562
#define NO_SCALE         0
#define SCALE_BY_N       1
#define SCALE_BY_SQRT_N  2
#define CENTER           1
#define NO_CENTER        0
#define FORWARD          1
#define INVERSE          0
/* -library_includes_end */


/****************************************************************
*
* Routine Name: ldhart1d - library call for dhart1d
*
* Purpose:
*    
*    Performs  the  1D  forward  Hartley  transform  of  a   data
*    sequence.
*    
*    

* Input:
*    
*    image          contains the signal data to be transformed.
*    
*    trans_dir      specifies the transform  direction.   0  indicates
*                   the  inverse  transform,  1  indicates the forward
*                   transform.  This is only useful when center_opt is
*                   set to 1 (center).
*    
*    center_opt     centering choice: 1 for centered, 0 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.
*    
*    process_dir    processing direction: 0 for vector processing or 1
*                   for band oriented processing.
*    
*    

* Output:
*    
*    image          holds the result of the transform.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Jeremy Worley
*    
*    Jeremy Worley 25 Oct 1990 11:09 MDT
*              Changed internal function log2  to  k_log2  to  prevent
*              compile errors.
*    
*    Jeremy Worley 05 Mar 1992 10:36 MST
*              Fixed up numerous implicit declarations of functions in
*              ldhart1d()  and  hart1d().   Fixed  a couple of fprintf
*              calls that were incorrect. Changed  permute()'s  return
*              type to void.
*    
*    

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


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

/* -library_code */
{
char *program = "ldhart1d";
  double k_log2();
  float         **data,                 /* array returned by dload_vector */
                scale_factor;

  int           dimension,              /* number of points per data sequence */                logdim,                 /* log of dimension */
                num_vects,              /* number of data sequences */
                i, j,                   /* loop variables */
                lvrtoc(),
                ldshift(),
                powtwo(),
                hart1d(),
                dunload_vector();

  char          **dload_vector();

 /*
  * The dhart1d routine works only on real data
  */

  if (image->data_storage_type != VFF_TYP_FLOAT){
      (void)fprintf(stderr,"%s: invalid data type!\n",program);
      return(0);
  }

 /*
  * Shift data if its the inverse and its necessary.
  */
  if ( center_opt == CENTER && trans_dir == INVERSE)
  {
        if(process_dir==DSP_BAND)
           dimension = image->row_size * image->col_size;
        else
           dimension = image->num_data_bands;
        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 */

 /*
  * Get data into 2-dimensional array vector format.
  */
  if ( ( data = (float **)dload_vector(image, &num_vects, &dimension,
                process_dir)) == NULL )
  {
    (void)fprintf(stderr,"%s: dload_vector failed\n",program);
    return(0);
  }

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

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

 /*
  * Perform the forward transform.
  */

  logdim = (int)k_log2((double)dimension);

  for (i=0; i<num_vects; i++)
  {
    (void) hart1d ( data[i], dimension, logdim );
  }

 /*
  * 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<dimension; j++)
           {
             data[i][j] /= scale_factor;
           }
         }
         break;

    case SCALE_BY_SQRT_N:
         scale_factor = sqrt((double)dimension);
         for (i=0; i<num_vects; i++)
         {
           for (j=0; j<dimension; 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_FLOAT, 
        num_vects, dimension, process_dir) )
  {
    (void)fprintf(stderr,"%s: dunload_vector failed!!\n",program);
    return(0);
  }

 /*
  * Shift data if its the forward and its necessary.
  */
  if ( center_opt == CENTER && trans_dir == FORWARD)
  {
        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 */

  for (i=0; i<num_vects; i++)
  {
    free((float *)data[i]);
  }
  free ((float**)data);

  return(1);
}


/*
** the real thing
*/

int hart1d(data,n,m)
  float *data;
  int n,m;
{
  int n2,n4,n8,id,is,i,j,k,l11,l12,l13,l14,l21,l22,l23,l24;
  float t1,t2,t3,t4,t5,a,a3,e,cc1,cc3,ss1,ss3;
  void permute();

  n2 = 2*n;
  for(k=1;k<m;k++){
      is = 1;
      id = n2;
      n2 = n2/2;
      n4 = n2/4;
      n8 = n2/8;
      e  = TWOPI/n2;
      do{
         for(i=is;i<=n;i+=id){                          /* 40 */
             l12 = i + n4;
             l13 = l12 + n4;
             l14 = l13 + n4;
             t1 = data[l12 - 1] - data[l14 - 1];
             t2 = data[i - 1] + data[l13 - 1];
             t3 = data[l12 - 1] + data[l14 - 1];
             t4 = data[i - 1] - data[l13 - 1];
             data[l14 - 1] = t4 - t1;
             data[l13 - 1] = t4 + t1;
             data[l12 - 1] = t3;
             data[i - 1] = t2;
             if(n4!=1){                                 /* rest of loop */
                 l21 = i + n8;
                 l22 = l21 + n4;
                 l23 = l22 + n4;
                 l24 = l23 + n4;
                 t1 = data[l22 - 1];
                 t2 = data[l23 - 1];
                 t3 = data[l24 - 1];
                 data[l24 - 1] = (t1 - t3)*SQRT2;
                 data[l23 - 1] = (data[l21 - 1] - t2)*SQRT2;
                 data[l22 - 1] = t1 + t3;
                 data[l21 - 1] = data[l21 - 1] + t2;
                 a = e;

                 for(j=2;j<=n8;j++){
                     l11 = i + j - 1;
                     l12 = l11 + n4;
                     l13 = l12 + n4;
                     l14 = l13 + n4;
                     l21 = i + n4 - j + 1;
                     l22 = l21 + n4;
                     l23 = l22 + n4;
                     l24 = l23 + n4;
                     a3 = 3*a;
                     cc1 = cos(a);
                     ss1 = sin(a);
                     cc3 = cos(a3);
                     ss3 = sin(a3);
                     a = j*e;
                     t5 = data[l21 - 1] - data[l23 - 1];
                     t2 = data[l11 - 1] - data[l13 - 1];
                     t1 = t2 + t5;
                     t2 = t2 - t5;
                     t5 = data[l22 - 1] - data[l24 - 1];
                     t4 = data[l14 - 1] - data[l12 - 1];
                     t3 = t4 + t5;
                     t4 = t4 - t5;
                     data[l11 - 1] = data[l11 - 1] + data[l13 - 1];
                     data[l12 - 1] = data[l12 - 1] + data[l14 - 1];
                     data[l21 - 1] = data[l21 - 1] + data[l23 - 1];
                     data[l22 - 1] = data[l22 - 1] + data[l24 - 1];
                     data[l13 - 1] = t1*cc1 + t3*ss1;
                     data[l14 - 1] = t2*cc3 - t4*ss3;
                     data[l23 - 1] = t1*ss1 - t3*cc1;
                     data[l24 - 1] = t2*ss3 + t4*cc3;
                 } /* end for j */
              } /* end if */
          } /* end for i */
   
          is = 2*id - n2 + 1;
          id = 4*id;

      }while(is<n); 

  } /* end for k */
 
/*
** last stage of transform 
*/

  is = 1;
  id = 4;
  do {
      for(i=is;i<=n;i+=id){
          t1 = data[i - 1];
          data[i - 1] = t1 + data[i];
          data[i]     = t1 - data[i];
      }
      is = 2*id - 1;
      id = 4*id;
  }while(is<n);

  permute(data,n);

  return(1);
}

/*
** permutation of data
*/

void permute(data,n)
  float *data;
  int n;
{
  int j,i,m;

        j = 1;
        for (i=1;i<n;i++)  {
                if (j > i)  {
                        FSWAP(data[j-1],data[i-1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m)  {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
}

double k_log2(x)
  double x;
{
  return(log10(x)*3.32192809489);
}

/* -library_code_end */
