Commit 82c718f2 authored by David Rohr's avatar David Rohr
Browse files

use improved FFTs for real and hermitian input

parent da3cdbbc
...@@ -240,9 +240,8 @@ int bioem::run() ...@@ -240,9 +240,8 @@ int bioem::run()
myfloat_t sumCONV,sumsquareCONV; myfloat_t sumCONV,sumsquareCONV;
//allocating fftw_complex vector //allocating fftw_complex vector
proj_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); proj_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
conv_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t)*param.param_device.NumberPixels*param.param_device.NumberPixels); conv_mapFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
HighResTimer timer; HighResTimer timer;
...@@ -378,9 +377,10 @@ int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myf ...@@ -378,9 +377,10 @@ int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myf
{ {
#pragma omp parallel #pragma omp parallel
{ {
mycomplex_t *localCCT, *lCC; mycomplex_t *localCCT;
localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); myfloat_t *lCC;
lCC= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); localCCT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
lCC= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
const int num_threads = omp_get_num_threads(); const int num_threads = omp_get_num_threads();
const int thread_id = omp_get_thread_num(); const int thread_id = omp_get_thread_num();
...@@ -400,40 +400,40 @@ int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myf ...@@ -400,40 +400,40 @@ int bioem::compareRefMaps2(int iOrient, int iConv, mycomplex_t* localConvFFT,myf
} }
/////////////NEW ROUTINE //////////////// /////////////NEW ROUTINE ////////////////
inline int bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC) inline int bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC)
{ {
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
{ {
for(int j=0; j < param.param_device.NumberPixels ; j++ ) for(int j=0; j < param.param_device.NumberFFTPixels1D ; j++ )
{ {
localCCT[i*param.param_device.NumberPixels+j][0]=localConvFFT[i*param.param_device.NumberPixels+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]+localConvFFT[i*param.param_device.NumberPixels+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1]; localCCT[i*param.param_device.NumberFFTPixels1D+j][0]=localConvFFT[i*param.param_device.NumberFFTPixels1D+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]+localConvFFT[i*param.param_device.NumberFFTPixels1D+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
localCCT[i*param.param_device.NumberPixels+j][1]=localConvFFT[i*param.param_device.NumberPixels+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]-localConvFFT[i*param.param_device.NumberPixels+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1]; localCCT[i*param.param_device.NumberFFTPixels1D+j][1]=localConvFFT[i*param.param_device.NumberFFTPixels1D+j][1]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]-localConvFFT[i*param.param_device.NumberFFTPixels1D+j][0]*RefMap.RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
} }
} }
myfftw_execute_dft(param.fft_plan_c2c_backward,localCCT,lCC); myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,localCCT,lCC);
// Storing CORRELATIONS FOR CORRESPONDING DISPLACEMENTS & Normalizing after Backward FFT // Storing CORRELATIONS FOR CORRESPONDING DISPLACEMENTS & Normalizing after Backward FFT
for (int cent_x = 0; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter) for (int cent_x = 0; cent_x <= param.param_device.maxDisplaceCenter; cent_x=cent_x+param.param_device.GridSpaceCenter)
{ {
for (int cent_y = 0; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = 0; cent_y <= param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y][0]/ (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels), cent_x, cent_y); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y]/ (myfloat_t) (param.param_device.NumberPixels * param.param_device.NumberPixels), cent_x, cent_y);
} }
for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y < param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y < param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y][0]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), cent_x, param.param_device.NumberPixels-cent_y); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), cent_x, param.param_device.NumberPixels-cent_y);
} }
} }
for (int cent_x = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_x < param.param_device.NumberPixels; cent_x=cent_x+param.param_device.GridSpaceCenter) for (int cent_x = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_x < param.param_device.NumberPixels; cent_x=cent_x+param.param_device.GridSpaceCenter)
{ {
for (int cent_y = 0; cent_y < param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = 0; cent_y < param.param_device.maxDisplaceCenter; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y][0]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, cent_y); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, cent_y);
} }
for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y <= param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter) for (int cent_y = param.param_device.NumberPixels-param.param_device.maxDisplaceCenter; cent_y <= param.param_device.NumberPixels; cent_y=cent_y+param.param_device.GridSpaceCenter)
{ {
calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y][0]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, param.param_device.NumberPixels-cent_y); calProb(iRefMap, iOrient, iConv, sumC, sumsquareC, (myfloat_t) lCC[cent_x*param.param_device.NumberPixels+cent_y]/ (myfloat_t) (param.param_device.NumberPixels*param.param_device.NumberPixels), param.param_device.NumberPixels-cent_x, param.param_device.NumberPixels-cent_y);
} }
} }
...@@ -499,9 +499,9 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT) ...@@ -499,9 +499,9 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
myfloat3_t RotatedPointsModel[Model.nPointsModel]; myfloat3_t RotatedPointsModel[Model.nPointsModel];
myfloat_t rotmat[3][3]; myfloat_t rotmat[3][3];
myfloat_t alpha, gam,beta; myfloat_t alpha, gam,beta;
mycomplex_t* localproj; myfloat_t* localproj;
localproj= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); localproj= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
memset(localproj,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj)); memset(localproj,0,param.param_device.NumberPixels*param.param_device.NumberPixels*sizeof(*localproj));
alpha=param.angles[iMap].pos[0]; alpha=param.angles[iMap].pos[0];
...@@ -549,7 +549,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT) ...@@ -549,7 +549,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
i=floor(RotatedPointsModel[n].pos[0]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f); i=floor(RotatedPointsModel[n].pos[0]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f); j=floor(RotatedPointsModel[n].pos[1]/param.pixelSize+ (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
localproj[i*param.param_device.NumberPixels+j][0]+=Model.densityPointsModel[n]/Model.NormDen; localproj[i*param.param_device.NumberPixels+j]+=Model.densityPointsModel[n]/Model.NormDen;
} }
/**** Output Just to check****/ /**** Output Just to check****/
...@@ -562,7 +562,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT) ...@@ -562,7 +562,7 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n"; myexamplemap << "ANGLES " << alpha << " " << beta << " " << gam << "\n";
for(int k=0; k<param.param_device.NumberPixels; k++) for(int k=0; k<param.param_device.NumberPixels; k++)
{ {
for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j][0]; for(int j=0; j<param.param_device.NumberPixels; j++) myexamplemap << "\nMAP " << k << " " << j<< " " <<localproj[k*param.param_device.NumberPixels+j];
} }
myexamplemap << " \n"; myexamplemap << " \n";
for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2]; for(int n=0; n< Model.nPointsModel; n++)myexampleRot << "\nCOOR " << RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2];
...@@ -572,12 +572,12 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT) ...@@ -572,12 +572,12 @@ int bioem::createProjection(int iMap,mycomplex_t* mapFFT)
/***** Converting projection to Fourier Space for Convolution later with kernel****/ /***** Converting projection to Fourier Space for Convolution later with kernel****/
/********** Omp Critical is necessary with FFTW*******/ /********** Omp Critical is necessary with FFTW*******/
myfftw_execute_dft(param.fft_plan_c2c_forward,localproj,mapFFT); myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localproj,mapFFT);
return(0); return(0);
} }
int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv,mycomplex_t* localmultFFT,myfloat_t& sumC,myfloat_t& sumsquareC) int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,bioem_map& Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC)
{ {
/**************************************************************************************/ /**************************************************************************************/
/**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier ********** /**** BioEM Create Convoluted Projection Map routine, multiplies in Fourier **********
...@@ -585,31 +585,35 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b ...@@ -585,31 +585,35 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
*************** and Backtransforming it to real Space ******************************/ *************** and Backtransforming it to real Space ******************************/
/**************************************************************************************/ /**************************************************************************************/
mycomplex_t* localconvFFT; myfloat_t* localconvFFT;
localconvFFT= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t)*param.param_device.NumberPixels*param.param_device.NumberPixels); localconvFFT= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t)*param.param_device.NumberPixels*param.param_device.NumberPixels);
mycomplex_t* tmp;
tmp = (mycomplex_t*) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
/**** Multiplying FFTmap with corresponding kernel ****/ /**** Multiplying FFTmap with corresponding kernel ****/
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
{ {
for(int j=0; j < param.param_device.NumberPixels ; j++ ) for(int j=0; j < param.param_device.NumberFFTPixels1D ; j++ )
{ //Projection*CONJ(KERNEL) { //Projection*CONJ(KERNEL)
localmultFFT[i*param.param_device.NumberPixels+j][0]=lproj[i*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0]+lproj[i*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1]; localmultFFT[i*param.param_device.NumberFFTPixels1D+j][0]=lproj[i*param.param_device.NumberFFTPixels1D+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]+lproj[i*param.param_device.NumberFFTPixels1D+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
localmultFFT[i*param.param_device.NumberPixels+j][1]=lproj[i*param.param_device.NumberPixels+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0]-lproj[i*param.param_device.NumberPixels+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1]; localmultFFT[i*param.param_device.NumberFFTPixels1D+j][1]=lproj[i*param.param_device.NumberFFTPixels1D+j][1]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]-lproj[i*param.param_device.NumberFFTPixels1D+j][0]*param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1];
// cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][0] << " " <<param.refCTF[iConv].cpoints[i*param.param_device.NumberPixels+j][1] <<" " <<lproj[i*param.param_device.NumberPixels+j][0] <<" " <<lproj[i*param.param_device.NumberPixels+j][1] << "\n"; // cout << "GG " << i << " " << j << " " << param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][0] << " " <<param.refCTF[iConv].cpoints[i*param.param_device.NumberFFTPixels1D+j][1] <<" " <<lproj[i*param.param_device.NumberFFTPixels1D+j][0] <<" " <<lproj[i*param.param_device.NumberFFTPixels1D+j][1] << "\n";
} }
} }
//FFTW_C2R will destroy the input array, so we have to work on a copy here
memcpy(tmp, localmultFFT, sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
/**** Bringing convoluted Map to real Space ****/ /**** Bringing convoluted Map to real Space ****/
myfftw_execute_dft(param.fft_plan_c2c_backward,localmultFFT,localconvFFT); myfftw_execute_dft_c2r(param.fft_plan_c2r_backward,tmp,localconvFFT);
/****Asigning convolution fftw_complex to bioem_map ****/ /****Asigning convolution fftw_complex to bioem_map ****/
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
{ {
for(int j=0; j < param.param_device.NumberPixels ; j++ ) for(int j=0; j < param.param_device.NumberPixels ; j++ )
{ {
Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j][0]; Mapconv.points[i][j]=localconvFFT[i*param.param_device.NumberPixels+j];
} }
} }
...@@ -620,8 +624,8 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b ...@@ -620,8 +624,8 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
{ {
for(int j=0; j < param.param_device.NumberPixels ; j++ ) for(int j=0; j < param.param_device.NumberPixels ; j++ )
{ {
sumC+=localconvFFT[i*param.param_device.NumberPixels+j][0]; sumC+=localconvFFT[i*param.param_device.NumberPixels+j];
sumsquareC+=localconvFFT[i*param.param_device.NumberPixels+j][0]*localconvFFT[i*param.param_device.NumberPixels+j][0]; sumsquareC+=localconvFFT[i*param.param_device.NumberPixels+j]*localconvFFT[i*param.param_device.NumberPixels+j];
} }
} }
/*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***/ /*** The DTF gives an unnormalized value so have to divded by the total number of pixels in Fourier ***/
...@@ -631,6 +635,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b ...@@ -631,6 +635,7 @@ int bioem::createConvolutedProjectionMap(int iMap,int iConv,mycomplex_t* lproj,b
/**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/ /**** Freeing fftw_complex created (dont know if omp critical is necessary) ****/
myfftw_free(localconvFFT); myfftw_free(localconvFFT);
myfftw_free(tmp);
return(0); return(0);
} }
......
...@@ -30,7 +30,7 @@ public: ...@@ -30,7 +30,7 @@ public:
int createProjection(int iMap, mycomplex_t* map); int createProjection(int iMap, mycomplex_t* map);
int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare); int calcross_cor(bioem_map& localmap,myfloat_t& sum,myfloat_t& sumsquare);
int calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, int value, int disx, int disy); int calProb(int iRefMap,int iOrient, int iConv,myfloat_t sumC,myfloat_t sumsquareC, int value, int disx, int disy);
int calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,mycomplex_t* lCC); int calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC,myfloat_t sumsquareC, mycomplex_t* localConvFFT,mycomplex_t* localCCT,myfloat_t* lCC);
bioem_Probability* pProb; bioem_Probability* pProb;
......
...@@ -10,7 +10,11 @@ typedef float myfloat_t; ...@@ -10,7 +10,11 @@ typedef float myfloat_t;
#define myfftw_destroy_plan fftwf_destroy_plan #define myfftw_destroy_plan fftwf_destroy_plan
#define myfftw_execute fftwf_execute #define myfftw_execute fftwf_execute
#define myfftw_execute_dft fftwf_execute_dft #define myfftw_execute_dft fftwf_execute_dft
#define myfftw_execute_dft_r2c fftwf_execute_dft_r2c
#define myfftw_execute_dft_c2r fftwf_execute_dft_c2r
#define myfftw_plan_dft_2d fftwf_plan_dft_2d #define myfftw_plan_dft_2d fftwf_plan_dft_2d
#define myfftw_plan_dft_r2c_2d fftwf_plan_dft_r2c_2d
#define myfftw_plan_dft_c2r_2d fftwf_plan_dft_c2r_2d
#define myfftw_plan fftwf_plan #define myfftw_plan fftwf_plan
#else #else
typedef double myfloat_t; typedef double myfloat_t;
...@@ -19,7 +23,11 @@ typedef double myfloat_t; ...@@ -19,7 +23,11 @@ typedef double myfloat_t;
#define myfftw_destroy_plan fftw_destroy_plan #define myfftw_destroy_plan fftw_destroy_plan
#define myfftw_execute fftw_execute #define myfftw_execute fftw_execute
#define myfftw_execute_dft fftw_execute_dft #define myfftw_execute_dft fftw_execute_dft
#define myfftw_execute_dft_r2c fftw_execute_dft_r2c
#define myfftw_execute_dft_c2r fftw_execute_dft_c2r
#define myfftw_plan_dft_2d fftw_plan_dft_2d #define myfftw_plan_dft_2d fftw_plan_dft_2d
#define myfftw_plan_dft_r2c_2d fftw_plan_dft_r2c_2d
#define myfftw_plan_dft_c2r_2d fftw_plan_dft_c2r_2d
#define myfftw_plan fftw_plan #define myfftw_plan fftw_plan
#endif #endif
typedef myfloat_t mycomplex_t[2]; typedef myfloat_t mycomplex_t[2];
......
...@@ -17,7 +17,6 @@ class bioem_map_forFFT ...@@ -17,7 +17,6 @@ class bioem_map_forFFT
{ {
public: public:
mycomplex_t cpoints[BIOEM_MAP_SIZE_X*BIOEM_MAP_SIZE_Y]; mycomplex_t cpoints[BIOEM_MAP_SIZE_X*BIOEM_MAP_SIZE_Y];
}; };
class bioem_convolutedMap class bioem_convolutedMap
...@@ -49,7 +48,6 @@ public: ...@@ -49,7 +48,6 @@ public:
__host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);} __host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);}
}; };
class bioem_RefMap_Mod class bioem_RefMap_Mod
{ {
public: public:
......
...@@ -15,6 +15,7 @@ public: ...@@ -15,6 +15,7 @@ public:
int maxDisplaceCenter; int maxDisplaceCenter;
int GridSpaceCenter; int GridSpaceCenter;
int NumberPixels; int NumberPixels;
int NumberFFTPixels1D;
int NtotDist; int NtotDist;
myfloat_t Ntotpi; myfloat_t Ntotpi;
myfloat_t volu; myfloat_t volu;
......
...@@ -9,7 +9,6 @@ ...@@ -9,7 +9,6 @@
#include "map.h" #include "map.h"
#include "param.h" #include "param.h"
using namespace std; using namespace std;
int bioem_RefMap::readRefMaps(bioem_param& param) int bioem_RefMap::readRefMaps(bioem_param& param)
...@@ -134,12 +133,14 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param) ...@@ -134,12 +133,14 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
/********** Routine that pre-calculates Kernels for Convolution **********************/ /********** Routine that pre-calculates Kernels for Convolution **********************/
/************************************************************************************/ /************************************************************************************/
mycomplex_t* localMap; myfloat_t* localMap;
localMap= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
localMap= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels);
RefMapFFT = new bioem_map_forFFT[ntotRefMap]; RefMapFFT = new bioem_map_forFFT[ntotRefMap];
mycomplex_t* localout; mycomplex_t* localout;
localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberPixels); localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param.param_device.NumberPixels*param.param_device.NumberFFTPixels1D);
for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++) for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
{ {
...@@ -148,21 +149,20 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param) ...@@ -148,21 +149,20 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{ {
for(int j=0; j< param.param_device.NumberPixels; j++) for(int j=0; j< param.param_device.NumberPixels; j++)
{ {
localMap[i*param.param_device.NumberPixels+j][0]=Ref[iRefMap].points[i][j]; localMap[i*param.param_device.NumberPixels+j]=Ref[iRefMap].points[i][j];
localMap[i*param.param_device.NumberPixels+j][1]=0.0;
} }
} }
//Calling FFT_Forward //Calling FFT_Forward
myfftw_execute_dft(param.fft_plan_c2c_forward,localMap,localout); myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localMap,localout);
// Normalizing and saving the Reference CTFs // Normalizing and saving the Reference CTFs
for(int i=0; i < param.param_device.NumberPixels ; i++ ) for(int i=0; i < param.param_device.NumberPixels ; i++ )
{ {
for(int j=0; j < param.param_device.NumberPixels ; j++ ) for(int j=0; j < param.param_device.NumberFFTPixels1D ; j++ )
{ {
RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][0]=localout[i*param.param_device.NumberPixels+j][0]; RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][0]=localout[i*param.param_device.NumberFFTPixels1D+j][0];
RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberPixels+j][1]=localout[i*param.param_device.NumberPixels+j][1]; RefMapFFT[iRefMap].cpoints[i*param.param_device.NumberFFTPixels1D+j][1]=localout[i*param.param_device.NumberFFTPixels1D+j][1];
} }
} }
...@@ -172,6 +172,5 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param) ...@@ -172,6 +172,5 @@ int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
myfftw_free(localMap); myfftw_free(localMap);
myfftw_free(localout); myfftw_free(localout);
return(0); return(0);
} }
...@@ -65,7 +65,6 @@ int bioem_model::readModel() ...@@ -65,7 +65,6 @@ int bioem_model::readModel()
float z = 0.0; float z = 0.0;
char tmp[6] = {' '}; char tmp[6] = {' '};
// parse name // parse name
strncpy(tmp,line+12,4); strncpy(tmp,line+12,4);
sscanf (tmp,"%s",name); sscanf (tmp,"%s",name);
...@@ -172,8 +171,6 @@ int bioem_model::readModel() ...@@ -172,8 +171,6 @@ int bioem_model::readModel()
return(0); return(0);
} }
myfloat_t bioem_model::getAminoAcidRad(char *name) myfloat_t bioem_model::getAminoAcidRad(char *name)
{ {
/*************** Function that gets the radius for each amino acid ****************/ /*************** Function that gets the radius for each amino acid ****************/
...@@ -209,7 +206,6 @@ myfloat_t bioem_model::getAminoAcidRad(char *name) ...@@ -209,7 +206,6 @@ myfloat_t bioem_model::getAminoAcidRad(char *name)
} }
myfloat_t bioem_model::getAminoAcidDensity(char *name) myfloat_t bioem_model::getAminoAcidDensity(char *name)
{ {
/*************** Function that gets the number of electrons for each amino acid ****************/ /*************** Function that gets the number of electrons for each amino acid ****************/
...@@ -242,6 +238,5 @@ myfloat_t bioem_model::getAminoAcidDensity(char *name) ...@@ -242,6 +238,5 @@ myfloat_t bioem_model::getAminoAcidDensity(char *name)
exit(0); exit(0);
} }
return iaa; return iaa;
} }
...@@ -15,6 +15,7 @@ bioem_param::bioem_param() ...@@ -15,6 +15,7 @@ bioem_param::bioem_param()
{ {
//Number of Pixels //Number of Pixels
param_device.NumberPixels=0; param_device.NumberPixels=0;
param_device.NumberFFTPixels1D = 0;
// Euler angle grid spacing // Euler angle grid spacing
angleGridPointsAlpha = 0; angleGridPointsAlpha = 0;
angleGridPointsBeta = 0; angleGridPointsBeta = 0;
...@@ -175,6 +176,7 @@ int bioem_param::readParameters() ...@@ -175,6 +176,7 @@ int bioem_param::readParameters()
} }
input.close(); input.close();
param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
cout << " +++++++++++++++++++++++++++++++++++++++++ \n"; cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
cout << "Preparing FFTs\n"; cout << "Preparing FFTs\n";
...@@ -183,10 +185,12 @@ int bioem_param::readParameters() ...@@ -183,10 +185,12 @@ int bioem_param::readParameters()
tmp_map = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels); tmp_map = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels);
tmp_map2 = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels); tmp_map2 = (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels * param_device.NumberPixels);
fft_plan_c2c_forward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_FORWARD, FFTW_MEASURE); fft_plan_c2c_forward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_FORWARD, FFTW_MEASURE | FFTW_DESTROY_INPUT);
fft_plan_c2c_backward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_BACKWARD, FFTW_MEASURE); fft_plan_c2c_backward = myfftw_plan_dft_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, tmp_map2, FFTW_BACKWARD, FFTW_MEASURE | FFTW_DESTROY_INPUT);
fft_plan_r2c_forward = myfftw_plan_dft_r2c_2d(param_device.NumberPixels, param_device.NumberPixels, (myfloat_t*) tmp_map, tmp_map2, FFTW_MEASURE | FFTW_DESTROY_INPUT);
fft_plan_c2r_backward = myfftw_plan_dft_c2r_2d(param_device.NumberPixels, param_device.NumberPixels, tmp_map, (myfloat_t*) tmp_map2, FFTW_MEASURE | FFTW_DESTROY_INPUT);
if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0) if (fft_plan_c2c_forward == 0 || fft_plan_c2c_backward == 0 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
{ {
cout << "Error planing FFTs\n"; cout << "Error planing FFTs\n";
exit(1); exit(1);
...@@ -255,7 +259,6 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS ...@@ -255,7 +259,6 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
} }
int bioem_param::CalculateRefCTF() int bioem_param::CalculateRefCTF()
{ {
/**************************************************************************************/ /**************************************************************************************/
...@@ -263,14 +266,13 @@ int bioem_param::CalculateRefCTF() ...@@ -263,14 +266,13 @@ int bioem_param::CalculateRefCTF()
/************************************************************************************/ /************************************************************************************/
myfloat_t amp,env,phase,ctf,radsq; myfloat_t amp,env,phase,ctf,radsq;
mycomplex_t* localCTF; myfloat_t* localCTF;
int nctfmax= param_device.NumberPixels / 2; int nctfmax= param_device.NumberPixels / 2;
int n=0; int n=0;
localCTF= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param_device.NumberPixels*param_device.NumberPixels); localCTF= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels*param_device.NumberPixels);
refCTF = new bioem_map_forFFT[MAX_REF_CTF]; refCTF = new bioem_map_forFFT[MAX_REF_CTF];
for (int iamp = 0; iamp < numberGridPointsCTF_amp ; iamp++) //Loop over amplitud for (int iamp = 0; iamp < numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
{ {
amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp; amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
...@@ -283,7 +285,7 @@ int bioem_param::CalculateRefCTF() ...@@ -283,7 +285,7 @@ int bioem_param::CalculateRefCTF()
{ {
env= (myfloat_t) ienv * gridEnvelop + startGridEnvelop; env= (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
memset(localCTF,0,param_device.NumberPixels*param_device.NumberPixels*sizeof(mycomplex_t)); memset(localCTF,0,param_device.NumberPixels*param_device.NumberPixels*sizeof(myfloat_t));
//Assigning CTF values //Assigning CTF values
for(int i=0; i< nctfmax; i++) for(int i=0; i< nctfmax; i++)
...@@ -293,25 +295,25 @@ int bioem_param::CalculateRefCTF() ...@@ -293,25 +295,25 @@ int bioem_param::CalculateRefCTF()
radsq=(myfloat_t) (i*i+j*j) * pixelSize * pixelSize; radsq=(myfloat_t) (i*i+j*j) * pixelSize * pixelSize;
ctf=exp(-radsq*env/2.0)*(amp*cos(radsq*phase/2.0)-sqrt((1-amp*amp))*sin(radsq*phase/2.0)); ctf=exp(-radsq*env/2.0)*(amp*cos(radsq*phase/2.0)-sqrt((1-amp*amp))*sin(radsq*phase/2.0));
localCTF[i*param_device.NumberPixels+j][0]=(myfloat_t