Commit 470535fe authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Proper Correspondance Euler Angles and Quaternions

parent f34833f5
Pipeline #2681 skipped
...@@ -1065,16 +1065,21 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1065,16 +1065,21 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
//*************** Quaternions **************************** //*************** Quaternions ****************************
if(param.doquater){ if(param.doquater){
myfloat_t quater[4],cq[4],temp[4]; myfloat_t quater[4];
myfloat_t a[4]; // myfloat_t cq[4],temp[4];
// myfloat_t a[4];
//quaternion //quaternion
quater[1]=param.angles[iMap].pos[1]; /* quater[1]=param.angles[iMap].pos[1];
quater[2]=param.angles[iMap].pos[2]; quater[2]=param.angles[iMap].pos[2];
quater[3]=param.angles[iMap].pos[0]; quater[3]=param.angles[iMap].pos[0];
quater[0]=param.angles[iMap].quat4; quater[0]=param.angles[iMap].quat4;*/
quater[0]=param.angles[iMap].pos[0];
quater[1]=param.angles[iMap].pos[1];
quater[2]=param.angles[iMap].pos[2];
quater[3]=param.angles[iMap].quat4;
// its conjugate // its conjugate
cq[0]=quater[0]; /* cq[0]=quater[0];
cq[1]=-quater[1]; cq[1]=-quater[1];
cq[2]=-quater[2]; cq[2]=-quater[2];
cq[3]=-quater[3]; cq[3]=-quater[3];
...@@ -1097,12 +1102,25 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1097,12 +1102,25 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
// First coordinate is zero Going back to real space // 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[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[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]; RotatedPointsModel[n].pos[2] = quater[0]*temp[3]+quater[1]*temp[2]-quater[2]*temp[1]+quater[3]*temp[0];
}*/
} //Rotation Matrix for Quaterions (wikipeda)
rotmat[0][0] = 1- 2 * quater[1] * quater[1] - 2 * quater[2] * quater[2];
rotmat[1][0] = 2 * ( quater[0] * quater[1] - quater[2] * quater[3]);
rotmat[2][0] = 2 * ( quater[0] * quater[2] + quater[1] * quater[3]);
rotmat[0][1] = 2 * ( quater[0] * quater[1] + quater[2] * quater[3]);
rotmat[1][1] = 1- 2 * quater[0] * quater[0] - 2 * quater[2] * quater[2];
rotmat[2][1] = 2 * ( quater[1] * quater[2] - quater[0] * quater[3]);
rotmat[0][2] = 2 * ( quater[0] * quater[2] - quater[1] * quater[3]);
rotmat[1][2] = 2 * ( quater[1] * quater[2] + quater[0] * quater[3]);
rotmat[2][2] = 1- 2 * quater[0] * quater[0] - 2 * quater[1] * quater[1];
/*Debug for(int i =0; i <3 ; i++){
for(int j =0; j <3 ; j++){cout << "MAT "<< i << " " << j << " " << rotmat[i][j] << "\n";
} } */
} else{ } else{
//*************** Euler Angles**************************** //*************** Euler Angles****************************
...@@ -1116,6 +1134,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1116,6 +1134,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n"; cout << "Id " << omp_get_thread_num() << " Angs: " << alpha << " " << beta << " " << gam << "\n";
#endif #endif
// ********** Creat Rotation with pre-defiend grid of orientations********** // ********** Creat Rotation with pre-defiend grid of orientations**********
// Same notation as in Goldstein and Mathematica
rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam); rotmat[0][0] = cos(gam) * cos(alpha) - cos(beta) * sin(alpha) * sin(gam);
rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam); rotmat[0][1] = cos(gam) * sin(alpha) + cos(beta) * cos(alpha) * sin(gam);
rotmat[0][2] = sin(gam) * sin(beta); rotmat[0][2] = sin(gam) * sin(beta);
...@@ -1126,7 +1145,10 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1126,7 +1145,10 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
rotmat[2][1] = -sin(beta) * cos(alpha); rotmat[2][1] = -sin(beta) * cos(alpha);
rotmat[2][2] = cos(beta); rotmat[2][2] = cos(beta);
for(int n = 0; n < Model.nPointsModel; n++) }
// The rotation matrix is calculated either for the quaternions or for the euler angles
for(int n = 0; n < Model.nPointsModel; n++)
{ {
RotatedPointsModel[n].pos[0] = 0.0; RotatedPointsModel[n].pos[0] = 0.0;
RotatedPointsModel[n].pos[1] = 0.0; RotatedPointsModel[n].pos[1] = 0.0;
...@@ -1143,7 +1165,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT) ...@@ -1143,7 +1165,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
} }
} }
}
if(param.printrotmod) { 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"; 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";
......
...@@ -918,7 +918,8 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN ...@@ -918,7 +918,8 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
angles[n].pos[2] = (myfloat_t) g; angles[n].pos[2] = (myfloat_t) g;
angles[n].quat4 =0.0;//just to be sure */ angles[n].quat4 =0.0;//just to be sure */
#ifdef DEBUG #ifdef DEBUG
if(yespriorAngles) cout << "check orient: " << n << " " << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n]<< "\n"; // if(yespriorAngles)
cout << "check orient: " << n << " " << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior:\n ";// << angprior[n]<< "\n";
#endif #endif
} }
n++; n++;
...@@ -1075,11 +1076,12 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN ...@@ -1075,11 +1076,12 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
strncpy (tmpVals, line + 48, 12); strncpy (tmpVals, line + 48, 12);
sscanf (tmpVals, "%f", &pp); sscanf (tmpVals, "%f", &pp);
if(pp <0.0000001)cout << "Sure you're input is correct? Very small prior.\n"; if(pp <0.0000001)cout << "Sure you're input is correct? Very small prior.\n";
angprior[n] = pp; angprior[n] = pp;}
#ifdef DEBUG #ifdef DEBUG
if(yespriorAngles) cout << "check orient: " << n << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angprior[n] << "\n"; // if(yespriorAngles)
cout << "check orient: " << n << " " << angles[n].pos[0] << " " << angles[n].pos[1] << " " << angles[n].pos[2] << " prior: " << angles[n].quat4 << "\n";
#endif #endif
}
} }
n++; n++;
if(ntotquat+1 < n) { if(ntotquat+1 < n) {
......
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