newfct.c

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

	    levelsize /= 2;
	    modptr = get_mods(levelsize);
	    highsize *= 2;
	    lowsize *= 2;
	    p2 = lowsize - 1;
	    p3 = 1;
	    
	    if ((loopcntr % 2) == 0)
	      {
		lowptr = lowpoly;
		highptr = highpoly;
	      }
	    else
	      {
		lowptr = highpoly;
		highptr = lowpoly;
	      }
	    loopcntr++;
	  }
      }
    
    /* find where the top level polynomials are */
    if (!(loopcntr % 2))
      temp = highpoly;
    else
      temp = lowpoly;
    
    if (k < n) 
      PartitionAdd(temp,result,n,k);
    else
      memcpy(result, temp, sizeof(double) * n);
    
    /* normalize */
    dn2 = 2.0/((double) n);
    result[0] /= ((double) n);
    for(i=1; i<k; i++)
      result[i] *= dn2;
        
    /* later */
}

/************************************************************************/
/* 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 interpoltaed 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
   lowpoly, highpoly - lowhigh arrays of length n / 2

 */

/*** just like kFCT EXCEPT will use structures ***/

void kFCTX( double *data,
	    double *result,
	    double *workspace,
	    int n,
	    int k,
	    int permflag,
	    struct lowhigh *lowpoly,
	    struct lowhigh *highpoly)
{
  double dn2;
  int p2, p3;
  register int i, j;
  int levelsize, highsize, lowsize, loopcntr;

  int dummy_int, dummy_int2;

  double e0, e1, f0, f1;
  double e0x, e1x, f0x, f1x;

  struct lowhigh *lowptr, *highptr, *temp;
  const double *modptrlow, *modptrhigh;
  register double modptrlow0, modptrhigh0;
  struct lowhigh *tmphigh;
  struct lowhigh *tmplowa, *tmplowb;

  lowptr = lowpoly;
  highptr = highpoly;
  
  levelsize = n;

  modptrlow = get_mods(levelsize);
  modptrhigh = modptrlow + (levelsize/2);
  
  if (permflag == 1)
    {
      OURpermuteX(data,lowpoly,n);
    }
  else
    {
      dummy_int = n / 2;
      for(i = 0; i < dummy_int ; i++)
	{
	  lowpoly[i].low = data[i];
	  lowpoly[i].high = data[i + dummy_int];
	}
    }
  
  /* do first level of interpolation  */
  dummy_int = n / 2;

  for ( i = 0 ; i < dummy_int ; i += 2)
    {
      (*highptr).low = (*lowptr).low + (*(lowptr+1)).low;
      (*highptr).high = (*lowptr).high + (*(lowptr+1)).high;
      highptr++;
    
      (*highptr).low = -(*modptrlow) *
	( (*lowptr).low - (*(lowptr+1)).low) ;
      (*highptr).high = -(*modptrhigh) *
	( (*lowptr).high - (*(lowptr+1)).high) ;
      highptr++;
      modptrlow += 2; modptrhigh += 2;
      lowptr += 2;

    }

  /* lowptr = highpoly; */
  lowptr = highpoly;
  highptr = lowpoly;
  loopcntr = 0;
  
  /* now do higher-order interpolations */
  if (k > 2) {
    
    /* reset pointers to logical polynomial workspace */
    levelsize /= 2;
    modptrlow = get_mods(levelsize);
    modptrhigh = modptrlow + (levelsize/2);
    highsize = 4;
    lowsize = 2;
    
    /* now set counter pointers to appropriate matrix locations */
    p2 = lowsize - 1;
    p3 = 1;
    
    /* main loop */
    while (highsize <= k) {
      dummy_int = n / 2;      

      for (j=0; j< dummy_int; j+=highsize) {
	/* do loworder coeffs first */

	tmplowa = lowptr;
	tmphigh = highptr;
	i = 0;
	do
	  {

	    tmphigh->low =
	      tmplowa->low + (tmplowa+lowsize)->low;
	    (tmphigh) ->high =
	      (tmplowa) ->high + (tmplowa+lowsize)->high;

	    tmplowa ++;
	    tmphigh ++;

	    tmphigh->low =
	      tmplowa->low + (tmplowa+lowsize)->low;
	    (tmphigh) ->high =
	      (tmplowa) ->high + (tmplowa+lowsize)->high;

	    tmplowa ++; tmphigh ++;
	    i += 2;

	  } while (i < lowsize );

	highptr += lowsize;
	
	/* now do higher order coeffs */
	
	modptrlow0 = *modptrlow;
	highptr->low = modptrlow0 * ((lowptr+lowsize)->low) -
	  modptrlow0 * (lowptr->low) ;

	modptrhigh0 = *modptrhigh;
	highptr->high = modptrhigh0 * ((lowptr+lowsize)->high) -
	  modptrhigh0 * (lowptr->high);

	modptrlow0 *= 2.0;
	modptrhigh0 *= 2.0;

	tmplowa = lowptr+p3;
	tmplowb = lowptr+p2;
	tmphigh = highptr+1;

	e0 = (tmplowa + lowsize)->low - tmplowa->low;
	e1 = (tmplowb + lowsize)->low + tmplowb->low;
	tmphigh->low = modptrlow0 * e0 - e1;
	
	f0 = (tmplowa + lowsize)->high - tmplowa->high;
	f1 = (tmplowb + lowsize)->high + tmplowb->high;
	tmphigh->high = modptrhigh0 * f0 - f1;
	
	tmplowa++; tmplowb--; tmphigh++;

	dummy_int2 = lowsize - 1;

	for( i = 1 ; i < dummy_int2 ; i += 2 )
	  {
	    e0 = (tmplowa + lowsize)->low - tmplowa->low;
	    e0x = (tmplowa + lowsize)->high - (tmplowa)->high;
	    e1 = (tmplowb + lowsize)->low + tmplowb->low;
	    e1x = (tmplowb + lowsize)->high + (tmplowb)->high;

	    tmphigh->low = modptrlow0 * e0 - e1;
	    (tmphigh) ->high = modptrhigh0 * e0x - e1x;

	    tmphigh++;
	    tmplowa++; tmplowb--;

	    f0 = (tmplowa + lowsize)->low - tmplowa->low;
	    f0x = (tmplowa + lowsize)->high - (tmplowa)->high;

	    f1 = (tmplowb + lowsize)->low + tmplowb->low;
	    f1x = (tmplowb + lowsize)->high + (tmplowb)->high;

	    tmphigh->low = modptrlow0 * f0 - f1;
	    (tmphigh)->high = modptrhigh0 * f0x - f1x;

	    tmplowa ++ ; tmplowb --; tmphigh ++;
	  }

	highptr += lowsize;
	lowptr += highsize;
	modptrlow += 2;
	modptrhigh += 2;
	p2 = lowsize - 1;
	p3 = 1;

      }
      
      levelsize /= 2;
      modptrlow = get_mods(levelsize);
      modptrhigh = modptrlow + (levelsize/2);
      highsize *= 2;
      lowsize *= 2;
      p2 = lowsize - 1;
      p3 = 1;
      
      if ((loopcntr % 2) == 0) {
	lowptr = lowpoly;
	highptr = highpoly;
      }
      else {
	lowptr = highpoly;
	highptr = lowpoly;
      }
      loopcntr++;
    }
  }
  
  /* find where the top level polynomials are */
  if (!(loopcntr % 2))
    temp = highpoly;
  else
    temp = lowpoly;
  
  if (k < n) 
    PartitionAddX(temp,result,n,k);
  else
    for (i=0; i<n; i++){
      result[i] = temp[i].low;
      result[i+(n/2)] = temp[i].high;
    }

  /* normalize  */
  dn2 = 2.0/((double) n);
  result[0] /= ((double) n);
  for(i=1; i<k; i++)
    result[i] *= dn2;

  /* later */
}


#endif   /* FFTPACK */


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

  Now for the inverse DCT routine, ExpIFCT, which is used regardless
  of FFTPACK being defined or not.

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

/***********************************************************************
  Chebyshev division program.  This is specifically coded for dividing 
  a chebyshev series of degree (n-1) by a modulus of the form
  T_n/2 + m(0)T_0, which returns a remainder sequence of degree
  (n/2) - 1.  Expects dividend coefficient sequence to be ordered from high
  degree to low degree.  Modulus should be a double value, m(0).  size is
  just n, the degree + 1 of the dividend.  The remainder is returned 
  with coefficients ordered from high degree to low degree.
  
  This function is used in the inverse DCT routine ExpIFCT.
  

  dividend - double array of size n
  remres - double array of size n/2
  m0 - supermodulus zero-degree coefficient value 
  
  *******************************************************************/

/****
  THE ORIGINAL ... easiest to understand how the algorithm actually
  works ... the code uses more efficient implementations of this
  original version.

  ****/

/*
static void ChebyDivRemOrig(double *dividend,
			    double m0,
			    double *remres,
			    int n)
{
  int rlim, dd;

  rlim = n/2;

  memcpy(remres, dividend + rlim, sizeof(double) * rlim );
  for (dd=0; dd<rlim-1; dd++) {
    remres[rlim-1-dd-1] -= dividend[dd];
    remres[dd] -= (2.0 * m0 * dividend[dd]);
  }
  remres[rlim-1] -= m0 * dividend[rlim-1];
}
*/

/***

  The version that is actually used in the code.

  ***/

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

static void ChebyDivRem(double *dividend,
			double m0,
			double *remres,
			int n)
{
  int rlim, dd;
  int tmpfloor;
  double tm0;
  int halfrlim;

  rlim = n/2;
  tm0 = 2.0 * m0;
  tmpfloor = (rlim - 1)/2;
  halfrlim = rlim/2;

  for(dd = 0; dd < tmpfloor; dd++)
    {
      remres[dd] =
	-tm0 * dividend[dd] - dividend[rlim-dd-2] + dividend[rlim+dd];
      remres[dd+halfrlim] =
	-dividend[halfrlim-2-dd] - tm0 * dividend[dd+halfrlim]
	+ dividend[rlim+dd+halfrlim];
    }
  remres[tmpfloor] = dividend[tmpfloor+rlim] -
    dividend[tmpfloor] - (tm0 * dividend[tmpfloor]);

  remres[rlim-1] = dividend[n-1] - (m0 * dividend[rlim-1]);

}

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

/*****************************************************************
  Fast Expanded Inverse Chebyshev transform

  Expects input data to be a Chebyshev series, with low-order 
  coefficients occuring first.
  Evaluates a Chebyshev series of length k at n data points,
  where n >= k.  Assumes power of 2

  data - double array of length k
  result - double array of length n
  workspace - double array of length (2*n)
  n - number of samples (evaluation points) desired
  k - number of coefficients 
  permflag = if 0, then don't permute the result.  If 1,
             permute data using OUR permutation

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

void ExpIFCT( double *data,
	      double *result,
	      double *workspace,
	      int n,
	      int k,
	      int permflag)
{
  int j, loopcntr;
  int levelsize, divsize;
  double *dividend, *remres;
  const double *modptr;
  double *divptr, *remptr;
  double divptr0, divptr1, divptr2, divptr3;
  double tmpmod, tmpmod2;

  /* assign workspace locations */
  remres = workspace;
  dividend = workspace + n;
  
  /* first reverse the Chebyshev series to order from high to low */
  reverse_it(data, dividend, k);
  
  /* compute initial values of variables */
  /* levelsize is the number of supermoduli at first level of division */
  /* divsize is the size of each divisor */
  /* modptr points to list of supermods at this level */
  
  levelsize = 2 * (n/k);
  divsize = k;
  modptr = get_mods(levelsize);
  divptr = dividend;
  remptr = remres;
  loopcntr = 0;
  
  /* now do divisions */
  for(j = 0; j < levelsize; j += 2)
    {
      ChebyDivRem(divptr,modptr[j],remptr,divsize);
      remptr += (divsize>>1);
      ChebyDivRem(divptr,modptr[j+1],remptr,divsize);
      remptr += (divsize>>1);
    }
  
  /* now reset everything */
  divptr = remres;
  remptr = dividend;
 
  loopcntr++;
  divsize /= 2;
  levelsize *= 2;
  modptr = get_mods(levelsize);

  /*** now loopcntr = 1 ***/
  while (divsize > 4)
    {

      for(j = 0 ; j < levelsize; j += 2)
	{
	  ChebyDivRem(divptr,modptr[j],remptr,divsize);
	  remptr += (divsize>>1);
	  ChebyDivRem(divptr,modptr[j+1],remptr,divsize);
	  remptr += (divsize>>1);
	  divptr += divsize;
	}
      
      /* now reset everything */
      if ( (loopcntr % 2) == 0 )
	{
	  divptr = remres;
	  remptr = dividend;
	}
      else 
	{
	  divptr = dividend;
	  remptr = remres;
	}
      
      loopcntr++;
      divsize /= 2;
      levelsize *= 2;
      modptr = get_mods(levelsize);
    }


  /***** in what follows, am using the fact
    that modptr[i+1] = -modptr[i] for i even ****/
  
  if(divsize == 4)
    {
      for (j=0; j<levelsize; j+=2)
	{

	  divptr0 = divptr[0]; divptr1 = divptr[1];
	  divptr2 = divptr[2]; divptr3 = divptr[3];
	  tmpmod = modptr[j];
	  tmpmod2 = 2.0 * modptr[j];
	  remptr[0] = divptr2 - divptr0 - (tmpmod2 * divptr0);
	  remptr[1] = (divptr3 - (tmpmod * divptr1));
	  remptr[2] = divptr2 - divptr0 + (tmpmod2 * divptr0);
	  remptr[3] = (divptr3 + (tmpmod * divptr1));
	  
	  remptr += 4;
	  divptr += 4; 
	}
      
      /* now reset everything */
      if ( (loopcntr % 2) == 0)
	{
	  divptr = remres;
	  remptr = dividend;
	}
      else
	{
	  divptr = dividend;
	  remptr = remres;
	}      
      loopcntr++;
      divsize /= 2;
      levelsize *= 2;
      modptr = get_mods(levelsize);
    }
  
  if(divsize == 2)
    {
      for( j = 0; j < levelsize; j += 2)
	{
	  divptr0 = divptr[0]; divptr1 = divptr[1];
	  
	  remptr[0] = divptr1 - (modptr[j] * divptr0);
	  remptr[1] = divptr1 + (modptr[j] * divptr0);
	  
	  remptr += 2;
	  divptr += divsize;
	}
      
      /* now reset everything */
      if ( (loopcntr % 2) == 0 )
	{
	  divptr = remres;
	  remptr = dividend;
	}
      else
	{
	  divptr = dividend;
	  remptr = remres;
	}
    }
  
  /******** *********/
  /* divptr now points to result - copy into result space 
     after checking permute flag */
  if ( permflag == 0 ) {
    memcpy(result, divptr, sizeof(double) * n);
  }
  else 
    OURpermute(divptr,result,n);
  
  /* later */
  
}

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