diff --git a/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision b/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision
index 42119b0a2737a832cd18594f6a9f25d6840e2cf8..89bcf4e7ec6571f78c536d52a8e09a93450d5347 100644
--- a/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision
+++ b/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision
@@ -8,9 +8,13 @@ USE_QUATERNIONS
##### Constrast transfer integration: #######
CTF_B_ENV 10.0 500. 5
CTF_DEFOCUS 0.5 4.5 8
-CTF_AMPLITUDE 0.01 0.401 4
+CTF_AMPLITUDE 0.01 0.601 5
SIGMA_PRIOR_B_CTF 50.
+SIGMA_PRIOR_DEFOCUS 0.4
+PRIOR_DEFOCUS_CENTER 2.7
+
+
##### Center displacement: #######
DISPLACE_CENTER 40 1
diff --git a/Tutorial_BioEM/Param_ProRun b/Tutorial_BioEM/Param_ProRun
index f359018b9b7165e56588777f9818110174a75454..7ba177358b0fcac31d07d13d8543aaea4bc96752 100644
--- a/Tutorial_BioEM/Param_ProRun
+++ b/Tutorial_BioEM/Param_ProRun
@@ -13,5 +13,10 @@ CTF_AMPLITUDE 0.01 0.401 4.
## Gaussian width of Prior of b-parameter
SIGMA_PRIOR_B_CTF 50.
+##Prior defocus
+##Usefull if defocus from micrographs is known, e.g. at 2.8
+SIGMA_PRIOR_DEFOCUS 0.4
+PRIOR_DEFOCUS_CENTER 2.8
+
##### Center displacement: #######
DISPLACE_CENTER 40 1
diff --git a/bioem.cpp b/bioem.cpp
index 7cfd6703dec712fc03580e8012bb7da20d0649a8..dcc5151af32986a90647bd2d82407ca76f102684 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -135,15 +135,15 @@ int bioem::configure(int ac, char* av[])
po::options_description desc("Command line inputs");
desc.add_options()
("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
- ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of paricles file")
+ ("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of particle-image file")
("Inputfile", po::value<std::string>(), "if BioEM (Mandatory) Name of input parameter file")
- ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map (file nec.). NO BioEM (!)")
- ("ReadOrientation", po::value< std::string>(), "(Optional) Read orientation list instead of uniform grid (file nec.)")
+ ("PrintBestCalMap", po::value< std::string>(), "(Optional) Only print best calculated map. NO BioEM (!)")
+ ("ReadOrientation", po::value< std::string>(), "(Optional) Read file name containing orientations")
("ReadPDB", "(Optional) If reading model file in PDB format")
("ReadMRC", "(Optional) If reading particle file in MRC format")
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
- ("DumpMaps", "(Optional) Dump maps after they were red from maps file")
- ("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
+ ("DumpMaps", "(Optional) Dump maps after they were read from particle-image file")
+ ("LoadMapDump", "(Optional) Read Maps from dump option")
("OutputFile", po::value< std::string>(), "(Optional) For changing the outputfile name")
("help", "(Optional) Produce help message")
;
@@ -1156,6 +1156,8 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
int irad;
myfloat_t dist, rad2;
+ myfloat_t tempden=0.0;
+
for(int n = 0; n < Model.nPointsModel; n++)
{
if(Model.points[n].radius <= param.pixelSize){
@@ -1170,7 +1172,8 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
if(not param.ignorepointsout)exit(1);
}
- localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;//Model.NormDen;
+ localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;
+ tempden += Model.points[n].density;
// exit(1);
}else{
@@ -1204,6 +1207,8 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
{
localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
* 3 / (4 * M_PI * Model.points[n].radius * rad2 );
+ tempden += param.pixelSize * param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
+ * 3 / (4 * M_PI * Model.points[n].radius * rad2 );
}
}
}
@@ -1211,7 +1216,18 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
}
+ // To avoid numerical mismatch in projection errors we normalize by the initial density
+
+ myfloat_t ratioDen;
+
+ ratioDen = Model.NormDen / tempden ;
+ for(int i = 0; i < param.param_device.NumberPixels ; i++){
+ for(int j = 0; j < param.param_device.NumberPixels ; j++){
+ localproj[ i * param.param_device.NumberPixels + j] *= ratioDen;
+ }
+ }
+// cout << "HERE " << tempden << " " << Model.NormDen << " \n";
// **** Output Just to check****
#ifdef DEBUG
diff --git a/model.cpp b/model.cpp
index 252d579e262c65dd63fb4f3212ac3c0d4801c509..aa5dc8db4a2d420c4f5f11cc0a4d1fe1cd4ca2a7 100644
--- a/model.cpp
+++ b/model.cpp
@@ -208,6 +208,7 @@ int bioem_model::readModel(bioem_param& param, const char* filemodel)
}
points = (bioem_model_point*) reallocchk(points, sizeof(bioem_model_point) * nPointsModel);
cout << "Total Number of Voxels " << nPointsModel ;
+ cout << "\nTotal Number of Electrons " << NormDen ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
exampleReadCoor.close();
diff --git a/param.cpp b/param.cpp
index 9158a3200a0f32ecb843288eb97c9b3a9885c462..f037002c24a390aaafc07a4ac7045ae319bd9653 100644
--- a/param.cpp
+++ b/param.cpp
@@ -92,6 +92,7 @@ int bioem_param::readParameters(const char* fileinput)
usepsf=false;
yespriorAngles=false;
ignorepointsout=false;
+ printrotmod=false;
NotUn_angles=0;
priorMod=1; //Default
@@ -411,6 +412,11 @@ int bioem_param::readParameters(const char* fileinput)
ignorepointsout=true;
cout << "Ignoring model points outside the map\n" ;
}
+ else if (strcmp(token, "PRINT_ROTATED_MODELS") == 0)//Number of Euler angle tripplets in non uniform Euler angle sampling
+ {
+ printrotmod=true;
+ cout << "PRINTING out rotatted models (best for debugging)\n";
+ }
}
@@ -435,6 +441,12 @@ int bioem_param::readParameters(const char* fileinput)
if(param_device.writeAngles){cout << "Calculate Cross-cor and write prob angles are mutualy exclusive options\n"; exit(1);}
}
+
+ cout << "To verify input of Priors:\n" ;
+ cout << "Sigma Prior B-Env: " << param_device.sigmaPriorbctf << "\n";
+ cout << "Sigma Prior Defocus: " << param_device.sigmaPriordefo << "\n";
+ cout << "Center Prior Defocus: " <<param_device.Priordefcent <<"\n";
+
// PSF or CTF Checks and asigments
if(usepsf){
if( not ( yesPSFpha ) ){ cout << "**** INPUT MISSING: Please provide Grid PSF PHASE \n" ; exit (1);};