00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include "transform.h"
00037 #include "util.h"
00038 #include "emobject.h"
00039 #include <cctype>
00040 #include <cstring>
00041 #include "symmetry.h"
00042 using namespace EMAN;
00043
00044 #ifdef WIN32
00045 #ifndef M_PI
00046 #define M_PI 3.14159265358979323846
00047 #endif
00048 #endif
00049
00050 #include <algorithm>
00051
00052 #include <gsl_matrix.h>
00053 #include <gsl_blas.h>
00054 #include <gsl_linalg.h>
00055
00056 #include <ostream>
00057 using std::ostream_iterator;
00058
00059
00060
00061 const float Transform::ERR_LIMIT = 0.000001f;
00062
00063 vector<string> Transform::permissable_2d_not_rot;
00064 vector<string> Transform::permissable_3d_not_rot;
00065 map<string,vector<string> > Transform::permissable_rot_keys;
00066
00067 Transform Transform::icos_5_to_2() {
00068 Transform t;
00069 Dict d;
00070 d["type"] = "eman";
00071 d["phi"] = 0;
00072 d["az"] = 90.0f;
00073 d["alt"] = 31.717474;
00085 t.set_rotation(d);
00086 return t;
00087 }
00088
00089 Transform Transform::tet_3_to_2() {
00090 Transform t;
00091 Dict d;
00092 d["type"] = "eman";
00093 d["phi"] = 45.0f;
00094 d["az"] = 0.0f;
00095 d["alt"] = 54.73561f;
00099 t.set_rotation(d);
00100 return t;
00101 }
00102
00103
00104 Transform::Transform()
00105 {
00106 to_identity();
00107 }
00108
00109 Transform::Transform( const Transform& that )
00110 {
00111 *this = that;
00112 }
00113
00114 Transform& Transform::operator=(const Transform& that ) {
00115 if (this != &that ) {
00116 memcpy(matrix,that.matrix,12*sizeof(float));
00117
00118 }
00119 return *this;
00120 }
00121
00122 bool Transform::operator==(const Transform& rhs) const{
00123 if (memcmp(this->matrix, rhs.matrix, 3*4*sizeof(float)) == 0) {
00124 return true;
00125 }
00126 else {
00127 return false;
00128 }
00129 }
00130
00131 bool Transform::operator!=(const Transform& rhs) const{
00132 return !(operator==(rhs));
00133 }
00134
00135 Transform::Transform(const Dict& d) {
00136 to_identity();
00137 set_params(d);
00138 }
00139
00140
00141 Transform::Transform(const float array[12]) {
00142 memcpy(matrix,array,12*sizeof(float));
00143 }
00144
00145 Transform::Transform(const vector<float> array)
00146 {
00147 set_matrix(array);
00148 }
00149
00150 void Transform::set_matrix(const vector<float>& v)
00151 {
00152 if (v.size() != 12 ) throw InvalidParameterException("The construction array must be of size 12");
00153
00154 for(int i=0; i<3; ++i) {
00155 for(int j=0; j<4; ++j) {
00156 matrix[i][j] = v[i*4+j];
00157 }
00158 }
00159 }
00160
00161 void Transform::copy_matrix_into_array(float* const array) const {
00162
00163 int idx = 0;
00164 for(int i=0; i<3; ++i) {
00165 for(int j=0; j<4; ++j) {
00166 array[idx] = matrix[i][j];
00167 idx ++;
00168 }
00169 }
00170 }
00171
00172 vector<float> Transform::get_matrix() const
00173 {
00174 vector<float> ret(12);
00175 for(int i=0; i<3; ++i) {
00176 for(int j=0; j<4; ++j) {
00177 ret[i*4+j] = matrix[i][j];
00178 }
00179 }
00180 return ret;
00181 }
00182
00183 vector<float> Transform::get_matrix_4x4() const
00184 {
00185 vector<float> ret(16);
00186 for(int i=0; i<3; ++i) {
00187 for(int j=0; j<4; ++j) {
00188 ret[i*4+j] = matrix[i][j];
00189 }
00190 }
00191 ret[12] = 0.0;
00192 ret[13] = 0.0;
00193 ret[14] = 0.0;
00194 ret[15] = 1.0;
00195
00196 return ret;
00197 }
00198 void Transform::to_identity()
00199 {
00200
00201 for(int i=0; i<3; ++i) {
00202 for(int j=0; j<4; ++j) {
00203 if(i==j) {
00204 matrix[i][j] = 1;
00205 }
00206 else {
00207 matrix[i][j] = 0;
00208 }
00209 }
00210 }
00211 }
00212
00213 bool Transform::is_identity() const {
00214 for(int i=0; i<3; ++i) {
00215 for(int j=0; j<4; ++j) {
00216 float c = matrix[i][j];
00217 Util::apply_precision(c,ERR_LIMIT);
00218 if(i==j) {
00219 if (c != 1.0) return false;
00220 }
00221 else {
00222 if (c != 0.0) return false;
00223 }
00224 }
00225 }
00226 return true;
00227 }
00228
00229 bool Transform::is_rot_identity() const {
00230 for(int i=0; i<3; ++i) {
00231 for(int j=0; j<3; ++j) {
00232 float c = matrix[i][j];
00233 Util::apply_precision(c,ERR_LIMIT);
00234 if(i==j) {
00235 if (c != 1.0) return false;
00236 }
00237 else {
00238 if (c != 0.0) return false;
00239 }
00240 }
00241 }
00242 return true;
00243 }
00244
00245 void Transform::set_params(const Dict& d) {
00246 detect_problem_keys(d);
00247
00248 if (d.has_key_ci("type") ) set_rotation(d);
00249
00250 if (d.has_key_ci("scale")) {
00251 float scale = static_cast<float>(d.get_ci("scale"));
00252 set_scale(scale);
00253 }
00254
00255 float dx=0,dy=0,dz=0;
00256
00257 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
00258 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
00259 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
00260
00261 if ( dx != 0.0 || dy != 0.0 || dz != 0.0 ) {
00262 set_trans(dx,dy,dz);
00263 }
00264
00265 if (d.has_key_ci("mirror")) {
00266 EMObject e = d.get_ci("mirror");
00267 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
00268 throw InvalidParameterException("Error, mirror must be a bool or an int");
00269
00270 bool mirror = static_cast<bool>(e);
00271 set_mirror(mirror);
00272 }
00273 }
00274
00275
00276 void Transform::init_permissable_keys()
00277 {
00278
00279 permissable_2d_not_rot.push_back("tx");
00280 permissable_2d_not_rot.push_back("ty");
00281 permissable_2d_not_rot.push_back("scale");
00282 permissable_2d_not_rot.push_back("mirror");
00283 permissable_2d_not_rot.push_back("type");
00284
00285 permissable_3d_not_rot.push_back("tx");
00286 permissable_3d_not_rot.push_back("ty");
00287 permissable_3d_not_rot.push_back("tz");
00288 permissable_3d_not_rot.push_back("scale");
00289 permissable_3d_not_rot.push_back("mirror");
00290 permissable_3d_not_rot.push_back("type");
00291
00292 vector<string> tmp;
00293 tmp.push_back("alpha");
00294 permissable_rot_keys["2d"] = tmp;
00295
00296 tmp.clear();
00297 tmp.push_back("alt");
00298 tmp.push_back("az");
00299 tmp.push_back("phi");
00300 permissable_rot_keys["eman"] = tmp;
00301
00302 tmp.clear();
00303 tmp.push_back("psi");
00304 tmp.push_back("theta");
00305 tmp.push_back("phi");
00306 permissable_rot_keys["spider"] = tmp;
00307
00308 tmp.clear();
00309 tmp.push_back("alpha");
00310 tmp.push_back("beta");
00311 tmp.push_back("gamma");
00312 permissable_rot_keys["imagic"] = tmp;
00313
00314 tmp.clear();
00315 tmp.push_back("ztilt");
00316 tmp.push_back("xtilt");
00317 tmp.push_back("ytilt");
00318 permissable_rot_keys["xyz"] = tmp;
00319
00320 tmp.clear();
00321 tmp.push_back("phi");
00322 tmp.push_back("theta");
00323 tmp.push_back("omega");
00324 permissable_rot_keys["mrc"] = tmp;
00325
00326 tmp.clear();
00327 tmp.push_back("e0");
00328 tmp.push_back("e1");
00329 tmp.push_back("e2");
00330 tmp.push_back("e3");
00331 permissable_rot_keys["quaternion"] = tmp;
00332
00333 tmp.clear();
00334 tmp.push_back("n1");
00335 tmp.push_back("n2");
00336 tmp.push_back("n3");
00337 tmp.push_back("Omega");
00338 permissable_rot_keys["spin"] = tmp;
00339
00340 tmp.clear();
00341 tmp.push_back("n1");
00342 tmp.push_back("n2");
00343 tmp.push_back("n3");
00344 tmp.push_back("q");
00345 permissable_rot_keys["sgirot"] = tmp;
00346
00347 tmp.clear();
00348 tmp.push_back("m11");
00349 tmp.push_back("m12");
00350 tmp.push_back("m13");
00351 tmp.push_back("m21");
00352 tmp.push_back("m22");
00353 tmp.push_back("m23");
00354 tmp.push_back("m31");
00355 tmp.push_back("m32");
00356 tmp.push_back("m33");
00357 permissable_rot_keys["matrix"] = tmp;
00358 }
00359
00360
00361 void Transform::detect_problem_keys(const Dict& d) {
00362 if (permissable_rot_keys.size() == 0 ) {
00363 init_permissable_keys();
00364 }
00365
00366 vector<string> verification;
00367 vector<string> problem_keys;
00368 bool is_2d = false;
00369 if (d.has_key_ci("type") ) {
00370 string type = Util::str_to_lower((string)d["type"]);
00371 bool problem = false;
00372 if (permissable_rot_keys.find(type) == permissable_rot_keys.end() ) {
00373 problem_keys.push_back(type);
00374 problem = true;
00375 }
00376 if ( !problem ) {
00377 vector<string> perm = permissable_rot_keys[type];
00378 std::copy(perm.begin(),perm.end(),back_inserter(verification));
00379
00380 if ( type == "2d" ) {
00381 is_2d = true;
00382 std::copy(permissable_2d_not_rot.begin(),permissable_2d_not_rot.end(),back_inserter(verification));
00383 }
00384 }
00385 }
00386 if ( !is_2d ) {
00387 std::copy(permissable_3d_not_rot.begin(),permissable_3d_not_rot.end(),back_inserter(verification));
00388 }
00389
00390 for (Dict::const_iterator it = d.begin(); it != d.end(); ++it) {
00391 if ( std::find(verification.begin(),verification.end(), it->first) == verification.end() ) {
00392 problem_keys.push_back(it->first);
00393 }
00394 }
00395
00396 if (problem_keys.size() != 0 ) {
00397 string error;
00398 if (problem_keys.size() == 1) {
00399 error = "Transform Error: The \"" +problem_keys[0]+ "\" key is unsupported";
00400 } else {
00401 error = "Transform Error: The ";
00402 for(vector<string>::const_iterator cit = problem_keys.begin(); cit != problem_keys.end(); ++cit ) {
00403 if ( cit != problem_keys.begin() ) {
00404 if (cit == (problem_keys.end() -1) ) error += " and ";
00405 else error += ", ";
00406 }
00407 error += "\"";
00408 error += *cit;
00409 error += "\"";
00410 }
00411 error += " keys are unsupported";
00412 }
00413 throw InvalidParameterException(error);
00414 }
00415 }
00416
00417 void Transform::set_params_inverse(const Dict& d) {
00418 detect_problem_keys(d);
00419
00420 if (d.has_key_ci("type") ) set_rotation(d);
00421
00422 float dx=0,dy=0,dz=0;
00423 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
00424 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
00425 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
00426
00427 if ( (dx != 0.0 || dy != 0.0 || dz != 0.0) && d.has_key_ci("type") ) {
00428 Transform pre_trans;
00429 pre_trans.set_trans(dx,dy,dz);
00430
00431 Transform tmp;
00432 tmp.set_rotation(d);
00433
00434 if (d.has_key_ci("scale")) {
00435 float scale = static_cast<float>(d.get_ci("scale"));
00436 tmp.set_scale(scale);
00437 }
00438
00439 Transform solution_trans = tmp*pre_trans;
00440
00441 if (d.has_key_ci("scale")) {
00442 Transform tmp;
00443 float scale = static_cast<float>(d.get_ci("scale"));
00444 tmp.set_scale(scale);
00445 solution_trans = solution_trans*tmp;
00446 }
00447
00448 tmp = Transform();
00449 tmp.set_rotation(d);
00450 solution_trans = solution_trans*tmp;
00451 set_trans(solution_trans.get_trans());
00452 }
00453
00454 if (d.has_key_ci("scale")) {
00455 float scale = static_cast<float>(d.get_ci("scale"));
00456 set_scale(scale);
00457 }
00458
00459 if (d.has_key_ci("mirror")) {
00460 EMObject e = d.get_ci("mirror");
00461 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
00462 throw InvalidParameterException("Error, mirror must be a bool or an int");
00463
00464 bool mirror = static_cast<bool>(e);
00465 set_mirror(mirror);
00466 }
00467 invert();
00468 }
00469
00470
00471 Dict Transform::get_params(const string& euler_type) const {
00472 Dict params = get_rotation(euler_type);
00473
00474 Vec3f v = get_trans();
00475 params["tx"] = v[0]; params["ty"] = v[1];
00476
00477 string type = Util::str_to_lower(euler_type);
00478 if ( type != "2d") params["tz"] = v[2];
00479
00480 float scale = get_scale();
00481 params["scale"] = scale;
00482
00483 bool mirror = get_mirror();
00484 params["mirror"] = mirror;
00485
00486 return params;
00487 }
00488
00489
00490
00491 Dict Transform::get_params_inverse(const string& euler_type) const {
00492 Transform inv(inverse());
00493
00494 Dict params = inv.get_rotation(euler_type);
00495 Vec3f v = inv.get_pre_trans();
00496 params["tx"] = v[0]; params["ty"] = v[1];
00497
00498 string type = Util::str_to_lower(euler_type);
00499 if ( type != "2d") params["tz"] = v[2];
00500
00501 float scale = inv.get_scale();
00502 params["scale"] = scale;
00503
00504 bool mirror = inv.get_mirror();
00505 params["mirror"] = mirror;
00506
00507 return params;
00508 }
00509
00510
00511 void Transform::set_rotation(const Dict& rotation)
00512 {
00513 detect_problem_keys(rotation);
00514 string euler_type;
00515
00516 if (!rotation.has_key_ci("type") ){
00517 throw InvalidParameterException("argument dictionary does not contain the type key");
00518 }
00519
00520 euler_type = static_cast<string>(rotation.get_ci("type"));
00521
00522
00523 double e0=0; double e1=0; double e2=0; double e3=0;
00524 double Omega=0;
00525 double az = 0;
00526 double alt = 0;
00527 double phi = 0;
00528 double cxtilt = 0;
00529 double sxtilt = 0;
00530 double cytilt = 0;
00531 double sytilt = 0;
00532 double cztilt = 0;
00533 double sztilt = 0;
00534 bool is_quaternion = 0;
00535 bool is_matrix = 0;
00536 bool is_xyz = 0;
00537
00538 bool x_mirror;
00539 float scale;
00540
00541 get_scale_and_mirror(scale,x_mirror);
00542 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
00543
00544 string type = Util::str_to_lower(euler_type);
00545 if (type == "2d") {
00546 assert_valid_2d();
00547 az = 0;
00548 alt = 0;
00549 phi = (double)rotation["alpha"] ;
00550 } else if ( type == "eman" ) {
00551
00552 az = (double)rotation["az"] ;
00553 alt = (double)rotation["alt"] ;
00554 phi = (double)rotation["phi"] ;
00555 } else if ( type == "imagic" ) {
00556
00557 az = (double)rotation["alpha"] ;
00558 alt = (double)rotation["beta"] ;
00559 phi = (double)rotation["gamma"] ;
00560 } else if ( type == "spider" ) {
00561
00562 az = (double)rotation["phi"] + 90.0;
00563 alt = (double)rotation["theta"] ;
00564 phi = (double)rotation["psi"] - 90.0;
00565 } else if ( type == "xyz" ) {
00566
00567 is_xyz = 1;
00568 cxtilt = cos(EMConsts::deg2rad*(double)rotation["xtilt"]);
00569 sxtilt = sin(EMConsts::deg2rad*(double)rotation["xtilt"]);
00570 cytilt = cos(EMConsts::deg2rad*(double)rotation["ytilt"]);
00571 sytilt = sin(EMConsts::deg2rad*(double)rotation["ytilt"]);
00572 cztilt = cos(EMConsts::deg2rad*(double)rotation["ztilt"]);
00573 sztilt = sin(EMConsts::deg2rad*(double)rotation["ztilt"]);
00574 } else if ( type == "mrc" ) {
00575
00576 az = (double)rotation["phi"] + 90.0f ;
00577 alt = (double)rotation["theta"] ;
00578 phi = (double)rotation["omega"] - 90.0f ;
00579 } else if ( type == "quaternion" ) {
00580
00581 is_quaternion = 1;
00582 e0 = (double)rotation["e0"];
00583 e1 = (double)rotation["e1"];
00584 e2 = (double)rotation["e2"];
00585 e3 = (double)rotation["e3"];
00586 } else if ( type == "spin" ) {
00587
00588 is_quaternion = 1;
00589 Omega = (double)rotation["Omega"];
00590 e0 = cos(Omega*EMConsts::deg2rad/2.0);
00591 e1 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"];
00592 e2 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"];
00593 e3 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"];
00594 } else if ( type == "sgirot" ) {
00595
00596 is_quaternion = 1;
00597 Omega = (double)rotation["q"] ;
00598 e0 = cos(Omega*EMConsts::deg2rad/2.0);
00599 e1 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"];
00600 e2 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"];
00601 e3 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"];
00602 } else if ( type == "matrix" ) {
00603 is_matrix = 1;
00604 matrix[0][0] = (float)rotation["m11"];
00605 matrix[0][1] = (float)rotation["m12"];
00606 matrix[0][2] = (float)rotation["m13"];
00607 matrix[1][0] = (float)rotation["m21"];
00608 matrix[1][1] = (float)rotation["m22"];
00609 matrix[1][2] = (float)rotation["m23"];
00610 matrix[2][0] = (float)rotation["m31"];
00611 matrix[2][1] = (float)rotation["m32"];
00612 matrix[2][2] = (float)rotation["m33"];
00613 } else {
00614
00615 throw InvalidStringException(euler_type, "unknown Euler Type");
00616 }
00617
00618 double azp = az*EMConsts::deg2rad;
00619 double altp = alt*EMConsts::deg2rad;
00620 double phip = phi*EMConsts::deg2rad;
00621
00622 if (!is_quaternion && !is_matrix && !is_xyz) {
00623 matrix[0][0] = (float)(cos(phip)*cos(azp) - cos(altp)*sin(azp)*sin(phip));
00624 matrix[0][1] = (float)(cos(phip)*sin(azp) + cos(altp)*cos(azp)*sin(phip));
00625 matrix[0][2] = (float)(sin(altp)*sin(phip));
00626 matrix[1][0] = (float)(-sin(phip)*cos(azp) - cos(altp)*sin(azp)*cos(phip));
00627 matrix[1][1] = (float)(-sin(phip)*sin(azp) + cos(altp)*cos(azp)*cos(phip));
00628 matrix[1][2] = (float)(sin(altp)*cos(phip));
00629 matrix[2][0] = (float)(sin(altp)*sin(azp));
00630 matrix[2][1] = (float)(-sin(altp)*cos(azp));
00631 matrix[2][2] = (float)cos(altp);
00632 }
00633 if (is_quaternion){
00634 matrix[0][0] = (float)(e0 * e0 + e1 * e1 - e2 * e2 - e3 * e3);
00635 matrix[0][1] = (float)(2.0f * (e1 * e2 + e0 * e3));
00636 matrix[0][2] = (float)(2.0f * (e1 * e3 - e0 * e2));
00637 matrix[1][0] = (float)(2.0f * (e2 * e1 - e0 * e3));
00638 matrix[1][1] = (float)(e0 * e0 - e1 * e1 + e2 * e2 - e3 * e3);
00639 matrix[1][2] = (float)(2.0f * (e2 * e3 + e0 * e1));
00640 matrix[2][0] = (float)(2.0f * (e3 * e1 + e0 * e2));
00641 matrix[2][1] = (float)(2.0f * (e3 * e2 - e0 * e1));
00642 matrix[2][2] = (float)(e0 * e0 - e1 * e1 - e2 * e2 + e3 * e3);
00643
00644 }
00645 if (is_xyz){
00646 matrix[0][0] = (float)(cytilt*cztilt);
00647 matrix[0][1] = (float)(cxtilt*sztilt+sxtilt*sytilt*cztilt);
00648 matrix[0][2] = (float)(sxtilt*sztilt-cxtilt*sytilt*cztilt);
00649 matrix[1][0] = (float)(-cytilt*sztilt);
00650 matrix[1][1] = (float)(cxtilt*cztilt-sxtilt*sytilt*sztilt);
00651 matrix[1][2] = (float)(sxtilt*cztilt+cxtilt*sytilt*sztilt);
00652 matrix[2][0] = (float)(sytilt);
00653 matrix[2][1] = (float)(-sxtilt*cytilt);
00654 matrix[2][2] = (float)(cxtilt*cytilt);
00655 }
00656
00657
00658 if (scale != 1.0f) {
00659 for(int i=0; i<3; ++i) {
00660 for(int j=0; j<3; ++j) {
00661 matrix[i][j] *= scale;
00662 }
00663 }
00664 }
00665
00666
00667 if ( x_mirror ) {
00668 for(int j=0; j<3; ++j) {
00669 matrix[0][j] *= -1.0f;
00670 }
00671 }
00672 }
00673
00674 Transform Transform::get_rotation_transform() const
00675 {
00676 Transform ret(*this);
00677 ret.set_scale(1.0);
00678 ret.set_mirror(false);
00679 ret.set_trans(0,0,0);
00680
00681 return ret;
00682 }
00683
00684 void Transform::set_rotation(const Vec3f & v)
00685 {
00686 if ( v[0] == 0 && v[1] == 0 && v[2] == 0 )
00687 throw UnexpectedBehaviorException("Can't set rotation for the null vector");
00688
00689 Vec3f v1(v);
00690 v1.normalize();
00691
00692 double theta = acos(v1[2]);
00693 double psi = atan2(v1[1],-v1[0]);
00694
00695 Dict d;
00696 d["theta"] = (double)EMConsts::rad2deg*theta;
00697 d["psi"] = (double)EMConsts::rad2deg*psi;
00698 d["phi"] = (double)0.0;
00699 d["type"] = "spider";
00700
00701 set_rotation(d);
00702
00703
00704 }
00705
00706 void Transform::rotate_origin(const Transform& by)
00707 {
00708 vector<float> multmatrix = by.get_matrix();
00709
00710 Transform result;
00711 for (int i=0; i<3; i++) {
00712 for (int j=0; j<3; j++) {
00713 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
00714 }
00715 }
00716
00717 for (int i=0; i<3; i++) {
00718 for (int j=0; j<3; j++) {
00719 matrix[i][j] = result[i][j];
00720 }
00721 }
00722 }
00723
00724 void Transform::rotate_origin_newBasis(const Transform& tcs, const float& Omega, const float& n1, const float& n2, const float& n3)
00725 {
00726
00727 Transform tcsinv = Transform(tcs);
00728 tcsinv.set_trans(0.0, 0.0, 0.0);
00729 tcsinv.set_scale(1.0);
00730 tcsinv.invert();
00731
00732
00733 Transform temp = Transform();
00734 temp.set_trans(n1, n2, n3);
00735 Transform cc = tcsinv*temp;
00736 Vec3f cctrans = cc.get_trans();
00737
00738
00739 Dict spinrot = Dict();
00740 spinrot["type"] = "spin";
00741 spinrot["Omega"] = Omega;
00742 spinrot["n1"] = cctrans[0];
00743 spinrot["n2"] = cctrans[1];
00744 spinrot["n3"] = cctrans[2];
00745 Transform rightrot = Transform(spinrot);
00746 rotate_origin(rightrot);
00747 }
00748
00749 void Transform::rotate(const Transform& by)
00750 {
00751 vector<float> multmatrix = by.get_matrix();
00752
00753 Transform result;
00754 for (int i=0; i<3; i++) {
00755 for (int j=0; j<4; j++) {
00756 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
00757 }
00758 }
00759
00760 for (int i=0; i<3; i++) {
00761 for (int j=0; j<4; j++) {
00762 matrix[i][j] = result[i][j];
00763 }
00764 }
00765 }
00766
00767 Transform Transform::negate() const
00768 {
00769 Transform t(*this);
00770 for(unsigned int i = 0; i < 3; ++i) {
00771 for(unsigned int j = 0; j < 4; ++j) {
00772 t.set(i,j,t[i][j]*-1);
00773 }
00774 }
00775 return t;
00776 }
00777
00778 Transform Transform::get_hflip_transform() const {
00779
00780 Dict rot = get_rotation("eman");
00781 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
00782 rot["phi"] = 180.0f - static_cast<float>(rot["phi"]);
00783
00784 Transform ret(*this);
00785 ret.set_rotation(rot);
00786
00787 Vec3f trans = get_trans();
00788 trans[0] = -trans[0];
00789 ret.set_trans(trans);
00790
00791
00792
00793 return ret;
00794 }
00795
00796 Transform Transform::get_vflip_transform() const {
00797
00798 Dict rot = get_rotation("eman");
00799 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
00800 rot["phi"] = - static_cast<float>(rot["phi"]);
00801
00802 Transform ret(*this);
00803 ret.set_rotation(rot);
00804
00805 Vec3f trans = get_trans();
00806 trans[1] = -trans[1];
00807 ret.set_trans(trans);
00808
00809 return ret;
00810 }
00811
00812 Dict Transform::get_rotation(const string& euler_type) const
00813 {
00814 Dict result;
00815
00816
00817 float scale;
00818 bool x_mirror;
00819 get_scale_and_mirror(scale,x_mirror);
00820 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
00821
00822 double cosalt = matrix[2][2]/scale;
00823 double x_mirror_scale = (x_mirror ? -1.0f : 1.0f);
00824 double inv_scale = 1.0f/scale;
00825
00826 double az = 0;
00827 double alt = 0;
00828 double phi = 0;
00829 double phiS = 0;
00830 double psiS = 0;
00831
00832
00833
00834 if (cosalt >= 1) {
00835 alt = 0;
00836 az = 0;
00837 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00838 } else if (cosalt <= -1) {
00839 alt = 180;
00840 az = 0;
00841 phi = (double)EMConsts::rad2deg * atan2(-x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00842 } else {
00843
00844 az = (double)EMConsts::rad2deg * atan2(scale*matrix[2][0], -scale*matrix[2][1]);
00845
00846 if (matrix[2][2]==0.0)
00847 alt = 90.0;
00848 else
00849 alt = (double)EMConsts::rad2deg * atan(sqrt((double)matrix[2][0]*matrix[2][0]+(double)matrix[2][1]*matrix[2][1])/fabs(matrix[2][2]));
00850
00851 if (matrix[2][2] * scale < 0)
00852 alt = 180.0f-alt;
00853
00854 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*(double)matrix[0][2], (double)matrix[1][2]);
00855
00856 }
00857
00858 phi = phi-360.0*floor(phi/360.0);
00859 az = az -360.0*floor(az/360.0);
00860
00861
00862 if (cosalt >= 1) {
00863 phiS = 0;
00864 psiS = phi;
00865 } else if (cosalt <= -1) {
00866 phiS = 0;
00867 psiS = phi + 180.0;
00868 } else {
00869 phiS = az - 90.0;
00870 psiS = phi + 90.0;
00871 }
00872
00873 phiS = phiS-360.0*floor(phiS/360.0);
00874 psiS = psiS-360.0*floor(psiS/360.0);
00875
00876
00877
00878 double nphi = (az-phi)/2.0;
00879
00880 double cosOover2 = cos((az+phi)*EMConsts::deg2rad/2.0) * cos(alt*EMConsts::deg2rad/2.0);
00881 double sinOover2 = sqrt(1 -cosOover2*cosOover2);
00882 double cosnTheta = sin((az+phi)*EMConsts::deg2rad/2.0) * cos(alt*EMConsts::deg2rad/2.0) / sqrt(1-cosOover2*cosOover2);
00883 double sinnTheta = sqrt(1-cosnTheta*cosnTheta);
00884 double n1 = sinnTheta*cos(nphi*EMConsts::deg2rad);
00885 double n2 = sinnTheta*sin(nphi*EMConsts::deg2rad);
00886 double n3 = cosnTheta;
00887 double xtilt = 0;
00888 double ytilt = 0;
00889 double ztilt = 0;
00890
00891
00892 if (cosOover2<0) {
00893 cosOover2*=-1; n1 *=-1; n2*=-1; n3*=-1;
00894 }
00895
00896 string type = Util::str_to_lower(euler_type);
00897
00898 result["type"] = type;
00899 if (type == "2d") {
00900 assert_valid_2d();
00901 result["alpha"] = phi;
00902 } else if (type == "eman") {
00903
00904 result["az"] = az;
00905 result["alt"] = alt;
00906 result["phi"] = phi;
00907 } else if (type == "imagic") {
00908
00909 result["alpha"] = az;
00910 result["beta"] = alt;
00911 result["gamma"] = phi;
00912 } else if (type == "spider") {
00913
00914 result["phi"] = phiS;
00915 result["theta"] = alt;
00916 result["psi"] = psiS;
00917 } else if (type == "mrc") {
00918
00919 result["phi"] = phiS;
00920 result["theta"] = alt;
00921 result["omega"] = psiS;
00922 } else if (type == "xyz") {
00923
00924 xtilt = atan2(-sin(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt),cos(EMConsts::deg2rad*alt));
00925 ytilt = asin( cos(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt));
00926 ztilt = psiS*EMConsts::deg2rad - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
00927
00928 xtilt *= EMConsts::rad2deg; ytilt *= EMConsts::rad2deg; ztilt *= EMConsts::rad2deg;
00929 xtilt = xtilt-360*.0*floor((xtilt+180.0)/360.0);
00930 ytilt = ytilt-360*.0*floor((ytilt+180.0)/360.0);
00931 ztilt = ztilt-360*.0*floor((ztilt+180.0)/360.0);
00932
00933 result["xtilt"] = xtilt;
00934 result["ytilt"] = ytilt;
00935 result["ztilt"] = ztilt;
00936 } else if (type == "quaternion") {
00937
00938 result["e0"] = cosOover2 ;
00939 result["e1"] = sinOover2 * n1 ;
00940 result["e2"] = sinOover2 * n2;
00941 result["e3"] = sinOover2 * n3;
00942 } else if (type == "spin") {
00943
00944 result["Omega"] = 2.0*EMConsts::rad2deg * acos(cosOover2);
00945 result["n1"] = n1;
00946 result["n2"] = n2;
00947 result["n3"] = n3;
00948 } else if (type == "sgirot") {
00949
00950 result["q"] = 2.0*EMConsts::rad2deg * acos(cosOover2);
00951 result["n1"] = n1;
00952 result["n2"] = n2;
00953 result["n3"] = n3;
00954 } else if (type == "matrix") {
00955
00956 result["m11"] = x_mirror_scale*matrix[0][0]*inv_scale;
00957 result["m12"] = x_mirror_scale*matrix[0][1]*inv_scale;
00958 result["m13"] = x_mirror_scale*matrix[0][2]*inv_scale;
00959 result["m21"] = matrix[1][0]*inv_scale;
00960 result["m22"] = matrix[1][1]*inv_scale;
00961 result["m23"] = matrix[1][2]*inv_scale;
00962 result["m31"] = matrix[2][0]*inv_scale;
00963 result["m32"] = matrix[2][1]*inv_scale;
00964 result["m33"] = matrix[2][2]*inv_scale;
00965 } else {
00966 throw InvalidStringException(euler_type, "unknown Euler Type");
00967 }
00968
00969 return result;
00970 }
00971
00972 void Transform::set_trans(const float& x, const float& y, const float& z)
00973 {
00974 bool x_mirror = get_mirror();
00975
00976 if (x_mirror) matrix[0][3] = -x;
00977 else matrix[0][3] = x;
00978 matrix[1][3] = y;
00979 matrix[2][3] = z;
00980 }
00981
00982 Vec3f Transform::get_trans() const
00983 {
00984
00985 bool x_mirror = get_mirror();
00986 Vec3f v;
00987 if (x_mirror) v[0] = -matrix[0][3];
00988 else v[0] = matrix[0][3];
00989 v[1] = matrix[1][3];
00990 v[2] = matrix[2][3];
00991
00992 Util::apply_precision(v[0],ERR_LIMIT);
00993 Util::apply_precision(v[1],ERR_LIMIT);
00994 Util::apply_precision(v[2],ERR_LIMIT);
00995
00996 return v;
00997 }
00998
00999 void Transform::translate(const float& tx, const float& ty, const float& tz)
01000 {
01001 bool x_mirror = get_mirror();
01002 if (x_mirror) matrix[0][3] = -matrix[0][3] + tx;
01003 else matrix[0][3] = matrix[0][3] + tx;
01004 matrix[1][3] = matrix[1][3] + ty;
01005 matrix[2][3] = matrix[2][3] + tz;
01006 }
01007
01008 void Transform::translate_newBasis(const Transform& tcs, const float& tx, const float& ty, const float& tz)
01009 {
01010
01011 Transform tcsinv = Transform(tcs);
01012 tcsinv.set_trans(0.0, 0.0, 0.0);
01013 tcsinv.invert();
01014
01015
01016 Transform temp = Transform();
01017 temp.set_trans(tx, ty, tz);
01018 Transform nb_trans = tcsinv*temp;
01019
01020 translate(nb_trans.get_trans());
01021
01022 }
01023
01024 Vec2f Transform::get_trans_2d() const
01025 {
01026 bool x_mirror = get_mirror();
01027 Vec2f v;
01028 if (x_mirror) v[0] = -matrix[0][3];
01029 else v[0] = matrix[0][3];
01030 v[1] = matrix[1][3];
01031 return v;
01032 }
01033
01034
01035
01036 Vec3f Transform::get_pre_trans() const
01037 {
01038 Transform T(*this);
01039 T.set_trans(0,0,0);
01040 T.invert();
01041
01042 Transform soln = T*(*this);
01043
01044 return soln.get_trans();
01045 }
01046
01047 Vec2f Transform::get_pre_trans_2d() const
01048 {
01049 Transform T(*this);
01050 T.set_trans(0,0,0);
01051 T.invert();
01052
01053 Transform soln = T*(*this);
01054
01055 return soln.get_trans_2d();
01056 }
01057
01058
01059 void Transform::set_scale(const float& new_scale) {
01060 if (new_scale <= 0) {
01061 throw InvalidValueException(new_scale,"The scale factor in a Transform object must be positive and non zero");
01062 }
01063
01064
01065
01066 float old_scale = get_scale();
01067
01068 float n_scale = new_scale;
01069 Util::apply_precision(n_scale,ERR_LIMIT);
01070
01071 float corrected_scale = n_scale/old_scale;
01072 if ( corrected_scale != 1.0 ) {
01073 for(int i = 0; i < 3; ++i ) {
01074 for(int j = 0; j < 3; ++j ) {
01075 matrix[i][j] *= corrected_scale;
01076 }
01077 }
01078 }
01079 }
01080
01081 float Transform::get_scale() const {
01082 float determinant = get_determinant();
01083 if (determinant < 0 ) determinant *= -1;
01084
01085 float scale = std::pow(determinant,1.0f/3.0f);
01086 int int_scale = static_cast<int>(scale);
01087 float scale_residual = scale-static_cast<float>(int_scale);
01088 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01089
01090 Util::apply_precision(scale, ERR_LIMIT);
01091
01092 return scale;
01093 }
01094
01095 void Transform::scale(const float& scale)
01096 {
01097 float determinant = get_determinant();
01098 if (determinant < 0) determinant *= -1.0f;
01099 float newscale = std::pow(determinant,1.0f/3.0f) + scale;
01100
01101 if(newscale > 0.01)
01102 {
01103
01104 set_scale(newscale);
01105 }
01106 }
01107
01108 void print_matrix(gsl_matrix* M, unsigned int r, unsigned int c, const string& message ) {
01109 cout << "Message is " << message << endl;
01110 for ( unsigned int i = 0; i < r; ++i )
01111 {
01112 for ( unsigned int j = 0; j < c; ++j )
01113 {
01114 cout << gsl_matrix_get(M,i,j) << " ";
01115 }
01116 cout << endl;
01117 }
01118 }
01119
01120 void Transform::orthogonalize()
01121 {
01122 float scale;
01123 bool x_mirror;
01124 get_scale_and_mirror(scale,x_mirror);
01125 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
01126 double inv_scale = 1.0/static_cast<double>(scale);
01127 double mirror_scale = (x_mirror == true ? -1.0:1.0);
01128
01129 gsl_matrix * R = gsl_matrix_calloc(3,3);
01130 for ( unsigned int i = 0; i < 3; ++i )
01131 {
01132 for ( unsigned int j = 0; j < 3; ++j )
01133 {
01134 if (i == 0 && mirror_scale != 1.0 ) {
01135 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*mirror_scale*inv_scale );
01136 }
01137 else {
01138 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*inv_scale );
01139 }
01140 }
01141 }
01142
01143 gsl_matrix * V = gsl_matrix_calloc(3,3);
01144 gsl_vector * S = gsl_vector_calloc(3);
01145 gsl_vector * work = gsl_vector_calloc(3);
01146 gsl_linalg_SV_decomp (R, V, S, work);
01147
01148 gsl_matrix * Soln = gsl_matrix_calloc(3,3);
01149 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, R, V, 0.0, Soln);
01150
01151 for ( unsigned int i = 0; i < 3; ++i )
01152 {
01153 for ( unsigned int j = 0; j < 3; ++j )
01154 {
01155 matrix[i][j] = static_cast<float>( gsl_matrix_get(Soln,i,j) );
01156 }
01157 }
01158
01159
01160 if (scale != 1.0f) {
01161 for(int i=0; i<3; ++i) {
01162 for(int j=0; j<3; ++j) {
01163 matrix[i][j] *= scale;
01164 }
01165 }
01166 }
01167
01168
01169 if ( x_mirror ) {
01170 for(int j=0; j<3; ++j) {
01171 matrix[0][j] *= -1.0f;
01172 }
01173 }
01174
01175 gsl_matrix_free(V); gsl_matrix_free(R); gsl_matrix_free(Soln);
01176 gsl_vector_free(S); gsl_vector_free(work);
01177 }
01178
01179 void Transform::set_mirror(const bool x_mirror ) {
01180
01181 bool old_x_mirror = get_mirror();
01182 if (old_x_mirror == x_mirror) return;
01183 else {
01184
01185 for (int j = 0; j < 4; ++j ) {
01186 matrix[0][j] *= -1;
01187 }
01188 }
01189 }
01190
01191 bool Transform::get_mirror() const {
01192 float determinant = get_determinant();
01193
01194 bool x_mirror = false;
01195 if ( determinant < 0 ) x_mirror = true;
01196
01197 return x_mirror;
01198
01199 }
01200
01201 void Transform::get_scale_and_mirror(float& scale, bool& x_mirror) const {
01202
01203 float determinant = get_determinant();
01204 x_mirror = false;
01205 if ( determinant < 0 ) {
01206 x_mirror = true;
01207 determinant *= -1;
01208 }
01209 if (determinant != 1 ) {
01210 scale = std::pow(determinant,1.0f/3.0f);
01211 int int_scale = static_cast<int>(scale);
01212 float scale_residual = scale-static_cast<float>(int_scale);
01213 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01214 }
01215 else scale = 1;
01216
01217 Util::apply_precision(scale,ERR_LIMIT);
01218 }
01219
01220 float Transform::get_determinant() const
01221 {
01222 float det;
01223 double det2;
01224 det2 = matrix[0][0]*((double)matrix[1][1]*matrix[2][2]-(double)matrix[2][1]*matrix[1][2]);
01225 det2 -= matrix[0][1]*((double)matrix[1][0]*matrix[2][2]-(double)matrix[2][0]*matrix[1][2]);
01226 det2 += matrix[0][2]*((double)matrix[1][0]*matrix[2][1]-(double)matrix[2][0]*matrix[1][1]);
01227
01228 det = (float)det2;
01229 Util::apply_precision(det,ERR_LIMIT);
01230
01231 return det;
01232 }
01233
01234 void Transform::invert() {
01235
01236 double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
01237 double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
01238 double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
01239 double v0 = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
01240
01241 double cof00 = m11*m22-m12*m21;
01242 double cof11 = m22*m00-m20*m02;
01243 double cof22 = m00*m11-m01*m10;
01244 double cof01 = m10*m22-m20*m12;
01245 double cof02 = m10*m21-m20*m11;
01246 double cof12 = m00*m21-m01*m20;
01247 double cof10 = m01*m22-m02*m21;
01248 double cof20 = m01*m12-m02*m11;
01249 double cof21 = m00*m12-m10*m02;
01250
01251 double det = m00* cof00 + m02* cof02 -m01*cof01;
01252
01253 matrix[0][0] = (float)(cof00/det);
01254 matrix[0][1] = - (float)(cof10/det);
01255 matrix[0][2] = (float)(cof20/det);
01256 matrix[1][0] = - (float)(cof01/det);
01257 matrix[1][1] = (float)(cof11/det);
01258 matrix[1][2] = - (float)(cof21/det);
01259 matrix[2][0] = (float)(cof02/det);
01260 matrix[2][1] = - (float)(cof12/det);
01261 matrix[2][2] = (float)(cof22/det);
01262
01263 matrix[0][3] = (float)((- cof00*v0 + cof10*v1 - cof20*v2)/det);
01264 matrix[1][3] = (float)(( cof01*v0 - cof11*v1 + cof21*v2)/det);
01265 matrix[2][3] = (float)((- cof02*v0 + cof12*v1 - cof22*v2)/det);
01266 }
01267
01268 Transform Transform::inverse() const {
01269 Transform t(*this);
01270 t.invert();
01271 return t;
01272 }
01273
01274 void Transform::transpose_inplace() {
01275 float tempij;
01276 for (int i = 0; i < 3; i++) {
01277 for (int j = 0; j < i; j++) {
01278 if (i != j) {
01279 tempij= matrix[i][j];
01280 matrix[i][j] = matrix[j][i];
01281 matrix[j][i] = tempij;
01282 }
01283 }
01284 }
01285 }
01286
01287 Transform Transform::transpose() const {
01288 Transform t(*this);
01289 t.transpose_inplace();
01290 return t;
01291 }
01292
01293
01294 Transform EMAN::operator*(const Transform & M2, const Transform & M1)
01295 {
01296 Transform result;
01297 for (int i=0; i<3; i++) {
01298 for (int j=0; j<4; j++) {
01299 result[i][j] = M2[i][0] * M1[0][j] + M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
01300 }
01301 result[i][3] += M2[i][3];
01302 }
01303
01304 return result;
01305 }
01306
01307 void Transform::assert_valid_2d() const {
01308 int rotation_error = 0;
01309 int translation_error = 0;
01310 if (fabs(matrix[2][0]) > ERR_LIMIT) rotation_error++;
01311 if (fabs(matrix[2][1]) > ERR_LIMIT) rotation_error++;
01312 if (fabs(matrix[2][3]) > ERR_LIMIT) translation_error++;
01313 if (fabs(matrix[0][2]) > ERR_LIMIT) rotation_error++;
01314 if (fabs(matrix[1][2]) > ERR_LIMIT) rotation_error++;
01315
01316 if ( translation_error && rotation_error ) {
01317 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and 3D translations. This object can not be considered 2D");
01318 } else if ( translation_error ) {
01319 throw UnexpectedBehaviorException("Error, the internal matrix contains a non zero z component for a 3D translation. This object can not be considered 2D");
01320 }
01321 else if ( rotation_error ) {
01322 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and this object can not be considered 2D");
01323 }
01324
01325 }
01326
01327
01328
01329 Transform Transform::get_sym(const string & sym_name, int n) const
01330 {
01331 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01332 Transform ret;
01333 ret = (*this) * sym->get_sym(n);
01334 delete sym;
01335 return ret;
01336 }
01337
01338 vector<Transform > Transform::get_sym_proj(const string & sym_name) const
01339 {
01340 vector<Transform> ret;
01341 Transform t;
01342 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01343 int nsym = sym->get_nsym();
01344 int n = nsym;
01345
01346 if ((sym_name[0] == 'c' || sym_name[0] == 'd' ) && fabs(matrix[2][2]) < 1.e-6){
01347
01348 Dict d1,d2;
01349 d2["theta"] = (double)90.0;
01350 d2["psi"] = (double)0.0;
01351 d2["phi"] = (double)0.0;
01352 d2["type"] = "spider";
01353 d1 = this->get_rotation("spider");
01354
01355 if (sym_name[0] == 'c') {
01356 if( nsym%2 == 0) n = nsym/2;
01357
01358 for (int k=0;k<n;k++) {
01359 d2["phi"] = (double)d1["phi"] + k*double(360.0)/ nsym;
01360 d2["psi"] = d1["psi"];
01361 t.set_rotation(d2);
01362 ret.push_back( t );
01363 }
01364
01365 }
01366 else {
01367 nsym = nsym/2;
01368
01369 if (nsym%2 == 0) {
01370 n = nsym;
01371 float cos_phi = cos( EMConsts::deg2rad*360.0/2/nsym );
01372
01373 for (int k=0;k<n;k++){
01374
01375 if(k%2==0) {
01376
01377 d2["phi"] = (double)d1["phi"] + k/2*double(360.0)/ nsym;
01378 d2["psi"] = d1["psi"];
01379 t.set_rotation(d2);
01380 ret.push_back( t );
01381 }
01382 else {
01383
01384 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01385
01386 d2["phi"] = k/2*double(360.0)/ nsym +180 - (double)d1["phi"];
01387 d2["psi"] = (double)d1["psi"] + 180;
01388 t.set_rotation(d2);
01389 ret.push_back( t );
01390 }
01391 }
01392
01393 }
01394 }
01395
01396
01397
01398 else {
01399 n = nsym*2;
01400 float cos_phi = cos( EMConsts::deg2rad*360.0/4/nsym );
01401 for (int k=0;k<n;k++){
01402
01403 if(k%4==0) {
01404
01405 d2["phi"] = (double)d1["phi"] + k/4*360.0/ nsym;
01406 d2["psi"] = (double)d1["psi"];
01407 t.set_rotation(d2);
01408 ret.push_back( t );
01409 }
01410 else if( k%4 ==1) {
01411 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01412
01413 d2["phi"] = k/4*360.0/nsym + 360.0/2/nsym+180 - (double)d1["phi"];
01414 d2["psi"] = (double)d1["psi"] + 180;
01415 t.set_rotation(d2);
01416 ret.push_back( t );
01417 }
01418
01419 }
01420
01421 else if( k%4 ==2) {
01422
01423 d2["phi"] = k/4*360.0/ nsym+360.0/2/nsym+180 + (double)d1["phi"];
01424 d2["psi"] = (double)d1["psi"];
01425 t.set_rotation(d2);
01426 ret.push_back( t );
01427
01428 }
01429
01430 else if( k%4 ==3) {
01431 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ) {
01432 d2["phi"] = k/4*360.0/nsym+ 2.0*360.0/2/nsym - (double)d1["phi"];
01433 d2["psi"] = (double)d1["psi"] + 180;
01434 t.set_rotation(d2);
01435 ret.push_back( t );
01436 }
01437 }
01438
01439 }
01440 }
01441
01442 }
01443
01444 }
01445 else {
01446 for (int k=0;k<nsym;k++) {
01447 t = sym->get_sym(k);
01448 ret.push_back( (*this) * t );
01449 }
01450
01451 }
01452 delete sym;
01453 return ret;
01454 }
01455
01456
01457 int Transform::get_nsym(const string & sym_name)
01458 {
01459 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01460 int nsym = sym->get_nsym();
01461 delete sym;
01462 return nsym;
01463 }
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01483
01484
01485
01486
01487
01488
01489
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01521
01522
01523
01524
01525
01526
01527
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01552
01553
01554
01555
01556
01557
01558
01559
01560
01561
01562
01563
01564
01565
01566
01567
01568
01569
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01591
01592
01593
01594
01595
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01611
01614
01615
01616
01617
01618
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01647
01648
01649
01650
01651
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730
01731
01732
01733
01734
01735
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01781
01782
01783
01784
01785
01786
01787
01788
01789
01791
01792
01793
01794
01795
01796
01797
01798
01799
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815
01816
01817
01818
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
01983
01984
01985
01986
01987
01988
01989
01990
01991
01992
01993
01994
01995
01996
01997
01998
01999
02000
02001
02002
02003
02004
02005
02006
02007
02008
02009
02010
02011
02012
02013
02014
02015
02016
02017
02018
02019
02020
02021
02022
02023
02024
02025
02026
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067
02068
02069
02070
02071
02072
02073
02074
02075
02076
02077
02078
02079
02080
02081
02082
02083
02084
02085
02086
02087
02088
02089
02090
02091
02092
02093
02094
02095
02096
02097
02098
02099
02100
02101
02102
02103
02104
02105
02106
02107
02108
02109
02110
02111
02112
02113
02114
02115
02116
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02129
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02147
02148
02149
02150
02151
02152
02153
02155
02157
02158
02160
02161
02162
02163
02164
02165
02166
02167
02168
02169
02170
02171
02172
02173
02174
02175
02176
02177
02178
02179
02180
02181
02182
02183
02184
02185
02186
02187
02188
02189
02190
02191
02192
02193
02194
02195
02196
02197
02198
02199
02200
02201
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02235
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245
02246
02247
02248
02249
02250
02251
02252
02253
02254
02255
02256
02257
02258
02259
02260
02261
02262
02263
02264
02265
02266
02267
02268
02269
02270
02271
02272
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301
02302
02303
02304
02305
02306
02307
02308
02309
02310
02311
02312
02313
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341
02342
02343
02344
02345
02346
02347
02348
02355
02356
02357
02358
02359
02360
02361
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384
02385
02386
02387
02388
02389
02390
02391
02392
02393
02394
02395
02396
02397
02398
02399
02400
02401
02402
02403
02404
02405
02406
02407
02408
02409
02410
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02424
02425
02426
02427
02428
02431
02432
02433
02434
02435
02436
02437
02438
02439
02440
02441
02442
02443
02444
02445
02446
02447
02448
02449
02450
02451
02452
02453
02454
02455
02456
02457
02458
02459
02460
02461
02462
02463
02464
02465
02466
02467
02468
02469
02470
02471
02472
02473
02474
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509
02510
02511
02512
02513
02514
02515
02516
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541
02542
02543
02544
02545
02546
02547
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566
02567
02568
02569
02570
02571
02572
02573
02574
02575
02576
02577
02578
02579
02580
02581
02582
02583
02584
02585
02586
02587
02588
02589
02590
02591
02592
02593
02594
02595
02596
02597
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632
02633
02634
02635
02636
02637
02638
02639
02640
02641
02642