cospmls.cpp

普林斯顿开发的快速球面调和变换算法

// cospmls.cpp: implementation of the cospmls class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "cospmls.h"
#include <math.h>
#include <string.h>   /* to declare memcpy */
#include <malloc.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

cospmls::cospmls()
{

}

cospmls::~cospmls()
{

}

int cospmls::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 ;
   }

}

int cospmls::TableSize(int m, int bw)
{
	return (((bw/2) * ((bw/2) + 1)) -
	    ((m/2)*((m/2)+1)) - ((bw/2) * (m % 2)));

}

int cospmls::Spharmonic_TableSize(int bw)
{
	  int m, sum;
  
  if (bw > 512)
    {
      sum = 0;
      
      for (m=0; m<bw; m++)
	  sum += this->TableSize(m,bw);
      
      return sum;
    }
  else   //适合我们的设置情况
    {
      return (
	      (((4*(bw*bw*bw)) + (6*(bw*bw)) - (8*bw))/24)
	      + bw
	      );
    }

}

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

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

  return sum;

}

int cospmls::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;

}

int cospmls::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;

}

void cospmls::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;
    double *CoswSave;

    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 = newf.precomp_dct( bw );

#endif
    
    /* main loop */

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

    /* now get the evaluation nodes */
    prim.EvalPts(bw,x_i);
    prim.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 
      prim.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
    newf.kFCT(prev, cosres, cosworkspace, bw, bw, 1);
#else
    memcpy(cosres, prev, sizeof(double) * bw);
    newf.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++) {
	prim.vec_mul(prim.L2_cn(m,m+i),prevprev,temp1,bw);
	prim.vec_pt_mul(prev, x_i, temp2, bw);
	prim.vec_mul(prim.L2_an(m,m+i), temp2, temp3, bw);
	prim.vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */
	/* compute cosine transform */
#ifndef FFTPACK	
	newf.kFCT(temp4, cosres, cosworkspace, bw, bw, 1);
#else
	memcpy(cosres, temp4, sizeof(double) * bw);
	newf.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

}

void cospmls::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;
    double *CoswSave;

    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 = newf.precomp_dct( bw );

#endif

    /* main loop */

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

    /* now get the evaluation nodes */
    prim.EvalPts(bw,x_i);
    prim.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 
      prim.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
    newf.kFCT(prev, cosres, cosworkspace, bw, bw, 1);
#else
    memcpy(cosres, prev, sizeof(double) * bw);
    newf.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 ++) {
	prim.vec_mul(prim.L2_cn(m,m+i),prevprev,temp1,bw);
	prim.vec_pt_mul(prev, x_i, temp2, bw);
	prim.vec_mul(prim.L2_an(m,m+i), temp2, temp3, bw);
	prim.vec_add(temp3, temp1, temp4, bw); /* temp4 now contains P(m,m+i+1) */
	/* compute cosine transform */
#ifndef FFTPACK	
	newf.kFCT(temp4, cosres, cosworkspace, bw, bw, 1);
#else
	memcpy(cosres, temp4, sizeof(double) * bw);
	newf.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


}

int cospmls::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);
    }


}

int cospmls::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));
    }

}

void cospmls::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.
  */
  
  double *trans_tableptr, *tableptr;
  int i, row, rowsize, stride, offset, costable_offset;