Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

transform.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "transform.h"
00037 #include "util.h"
00038 #include "emobject.h"
00039 #include <cctype> // for std::tolower
00040 #include <cstring>  // for memcpy
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> // for std::transform
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 //const float Transform3D::ERR_LIMIT = 0.000001f;
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; // 5 fold to a 3 fold to a 2 fold
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; // 3 fold to a 2 fold
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 //              transform_type = that.transform_type;
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 //      transform_type = UNKNOWN;
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"));// Warning, will throw
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         // Get these before anything changes so we can apply them again after the rotation is set
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 //              validate_and_set_type(THREED);
00531                 az  = (double)rotation["az"] ;
00532                 alt = (double)rotation["alt"]  ;
00533                 phi = (double)rotation["phi"] ;
00534         } else if ( type == "imagic" ) {
00535 //              validate_and_set_type(THREED);
00536                 az  = (double)rotation["alpha"] ;
00537                 alt = (double)rotation["beta"]  ;
00538                 phi = (double)rotation["gamma"] ;
00539         } else if ( type == "spider" ) {
00540 //              validate_and_set_type(THREED);
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 //              validate_and_set_type(THREED);
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 //              validate_and_set_type(THREED);
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 //              validate_and_set_type(THREED);
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 //              validate_and_set_type(THREED);
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 //              validate_and_set_type(THREED);
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 //              transform_type = UNKNOWN;
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                 // keep in mind matrix[0][2] is M13 gives an e0 e2 piece, etc
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         // Apply scale if it existed previously
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         // Apply post x mirroring if it was applied previously
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         //ret.orthogonalize(); // ?
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]); // in radians
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         // First Multiply and put the result in a temp matrix
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         //Then put the result from the tmep matrix in the original one
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         // First Multiply and put the result in a temp matrix
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         //Then put the result from the tmep matrix in the original one
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); // Is the identity
00739         ret.set_rotation(rot);
00740 
00741         Vec3f trans = get_trans();
00742         trans[0] = -trans[0];
00743         ret.set_trans(trans);
00744 
00745 //      ret.set_mirror(self.get_mirror());
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         //float max = 1 - ERR_LIMIT;
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;  // like az  (but in SPIDER ZYZ)
00784         double psiS = 0;  // like phi (but in SPIDER ZYZ)
00785 
00786         // get alt, az, phi in EMAN convention
00787 
00788         if (cosalt >= 1) {  // that is, alt close to 0
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) {  // that is, alt close to 180
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 {   // for non exceptional cases:  0 < alt < 180
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         } // ends separate cases: alt close to 0, 180, or neither
00811 
00812         phi = phi-360.0*floor(phi/360.0);
00813         az  = az -360.0*floor(az/360.0);
00814 
00815 //  get phiS, psiS (SPIDER)
00816         if (cosalt >= 1) {  // that is, alt close to 0
00817                 phiS = 0;
00818                 psiS = phi;
00819         } else if (cosalt <= -1) {  // that is, alt close to 180
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 //   do some quaternionic stuff here
00831 
00832         double nphi = (az-phi)/2.0;
00833     // The next is also e0
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 //              assert_consistent_type(THREED);
00858                 result["az"]  = az;
00859                 result["alt"] = alt;
00860                 result["phi"] = phi;
00861         } else if (type == "imagic") {
00862 //              assert_consistent_type(THREED);
00863                 result["alpha"] = az;
00864                 result["beta"]  = alt;
00865                 result["gamma"] = phi;
00866         } else if (type == "spider") {
00867 //              assert_consistent_type(THREED);
00868                 result["phi"]   = phiS;  // The first Euler like az
00869                 result["theta"] = alt;
00870                 result["psi"]   = psiS;
00871         } else if (type == "mrc") {
00872 //              assert_consistent_type(THREED);
00873                 result["phi"]   = phiS;
00874                 result["theta"] = alt;
00875                 result["omega"] = psiS;
00876         } else if (type == "xyz") {               // need to double-check these 3 equations ********
00877 //              assert_consistent_type(THREED);
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);  //already in range [-90,90] but anyway...
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 //              assert_consistent_type(THREED);
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 //              assert_consistent_type(THREED);
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 //              assert_consistent_type(THREED);
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 //              assert_consistent_type(THREED);
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         // No type asserted
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 //      soln.printme();
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 //      soln.printme();
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         // Transform = MTSR (Mirroring, Translation, Scaling, Rotate)
01002         // So changing the scale boils down to this....
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); // If scale ~ 0 things blowup, so we need a little fudge factor
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); // Now R is U of the SVD R = USV^T
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         // Apply scale if it existed previously
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         // Apply post x mirroring if it was applied previouslys
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; // The user is setting the same value
01116         else {
01117                 // Toggle the mirroring operation
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)     // YYY
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 //      if (fabs(matrix[2][2]-1.0) >ERR_LIMIT) rotation_error++; 
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                                                         //cout<<"jumped into"<<endl;
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 //Transform3D::Transform3D()  //    C1
01400 //{
01401 //      init();
01402 //}
01403 //
01404 //Transform3D::Transform3D( const Transform3D& rhs )
01405 //{
01406 //    for( int i=0; i < 4; ++i )
01407 //    {
01408 //        for( int j=0; j < 4; ++j )
01409 //      {
01410 //          matrix[i][j] = rhs.matrix[i][j];
01411 //      }
01412 //    }
01413 //}
01414 //
01416 //Transform3D::Transform3D(const float& az, const float& alt, const float& phi)
01417 //{
01418 //      init();
01419 //      set_rotation(az,alt,phi);
01420 //}
01421 //
01422 //
01424 //Transform3D::Transform3D(const float& az, const float& alt, const float& phi, const Vec3f& posttrans )
01425 //{
01426 //      init(); // This is called in set_rotation
01427 //      set_rotation(az,alt,phi);
01428 //      set_posttrans(posttrans);
01429 //}
01430 //
01431 //Transform3D::Transform3D(const float& m11, const float& m12, const float& m13,
01432 //                                               const float& m21, const float& m22, const float& m23,
01433 //                                               const float& m31, const float& m32, const float& m33)
01434 //{
01435 //      init();
01436 //      set_rotation(m11,m12,m13,m21,m22,m23,m31,m32,m33);
01437 //}
01438 //
01440 //Transform3D::Transform3D(EulerType euler_type, const float& a1, const float& a2, const float& a3)
01441 //{
01442 //      init();
01443 //      set_rotation(euler_type,a1,a2,a3);
01444 //}
01445 //
01446 //Transform3D::Transform3D(EulerType euler_type, const float& a1, const float& a2, const float& a3, const float& a4)
01447 //{
01448 //      init();
01449 //      set_rotation(euler_type,a1,a2,a3,a4);
01450 //}
01451 //
01452 //
01454 //Transform3D::Transform3D(EulerType euler_type, const Dict& rotation)  //YYY
01455 //{
01456 //      init();
01457 //      set_rotation(euler_type,rotation);
01458 //}
01459 //
01460 //
01462 //
01463 //Transform3D::Transform3D(  const Vec3f& pretrans,  const float& az, const float& alt, const float& phi, const Vec3f& posttrans )  //YYY  by default EMAN
01464 //{
01465 //      init();
01466 //      set_pretrans(pretrans);
01467 //      set_rotation(az,alt,phi);
01468 //      set_posttrans(posttrans);
01469 //}
01470 //
01471 //
01472 //
01473 //
01474 //Transform3D::~Transform3D()
01475 //{
01476 //}
01477 //
01478 //
01479 //
01480 //void Transform3D::to_identity()
01481 //{
01485 //
01486 //      for(int i=0; i<4; ++i) {
01487 //              for(int j=0; j<4; ++j) {
01488 //                      if(i==j) {
01489 //                              matrix[i][j] = 1;
01490 //                      }
01491 //                      else {
01492 //                              matrix[i][j] = 0;
01493 //                      }
01494 //              }
01495 //      }
01496 //      post_x_mirror = false;
01497 //      set_center(Vec3f(0,0,0));
01498 //}
01499 //
01500 //
01501 //
01502 //bool Transform3D::is_identity()  // YYY
01503 //{
01504 //      for (int i=0; i<4; i++) {
01505 //              for (int j=0; j<4; j++) {
01506 //                      if (i==j && matrix[i][j]!=1.0) return 0;
01507 //                      if (i!=j && matrix[i][j]!=0.0) return 0;
01508 //              }
01509 //      }
01510 //      return 1;
01511 //}
01512 //
01513 //
01514 //void Transform3D::set_center(const Vec3f & center) //YYN
01515 //{
01516 //      set_pretrans( Vec3f(0,0,0)-center);
01517 //      for (int i = 0; i < 3; i++) {
01518 //              matrix[i][3]=center[i];
01519 //      }
01520 //}
01521 //
01524 //void Transform3D::init()  // M1
01525 //{
01526 //      to_identity();
01527 //}
01528 //
01530 //
01531 //void Transform3D::set_pretrans(const float& dx, const float& dy, const float& dz) // YYY
01532 //{    set_pretrans( Vec3f(dx,dy,dz)); }
01533 //
01534 //
01535 //void Transform3D::set_pretrans(const float& dx, const float& dy) // YYY
01536 //{    set_pretrans( Vec3f(dx,dy,0)); }
01537 //
01538 //void Transform3D::set_pretrans(const Vec2f& pretrans) // YYY
01539 //{    set_pretrans( Vec3f(pretrans[0],pretrans[1],0)); }
01540 //
01541 //void Transform3D::set_pretrans(const Vec3f & preT)  // flag=1 means keep the old value of total trans
01542 //{
01543 //              int flag=0;
01544 //
01547 //    if (flag==0){
01548 //              matrix[0][3] = matrix[3][0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2]  ;
01549 //              matrix[1][3] = matrix[3][1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2]  ;
01550 //              matrix[2][3] = matrix[3][2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2]  ;
01551 //      }
01553 //    if (flag==1){
01554 //              matrix[3][0] = matrix[0][3] - (matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2])  ;
01555 //              matrix[3][1] = matrix[1][3] - (matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2])  ;
01556 //              matrix[3][2] = matrix[2][3] - (matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2])  ;
01557 //      }
01558 //}
01559 //
01560 //
01561 //void Transform3D::set_posttrans(const float& dx, const float& dy, const float& dz) // YYY
01562 //{    set_posttrans( Vec3f(dx,dy,dz)); }
01563 //
01564 //
01565 //void Transform3D::set_posttrans(const float& dx, const float& dy) // YYY
01566 //{    set_posttrans( Vec3f(dx,dy,0)); }
01567 //
01568 //void Transform3D::set_posttrans(const Vec2f& posttrans) // YYY
01569 //{    set_pretrans( Vec3f(posttrans[0],posttrans[1],0)); }
01570 //
01571 //void Transform3D::set_posttrans(const Vec3f & posttrans) // flag=1 means keep the old value of total trans
01572 //{
01573 //      int flag=0;
01574 //    Vec3f preT   = get_pretrans(0) ;
01575 //      for (int i = 0; i < 3; i++) {
01576 //              matrix[3][i] = posttrans[i];
01577 //      }
01580 //      if (flag==0) {
01581 //              matrix[0][3] = matrix[3][0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2]  ;
01582 //              matrix[1][3] = matrix[3][1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2]  ;
01583 //              matrix[2][3] = matrix[3][2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2]  ;
01584 //      }
01586 //      if (flag==1) { // Don't do anything
01587 //      }
01588 //}
01589 //
01590 //
01591 //
01592 //
01593 //void Transform3D::apply_scale(const float& scale)    // YYY
01594 //{
01595 //      for (int i = 0; i < 3; i++) {
01596 //              for (int j = 0; j < 4; j++) {
01597 //                      matrix[i][j] *= scale;
01598 //              }
01599 //      }
01600 //      for (int j = 0; j < 3; j++) {
01601 //              matrix[3][j] *= scale;
01602 //      }
01603 //}
01604 //
01605 //void Transform3D::orthogonalize()  // YYY
01606 //{
01607 //      //EulerType EMAN;
01608 //      float scale = get_scale() ;
01609 //      float inverseScale= 1/scale ;
01610 //      apply_scale(inverseScale);
01613 //}
01614 //
01615 //
01616 //void Transform3D::transpose()  // YYY
01617 //{
01618 //      float tempij;
01619 //      for (int i = 0; i < 3; i++) {
01620 //              for (int j = 0; j < i; j++) {
01621 //                      tempij= matrix[i][j];
01622 //                      matrix[i][j] = matrix[j][i];
01623 //                      matrix[j][i] = tempij;
01624 //              }
01625 //      }
01626 //}
01627 //
01628 //void Transform3D::set_scale(const float& scale)    // YYY
01629 //{
01630 //      float OldScale= get_scale();
01631 //      float Scale2Apply = scale/OldScale;
01632 //      apply_scale(Scale2Apply);
01633 //}
01634 //
01635 //float Transform3D::get_mag() const //
01636 //{
01637 //      EulerType eulertype= SPIN ;
01638 //      Dict AA= get_rotation(eulertype);
01639 //      return AA["Omega"];
01640 //}
01641 //
01642 //Vec3f Transform3D::get_finger() const //
01643 //{
01644 //      EulerType eulertype= SPIN ;
01645 //      Dict AA= get_rotation(eulertype);
01646 //      return Vec3f(AA["n1"],AA["n2"],AA["n3"]);
01647 //}
01648 //
01649 //Vec3f Transform3D::get_posttrans(int flag) const    //
01650 //{
01651 //      if (flag==0){
01652 //              return Vec3f(matrix[3][0], matrix[3][1], matrix[3][2]);
01653 //      }
01654 //      // otherwise as if all the translation was post
01655 //      return Vec3f(matrix[0][3], matrix[1][3], matrix[2][3]);
01656 //}
01657 //
01658 //Vec3f Transform3D::get_total_posttrans() const {
01659 //      return get_posttrans(1);
01660 //}
01661 //
01662 //Vec3f Transform3D::get_total_pretrans() const {
01663 //      return get_pretrans(1);
01664 //}
01665 //
01666 //
01667 //Vec3f Transform3D::get_pretrans(int flag) const    // Fix Me
01668 //{
01670 //
01671 //      Vec3f pretrans;
01672 //      Vec3f posttrans(matrix[3][0], matrix[3][1], matrix[3][2]);
01673 //      Vec3f tottrans(matrix[0][3], matrix[1][3], matrix[2][3]);
01674 //      Vec3f totminuspost;
01675 //
01676 //      totminuspost = tottrans;
01677 //      if (flag==0) {
01678 //              totminuspost = tottrans-posttrans;
01679 //      }
01680 //
01681 //      Transform3D Rinv = inverse();
01682 //      for (int i=0; i<3; i++) {
01683 //                float ptnow=0;
01684 //              for (int j=0; j<3; j++) {
01685 //                      ptnow +=   Rinv.matrix[i][j]* totminuspost[j] ;
01686 //              }
01687 //              pretrans.set_value_at(i,ptnow) ;  //
01688 //      }
01689 //      return pretrans;
01690 //}
01691 //
01692 //
01693 // Vec3f Transform3D::get_center() const  // YYY
01694 // {
01695 //      return Vec3f();
01696 // }
01697 //
01698 //
01699 //
01700 //Vec3f Transform3D::get_matrix3_col(int i) const     // YYY
01701 //{
01702 //      return Vec3f(matrix[0][i], matrix[1][i], matrix[2][i]);
01703 //}
01704 //
01705 //
01706 //Vec3f Transform3D::get_matrix3_row(int i) const     // YYY
01707 //{
01708 //      return Vec3f(matrix[i][0], matrix[i][1], matrix[i][2]);
01709 //}
01710 //
01711 //Vec3f Transform3D::transform(const Vec3f & v3f) const     // YYY
01712 //{
01714 //      float x = matrix[0][0] * v3f[0] + matrix[0][1] * v3f[1] + matrix[0][2] * v3f[2] + matrix[0][3] ;
01715 //      float y = matrix[1][0] * v3f[0] + matrix[1][1] * v3f[1] + matrix[1][2] * v3f[2] + matrix[1][3] ;
01716 //      float z = matrix[2][0] * v3f[0] + matrix[2][1] * v3f[1] + matrix[2][2] * v3f[2] + matrix[2][3] ;
01717 //      return Vec3f(x, y, z);
01718 //}
01719 //
01720 //
01721 //Vec3f Transform3D::rotate(const Vec3f & v3f) const     // YYY
01722 //{
01724 //      float x = matrix[0][0] * v3f[0] + matrix[0][1] * v3f[1] + matrix[0][2] * v3f[2]  ;
01725 //      float y = matrix[1][0] * v3f[0] + matrix[1][1] * v3f[1] + matrix[1][2] * v3f[2]  ;
01726 //      float z = matrix[2][0] * v3f[0] + matrix[2][1] * v3f[1] + matrix[2][2] * v3f[2]  ;
01727 //      return Vec3f(x, y, z);
01728 //}
01729 //
01730 //
01731 //Transform3D EMAN::operator*(const Transform3D & M2, const Transform3D & M1)     // YYY
01732 //{
01735 //      Transform3D resultant;
01736 //      for (int i=0; i<3; i++) {
01737 //              for (int j=0; j<4; j++) {
01738 //                      resultant[i][j] = M2[i][0] * M1[0][j] +  M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
01739 //              }
01740 //              resultant[i][3] += M2[i][3];  // add on the new translation (not included above)
01741 //      }
01742 //
01743 //      for (int j=0; j<3; j++) {
01744 //              resultant[3][j] = M2[3][j];
01745 //      }
01746 //
01747 //      return resultant; // This will have the post_trans of M2
01748 //}
01749 
01750 /*             Here starts the pure rotation stuff */
01751 
01821 //
01822 //void Transform3D::set_rotation(const float& az, const float& alt, const float& phi )
01823 //{
01824 //      EulerType euler_type=EMAN;
01825 //      Dict rot;
01826 //      rot["az"]  = az;
01827 //      rot["alt"] = alt;
01828 //      rot["phi"] = phi;
01829 //      set_rotation(euler_type, rot);
01830 //}
01831 //
01833 //void Transform3D::set_rotation(EulerType euler_type, const float& a1, const float& a2, const float& a3)
01834 //{
01835 //      init();
01836 //      Dict rot;
01837 //      switch(euler_type) {
01838 //              case EMAN:
01839 //                      rot["az"]  = a1;
01840 //                      rot["alt"] = a2;
01841 //                      rot["phi"] = a3;
01842 //                      break;
01843 //              case SPIDER:
01844 //                      rot["phi"]   = a1;
01845 //                      rot["theta"] = a2;
01846 //                      rot["psi"]   = a3;
01847 //                      break;
01848 //              case IMAGIC:
01849 //                      rot["alpha"]   = a1;
01850 //                      rot["beta"] = a2;
01851 //                      rot["gamma"]   = a3;
01852 //                      break;
01853 //              case MRC:
01854 //                      rot["phi"]   = a1;
01855 //                      rot["theta"] = a2;
01856 //                      rot["omega"]   = a3;
01857 //                      break;
01858 //              case XYZ:
01859 //                      rot["xtilt"]   = a1;
01860 //                      rot["ytilt"] = a2;
01861 //                      rot["ztilt"]   = a3;
01862 //                      break;
01863 //              default:
01864 //              throw InvalidValueException(euler_type, "cannot instantiate this Euler Type");
01865 //      }  // ends switch euler_type
01866 //      set_rotation(euler_type, rot);
01867 //}
01868 //
01870 //void Transform3D::set_rotation(EulerType euler_type, const float& a1, const float& a2, const float& a3, const float& a4)
01871 //{
01872 //      init();
01873 //      Dict rot;
01874 //      switch(euler_type) {
01875 //              case QUATERNION:
01876 //                      rot["e0"]  = a1;
01877 //                      rot["e1"] = a2;
01878 //                      rot["e2"] = a3;
01879 //                      rot["e3"] = a4;
01880 //                      break;
01881 //              case SGIROT:
01882 //                      rot["q"]  = a1;
01883 //                      rot["n1"] = a2;
01884 //                      rot["n2"] = a3;
01885 //                      rot["n3"] = a4;
01886 //              case SPIN:
01887 //                      rot["Omega"]  = a1;
01888 //                      rot["n1"] = a2;
01889 //                      rot["n2"] = a3;
01890 //                      rot["n3"] = a4;
01891 //                      break;
01892 //              default:
01893 //                      throw InvalidValueException(euler_type, "cannot instantiate this Euler Type");
01894 //      }  // ends switch euler_type
01895 //      set_rotation(euler_type, rot);
01896 //}
01897 //
01898 //void Transform3D::set_rotation(EulerType euler_type, const Dict& rotation)
01899 //{
01900 //      double e0  = 0; double e1=0; double e2=0; double e3=0;
01901 //      double Omega=0;
01902 //      double az  = 0;
01903 //      double alt = 0;
01904 //      double phi = 0;
01905 //      double cxtilt = 0;
01906 //      double sxtilt = 0;
01907 //      double cytilt = 0;
01908 //      double sytilt = 0;
01909 //      bool is_quaternion = 0;
01910 //      bool is_matrix = 0;
01911 //
01912 //      switch(euler_type) {
01913 //      case EMAN:
01914 //              az  = (double)rotation["az"] ;
01915 //              alt = (double)rotation["alt"]  ;
01916 //              phi = (double)rotation["phi"] ;
01917 //              break;
01918 //      case IMAGIC:
01919 //              az  = (double)rotation["alpha"] ;
01920 //              alt = (double)rotation["beta"]  ;
01921 //              phi = (double)rotation["gamma"] ;
01922 //              break;
01923 //
01924 //      case SPIDER:
01925 //              az =  (double)rotation["phi"]    + 90.0;
01926 //              alt = (double)rotation["theta"] ;
01927 //              phi = (double)rotation["psi"]    - 90.0;
01928 //              break;
01929 //
01930 //      case XYZ:
01931 //              cxtilt = cos( (M_PI/180.0f)*(double)rotation["xtilt"]);
01932 //              sxtilt = sin( (M_PI/180.0f)*(double)rotation["xtilt"]);
01933 //              cytilt = cos( (M_PI/180.0f)*(double)rotation["ytilt"]);
01934 //              sytilt = sin( (M_PI/180.0f)*(double)rotation["ytilt"]);
01935 //              az =  (180.0f/M_PI)*atan2(-cytilt*sxtilt,sytilt)   + 90.0f ;
01936 //              alt = (180.0f/M_PI)*acos(cytilt*cxtilt)  ;
01937 //              phi = (double)rotation["ztilt"] +(180.0f/M_PI)*atan2(sxtilt,cxtilt*sytilt)   - 90.0f ;
01938 //              break;
01939 //
01940 //      case MRC:
01941 //              az  = (double)rotation["phi"]   + 90.0f ;
01942 //              alt = (double)rotation["theta"] ;
01943 //              phi = (double)rotation["omega"] - 90.0f ;
01944 //              break;
01945 //
01946 //      case QUATERNION:
01947 //              is_quaternion = 1;
01948 //              e0 = (double)rotation["e0"];
01949 //              e1 = (double)rotation["e1"];
01950 //              e2 = (double)rotation["e2"];
01951 //              e3 = (double)rotation["e3"];
01952 //              break;
01953 //
01954 //      case SPIN:
01955 //              is_quaternion = 1;
01956 //              Omega = (double)rotation["Omega"];
01957 //              e0 = cos(Omega*M_PI/360.0f);
01958 //              e1 = sin(Omega*M_PI/360.0f)* (double)rotation["n1"];
01959 //              e2 = sin(Omega*M_PI/360.0f)* (double)rotation["n2"];
01960 //              e3 = sin(Omega*M_PI/360.0f)* (double)rotation["n3"];
01961 //              break;
01962 //
01963 //      case SGIROT:
01964 //              is_quaternion = 1;
01965 //              Omega = (double)rotation["q"]  ;
01966 //              e0 = cos(Omega*M_PI/360.0f);
01967 //              e1 = sin(Omega*M_PI/360.0f)* (double)rotation["n1"];
01968 //              e2 = sin(Omega*M_PI/360.0f)* (double)rotation["n2"];
01969 //              e3 = sin(Omega*M_PI/360.0f)* (double)rotation["n3"];
01970 //              break;
01971 //
01972 //      case MATRIX:
01973 //              is_matrix = 1;
01974 //              matrix[0][0] = (float)rotation["m11"]  ;
01975 //              matrix[0][1] = (float)rotation["m12"]  ;
01976 //              matrix[0][2] = (float)rotation["m13"]  ;
01977 //              matrix[1][0] = (float)rotation["m21"]  ;
01978 //              matrix[1][1] = (float)rotation["m22"]  ;
01979 //              matrix[1][2] = (float)rotation["m23"]  ;
01980 //              matrix[2][0] = (float)rotation["m31"]  ;
01981 //              matrix[2][1] = (float)rotation["m32"]  ;
01982 //              matrix[2][2] = (float)rotation["m33"]  ;
01983 //              break;
01984 //
01985 //      default:
01986 //              throw InvalidValueException(euler_type, "unknown Euler Type");
01987 //      }  // ends switch euler_type
01988 //
01989 //
01990 //      Vec3f postT  = get_posttrans( ) ;
01991 //      Vec3f preT   = get_pretrans( ) ;
01992 //
01993 //
01994 //      double azp  = fmod(az,360.0)*M_PI/180.0;
01995 //      double altp  = alt*M_PI/180.0;
01996 //      double phip = fmod(phi,360.0)*M_PI/180.0;
01997 //
01998 //      if (!is_quaternion && !is_matrix) {
01999 //              matrix[0][0] =  (float)(cos(phip)*cos(azp) - cos(altp)*sin(azp)*sin(phip));
02000 //              matrix[0][1] =  (float)(cos(phip)*sin(azp) + cos(altp)*cos(azp)*sin(phip));
02001 //              matrix[0][2] =  (float)(sin(altp)*sin(phip));
02002 //              matrix[1][0] =  (float)(-sin(phip)*cos(azp) - cos(altp)*sin(azp)*cos(phip));
02003 //              matrix[1][1] =  (float)(-sin(phip)*sin(azp) + cos(altp)*cos(azp)*cos(phip));
02004 //              matrix[1][2] =  (float)(sin(altp)*cos(phip));
02005 //              matrix[2][0] =  (float)(sin(altp)*sin(azp));
02006 //              matrix[2][1] =  (float)(-sin(altp)*cos(azp));
02007 //              matrix[2][2] =  (float)cos(altp);
02008 //      }
02009 //      if (is_quaternion){
02010 //              matrix[0][0] = (float)(e0 * e0 + e1 * e1 - e2 * e2 - e3 * e3);
02011 //              matrix[0][1] = (float)(2.0f * (e1 * e2 + e0 * e3));
02012 //              matrix[0][2] = (float)(2.0f * (e1 * e3 - e0 * e2));
02013 //              matrix[1][0] = (float)(2.0f * (e2 * e1 - e0 * e3));
02014 //              matrix[1][1] = (float)(e0 * e0 - e1 * e1 + e2 * e2 - e3 * e3);
02015 //              matrix[1][2] = (float)(2.0f * (e2 * e3 + e0 * e1));
02016 //              matrix[2][0] = (float)(2.0f * (e3 * e1 + e0 * e2));
02017 //              matrix[2][1] = (float)(2.0f * (e3 * e2 - e0 * e1));
02018 //              matrix[2][2] = (float)(e0 * e0 - e1 * e1 - e2 * e2 + e3 * e3);
02019 //              // keep in mind matrix[0][2] is M13 gives an e0 e2 piece, etc
02020 //      }
02021 //      // Now do post and pretrans: vfinal = vpost + R vpre;
02022 //
02023 //      matrix[0][3] = postT[0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2]  ;
02024 //      matrix[1][3] = postT[1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2]  ;
02025 //      matrix[2][3] = postT[2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2]  ;
02026 //}
02027 //
02028 //
02029 //void Transform3D::set_rotation(const float& m11, const float& m12, const float& m13,
02030 //                                                         const float& m21, const float& m22, const float& m23,
02031 //                const float& m31, const float& m32, const float& m33)
02032 //{
02033 //      EulerType euler_type = MATRIX;
02034 //      Dict rot;
02035 //      rot["m11"]  = m11;
02036 //      rot["m12"]  = m12;
02037 //      rot["m13"]  = m13;
02038 //      rot["m21"]  = m21;
02039 //      rot["m22"]  = m22;
02040 //      rot["m23"]  = m23;
02041 //      rot["m31"]  = m31;
02042 //      rot["m32"]  = m32;
02043 //      rot["m33"]  = m33;
02044 //      set_rotation(euler_type, rot);  // Or should it be &rot ?
02045 //}
02046 //
02047 //void Transform3D::set_rotation(const Vec3f & eahat, const Vec3f & ebhat,
02048 //                                    const Vec3f & eAhat, const Vec3f & eBhat)
02049 //{// this rotation rotates unit vectors a,b into A,B;
02051 //      Vec3f eahatcp(eahat);
02052 //      Vec3f ebhatcp(ebhat);
02053 //      Vec3f eAhatcp(eAhat);
02054 //      Vec3f eBhatcp(eBhat);
02055 //
02056 //      eahatcp.normalize();
02057 //      ebhatcp.normalize();
02058 //      eAhatcp.normalize();
02059 //      eBhatcp.normalize();
02060 //
02061 //      Vec3f aMinusA(eahatcp);
02062 //      aMinusA  -= eAhatcp;
02063 //      Vec3f bMinusB(ebhatcp);
02064 //      bMinusB  -= eBhatcp;
02065 //
02066 //
02067 //      Vec3f  nhat;
02068 //      float aAlength = aMinusA.length();
02069 //      float bBlength = bMinusB.length();
02070 //      if (aAlength==0){
02071 //              nhat=eahatcp;
02072 //      }else if (bBlength==0){
02073 //              nhat=ebhatcp;
02074 //      }else{
02075 //              nhat= aMinusA.cross(bMinusB);
02076 //              nhat.normalize();
02077 //      }
02078 //
02080 //
02081 //      Vec3f neahat  = eahatcp.cross(nhat);
02082 //      Vec3f nebhat  = ebhatcp.cross(nhat);
02083 //      Vec3f neAhat  = eAhatcp.cross(nhat);
02084 //      Vec3f neBhat  = eBhatcp.cross(nhat);
02085 //
02086 //      double cosOmegaA = (neahat.dot(neAhat))  / (neahat.dot(neahat));
02088 //      double sinOmegaA = (neahat.dot(eAhatcp)) / (neahat.dot(neahat));
02090 //
02091 //      double OmegaA = atan2(sinOmegaA,cosOmegaA);
02093 //
02094 //      EulerType euler_type=SPIN;
02095 //      Dict rotation;
02096 //      rotation["n1"]= nhat[0];
02097 //      rotation["n2"]= nhat[1];
02098 //      rotation["n3"]= nhat[2];
02099 //      rotation["Omega"] = (float)(OmegaA*180.0/M_PI);
02100 //      set_rotation(euler_type,  rotation);
02101 //}
02102 //
02103 //
02104 //float Transform3D::get_scale() const     // YYY
02105 //{
02106 //      // Assumes uniform scaling, calculation uses Z only.
02107 //      float scale =0;
02108 //      for (int i=0; i<3; i++) {
02109 //              for (int j=0; j<3; j++) {
02110 //                      scale = scale + matrix[i][j]*matrix[i][j];
02111 //              }
02112 //      }
02113 //
02114 //      return sqrt(scale/3);
02115 //}
02116 //
02117 //
02118 //
02119 //Dict Transform3D::get_rotation(EulerType euler_type) const
02120 //{
02121 //      Dict result;
02122 //
02123 //      double max = 1 - ERR_LIMIT;
02124 //      double sca=get_scale();
02125 //      double cosalt=matrix[2][2]/sca;
02126 //
02127 //
02128 //      double az=0;
02129 //      double alt = 0;
02130 //      double phi=0;
02131 //      double phiS = 0; // like az     (but in SPIDER ZXZ)
02132 //      double psiS =0;  // like phi  (but in SPIDER ZYZ)
02133 //
02134 //
02136 //
02137 //      if (cosalt > max) {  // that is, alt close to 0
02138 //              alt = 0;
02139 //              az=0;
02140 //              phi = (double)EMConsts::rad2deg * atan2(matrix[0][1], matrix[0][0]);
02141 //      }
02142 //      else if (cosalt < -max) { // alt close to pi
02143 //              alt = 180;
02144 //              az=0;
02145 //              phi=360.0f-(double)EMConsts::rad2deg * atan2(matrix[0][1], matrix[0][0]);
02146 //      }
02147 //      else {
02148 //              alt = (double)EMConsts::rad2deg * acos(cosalt);
02149 //              az  = 360.0f+(double)EMConsts::rad2deg * atan2(matrix[2][0], -matrix[2][1]);
02150 //              phi = 360.0f+(double)EMConsts::rad2deg * atan2(matrix[0][2], matrix[1][2]);
02151 //      }
02152 //      az =fmod(az+180.0,360.0)-180.0;
02153 //      phi=fmod(phi+180.0,360.0)-180.0;
02154 //
02156 //      if (fabs(cosalt) > max) {  // that is, alt close to 0
02157 //              phiS=0;
02158 //              psiS = az+phi;
02159 //      }
02160 //      else {
02161 //              phiS = az   - 90.0;
02162 //              psiS = phi  + 90.0;
02163 //      }
02164 //      phiS = fmod((phiS   + 360.0 ), 360.0) ;
02165 //      psiS = fmod((psiS   + 360.0 ), 360.0) ;
02166 //
02168 //
02169 //      double nphi = (az-phi)/2.0;
02170 //    // The next is also e0
02171 //      double cosOover2 = (cos((az+phi)*M_PI/360) * cos(alt*M_PI/360)) ;
02172 //      double sinOover2 = sqrt(1 -cosOover2*cosOover2);
02173 //      double cosnTheta = sin((az+phi)*M_PI/360) * cos(alt*M_PI/360) / sqrt(1-cosOover2*cosOover2) ;
02174 //      double sinnTheta = sqrt(1-cosnTheta*cosnTheta);
02175 //      double n1 = sinnTheta*cos(nphi*M_PI/180);
02176 //      double n2 = sinnTheta*sin(nphi*M_PI/180);
02177 //      double n3 = cosnTheta;
02178 //      double xtilt = 0;
02179 //      double ytilt = 0;
02180 //      double ztilt = 0;
02181 //
02182 //
02183 //      if (cosOover2<0) {
02184 //              cosOover2*=-1; n1 *=-1; n2*=-1; n3*=-1;
02185 //      }
02186 //
02187 //
02188 //      switch (euler_type) {
02189 //      case EMAN:
02190 //              result["az"]  = az;
02191 //              result["alt"] = alt;
02192 //              result["phi"] = phi;
02193 //              break;
02194 //
02195 //      case IMAGIC:
02196 //              result["alpha"] = az;
02197 //              result["beta"] = alt;
02198 //              result["gamma"] = phi;
02199 //              break;
02200 //
02201 //      case SPIDER:
02202 //              result["phi"]   = phiS;  // The first Euler like az
02203 //              result["theta"] = alt;
02204 //              result["psi"]   = psiS;
02205 //              break;
02206 //
02207 //      case MRC:
02208 //              result["phi"]   = phiS;
02209 //              result["theta"] = alt;
02210 //              result["omega"] = psiS;
02211 //              break;
02212 //
02213 //      case XYZ:
02214 //              xtilt = atan2(-sin((M_PI/180.0f)*phiS)*sin((M_PI/180.0f)*alt),cos((M_PI/180.0f)*alt));
02215 //              ytilt = asin(  cos((M_PI/180.0f)*phiS)*sin((M_PI/180.0f)*alt));
02216 //              ztilt = psiS*M_PI/180.0f - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
02217 //
02218 //              xtilt=fmod(xtilt*180/M_PI+540.0,360.0) -180.0;
02219 //              ztilt=fmod(ztilt*180/M_PI+540.0,360.0) -180.0;
02220 //
02221 //              result["xtilt"]  = xtilt;
02222 //              result["ytilt"]  = ytilt*180/M_PI;
02223 //              result["ztilt"]  = ztilt;
02224 //              break;
02225 //
02226 //      case QUATERNION:
02227 //              result["e0"] = cosOover2 ;
02228 //              result["e1"] = sinOover2 * n1 ;
02229 //              result["e2"] = sinOover2 * n2;
02230 //              result["e3"] = sinOover2 * n3;
02231 //              break;
02232 //
02233 //      case SPIN:
02234 //              result["Omega"] =360.0f* acos(cosOover2)/ M_PI ;
02235 //              result["n1"] = n1;
02236 //              result["n2"] = n2;
02237 //              result["n3"] = n3;
02238 //              break;
02239 //
02240 //      case SGIROT:
02241 //              result["q"] = 360.0f*acos(cosOover2)/M_PI ;
02242 //              result["n1"] = n1;
02243 //              result["n2"] = n2;
02244 //              result["n3"] = n3;
02245 //              break;
02246 //
02247 //      case MATRIX:
02248 //              result["m11"] = matrix[0][0] ;
02249 //              result["m12"] = matrix[0][1] ;
02250 //              result["m13"] = matrix[0][2] ;
02251 //              result["m21"] = matrix[1][0] ;
02252 //              result["m22"] = matrix[1][1] ;
02253 //              result["m23"] = matrix[1][2] ;
02254 //              result["m31"] = matrix[2][0] ;
02255 //              result["m32"] = matrix[2][1] ;
02256 //              result["m33"] = matrix[2][2] ;
02257 //              break;
02258 //
02259 //      default:
02260 //              throw InvalidValueException(euler_type, "unknown Euler Type");
02261 //      }
02262 //
02263 //      return result;
02264 //}
02265 //
02266 //Transform3D Transform3D::inverseUsingAngs() const    //   YYN need to test it for sure
02267 //{
02268 //      // First Find the scale
02269 //      EulerType eE=EMAN;
02270 //
02271 //
02272 //      float scale   = get_scale();
02273 //      Vec3f preT   = get_pretrans( ) ;
02274 //      Vec3f postT   = get_posttrans( ) ;
02275 //      Dict angs     = get_rotation(eE);
02276 //      Dict invAngs  ;
02277 //
02278 //      invAngs["phi"]   = 180.0f - (float) angs["az"] ;
02279 //      invAngs["az"]    = 180.0f - (float) angs["phi"] ;
02280 //      invAngs["alt"]   = angs["alt"] ;
02281 //
02288 //
02289 //      float inverseScale= 1/scale ;
02290 //
02291 //      Transform3D invM;
02292 //
02293 //      invM.set_rotation(EMAN, invAngs);
02294 //      invM.apply_scale(inverseScale);
02295 //      invM.set_pretrans(-postT );
02296 //      invM.set_posttrans(-preT );
02297 //
02298 //
02299 //      return invM;
02300 //
02301 //}
02302 //
02303 //Transform3D Transform3D::inverse() const    //   YYN need to test it for sure
02304 //{
02305 //      // This assumes the matrix is 4 by 4 and the last row reads [0 0 0 1]
02306 //
02307 //      double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
02308 //      double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
02309 //      double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
02310 //      double v0  = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
02311 //
02312 //    double cof00 = m11*m22-m12*m21;
02313 //    double cof11 = m22*m00-m20*m02;
02314 //    double cof22 = m00*m11-m01*m10;
02315 //    double cof01 = m10*m22-m20*m12;
02316 //    double cof02 = m10*m21-m20*m11;
02317 //    double cof12 = m00*m21-m01*m20;
02318 //    double cof10 = m01*m22-m02*m21;
02319 //    double cof20 = m01*m12-m02*m11;
02320 //    double cof21 = m00*m12-m10*m02;
02321 //
02322 //    double Det = m00* cof00 + m02* cof02 -m01*cof01;
02323 //
02324 //    Transform3D invM;
02325 //
02326 //    invM.matrix[0][0] =  (float)(cof00/Det);
02327 //    invM.matrix[0][1] =  (float)(-cof10/Det);
02328 //    invM.matrix[0][2] =  (float)(cof20/Det);
02329 //    invM.matrix[1][0] =  (float)(-cof01/Det);
02330 //    invM.matrix[1][1] =  (float)(cof11/Det);
02331 //    invM.matrix[1][2] =  (float)(-cof21/Det);
02332 //    invM.matrix[2][0] =  (float)(cof02/Det);
02333 //    invM.matrix[2][1] =  (float)(-cof12/Det);
02334 //    invM.matrix[2][2] =  (float)(cof22/Det);
02335 //
02336 //    invM.matrix[0][3] =  (float)((-cof00*v0 + cof10*v1 - cof20*v2)/Det);
02337 //    invM.matrix[1][3] =  (float)(( cof01*v0 - cof11*v1 + cof21*v2)/Det);
02338 //    invM.matrix[2][3] =  (float)((-cof02*v0 + cof12*v1 - cof22*v2)/Det);
02339 //
02340 //      Vec3f postT   = get_posttrans( ) ;
02341 //      Vec3f invMpre   = - postT;
02342 //      Vec3f invMpost   ;
02343 //      for ( int i = 0; i < 3; i++) {
02344 //              invMpost[i] = invM.matrix[i][3];
02345 //              for ( int j = 0; j < 3; j++) {
02346 //                      invMpost[i] += - invM.matrix[i][j]*invMpre[j];
02347 //              }
02348 //              invM.matrix[3][i] = invMpost[i];
02349 //      }
02350 //
02351 //      return invM;
02352 //}
02353 //
02354 //
02355 //
02357 //
02358 //Transform3D Transform3D::get_sym(const string & symname, int n) const
02359 //{
02360 //      int nsym = get_nsym(symname);
02361 //
02364 //
02365 //      // see www.math.utah.edu/~alfeld/math/polyhedra/polyhedra.html for pictures
02366 //      // By default we will put largest symmetry along z-axis.
02367 //
02368 //      // Each Platonic Solid has 2E symmetry elements.
02369 //
02370 //
02371 //      // An icosahedron has   m=5, n=3, F=20 E=30=nF/2, V=12=nF/m,since vertices shared by 5 triangles;
02372 //      // It is composed of 20 triangles. E=3*20/2;
02373 //
02374 //
02375 //      // An dodecahedron has m=3, n=5   F=12 E=30  V=20
02376 //      // It is composed of 12 pentagons. E=5*12/2;   V= 5*12/3, since vertices shared by 3 pentagons;
02377 //
02378 //
02379 //
02380 //    // The ICOS symmetry group has the face along z-axis
02381 //
02382 //      float lvl0=0;                             //  there is one pentagon on top; five-fold along z
02383 //      float lvl1= 63.4349f; // that is atan(2)  // there are 5 pentagons with centers at this height (angle)
02384 //      float lvl2=116.5651f; //that is 180-lvl1  // there are 5 pentagons with centers at this height (angle)
02385 //      float lvl3=180.f;                           // there is one pentagon on the bottom
02386 //             // Notice that 63.439 is the angle between two faces of the dual object
02387 //
02388 //      static double ICOS[180] = { // This is with a pentagon normal to z
02389 //                0,lvl0,0,    0,lvl0,288,   0,lvl0,216,   0,lvl0,144,  0,lvl0,72,
02390 //                0,lvl1,36,   0,lvl1,324,   0,lvl1,252,   0,lvl1,180,  0,lvl1,108,
02391 //               72,lvl1,36,  72,lvl1,324,  72,lvl1,252,  72,lvl1,180,  72,lvl1,108,
02392 //              144,lvl1,36, 144,lvl1,324, 144,lvl1,252, 144,lvl1,180, 144,lvl1,108,
02393 //              216,lvl1,36, 216,lvl1,324, 216,lvl1,252, 216,lvl1,180, 216,lvl1,108,
02394 //              288,lvl1,36, 288,lvl1,324, 288,lvl1,252, 288,lvl1,180, 288,lvl1,108,
02395 //               36,lvl2,0,   36,lvl2,288,  36,lvl2,216,  36,lvl2,144,  36,lvl2,72,
02396 //              108,lvl2,0,  108,lvl2,288, 108,lvl2,216, 108,lvl2,144, 108,lvl2,72,
02397 //              180,lvl2,0,  180,lvl2,288, 180,lvl2,216, 180,lvl2,144, 180,lvl2,72,
02398 //              252,lvl2,0,  252,lvl2,288, 252,lvl2,216, 252,lvl2,144, 252,lvl2,72,
02399 //              324,lvl2,0,  324,lvl2,288, 324,lvl2,216, 324,lvl2,144, 324,lvl2,72,
02400 //                0,lvl3,0,    0,lvl3,288,   0,lvl3,216,   0,lvl3,144,   0,lvl3,72
02401 //      };
02402 //
02403 //
02404 //      // A cube has   m=3, n=4, F=6 E=12=nF/2, V=8=nF/m,since vertices shared by 3 squares;
02405 //      // It is composed of 6 squares.
02406 //
02407 //
02408 //      // An octahedron has   m=4, n=3, F=8 E=12=nF/2, V=6=nF/m,since vertices shared by 4 triangles;
02409 //      // It is composed of 8 triangles.
02410 //
02411 //    // We have placed the OCT symmetry group with a face along the z-axis
02412 //        lvl0=0;
02413 //      lvl1=90;
02414 //      lvl2=180;
02415 //
02416 //      static float OCT[72] = {// This is with a face of a cube along z
02417 //                    0,lvl0,0,   0,lvl0,90,    0,lvl0,180,    0,lvl0,270,
02418 //                    0,lvl1,0,   0,lvl1,90,    0,lvl1,180,    0,lvl1,270,
02419 //                   90,lvl1,0,  90,lvl1,90,   90,lvl1,180,   90,lvl1,270,
02420 //                  180,lvl1,0, 180,lvl1,90,  180,lvl1,180,  180,lvl1,270,
02421 //                  270,lvl1,0, 270,lvl1,90,  270,lvl1,180,  270,lvl1,270,
02422 //                    0,lvl2,0,   0,lvl2,90,    0,lvl2,180,    0,lvl2,270
02423 //      };
02424 //      // B^4=A^3=1;  BABA=1; implies   AA=BAB, ABA=B^3 , AB^2A = BBBABBB and
02425 //      //   20 words with at most a single A
02426 //    //   1 B BB BBB A  BA AB BBA BAB ABB BBBA BBAB BABB ABBB BBBAB BBABB BABBB
02427 //    //                        BBBABB BBABBB BBBABBB
02428 //     // also     ABBBA is distinct yields 4 more words
02429 //     //    ABBBA   BABBBA BBABBBA BBBABBBA
02430 //     // for a total of 24 words
02431 //     // Note A BBB A BBB A  reduces to BBABB
02432 //     //  and  B A BBB A is the same as A BBB A BBB etc.
02433 //
02434 //    // The TET symmetry group has a face along the z-axis
02435 //    // It has n=m=3; F=4, E=6=nF/2, V=4=nF/m
02436 //        lvl0=0;         // There is a face along z
02437 //      lvl1=109.4712f;  //  that is acos(-1/3)  // There  are 3 faces at this angle
02438 //
02439 //      static float TET[36] = {// This is with the face along z
02440 //            0,lvl0,0,   0,lvl0,120,    0,lvl0,240,
02441 //            0,lvl1,60,   0,lvl1,180,    0,lvl1,300,
02442 //          120,lvl1,60, 120,lvl1,180,  120,lvl1,300,
02443 //          240,lvl1,60, 240,lvl1,180,  240,lvl1,300
02444 //      };
02445 //      // B^3=A^3=1;  BABA=1; implies   A^2=BAB, ABA=B^2 , AB^2A = B^2AB^2 and
02446 //      //   12 words with at most a single A
02447 //    //   1 B BB  A  BA AB BBA BAB ABB BBAB BABB BBABB
02448 //    // at most one A is necessary
02449 //
02450 //      Transform3D ret;
02451 //      SymType type = get_sym_type(symname);
02452 //
02453 //      switch (type) {
02454 //      case CSYM:
02455 //              ret.set_rotation( n * 360.0f / nsym, 0, 0);
02456 //              break;
02457 //      case DSYM:
02458 //              if (n >= nsym / 2) {
02459 //                      ret.set_rotation((n - nsym/2) * 360.0f / (nsym / 2),180.0f, 0);
02460 //              }
02461 //              else {
02462 //                      ret.set_rotation( n * 360.0f / (nsym / 2),0, 0);
02463 //              }
02464 //              break;
02465 //      case ICOS_SYM:
02466 //              ret.set_rotation((float)ICOS[n * 3 ],
02467 //                               (float)ICOS[n * 3 + 1],
02468 //                               (float)ICOS[n * 3 + 2] );
02469 //              break;
02470 //      case OCT_SYM:
02471 //              ret.set_rotation((float)OCT[n * 3],
02472 //                               (float)OCT[n * 3 + 1],
02473 //                               (float)OCT[n * 3 + 2] );
02474 //              break;
02475 //      case TET_SYM:
02476 //              ret.set_rotation((float)TET[n * 3 ],
02477 //                               (float)TET[n * 3 + 1] ,
02478 //                               (float)TET[n * 3 + 2] );
02479 //              break;
02480 //      case ISYM:
02481 //              ret.set_rotation(0, 0, 0);
02482 //              break;
02483 //      default:
02484 //              throw InvalidValueException(type, symname);
02485 //      }
02486 //
02487 //      ret = (*this) * ret;
02488 //
02489 //      return ret;
02490 //}
02491 //
02492 //int Transform3D::get_nsym(const string & name)
02493 //{
02494 //      string symname = name;
02495 //
02496 //      for (size_t i = 0; i < name.size(); i++) {
02497 //              if (isalpha(name[i])) {
02498 //                      symname[i] = (char)tolower(name[i]);
02499 //              }
02500 //      }
02501 //
02502 //      SymType type = get_sym_type(symname);
02503 //      int nsym = 0;
02504 //
02505 //      switch (type) {
02506 //      case CSYM:
02507 //              nsym = atoi(symname.c_str() + 1);
02508 //              break;
02509 //      case DSYM:
02510 //              nsym = atoi(symname.c_str() + 1) * 2;
02511 //              break;
02512 //      case ICOS_SYM:
02513 //              nsym = 60;
02514 //              break;
02515 //      case OCT_SYM:
02516 //              nsym = 24;
02517 //              break;
02518 //      case TET_SYM:
02519 //              nsym = 12;
02520 //              break;
02521 //      case ISYM:
02522 //              nsym = 1;
02523 //              break;
02524 //      case UNKNOWN_SYM:
02525 //      default:
02526 //              throw InvalidValueException(type, name);
02527 //      }
02528 //      return nsym;
02529 //}
02530 //
02531 //
02532 //
02533 //Transform3D::SymType Transform3D::get_sym_type(const string & name)
02534 //{
02535 //      SymType t = UNKNOWN_SYM;
02536 //
02537 //      if (name[0] == 'c') {
02538 //              t = CSYM;
02539 //      }
02540 //      else if (name[0] == 'd') {
02541 //              t = DSYM;
02542 //      }
02543 //      else if (name == "icos") {
02544 //              t = ICOS_SYM;
02545 //      }
02546 //      else if (name == "oct") {
02547 //              t = OCT_SYM;
02548 //      }
02549 //      else if (name == "tet") {
02550 //              t = TET_SYM;
02551 //      }
02552 //      else if (name == "i" || name == "") {
02553 //              t = ISYM;
02554 //      }
02555 //      return t;
02556 //}
02557 //
02558 //vector<Transform3D*>
02559 //Transform3D::angles2tfvec(EulerType eulertype, const vector<float> ang) {
02560 //      int nangles = ang.size() / 3;
02561 //      vector<Transform3D*> tfvec;
02562 //      for (int i = 0; i < nangles; i++) {
02563 //              tfvec.push_back(new Transform3D(eulertype,ang[3*i],ang[3*i+1],ang[3*i+2]));
02564 //      }
02565 //      return tfvec;
02566 //}
02567 //
02568 
02569 
02570 
02571 /* vim: set ts=4 noet: */
02572 
02573 
02574 /*    Rotation stuff */
02575 

Generated on Tue Jul 12 13:49:04 2011 for EMAN2 by  doxygen 1.3.9.1