cospmls.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.
  
  ************************************************************************
  ************************************************************************/


/* source code for generating cosine transforms 
   of Pml and Gml functions */

#include <math.h>
#include <string.h>   /* to declare memcpy */
#include <stdlib.h>
#include "newFCT.h"
#include "primitive.h"


#ifdef FFTPACK
#include "fftpack.h"
#endif


#ifndef PI
#define PI 3.14159265358979
#endif

/************************************************************************/
/* look up table for the power of 2 which is >= i */

int Power2Ceiling(int i)
{
  int pow_2;
  
  pow_2 = 1;

  if (i == 0)
    return 0;
  else
    {
      while ( pow_2 < i )
	pow_2 *= 2;
      return pow_2 ;
    }
}

/************************************************************************/
/* utility functions for table management */
/************************************************************************/
/* Computes the number of non-zero entries in a table containing
   cosine series coefficients of all the P(m,l) or G(m,l) functions
   necessary for computing the seminaive transform for a given
   bandwidth bw and order m.  Works specifically for tables
   generated by CosPmlTableGen() 
*/

int TableSize(int m,
	      int bw)
{
    return (((bw/2) * ((bw/2) + 1)) -
	    ((m/2)*((m/2)+1)) - ((bw/2) * (m % 2)));
}
/************************************************************************/
/* Spharmonic_TableSize(bw) returns an integer value for
   the amount of space necessary to fill out an entire spharmonic
   table.  Note that in the above TableSize() formula, 
   you need to sum this formula over m as m ranges from 0 to
   (bw-1).  The critical closed form that you need is that

   \sum_{k=0}^n = \frac{(n(n+1)(2n+1)}{6}

   You also need to account for integer division.
   From this you should derive an upper bound on the
   amount of space. 

   Some notes - because of integer division, you need to account for
   a fudge factor - this is the additional " + bw" at the
   end.  This gaurantees that you will always have slightly more
   space than you need, which is clearly better than underestimating!
   Also, if bw > 512, the closed form
   fails because of the bw*bw*bw term (at least on Sun Sparcstations)
   so the loop computation is used instead.

   Also, the transpose is exactly the same size, obviously.

*/

int Spharmonic_TableSize(int bw)
{
  int m, sum;
  
  if (bw > 512)
    {
      sum = 0;
      
      for (m=0; m<bw; m++)
	  sum += TableSize(m,bw);
      
      return sum;
    }
  else
    {
      return (
	      (((4*(bw*bw*bw)) + (6*(bw*bw)) - (8*bw))/24)
	      + bw
	      );
    }
}

/************************************************************************/
/* Reduced_Spharmonic_TableSize(bw,m) returns an integer value for
   the amount of space necessary to fill out a spharmonic table
   if interesting in using it only for orders up to (but NOT
   including) order m.
   This will be used in the hybrid algorithm's call of the
   semi-naive algorithm (which won't need the full table ... hopefully
   this'll cut down on the memory usage).

   Also, the transpose is exactly the same size, obviously.

   This is a "reduced" version of Spharmonic_TableSize(m).

*/

int Reduced_SpharmonicTableSize(int bw,
				int m)
{
  
  int i, sum;

  sum = 0;
  
  for (i=0; i<m; i++)
      sum += TableSize(i,bw);

  return sum;
}


/************************************************************************/
/* For an array containing cosine series coefficients of Pml or Gml
   functions, computes the location of the first coefficient of Pml.
   This supersedes the TableOffset() function.
   Assumes table is generated by CosPmlTableGen()
*/

int NewTableOffset(int m,
		   int l)
{
    int offset;
    int tm, tl;
    
    if ((m % 2) == 1)
      {
	tl = l-1;
	tm = m-1;
      }
    else
      {
	tl = l;
	tm = m;
      }

    offset = ((tl/2)*((tl/2)+1)) - ((tm/2)*((tm/2)+1));
    if (tl % 2 == 1)
      offset += (tl/2)+1;

    return offset;
}
/************************************************************************/
/* Computes the offset in the table generated by CosPmlTableGen for
   the cosine series of Pml OR GML.  NOTE THAT M MUST BE EVEN - 
   MUST BE HANDLED BY CALLER 

   Superseded by NewTableOffset()
*/

int TableOffset(int m,
		int l)
{
    int offset;

    offset = ((l/2)*((l/2)+1)) - ((m/2)*((m/2)+1));
    if (l % 2 == 1)
      offset += (l/2)+1;

    return offset;
}
/************************************************************************/
/* generate all of the cosine series for Pmls or Gmls for a specified 
   value of m.  Note especially that since series are zero-striped,
   all zeroes have been removed.  

   tablespace points to a double array of size TableSize(m,bw);

   Workspace needs to be
   16 * bw 

   Let P(m,l,j) represent the jth coefficient of the
   cosine series representation of Pml.  The array
   stuffed into tablespace is organized as follows:

   P(m,m,0)    P(m,m,2)  ...  P(m,m,m)
   P(m,m+1,1)  P(m,m+1,3)...  P(m,m+1,m+1)
   P(m,m+2,0)  P(m,m+2,2) ... P(m,m+2,m+2)

   etc.  Appropriate modifications are made for m odd (Gml functions).

*/

void CosPmlTableGen(int bw,
		    int m,
		    double *tablespace,
		    double *workspace)
{
    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;
    double *tableptr, *cosres, *cosworkspace;
    int i, j, k;

    prevprev = workspace;
    prev = prevprev + bw;
    temp1 = prev + bw;
    temp2 = temp1 + bw;
    temp3 = temp2 + bw;
    temp4 = temp3 + bw;
    x_i = temp4 + bw;
    eval_args = x_i + bw;
    cosres = eval_args + bw;
    cosworkspace = cosres + bw;

    tableptr = tablespace;

#ifdef FFTPACK

    CoswSave = precomp_dct( bw );

#endif
    
    /* main loop */

    /* Set the initial number of evaluation points to appropriate
       amount */

    /* now get the evaluation nodes */
    EvalPts(bw,x_i);
    ArcCosEvalPts(bw,eval_args);
    
    /* set initial values of first two Pmls */
    for (i=0; i<bw; i++) 
      prevprev[i] = 0.0;
    if (m == 0) {
	for (i=0; i<bw; i++) {
	  prev[i] = 1.0 ;
	}
    }
    else 
      Pmm_L2(m, eval_args, bw, prev);

    
    if ((m % 2) == 1) { /* need to divide out sin x */
	for (i=0; i<bw; i++)
	  prev[i] /= sin(eval_args[i]);
    }
  
    /* set k to highest degree coefficient */
    if ((m % 2) == 0)
      k = m;
    else
      k = m-1;	

    /* now compute cosine transform */
#ifndef FFTPACK
    kFCT(prev, cosres, cosworkspace, bw, bw, 1);
#else
    memcpy(cosres, prev, sizeof(double) * bw);
    DCTf( cosres, bw, bw, CoswSave );
#endif

    for (i=0; i<=k; i+=2)
      tableptr[i/2] = cosres[i];

    /* update tableptr */
    tableptr += k/2+1;

    /* now generate remaining pmls  */

    for (i=0; i<bw-m-1; i++) {
	vec_mul(L2_cn(m,m+i),prevprev,temp1,bw);
	vec_pt_mul(prev, x_i, temp2, bw);
	vec_mul(L2_an(m,m+i), temp2, temp3, bw);
	vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */
	/* compute cosine transform */
#ifndef FFTPACK	
	kFCT(temp4, cosres, cosworkspace, bw, bw, 1);
#else
	memcpy(cosres, temp4, sizeof(double) * bw);
	DCTf( cosres, bw, bw, CoswSave );
#endif

	/* update degree counter */
	k++;

	/* now put decimated result into table */
	if ((i % 2) == 1) {
	    for (j=0; j<=k; j+=2)
	      tableptr[j/2] = cosres[j];
	}
	else {
	    for (j=1; j<=k; j+=2)
	      tableptr[j/2] = cosres[j];
	}
	/* update tableptr */
	tableptr += k/2+1;


	/* now update Pi and P(i+1) */
	memcpy(prevprev, prev, sizeof(double) * bw);
	memcpy(prev, temp4, sizeof(double) * bw);
    }

#ifdef FFTPACK
    free( CoswSave );
#endif


}

/****

  just like above except does not compute all the
  legendres for a given m, bw that semi-naive would
  use, just those that semi-naive IN THE FLT_HYBRID CODE
  would use

  lim = degree of FIRST coefficient in the flt_hybrid
        code that will NOT be computed using semi-naive

  lim is defined in this strange way so that this
  function will produce the SAME array as
  CosPmlTableGen if lim = bw.

  tablespace points to a double array of size cos_size
  as computed in SplitLocales().

  *****/


void CosPmlTableGenLim( int bw, 
			int m,
			int lim,
			double *tablespace,
			double *workspace)
{
    double *prev, *prevprev, *temp1, *temp2, *temp3, *temp4, *x_i, *eval_args;
    double *tableptr, *cosres, *cosworkspace;
    int i, j, k;

    prevprev = workspace;
    prev = prevprev + bw;
    temp1 = prev + bw;
    temp2 = temp1 + bw;
    temp3 = temp2 + bw;
    temp4 = temp3 + bw;
    x_i = temp4 + bw;
    eval_args = x_i + bw;
    cosres = eval_args + bw;
    cosworkspace = cosres + bw;

    tableptr = tablespace;

#ifdef FFTPACK

    CoswSave = precomp_dct( bw );

#endif

    /* main loop */

    /* Set the initial number of evaluation points to appropriate
       amount */

    /* now get the evaluation nodes */
    EvalPts(bw,x_i);
    ArcCosEvalPts(bw,eval_args);
    
    /* set initial values of first two Pmls */
    for (i=0; i<bw; i++) 
      prevprev[i] = 0.0;
    if (m == 0) {
	for (i=0; i<bw; i++) {
	    prev[i] = 1.0;
	}
    }
    else 
      Pmm_L2(m, eval_args, bw, prev);

    
    if ((m % 2) == 1) { /* need to divide out sin x */
	for (i=0; i<bw; i++)
	  prev[i] /= sin(eval_args[i]);
    }
  
    /* set k to highest degree coefficient */
    if ((m % 2) == 0)
      k = m;
    else
      k = m-1;	

    /* now compute cosine transform */
#ifndef FFTPACK
    kFCT(prev, cosres, cosworkspace, bw, bw, 1);
#else
    memcpy(cosres, prev, sizeof(double) * bw);
    DCTf( cosres, bw, bw, CoswSave );
#endif

    for (i=0; i<=k; i+=2)
      tableptr[i/2] = cosres[i];

    /* update tableptr */
    tableptr += k/2+1;

    /* now generate remaining pmls  */
    for (i = 0 ; i < lim - m - 1 ; i ++) {
	vec_mul(L2_cn(m,m+i),prevprev,temp1,bw);
	vec_pt_mul(prev, x_i, temp2, bw);
	vec_mul(L2_an(m,m+i), temp2, temp3, bw);
	vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */
	/* compute cosine transform */
#ifndef FFTPACK	
	kFCT(temp4, cosres, cosworkspace, bw, bw, 1);
#else
	memcpy(cosres, temp4, sizeof(double) * bw);
	DCTf( cosres, bw, bw, CoswSave );
#endif

	/* update degree counter */
	k++;

	/* now put decimated result into table */
	if ((i % 2) == 1) {
	    for (j=0; j<=k; j+=2)
	      tableptr[j/2] = cosres[j];
	}
	else {
	    for (j=1; j<=k; j+=2)
	      tableptr[j/2] = cosres[j];
	}
	/* update tableptr */
	tableptr += k/2+1;

	/* now update Pi and P(i+1) */
	memcpy(prevprev, prev, sizeof(double) * bw);
	memcpy(prev, temp4, sizeof(double) * bw);

    }

#ifdef FFTPACK
    free( CoswSave );
#endif

}



/************************************************************************/
/* RowSize returns the number of non-zero coefficients in a row of the
   cospmltable if were really in matrix form.  Helpful in transpose
   computations.  It is helpful to think of the parameter l as
   the row of the corresponding matrix.
*/

int RowSize(int m,
	    int l)
{
  if (l < m)
    return 0;
  else
    {
      if ((m % 2) == 0)
	return ((l/2)+1);
      else
	return (((l-1)/2)+1);
    }
}
/************************************************************************/
/* Transposed row size returns the number of non-zero coefficients
   in the transposition of the matrix representing a cospmltable.
   Used for generating arrays for inverse seminaive transform.
   Unlike RowSize, need to know the bandwidth bw.  Also, in
   the cospml array, the first m+1 rows are empty, but in
   the transpose, all rows have non-zero entries, and the first
   m+1 columns are empty.  So the input parameters are a bit different
   in the you need to specify the row you want.

*/

int Transpose_RowSize(int row,
		      int m,
		      int bw)
{

  if (row >= bw)
    return 0;
  else if ((m % 2) == 0)
    {
      if (row <= m)
	return ( ((bw-m)/2) );
      else
	return ( ((bw-row-1)/2) + 1);
    }
  else
    {
      if (row == (bw-1))
	return 0;
      else if (row >= m)
	return (Transpose_RowSize(row+1,m-1,bw));
      else    /*  (row < m)  */
	return (Transpose_RowSize(row+1,m-1,bw) - (row % 2));
    }
}

/************************************************************************/
/* Inverse transform is transposition of forward transform.
   Thus, need to provide transposed version of table
   returned by CosPmlTableGen.  This function does that
   by taking as input a cos_pml_table for a particular value
   of bw and m, and loads the result as a
   transposed, decimated version of it for use by an inverse 
   seminaive transform computation.

   result needs to be of size TableSize(m,bw)

*/

void Transpose_CosPmlTableGen(int bw,
			      int m,
			      double *cos_pml_table,
			      double *result)
{
  /* recall that cospml_table has had all the zeroes
     stripped out, and that if m is odd, then it is
     really a Gml function, which affects indexing a bit.