spherical.cpp

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

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

#include "stdafx.h"
#include "Spherical.h"
#include <stdlib.h>
#include <stdlib.h>
#include <time.h>
#include "FST_semi_memo.h"
#include "MathFace.h"
#include "cospmls.h"
#include "FFT.h"
#include "newFCT.h"
#include "primitive_FST.h"
#define max(A, B) ((A) > (B) ? (A) : (B))
#define M_PI 3.14159265358979323846

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

Spherical::Spherical(int a, int b)
{

	R=a;
	B=b;
   
//	rsample = reinterpret_cast<double **>(new double[R * sizeof(double *)]); 
	rsample=(double **)malloc(R*sizeof(double *));
	isample=(double **)malloc(R*sizeof(double *));
	rcoeff=(double **)malloc(R*sizeof(double *));
	icoeff=(double **)malloc(R*sizeof(double *));
	feature=(double **)malloc(R*sizeof(double *));
    //为较高的维分配空间（也可以理解为二维数组的“行”）
    for(int i=0;i<R;i++)
	{
        rsample[i]=(double *)malloc(4*B*B*sizeof(double));
        isample[i]=(double *)malloc(4*B*B*sizeof(double));
		rcoeff[i]=(double *)malloc(B*B*sizeof(double));
		icoeff[i]=(double *)malloc(B*B*sizeof(double));
		feature[i]=(double *)malloc(B*sizeof(double));
	}
	//用循环为较低的维分配空间（也可以理解为二维数组的“列”）
    int j=0;
	int k;
	for(i=0;i<R;i++)
	{
		for(j=0;j<4*B*B;j++)
		{
			rsample[i][j]=0;
			isample[i][j]=0;
		}
		for(k=0;k<B*B;k++)
		{
			rcoeff[i][k]=0;
			icoeff[i][k]=0;
		}
	}



}

Spherical::~Spherical()
{
	
	free(rsample);
	free(isample);
	free(rcoeff);
	free(icoeff);
	free(feature);

}

void Spherical::Read(double **vox)
{

	int i,j;
	for(i=0;i<R;i++)
	{
		for(j=0;j<4*B*B;j++)
			rsample[i][j]=vox[i][j];
	}

	

}

void Spherical::CountA()
{
	srand( (unsigned)time( NULL ) );

	primitive_FST prif1;
	cospmls cosp1;
	FST_semi_memo fst1;
	newFCT newf1;

	int bw, size;
	register int i, j;
	double *seminaive_naive_tablespace, *trans_seminaive_naive_tablespace,
         *workspace;
	double **seminaive_naive_table, **trans_seminaive_naive_table;
	int l, m, dummy;
	double dumx, dumy, fudge;
	int cutoff;

	double realtmp, imagtmp,origmag, tmpmag;

	double *wSave, *CoswSave;
  /* for fftpack */
  bw=B;

    cutoff = bw ;

  size = 2*bw;   

#ifdef FFTPACK
  fprintf(stdout,"precomputing for fft, dct\n");

  wSave = precomp_fft( size );

  CoswSave = newf1.precomp_dct( size );

#endif

  seminaive_naive_tablespace =
    (double *) malloc(sizeof(double) *		      (cosp1.Reduced_Naive_TableSize(bw,cutoff) +
		       cosp1.Reduced_SpharmonicTableSize(bw,cutoff)));

  trans_seminaive_naive_tablespace =
    (double *) malloc(sizeof(double) *
		      (cosp1.Reduced_Naive_TableSize(bw,cutoff) +
		       cosp1.Reduced_SpharmonicTableSize(bw,cutoff)));
  
#ifndef FFTPACK  
  workspace = (double *) malloc(sizeof(double) * 
				((8 * (bw*bw)) + 
				 (16 * bw)));

#else
  workspace = (double *) malloc(sizeof(double) * 
				((8 * (bw*bw)) + 
				 (16 * bw)));
#endif
  
int k=0;
for(k=0;k<R;k++)
{
  #ifndef FFTPACK
  if ( (rsample[k] == NULL) || (isample[k] == NULL) ||
       (rcoeff[k] == NULL) || (icoeff[k] == NULL) ||
       (seminaive_naive_tablespace == NULL) ||
       (trans_seminaive_naive_tablespace == NULL) ||
       (workspace == NULL) )
    {
      perror("Error in allocating memory");
      exit( 1 ) ;
    }
#else
  if ( (rsample[k] == NULL) ||
       (rcoeff[k] == NULL) || (icoeff[k] == NULL) ||
       (seminaive_naive_tablespace == NULL) ||
       (trans_seminaive_naive_tablespace == NULL) ||
       (workspace == NULL) )
    {
      perror("Error in allocating memory");
      exit( 1 ) ;
    }
#endif

  fprintf(stdout,"Generating seminaive_naive tables...\n");
  seminaive_naive_table = cosp1.SemiNaive_Naive_Pml_Table(bw, cutoff,
						    seminaive_naive_tablespace,
						    workspace);

  #ifndef FFTPACK
  	
    fst1.FST_semi_memo1(rsample[k], isample[k],
		  rcoeff[k], icoeff[k],
		  size,
		  seminaive_naive_table,
		  workspace,
		  1,
		  cutoff);
	
#else
    fst1.FST_semi_memo2(rsample[k],
		  rcoeff[k], icoeff[k],
		  size,
		  seminaive_naive_table,
		  workspace,
		  wSave,
		  CoswSave,
		  cutoff);
#endif
	
}


  free(workspace);
  free(trans_seminaive_naive_tablespace);
  free(seminaive_naive_tablespace);
#ifdef FFTPACK
  free(CoswSave);
  free(wSave);
#endif
 
 CountFea();	



}


void Spherical::CountFea()
{
	int i,j,l;
	double mf;
	int dummy;
	primitive_FST prif2;
	for(i=0;i<R;i++)
	{
		
		for(j=0;j<B;j++)
		{
			mf=0;
			for(l = -j ; l<= j ; l++ )
			{
				dummy = prif2.seanindex(l,j,B);	
				mf=mf+(rcoeff[i][dummy])*(rcoeff[i][dummy])+(icoeff[i][dummy])*(icoeff[i][dummy]);
			}
			feature[i][j]=sqrt(mf);

		}
	}
}


double* Spherical::precomp_fft(int size)
{
	double *wSave;
  /* allocate space *///分配空间
  wSave = (double *) malloc( sizeof(double) * (2 * size + 15) );

  /* now precompute data */
  FFT fft2; 
  fft2.rffti_( &size, wSave );
  return wSave;

}