diff --git a/Tutorial_BioEM/Param_Input b/Tutorial_BioEM/Param_Input
index ced303c77163c9434b08a88135c527fcaeea4932..d903911b96befe0dc950feee198bdbcf135c59df 100644
--- a/Tutorial_BioEM/Param_Input
+++ b/Tutorial_BioEM/Param_Input
@@ -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
diff --git a/Tutorial_BioEM/Param_Input_Quat b/Tutorial_BioEM/Param_Input_Quat
index 92f5f095f16127a8e4f93e8b71bab4179067110d..0c1403e97552c2ac23c7909a955b11d47a987d60 100644
--- a/Tutorial_BioEM/Param_Input_Quat
+++ b/Tutorial_BioEM/Param_Input_Quat
@@ -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
diff --git a/Tutorial_BioEM/Param_Input_Quaternions b/Tutorial_BioEM/Param_Input_Quaternions
deleted file mode 100644
index 217bc99ee931b06a481d01a34e1ef1a66b7bfb12..0000000000000000000000000000000000000000
--- a/Tutorial_BioEM/Param_Input_Quaternions
+++ /dev/null
@@ -1,16 +0,0 @@
-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
diff --git a/Tutorial_BioEM/Param_ProRun b/Tutorial_BioEM/Param_ProRun
new file mode 100644
index 0000000000000000000000000000000000000000..b72a1d2dae5b8d3083b4f8e23ad444081f55767a
--- /dev/null
+++ b/Tutorial_BioEM/Param_ProRun
@@ -0,0 +1,17 @@
+##### 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
diff --git a/Tutorial_BioEM/Quat_list_Small b/Tutorial_BioEM/Quat_list_Small
new file mode 100644
index 0000000000000000000000000000000000000000..bd436cf2685b22b9ad395067c377f14e96e77efc
--- /dev/null
+++ b/Tutorial_BioEM/Quat_list_Small
@@ -0,0 +1,11 @@
+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
diff --git a/bioem.cpp b/bioem.cpp
index 1b788a500b396abfdb8b7c3e17101f7076e21655..81e1eda22f407eaa07819b0fd1f7b607d6bc41ae 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -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);
}
diff --git a/bioem_cuda.cu b/bioem_cuda.cu
index 00f3f6d62993515638352495b1328bef07c8e363..efc7d256dc68674067db90dae6f52a452b672e94 100644
--- a/bioem_cuda.cu
+++ b/bioem_cuda.cu
@@ -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)
{
diff --git a/include/bioem.h b/include/bioem.h
index e17180258a5c4e408766bffd0046dfb500ea6045..d5b7c8e32a516ec8081df71442093f534fc9c044 100644
--- a/include/bioem.h
+++ b/include/bioem.h
@@ -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;
diff --git a/include/bioem_cuda_internal.h b/include/bioem_cuda_internal.h
index 9aaa689d7aac7d445008cae1f2aaf2e6f5508ea1..cc90d9aa7a7b176b83dff7fd8229809bf229cd0b 100644
--- a/include/bioem_cuda_internal.h
+++ b/include/bioem_cuda_internal.h
@@ -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);
diff --git a/include/param.h b/include/param.h
index 2f37607224ad01eb233d298909dea5d5d5ef4c77..52496ac0b6389f2a2def552601b489f34d7e7da0 100644
--- a/include/param.h
+++ b/include/param.h
@@ -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;
diff --git a/map.cpp b/map.cpp
index 6d0caf5e0242f3e64e7aefcefe91de6e30dd263e..198a9a34652d65efea0f433ee3fabf6534ee7155 100644
--- a/map.cpp
+++ b/map.cpp
@@ -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";
}
diff --git a/param.cpp b/param.cpp
index 9a0679830f8b897c51b45feaf7a451be9891eade..3b940d8933c0238fb3a5c73c01b3ecea6974a84c 100644
--- a/param.cpp
+++ b/param.cpp
@@ -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++)
- {
- 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)) ;
- normctf += 4.0 * (myfloat_t) ctf ;
- }
- }
+ normctf=0.0;
+
+ for(int i = 0; i < param_device.NumberPixels; i++)
+ {
+ for(int j = 0; j < param_device.NumberPixels; j++)
+ {
+ int ri=0,rj=0;
- //Assigning PSF 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)) / normctf ;
+ //Calculating the distance from the periodic center at 0,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;
- }
- }
+ 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)) ;
+
+ 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";
+ }
+ }
+
+ //Normalization
+ for(int i = 0; i < param_device.NumberPixels; i++)
+ {
+ for(int j = 0; j < param_device.NumberPixels; j++)
+ {
+ 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);
}