newfct.c

基于java的3d开发库。对坐java3d的朋友有很大的帮助。

/***************************************************************************
  **************************************************************************
  
                Spherical Harmonic Transform Kit 2.7
  
  
   Contact: Peter Kostelec
            [e-mail]
  
  
   Copyright 1997-2003  Sean Moore, Dennis Healy,
                        Dan Rockmore, Peter Kostelec
  
  
   Copyright 2004  Peter Kostelec, Dan Rockmore


     SpharmonicKit is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
     (at your option) any later version.
  
     SpharmonicKit is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.
  
     You should have received a copy of the GNU General Public License
     along with this program; if not, write to the Free Software
     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  
  
   Commercial use is absolutely prohibited.
  
   See the accompanying LICENSE file for details.
  
  ************************************************************************
  ************************************************************************/

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

  If FFTPACK *is* defined, use the FFTPACK-based dcts. Small
  C-interfaces have been written, to act as front ends to the
  functions which live in the fftpack library that I'll be linking
  to.

  NOTE THAT even if FFTPACK is defined, I still need the inverse
  DCT routine ExpIFCT. That is because, in the inverse discrete
  Legendre transform, I need to take N coefficients and inverse
  DCT them to produce 2*N samples. I don't know how to do that
  with FFTPACK's inverse dct. If I have to do an inverse DCT whose
  input and output are the same lengths, *then* I can use FFTPACK's
  routine.

  The "standard" (i.e. non-FFTPACK) forward and inverse dct
  routines defined in this file are

  forward dct: kFCT, kFCTX
  inverse dct: ExpIFCT

  When using FFTPACK, the routines are

  forward dct: DCTf    ("f" for forward)
  inverse dct: DCTb    ("b" for backward)

  If you have your own optimized dct routines, just modify
  the definitions of DCTf and DCTb accordingly (and still
  compile with the -DFFTPACK flag and change the library
  you're linking to in the Makefile).

  NOTE: DCTf and DCTb use slightly modified versions of
        the fftpack functions cosqb ( in DCTf ) and cosqf (in DCTb ).
	Before compiling the fftpack library, these functions
	were modified in order to scale the coefficients the same way
	as they are in kFCT, kFCTX, and ExpIFCT. Originally,
	the transforms cosqb and cosqf were not orthogonal.
	THEY ARE NOW.

	Exactly how was the scaling modified in the original,
	unmodified fftpack routines cosqb, cosqf ?

	Let X be an input array of length N. Let Y be the output
	of kFCT (or kFCTX - they're identical), Z be the
	output of the original, unmodified fftpack routine cosqb.
	Let SCALE = 1 / (2 * N) .Then (using C and not Fortran indexing)

	Y[0] = SCALE * Z[0] / 2
	Y[i] = SCALE * Z[i]    for i = 1, 2, ..., N-1

	The original fftpack library routine cosqb was further
	modified to return only the first p <= N many
	coefficients. This required modifying the original,
	unmodified fftpack-source file cosqb.f .


	Let X be an input array of length N. Let Y be the output
	of ExpIFCT. Now let me multiply the first coefficient of
	X by 2 ( so X[0] *= 2 ) and plug this modified X into
	the original, unmodified fftpack-routine cosqf.
	Let Z be the resulting output of cosqf. Let SCALE = 1/2.
	Then (using C and not Fortran indexing)

	Y[i] = SCALE * Z[i]    for i = 0, 1, 2, ..., N-1


	All this scaling is now done within the source files
	cosqf1.f and cosqb1.f . It's a little hokey, but it
	works.


  IMPORTANT NOTE: DCTf and DCTb both write the output
                  over the input array!!!

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

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

  Some general background on the routines kFCT, kFCTX and ExpIFCT ...


  This is new C cource code for doing fast cosine transforms, and
  it contains variations for speeding up FCTs and convolutions
  for certain structures appearing in Legendre transforms.
  Many of these commands have both a real-valued version and
  a complex-valued version to save time when performing 
  spherical harmonic transforms.
  It also is written with memory management in mind.  
  
  This cosine transform is based on the Steidl and Tasche algorithm
  with some modifications.  See Sean's Ph.D. thesis for refs.
  

  NOTE: I will try to explicitly unwrap the ChebyDivRem routine
  as called in ExpIFCT: the while-loop will be for divsize > 4.
  Then I'll explicitly type in for the cases of divsize = 4 and
  divsize = 2. This, I believe, significantly reduces the overhead
  involved in calling ChebyDivRem so often. Thank you Doug and Sumit !


  Note about structures: some functions use the lowhigh structure,
  defined to be
  
  struct lowhigh{ double low; double high; }
  
  The reason for this is to "interweave" elements. That is, if
  I have a vector A of length 2*n, then I'll create an array with
  n-many lowhigh elements which will look like
  
  {{A[0],A[n]},{A[1],A[n+1]},{A[2],A[n+2]},...,{A[n-1],A[2*n]}}
  
  Having the elements arranged this way increases in a small
  (but real) way the efficiency of the forward fast cosine transform.
  So, for example, kFCTX is the "structure" version of kFCT. The
  structure is well suited for the transform, given how this
  particular algorithm is implemented. Basically, I believe,
  using this structure cuts down on cache-thrashing.
  
  General note about flags - throughout this code flags are often
  used as input arguments to tell the routine something about
  the input information.  The convention adopted throughout is
  that when a flag is set to 0, that means that no action needs
  to be performed on the pertinent data, and 1 means that some action
  must be taken.
  

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



#include <math.h>
#include <string.h>

#include "OURperms.h"
#include "OURmods.h"
#ifdef FFTPACK
#include "fftpack.h"
#endif
#include "newFCT.h"


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

  Utility functions for kFCT, kFCTX and ExpIFCT

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

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

  reverse the elements of a (double) vector of length n
  i.e. if a = [0 1 2 3] then what's returned is
          b = [3 2 1 0]

	  input -> vector to reverse
	  n     -> length of input (must be even !)
	  output -> where to store result

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

static void reverse_it(double *input,
		       double *output,
		       int n)
{
  int i, len;

  len = n / 2;

  if (len >= 4)
    {
      for(i = 0; i < (len % 4); ++i)
	{
	  output[i] = input[n - 1 - i];
	  output[n - 1 - i] = input[i];
	}
      for( ; i < len ; i += 4)
	{
	  output[i] = input[n - i - 1];
	  output[i + 1] = input[n - i - 2];
	  output[i + 2] = input[n - i - 3];
	  output[i + 3] = input[n - i - 4];

	  output[n - i - 1] = input[i];
	  output[n - i - 2] = input[i + 1];
	  output[n - i - 3] = input[i + 2];
	  output[n - i - 4] = input[i + 3];
	}
    }
  else if (len == 2)
    {
      output[0] = input[n - 1];
      output[1] = input[n - 2];
      output[n - 1] = input[0];
      output[n - 2] = input[1];
    }
  else if (len == 3)
    {
      output[0] = input[n - 1];
      output[1] = input[n - 2];
      output[2] = input[n - 3];
      output[n - 1] = input[0];
      output[n - 2] = input[1];
      output[n - 3] = input[2];
    }
  else
    {
      output[0] = input[n - 1];
      output[n - 1] = input[0];
    }

}

/************************************************************************/
/* permutation management code */
/* OUR permutations - a self-inverting permutation used for
   performing fast cosine transforms.  See Sean's Ph.D. thesis */
/************************************************************************/
/* permutes data and puts the result in result */
/* data and result should be arrays of size n */

static void OURpermute(double *data,
		       double *result,
		       int n)
{
  int i;
  const int *pptr;

  pptr = get_perm(n);
  
  for (i=0; i<n; i++)
    result[i] = data[pptr[i]];
  
}
/************************************************************************/
/* permutes data and puts the result in result */
/* data is size n, result array of structures of size n/2 */
/* just like above except designed to work with structures */

static void OURpermuteX(double *data,
			struct lowhigh *result,
			int n)
{
  register int i, dummy;
  const int *pptr;

  dummy = n / 2;

  pptr = get_perm(n);

  for ( i = 0 ; i < dummy; i ++)
    {
      result->low = data[*pptr];
      result->high = data[*(pptr+dummy)];
      result++; pptr++;
    }
}

/************************************************************************/
/* used by kFCT to add up polynomial coefficients once k level is reached */
/* data is double array of size n, result is double array of size k */

static void PartitionAdd(double *data,
			 double *result,
			 int n,
			 int k)
{
  int i, j;
  double tmp;

  for(i = 0; i < k ; i++)
    {
      tmp = 0.0;
      for(j = i; j < n; j += k)
	tmp += data[j];
      result[i] = tmp;
    }
}

/************************************************************************/
/* used by kFCT to add up polynomial coefficients once k level is reached */
/* data is double array of size n, result is double array of size k */
/* just like above except used for lowhigh structures ...
 ASSUMING k = (n/2) */

static void PartitionAddX(struct lowhigh *data,
			  double *result,
			  int n,
			  int k)
{
  int i, j;
  double tmp0, tmp1;
  int dummy;

  dummy = n / 2;

  for (i = 0 ; i < k ; i++)
    {
      tmp0 = 0.0 ; tmp1 = 0.0;
      for(j = i ; j < dummy ; j += k)
	{
	  tmp0 += data[j].low;
	  tmp1 += data[j].high;
	}
      result[i] = ( tmp0 + tmp1 );
    }
}
/************************************************************************/

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

  Ok, if FFTPACK *is* defined, use DCTf, DCTb. Otherwise use
  kFCT, kFCTX

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

#ifdef FFTPACK


/****

  If will use the fftpack routines cosqf and cosqb, first
  have to precompute an array of values that both routines
  use. This interface function will call the fftpack routine
  
       cosqi

  This function will allocate memory to store the precomputed
  values and then will return a double pointer to that array.

  NOTE: The memory allocated in this function is freed in
        the routine that calls precomp_dct!!!

  
  size = length of sequence to be transformed
  
  wSave = the array created by precomp_dct, of length
          3 * size + 15, which contains the precomputed
	  values. precomp_dct will return a double pointer
	  to this array.


  ****/

double *precomp_dct( int size )
{
  double *wSave;

  /* allocate space */
  wSave = (double *) malloc( sizeof(double) * (3 * size + 15) );

  /* now precompute data */
  cosqi_( &size, wSave );

  return wSave;

}


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

  This is the forward DCT routine DCTf. It uses the fftpack library
  function cosqb (NOT cosqf).

  array = double array of samples, of length n

  p = how many coefficients desired.

  CoswSave = array of precomputed data that cosqb requires (produced
             by a call to the function cosqi ).

  OUTPUT IS WRITTEN OVER INPUT !!!

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

void DCTf ( double *array,
	    int n,
	    int p,
	    double *CoswSave )
{
  
  cosqb_( &n, &p, array, CoswSave );

}

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

  This is the inverse DCT routine DCTb. It uses the fftpack library
  function cosqf (NOT cosqb).

  array = double array of coefficients, of length n;
          n many samples will be produced

  CoswSave = array of precomputed data that cosqf requires (produced
             by a call to the function cosqi ).

  
  OUTPUT IS WRITTEN OVER INPUT !!!

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

void DCTb( double *array ,
	   int n,
	   double *CoswSave )
{
  cosqf_( &n , array , CoswSave );

}

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

#else   /* FFTPACK is *not* defined, use the "local" functions kFCT, kFCTX */

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

/************************************************************************/
/* OK - here is the forward transform code.  The main modification here
   is that given n samples of data, kFCT will compute only the first
   k <= n coefficients, where n and k are assumed to be powers of 2.
   The basic idea is that the algorithm keeps interpolating a polynomial
   of degree 2d-1 from two polynomials of degree d-1.  Initially, lowpoly
   contains polynomials of degree 0 - the data samples.  From these
   values, polynomials of degree 1 are interpolated and stored in highpoly.
   This transform is (essentially) the transpose of the matrix formulation
   of ExpIFCT.  Series is returned with coefficients ordered from low
   degree to high degree

   data - double array of size n
   result - double array of size k
   workspace - double array of size 2*n
   n - number of samples
   k = number of coefficients desired 
   permflag - 0 if data in OUR permuted order, 1 if data needs to 
              be permuted */

void kFCT(double *data,
	  double *result,
	  double *workspace,
	  int n,
	  int k,
	  int permflag)
{
    double *lowpoly, *highpoly, *lowptr, *highptr;
    const double *modptr;
    double *temp, dn2;
    int p2, p3;
    unsigned int i, j, levelsize, highsize, lowsize, loopcntr;
    double modptr0;
    double e0, e1, f0, f1;

    lowpoly = workspace; 
    lowptr = lowpoly;
    highpoly = workspace + n; 
    highptr = highpoly;
    
    levelsize = n;
    modptr = get_mods(levelsize);

    if (permflag == 1)
      OURpermute(data,lowpoly,n);
    else {
      memcpy(lowpoly, data, sizeof(double) * n);
    }

    /***** in what follows, am using the fact
      that modptr[i+1] = -modptr[i] for i even ****/

    /* do first level of interpolation  */
    for (i=0; i<n; i+=2) {
	*highptr++ = *lowptr + *(lowptr+1);
	*highptr++ = *modptr * (*(lowptr+1) - *lowptr);
	modptr+=2;lowptr+=2;
    }

    lowptr = highpoly;
    highptr = lowpoly;
    loopcntr = 0;

    /* now do higher-order interpolations */
    if (k > 2)
      {
	/* reset pointers to logical polynomial workspace */	
	levelsize /= 2;
	modptr = get_mods(levelsize);
	highsize = 4;
	lowsize = 2;

	/* now set counter pointers to appropriate matrix locations */
	p2 = lowsize - 1;
	p3 = 1;
	
	/* main loop */
	while (highsize <= k)
	  {
	    for (j=0; j<n; j+=highsize)
	      {
		/* do loworder coeffs first */
		for (i=0; i<lowsize; i++)
		  *highptr++ = lowptr[i] + lowptr[i+lowsize];

		/* now do higher order coeffs */
		
		modptr0 = modptr[0];
		*highptr++ = modptr0 * (lowptr[lowsize] - lowptr[0]);
		modptr0 *= 2.0;
		
		e0 = modptr0 * ( lowptr[p3+lowsize] - lowptr[p3] );
		e1 = lowptr[p2+lowsize] + lowptr[p2];
		*highptr++ = e0 - e1;
		p2--;
		p3++;
		
		for ( i = 0 ; i<(lowsize-2); i += 2)
		  {
		    e0 = -lowptr[p3] + lowptr[p3+lowsize];
		    f0 = -lowptr[p3+1] + lowptr[p3+lowsize+1];
		    f1 = lowptr[p2-1] + lowptr[p2+lowsize-1];
		    e1 = lowptr[p2] + lowptr[p2+lowsize];		    
		    *highptr++ = modptr0 * e0 - e1;
		    *highptr++ = modptr0 * f0 - f1;
		    p2 -= 2;
		    p3 += 2;
		  }

		lowptr += highsize;
		modptr += 2;
		p2 = lowsize - 1;
		p3 = 1;
	      }