Skip to content

GitLab

B
BioEM
Commit 797ab249 authored by Pilar Cossio's avatar Pilar Cossio
Browse files
Options

No B-factor, Prior B-envelop, PSF correction

parent 048ebbce
Show whitespace changes
Inline Side-by-side
Showing with 170 additions and 87 deletions
Tutorial_BioEM/Param_Input
View file @ 797ab249
......@@ -7,7 +7,7 @@ GRIDPOINTS_ALPHA 2
GRIDPOINTS_BETA 2
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 2
CTF_B_ENV 100.0 300.5 2
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
......
Tutorial_BioEM/Param_Input_Quat
View file @ 797ab249
......@@ -4,10 +4,9 @@ PIXEL_SIZE 1.77
##### Quaterion grid points: #######
USE_QUATERNIONS
#GRIDPOINTS_QUATERNION 5
##### Constrast transfer integration: #######
CTF_B_FACTOR 100.0 300.5 8
CTF_B_ENV 100.0 300.5 2
CTF_DEFOCUS 1.0 6.0 5
CTF_AMPLITUDE 0.1 0.1 1
......
Tutorial_BioEM/Param_Input_Quaternions deleted 100644 → 0
View file @ 048ebbce
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
QUATERNIONS 4
GRIDPOINTS_ENVELOPE 4
START_ENVELOPE 0.0002
END_ENVELOPE 0.0012
GRIDPOINTS_PSF_PHASE 4
START_PSF_PHASE 0.004
END_PSF_PHASE 0.016
GRIDPOINTS_PSF_AMP 1
START_PSF_AMP 1.
END_PSF_AMP 1.
MAX_D_CENTER 25
PIXEL_GRID_CENTER 2
PROJECT_RADIUS
WRITE_PROB_ANGLES
Tutorial_BioEM/Param_ProRun 0 → 100644
View file @ 797ab249
##### Micrograph Parameters #######
NUMBER_PIXELS 220
PIXEL_SIZE 1.77
##### Using quaterions #######
USE_QUATERNIONS
##### Constrast transfer integration: #######
CTF_B_ENV 10.0 500. 5.
CTF_DEFOCUS 1.0 6.0 10.
CTF_AMPLITUDE 0.1 0.9 5.
## Gaussian width of Prior of b-parameter
SIGMA_PRIOR_B_CTF 50.
##### Center displacement: #######
DISPLACE_CENTER 40 1
Tutorial_BioEM/Quat_list_Small 0 → 100644
View file @ 797ab249
10
0.12635612 0.01592012 0.99114512 0.03758812
0.12337412 0.03159012 0.97867212 0.16119012
0.11845812 0.04676412 0.95084612 0.28226312
0.11168312 0.06120512 0.90810412 0.39890812
0.10315612 0.07468612 0.85111612 0.50929612
0.09301112 0.08699512 0.78077612 0.61169412
0.08140712 0.09793912 0.69818812 0.70449612
0.06852512 0.10734712 0.60464712 0.78624612
0.05456912 0.11507112 0.50162112 0.85566212
0.03975712 0.12099012 0.39072512 0.91165512
bioem.cpp
View file @ 797ab249
......@@ -439,6 +439,7 @@ int bioem::run()
if (mpi_rank == 0) printf("\tInitializing Probabilities\n");
// Inizialzing Probabilites to zero and constant to -Infinity
for (int iRefMap = 0; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
......@@ -479,6 +480,9 @@ int bioem::run()
if(!FFTAlgo){cout << "Cross correlation calculation must be with enviormental variable FFTALGO=1\n"; exit(1);}
}
}
if(!FFTAlgo){cout << "Remark: Not using FFT algorithm. Not using Prior in B-Env.";}
// **************************************************************************************
deviceStartRun();
......@@ -537,7 +541,13 @@ int bioem::run()
// ***************************************************************************************
// *** Comparing each calculated convoluted map with all experimental maps ***
if (DebugOutput >= 2) timer.ResetStart();
compareRefMaps(iOrient, iConv, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
myfloat_t amp,pha,env;
amp=param.CtfParam[iConv].pos[0];
pha=param.CtfParam[iConv].pos[1];
env=param.CtfParam[iConv].pos[2];
compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, conv_mapFFT, sumCONV, sumsquareCONV);
if (DebugOutput >= 2)
{
......@@ -724,15 +734,15 @@ int bioem::run()
if(!param.doquater){
if(param.usepsf){
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";}else{
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - CTF amp - CTF defocus - CTF B-factor - center x - center y - normalization - offsett \n";}
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - alpha[rad] - beta[rad] - gamma[rad] - CTF amp - CTF defocus - CTF B-Env - center x - center y - normalization - offsett \n";}
}else {
if(param.usepsf){
// if( localcc[rx * param.param_device.NumberPixels + ry] <
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - PSF amp - PSF phase - PSF envelope - center x - center y - normalization - offsett \n";
}else{
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - CTF amp - CTF defocus - CTF B-factor - center x - center y - normalization - offsett \n";
outputProbFile << "Notation= RefMap: MapNumber ; Maximizing Param: MaxLogProb - q1 - q2 - q3 - CTF amp - CTF defocus - CTF B-Env - center x - center y - normalization - offsett \n";
}}
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-factor (B ref. Penzeck 2010)\n";
if(param.writeCTF) outputProbFile << " RefMap: MapNumber ; CTFMaxParm: defocus - b-Env (B ref. Penzeck 2010)\n";
outputProbFile <<"************************* HEADER:: NOTATION *******************************************\n\n";
// Loop over reference maps
......@@ -779,7 +789,7 @@ int bioem::run()
// outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel << " [micro-m]; ";
// outputProbFile << param.CtfParam[pProbMap.max.max_prob_conv].pos[1] << param.CtfParam[pProbMap.max.max_prob_conv].pos[2] << denomi ;
outputProbFile << 2*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[1]/denomi/param.elecwavel*0.0001 << " [micro-m] ";
outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi*2. << " [A²] \n";
outputProbFile << 4*M_PI*M_PI*param.CtfParam[pProbMap.max.max_prob_conv].pos[2]/denomi << " [A²] \n";
}
// Writing Individual Angle probabilities
......@@ -905,7 +915,7 @@ int bioem::run()
return(0);
}
int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
int bioem::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
{
//***************************************************************************************
......@@ -917,7 +927,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc
for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
{
const int num = omp_get_thread_num();
calculateCCFFT(iRefMap, iOrient, iConv, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
calculateCCFFT(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC, localmultFFT, param.fft_scratch_complex[num], param.fft_scratch_real[num]);
}
}
else
......@@ -932,7 +942,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc
return(0);
}
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
{
//***************************************************************************************
//***** Calculating cross correlation in FFTALGOrithm *****
......@@ -946,7 +956,7 @@ inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t
myfftw_execute_dft_c2r(param.fft_plan_c2r_backward, localCCT, lCC);
doRefMapFFT(iRefMap, iOrient, iConv, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, lCC, sumC, sumsquareC, pProb, param.param_device, RefMap);
#ifdef PILAR_DEBUG
if (param.param_device.writeCC)
......@@ -1084,7 +1094,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
{
if (DebugOutput >= 0) cout << "WARNING:::: Model Point out of Projection map: " << i << ", " << j << "\n";
// continue;
exit(1);
if(not param.ignorepointsout)exit(1);
}
localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;//Model.NormDen;
......@@ -1107,7 +1117,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
cout << "WARNING: Angle orient " << n << " " << param.angles[iMap].pos[0] << " " << param.angles[iMap].pos[1] << " " << param.angles[iMap].pos[2] << " out " << i << " " << j << "\n";
cout << "WARNING: MPI rank " << mpi_rank <<"\n";
// continue;
exit(1);
if(not param.ignorepointsout)exit(1);
}
......
bioem_cuda.cu
View file @ 797ab249
......@@ -141,12 +141,12 @@ __global__ void multComplexMap(const mycomplex_t* convmap, const mycomplex_t* re
}
}
__global__ void cuDoRefMapsFFT(const int iOrient, const int iConv, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability pProb, const bioem_param_device param, const bioem_RefMap RefMap, const int maxRef, const int Offset)
__global__ void cuDoRefMapsFFT(const int iOrient, const int iConv, const myfloat_t amp, const myfloat_t pha, const myfloat_t env, const myfloat_t* lCC, const myfloat_t sumC, const myfloat_t sumsquareC, bioem_Probability pProb, const bioem_param_device param, const bioem_RefMap RefMap, const int maxRef, const int Offset)
{
if (myBlockIdxX * myBlockDimX + myThreadIdxX >= maxRef) return;
const int iRefMap = myBlockIdxX * myBlockDimX + myThreadIdxX + Offset;
const myfloat_t* mylCC = &lCC[(myBlockIdxX * myBlockDimX + myThreadIdxX) * param.NumberPixels * param.NumberPixels];
doRefMapFFT(iRefMap, iOrient, iConv, mylCC, sumC, sumsquareC, pProb, param, RefMap);
doRefMapFFT(iRefMap, iOrient, iConv, amp, pha, env, mylCC, sumC, sumsquareC, pProb, param, RefMap);
}
template <class T> static inline T divup(T num, T divider) {return((num + divider - 1) / divider);}
......@@ -162,7 +162,7 @@ static inline int ilog2 (int value)
static inline int ilog2(int value) {return 31 - __builtin_clz(value);}
#endif
int bioem_cuda::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
int bioem_cuda::compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
{
if (startMap)
{
......@@ -199,7 +199,7 @@ int bioem_cuda::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map
cout << "Error running CUFFT " << cufftGetErrorStrung(err) << "\n";
exit(1);
}
cuDoRefMapsFFT<<<divup(num, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream[j & 1]>>>(iOrient, iConv, pFFTtmp[j & 1], sumC, sumsquareC, pProb_device, param.param_device, *gpumap, num, i);
cuDoRefMapsFFT<<<divup(num, CUDA_THREAD_COUNT), CUDA_THREAD_COUNT, 0, cudaStream[j & 1]>>>(iOrient, iConv, amp, pha, env, pFFTtmp[j & 1], sumC, sumsquareC, pProb_device, param.param_device, *gpumap, num, i);
}
checkCudaErrors(cudaGetLastError());
if (GPUDualStream)
......@@ -260,7 +260,7 @@ int bioem_cuda::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map
}
if (GPUWorkload < 100)
{
bioem::compareRefMaps(iOrient, iConv, conv_map, localmultFFT, sumC, sumsquareC, maxRef);
bioem::compareRefMaps(iOrient, iConv, amp, pha, env, conv_map, localmultFFT, sumC, sumsquareC, maxRef);
}
if (GPUAsync)
{
......
include/bioem.h
View file @ 797ab249
......@@ -25,14 +25,14 @@ public:
int doProjections(int iMap);
int createConvolutedProjectionMap(int iOreint, int iMap, mycomplex_t* lproj, myfloat_t* Mapconv, mycomplex_t* localmultFFT, myfloat_t& sumC, myfloat_t& sumsquareC);
virtual int compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
virtual int compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
virtual void* malloc_device_host(size_t size);
virtual void free_device_host(void* ptr);
int createProjection(int iMap, mycomplex_t* map);
int calcross_cor(myfloat_t* localmap, myfloat_t& sum, myfloat_t& sumsquare);
void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC);
void calculateCCFFT(int iMap, int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC);
bioem_Probability pProb;
......
include/bioem_cuda_internal.h
View file @ 797ab249
......@@ -17,7 +17,7 @@ public:
bioem_cuda();
virtual ~bioem_cuda();
virtual int compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
virtual int compareRefMaps(int iOrient, int iConv, myfloat_t amp, myfloat_t pha, myfloat_t env, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap = 0);
virtual void* malloc_device_host(size_t size);
virtual void free_device_host(void* ptr);
......
include/param.h
View file @ 797ab249
......@@ -21,6 +21,8 @@ public:
myfloat_t Ntotpi;
myfloat_t volu;
myfloat_t sigmaPriorbctf;
// If to write Probabilities of Angles from Model
bool writeAngles;
bool writeCC;
......@@ -28,6 +30,7 @@ public:
bool debugterm;
int CCdisplace;
bool CCwithBayes;
bool tousepsf;
};
......@@ -54,6 +57,7 @@ public:
bool notsqure;
bool notnormmap;
bool usepsf;
bool ignorepointsout;
myfloat_t elecwavel;
......
map.cpp
View file @ 797ab249
......@@ -54,6 +54,8 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
if(readMultMRC)
{
cout << "Opening File with MRC list names: " << filemap << "\n";
ifstream input(filemap);
if (!input.good())
......@@ -87,13 +89,9 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
{
indifile=strline.c_str();
size_t foundpos= strline.find(".mrc");
size_t foundpos= strline.find("mrc");
size_t endpos = strline.find_last_not_of(" \t");
if(foundpos > endpos){
cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n";
cout << "Warining:::: Are you sure you want to read an MRC? \n";
}
//Reading Multiple MRC
read_MRC(indifile,param);
......@@ -110,11 +108,11 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
string strfilename(filemap);
size_t foundpos= strfilename.find(".mrc");
size_t foundpos= strfilename.find("mrc");
size_t endpos = strfilename.find_last_not_of(" \t");
if(foundpos > endpos){
cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n";
cout << "Warining:::: mrc extension NOT dectected in file name::" << filemap <<" \n";
cout << "Warining:::: Are you sure you want to read an MRC? \n";
}
......
param.cpp
View file @ 797ab249
......@@ -75,9 +75,10 @@ int bioem_param::readParameters(const char* fileinput)
param_device.flipped=false;
param_device.debugterm=false;
param_device.writeCC=false;
param_device.tousepsf=false;
param_device.CCwithBayes=true;
writeCTF=false;
elecwavel=0.0220;
elecwavel=0.019866;
ignoreCCoff=false;
doquater= false;
nocentermass=false;
......@@ -87,6 +88,7 @@ int bioem_param::readParameters(const char* fileinput)
notnormmap=false;
usepsf=false;
yespriorAngles=false;
ignorepointsout=false;
//Storing angle file name if existing.
//f(notuniformangles)inanglef=std::string(fileangles);
......@@ -94,6 +96,7 @@ int bioem_param::readParameters(const char* fileinput)
priorMod=1; //Default
shiftX=0;
shiftY=0;
param_device.sigmaPriorbctf=100.;
ifstream input(fileinput);
......@@ -175,18 +178,18 @@ int bioem_param::readParameters(const char* fileinput)
doquater= true;
}
//CTF PARAMETERS
else if (strcmp(token, "CTF_B_FACTOR") == 0)
else if (strcmp(token, "CTF_B_ENV") == 0)
{
token = strtok(NULL, " ");
startBfactor = atof(token);
if (startBfactor < 0 ) { cout << "*** Error: Negative START B factor "; exit(1);}
if (startBfactor < 0 ) { cout << "*** Error: Negative START B Env "; exit(1);}
token = strtok(NULL, " ");
endBfactor = atof(token);
if (endBfactor < 0 ) { cout << "*** Error: Negative END B factor "; exit(1);}
if (endBfactor < 0 ) { cout << "*** Error: Negative END B Env "; exit(1);}
token = strtok(NULL, " ");
numberGridPointsEnvelop = int(atoi(token));
if (numberGridPointsEnvelop < 0 ) { cout << "*** Error: Negative Number of Grid points Bfactor "; exit(1);}
cout << "Grid CTF B-FACTOR: " << startBfactor << " " << endBfactor <<" " << numberGridPointsEnvelop<< "\n";
if (numberGridPointsEnvelop < 0 ) { cout << "*** Error: Negative Number of Grid points BEnv "; exit(1);}
cout << "Grid CTF B-ENV: " << startBfactor << " " << endBfactor <<" " << numberGridPointsEnvelop<< "\n";
if(startBfactor > endBfactor){ cout << "Error: Grid ill defined END > START\n"; exit(1);};
yesBFact = true;
}
......@@ -235,6 +238,7 @@ int bioem_param::readParameters(const char* fileinput)
else if (strcmp(token, "USE_PSF") == 0)
{
usepsf=true;
param_device.tousepsf=true;
cout << "IMPORTANT: Using Point Spread Function. Thus, all parameters are in Real Space. \n";
}
else if (strcmp(token,"PSF_AMPLITUDE")==0)
......@@ -387,6 +391,17 @@ int bioem_param::readParameters(const char* fileinput)
shiftY=atoi(token);
cout << "Shifting initial model Y by "<< shiftY << "\n" ;
}
else if (strcmp(token, "SIGMA_PRIOR_B_CTF") == 0)
{
token = strtok(NULL, " ");
param_device.sigmaPriorbctf=atof(token);
cout << "Chainging Gaussian width in Prior of Envelope b parameter" << param_device.sigmaPriorbctf << "\n";
}
else if (strcmp(token, "IGNORE_POINTSOUT") == 0)
{
ignorepointsout=true;
cout << "Ignoring model points outside the map\n" ;
}
}
......@@ -418,18 +433,18 @@ int bioem_param::readParameters(const char* fileinput)
if( not ( yesAMP ) ){ cout << "**** INPUT MISSING: Please provide Grid PSF AMPLITUD \n" ; exit (1);};
} else {
cout << "**Note:: Calculation using CTF values (not PSF). If this is not correct then key word: USE_PSF missing in inputfile**\n";
if( not ( yesBFact ) ){ cout << "**** INPUT MISSING: Please provide Grid CTF B-factor \n" ; exit (1);};
if( not ( yesBFact ) ){ cout << "**** INPUT MISSING: Please provide Grid CTF B-ENV \n" ; exit (1);};
if( not ( yesDefocus ) ){ cout << "**** INPUT MISSING: Please provide Grid CTF Defocus \n" ; exit (1);};
if( not ( yesAMP ) ){ cout << "**** INPUT MISSING: Please provide Grid CTF AMPLITUD \n" ; exit (1);};
// Asigning values of phase according to defocus
startGridCTF_phase= startDefocus * M_PI * 2.f * 10000 * elecwavel ;
endGridCTF_phase= endDefocus * M_PI * 2.f * 10000 * elecwavel ;
//Asigning values of envelope according to b-factor
startGridEnvelop = startBfactor / 2.f;
endGridEnvelop = endBfactor / 2.f;
//Asigning values of envelope according to b-envelope no b-factor
startGridEnvelop = startBfactor ;// 2.f;
endGridEnvelop = endBfactor ; // / 2.f;
}
if(elecwavel==0.0220)cout << "Using default electron wave length: 0.0220 (A)\n";
if(elecwavel==0.019688)cout << "Using default electron wave length: 0.019688 (A) of 300kV microscope\n";
param_device.NumberFFTPixels1D = param_device.NumberPixels / 2 + 1;
FFTMapSize = param_device.NumberPixels * param_device.NumberFFTPixels1D;
......@@ -464,13 +479,14 @@ int bioem_param::forprintBest(const char* fileinput)
param_device.writeCC=false;
param_device.CCwithBayes=true;
writeCTF=false;
elecwavel=0.022;
elecwavel=0.019866;
ignoreCCoff=false;
doquater= false;
nocentermass=false;
printrotmod=false;
readquatlist=false;
doaaradius=true;
doquater=true;
shiftX=0;
shiftY=0;
......@@ -537,6 +553,30 @@ int bioem_param::forprintBest(const char* fileinput)
angles[0].pos[2] = atof(token);
cout << "Best Gamma " << angles[0].pos[2] << "\n";
}
else if (strcmp(token, "USE_QUATERNIONS") == 0)
// else if (token=="USE_QUATERNIONS")
{
cout << "Orientations with Quaternions. \n";
doquater= true;
}
else if (strcmp(token, "BEST_Q1") == 0)
{
token = strtok(NULL, " ");
angles[0].pos[0] = atof(token);
cout << "Best q1 " << angles[0].pos[0] << "\n";
}
else if (strcmp(token, "BEST_Q2") == 0)
{
token = strtok(NULL, " ");
angles[0].pos[1] = atof(token);
cout << "Best q2 " << angles[0].pos[1] << "\n";
}
else if (strcmp(token, "BEST_Q3") == 0)
{
token = strtok(NULL, " ");
angles[0].pos[2] = atof(token);
cout << "Best Q3 " << angles[0].pos[2] << "\n";
}
else if (strcmp(token, "USE_PSF") == 0)
{
usepsf=true;
......@@ -562,12 +602,12 @@ int bioem_param::forprintBest(const char* fileinput)
if(startGridCTF_amp <0){cout << "Error Negative Amplitud\n";exit(1);}
cout << "Best Amplitud PSF " << startGridCTF_amp << "\n";
}
else if (strcmp(token, "BEST_CTF_B_FACTOR") == 0)
else if (strcmp(token, "BEST_CTF_B_ENV") == 0)
{
token = strtok(NULL, " ");
startGridEnvelop = atof(token) / 2.f;
if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START B-factor "; exit(1);}
cout << "Best B- factor " << startGridEnvelop << "\n";
startGridEnvelop = atof(token);// / 2.f;
if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START B-Env "; exit(1);}
cout << "Best B- Env " << startGridEnvelop << "\n";
ctfparam=true;
}
else if (strcmp(token,"BEST_CTF_DEFOCUS")==0)
......@@ -639,9 +679,10 @@ int bioem_param::forprintBest(const char* fileinput)
cout << "Shifting initial model Y by "<< shiftY << "\n" ;
}
}
if(doquater)angles[0].quat4=sqrt(1.f-angles[0].pos[0]*angles[0].pos[0]-angles[0].pos[1]*angles[0].pos[1]-angles[0].pos[2]*angles[0].pos[2]);
input.close();
if(usepsf && ctfparam){
......@@ -1011,6 +1052,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
if(n< ntotquat){
myfloat_t q1,q2,q3,q4,pp;
q1=-99999; q2=-99999;q3=-99999;q4=-99999;
char tmpVals[60] = {0};
strncpy (tmpVals, line, 12);
......@@ -1030,6 +1072,12 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
angles[n].pos[2] = q3;
angles[n].quat4 = q4;
if(q1<-1 || q1 > 1){ cout << "Error reading quaterions from list. Value out of range " << q1 << " row " << n << "\n"; exit(1);};
if(q2<-1 || q2 > 1){ cout << "Error reading quaterions from list. Value out of range " << q2 << " row " << n << "\n"; exit(1);};
if(q3<-1 || q3 > 1){ cout << "Error reading quaterions from list. Value out of range " << q3 << " row " << n << "\n"; exit(1);};
if(q4<-1 || q4 > 1){ cout << "Error reading quaterions from list. Value out of range " << q4 << " row " << n << "\n"; exit(1);};
if(yespriorAngles){
strncpy (tmpVals, line + 48, 12);
sscanf (tmpVals, "%f", &pp);
......@@ -1147,31 +1195,41 @@ int bioem_param::CalculateRefCTF()
mycomplex_t* curRef = &refCTF[n * FFTMapSize];
if(usepsf){
// Calculating in real space Point spread function
for(int i = 0; i < nctfmax; i++)
normctf=0.0;
for(int i = 0; i < param_device.NumberPixels; i++)
{
for(int j = 0; j < nctfmax; j++)
for(int j = 0; j < param_device.NumberPixels; j++)
{
radsq = (myfloat_t) (i * i + j * j) * pixelSize * pixelSize;
int ri=0,rj=0;
//Calculating the distance from the periodic center at 0,0
if(i<nctfmax+1){ ri=i; }else{ ri=param_device.NumberPixels-i;};
if(j<nctfmax+1){ rj=j; }else{ rj=param_device.NumberPixels-j;};
radsq = (myfloat_t) ((ri) * (ri) + (rj) *(rj)) * pixelSize * pixelSize;
ctf = exp(-radsq * env / 2.0) * (- amp * cos(radsq * phase / 2.0) - sqrt((1 - amp * amp)) * sin(radsq * phase / 2.0)) ;
normctf += 4.0 * (myfloat_t) ctf ;
localCTF[i * param_device.NumberPixels + j] = (myfloat_t) ctf;
normctf += localCTF[i * param_device.NumberPixels + j];
// cout << "TT " << i << " " << j << " " << localCTF[i * param_device.NumberPixels + j] << "\n";
}
}
//Assigning PSF values
for(int i = 0; i < nctfmax; i++)
//Normalization
for(int i = 0; i < param_device.NumberPixels; i++)
{
for(int j = 0; j < nctfmax; j++)
for(int j = 0; j < param_device.NumberPixels; 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)) / normctf ;
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;
localCTF[i * param_device.NumberPixels + j]= localCTF[i * param_device.NumberPixels + j]/normctf;
}
}
//Calling FFT_Forward
myfftw_execute_dft_r2c(fft_plan_r2c_forward, localCTF, localout);
......@@ -1239,9 +1297,10 @@ int bioem_param::CalculateRefCTF()
// ********** Calculating normalized volumen element *********
if(!printModel){
param_device.volu = voluang * (myfloat_t) param_device.GridSpaceCenter * pixelSize * (myfloat_t) param_device.GridSpaceCenter * pixelSize
* gridCTF_phase * gridCTF_amp * gridEnvelop / ((2.f * (myfloat_t) param_device.maxDisplaceCenter+1.)) / (2.f * (myfloat_t) (param_device.maxDisplaceCenter + 1.)) / ((myfloat_t) numberGridPointsCTF_amp * gridCTF_amp )
/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop );
/ ((myfloat_t) numberGridPointsCTF_phase * gridCTF_phase) / ((myfloat_t) numberGridPointsEnvelop * gridEnvelop ) / sqrt(2 * M_PI ) / param_device.sigmaPriorbctf ;
// cout << "VOLU " << param_device.volu << " " << gridCTF_amp << "\n";
// *** Number of total pixels***
......@@ -1249,6 +1308,7 @@ int bioem_param::CalculateRefCTF()
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);
}
nTotCC = (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1) * (int) ((myfloat_t) param_device.NumberPixels / (myfloat_t) param_device.CCdisplace + 1);
return(0);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment