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