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.0f;
00072 d["az"] = 270.0f;
00073 d["alt"] = 58.282523f;
00079 t.set_rotation(d);
00080 return t;
00081 }
00082
00083 Transform Transform::tet_3_to_2() {
00084 Transform t;
00085 Dict d;
00086 d["type"] = "eman";
00087 d["phi"] = 45.0f;
00088 d["az"] = 0.0f;
00089 d["alt"] = 54.73561f;
00093 t.set_rotation(d);
00094 return t;
00095 }
00096
00097
00098 Transform::Transform()
00099 {
00100 to_identity();
00101 }
00102
00103 Transform::Transform( const Transform& that )
00104 {
00105 *this = that;
00106 }
00107
00108 Transform& Transform::operator=(const Transform& that ) {
00109 if (this != &that ) {
00110 memcpy(matrix,that.matrix,12*sizeof(float));
00111
00112 }
00113 return *this;
00114 }
00115
00116 bool Transform::operator==(const Transform& rhs) const{
00117 if (memcmp(this->matrix, rhs.matrix, 3*4*sizeof(float)) == 0) {
00118 return true;
00119 }
00120 else {
00121 return false;
00122 }
00123 }
00124
00125 bool Transform::operator!=(const Transform& rhs) const{
00126 return !(operator==(rhs));
00127 }
00128
00129 Transform::Transform(const Dict& d) {
00130 to_identity();
00131 set_params(d);
00132 }
00133
00134
00135 Transform::Transform(const float array[12]) {
00136 memcpy(matrix,array,12*sizeof(float));
00137 }
00138
00139 Transform::Transform(const vector<float> array)
00140 {
00141 set_matrix(array);
00142 }
00143
00144 void Transform::set_matrix(const vector<float>& v)
00145 {
00146 if (v.size() != 12 ) throw InvalidParameterException("The construction array must be of size 12");
00147
00148 for(int i=0; i<3; ++i) {
00149 for(int j=0; j<4; ++j) {
00150 matrix[i][j] = v[i*4+j];
00151 }
00152 }
00153 }
00154
00155 void Transform::copy_matrix_into_array(float* const array) const {
00156
00157 int idx = 0;
00158 for(int i=0; i<3; ++i) {
00159 for(int j=0; j<4; ++j) {
00160 array[idx] = matrix[i][j];
00161 idx ++;
00162 }
00163 }
00164 }
00165
00166 vector<float> Transform::get_matrix() const
00167 {
00168 vector<float> ret(12);
00169 for(int i=0; i<3; ++i) {
00170 for(int j=0; j<4; ++j) {
00171 ret[i*4+j] = matrix[i][j];
00172 }
00173 }
00174 return ret;
00175 }
00176
00177 vector<float> Transform::get_matrix_4x4() const
00178 {
00179 vector<float> ret(16);
00180 for(int i=0; i<3; ++i) {
00181 for(int j=0; j<4; ++j) {
00182 ret[i*4+j] = matrix[i][j];
00183 }
00184 }
00185 ret[12] = 0.0;
00186 ret[13] = 0.0;
00187 ret[14] = 0.0;
00188 ret[15] = 1.0;
00189
00190 return ret;
00191 }
00192 void Transform::to_identity()
00193 {
00194
00195 for(int i=0; i<3; ++i) {
00196 for(int j=0; j<4; ++j) {
00197 if(i==j) {
00198 matrix[i][j] = 1;
00199 }
00200 else {
00201 matrix[i][j] = 0;
00202 }
00203 }
00204 }
00205 }
00206
00207 bool Transform::is_identity() const {
00208 for(int i=0; i<3; ++i) {
00209 for(int j=0; j<4; ++j) {
00210 float c = matrix[i][j];
00211 Util::apply_precision(c,ERR_LIMIT);
00212 if(i==j) {
00213 if (c != 1.0) return false;
00214 }
00215 else {
00216 if (c != 0.0) return false;
00217 }
00218 }
00219 }
00220 return true;
00221 }
00222
00223
00224 void Transform::set_params(const Dict& d) {
00225 detect_problem_keys(d);
00226
00227 if (d.has_key_ci("type") ) set_rotation(d);
00228
00229 if (d.has_key_ci("scale")) {
00230 float scale = static_cast<float>(d.get_ci("scale"));
00231 set_scale(scale);
00232 }
00233
00234 float dx=0,dy=0,dz=0;
00235
00236 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
00237 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
00238 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
00239
00240 if ( dx != 0.0 || dy != 0.0 || dz != 0.0 ) {
00241 set_trans(dx,dy,dz);
00242 }
00243
00244 if (d.has_key_ci("mirror")) {
00245 EMObject e = d.get_ci("mirror");
00246 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
00247 throw InvalidParameterException("Error, mirror must be a bool or an int");
00248
00249 bool mirror = static_cast<bool>(e);
00250 set_mirror(mirror);
00251 }
00252 }
00253
00254
00255 void Transform::init_permissable_keys()
00256 {
00257
00258 permissable_2d_not_rot.push_back("tx");
00259 permissable_2d_not_rot.push_back("ty");
00260 permissable_2d_not_rot.push_back("scale");
00261 permissable_2d_not_rot.push_back("mirror");
00262 permissable_2d_not_rot.push_back("type");
00263
00264 permissable_3d_not_rot.push_back("tx");
00265 permissable_3d_not_rot.push_back("ty");
00266 permissable_3d_not_rot.push_back("tz");
00267 permissable_3d_not_rot.push_back("scale");
00268 permissable_3d_not_rot.push_back("mirror");
00269 permissable_3d_not_rot.push_back("type");
00270
00271 vector<string> tmp;
00272 tmp.push_back("alpha");
00273 permissable_rot_keys["2d"] = tmp;
00274
00275 tmp.clear();
00276 tmp.push_back("alt");
00277 tmp.push_back("az");
00278 tmp.push_back("phi");
00279 permissable_rot_keys["eman"] = tmp;
00280
00281 tmp.clear();
00282 tmp.push_back("psi");
00283 tmp.push_back("theta");
00284 tmp.push_back("phi");
00285 permissable_rot_keys["spider"] = tmp;
00286
00287 tmp.clear();
00288 tmp.push_back("alpha");
00289 tmp.push_back("beta");
00290 tmp.push_back("gamma");
00291 permissable_rot_keys["imagic"] = tmp;
00292
00293 tmp.clear();
00294 tmp.push_back("ztilt");
00295 tmp.push_back("xtilt");
00296 tmp.push_back("ytilt");
00297 permissable_rot_keys["xyz"] = tmp;
00298
00299 tmp.clear();
00300 tmp.push_back("phi");
00301 tmp.push_back("theta");
00302 tmp.push_back("omega");
00303 permissable_rot_keys["mrc"] = tmp;
00304
00305 tmp.clear();
00306 tmp.push_back("e0");
00307 tmp.push_back("e1");
00308 tmp.push_back("e2");
00309 tmp.push_back("e3");
00310 permissable_rot_keys["quaternion"] = tmp;
00311
00312 tmp.clear();
00313 tmp.push_back("n1");
00314 tmp.push_back("n2");
00315 tmp.push_back("n3");
00316 tmp.push_back("Omega");
00317 permissable_rot_keys["spin"] = tmp;
00318
00319 tmp.clear();
00320 tmp.push_back("n1");
00321 tmp.push_back("n2");
00322 tmp.push_back("n3");
00323 tmp.push_back("q");
00324 permissable_rot_keys["sgirot"] = tmp;
00325
00326 tmp.clear();
00327 tmp.push_back("m11");
00328 tmp.push_back("m12");
00329 tmp.push_back("m13");
00330 tmp.push_back("m21");
00331 tmp.push_back("m22");
00332 tmp.push_back("m23");
00333 tmp.push_back("m31");
00334 tmp.push_back("m32");
00335 tmp.push_back("m33");
00336 permissable_rot_keys["matrix"] = tmp;
00337 }
00338
00339
00340 void Transform::detect_problem_keys(const Dict& d) {
00341 if (permissable_rot_keys.size() == 0 ) {
00342 init_permissable_keys();
00343 }
00344
00345 vector<string> verification;
00346 vector<string> problem_keys;
00347 bool is_2d = false;
00348 if (d.has_key_ci("type") ) {
00349 string type = Util::str_to_lower((string)d["type"]);
00350 bool problem = false;
00351 if (permissable_rot_keys.find(type) == permissable_rot_keys.end() ) {
00352 problem_keys.push_back(type);
00353 problem = true;
00354 }
00355 if ( !problem ) {
00356 vector<string> perm = permissable_rot_keys[type];
00357 std::copy(perm.begin(),perm.end(),back_inserter(verification));
00358
00359 if ( type == "2d" ) {
00360 is_2d = true;
00361 std::copy(permissable_2d_not_rot.begin(),permissable_2d_not_rot.end(),back_inserter(verification));
00362 }
00363 }
00364 }
00365 if ( !is_2d ) {
00366 std::copy(permissable_3d_not_rot.begin(),permissable_3d_not_rot.end(),back_inserter(verification));
00367 }
00368
00369 for (Dict::const_iterator it = d.begin(); it != d.end(); ++it) {
00370 if ( std::find(verification.begin(),verification.end(), it->first) == verification.end() ) {
00371 problem_keys.push_back(it->first);
00372 }
00373 }
00374
00375 if (problem_keys.size() != 0 ) {
00376 string error;
00377 if (problem_keys.size() == 1) {
00378 error = "Transform Error: The \"" +problem_keys[0]+ "\" key is unsupported";
00379 } else {
00380 error = "Transform Error: The ";
00381 for(vector<string>::const_iterator cit = problem_keys.begin(); cit != problem_keys.end(); ++cit ) {
00382 if ( cit != problem_keys.begin() ) {
00383 if (cit == (problem_keys.end() -1) ) error += " and ";
00384 else error += ", ";
00385 }
00386 error += "\"";
00387 error += *cit;
00388 error += "\"";
00389 }
00390 error += " keys are unsupported";
00391 }
00392 throw InvalidParameterException(error);
00393 }
00394 }
00395
00396 void Transform::set_params_inverse(const Dict& d) {
00397 detect_problem_keys(d);
00398
00399 if (d.has_key_ci("type") ) set_rotation(d);
00400
00401 float dx=0,dy=0,dz=0;
00402 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
00403 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
00404 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
00405
00406 if ( (dx != 0.0 || dy != 0.0 || dz != 0.0) && d.has_key_ci("type") ) {
00407 Transform pre_trans;
00408 pre_trans.set_trans(dx,dy,dz);
00409
00410 Transform tmp;
00411 tmp.set_rotation(d);
00412
00413 if (d.has_key_ci("scale")) {
00414 float scale = static_cast<float>(d.get_ci("scale"));
00415 tmp.set_scale(scale);
00416 }
00417
00418 Transform solution_trans = tmp*pre_trans;
00419
00420 if (d.has_key_ci("scale")) {
00421 Transform tmp;
00422 float scale = static_cast<float>(d.get_ci("scale"));
00423 tmp.set_scale(scale);
00424 solution_trans = solution_trans*tmp;
00425 }
00426
00427 tmp = Transform();
00428 tmp.set_rotation(d);
00429 solution_trans = solution_trans*tmp;
00430 set_trans(solution_trans.get_trans());
00431 }
00432
00433 if (d.has_key_ci("scale")) {
00434 float scale = static_cast<float>(d.get_ci("scale"));
00435 set_scale(scale);
00436 }
00437
00438 if (d.has_key_ci("mirror")) {
00439 EMObject e = d.get_ci("mirror");
00440 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
00441 throw InvalidParameterException("Error, mirror must be a bool or an int");
00442
00443 bool mirror = static_cast<bool>(e);
00444 set_mirror(mirror);
00445 }
00446 invert();
00447 }
00448
00449
00450 Dict Transform::get_params(const string& euler_type) const {
00451 Dict params = get_rotation(euler_type);
00452
00453 Vec3f v = get_trans();
00454 params["tx"] = v[0]; params["ty"] = v[1];
00455
00456 string type = Util::str_to_lower(euler_type);
00457 if ( type != "2d") params["tz"] = v[2];
00458
00459 float scale = get_scale();
00460 params["scale"] = scale;
00461
00462 bool mirror = get_mirror();
00463 params["mirror"] = mirror;
00464
00465 return params;
00466 }
00467
00468
00469
00470 Dict Transform::get_params_inverse(const string& euler_type) const {
00471 Transform inv(inverse());
00472
00473 Dict params = inv.get_rotation(euler_type);
00474 Vec3f v = inv.get_pre_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 = inv.get_scale();
00481 params["scale"] = scale;
00482
00483 bool mirror = inv.get_mirror();
00484 params["mirror"] = mirror;
00485
00486 return params;
00487 }
00488
00489
00490 void Transform::set_rotation(const Dict& rotation)
00491 {
00492 detect_problem_keys(rotation);
00493 string euler_type;
00494
00495 if (!rotation.has_key_ci("type") ){
00496 throw InvalidParameterException("argument dictionary does not contain the type key");
00497 }
00498
00499 euler_type = static_cast<string>(rotation.get_ci("type"));
00500
00501
00502 double e0=0; double e1=0; double e2=0; double e3=0;
00503 double Omega=0;
00504 double az = 0;
00505 double alt = 0;
00506 double phi = 0;
00507 double cxtilt = 0;
00508 double sxtilt = 0;
00509 double cytilt = 0;
00510 double sytilt = 0;
00511 double cztilt = 0;
00512 double sztilt = 0;
00513 bool is_quaternion = 0;
00514 bool is_matrix = 0;
00515 bool is_xyz = 0;
00516
00517 bool x_mirror;
00518 float scale;
00519
00520 get_scale_and_mirror(scale,x_mirror);
00521 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
00522
00523 string type = Util::str_to_lower(euler_type);
00524 if (type == "2d") {
00525 assert_valid_2d();
00526 az = 0;
00527 alt = 0;
00528 phi = (double)rotation["alpha"] ;
00529 } else if ( type == "eman" ) {
00530
00531 az = (double)rotation["az"] ;
00532 alt = (double)rotation["alt"] ;
00533 phi = (double)rotation["phi"] ;
00534 } else if ( type == "imagic" ) {
00535
00536 az = (double)rotation["alpha"] ;
00537 alt = (double)rotation["beta"] ;
00538 phi = (double)rotation["gamma"] ;
00539 } else if ( type == "spider" ) {
00540
00541 az = (double)rotation["phi"] + 90.0;
00542 alt = (double)rotation["theta"] ;
00543 phi = (double)rotation["psi"] - 90.0;
00544 } else if ( type == "xyz" ) {
00545
00546 is_xyz = 1;
00547 cxtilt = cos(EMConsts::deg2rad*(double)rotation["xtilt"]);
00548 sxtilt = sin(EMConsts::deg2rad*(double)rotation["xtilt"]);
00549 cytilt = cos(EMConsts::deg2rad*(double)rotation["ytilt"]);
00550 sytilt = sin(EMConsts::deg2rad*(double)rotation["ytilt"]);
00551 cztilt = cos(EMConsts::deg2rad*(double)rotation["ztilt"]);
00552 sztilt = sin(EMConsts::deg2rad*(double)rotation["ztilt"]);
00553 } else if ( type == "mrc" ) {
00554
00555 az = (double)rotation["phi"] + 90.0f ;
00556 alt = (double)rotation["theta"] ;
00557 phi = (double)rotation["omega"] - 90.0f ;
00558 } else if ( type == "quaternion" ) {
00559
00560 is_quaternion = 1;
00561 e0 = (double)rotation["e0"];
00562 e1 = (double)rotation["e1"];
00563 e2 = (double)rotation["e2"];
00564 e3 = (double)rotation["e3"];
00565 } else if ( type == "spin" ) {
00566
00567 is_quaternion = 1;
00568 Omega = (double)rotation["Omega"];
00569 e0 = cos(Omega*EMConsts::deg2rad/2.0);
00570 e1 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"];
00571 e2 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"];
00572 e3 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"];
00573 } else if ( type == "sgirot" ) {
00574
00575 is_quaternion = 1;
00576 Omega = (double)rotation["q"] ;
00577 e0 = cos(Omega*EMConsts::deg2rad/2.0);
00578 e1 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"];
00579 e2 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"];
00580 e3 = sin(Omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"];
00581 } else if ( type == "matrix" ) {
00582 is_matrix = 1;
00583 matrix[0][0] = (float)rotation["m11"];
00584 matrix[0][1] = (float)rotation["m12"];
00585 matrix[0][2] = (float)rotation["m13"];
00586 matrix[1][0] = (float)rotation["m21"];
00587 matrix[1][1] = (float)rotation["m22"];
00588 matrix[1][2] = (float)rotation["m23"];
00589 matrix[2][0] = (float)rotation["m31"];
00590 matrix[2][1] = (float)rotation["m32"];
00591 matrix[2][2] = (float)rotation["m33"];
00592 } else {
00593
00594 throw InvalidStringException(euler_type, "unknown Euler Type");
00595 }
00596
00597 double azp = az*EMConsts::deg2rad;
00598 double altp = alt*EMConsts::deg2rad;
00599 double phip = phi*EMConsts::deg2rad;
00600
00601 if (!is_quaternion && !is_matrix && !is_xyz) {
00602 matrix[0][0] = (float)(cos(phip)*cos(azp) - cos(altp)*sin(azp)*sin(phip));
00603 matrix[0][1] = (float)(cos(phip)*sin(azp) + cos(altp)*cos(azp)*sin(phip));
00604 matrix[0][2] = (float)(sin(altp)*sin(phip));
00605 matrix[1][0] = (float)(-sin(phip)*cos(azp) - cos(altp)*sin(azp)*cos(phip));
00606 matrix[1][1] = (float)(-sin(phip)*sin(azp) + cos(altp)*cos(azp)*cos(phip));
00607 matrix[1][2] = (float)(sin(altp)*cos(phip));
00608 matrix[2][0] = (float)(sin(altp)*sin(azp));
00609 matrix[2][1] = (float)(-sin(altp)*cos(azp));
00610 matrix[2][2] = (float)cos(altp);
00611 }
00612 if (is_quaternion){
00613 matrix[0][0] = (float)(e0 * e0 + e1 * e1 - e2 * e2 - e3 * e3);
00614 matrix[0][1] = (float)(2.0f * (e1 * e2 + e0 * e3));
00615 matrix[0][2] = (float)(2.0f * (e1 * e3 - e0 * e2));
00616 matrix[1][0] = (float)(2.0f * (e2 * e1 - e0 * e3));
00617 matrix[1][1] = (float)(e0 * e0 - e1 * e1 + e2 * e2 - e3 * e3);
00618 matrix[1][2] = (float)(2.0f * (e2 * e3 + e0 * e1));
00619 matrix[2][0] = (float)(2.0f * (e3 * e1 + e0 * e2));
00620 matrix[2][1] = (float)(2.0f * (e3 * e2 - e0 * e1));
00621 matrix[2][2] = (float)(e0 * e0 - e1 * e1 - e2 * e2 + e3 * e3);
00622
00623 }
00624 if (is_xyz){
00625 matrix[0][0] = (float)(cytilt*cztilt);
00626 matrix[0][1] = (float)(cxtilt*sztilt+sxtilt*sytilt*cztilt);
00627 matrix[0][2] = (float)(sxtilt*sztilt-cxtilt*sytilt*cztilt);
00628 matrix[1][0] = (float)(-cytilt*sztilt);
00629 matrix[1][1] = (float)(cxtilt*cztilt-sxtilt*sytilt*sztilt);
00630 matrix[1][2] = (float)(sxtilt*cztilt+cxtilt*sytilt*sztilt);
00631 matrix[2][0] = (float)(sytilt);
00632 matrix[2][1] = (float)(-sxtilt*cytilt);
00633 matrix[2][2] = (float)(cxtilt*cytilt);
00634 }
00635
00636
00637 if (scale != 1.0f) {
00638 for(int i=0; i<3; ++i) {
00639 for(int j=0; j<3; ++j) {
00640 matrix[i][j] *= scale;
00641 }
00642 }
00643 }
00644
00645
00646 if ( x_mirror ) {
00647 for(int j=0; j<3; ++j) {
00648 matrix[0][j] *= -1.0f;
00649 }
00650 }
00651 }
00652
00653 Transform Transform::get_rotation_transform() const
00654 {
00655 Transform ret(*this);
00656 ret.set_scale(1.0);
00657 ret.set_mirror(false);
00658 ret.set_trans(0,0,0);
00659
00660 return ret;
00661 }
00662
00663 void Transform::set_rotation(const Vec3f & v)
00664 {
00665 if ( v[0] == 0 && v[1] == 0 && v[2] == 0 )
00666 throw UnexpectedBehaviorException("Can't set rotation for the null vector");
00667
00668 Vec3f v1(v);
00669 v1.normalize();
00670
00671 double theta = acos(v1[2]);
00672 double psi = atan2(v1[1],-v1[0]);
00673
00674 Dict d;
00675 d["theta"] = (double)EMConsts::rad2deg*theta;
00676 d["psi"] = (double)EMConsts::rad2deg*psi;
00677 d["phi"] = (double)0.0;
00678 d["type"] = "spider";
00679
00680 set_rotation(d);
00681
00682
00683 }
00684
00685 void Transform::rotate_origin(const Transform& by)
00686 {
00687 vector<float> multmatrix = by.get_matrix();
00688
00689 Transform result;
00690 for (int i=0; i<3; i++) {
00691 for (int j=0; j<3; j++) {
00692 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
00693 }
00694 }
00695
00696 for (int i=0; i<3; i++) {
00697 for (int j=0; j<3; j++) {
00698 matrix[i][j] = result[i][j];
00699 }
00700 }
00701 }
00702
00703 void Transform::rotate(const Transform& by)
00704 {
00705 vector<float> multmatrix = by.get_matrix();
00706
00707 Transform result;
00708 for (int i=0; i<3; i++) {
00709 for (int j=0; j<4; j++) {
00710 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
00711 }
00712 }
00713
00714 for (int i=0; i<3; i++) {
00715 for (int j=0; j<4; j++) {
00716 matrix[i][j] = result[i][j];
00717 }
00718 }
00719 }
00720
00721 Transform Transform::negate() const
00722 {
00723 Transform t(*this);
00724 for(unsigned int i = 0; i < 3; ++i) {
00725 for(unsigned int j = 0; j < 4; ++j) {
00726 t.set(i,j,t[i][j]*-1);
00727 }
00728 }
00729 return t;
00730 }
00731
00732 Transform Transform::get_hflip_transform() const {
00733
00734 Dict rot = get_rotation("eman");
00735 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
00736 rot["phi"] = 180.0f - static_cast<float>(rot["phi"]);
00737
00738 Transform ret(*this);
00739 ret.set_rotation(rot);
00740
00741 Vec3f trans = get_trans();
00742 trans[0] = -trans[0];
00743 ret.set_trans(trans);
00744
00745
00746
00747 return ret;
00748 }
00749
00750 Transform Transform::get_vflip_transform() const {
00751
00752 Dict rot = get_rotation("eman");
00753 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
00754 rot["phi"] = - static_cast<float>(rot["phi"]);
00755
00756 Transform ret(*this);
00757 ret.set_rotation(rot);
00758
00759 Vec3f trans = get_trans();
00760 trans[1] = -trans[1];
00761 ret.set_trans(trans);
00762
00763 return ret;
00764 }
00765
00766 Dict Transform::get_rotation(const string& euler_type) const
00767 {
00768 Dict result;
00769
00770
00771 float scale;
00772 bool x_mirror;
00773 get_scale_and_mirror(scale,x_mirror);
00774 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
00775
00776 double cosalt = matrix[2][2]/scale;
00777 double x_mirror_scale = (x_mirror ? -1.0f : 1.0f);
00778 double inv_scale = 1.0f/scale;
00779
00780 double az = 0;
00781 double alt = 0;
00782 double phi = 0;
00783 double phiS = 0;
00784 double psiS = 0;
00785
00786
00787
00788 if (cosalt >= 1) {
00789 alt = 0;
00790 az = 0;
00791 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00792 } else if (cosalt <= -1) {
00793 alt = 180;
00794 az = 0;
00795 phi = (double)EMConsts::rad2deg * atan2(-x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
00796 } else {
00797
00798 az = (double)EMConsts::rad2deg * atan2(scale*matrix[2][0], -scale*matrix[2][1]);
00799
00800 if (matrix[2][2]==0.0)
00801 alt = 90.0;
00802 else
00803 alt = (double)EMConsts::rad2deg * atan(sqrt((double)matrix[2][0]*matrix[2][0]+(double)matrix[2][1]*matrix[2][1])/fabs(matrix[2][2]));
00804
00805 if (matrix[2][2] * scale < 0)
00806 alt = 180.0f-alt;
00807
00808 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*(double)matrix[0][2], (double)matrix[1][2]);
00809
00810 }
00811
00812 phi = phi-360.0*floor(phi/360.0);
00813 az = az -360.0*floor(az/360.0);
00814
00815
00816 if (cosalt >= 1) {
00817 phiS = 0;
00818 psiS = phi;
00819 } else if (cosalt <= -1) {
00820 phiS = 0;
00821 psiS = phi + 180.0;
00822 } else {
00823 phiS = az - 90.0;
00824 psiS = phi + 90.0;
00825 }
00826
00827 phiS = phiS-360.0*floor(phiS/360.0);
00828 psiS = psiS-360.0*floor(psiS/360.0);
00829
00830
00831
00832 double nphi = (az-phi)/2.0;
00833
00834 double cosOover2 = cos((az+phi)*EMConsts::deg2rad/2.0) * cos(alt*EMConsts::deg2rad/2.0);
00835 double sinOover2 = sqrt(1 -cosOover2*cosOover2);
00836 double cosnTheta = sin((az+phi)*EMConsts::deg2rad/2.0) * cos(alt*EMConsts::deg2rad/2.0) / sqrt(1-cosOover2*cosOover2);
00837 double sinnTheta = sqrt(1-cosnTheta*cosnTheta);
00838 double n1 = sinnTheta*cos(nphi*EMConsts::deg2rad);
00839 double n2 = sinnTheta*sin(nphi*EMConsts::deg2rad);
00840 double n3 = cosnTheta;
00841 double xtilt = 0;
00842 double ytilt = 0;
00843 double ztilt = 0;
00844
00845
00846 if (cosOover2<0) {
00847 cosOover2*=-1; n1 *=-1; n2*=-1; n3*=-1;
00848 }
00849
00850 string type = Util::str_to_lower(euler_type);
00851
00852 result["type"] = type;
00853 if (type == "2d") {
00854 assert_valid_2d();
00855 result["alpha"] = phi;
00856 } else if (type == "eman") {
00857
00858 result["az"] = az;
00859 result["alt"] = alt;
00860 result["phi"] = phi;
00861 } else if (type == "imagic") {
00862
00863 result["alpha"] = az;
00864 result["beta"] = alt;
00865 result["gamma"] = phi;
00866 } else if (type == "spider") {
00867
00868 result["phi"] = phiS;
00869 result["theta"] = alt;
00870 result["psi"] = psiS;
00871 } else if (type == "mrc") {
00872
00873 result["phi"] = phiS;
00874 result["theta"] = alt;
00875 result["omega"] = psiS;
00876 } else if (type == "xyz") {
00877
00878 xtilt = atan2(-sin(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt),cos(EMConsts::deg2rad*alt));
00879 ytilt = asin( cos(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt));
00880 ztilt = psiS*EMConsts::deg2rad - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
00881
00882 xtilt *= EMConsts::rad2deg; ytilt *= EMConsts::rad2deg; ztilt *= EMConsts::rad2deg;
00883 xtilt = xtilt-360*.0*floor((xtilt+180.0)/360.0);
00884 ytilt = ytilt-360*.0*floor((ytilt+180.0)/360.0);
00885 ztilt = ztilt-360*.0*floor((ztilt+180.0)/360.0);
00886
00887 result["xtilt"] = xtilt;
00888 result["ytilt"] = ytilt;
00889 result["ztilt"] = ztilt;
00890 } else if (type == "quaternion") {
00891
00892 result["e0"] = cosOover2 ;
00893 result["e1"] = sinOover2 * n1 ;
00894 result["e2"] = sinOover2 * n2;
00895 result["e3"] = sinOover2 * n3;
00896 } else if (type == "spin") {
00897
00898 result["Omega"] = 2.0*EMConsts::rad2deg * acos(cosOover2);
00899 result["n1"] = n1;
00900 result["n2"] = n2;
00901 result["n3"] = n3;
00902 } else if (type == "sgirot") {
00903
00904 result["q"] = 2.0*EMConsts::rad2deg * acos(cosOover2);
00905 result["n1"] = n1;
00906 result["n2"] = n2;
00907 result["n3"] = n3;
00908 } else if (type == "matrix") {
00909
00910 result["m11"] = x_mirror_scale*matrix[0][0]*inv_scale;
00911 result["m12"] = x_mirror_scale*matrix[0][1]*inv_scale;
00912 result["m13"] = x_mirror_scale*matrix[0][2]*inv_scale;
00913 result["m21"] = matrix[1][0]*inv_scale;
00914 result["m22"] = matrix[1][1]*inv_scale;
00915 result["m23"] = matrix[1][2]*inv_scale;
00916 result["m31"] = matrix[2][0]*inv_scale;
00917 result["m32"] = matrix[2][1]*inv_scale;
00918 result["m33"] = matrix[2][2]*inv_scale;
00919 } else {
00920 throw InvalidStringException(euler_type, "unknown Euler Type");
00921 }
00922
00923 return result;
00924 }
00925
00926 void Transform::set_trans(const float& x, const float& y, const float& z)
00927 {
00928 bool x_mirror = get_mirror();
00929
00930 if (x_mirror) matrix[0][3] = -x;
00931 else matrix[0][3] = x;
00932 matrix[1][3] = y;
00933 matrix[2][3] = z;
00934 }
00935
00936 Vec3f Transform::get_trans() const
00937 {
00938
00939 bool x_mirror = get_mirror();
00940 Vec3f v;
00941 if (x_mirror) v[0] = -matrix[0][3];
00942 else v[0] = matrix[0][3];
00943 v[1] = matrix[1][3];
00944 v[2] = matrix[2][3];
00945
00946 Util::apply_precision(v[0],ERR_LIMIT);
00947 Util::apply_precision(v[1],ERR_LIMIT);
00948 Util::apply_precision(v[2],ERR_LIMIT);
00949
00950 return v;
00951 }
00952
00953 void Transform::translate(const float& tx, const float& ty, const float& tz)
00954 {
00955 bool x_mirror = get_mirror();
00956 if (x_mirror) matrix[0][3] = -matrix[0][3] + tx;
00957 else matrix[0][3] = matrix[0][3] + tx;
00958 matrix[1][3] = matrix[1][3] + ty;
00959 matrix[2][3] = matrix[2][3] + tz;
00960 }
00961
00962 Vec2f Transform::get_trans_2d() const
00963 {
00964 bool x_mirror = get_mirror();
00965 Vec2f v;
00966 if (x_mirror) v[0] = -matrix[0][3];
00967 else v[0] = matrix[0][3];
00968 v[1] = matrix[1][3];
00969 return v;
00970 }
00971
00972
00973
00974 Vec3f Transform::get_pre_trans() const
00975 {
00976 Transform T(*this);
00977 T.set_trans(0,0,0);
00978 T.invert();
00979
00980 Transform soln = T*(*this);
00981
00982 return soln.get_trans();
00983 }
00984
00985 Vec2f Transform::get_pre_trans_2d() const
00986 {
00987 Transform T(*this);
00988 T.set_trans(0,0,0);
00989 T.invert();
00990
00991 Transform soln = T*(*this);
00992
00993 return soln.get_trans_2d();
00994 }
00995
00996
00997 void Transform::set_scale(const float& new_scale) {
00998 if (new_scale <= 0) {
00999 throw InvalidValueException(new_scale,"The scale factor in a Transform object must be positive and non zero");
01000 }
01001
01002
01003
01004 float old_scale = get_scale();
01005
01006 float n_scale = new_scale;
01007 Util::apply_precision(n_scale,ERR_LIMIT);
01008
01009 float corrected_scale = n_scale/old_scale;
01010 if ( corrected_scale != 1.0 ) {
01011 for(int i = 0; i < 3; ++i ) {
01012 for(int j = 0; j < 3; ++j ) {
01013 matrix[i][j] *= corrected_scale;
01014 }
01015 }
01016 }
01017 }
01018
01019 float Transform::get_scale() const {
01020 float determinant = get_determinant();
01021 if (determinant < 0 ) determinant *= -1;
01022
01023 float scale = std::pow(determinant,1.0f/3.0f);
01024 int int_scale = static_cast<int>(scale);
01025 float scale_residual = scale-static_cast<float>(int_scale);
01026 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01027
01028 Util::apply_precision(scale, ERR_LIMIT);
01029
01030 return scale;
01031 }
01032
01033 void Transform::scale(const float& scale)
01034 {
01035 float determinant = get_determinant();
01036 if (determinant < 0) determinant *= -1.0f;
01037 float newscale = std::pow(determinant,1.0f/3.0f) + scale;
01038 if(newscale > 0.0001) set_scale(newscale);
01039 }
01040
01041 void print_matrix(gsl_matrix* M, unsigned int r, unsigned int c, const string& message ) {
01042 cout << "Message is " << message << endl;
01043 for ( unsigned int i = 0; i < r; ++i )
01044 {
01045 for ( unsigned int j = 0; j < c; ++j )
01046 {
01047 cout << gsl_matrix_get(M,i,j) << " ";
01048 }
01049 cout << endl;
01050 }
01051 }
01052
01053 void Transform::orthogonalize()
01054 {
01055 float scale;
01056 bool x_mirror;
01057 get_scale_and_mirror(scale,x_mirror);
01058 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
01059 double inv_scale = 1.0/static_cast<double>(scale);
01060 double mirror_scale = (x_mirror == true ? -1.0:1.0);
01061
01062 gsl_matrix * R = gsl_matrix_calloc(3,3);
01063 for ( unsigned int i = 0; i < 3; ++i )
01064 {
01065 for ( unsigned int j = 0; j < 3; ++j )
01066 {
01067 if (i == 0 && mirror_scale != 1.0 ) {
01068 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*mirror_scale*inv_scale );
01069 }
01070 else {
01071 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*inv_scale );
01072 }
01073 }
01074 }
01075
01076 gsl_matrix * V = gsl_matrix_calloc(3,3);
01077 gsl_vector * S = gsl_vector_calloc(3);
01078 gsl_vector * work = gsl_vector_calloc(3);
01079 gsl_linalg_SV_decomp (R, V, S, work);
01080
01081 gsl_matrix * Soln = gsl_matrix_calloc(3,3);
01082 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, R, V, 0.0, Soln);
01083
01084 for ( unsigned int i = 0; i < 3; ++i )
01085 {
01086 for ( unsigned int j = 0; j < 3; ++j )
01087 {
01088 matrix[i][j] = static_cast<float>( gsl_matrix_get(Soln,i,j) );
01089 }
01090 }
01091
01092
01093 if (scale != 1.0f) {
01094 for(int i=0; i<3; ++i) {
01095 for(int j=0; j<3; ++j) {
01096 matrix[i][j] *= scale;
01097 }
01098 }
01099 }
01100
01101
01102 if ( x_mirror ) {
01103 for(int j=0; j<3; ++j) {
01104 matrix[0][j] *= -1.0f;
01105 }
01106 }
01107
01108 gsl_matrix_free(V); gsl_matrix_free(R); gsl_matrix_free(Soln);
01109 gsl_vector_free(S); gsl_vector_free(work);
01110 }
01111
01112 void Transform::set_mirror(const bool x_mirror ) {
01113
01114 bool old_x_mirror = get_mirror();
01115 if (old_x_mirror == x_mirror) return;
01116 else {
01117
01118 for (int j = 0; j < 4; ++j ) {
01119 matrix[0][j] *= -1;
01120 }
01121 }
01122 }
01123
01124 bool Transform::get_mirror() const {
01125 float determinant = get_determinant();
01126
01127 bool x_mirror = false;
01128 if ( determinant < 0 ) x_mirror = true;
01129
01130 return x_mirror;
01131
01132 }
01133
01134 void Transform::get_scale_and_mirror(float& scale, bool& x_mirror) const {
01135
01136 float determinant = get_determinant();
01137 x_mirror = false;
01138 if ( determinant < 0 ) {
01139 x_mirror = true;
01140 determinant *= -1;
01141 }
01142 if (determinant != 1 ) {
01143 scale = std::pow(determinant,1.0f/3.0f);
01144 int int_scale = static_cast<int>(scale);
01145 float scale_residual = scale-static_cast<float>(int_scale);
01146 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
01147 }
01148 else scale = 1;
01149
01150 Util::apply_precision(scale,ERR_LIMIT);
01151 }
01152
01153 float Transform::get_determinant() const
01154 {
01155 float det;
01156 double det2;
01157 det2 = matrix[0][0]*((double)matrix[1][1]*matrix[2][2]-(double)matrix[2][1]*matrix[1][2]);
01158 det2 -= matrix[0][1]*((double)matrix[1][0]*matrix[2][2]-(double)matrix[2][0]*matrix[1][2]);
01159 det2 += matrix[0][2]*((double)matrix[1][0]*matrix[2][1]-(double)matrix[2][0]*matrix[1][1]);
01160
01161 det = (float)det2;
01162 Util::apply_precision(det,ERR_LIMIT);
01163
01164 return det;
01165 }
01166
01167 void Transform::invert() {
01168
01169 double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
01170 double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
01171 double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
01172 double v0 = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
01173
01174 double cof00 = m11*m22-m12*m21;
01175 double cof11 = m22*m00-m20*m02;
01176 double cof22 = m00*m11-m01*m10;
01177 double cof01 = m10*m22-m20*m12;
01178 double cof02 = m10*m21-m20*m11;
01179 double cof12 = m00*m21-m01*m20;
01180 double cof10 = m01*m22-m02*m21;
01181 double cof20 = m01*m12-m02*m11;
01182 double cof21 = m00*m12-m10*m02;
01183
01184 double det = m00* cof00 + m02* cof02 -m01*cof01;
01185
01186 matrix[0][0] = (float)(cof00/det);
01187 matrix[0][1] = - (float)(cof10/det);
01188 matrix[0][2] = (float)(cof20/det);
01189 matrix[1][0] = - (float)(cof01/det);
01190 matrix[1][1] = (float)(cof11/det);
01191 matrix[1][2] = - (float)(cof21/det);
01192 matrix[2][0] = (float)(cof02/det);
01193 matrix[2][1] = - (float)(cof12/det);
01194 matrix[2][2] = (float)(cof22/det);
01195
01196 matrix[0][3] = (float)((- cof00*v0 + cof10*v1 - cof20*v2)/det);
01197 matrix[1][3] = (float)(( cof01*v0 - cof11*v1 + cof21*v2)/det);
01198 matrix[2][3] = (float)((- cof02*v0 + cof12*v1 - cof22*v2)/det);
01199 }
01200
01201 Transform Transform::inverse() const {
01202 Transform t(*this);
01203 t.invert();
01204 return t;
01205 }
01206
01207 void Transform::transpose_inplace() {
01208 float tempij;
01209 for (int i = 0; i < 3; i++) {
01210 for (int j = 0; j < i; j++) {
01211 if (i != j) {
01212 tempij= matrix[i][j];
01213 matrix[i][j] = matrix[j][i];
01214 matrix[j][i] = tempij;
01215 }
01216 }
01217 }
01218 }
01219
01220 Transform Transform::transpose() const {
01221 Transform t(*this);
01222 t.transpose_inplace();
01223 return t;
01224 }
01225
01226
01227 Transform EMAN::operator*(const Transform & M2, const Transform & M1)
01228 {
01229 Transform result;
01230 for (int i=0; i<3; i++) {
01231 for (int j=0; j<4; j++) {
01232 result[i][j] = M2[i][0] * M1[0][j] + M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
01233 }
01234 result[i][3] += M2[i][3];
01235 }
01236
01237 return result;
01238 }
01239
01240 void Transform::assert_valid_2d() const {
01241 int rotation_error = 0;
01242 int translation_error = 0;
01243 if (fabs(matrix[2][0]) > ERR_LIMIT) rotation_error++;
01244 if (fabs(matrix[2][1]) > ERR_LIMIT) rotation_error++;
01245 if (fabs(matrix[2][3]) > ERR_LIMIT) translation_error++;
01246 if (fabs(matrix[0][2]) > ERR_LIMIT) rotation_error++;
01247 if (fabs(matrix[1][2]) > ERR_LIMIT) rotation_error++;
01248
01249 if ( translation_error && rotation_error ) {
01250 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and 3D translations. This object can not be considered 2D");
01251 } else if ( translation_error ) {
01252 throw UnexpectedBehaviorException("Error, the internal matrix contains a non zero z component for a 3D translation. This object can not be considered 2D");
01253 }
01254 else if ( rotation_error ) {
01255 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and this object can not be considered 2D");
01256 }
01257
01258 }
01259
01260
01261
01262 Transform Transform::get_sym(const string & sym_name, int n) const
01263 {
01264 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01265 Transform ret;
01266 ret = (*this) * sym->get_sym(n);
01267 delete sym;
01268 return ret;
01269 }
01270
01271 vector<Transform > Transform::get_sym_proj(const string & sym_name) const
01272 {
01273 vector<Transform> ret;
01274 Transform t;
01275 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01276 int nsym = sym->get_nsym();
01277 int n = nsym;
01278
01279 if ((sym_name[0] == 'c' || sym_name[0] == 'd' ) && fabs(matrix[2][2]) < 1.e-6){
01280
01281 Dict d1,d2;
01282 d2["theta"] = (double)90.0;
01283 d2["psi"] = (double)0.0;
01284 d2["phi"] = (double)0.0;
01285 d2["type"] = "spider";
01286 d1 = this->get_rotation("spider");
01287
01288 if (sym_name[0] == 'c') {
01289 if( nsym%2 == 0) n = nsym/2;
01290
01291 for (int k=0;k<n;k++) {
01292 d2["phi"] = (double)d1["phi"] + k*double(360.0)/ nsym;
01293 d2["psi"] = d1["psi"];
01294 t.set_rotation(d2);
01295 ret.push_back( t );
01296 }
01297
01298 }
01299 else {
01300 nsym = nsym/2;
01301
01302 if (nsym%2 == 0) {
01303 n = nsym;
01304 float cos_phi = cos( EMConsts::deg2rad*360.0/2/nsym );
01305
01306 for (int k=0;k<n;k++){
01307
01308 if(k%2==0) {
01309
01310 d2["phi"] = (double)d1["phi"] + k/2*double(360.0)/ nsym;
01311 d2["psi"] = d1["psi"];
01312 t.set_rotation(d2);
01313 ret.push_back( t );
01314 }
01315 else {
01316
01317 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01318
01319 d2["phi"] = k/2*double(360.0)/ nsym +180 - (double)d1["phi"];
01320 d2["psi"] = (double)d1["psi"] + 180;
01321 t.set_rotation(d2);
01322 ret.push_back( t );
01323 }
01324 }
01325
01326 }
01327 }
01328
01329
01330
01331 else {
01332 n = nsym*2;
01333 float cos_phi = cos( EMConsts::deg2rad*360.0/4/nsym );
01334 for (int k=0;k<n;k++){
01335
01336 if(k%4==0) {
01337
01338 d2["phi"] = (double)d1["phi"] + k/4*360.0/ nsym;
01339 d2["psi"] = (double)d1["psi"];
01340 t.set_rotation(d2);
01341 ret.push_back( t );
01342 }
01343 else if( k%4 ==1) {
01344 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ){
01345
01346 d2["phi"] = k/4*360.0/nsym + 360.0/2/nsym+180 - (double)d1["phi"];
01347 d2["psi"] = (double)d1["psi"] + 180;
01348 t.set_rotation(d2);
01349 ret.push_back( t );
01350 }
01351
01352 }
01353
01354 else if( k%4 ==2) {
01355
01356 d2["phi"] = k/4*360.0/ nsym+360.0/2/nsym+180 + (double)d1["phi"];
01357 d2["psi"] = (double)d1["psi"];
01358 t.set_rotation(d2);
01359 ret.push_back( t );
01360
01361 }
01362
01363 else if( k%4 ==3) {
01364 if( ( fabs(1.0-matrix[2][0])>1.0e-6 )&& fabs( matrix[2][0]-cos_phi)>1.0e-6 ) {
01365 d2["phi"] = k/4*360.0/nsym+ 2.0*360.0/2/nsym - (double)d1["phi"];
01366 d2["psi"] = (double)d1["psi"] + 180;
01367 t.set_rotation(d2);
01368 ret.push_back( t );
01369 }
01370 }
01371
01372 }
01373 }
01374
01375 }
01376
01377 }
01378 else {
01379 for (int k=0;k<nsym;k++) {
01380 t = sym->get_sym(k);
01381 ret.push_back( (*this) * t );
01382 }
01383
01384 }
01385 delete sym;
01386 return ret;
01387 }
01388
01389
01390 int Transform::get_nsym(const string & sym_name)
01391 {
01392 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
01393 int nsym = sym->get_nsym();
01394 delete sym;
01395 return nsym;
01396 }
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01416
01417
01418
01419
01420
01421
01422
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01454
01455
01456
01457
01458
01459
01460
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01524
01525
01526
01527
01528
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01547
01548
01549
01550
01551
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
01580
01581
01582
01583
01584
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610
01613
01614
01615
01616
01617
01618
01619
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
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
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
01714
01715
01716
01717
01718
01719
01720
01721
01722
01724
01725
01726
01727
01728
01729
01730
01731
01732
01735
01736
01737
01738
01739
01740
01741
01742
01743
01744
01745
01746
01747
01748
01749
01750
01751
01821
01822
01823
01824
01825
01826
01827
01828
01829
01830
01831
01833
01834
01835
01836
01837
01838
01839
01840
01841
01842
01843
01844
01845
01846
01847
01848
01849
01850
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892
01893
01894
01895
01896
01897
01898
01899
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
01936
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
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
02080
02081
02082
02083
02084
02085
02086
02088
02090
02091
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
02129
02130
02131
02132
02133
02134
02136
02137
02138
02139
02140
02141
02142
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02156
02157
02158
02159
02160
02161
02162
02163
02164
02165
02166
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
02202
02203
02204
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217
02218
02219
02220
02221
02222
02223
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234
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
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
02357
02358
02359
02360
02361
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
02423
02424
02425
02426
02427
02428
02429
02430
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