Skip to content

GitLab

Commit f83b2d54 authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Quaternions and input organization

parent da09ec2b
Show whitespace changes
Inline Side-by-side
Showing with 264 additions and 25 deletions
bioem.cpp
View file @ f83b2d54
......@@ -129,7 +129,7 @@ int bioem::configure(int ac, char* av[])
("Particlesfile", po::value< std::string>(), "if BioEM (Mandatory) Name of paricles 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 (!)")
("ReadEulerAngles", po::value< std::string>(), "(Optional) Read Euler angle list instead of uniform grid (file nec.)")
("ReadOrientation", po::value< std::string>(), "(Optional) Read orientation list instead of uniform grid (file nec.)")
("ReadPDB", "(Optional) If reading model file in PDB format")
("ReadMRC", "(Optional) If reading particle file in MRC format")
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
......@@ -146,7 +146,7 @@ int bioem::configure(int ac, char* av[])
p.add("ReadPDB", -1);
p.add("ReadMRC", -1);
p.add("ReadMultipleMRC", -1);
p.add("ReadEulerAngles",-1);
p.add("ReadOrientation",-1);
p.add("PrintBestCalMap",-1);
p.add("DumpMaps", -1);
p.add("LoadMapDump", -1);
......@@ -183,12 +183,16 @@ int bioem::configure(int ac, char* av[])
cout << "Reading model file in PDB format.\n";
Model.readPDB = true;
}
if (vm.count("ReadEulerAngles"))
if (vm.count("ReadOrientation"))
{
cout << "Reading Euler Angles from file: "
<< vm["ReadEulerAngles"].as< std::string >() << "\n";
cout << "Note: Format should be alpha (12.6f) beta (12.6f) gamma (12.6f)\n";
Inputanglefile = vm["ReadEulerAngles"].as< std::string >();
cout << "Reading Orientation from file: "
<< vm["ReadOrientation"].as< std::string >() << "\n";
cout << "Important! if using Quaternions, include \n";
cout << "QUATERNIONS keyword in INPUT PARAMETER FILE\n";
cout << "First row in file should be the total number of orientations (int)\n";
cout << "Euler angle format should be alpha (12.6f) beta (12.6f) gamma (12.6f)\n";
cout << "Quaternion format q1 (12.6f) q2 (12.6f) q3 (12.6f) q4 (12.6f)\n";
Inputanglefile = vm["ReadOrientation"].as< std::string >();
param.notuniformangles=true;
}
if (vm.count("PrintBestCalMap"))
......@@ -265,7 +269,7 @@ int bioem::configure(int ac, char* av[])
}
// ********************* Reading Model Input ******************************
Model.readModel(modelfile.c_str());
Model.readModel(param, modelfile.c_str());
cout << "Remark: look at file COORDREAD to confirm that the Model coordinates are correct\n";
......@@ -723,7 +727,7 @@ int bioem::run()
{
bioem_Probability_angle& pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << " " << log(param.param_device.volu) << "\n";
angProbfile << " " << iRefMap << " " << param.angles[iOrient].pos[0] << " " << param.angles[iOrient].pos[1] << " " << param.angles[iOrient].pos[2] << " " << log(pProbAngle.forAngles) + pProbAngle.ConstAngle + 0.5 * log(M_PI) + (1 - param.param_device.Ntotpi * 0.5)*(log(2 * M_PI) + 1) + log(param.param_device.volu) << "\n";
}
}
......@@ -896,6 +900,48 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
localproj = param.fft_scratch_real[omp_get_thread_num()];
memset(localproj, 0, param.param_device.NumberPixels * param.param_device.NumberPixels * sizeof(*localproj));
if(param.doquater){
myfloat_t quater[4],cq[4],temp[4];
myfloat_t a[4];
//quaternion
quater[1]=param.angles[iMap].pos[1];
quater[2]=param.angles[iMap].pos[2];
quater[3]=param.angles[iMap].pos[0];
quater[0]=param.angles[iMap].quat4;
// its conjugate
cq[0]=quater[0];
cq[1]=-quater[1];
cq[2]=-quater[2];
cq[3]=-quater[3];
for(int n = 0; n < Model.nPointsModel; n++)
{
//Rotation with quaterion
// Initial coords
a[0]=0.;
a[1]=Model.points[n].point.pos[0];
a[2]=Model.points[n].point.pos[1];
a[3]=Model.points[n].point.pos[2];
temp[0]=a[0]*cq[0]-a[1]*cq[1]-a[2]*cq[2]-a[3]*cq[3];
temp[1]=a[0]*cq[1]+a[1]*cq[0]+a[2]*cq[3]-a[3]*cq[2];
temp[2]=a[0]*cq[2]-a[1]*cq[3]+a[2]*cq[0]+a[3]*cq[1];
temp[3]=a[0]*cq[3]+a[1]*cq[2]-a[2]*cq[1]+a[3]*cq[0];
// First coordinate is zero Going back to real space
RotatedPointsModel[n].pos[0] = quater[0]*temp[1]+quater[1]*temp[0]+quater[2]*temp[3]-quater[3]*temp[2];
RotatedPointsModel[n].pos[1] = quater[0]*temp[2]-quater[1]*temp[3]+quater[2]*temp[0]+quater[3]*temp[1];
RotatedPointsModel[n].pos[2] = quater[0]*temp[3]+quater[1]*temp[2]-quater[2]*temp[1]+quater[3]*temp[0];
}
} else{
// Doing Euler angles instead of Quaternions
alpha = param.angles[iMap].pos[0];
beta = param.angles[iMap].pos[1];
gam = param.angles[iMap].pos[2];
......@@ -930,6 +976,12 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
}
}
}
if(param.printrotmod) {
for(int n = 0; n < Model.nPointsModel; n++) cout << "ROTATED " << iMap << " " << n <<" "<< RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " << RotatedPointsModel[n].pos[2] << "\n";
}
int i, j;
//********** Projection with radius ***************
......
include/defs.h
View file @ f83b2d54
......@@ -49,6 +49,7 @@ typedef myfloat_t mycomplex_t[2];
struct myfloat3_t
{
myfloat_t pos[BIOEM_FLOAT_3_PHYSICAL_SIZE];
myfloat_t quat4;
};
#ifdef BIOEM_GPUCODE
......
include/model.h
View file @ f83b2d54
......@@ -2,6 +2,7 @@
#define BIOEM_MODEL_H
#include "defs.h"
#include "param.h"
#include <boost/concept_check.hpp>
class bioem_model
......@@ -18,7 +19,7 @@ public:
bioem_model();
~bioem_model();
int readModel(const char* filemodel);
int readModel(bioem_param& param, const char* filemodel);
bool readPDB;
......
include/param.h
View file @ f83b2d54
......@@ -45,6 +45,10 @@ public:
bool doaaradius;
bool writeCTF;
bool ignoreCCoff;
bool nocentermass;
bool printrotmod;
bool readquatlist;
myfloat_t elecwavel;
bioem_param_device param_device;
......@@ -61,12 +65,16 @@ public:
int angleGridPointsAlpha;
int angleGridPointsBeta;
int GridPointsQuatern;
bool doquater;
bool notuniformangles;
int NotUn_angles;
bool withnoise;
myfloat_t stnoise;
std::string inanglef;
std::string quatfile;
int numberGridPointsDisplaceCenter;
// Grid sampling for the convolution kernel
......
main.cpp
View file @ f83b2d54
......@@ -70,7 +70,7 @@ int main(int argc, char* argv[])
}
// ************ Configuration and Pre-calculating necessary objects *****************
if (mpi_rank == 0) printf("Configuring\n");
// if (mpi_rank == 0) printf("Configuring\n");
if (bio->configure(argc, argv) == 0)
{
......
map.cpp
View file @ f83b2d54
......@@ -175,7 +175,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
{
cout << "..." << nummap << "\n";
}
if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
if(lasti+1 != param.param_device.NumberPixels && lastj+1 != param.param_device.NumberPixels && nummap > 0)
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
exit(1);
......
param.cpp
View file @ f83b2d54
......@@ -73,7 +73,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
bool yesMDC = false ;
bool yesGCen = false ;
//Default not flipped micrographs
//Default VALUES
param_device.flipped=false;
param_device.debugterm=false;
param_device.writeCC=false;
......@@ -81,7 +81,11 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
writeCTF=false;
elecwavel=220;
ignoreCCoff=false;
//Assigning name of angles
doquater= false;
nocentermass=false;
printrotmod=false;
readquatlist=false;
doaaradius=true;
//Storing angle file name if existing.
inanglef=std::string(fileangles);
......@@ -143,6 +147,14 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
cout << "Grid points in Cosine ( beta ) " << angleGridPointsBeta << "\n";
yesGPbe= true;
}
else if (strcmp(token, "QUATERNIONS") == 0)
{
token = strtok(NULL, " ");
GridPointsQuatern = int(atoi(token));
doquater= true;
cout << "Gridpoints Quaternions " << GridPointsQuatern << "\n";
}
else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
{
......@@ -258,10 +270,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
ignoreCCoff = true;
cout << "Ignoring Cross-Correlation offset \n";
}
else if (strcmp(token, "PROJECT_RADIUS") == 0) //If projecting CA with amino-acid radius
else if (strcmp(token, "NO_PROJECT_RADIUS") == 0) //If projecting CA with amino-acid radius
{
doaaradius = true;
cout << "Projecting corresponding radius \n";
doaaradius = false;
cout << "Not Projecting corresponding radius \n";
}
else if (strcmp(token, "DEBUG_INDI_PROB_TERM") == 0)//writing out each term of the probability
{
......@@ -292,6 +304,16 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
cout << "Using default electron wave length: 220 (A)\n";
};
}
else if (strcmp(token, "NO_CENTEROFMASS") == 0)//Number of Euler angle tripplets in non uniform Euler angle sampling
{
nocentermass=true;
cout << "BE CAREFUL CENTER OF MASS IS NOT REMOVED \n Calculated images might be out of range \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";
}
}
......@@ -301,9 +323,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);};
if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);};
if(!notuniformangles){
if(!doquater){
if( not ( yesGPal) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);};
if( not ( yesGPbe )) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);};
}
}}
if( not ( yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);};
if( not ( yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);};
if( not ( yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);};
......@@ -316,7 +339,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles)
if( not ( yesMDC ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);};
if( not ( yesGCen ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);};
if( param_device.writeCC && param_device.CCdisplace < 1 ){ cout << "**** INPUT MISSING: Please provide CROSSCOR_DISPLACE \n" ; exit (1);};
if( !param_device.writeCC && param_device.CCwithBayes ){ cout << "**** INPUT MISSING: WRITE_CROSSCOR keyword \n"; exit(1);}
// if( !param_device.writeCC && param_device.CCwithBayes ){ cout << "**** INPUT MISSING: WRITE_CROSSCOR keyword \n"; exit(1);}
if( param_device.writeCC) {if(!param_device.CCwithBayes ){ cout << "Remark:: Not Using Bayesian method to store Cross-Correlation.\n Only Printing out Maximum\n";}
if(param_device.flipped){ cout << "Remark:: Micrographs are Flipped = Particles are white\n";} else { cout << "Remark:: Micrographs are NOT Flipped = Particles are dark\n";}
}
......@@ -475,9 +498,9 @@ int bioem_param::forprintBest(const char* fileinput)
withnoise=true;
cout << "Including noise with standard deviation " << stnoise << "\n";
}
else if (strcmp(token, "DO_AA_RADIUS") == 0) //If projecting CA with amino-acid radius
else if (strcmp(token, "NO_PROJECT_RADIUS") == 0) //If projecting CA with amino-acid radius
{
doaaradius = true;
doaaradius = false;
cout << "Projecting corresponding radius \n";
}
......@@ -571,7 +594,15 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
myfloat_t voluang;
// Euler angle GRID
if(!doquater){
if(!notuniformangles){
cout << "Calculating Grids in Euler Angles\n ";
myfloat_t grid_alpha, cos_grid_beta;
int n = 0;
......@@ -613,14 +644,26 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
char line[512] = {0};
int n=0;
// First line tels the number of rows
input.getline(line, 511);
char tmpVals[36] = {0};
strncpy (tmpVals, line, 12);
sscanf (tmpVals, "%d", &NotUn_angles);
cout << "Number of Euler angles " << NotUn_angles << "\n";
if(NotUn_angles<1) {
cout << "\nNot defined number of Euler angles in INPUT file:" << endl ;
cout << "Use key word: NOT_UNIFORM_TOTAL_ANGS\n";
// cout << "Use key word: NOT_UNIFORM_TOTAL_ANGS\n";
exit(1);
}
delete[] angles;
angles = new myfloat3_t[ NotUn_angles] ;
angles = new myfloat3_t[NotUn_angles] ;
while (!input.eof())
{
......@@ -656,6 +699,140 @@ int bioem_param::CalculateGridsParam() //TO DO FOR QUATERNIONS
voluang= 1./ (myfloat_t) NotUn_angles;
}
} else {
// Analysis with Quaternions
if(!notuniformangles){
cout << "Calculating Grids in Quaterions\n ";
if (GridPointsQuatern < 0 ) { cout << "*** Missing Gridpoints Quaternions \n after QUATERNIONS (int)\n (int)=Number of gridpoins per dimension"; exit(1);}
myfloat_t dgridq,q1,q2,q3;
int n=0;
dgridq=2.f/(myfloat_t) (GridPointsQuatern +1);
// loop to calculate the number ofpoints in the quaternion shpere rad < 1
for (int ialpha = 0; ialpha < GridPointsQuatern + 1 ; ialpha ++)
{
q1=(myfloat_t) ialpha * dgridq -1.f + 0.5 * dgridq;
for (int ibeta = 0; ibeta < GridPointsQuatern + 1 ; ibeta ++)
{
q2=(myfloat_t) ibeta * dgridq -1.f + 0.5 * dgridq;
for (int igamma = 0; igamma < GridPointsQuatern + 1; igamma ++)
{
q3= (myfloat_t) igamma * dgridq -1.f + 0.5 * dgridq;
if(q1*q1+q2*q2+q3*q3 <= 1.f)n=n+2;
}
}
}
//allocating angles
nTotGridAngles = n;
delete[] angles;
angles = new myfloat3_t[nTotGridAngles];
voluang= dgridq * dgridq * dgridq;
n=0;
// assigning values
for (int ialpha = 0; ialpha < GridPointsQuatern + 1; ialpha ++)
{
q1=(myfloat_t) ialpha * dgridq -1.f + 0.5 * dgridq;
for (int ibeta = 0; ibeta < GridPointsQuatern + 1; ibeta ++)
{
q2=(myfloat_t) ibeta * dgridq -1.f + 0.5 * dgridq;
for (int igamma = 0; igamma < GridPointsQuatern + 1 ; igamma ++)
{
q3= (myfloat_t) igamma * dgridq -1.f + 0.5 * dgridq;
if(q1*q1+q2*q2+q3*q3 <= 1.f){
angles[n].pos[0] = q1;
angles[n].pos[1] = q2;
angles[n].pos[2] = q3;
angles[n].quat4=sqrt(1.f-q1*q1-q2*q2-q3*q3);
n++;
//Adding the negative
angles[n].pos[0] = q1;
angles[n].pos[1] = q2;
angles[n].pos[2] = q3;
angles[n].quat4=-sqrt(1.f-q1*q1-q2*q2-q3*q3);
n++;
}
}
}
}
} else{
// ifstream input(quatfile);
ifstream input(inanglef.c_str());
// ifstream input("QUATERNION_LIST");
if (!input.good())
{
cout << "Problem with Quaterion List file " << quatfile.c_str() << " " << endl ;
exit(1);
}
char line[512] = {0};
int n=0;
// First line tels the number of rows
input.getline(line, 511);
int ntotquat;
char tmpVals[36] = {0};
strncpy (tmpVals, line, 12);
sscanf (tmpVals, "%d", &ntotquat);
cout << "Number of quaternions " << ntotquat << "\n";
delete[] angles;
angles = new myfloat3_t[ ntotquat] ;
while (!input.eof())
{
input.getline(line, 511);
if(n< ntotquat){
myfloat_t q1,q2,q3,q4;
char tmpVals[36] = {0};
strncpy (tmpVals, line, 12);
sscanf (tmpVals, "%f", &q1);
strncpy (tmpVals, line + 12, 12);
sscanf (tmpVals, "%f", &q2);
strncpy (tmpVals, line + 24, 12);
sscanf (tmpVals, "%f", &q3);
strncpy (tmpVals, line + 36, 12);
sscanf (tmpVals, "%f", &q4);
angles[n].pos[0] = q1;
angles[n].pos[1] = q2;
angles[n].pos[2] = q3;
angles[n].quat4 = q4;
//cout << "INTT2 " << n << " " << q1 << " " << q2 << " " << q3 << " " << q4 <<"\n";
}
n++;
if(ntotquat+1 < n) {
cout << "More quaternions than expected in header " << n << " instead of " << NotUn_angles << "\n";
exit(1);
}
}
nTotGridAngles = ntotquat;
voluang= 1./ (myfloat_t) ntotquat;
}
cout << "Analysis with Quaternions. Total number of quaternions " << nTotGridAngles << "\n";
}
// ********** Calculating normalized volumen element *********
......
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