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
00785 Transform ret(*this);
00786 ret.set_rotation(rot);
00787
00788 Vec3f trans = get_trans();
00789 trans[0] = -trans[0];
00790 ret.set_trans(trans);
00791
00792
00793
00794 return ret;
00795 }
00796
00797 Transform Transform::get_vflip_transform() const {
00798
00799 Dict rot = get_rotation("eman");
00800 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
00801 rot["phi"] = - static_cast<float>(rot["phi"]);
00802
00803
00804
00805 Transform ret(*this);
00806 ret.set_rotation(rot);
00807
00808 Vec3f trans = get_trans();
00809 trans[1] = -trans[1];
00810 ret.set_trans(trans);
00811
00812 return ret;
00813 }
00814
00815 Dict Transform::get_rotation(const string& euler_type) const
00816 {
00817 Dict result;
00818
00819
00820 float scale;
00821 bool x_mirror;
00822 get_scale_and_mirror(scale,x_mirror);
00823 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
00824
00825 double cosalt = matrix[2][2]/scale;
00826 double x_mirror_scale = (x_mirror ? -1.0f : 1.0f);
00827 double inv_scale = 1.0f/scale;
00828
00829 double az = 0;
00830 double alt = 0;
00831 double phi = 0;
00832 double phiS = 0;
00833 double psiS = 0;
00834
00835
00836
00837 if (cosalt >= 1) {
00838 alt = 0;
00839 az = 0;
00840 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00841 } else if (cosalt <= -1) {
00842 alt = 180;
00843 az = 0;
00844 phi = (double)EMConsts::rad2deg * atan2(-x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00845 } else {
00846
00847 az = (double)EMConsts::rad2deg * atan2(scale*matrix[2][0], -scale*matrix[2][1]);
00848
00849 if (matrix[2][2]==0.0)
00850 alt = 90.0;
00851 else
00852 alt = (double)EMConsts::rad2deg * atan(sqrt((double)matrix[2][0]*matrix[2][0]+(double)matrix[2][1]*matrix[2][1])/fabs(matrix[2][2]));
00853
00854 if (matrix[2][2] * scale < 0)
00855 alt = 180.0f-alt;
00856
00857 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*(double)matrix[0][2], (double)matrix[1][2]);
00858
00859 }
00860
00861 phi = phi-360.0*floor(phi/360.0);
00862 az = az -360.0*floor(az/360.0);
00863
00864
00865 if (cosalt >= 1) {
00866 phiS = 0;
00867 psiS = phi;
00868 } else if (cosalt <= -1) {
00869 phiS = 0;
00870 psiS = phi + 180.0;
00871 } else {
00872 phiS = az - 90.0;
00873 psiS = phi + 90.0;
00874 }
00875
00876 phiS = phiS-360.0*floor(phiS/360.0);
00877 psiS = psiS-360.0*floor(psiS/360.0);
00878
00879
00880 double xtilt = 0;
00881 double ytilt = 0;
00882 double ztilt = 0;
00883
00884
00885 string type = Util::str_to_lower(euler_type);
00886
00887 result["type"] = type;
00888 if (type == "2d") {
00889 assert_valid_2d();
00890 result["alpha"] = phi;
00891 } else if (type == "eman") {
00892
00893 result["az"] = az;
00894 result["alt"] = alt;
00895 result["phi"] = phi;
00896 } else if (type == "imagic") {
00897
00898 result["alpha"] = az;
00899 result["beta"] = alt;
00900 result["gamma"] = phi;
00901 } else if (type == "spider") {
00902
00903 result["phi"] = phiS;
00904 result["theta"] = alt;
00905 result["psi"] = psiS;
00906 } else if (type == "mrc") {
00907
00908 result["phi"] = phiS;
00909 result["theta"] = alt;
00910 result["omega"] = psiS;
00911 } else if (type == "xyz") {
00912
00913 xtilt = atan2(-sin(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt),cos(EMConsts::deg2rad*alt));
00914 ytilt = asin( cos(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt));
00915 ztilt = psiS*EMConsts::deg2rad - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
00916
00917 xtilt *= EMConsts::rad2deg; ytilt *= EMConsts::rad2deg; ztilt *= EMConsts::rad2deg;
00918 xtilt = xtilt-360*.0*floor((xtilt+180.0)/360.0);
00919 ytilt = ytilt-360*.0*floor((ytilt+180.0)/360.0);
00920 ztilt = ztilt-360*.0*floor((ztilt+180.0)/360.0);
00921
00922 result["xtilt"] = xtilt;
00923 result["ytilt"] = ytilt;
00924 result["ztilt"] = ztilt;
00925 } else if ((type == "quaternion") || (type == "spin") || (type == "sgirot")) {
00926
00927
00928
00929
00930
00931 double traceR = matrix[0][0]+matrix[1][1]+matrix[2][2];
00932 double cosomega = (traceR-1.0)/2.0;
00933 if (cosomega>1.0) cosomega=1.0;
00934 if (cosomega<-1.0) cosomega=-1.0;
00935
00936
00937
00938 double sinOover2= sqrt((1.0 -cosomega)/2.0);
00939 double cosOover2= sqrt(1.0 -sinOover2*sinOover2);
00940 double sinomega = 2* sinOover2*cosOover2;
00941 double n1 = 0; double n2 = 0; double n3 = 0;
00942 if (sinomega>0) {
00943 n1 = (matrix[1][2]-matrix[2][1])/2.0/sinomega ;
00944 n2 = (matrix[2][0]-matrix[0][2])/2.0/sinomega ;
00945 n3 = (matrix[0][1]-matrix[1][0])/2.0/sinomega ;
00946 }
00947
00948
00949
00950 if (type == "quaternion"){
00951 result["e0"] = cosOover2 ;
00952 result["e1"] = sinOover2 * n1 ;
00953 result["e2"] = sinOover2 * n2;
00954 result["e3"] = sinOover2 * n3;
00955 }
00956
00957 if (type == "spin"){
00958 result["omega"] = EMConsts::rad2deg * acos(cosomega);
00959 result["n1"] = n1;
00960 result["n2"] = n2;
00961 result["n3"] = n3;
00962 }
00963
00964 if (type == "sgirot"){
00965 result["q"] = EMConsts::rad2deg * acos(cosomega);
00966 result["n1"] = n1;
00967 result["n2"] = n2;
00968 result["n3"] = n3;
00969 }
00970
00971 } else if (type == "matrix") {
00972
00973 result["m11"] = x_mirror_scale*matrix[0][0]*inv_scale;
00974 result["m12"] = x_mirror_scale*matrix[0][1]*inv_scale;
00975 result["m13"] = x_mirror_scale*matrix[0][2]*inv_scale;
00976 result["m21"] = matrix[1][0]*inv_scale;
00977 result["m22"] = matrix[1][1]*inv_scale;
00978 result["m23"] = matrix[1][2]*inv_scale;
00979 result["m31"] = matrix[2][0]*inv_scale;
00980 result["m32"] = matrix[2][1]*inv_scale;
00981 result["m33"] = matrix[2][2]*inv_scale;
00982 } else {
00983 throw InvalidStringException(euler_type, "unknown Euler Type");
00984 }
00985
00986 return result;
00987 }
00988
00989 void Transform::set_trans(const float& x, const float& y, const float& z)
00990 {
00991 bool x_mirror = get_mirror();
00992
00993 if (x_mirror) matrix[0][3] = -x;
00994 else matrix[0][3] = x;
00995 matrix[1][3] = y;
00996 matrix[2][3] = z;
00997 }
00998
00999 Vec3f Transform::get_trans() const
01000 {
01001
01002 bool x_mirror = get_mirror();
01003 Vec3f v;
01004 if (x_mirror) v[0] = -matrix[0][3];
01005 else v[0] = matrix[0][3];
01006 v[1] = matrix[1][3];
01007 v[2] = matrix[2][3];
01008
01009 Util::apply_precision(v[0],ERR_LIMIT);
01010 Util::apply_precision(v[1],ERR_LIMIT);
01011 Util::apply_precision(v[2],ERR_LIMIT);
01012
01013 return v;
01014 }
01015
01016 void Transform::translate(const float& tx, const float& ty, const float& tz)
01017 {
01018 bool x_mirror = get_mirror();
01019 if (x_mirror) matrix[0][3] = -matrix[0][3] + tx;
01020 else matrix[0][3] = matrix[0][3] + tx;
01021 matrix[1][3] = matrix[1][3] + ty;
01022 matrix[2][3] = matrix[2][3] + tz;
01023 }
01024
01025 void Transform::translate_newBasis(const Transform& tcs, const float& tx, const float& ty, const float& tz)
01026 {
01027
01028 Transform tcsinv = Transform(tcs);
01029 tcsinv.set_trans(0.0, 0.0, 0.0);
01030 tcsinv.invert();
01031
01032
01033 Transform temp = Transform();
01034 temp.set_trans(tx, ty, tz);
01035 Transform nb_trans = tcsinv*temp;
01036
01037 translate(nb_trans.get_trans());
01038
01039 }
01040
01041 Vec2f Transform::get_trans_2d() const
01042 {
01043 bool x_mirror = get_mirror();
01044 Vec2f v;
01045 if (x_mirror) v[0] = -matrix[0][3];
01046 else v[0] = matrix[0][3];
01047 v[1] = matrix[1][3];
01048 return v;
01049 }
01050
01051
01052
01053 Vec3f Transform::get_pre_trans() const
01054 {
01055 Transform T(*this);
01056 T.set_trans(0,0,0);
01057 T.invert();
01058
01059 Transform soln = T*(*this);
01060
01061 return soln.get_trans();
01062 }
01063
01064 Vec2f Transform::get_pre_trans_2d() const
01065 {
01066 Transform T(*this);
01067 T.set_trans(0,0,0);
01068 T.invert();
01069
01070 Transform soln = T*(*this);
01071
01072 return soln.get_trans_2d();
01073 }
01074
01075
01076 void Transform::set_scale(const float& new_scale) {
01077 if (new_scale <= 0) {
01078 throw InvalidValueException(new_scale,"The scale factor in a Transform object must be positive and non zero");
01079 }
01080
01081
01082
01083 float old_scale = get_scale();
01084
01085 float n_scale = new_scale;
01086 Util::apply_precision(n_scale,ERR_LIMIT);
01087
01088 float corrected_scale = n_scale/old_scale;
01089 if ( corrected_scale != 1.0 ) {
01090 for(int i = 0; i < 3; ++i ) {
01091 for(int j = 0; j < 3; ++j ) {
01092 matrix[i][j] *= corrected_scale;
01093 }
01094 }
01095 }
01096 }
01097
01098 float Transform::get_scale() const {
01099 float determinant = get_determinant();
01100 if (determinant < 0 ) determinant *= -1;
01101
01102 float scale = std::pow(determinant,1.0f/3.0f);
01103 int int_scale = static_cast<int>(scale);
01104 float scale_residual = scale-static_cast<float>(int_scale);
01105 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01106
01107 Util::apply_precision(scale, ERR_LIMIT);
01108
01109 return scale;
01110 }
01111
01112 void Transform::scale(const float& scale)
01113 {
01114 float determinant = get_determinant();
01115 if (determinant < 0) determinant *= -1.0f;
01116 float newscale = std::pow(determinant,1.0f/3.0f) + scale;
01117 if(newscale > 0.0001) set_scale(newscale);
01118 }
01119
01120 void print_matrix(gsl_matrix* M, unsigned int r, unsigned int c, const string& message ) {
01121 cout << "Message is " << message << endl;
01122 for ( unsigned int i = 0; i < r; ++i )
01123 {
01124 for ( unsigned int j = 0; j < c; ++j )
01125 {
01126 cout << gsl_matrix_get(M,i,j) << " ";
01127 }
01128 cout << endl;
01129 }
01130 }
01131
01132 void Transform::orthogonalize()
01133 {
01134 float scale;
01135 bool x_mirror;
01136 get_scale_and_mirror(scale,x_mirror);
01137 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
01138 double inv_scale = 1.0/static_cast<double>(scale);
01139 double mirror_scale = (x_mirror == true ? -1.0:1.0);
01140
01141 gsl_matrix * R = gsl_matrix_calloc(3,3);
01142 for ( unsigned int i = 0; i < 3; ++i )
01143 {
01144 for ( unsigned int j = 0; j < 3; ++j )
01145 {
01146 if (i == 0 && mirror_scale != 1.0 ) {
01147 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*mirror_scale*inv_scale );
01148 }
01149 else {
01150 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*inv_scale );
01151 }
01152 }
01153 }
01154
01155 gsl_matrix * V = gsl_matrix_calloc(3,3);
01156 gsl_vector * S = gsl_vector_calloc(3);
01157 gsl_vector * work = gsl_vector_calloc(3);
01158 gsl_linalg_SV_decomp (R, V, S, work);
01159
01160 gsl_matrix * Soln = gsl_matrix_calloc(3,3);
01161 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, R, V, 0.0, Soln);
01162
01163 for ( unsigned int i = 0; i < 3; ++i )
01164 {
01165 for ( unsigned int j = 0; j < 3; ++j )
01166 {
01167 matrix[i][j] = static_cast<float>( gsl_matrix_get(Soln,i,j) );
01168 }
01169 }
01170
01171
01172 if (scale != 1.0f) {
01173 for(int i=0; i<3; ++i) {
01174 for(int j=0; j<3; ++j) {
01175 matrix[i][j] *= scale;
01176 }
01177 }
01178 }
01179
01180
01181 if ( x_mirror ) {
01182 for(int j=0; j<3; ++j) {
01183 matrix[0][j] *= -1.0f;
01184 }
01185 }
01186
01187 gsl_matrix_free(V); gsl_matrix_free(R); gsl_matrix_free(Soln);
01188 gsl_vector_free(S); gsl_vector_free(work);
01189 }
01190
01191 void Transform::set_mirror(const bool x_mirror ) {
01192
01193 bool old_x_mirror = get_mirror();
01194 if (old_x_mirror == x_mirror) return;
01195 else {
01196
01197 for (int j = 0; j < 4; ++j ) {
01198 matrix[0][j] *= -1;
01199 }
01200 }
01201 }
01202
01203 bool Transform::get_mirror() const {
01204 float determinant = get_determinant();
01205
01206 bool x_mirror = false;
01207 if ( determinant < 0 ) x_mirror = true;
01208
01209 return x_mirror;
01210
01211 }
01212
01213 void Transform::get_scale_and_mirror(float& scale, bool& x_mirror) const {
01214
01215 float determinant = get_determinant();
01216 x_mirror = false;
01217 if ( determinant < 0 ) {
01218 x_mirror = true;
01219 determinant *= -1;
01220 }
01221 if (determinant != 1 ) {
01222 scale = std::pow(determinant,1.0f/3.0f);
01223 int int_scale = static_cast<int>(scale);
01224 float scale_residual = scale-static_cast<float>(int_scale);
01225 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01226 }
01227 else scale = 1;
01228
01229 Util::apply_precision(scale,ERR_LIMIT);
01230 }
01231
01232 float Transform::get_determinant() const
01233 {
01234 float det;
01235 double det2;
01236 det2 = matrix[0][0]*((double)matrix[1][1]*matrix[2][2]-(double)matrix[2][1]*matrix[1][2]);
01237 det2 -= matrix[0][1]*((double)matrix[1][0]*matrix[2][2]-(double)matrix[2][0]*matrix[1][2]);
01238 det2 += matrix[0][2]*((double)matrix[1][0]*matrix[2][1]-(double)matrix[2][0]*matrix[1][1]);
01239
01240 det = (float)det2;
01241 Util::apply_precision(det,ERR_LIMIT);
01242
01243 return det;
01244 }
01245
01246 void Transform::invert() {
01247
01248 double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
01249 double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
01250 double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
01251 double v0 = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
01252
01253 double cof00 = m11*m22-m12*m21;
01254 double cof11 = m22*m00-m20*m02;
01255 double cof22 = m00*m11-m01*m10;
01256 double cof01 = m10*m22-m20*m12;
01257 double cof02 = m10*m21-m20*m11;
01258 double cof12 = m00*m21-m01*m20;
01259 double cof10 = m01*m22-m02*m21;
01260 double cof20 = m01*m12-m02*m11;
01261 double cof21 = m00*m12-m10*m02;
01262
01263 double det = m00* cof00 + m02* cof02 -m01*cof01;
01264
01265 matrix[0][0] = (float)(cof00/det);
01266 matrix[0][1] = - (float)(cof10/det);
01267 matrix[0][2] = (float)(cof20/det);
01268 matrix[1][0] = - (float)(cof01/det);
01269 matrix[1][1] = (float)(cof11/det);
01270 matrix[1][2] = - (float)(cof21/det);
01271 matrix[2][0] = (float)(cof02/det);
01272 matrix[2][1] = - (float)(cof12/det);
01273 matrix[2][2] = (float)(cof22/det);
01274
01275 matrix[0][3] = (float)((- cof00*v0 + cof10*v1 - cof20*v2)/det);
01276 matrix[1][3] = (float)(( cof01*v0 - cof11*v1 + cof21*v2)/det);
01277 matrix[2][3] = (float)((- cof02*v0 + cof12*v1 - cof22*v2)/det);
01278 }
01279
01280 Transform Transform::inverse() const {
01281 Transform t(*this);
01282 t.invert();
01283 return t;
01284 }
01285
01286 void Transform::transpose_inplace() {
01287 float tempij;
01288 for (int i = 0; i < 3; i++) {
01289 for (int j = 0; j < i; j++) {
01290 if (i != j) {
01291 tempij= matrix[i][j];
01292 matrix[i][j] = matrix[j][i];
01293 matrix[j][i] = tempij;
01294 }
01295 }
01296 }
01297 }
01298
01299 Transform Transform::transpose() const {
01300 Transform t(*this);
01301 t.transpose_inplace();
01302 return t;
01303 }
01304
01305
01306 Transform EMAN::operator*(const Transform & M2, const Transform & M1)
01307 {
01308 Transform result;
01309 for (int i=0; i<3; i++) {
01310 for (int j=0; j<4; j++) {
01311 result[i][j] = M2[i][0] * M1[0][j] + M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
01312 }
01313 result[i][3] += M2[i][3];
01314 }
01315
01316 return result;
01317 }
01318
01319 void Transform::assert_valid_2d() const {
01320 int rotation_error = 0;
01321 int translation_error = 0;
01322 if (fabs(matrix[2][0]) > ERR_LIMIT) rotation_error++;
01323 if (fabs(matrix[2][1]) > ERR_LIMIT) rotation_error++;
01324 if (fabs(matrix[2][3]) > ERR_LIMIT) translation_error++;
01325 if (fabs(matrix[0][2]) > ERR_LIMIT) rotation_error++;
01326 if (fabs(matrix[1][2]) > ERR_LIMIT) rotation_error++;
01327
01328 if ( translation_error && rotation_error ) {
01329 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and 3D translations. This object can not be considered 2D");
01330 } else if ( translation_error ) {
01331 throw UnexpectedBehaviorException("Error, the internal matrix contains a non zero z component for a 3D translation. This object can not be considered 2D");
01332 }
01333 else if ( rotation_error ) {
01334 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and this object can not be considered 2D");
01335 }
01336
01337 }
01338
01339
01340
01341 Transform Transform::get_sym(const string & sym_name, int n) const
01342 {
01343 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01344 Transform ret;
01345 ret = (*this) * sym->get_sym(n);
01346 delete sym;
01347 return ret;
01348 }
01349
01350 vector<Transform > Transform::get_sym_proj(const string & sym_name) const
01351 {
01352 vector<Transform> ret;
01353 Transform t;
01354 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01355 int nsym = sym->get_nsym();
01356 int n = nsym;
01357
01358 if ((sym_name[0] == 'c' || sym_name[0] == 'd' ) && fabs(matrix[2][2]) < 1.e-6){
01359
01360 Dict d1,d2;
01361 d2["theta"] = (double)90.0;
01362 d2["psi"] = (double)0.0;
01363 d2["phi"] = (double)0.0;
01364 d2["type"] = "spider";
01365 d1 = this->get_rotation("spider");
01366
01367 if (sym_name[0] == 'c') {
01368 if( nsym%2 == 0) n = nsym/2;
01369
01370 for (int k=0;k<n;k++) {
01371 d2["phi"] = (double)d1["phi"] + k*double(360.0)/ nsym;
01372 d2["psi"] = d1["psi"];
01373 t.set_rotation(d2);
01374 ret.push_back( t );
01375 }
01376
01377 }
01378 else {
01379 nsym = nsym/2;
01380
01381 if (nsym%2 == 0) {
01382 n = nsym;
01383 float cos_phi = cos( EMConsts::deg2rad*360.0/2/nsym );
01384
01385 for (int k=0;k<n;k++){
01386
01387 if(k%2==0) {
01388
01389 d2["phi"] = (double)d1["phi"] + k/2*double(360.0)/ nsym;
01390 d2["psi"] = d1["psi"];
01391 t.set_rotation(d2);
01392 ret.push_back( t );
01393 }
01394 else {
01395
01396 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01397
01398 d2["phi"] = k/2*double(360.0)/ nsym +180 - (double)d1["phi"];
01399 d2["psi"] = (double)d1["psi"] + 180;
01400 t.set_rotation(d2);
01401 ret.push_back( t );
01402 }
01403 }
01404
01405 }
01406 }
01407
01408
01409
01410 else {
01411 n = nsym*2;
01412 float cos_phi = cos( EMConsts::deg2rad*360.0/4/nsym );
01413 for (int k=0;k<n;k++){
01414
01415 if(k%4==0) {
01416
01417 d2["phi"] = (double)d1["phi"] + k/4*360.0/ nsym;
01418 d2["psi"] = (double)d1["psi"];
01419 t.set_rotation(d2);
01420 ret.push_back( t );
01421 }
01422 else if( k%4 ==1) {
01423 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01424
01425 d2["phi"] = k/4*360.0/nsym + 360.0/2/nsym+180 - (double)d1["phi"];
01426 d2["psi"] = (double)d1["psi"] + 180;
01427 t.set_rotation(d2);
01428 ret.push_back( t );
01429 }
01430
01431 }
01432
01433 else if( k%4 ==2) {
01434
01435 d2["phi"] = k/4*360.0/ nsym+360.0/2/nsym+180 + (double)d1["phi"];
01436 d2["psi"] = (double)d1["psi"];
01437 t.set_rotation(d2);
01438 ret.push_back( t );
01439
01440 }
01441
01442 else if( k%4 ==3) {
01443 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ) {
01444 d2["phi"] = k/4*360.0/nsym+ 2.0*360.0/2/nsym - (double)d1["phi"];
01445 d2["psi"] = (double)d1["psi"] + 180;
01446 t.set_rotation(d2);
01447 ret.push_back( t );
01448 }
01449 }
01450
01451 }
01452 }
01453
01454 }
01455
01456 }
01457 else {
01458 for (int k=0;k<nsym;k++) {
01459 t = sym->get_sym(k);
01460 ret.push_back( (*this) * t );
01461 }
01462
01463 }
01464 delete sym;
01465 return ret;
01466 }
01467
01468
01469 int Transform::get_nsym(const string & sym_name)
01470 {
01471 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01472 int nsym = sym->get_nsym();
01473 delete sym;
01474 return nsym;
01475 }
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01495
01496
01497
01498
01499
01500
01501
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01533
01534
01535
01536
01537
01538
01539
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550
01551
01552
01553
01554
01555
01556
01557
01558
01559
01560
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
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01603
01604
01605
01606
01607
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01626
01627
01628
01629
01630
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01659
01660
01661
01662
01663
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
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
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
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
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01793
01794
01795
01796
01797
01798
01799
01800
01801
01803
01804
01805
01806
01807
01808
01809
01810
01811
01814
01815
01816
01817
01818
01819
01820
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01912
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924
01925
01926
01927
01928
01929
01930
01931
01932
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
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
02117
02118
02119
02120
02121
02122
02123
02124
02125
02126
02127
02128
02130
02131
02132
02133
02134
02135
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02159
02160
02161
02162
02163
02164
02165
02167
02169
02170
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
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02235
02236
02237
02238
02239
02240
02241
02242
02243
02244
02245
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
02349
02350
02351
02352
02353
02354
02355
02356
02357
02358
02359
02360
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
02423
02424
02425
02426
02427
02428
02429
02430
02431
02432
02433
02434
02436
02437
02438
02439
02440
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
02643
02644
02645
02646
02647
02648
02649
02650
02651
02652
02653
02654