param.cpp

#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;
}