-
David Rohr authoredDavid Rohr authored
param.cpp 11.17 KiB
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <cstring>
#include <math.h>
#include <fftw3.h>
#include "param.h"
#include "map.h"
using namespace std;
bioem_param::bioem_param()
{
//Number of Pixels
param_device.NumberPixels=0;
param_device.NumberFFTPixels1D = 0;
// Euler angle grid spacing
angleGridPointsAlpha = 0;
angleGridPointsBeta = 0;
//Envelop function paramters
numberGridPointsEnvelop = 0;
//Contrast transfer function paramters
numberGridPointsCTF_amp = 0;
numberGridPointsCTF_phase = 0;
/****center displacement paramters Equal in both directions***/
param_device.maxDisplaceCenter = 0;
numberGridPointsDisplaceCenter = 0;
fft_plans_created = 0;
}
int bioem_param::readParameters()
{
/**************************************************************************************/
/***************************** Reading Input Parameters ******************************/
/**************************************************************************************/
ifstream input(fileinput);
if (!input.good())
{
cout << "Failed to open file: " << fileinput << "\n";
exit(0);
}
char line[512] = {0};
char saveline[512];
cout << "\n +++++++++++++++++++++++++++++++++++++++++ \n";
cout << "\n READING PARAMETERS \n\n";
cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
while (!input.eof())
{
input.getline(line,512);
strcpy(saveline,line);
char *token = strtok(line," ");
if (token==NULL || line[0] == '#' || strlen(token)==0)
{
// comment or blank line
}
else if (strcmp(token,"PIXEL_SIZE")==0)
{
token = strtok(NULL," ");
pixelSize=atof(token);
cout << "Pixel Sixe " << pixelSize << "\n";
}
else if (strcmp(token,"NUMBER_PIXELS")==0)
{
token = strtok(NULL," ");
param_device.NumberPixels=int(atoi(token));
cout << "Number of Pixels " << param_device.NumberPixels << "\n";
}
else if (strcmp(token,"GRIDPOINTS_ALPHA")==0)
{
token = strtok(NULL," ");
angleGridPointsAlpha=int(atoi(token));
cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
}
else if (strcmp(token,"GRIDPOINTS_BETA")==0)
{
token = strtok(NULL," ");
angleGridPointsBeta=int(atoi(token));
cout << "Grid points beta " << angleGridPointsBeta << "\n";
}
else if (strcmp(token,"GRIDPOINTS_ENVELOPE")==0)
{
token = strtok(NULL," ");
numberGridPointsEnvelop=int(atoi(token));
cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
}
else if (strcmp(token,"START_ENVELOPE")==0)
{
token = strtok(NULL," ");
startGridEnvelop=atof(token);
cout << "Start Envelope " << startGridEnvelop << "\n";
}
else if (strcmp(token,"GRIDSPACE_ENVELOPE")==0)
{
token = strtok(NULL," ");
gridEnvelop=atof(token);
cout << "Grid spacing Envelope " << gridEnvelop << "\n";
}
else if (strcmp(token,"GRIDPOINTS_CTF_PHASE")==0)
{
token = strtok(NULL," ");
numberGridPointsCTF_phase=int(atoi(token));
cout << "Grid points CTF " << numberGridPointsCTF_phase << "\n";
}
else if (strcmp(token,"START_CTF_PHASE")==0)
{
token = strtok(NULL," ");
startGridCTF_phase=atof(token);
cout << "Start CTF " << startGridCTF_phase << "\n";
}
else if (strcmp(token,"GRIDSPACE_CTF_PHASE")==0)
{
token = strtok(NULL," ");
gridCTF_phase=atof(token);
cout << "Grid Space CTF " << gridCTF_phase << "\n";
} else if (strcmp(token,"GRIDPOINTS_CTF_AMP")==0)
{
token = strtok(NULL," ");
numberGridPointsCTF_amp=int(atoi(token));
cout << "Grid points CTF " << numberGridPointsCTF_amp << "\n";
}
else if (strcmp(token,"START_CTF_AMP")==0)
{
token = strtok(NULL," ");
startGridCTF_amp=atof(token);
cout << "Start CTF " << startGridCTF_amp << "\n";
}
else if (strcmp(token,"GRIDSPACE_CTF_AMP")==0)
{
token = strtok(NULL," ");
gridCTF_amp=atof(token);
cout << "Grid Space CTF " << gridCTF_amp << "\n";
}
else if (strcmp(token,"MAX_D_CENTER")==0)
{
token = strtok(NULL," ");
param_device.maxDisplaceCenter=int(atoi(token));
cout << "Maximum displacement Center " << param_device.maxDisplaceCenter << "\n";
}
else if (strcmp(token,"PIXEL_GRID_CENTER")==0)
{
token = strtok(NULL," ");
param_device.GridSpaceCenter=int(atoi(token));
cout << "Grid space displacement center " << param_device.GridSpaceCenter << "\n";
}
else if (strcmp(token,"WRITE_PROB_ANGLES")==0)
{
writeAngles=true;
cout << "Writing Probabilies of each angle \n";
}
}
input.close();
param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
cout << "Preparing FFTs\n";
releaseFFTPlans();
mycomplex_t *tmp_map, *tmp_map2;
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);
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 | 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 || fft_plan_r2c_forward == 0 || fft_plan_c2r_backward == 0)
{
cout << "Error planing FFTs\n";
exit(1);
}
myfftw_free(tmp_map);
myfftw_free(tmp_map2);
fft_plans_created = 1;
cout << " +++++++++++++++++++++++++++++++++++++++++ \n";
return(0);
}
void bioem_param::releaseFFTPlans()
{
if (fft_plans_created)
{
myfftw_destroy_plan(fft_plan_c2c_forward);
myfftw_destroy_plan(fft_plan_c2c_backward);
}
fft_plans_created = 0;
}
int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
{
/**************************************************************************************/
/**************** Routine that pre-calculates Euler angle grids **********************/
/************************************************************************************/
myfloat_t grid_alpha,cos_grid_beta;
int n=0;
//alpha and gamma are uniform in -PI,PI
grid_alpha=2.f * M_PI / (myfloat_t) angleGridPointsAlpha;
//cosine beta is uniform in -1,1
cos_grid_beta=2.f / (myfloat_t) angleGridPointsBeta;
// Euler Angle Array
for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
{
for (int ibeta = 0; ibeta < angleGridPointsBeta; ibeta ++)
{
for (int igamma = 0; igamma < angleGridPointsAlpha; igamma ++)
{
angles[n].pos[0]= (myfloat_t) ialpha * grid_alpha - M_PI + grid_alpha * 0.5f; //ALPHA centered in the middle
angles[n].pos[1]= acos((myfloat_t) ibeta * cos_grid_beta - 1 + cos_grid_beta * 0.5f); //BETA centered in the middle
angles[n].pos[2]= (myfloat_t) igamma * grid_alpha - M_PI + grid_alpha * 0.5f; //GAMMA centered in the middle
n++;
}
}
}
nTotGridAngles=n;
/********** Calculating normalized volumen element *********/
param_device.volu = grid_alpha * grid_alpha * cos_grid_beta * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
* gridCTF_phase * gridCTF_amp * gridEnvelop / (2.f * M_PI) / (2.f * M_PI) / 2.f / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / (2.f * (myfloat_t) param_device.maxDisplaceCenter) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp)
/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase + startGridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop + startGridEnvelop);
/*** Number of total pixels***/
param_device.Ntotpi= (myfloat_t) (param_device.NumberPixels * param_device.NumberPixels);
param_device.NtotDist=(2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1 ) * (2 * (int) (param_device.maxDisplaceCenter/param_device.GridSpaceCenter) + 1);
return(0);
}
int bioem_param::CalculateRefCTF()
{
/**************************************************************************************/
/********** Routine that pre-calculates Kernels for Convolution **********************/
/************************************************************************************/
myfloat_t amp,env,phase,ctf,radsq;
myfloat_t* localCTF;
int nctfmax= param_device.NumberPixels / 2;
int n=0;
localCTF= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param_device.NumberPixels*param_device.NumberPixels);
refCTF = new bioem_map_forFFT[MAX_REF_CTF];
for (int iamp = 0; iamp < numberGridPointsCTF_amp ; iamp++) //Loop over amplitud
{
amp = (myfloat_t) iamp * gridCTF_amp + startGridCTF_amp;
for (int iphase = 0; iphase < numberGridPointsCTF_phase ; iphase++)//Loop over phase
{
phase = (myfloat_t) iphase * gridCTF_phase + startGridCTF_phase;
for (int ienv = 0; ienv < numberGridPointsEnvelop ; ienv++)//Loop over envelope
{
env= (myfloat_t) ienv * gridEnvelop + startGridEnvelop;
memset(localCTF,0,param_device.NumberPixels*param_device.NumberPixels*sizeof(myfloat_t));
//Assigning CTF values
for(int i=0; i< nctfmax; i++)
{
for(int j=0; j< nctfmax; j++)
{
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));
localCTF[i*param_device.NumberPixels+j]=(myfloat_t) ctf;
localCTF[i*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+j]=(myfloat_t) ctf;
localCTF[(param_device.NumberPixels-i-1)*param_device.NumberPixels+param_device.NumberPixels-j-1]=(myfloat_t) ctf;
}
}
//Creatting FFT_Forward of Kernel to store
mycomplex_t* localout;
localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) *param_device.NumberPixels*param_device.NumberFFTPixels1D);
//Calling FFT_Forward
myfftw_execute_dft_r2c(fft_plan_r2c_forward,localCTF,localout);
// Normalizing and saving the Reference CTFs
for(int i=0; i < param_device.NumberPixels ; i++ )
{
for(int j=0; j < param_device.NumberFFTPixels1D ; j++ )
{
refCTF[n].cpoints[i*param_device.NumberFFTPixels1D+j][0]=localout[i*param_device.NumberFFTPixels1D+j][0];
refCTF[n].cpoints[i*param_device.NumberFFTPixels1D+j][1]=localout[i*param_device.NumberFFTPixels1D+j][1];
}
}
CtfParam[n].pos[0]=amp;
CtfParam[n].pos[1]=phase;
CtfParam[n].pos[2]=env;
n++;
myfftw_free(localout);
if(n>MAX_REF_CTF)
{ cout << n << "PROBLEM WITH CTF KERNEL PARAMETERS AND MAX NUMBER ALLOWED\n";
exit(1);
}
}
}
}
myfftw_free(localCTF);
nTotCTFs=n;
return(0);
}
bioem_param::~bioem_param()
{
releaseFFTPlans();
param_device.NumberPixels=0;
angleGridPointsAlpha = 0;
angleGridPointsBeta = 0;
numberGridPointsEnvelop = 0;
numberGridPointsCTF_amp = 0;
numberGridPointsCTF_phase = 0;
param_device.maxDisplaceCenter = 0;
numberGridPointsDisplaceCenter = 0;
}