EMAN2
transform.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include "transform.h"
33#include "util.h"
34#include "emobject.h"
35#include <cctype> // for std::tolower
36#include <cstring> // for memcpy
37#include "symmetry.h"
38using namespace EMAN;
39
40#ifdef WIN32
41#ifndef M_PI
42#define M_PI 3.14159265358979323846
43#endif
44#endif
45
46#include <algorithm> // for std::transform
47
48#include <gsl/gsl_matrix.h>
49#include <gsl/gsl_blas.h>
50#include <gsl/gsl_linalg.h>
51
52#include <ostream>
53using std::ostream_iterator;
54
55//const float Transform3D::ERR_LIMIT = 0.000001f;
56
57const float Transform::ERR_LIMIT = 0.000001f;
58
59vector<string> Transform::permissable_2d_not_rot;
60vector<string> Transform::permissable_3d_not_rot;
61map<string,vector<string> > Transform::permissable_rot_keys;
62
63Transform Transform::icos_5_to_2() {
64 Transform t;
65 Dict d;
66 d["type"] = "eman";
67 d["phi"] = 0;
68 d["az"] = 90.0f;
69 d["alt"] = 31.717474; // 5 fold to the nearest 2-fold
81 t.set_rotation(d);
82 return t;
83}
84
85Transform Transform::tet_3_to_2() {
86 Transform t;
87 Dict d;
88 d["type"] = "eman";
89 d["phi"] = 45.0f;
90 d["az"] = 0.0f;
91 d["alt"] = 54.73561f; // 3 fold to a 2 fold
95 t.set_rotation(d);
96 return t;
97}
98
99
100Transform::Transform()
101{
102 to_identity();
103}
104
105Transform::Transform( const Transform& that )
106{
107 *this = that;
108}
109
110Transform& Transform::operator=(const Transform& that ) {
111 if (this != &that ) {
112// for (int i=0; i<3; i++)
113// for (int j=0; j<4; j++) matrix[i][j]=that.matrix[i][j];
114 memcpy(matrix,that.matrix,12*sizeof(float));
115 }
116 return *this;
117}
118
119bool Transform::operator==(const Transform& rhs) const{
120 if (memcmp(this->matrix, rhs.matrix, 3*4*sizeof(float)) == 0) {
121 return true;
122 }
123 else {
124 return false;
125 }
126}
127
128bool Transform::operator!=(const Transform& rhs) const{
129 return !(operator==(rhs));
130}
131
132Transform::Transform(const Dict& d) {
133 to_identity();
134 set_params(d);
135}
136
137
138Transform::Transform(const float array[12]) {
139 memcpy(matrix,array,12*sizeof(float));
140}
141
142Transform::Transform(const vector<float> array)
143{
144 set_matrix(array);
145}
146
147void Transform::set_matrix(const vector<float>& v)
148{
149 if (v.size() != 12 ) throw InvalidParameterException("The construction array must be of size 12");
150
151 for(int i=0; i<3; ++i) {
152 for(int j=0; j<4; ++j) {
153 matrix[i][j] = v[i*4+j];
154 }
155 }
156}
157
158void Transform::copy_matrix_into_array(float* const array) const {
159
160 int idx = 0;
161 for(int i=0; i<3; ++i) {
162 for(int j=0; j<4; ++j) {
163 array[idx] = matrix[i][j];
164 idx ++;
165 }
166 }
167}
168
169string Transform::get_matrix_string(int precision) {
170 if (precision<2||precision>9) precision=9;
171 char fmt[256];
172 char buf[256]; // sufficient for any potential result
173 sprintf(fmt,"[%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg,%%1.%0dg]",
174 precision,precision,precision,precision,precision,precision,precision,precision,precision,precision,precision,precision);
175 sprintf(buf,fmt,matrix[0][0],matrix[0][1],matrix[0][2],matrix[0][3],matrix[1][0],matrix[1][1],matrix[1][2],matrix[1][3],
176 matrix[2][0],matrix[2][1],matrix[2][2],matrix[2][3]);
177 return string(buf);
178}
179
180
181vector<float> Transform::get_matrix() const
182{
183 vector<float> ret(12);
184 for(int i=0; i<3; ++i) {
185 for(int j=0; j<4; ++j) {
186 ret[i*4+j] = matrix[i][j];
187 }
188 }
189 return ret;
190}
191
192vector<float> Transform::get_matrix_4x4() const
193{
194 vector<float> ret(16);
195 for(int i=0; i<3; ++i) {
196 for(int j=0; j<4; ++j) {
197 ret[i*4+j] = matrix[i][j];
198 }
199 }
200 ret[12] = 0.0;
201 ret[13] = 0.0;
202 ret[14] = 0.0;
203 ret[15] = 1.0;
204
205 return ret;
206}
207void Transform::to_identity()
208{
209// transform_type = UNKNOWN;
210 for(int i=0; i<3; ++i) {
211 for(int j=0; j<4; ++j) {
212 if(i==j) {
213 matrix[i][j] = 1;
214 }
215 else {
216 matrix[i][j] = 0;
217 }
218 }
219 }
220}
221
222bool Transform::is_identity() const {
223 for(int i=0; i<3; ++i) {
224 for(int j=0; j<4; ++j) {
225 float c = matrix[i][j];
226 Util::apply_precision(c,ERR_LIMIT);
227 if(i==j) {
228 if (c != 1.0) return false;
229 }
230 else {
231 if (c != 0.0) return false;
232 }
233 }
234 }
235 return true;
236}
237
238bool Transform::is_rot_identity() const {
239 for(int i=0; i<3; ++i) {
240 for(int j=0; j<3; ++j) {
241 float c = matrix[i][j];
242 Util::apply_precision(c,ERR_LIMIT);
243 if(i==j) {
244 if (c != 1.0) return false;
245 }
246 else {
247 if (c != 0.0) return false;
248 }
249 }
250 }
251 return true;
252}
253
254void Transform::set_params(const Dict& d) {
255 detect_problem_keys(d);
256
257 if (d.has_key_ci("type") ) set_rotation(d);
258
259 if (d.has_key_ci("scale")) {
260 float scale = static_cast<float>(d.get_ci("scale"));
261 set_scale(scale);
262 }
263
264 float dx=0,dy=0,dz=0;
265
266 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
267 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
268 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
269
270 if ( dx != 0.0 || dy != 0.0 || dz != 0.0 ) {
271 set_trans(dx,dy,dz);
272 }
273
274 if (d.has_key_ci("mirror")) {
275 EMObject e = d.get_ci("mirror");
276 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
277 throw InvalidParameterException("Error, mirror must be a bool or an int");
278
279 bool mirror = static_cast<bool>(e);
280 set_mirror(mirror);
281 }
282}
283
284
285void Transform::init_permissable_keys()
286{
287
288 permissable_2d_not_rot.push_back("tx");
289 permissable_2d_not_rot.push_back("ty");
290 permissable_2d_not_rot.push_back("scale");
291 permissable_2d_not_rot.push_back("mirror");
292 permissable_2d_not_rot.push_back("type");
293
294 permissable_3d_not_rot.push_back("tx");
295 permissable_3d_not_rot.push_back("ty");
296 permissable_3d_not_rot.push_back("tz");
297 permissable_3d_not_rot.push_back("scale");
298 permissable_3d_not_rot.push_back("mirror");
299 permissable_3d_not_rot.push_back("type");
300
301 vector<string> tmp;
302 tmp.push_back("alpha");
303 permissable_rot_keys["2d"] = tmp;
304
305 tmp.clear();
306 tmp.push_back("alt");
307 tmp.push_back("az");
308 tmp.push_back("phi");
309 permissable_rot_keys["eman"] = tmp;
310
311 tmp.clear();
312 tmp.push_back("psi");
313 tmp.push_back("theta");
314 tmp.push_back("phi");
315 permissable_rot_keys["spider"] = tmp;
316
317 tmp.clear();
318 tmp.push_back("alpha");
319 tmp.push_back("beta");
320 tmp.push_back("gamma");
321 permissable_rot_keys["imagic"] = tmp;
322
323 tmp.clear();
324 tmp.push_back("ztilt");
325 tmp.push_back("xtilt");
326 tmp.push_back("ytilt");
327 permissable_rot_keys["xyz"] = tmp;
328
329 tmp.clear();
330 tmp.push_back("phi");
331 tmp.push_back("theta");
332 tmp.push_back("omega");
333 permissable_rot_keys["mrc"] = tmp;
334
335 tmp.clear();
336 tmp.push_back("e0");
337 tmp.push_back("e1");
338 tmp.push_back("e2");
339 tmp.push_back("e3");
340 permissable_rot_keys["quaternion"] = tmp;
341
342 tmp.clear();
343 tmp.push_back("n1");
344 tmp.push_back("n2");
345 tmp.push_back("n3");
346 tmp.push_back("omega");
347 permissable_rot_keys["spin"] = tmp;
348
349 tmp.clear();
350 tmp.push_back("n1");
351 tmp.push_back("n2");
352 tmp.push_back("n3");
353 tmp.push_back("q");
354 permissable_rot_keys["sgirot"] = tmp;
355
356 tmp.clear();
357 tmp.push_back("m11");
358 tmp.push_back("m12");
359 tmp.push_back("m13");
360 tmp.push_back("m21");
361 tmp.push_back("m22");
362 tmp.push_back("m23");
363 tmp.push_back("m31");
364 tmp.push_back("m32");
365 tmp.push_back("m33");
366 permissable_rot_keys["matrix"] = tmp;
367}
368
369void Transform::detect_problem_keys(const Dict& d) {
370 if (permissable_rot_keys.size() == 0 ) {
371 init_permissable_keys();
372 }
373
374 vector<string> verification;
375 vector<string> problem_keys;
376 bool is_2d = false;
377 if (d.has_key_ci("type") ) {
378 string type = Util::str_to_lower((string)d["type"]);
379 bool problem = false;
380 if (permissable_rot_keys.find(type) == permissable_rot_keys.end() ) {
381 problem_keys.push_back(type);
382 problem = true;
383 }
384 if ( !problem ) {
385 vector<string> perm = permissable_rot_keys[type];
386 std::copy(perm.begin(),perm.end(),back_inserter(verification));
387
388 if ( type == "2d" ) {
389 is_2d = true;
390 std::copy(permissable_2d_not_rot.begin(),permissable_2d_not_rot.end(),back_inserter(verification));
391 }
392 }
393 }
394 if ( !is_2d ) {
395 std::copy(permissable_3d_not_rot.begin(),permissable_3d_not_rot.end(),back_inserter(verification));
396 }
397
398 for (Dict::const_iterator it = d.begin(); it != d.end(); ++it) {
399 if ( std::find(verification.begin(),verification.end(), it->first) == verification.end() ) {
400 problem_keys.push_back(it->first);
401 }
402 }
403
404 if (problem_keys.size() != 0 ) {
405 string error;
406 if (problem_keys.size() == 1) {
407 error = "Transform Error: The \"" +problem_keys[0]+ "\" key is unsupported";
408 } else {
409 error = "Transform Error: The ";
410 for(vector<string>::const_iterator cit = problem_keys.begin(); cit != problem_keys.end(); ++cit ) {
411 if ( cit != problem_keys.begin() ) {
412 if (cit == (problem_keys.end() -1) ) error += " and ";
413 else error += ", ";
414 }
415 error += "\"";
416 error += *cit;
417 error += "\"";
418 }
419 error += " keys are unsupported";
420 }
421 throw InvalidParameterException(error);
422 }
423}
424
425void Transform::set_params_inverse(const Dict& d) {
426 detect_problem_keys(d);
427
428 if (d.has_key_ci("type") ) set_rotation(d);
429
430 float dx=0,dy=0,dz=0;
431 if (d.has_key_ci("tx")) dx = static_cast<float>(d.get_ci("tx"));
432 if (d.has_key_ci("ty")) dy = static_cast<float>(d.get_ci("ty"));
433 if (d.has_key_ci("tz")) dz = static_cast<float>(d.get_ci("tz"));
434
435 if ( (dx != 0.0 || dy != 0.0 || dz != 0.0) && d.has_key_ci("type") ) {
436 Transform pre_trans;
437 pre_trans.set_trans(dx,dy,dz);
438
439 Transform tmp;
440 tmp.set_rotation(d);
441
442 if (d.has_key_ci("scale")) {
443 float scale = static_cast<float>(d.get_ci("scale"));
444 tmp.set_scale(scale);
445 }
446
447 Transform solution_trans = tmp*pre_trans;
448
449 if (d.has_key_ci("scale")) {
450 Transform tmp;
451 float scale = static_cast<float>(d.get_ci("scale"));
452 tmp.set_scale(scale);
453 solution_trans = solution_trans*tmp;
454 }
455
456 tmp = Transform();
457 tmp.set_rotation(d);
458 solution_trans = solution_trans*tmp;
459 set_trans(solution_trans.get_trans());
460 }
461
462 if (d.has_key_ci("scale")) {
463 float scale = static_cast<float>(d.get_ci("scale"));
464 set_scale(scale);
465 }
466
467 if (d.has_key_ci("mirror")) {
468 EMObject e = d.get_ci("mirror");
469 if ( (e.get_type() != EMObject::BOOL ) && (e.get_type() != EMObject::INT ) && (e.get_type() != EMObject::UNSIGNEDINT ) )
470 throw InvalidParameterException("Error, mirror must be a bool or an int");
471
472 bool mirror = static_cast<bool>(e);
473 set_mirror(mirror);
474 }
475 invert();
476}
477
478
479Dict Transform::get_params(const string& euler_type) const {
480 Dict params = get_rotation(euler_type);
481
482 Vec3f v = get_trans();
483 params["tx"] = v[0]; params["ty"] = v[1];
484
485 string type = Util::str_to_lower(euler_type);
486 if ( type != "2d") params["tz"] = v[2];
487
488 float scale = get_scale();
489 params["scale"] = scale;
490
491 bool mirror = get_mirror();
492 params["mirror"] = mirror;
493
494 return params;
495}
496
497
498
499Dict Transform::get_params_inverse(const string& euler_type) const {
500 Transform inv(inverse());
501
502 Dict params = inv.get_rotation(euler_type);
503 Vec3f v = inv.get_pre_trans();
504 params["tx"] = v[0]; params["ty"] = v[1];
505
506 string type = Util::str_to_lower(euler_type);
507 if ( type != "2d") params["tz"] = v[2];
508
509 float scale = inv.get_scale();
510 params["scale"] = scale;
511
512 bool mirror = inv.get_mirror();
513 params["mirror"] = mirror;
514
515 return params;
516}
517
518
519void Transform::set_rotation(const Dict& rotation)
520{
521 detect_problem_keys(rotation);
522 string euler_type;
523
524 if (!rotation.has_key_ci("type") ){
525 throw InvalidParameterException("argument dictionary does not contain the type key");
526 }
527
528 euler_type = static_cast<string>(rotation.get_ci("type"));// Warning, will throw
529
530
531 double e0=0; double e1=0; double e2=0; double e3=0;
532 double omega=0;
533 double az = 0;
534 double alt = 0;
535 double phi = 0;
536 double cxtilt = 0;
537 double sxtilt = 0;
538 double cytilt = 0;
539 double sytilt = 0;
540 double cztilt = 0;
541 double sztilt = 0;
542 bool is_quaternion = 0;
543 bool is_matrix = 0;
544 bool is_xyz = 0;
545
546 bool x_mirror;
547 float scale;
548 // Get these before anything changes so we can apply them again after the rotation is set
549 get_scale_and_mirror(scale,x_mirror);
550 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
551
552 string type = Util::str_to_lower(euler_type);
553 if (type == "2d") {
554 assert_valid_2d();
555 az = 0;
556 alt = 0;
557 phi = (double)rotation["alpha"] ;
558 } else if ( type == "eman" ) {
559// validate_and_set_type(THREED);
560 az = (double)rotation["az"] ;
561 alt = (double)rotation["alt"] ;
562 phi = (double)rotation["phi"] ;
563 } else if ( type == "imagic" ) {
564// validate_and_set_type(THREED);
565 az = (double)rotation["alpha"] ;
566 alt = (double)rotation["beta"] ;
567 phi = (double)rotation["gamma"] ;
568 } else if ( type == "spider" ) {
569// validate_and_set_type(THREED);
570 az = (double)rotation["phi"] + 90.0;
571 alt = (double)rotation["theta"] ;
572 phi = (double)rotation["psi"] - 90.0;
573 } else if ( type == "xyz" ) {
574// validate_and_set_type(THREED);
575 is_xyz = 1;
576 cxtilt = cos(EMConsts::deg2rad*(double)rotation["xtilt"]);
577 sxtilt = sin(EMConsts::deg2rad*(double)rotation["xtilt"]);
578 cytilt = cos(EMConsts::deg2rad*(double)rotation["ytilt"]);
579 sytilt = sin(EMConsts::deg2rad*(double)rotation["ytilt"]);
580 cztilt = cos(EMConsts::deg2rad*(double)rotation["ztilt"]);
581 sztilt = sin(EMConsts::deg2rad*(double)rotation["ztilt"]);
582 } else if ( type == "mrc" ) {
583// validate_and_set_type(THREED);
584 az = (double)rotation["phi"] + 90.0f ;
585 alt = (double)rotation["theta"] ;
586 phi = (double)rotation["omega"] - 90.0f ;
587 } else if ( type == "quaternion" ) {
588// validate_and_set_type(THREED);
589 is_quaternion = 1;
590 e0 = (double)rotation["e0"];
591 e1 = (double)rotation["e1"];
592 e2 = (double)rotation["e2"];
593 e3 = (double)rotation["e3"];
594 } else if ( type == "spin" ) {
595// validate_and_set_type(THREED);
596 is_quaternion = 1;
597 omega = (double)rotation["omega"];
598 double norm=Util::hypot3((double)rotation["n1"],(double)rotation["n2"],(double)rotation["n3"]);
599 if (norm==0.0) {
600 e0=1.0;
601 e1=e2=e3=0.0;
602 } else {
603 e0 = cos(omega*EMConsts::deg2rad/2.0);
604 e1 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"]/norm;
605 e2 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"]/norm;
606 e3 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"]/norm;
607 }
608 } else if ( type == "sgirot" ) {
609// validate_and_set_type(THREED);
610 is_quaternion = 1;
611 omega = (double)rotation["q"] ;
612 e0 = cos(omega*EMConsts::deg2rad/2.0);
613 e1 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n1"];
614 e2 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n2"];
615 e3 = sin(omega*EMConsts::deg2rad/2.0) * (double)rotation["n3"];
616 } else if ( type == "matrix" ) {
617 is_matrix = 1;
618 matrix[0][0] = (float)rotation["m11"];
619 matrix[0][1] = (float)rotation["m12"];
620 matrix[0][2] = (float)rotation["m13"];
621 matrix[1][0] = (float)rotation["m21"];
622 matrix[1][1] = (float)rotation["m22"];
623 matrix[1][2] = (float)rotation["m23"];
624 matrix[2][0] = (float)rotation["m31"];
625 matrix[2][1] = (float)rotation["m32"];
626 matrix[2][2] = (float)rotation["m33"];
627 } else {
628// transform_type = UNKNOWN;
629 throw InvalidStringException(euler_type, "unknown Euler Type");
630 }
631
632 double azp = az*EMConsts::deg2rad;
633 double altp = alt*EMConsts::deg2rad;
634 double phip = phi*EMConsts::deg2rad;
635
636 if (!is_quaternion && !is_matrix && !is_xyz) {
637 matrix[0][0] = (float)(cos(phip)*cos(azp) - cos(altp)*sin(azp)*sin(phip));
638 matrix[0][1] = (float)(cos(phip)*sin(azp) + cos(altp)*cos(azp)*sin(phip));
639 matrix[0][2] = (float)(sin(altp)*sin(phip));
640 matrix[1][0] = (float)(-sin(phip)*cos(azp) - cos(altp)*sin(azp)*cos(phip));
641 matrix[1][1] = (float)(-sin(phip)*sin(azp) + cos(altp)*cos(azp)*cos(phip));
642 matrix[1][2] = (float)(sin(altp)*cos(phip));
643 matrix[2][0] = (float)(sin(altp)*sin(azp));
644 matrix[2][1] = (float)(-sin(altp)*cos(azp));
645 matrix[2][2] = (float)cos(altp);
646 }
647 if (is_quaternion){
648 matrix[0][0] = (float)(e0 * e0 + e1 * e1 - e2 * e2 - e3 * e3);
649 matrix[0][1] = (float)(2.0f * (e1 * e2 + e0 * e3));
650 matrix[0][2] = (float)(2.0f * (e1 * e3 - e0 * e2));
651 matrix[1][0] = (float)(2.0f * (e2 * e1 - e0 * e3));
652 matrix[1][1] = (float)(e0 * e0 - e1 * e1 + e2 * e2 - e3 * e3);
653 matrix[1][2] = (float)(2.0f * (e2 * e3 + e0 * e1));
654 matrix[2][0] = (float)(2.0f * (e3 * e1 + e0 * e2));
655 matrix[2][1] = (float)(2.0f * (e3 * e2 - e0 * e1));
656 matrix[2][2] = (float)(e0 * e0 - e1 * e1 - e2 * e2 + e3 * e3);
657 // keep in mind matrix[0][2] is M13 gives an e0 e2 piece, etc
658 }
659 if (is_xyz){
660 matrix[0][0] = (float)(cytilt*cztilt);
661 matrix[0][1] = (float)(cxtilt*sztilt+sxtilt*sytilt*cztilt);
662 matrix[0][2] = (float)(sxtilt*sztilt-cxtilt*sytilt*cztilt);
663 matrix[1][0] = (float)(-cytilt*sztilt);
664 matrix[1][1] = (float)(cxtilt*cztilt-sxtilt*sytilt*sztilt);
665 matrix[1][2] = (float)(sxtilt*cztilt+cxtilt*sytilt*sztilt);
666 matrix[2][0] = (float)(sytilt);
667 matrix[2][1] = (float)(-sxtilt*cytilt);
668 matrix[2][2] = (float)(cxtilt*cytilt);
669 }
670
671 // Apply scale if it existed previously
672 if (scale != 1.0f) {
673 for(int i=0; i<3; ++i) {
674 for(int j=0; j<3; ++j) {
675 matrix[i][j] *= scale;
676 }
677 }
678 }
679
680 // Apply post x mirroring if it was applied previously
681 if ( x_mirror ) {
682 for(int j=0; j<3; ++j) {
683 matrix[0][j] *= -1.0f;
684 }
685 }
686}
687
688Transform Transform::get_rotation_transform() const
689{
690 Transform ret(*this);
691 ret.set_scale(1.0);
692 ret.set_mirror(false);
693 ret.set_trans(0,0,0);
694 //ret.orthogonalize(); // ?
695 return ret;
696}
697
698void Transform::set_rotation(const Vec3f & v)
699{
700 if ( v[0] == 0 && v[1] == 0 && v[2] == 0 )
701 throw UnexpectedBehaviorException("Can't set rotation for the null vector");
702
703 Vec3f v1(v);
704 v1.normalize();
705
706 double theta = acos(v1[2]); // in radians
707 double psi = atan2(v1[1],-v1[0]);
708
709 Dict d;
710 d["theta"] = (double)EMConsts::rad2deg*theta;
711 d["psi"] = (double)EMConsts::rad2deg*psi;
712 d["phi"] = (double)0.0;
713 d["type"] = "spider";
714
715 set_rotation(d);
716
717
718}
719
720void Transform::rotate_origin(const Transform& by)
721{
722 vector<float> multmatrix = by.get_matrix();
723 // First Multiply and put the result in a temp matrix
724 Transform result;
725 for (int i=0; i<3; i++) {
726 for (int j=0; j<3; j++) {
727 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
728 }
729 }
730 //Then put the result from the tmep matrix in the original one
731 for (int i=0; i<3; i++) {
732 for (int j=0; j<3; j++) {
733 matrix[i][j] = result[i][j];
734 }
735 }
736}
737
738void Transform::rotate_origin_newBasis(const Transform& tcs, const float& omega, const float& n1, const float& n2, const float& n3)
739{
740 //Get the rotational inverse
741 Transform tcsinv = Transform(tcs);
742 tcsinv.set_trans(0.0, 0.0, 0.0);
743 tcsinv.set_scale(1.0);
744 tcsinv.invert();
745
746 //Get the current rotation
747 Transform temp = Transform();
748 temp.set_trans(n1, n2, n3);
749 Transform cc = tcsinv*temp;
750 Vec3f cctrans = cc.get_trans();
751
752 //set the right rotation
753 Dict spinrot = Dict();
754 spinrot["type"] = "spin";
755 spinrot["omega"] = omega;
756 spinrot["n1"] = cctrans[0];
757 spinrot["n2"] = cctrans[1];
758 spinrot["n3"] = cctrans[2];
759 Transform rightrot = Transform(spinrot);
760 rotate_origin(rightrot);
761}
762
763void Transform::rotate(const Transform& by)
764{
765 vector<float> multmatrix = by.get_matrix();
766 // First Multiply and put the result in a temp matrix
767 Transform result;
768 for (int i=0; i<3; i++) {
769 for (int j=0; j<4; j++) {
770 result[i][j] = multmatrix[i*4]*matrix[0][j] + multmatrix[i*4+1]*matrix[1][j] + multmatrix[i*4+2]*matrix[2][j];
771 }
772 }
773 //Then put the result from the tmep matrix in the original one
774 for (int i=0; i<3; i++) {
775 for (int j=0; j<4; j++) {
776 matrix[i][j] = result[i][j];
777 }
778 }
779}
780
781Transform Transform::negate() const
782{
783 Transform t(*this);
784 for(unsigned int i = 0; i < 3; ++i) {
785 for(unsigned int j = 0; j < 4; ++j) {
786 t.set(i,j,t[i][j]*-1);
787 }
788 }
789 return t;
790}
791
792Transform Transform::get_hflip_transform() const {
793
794 Dict rot = get_rotation("eman");
795 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
796 rot["phi"] = 180.0f - static_cast<float>(rot["phi"]);
797 // This is the same as new_alt= 180-alt, new_phi=-phi, new_az=180+az
798
799 Transform ret(*this); // Is the identity
800 ret.set_rotation(rot);
801
802 Vec3f trans = get_trans();
803 trans[0] = -trans[0];
804 ret.set_trans(trans);
805
806// ret.set_mirror(self.get_mirror());
807
808 return ret;
809}
810
811Transform Transform::get_vflip_transform() const {
812
813 Dict rot = get_rotation("eman");
814 rot["alt"] = 180.0f + static_cast<float>(rot["alt"]);
815 rot["phi"] = - static_cast<float>(rot["phi"]);
816
817 // This is the same as new_alt= 180-alt, new_phi=180-phi, new_az=180+az
818
819 Transform ret(*this);
820 ret.set_rotation(rot);
821
822 Vec3f trans = get_trans();
823 trans[1] = -trans[1];
824 ret.set_trans(trans);
825
826 return ret;
827}
828
829Dict Transform::get_rotation(const string& euler_type) const
830{
831 Dict result;
832 //printf(" Hello from Transform \n");
833
834 //float max = 1 - ERR_LIMIT;
835 float scale;
836 bool x_mirror;
837 get_scale_and_mirror(scale,x_mirror);
838 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
839
840 double cosalt = matrix[2][2]/scale;
841 double sinalt = sqrt(matrix[2][0]*matrix[2][0]+matrix[2][1]*matrix[2][1])/scale; // PRB May 2021
842 double x_mirror_scale = (x_mirror ? -1.0f : 1.0f);
843 double inv_scale = 1.0f/scale;
844
845 double az = 0;
846 double alt = 0;
847 double phi = 0;
848 double phiS = 0; // like az (but in SPIDER ZYZ)
849 double psiS = 0; // like phi (but in SPIDER ZYZ)
850
851 // get alt, az, phi in EMAN convention
852
853 if (cosalt >= 1) { // that is, alt close to 0
854 alt = 0;
855 az = 0;
856 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
857 } else if (cosalt <= -1) { // that is, alt close to 180
858 alt = 180;
859 az = 0;
860 phi = (double)EMConsts::rad2deg * atan2(-x_mirror_scale*matrix[0][1], x_mirror_scale*matrix[0][0]);
861 } else { // for non exceptional cases: 0 < alt < 180
862
863 az = (double)EMConsts::rad2deg * atan2(scale*matrix[2][0], -scale*matrix[2][1]);
864
865 if (matrix[2][2]==0.0)
866 alt = 90.0;
867 else
868 alt = (double)EMConsts::rad2deg * atan(sqrt((double)matrix[2][0]*matrix[2][0]+(double)matrix[2][1]*matrix[2][1])/fabs(matrix[2][2]));
869
870 if (matrix[2][2] * scale < 0)
871 alt = 180.0f-alt;
872
873 phi = (double)EMConsts::rad2deg * atan2(x_mirror_scale*(double)matrix[0][2], (double)matrix[1][2]);
874
875 } // ends separate cases: alt close to 0, 180, or neither
876
877 //phi = phi-360.0*floor(phi/360.0);
878 //az = az -360.0*floor(az/360.0);
879
880// get phiS, psiS (SPIDER)
881 if (cosalt >= 1) { // that is, alt close to 0
882 phiS = 0;
883 psiS = phi;
884 } else if (cosalt <= -1) { // that is, alt close to 180
885 phiS = 0;
886 psiS = phi + 180.0;
887 } else {
888 phiS = az - 90.0;
889 psiS = phi + 90.0;
890 }
891
892 phiS = phiS-360.0*floor(phiS/360.0);
893 psiS = psiS-360.0*floor(psiS/360.0);
894
895// do some quaternionic stuff here
896 double xtilt = 0;
897 double ytilt = 0;
898 double ztilt = 0;
899
900
901 string type = Util::str_to_lower(euler_type);
902
903 result["type"] = type;
904 if (type == "2d") {
905 assert_valid_2d();
906 result["alpha"] = phi;
907 } else if (type == "eman") {
908// assert_consistent_type(THREED);
909 result["az"] = az;
910 result["alt"] = alt;
911 result["phi"] = phi;
912 } else if (type == "imagic") {
913// assert_consistent_type(THREED);
914 result["alpha"] = az;
915 result["beta"] = alt;
916 result["gamma"] = phi;
917 } else if (type == "spider") {
918// assert_consistent_type(THREED);
919 result["phi"] = phiS; // The first Euler like az
920 result["theta"] = alt;
921 result["psi"] = psiS;
922 } else if (type == "mrc") {
923// assert_consistent_type(THREED);
924 result["phi"] = phiS;
925 result["theta"] = alt;
926 result["omega"] = psiS;
927 } else if (type == "xyz") { // need to double-check these 3 equations ********
928// assert_consistent_type(THREED);
929 xtilt = atan2(-sin(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt),cos(EMConsts::deg2rad*alt));
930 ytilt = asin( cos(EMConsts::deg2rad*phiS)*sin(EMConsts::deg2rad*alt));
931 ztilt = psiS*EMConsts::deg2rad - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
932
933 xtilt *= EMConsts::rad2deg; ytilt *= EMConsts::rad2deg; ztilt *= EMConsts::rad2deg;
934 xtilt = xtilt-360*.0*floor((xtilt+180.0)/360.0);
935 ytilt = ytilt-360*.0*floor((ytilt+180.0)/360.0); //already in range [-90,90] but anyway...
936 ztilt = ztilt-360*.0*floor((ztilt+180.0)/360.0);
937
938 result["xtilt"] = xtilt;
939 result["ytilt"] = ytilt;
940 result["ztilt"] = ztilt;
941 } else if ((type == "quaternion") || (type == "spin") || (type == "sgirot")) {
942
943 // The cosOover2 is also e0
944// double nphi = (az-phi)/2.0;
945// double cosOover2 = cos((az+phi)*EMConsts::deg2rad/2.0) * cos(alt*EMConsts::deg2rad/2.0);
946// printf("%f %f %f",matrix[0][0],matrix[1][1],matrix[2][2]);
947 double traceR = matrix[0][0]+matrix[1][1]+matrix[2][2]; // This should be 1 + 2 cos omega
948 double cosomega = (traceR-1.0)/2.0;
949
950 double A0 =(matrix[1][2]-matrix[2][1])/2.0;
951 double A1 =(matrix[2][0]-matrix[0][2])/2.0;
952 double A2 =(matrix[0][1]-matrix[1][0])/2.0;
953 double sinomega = sqrt(A0*A0+A1*A1+A2*A2);
954
955 if (cosomega> 1.0) cosomega= 1.0;
956 if (cosomega<-1.0) cosomega=-1.0;
957
958 // matrix(x,y)-matrix(y,x) = 2 n_z sin(omega) etc
959 // trace matrix = 1 + 2 cos(omega)
960 double cosOover2= sqrt((1.0 +cosomega)/2.0);
961 //double sinOover2= sinomega/cosOover2/2.0;
962 double sinOover2= sqrt((1.0 -cosomega)/2.0);
963 //double sinomega = 2* sinOover2*cosOover2; // Not a great formula. See above
964
965 double n1 = A0/sinomega;
966 double n2 = A1/sinomega;
967 double n3 = A2/sinomega;
968
969
970
971 if (sinOover2==1) {// This will also mean sinomega=0, omega =pi,
972 n1 = sqrt((matrix[0][0]+1)/2.0) ;
973 n2 = sqrt((matrix[1][1]+1)/2.0) ;
974 n3 = sqrt((matrix[2][2]+1)/2.0) ;
975 if (n1 <= n2){ // n1 is the smallest
976 if (n2<=n3) {
977 n3 = sqrt((matrix[2][2]+1)/2.0) ; n1=matrix[0][2]/n3/2.0; n2=matrix[1][2]/n3/2.0; }
978 else { n2 = sqrt((matrix[1][1]+1)/2.0) ; n1=matrix[0][1]/n2/2.0; n3=matrix[2][1]/n2/2.0; }}
979 else { // n2 is the smallest
980 if (n1<=n3) {
981 n3 = sqrt((matrix[2][2]+1)/2.0) ; n1=matrix[0][2]/n3/2.0; n2=matrix[1][2]/n3/2.0;}
982 else { n1 = sqrt((matrix[0][0]+1)/2.0) ; n2=matrix[1][0]/n1/2.0; n3=matrix[2][0]/n1/2.0;}}
983 }
984
985
986 if (sinOover2==0) {// This will also mean omega =0,
987 n1 = 0;
988 n2 = 0;
989 n3 = 1;
990 }
991
992
993// printf("traceR=%lf,OneMinusCosomega=%lf,sinOover2=%lf,cosOover2=%lf,sinomega=%lf,cosomega=%lf,n3=%lf \n",traceR,1-cosomega,sinOover2,cosOover2,sinomega,cosomega,n3);
994
995 if (type == "quaternion"){
996 result["e0"] = cosOover2 ;
997 result["e1"] = sinOover2 * n1 ;
998 result["e2"] = sinOover2 * n2;
999 result["e3"] = sinOover2 * n3;
1000 }
1001
1002 if (type == "spin"){
1003 double Omega = EMConsts::rad2deg * asin(sinomega);
1004 if (cosomega< 0) { Omega = 180-Omega;}
1005 result["omega"] = Omega; //Changed by PRB
1006 result["n1"] = n1;
1007 result["n2"] = n2;
1008 result["n3"] = n3;
1009 }
1010
1011 if (type == "sgirot"){
1012 result["q"] = EMConsts::rad2deg * acos(cosomega);
1013 result["n1"] = n1;
1014 result["n2"] = n2;
1015 result["n3"] = n3;
1016 }
1017
1018 } else if (type == "matrix") {
1019// assert_consistent_type(THREED);
1020 result["m11"] = x_mirror_scale*matrix[0][0]*inv_scale;
1021 result["m12"] = x_mirror_scale*matrix[0][1]*inv_scale;
1022 result["m13"] = x_mirror_scale*matrix[0][2]*inv_scale;
1023 result["m21"] = matrix[1][0]*inv_scale;
1024 result["m22"] = matrix[1][1]*inv_scale;
1025 result["m23"] = matrix[1][2]*inv_scale;
1026 result["m31"] = matrix[2][0]*inv_scale;
1027 result["m32"] = matrix[2][1]*inv_scale;
1028 result["m33"] = matrix[2][2]*inv_scale;
1029 } else {
1030 throw InvalidStringException(euler_type, "unknown Euler Type");
1031 }
1032
1033 return result;
1034}
1035
1036void Transform::set_trans(const float& x, const float& y, const float& z)
1037{
1038 bool x_mirror = get_mirror();
1039
1040 if (x_mirror) matrix[0][3] = -x;
1041 else matrix[0][3] = x;
1042 matrix[1][3] = y;
1043 matrix[2][3] = z;
1044}
1045
1046Vec3f Transform::get_trans() const
1047{
1048 // No type asserted
1049 bool x_mirror = get_mirror();
1050 Vec3f v;
1051 if (x_mirror) v[0] = -matrix[0][3];
1052 else v[0] = matrix[0][3];
1053 v[1] = matrix[1][3];
1054 v[2] = matrix[2][3];
1055
1056 Util::apply_precision(v[0],ERR_LIMIT);
1057 Util::apply_precision(v[1],ERR_LIMIT);
1058 Util::apply_precision(v[2],ERR_LIMIT);
1059
1060 return v;
1061}
1062
1063void Transform::translate(const float& tx, const float& ty, const float& tz)
1064{
1065 bool x_mirror = get_mirror();
1066 if (x_mirror) matrix[0][3] = -matrix[0][3] + tx;
1067 else matrix[0][3] = matrix[0][3] + tx;
1068 matrix[1][3] = matrix[1][3] + ty;
1069 matrix[2][3] = matrix[2][3] + tz;
1070}
1071
1072void Transform::translate_newBasis(const Transform& tcs, const float& tx, const float& ty, const float& tz)
1073{
1074 //Get the rotational inverse
1075 Transform tcsinv = Transform(tcs);
1076 tcsinv.set_trans(0.0, 0.0, 0.0);
1077 tcsinv.invert();
1078
1079 //Now move the coordinate system
1080 Transform temp = Transform();
1081 temp.set_trans(tx, ty, tz);
1082 Transform nb_trans = tcsinv*temp;
1083
1084 translate(nb_trans.get_trans());
1085
1086}
1087
1088Vec2f Transform::get_trans_2d() const
1089{
1090 bool x_mirror = get_mirror();
1091 Vec2f v;
1092 if (x_mirror) v[0] = -matrix[0][3];
1093 else v[0] = matrix[0][3];
1094 v[1] = matrix[1][3];
1095 return v;
1096}
1097
1098
1099
1100Vec3f Transform::get_pre_trans() const
1101{
1102 Transform T(*this);
1103 T.set_trans(0,0,0);
1104 T.invert();
1105
1106 Transform soln = T*(*this);
1107// soln.printme();
1108 return soln.get_trans();
1109}
1110
1111Vec2f Transform::get_pre_trans_2d() const
1112{
1113 Transform T(*this);
1114 T.set_trans(0,0,0);
1115 T.invert();
1116
1117 Transform soln = T*(*this);
1118// soln.printme();
1119 return soln.get_trans_2d();
1120}
1121
1122
1123void Transform::set_scale(const float& new_scale) {
1124 if (new_scale <= 0) {
1125 throw InvalidValueException(new_scale,"The scale factor in a Transform object must be positive and non zero");
1126 }
1127 // Transform = MTSR (Mirroring, Translation, Scaling, Rotate)
1128 // So changing the scale boils down to this....
1129
1130 float old_scale = get_scale();
1131
1132 float n_scale = new_scale;
1133 Util::apply_precision(n_scale,ERR_LIMIT);
1134
1135 float corrected_scale = n_scale/old_scale;
1136 if ( corrected_scale != 1.0 ) {
1137 for(int i = 0; i < 3; ++i ) {
1138 for(int j = 0; j < 3; ++j ) {
1139 matrix[i][j] *= corrected_scale;
1140 }
1141 }
1142 }
1143}
1144
1145float Transform::get_scale() const {
1146 float determinant = get_determinant();
1147 if (determinant < 0 ) determinant *= -1;
1148
1149 float scale = std::pow(determinant,1.0f/3.0f);
1150 int int_scale = static_cast<int>(scale);
1151 float scale_residual = scale-static_cast<float>(int_scale);
1152 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
1153
1154 Util::apply_precision(scale, ERR_LIMIT);
1155
1156 return scale;
1157}
1158
1159void Transform::scale(const float& scale)
1160{
1161 float determinant = get_determinant();
1162 if (determinant < 0) determinant *= -1.0f;
1163 float newscale = std::pow(determinant,1.0f/3.0f) + scale;
1164 if(newscale > 0.0001) set_scale(newscale); // If scale ~ 0 things blowup, so we need a little fudge factor
1165}
1166
1167void print_matrix(gsl_matrix* M, unsigned int r, unsigned int c, const string& message ) {
1168 cout << "Message is " << message << endl;
1169 for ( unsigned int i = 0; i < r; ++i )
1170 {
1171 for ( unsigned int j = 0; j < c; ++j )
1172 {
1173 cout << gsl_matrix_get(M,i,j) << " ";
1174 }
1175 cout << endl;
1176 }
1177}
1178
1179void Transform::orthogonalize()
1180{
1181 float scale;
1182 bool x_mirror;
1183 get_scale_and_mirror(scale,x_mirror);
1184 if (scale == 0) throw UnexpectedBehaviorException("The determinant of the Transform is 0. This is unexpected.");
1185 double inv_scale = 1.0/static_cast<double>(scale);
1186 double mirror_scale = (x_mirror == true ? -1.0:1.0);
1187
1188 gsl_matrix * R = gsl_matrix_calloc(3,3);
1189 for ( unsigned int i = 0; i < 3; ++i )
1190 {
1191 for ( unsigned int j = 0; j < 3; ++j )
1192 {
1193 if (i == 0 && mirror_scale != 1.0 ) {
1194 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*mirror_scale*inv_scale );
1195 }
1196 else {
1197 gsl_matrix_set( R, i, j, static_cast<double>(matrix[i][j])*inv_scale );
1198 }
1199 }
1200 }
1201
1202 gsl_matrix * V = gsl_matrix_calloc(3,3);
1203 gsl_vector * S = gsl_vector_calloc(3);
1204 gsl_vector * work = gsl_vector_calloc(3);
1205 gsl_linalg_SV_decomp (R, V, S, work); // Now R is U of the SVD R = USV^T
1206
1207 gsl_matrix * Soln = gsl_matrix_calloc(3,3);
1208 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, R, V, 0.0, Soln);
1209
1210 for ( unsigned int i = 0; i < 3; ++i )
1211 {
1212 for ( unsigned int j = 0; j < 3; ++j )
1213 {
1214 matrix[i][j] = static_cast<float>( gsl_matrix_get(Soln,i,j) );
1215 }
1216 }
1217
1218 // Apply scale if it existed previously
1219 if (scale != 1.0f) {
1220 for(int i=0; i<3; ++i) {
1221 for(int j=0; j<3; ++j) {
1222 matrix[i][j] *= scale;
1223 }
1224 }
1225 }
1226
1227 // Apply post x mirroring if it was applied previouslys
1228 if ( x_mirror ) {
1229 for(int j=0; j<3; ++j) {
1230 matrix[0][j] *= -1.0f;
1231 }
1232 }
1233
1234 gsl_matrix_free(V); gsl_matrix_free(R); gsl_matrix_free(Soln);
1235 gsl_vector_free(S); gsl_vector_free(work);
1236}
1237
1238void Transform::set_mirror(const bool x_mirror ) {
1239
1240 bool old_x_mirror = get_mirror();
1241 if (old_x_mirror == x_mirror) return; // The user is setting the same value
1242 else {
1243 // Toggle the mirroring operation
1244 for (int j = 0; j < 4; ++j ) {
1245 matrix[0][j] *= -1;
1246 }
1247 }
1248}
1249
1250bool Transform::get_mirror() const {
1251 float determinant = get_determinant();
1252
1253 bool x_mirror = false;
1254 if ( determinant < 0 ) x_mirror = true;
1255
1256 return x_mirror;
1257
1258}
1259
1260void Transform::get_scale_and_mirror(float& scale, bool& x_mirror) const {
1261
1262 float determinant = get_determinant();
1263 x_mirror = false;
1264 if ( determinant < 0 ) {
1265 x_mirror = true;
1266 determinant *= -1;
1267 }
1268 if (determinant != 1 ) {
1269 scale = std::pow(determinant,1.0f/3.0f);
1270 int int_scale = static_cast<int>(scale);
1271 float scale_residual = scale-static_cast<float>(int_scale);
1272 if ( scale_residual < ERR_LIMIT ) { scale = static_cast<float>(int_scale); };
1273 }
1274 else scale = 1;
1275
1276 Util::apply_precision(scale,ERR_LIMIT);
1277}
1278
1279float Transform::get_determinant() const
1280{
1281 float det;
1282 double det2;
1283 det2 = matrix[0][0]*((double)matrix[1][1]*matrix[2][2]-(double)matrix[2][1]*matrix[1][2]);
1284 det2 -= matrix[0][1]*((double)matrix[1][0]*matrix[2][2]-(double)matrix[2][0]*matrix[1][2]);
1285 det2 += matrix[0][2]*((double)matrix[1][0]*matrix[2][1]-(double)matrix[2][0]*matrix[1][1]);
1286
1287 det = (float)det2;
1288 Util::apply_precision(det,ERR_LIMIT);
1289
1290 return det;
1291}
1292
1293void Transform::invert() {
1294
1295 double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
1296 double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
1297 double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
1298 double v0 = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
1299
1300 double cof00 = m11*m22-m12*m21;
1301 double cof11 = m22*m00-m20*m02;
1302 double cof22 = m00*m11-m01*m10;
1303 double cof01 = m10*m22-m20*m12;
1304 double cof02 = m10*m21-m20*m11;
1305 double cof12 = m00*m21-m01*m20;
1306 double cof10 = m01*m22-m02*m21;
1307 double cof20 = m01*m12-m02*m11;
1308 double cof21 = m00*m12-m10*m02;
1309
1310 double det = m00* cof00 + m02* cof02 -m01*cof01;
1311
1312 matrix[0][0] = (float)(cof00/det);
1313 matrix[0][1] = - (float)(cof10/det);
1314 matrix[0][2] = (float)(cof20/det);
1315 matrix[1][0] = - (float)(cof01/det);
1316 matrix[1][1] = (float)(cof11/det);
1317 matrix[1][2] = - (float)(cof21/det);
1318 matrix[2][0] = (float)(cof02/det);
1319 matrix[2][1] = - (float)(cof12/det);
1320 matrix[2][2] = (float)(cof22/det);
1321
1322 matrix[0][3] = (float)((- cof00*v0 + cof10*v1 - cof20*v2)/det);
1323 matrix[1][3] = (float)(( cof01*v0 - cof11*v1 + cof21*v2)/det);
1324 matrix[2][3] = (float)((- cof02*v0 + cof12*v1 - cof22*v2)/det);
1325}
1326
1327Transform Transform::inverse() const {
1328 Transform t(*this);
1329 t.invert();
1330 return t;
1331}
1332
1333void Transform::transpose_inplace() {
1334 float tempij;
1335 for (int i = 0; i < 3; i++) {
1336 for (int j = 0; j < i; j++) {
1337 if (i != j) {
1338 tempij= matrix[i][j];
1339 matrix[i][j] = matrix[j][i];
1340 matrix[j][i] = tempij;
1341 }
1342 }
1343 }
1344}
1345
1346Transform Transform::transpose() const {
1347 Transform t(*this);
1348 t.transpose_inplace();
1349 return t;
1350}
1351
1352
1353Transform EMAN::operator*(const Transform & M2, const Transform & M1) // YYY
1354{
1355 Transform result;
1356 for (int i=0; i<3; i++) {
1357 for (int j=0; j<4; j++) {
1358 result[i][j] = M2[i][0] * M1[0][j] + M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
1359 }
1360 result[i][3] += M2[i][3];
1361 }
1362
1363 return result;
1364}
1365
1366void Transform::assert_valid_2d() const {
1367 int rotation_error = 0;
1368 int translation_error = 0;
1369 if (fabs(matrix[2][0]) > ERR_LIMIT) rotation_error++;
1370 if (fabs(matrix[2][1]) > ERR_LIMIT) rotation_error++;
1371 if (fabs(matrix[2][3]) > ERR_LIMIT) translation_error++;
1372 if (fabs(matrix[0][2]) > ERR_LIMIT) rotation_error++;
1373 if (fabs(matrix[1][2]) > ERR_LIMIT) rotation_error++;
1374// if (fabs(matrix[2][2]-1.0) >ERR_LIMIT) rotation_error++;
1375 if (matrix[2][2] <=0) rotation_error++; // previous line commented out due to addition of scaling. This line will insure we don't have alt=180.0
1376 if ( translation_error && rotation_error ) {
1377 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and 3D translations. This object can not be considered 2D");
1378 } else if ( translation_error ) {
1379 throw UnexpectedBehaviorException("Error, the internal matrix contains a non zero z component for a 3D translation. This object can not be considered 2D");
1380 }
1381 else if ( rotation_error ) {
1382 throw UnexpectedBehaviorException("Error, the internal matrix contains 3D rotations and this object can not be considered 2D");
1383 }
1384
1385}
1386
1387
1388
1389Transform Transform::get_sym(const string & sym_name, int n) const
1390{
1391 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
1392 Transform ret;
1393 ret = (*this) * sym->get_sym(n);
1394 delete sym;
1395 return ret;
1396}
1397
1398Transform Transform::get_sym_sparx(const string & sym_name, int n) const
1399{
1400 Transform ret;
1401 vector<Transform> tut;
1402 tut = ret.get_sym_proj(sym_name);
1403 ret = (*this) * tut[n];
1404 return ret;
1405}
1406
1407vector<Transform > Transform::get_sym_proj(const string & sym_name) const
1408{
1409 vector<Transform> ret;
1410 int nsym;
1411 string lstr = Util::str_to_lower(sym_name);
1412 if( lstr == "tet" ) {
1413 nsym = 12;
1414 static float TET[108] = {
1415 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1416 -0.5f, 0.86602540378f, 0.0f, -0.86602540378f, -0.5f, 0.0f, 0.0f, 0.0f, 1.0f,
1417 -0.5f, -0.86602540378f, 0.0f, 0.86602540378f, -0.5f, -0.0f, 0.0f, 0.0f, 1.0f,
1418 -0.16666666667f, 0.86602540378f, -0.47140452079f, 0.28867513459f, 0.5f, 0.81649658093f, 0.94280904158f, 0.0f, -0.33333333333f,
1419 0.33333333333f, 0.0f, 0.94280904158f, 0.0f, -1.0f, 0.0f, 0.94280904158f, 0.0f, -0.33333333333f,
1420 -0.16666666667f, -0.86602540378f, -0.47140452079f, -0.28867513459f, 0.5f, -0.81649658093f, 0.94280904158f, 0.0f, -0.33333333333f,
1421 -0.66666666667f, -0.57735026919f, -0.47140452079f, -0.57735026919f, 0.0f, 0.81649658093f, -0.47140452079f, 0.81649658093f, -0.33333333333f,
1422 -0.16666666667f, 0.28867513459f, 0.94280904158f, 0.86602540378f, 0.5f, 0.0f, -0.47140452079f, 0.81649658093f, -0.33333333333f,
1423 0.83333333333f, 0.28867513459f, -0.47140452079f, -0.28867513459f, -0.5f, -0.81649658093f, -0.47140452079f, 0.81649658093f, -0.33333333333f,
1424 0.83333333333f, -0.28867513459f, -0.47140452079f, 0.28867513459f, -0.5f, 0.81649658093f, -0.47140452079f, -0.81649658093f, -0.33333333333f,
1425 -0.16666666667f, -0.28867513459f, 0.94280904158f, -0.86602540378f, 0.5f, 0.0f, -0.47140452079f, -0.81649658093f, -0.33333333333f,
1426 -0.66666666667f, 0.57735026919f, -0.47140452079f, 0.57735026919f, -0.0f, -0.81649658093f, -0.47140452079f, -0.81649658093f, -0.33333333333f
1427 };
1428
1429 for (int k=0;k<nsym;k++) {
1430 Transform t;
1431 for (int i=0; i<3; i++) {
1432 for (int j=0; j<3; j++) {
1433 t.matrix[i][j] = TET[9*k + i*3 +j];
1434 }
1435 }
1436 //vector<float> z = t.get_matrix();
1437 //for (int i=0; i<12; i++) cout<<z[i]<<endl;
1438 ret.push_back( (*this) * t );
1439 }
1440 } else if( lstr == "oct" ) {
1441 nsym = 24;
1442 static float TET[216] = {
1443 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1444 0.0f, 1.0f, 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1445 -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1446 0.0f, -1.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1447 0.0f, 0.0f, -1.0f, 0.0f, 1.0f, 0.0f, 1.0f, 0.0f, 0.0f,
1448 0.0f, 0.0f, -1.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f,
1449 0.0f, 0.0f, -1.0f, 0.0f, -1.0f, 0.0f, -1.0f, 0.0f, 0.0f,
1450 0.0f, 0.0f, -1.0f, 1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f,
1451 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 1.0f, 0.0f, 0.0f,
1452 -1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 1.0f, 0.0f,
1453 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, -1.0f, 0.0f, 0.0f,
1454 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, -1.0f, 0.0f,
1455 0.0f, 0.0f, 1.0f, 0.0f, -1.0f, 0.0f, 1.0f, 0.0f, 0.0f,
1456 0.0f, 0.0f, 1.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f,
1457 0.0f, 0.0f, 1.0f, 0.0f, 1.0f, 0.0f, -1.0f, 0.0f, 0.0f,
1458 0.0f, 0.0f, 1.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f,
1459 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 1.0f, 0.0f, 0.0f,
1460 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f, 1.0f, 0.0f,
1461 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, -1.0f, -1.0f, 0.0f, 0.0f,
1462 -1.0f, 0.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f, -1.0f, 0.0f,
1463 -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, -1.0f,
1464 0.0f, 1.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 0.0f, -1.0f,
1465 1.0f, 0.0f, 0.0f, 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, -1.0f,
1466 0.0f, -1.0f, 0.0f, -1.0f, 0.0f, 0.0f, 0.0f, 0.0f, -1.0f
1467 };
1468
1469 for (int k=0;k<nsym;k++) {
1470 Transform t;
1471 for (int i=0; i<3; i++) {
1472 for (int j=0; j<3; j++) {
1473 t.matrix[i][j] = TET[9*k + i*3 +j];
1474 }
1475 }
1476 //vector<float> z = t.get_matrix();
1477 //for (int i=0; i<12; i++) cout<<z[i]<<endl;
1478 ret.push_back( (*this) * t );
1479 }
1480 } else if( lstr == "icos" ) {
1481 nsym = 60;
1482 static float TET[540] = {
1483 1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, 1.0f,
1484 0.30901699437f, 0.9510565163f, 0.0f, -0.9510565163f, 0.30901699437f, 0.0f, 0.0f, 0.0f, 1.0f,
1485 -0.80901699437f, 0.58778525229f, 0.0f, -0.58778525229f, -0.80901699437f, 0.0f, 0.0f, 0.0f, 1.0f,
1486 -0.80901699437f, -0.58778525229f, 0.0f, 0.58778525229f, -0.80901699437f, 0.0f, 0.0f, 0.0f, 1.0f,
1487 0.30901699437f, -0.9510565163f, 0.0f, 0.9510565163f, 0.30901699437f, 0.0f, 0.0f, 0.0f, 1.0f,
1488 0.36180339887f, 0.58778525229f, -0.72360679775f, -0.26286555606f, 0.80901699437f, 0.52573111212f, 0.894427191f, 0.0f, 0.4472135955f,
1489 -0.13819660113f, 0.9510565163f, 0.27639320225f, -0.42532540418f, -0.30901699437f, 0.85065080835f, 0.894427191f, 0.0f, 0.4472135955f,
1490 -0.4472135955f, 0.0f, 0.894427191f, 0.0f, -1.0f, 0.0f, 0.894427191f, 0.0f, 0.4472135955f,
1491 -0.13819660113f, -0.9510565163f, 0.27639320225f, 0.42532540418f, -0.30901699437f, -0.85065080835f, 0.894427191f, 0.0f, 0.4472135955f,
1492 0.36180339887f, -0.58778525229f, -0.72360679775f, 0.26286555606f, 0.80901699437f, -0.52573111212f, 0.894427191f, 0.0f, 0.4472135955f,
1493 -0.4472135955f, 0.52573111212f, -0.72360679775f, -0.85065080835f, 0.0f, 0.52573111212f, 0.27639320225f, 0.85065080835f, 0.4472135955f,
1494 -0.9472135955f, 0.16245984812f, 0.27639320225f, 0.16245984812f, -0.5f, 0.85065080835f, 0.27639320225f, 0.85065080835f, 0.4472135955f,
1495 -0.13819660113f, -0.42532540418f, 0.894427191f, 0.9510565163f, -0.30901699437f, 0.0f, 0.27639320225f, 0.85065080835f, 0.4472135955f,
1496 0.86180339887f, -0.42532540418f, 0.27639320225f, 0.42532540418f, 0.30901699437f, -0.85065080835f, 0.27639320225f, 0.85065080835f, 0.4472135955f,
1497 0.67082039325f, 0.16245984812f, -0.72360679775f, -0.68819096024f, 0.5f, -0.52573111212f, 0.27639320225f, 0.85065080835f, 0.4472135955f,
1498 -0.63819660113f, -0.26286555606f, -0.72360679775f, -0.26286555606f, -0.80901699437f, 0.52573111212f, -0.72360679775f, 0.52573111212f, 0.4472135955f,
1499 -0.4472135955f, -0.85065080835f, 0.27639320225f, 0.52573111212f, 0.0f, 0.85065080835f, -0.72360679775f, 0.52573111212f, 0.4472135955f,
1500 0.36180339887f, -0.26286555606f, 0.894427191f, 0.58778525229f, 0.80901699437f, 0.0f, -0.72360679775f, 0.52573111212f, 0.4472135955f,
1501 0.67082039325f, 0.68819096024f, 0.27639320225f, -0.16245984812f, 0.5f, -0.85065080835f, -0.72360679775f, 0.52573111212f, 0.4472135955f,
1502 0.0527864045f, 0.68819096024f, -0.72360679775f, -0.68819096024f, -0.5f, -0.52573111212f, -0.72360679775f, 0.52573111212f, 0.4472135955f,
1503 0.0527864045f, -0.68819096024f, -0.72360679775f, 0.68819096024f, -0.5f, 0.52573111212f, -0.72360679775f, -0.52573111212f, 0.4472135955f,
1504 0.67082039325f, -0.68819096024f, 0.27639320225f, 0.16245984812f, 0.5f, 0.85065080835f, -0.72360679775f, -0.52573111212f, 0.4472135955f,
1505 0.36180339887f, 0.26286555606f, 0.894427191f, -0.58778525229f, 0.80901699437f, 0.0f, -0.72360679775f, -0.52573111212f, 0.4472135955f,
1506 -0.4472135955f, 0.85065080835f, 0.27639320225f, -0.52573111212f, 0.0f, -0.85065080835f, -0.72360679775f, -0.52573111212f, 0.4472135955f,
1507 -0.63819660113f, 0.26286555606f, -0.72360679775f, 0.26286555606f, -0.80901699437f, -0.52573111212f, -0.72360679775f, -0.52573111212f, 0.4472135955f,
1508 0.67082039325f, -0.16245984812f, -0.72360679775f, 0.68819096024f, 0.5f, 0.52573111212f, 0.27639320225f, -0.85065080835f, 0.4472135955f,
1509 0.86180339887f, 0.42532540418f, 0.27639320225f, -0.42532540418f, 0.30901699437f, 0.85065080835f, 0.27639320225f, -0.85065080835f, 0.4472135955f,
1510 -0.13819660113f, 0.42532540418f, 0.894427191f, -0.9510565163f, -0.30901699437f, 0.0f, 0.27639320225f, -0.85065080835f, 0.4472135955f,
1511 -0.9472135955f, -0.16245984812f, 0.27639320225f, -0.16245984812f, -0.5f, -0.85065080835f, 0.27639320225f, -0.85065080835f, 0.4472135955f,
1512 -0.4472135955f, -0.52573111212f, -0.72360679775f, 0.85065080835f, 0.0f, -0.52573111212f, 0.27639320225f, -0.85065080835f, 0.4472135955f,
1513 -0.36180339887f, -0.26286555606f, -0.894427191f, -0.58778525229f, 0.80901699437f, 0.0f, 0.72360679775f, 0.52573111212f, -0.4472135955f,
1514 -0.67082039325f, 0.68819096024f, -0.27639320225f, 0.16245984812f, 0.5f, 0.85065080835f, 0.72360679775f, 0.52573111212f, -0.4472135955f,
1515 -0.0527864045f, 0.68819096024f, 0.72360679775f, 0.68819096024f, -0.5f, 0.52573111212f, 0.72360679775f, 0.52573111212f, -0.4472135955f,
1516 0.63819660113f, -0.26286555606f, 0.72360679775f, 0.26286555606f, -0.80901699437f, -0.52573111212f, 0.72360679775f, 0.52573111212f, -0.4472135955f,
1517 0.4472135955f, -0.85065080835f, -0.27639320225f, -0.52573111212f, 0.0f, -0.85065080835f, 0.72360679775f, 0.52573111212f, -0.4472135955f,
1518 0.13819660113f, -0.42532540418f, -0.894427191f, -0.9510565163f, -0.30901699437f, 0.0f, -0.27639320225f, 0.85065080835f, -0.4472135955f,
1519 -0.86180339887f, -0.42532540418f, -0.27639320225f, -0.42532540418f, 0.30901699437f, 0.85065080835f, -0.27639320225f, 0.85065080835f, -0.4472135955f,
1520 -0.67082039325f, 0.16245984812f, 0.72360679775f, 0.68819096024f, 0.5f, 0.52573111212f, -0.27639320225f, 0.85065080835f, -0.4472135955f,
1521 0.4472135955f, 0.52573111212f, 0.72360679775f, 0.85065080835f, 0.0f, -0.52573111212f, -0.27639320225f, 0.85065080835f, -0.4472135955f,
1522 0.9472135955f, 0.16245984812f, -0.27639320225f, -0.16245984812f, -0.5f, -0.85065080835f, -0.27639320225f, 0.85065080835f, -0.4472135955f,
1523 0.4472135955f, 0.0f, -0.894427191f, 0.0f, -1.0f, 0.0f, -0.894427191f, 0.0f, -0.4472135955f,
1524 0.13819660113f, -0.9510565163f, -0.27639320225f, -0.42532540418f, -0.30901699437f, 0.85065080835f, -0.894427191f, 0.0f, -0.4472135955f,
1525 -0.36180339887f, -0.58778525229f, 0.72360679775f, -0.26286555606f, 0.80901699437f, 0.52573111212f, -0.894427191f, 0.0f, -0.4472135955f,
1526 -0.36180339887f, 0.58778525229f, 0.72360679775f, 0.26286555606f, 0.80901699437f, -0.52573111212f, -0.894427191f, 0.0f, -0.4472135955f,
1527 0.13819660113f, 0.9510565163f, -0.27639320225f, 0.42532540418f, -0.30901699437f, -0.85065080835f, -0.894427191f, 0.0f, -0.4472135955f,
1528 0.13819660113f, 0.42532540418f, -0.894427191f, 0.9510565163f, -0.30901699437f, 0.0f, -0.27639320225f, -0.85065080835f, -0.4472135955f,
1529 0.9472135955f, -0.16245984812f, -0.27639320225f, 0.16245984812f, -0.5f, 0.85065080835f, -0.27639320225f, -0.85065080835f, -0.4472135955f,
1530 0.4472135955f, -0.52573111212f, 0.72360679775f, -0.85065080835f, 0.0f, 0.52573111212f, -0.27639320225f, -0.85065080835f, -0.4472135955f,
1531 -0.67082039325f, -0.16245984812f, 0.72360679775f, -0.68819096024f, 0.5f, -0.52573111212f, -0.27639320225f, -0.85065080835f, -0.4472135955f,
1532 -0.86180339887f, 0.42532540418f, -0.27639320225f, 0.42532540418f, 0.30901699437f, -0.85065080835f, -0.27639320225f, -0.85065080835f, -0.4472135955f,
1533 -0.36180339887f, 0.26286555606f, -0.894427191f, 0.58778525229f, 0.80901699437f, 0.0f, 0.72360679775f, -0.52573111212f, -0.4472135955f,
1534 0.4472135955f, 0.85065080835f, -0.27639320225f, 0.52573111212f, 0.0f, 0.85065080835f, 0.72360679775f, -0.52573111212f, -0.4472135955f,
1535 0.63819660113f, 0.26286555606f, 0.72360679775f, -0.26286555606f, -0.80901699437f, 0.52573111212f, 0.72360679775f, -0.52573111212f, -0.4472135955f,
1536 -0.0527864045f, -0.68819096024f, 0.72360679775f, -0.68819096024f, -0.5f, -0.52573111212f, 0.72360679775f, -0.52573111212f, -0.4472135955f,
1537 -0.67082039325f, -0.68819096024f, -0.27639320225f, -0.16245984812f, 0.5f, -0.85065080835f, 0.72360679775f, -0.52573111212f, -0.4472135955f,
1538 -1.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f, 0.0f, 0.0f, -1.0f,
1539 -0.30901699437f, 0.9510565163f, 0.0f, 0.9510565163f, 0.30901699437f, 0.0f, 0.0f, 0.0f, -1.0f,
1540 0.80901699437f, 0.58778525229f, 0.0f, 0.58778525229f, -0.80901699437f, 0.0f, 0.0f, 0.0f, -1.0f,
1541 0.80901699437f, -0.58778525229f, 0.0f, -0.58778525229f, -0.80901699437f, 0.0f, 0.0f, 0.0f, -1.0f,
1542 -0.30901699437f, -0.9510565163f, 0.0f, -0.9510565163f, 0.30901699437f, 0.0f, 0.0f, 0.0f, -1.0f
1543 };
1544
1545 for (int k=0;k<nsym;k++) {
1546 Transform t;
1547 for (int i=0; i<3; i++) {
1548 for (int j=0; j<3; j++) {
1549 t.matrix[i][j] = TET[9*k + i*3 +j];
1550 }
1551 }
1552 //vector<float> z = t.get_matrix();
1553 //for (int i=0; i<12; i++) cout<<z[i]<<endl;
1554 ret.push_back( (*this) * t );
1555 }
1556 } else {
1557 Symmetry3D* sym = Factory<Symmetry3D>::get(lstr);
1558 nsym = sym->get_nsym();
1559 for (int k=0;k<nsym;k++) {
1560 Transform t;
1561 t = sym->get_sym(k);
1562 ret.push_back( (*this) * t );
1563 }
1564 delete sym;
1565 }
1566 return ret;
1567}
1568
1569
1570int Transform::get_nsym(const string & sym_name)
1571{
1572 Symmetry3D* sym = Factory<Symmetry3D>::get(sym_name);
1573 int nsym = sym->get_nsym();
1574 delete sym;
1575 return nsym;
1576}
1577
1578//
1579//Transform3D::Transform3D() // C1
1580//{
1581// init();
1582//}
1583//
1584//Transform3D::Transform3D( const Transform3D& rhs )
1585//{
1586// for( int i=0; i < 4; ++i )
1587// {
1588// for( int j=0; j < 4; ++j )
1589// {
1590// matrix[i][j] = rhs.matrix[i][j];
1591// }
1592// }
1593//}
1594//
1596//Transform3D::Transform3D(const float& az, const float& alt, const float& phi)
1597//{
1598// init();
1599// set_rotation(az,alt,phi);
1600//}
1601//
1602//
1604//Transform3D::Transform3D(const float& az, const float& alt, const float& phi, const Vec3f& posttrans )
1605//{
1606// init(); // This is called in set_rotation
1607// set_rotation(az,alt,phi);
1608// set_posttrans(posttrans);
1609//}
1610//
1611//Transform3D::Transform3D(const float& m11, const float& m12, const float& m13,
1612// const float& m21, const float& m22, const float& m23,
1613// const float& m31, const float& m32, const float& m33)
1614//{
1615// init();
1616// set_rotation(m11,m12,m13,m21,m22,m23,m31,m32,m33);
1617//}
1618//
1620//Transform3D::Transform3D(EulerType euler_type, const float& a1, const float& a2, const float& a3)
1621//{
1622// init();
1623// set_rotation(euler_type,a1,a2,a3);
1624//}
1625//
1626//Transform3D::Transform3D(EulerType euler_type, const float& a1, const float& a2, const float& a3, const float& a4)
1627//{
1628// init();
1629// set_rotation(euler_type,a1,a2,a3,a4);
1630//}
1631//
1632//
1634//Transform3D::Transform3D(EulerType euler_type, const Dict& rotation) //YYY
1635//{
1636// init();
1637// set_rotation(euler_type,rotation);
1638//}
1639//
1640//
1642//
1643//Transform3D::Transform3D( const Vec3f& pretrans, const float& az, const float& alt, const float& phi, const Vec3f& posttrans ) //YYY by default EMAN
1644//{
1645// init();
1646// set_pretrans(pretrans);
1647// set_rotation(az,alt,phi);
1648// set_posttrans(posttrans);
1649//}
1650//
1651//
1652//
1653//
1654//Transform3D::~Transform3D()
1655//{
1656//}
1657//
1658//
1659//
1660//void Transform3D::to_identity()
1661//{
1665//
1666// for(int i=0; i<4; ++i) {
1667// for(int j=0; j<4; ++j) {
1668// if(i==j) {
1669// matrix[i][j] = 1;
1670// }
1671// else {
1672// matrix[i][j] = 0;
1673// }
1674// }
1675// }
1676// post_x_mirror = false;
1677// set_center(Vec3f(0,0,0));
1678//}
1679//
1680//
1681//
1682//bool Transform3D::is_identity() // YYY
1683//{
1684// for (int i=0; i<4; i++) {
1685// for (int j=0; j<4; j++) {
1686// if (i==j && matrix[i][j]!=1.0) return 0;
1687// if (i!=j && matrix[i][j]!=0.0) return 0;
1688// }
1689// }
1690// return 1;
1691//}
1692//
1693//
1694//void Transform3D::set_center(const Vec3f & center) //YYN
1695//{
1696// set_pretrans( Vec3f(0,0,0)-center);
1697// for (int i = 0; i < 3; i++) {
1698// matrix[i][3]=center[i];
1699// }
1700//}
1701//
1704//void Transform3D::init() // M1
1705//{
1706// to_identity();
1707//}
1708//
1710//
1711//void Transform3D::set_pretrans(const float& dx, const float& dy, const float& dz) // YYY
1712//{ set_pretrans( Vec3f(dx,dy,dz)); }
1713//
1714//
1715//void Transform3D::set_pretrans(const float& dx, const float& dy) // YYY
1716//{ set_pretrans( Vec3f(dx,dy,0)); }
1717//
1718//void Transform3D::set_pretrans(const Vec2f& pretrans) // YYY
1719//{ set_pretrans( Vec3f(pretrans[0],pretrans[1],0)); }
1720//
1721//void Transform3D::set_pretrans(const Vec3f & preT) // flag=1 means keep the old value of total trans
1722//{
1723// int flag=0;
1724//
1727// if (flag==0){
1728// matrix[0][3] = matrix[3][0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2] ;
1729// matrix[1][3] = matrix[3][1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2] ;
1730// matrix[2][3] = matrix[3][2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2] ;
1731// }
1733// if (flag==1){
1734// matrix[3][0] = matrix[0][3] - (matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2]) ;
1735// matrix[3][1] = matrix[1][3] - (matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2]) ;
1736// matrix[3][2] = matrix[2][3] - (matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2]) ;
1737// }
1738//}
1739//
1740//
1741//void Transform3D::set_posttrans(const float& dx, const float& dy, const float& dz) // YYY
1742//{ set_posttrans( Vec3f(dx,dy,dz)); }
1743//
1744//
1745//void Transform3D::set_posttrans(const float& dx, const float& dy) // YYY
1746//{ set_posttrans( Vec3f(dx,dy,0)); }
1747//
1748//void Transform3D::set_posttrans(const Vec2f& posttrans) // YYY
1749//{ set_pretrans( Vec3f(posttrans[0],posttrans[1],0)); }
1750//
1751//void Transform3D::set_posttrans(const Vec3f & posttrans) // flag=1 means keep the old value of total trans
1752//{
1753// int flag=0;
1754// Vec3f preT = get_pretrans(0) ;
1755// for (int i = 0; i < 3; i++) {
1756// matrix[3][i] = posttrans[i];
1757// }
1760// if (flag==0) {
1761// matrix[0][3] = matrix[3][0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2] ;
1762// matrix[1][3] = matrix[3][1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2] ;
1763// matrix[2][3] = matrix[3][2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2] ;
1764// }
1766// if (flag==1) { // Don't do anything
1767// }
1768//}
1769//
1770//
1771//
1772//
1773//void Transform3D::apply_scale(const float& scale) // YYY
1774//{
1775// for (int i = 0; i < 3; i++) {
1776// for (int j = 0; j < 4; j++) {
1777// matrix[i][j] *= scale;
1778// }
1779// }
1780// for (int j = 0; j < 3; j++) {
1781// matrix[3][j] *= scale;
1782// }
1783//}
1784//
1785//void Transform3D::orthogonalize() // YYY
1786//{
1787// //EulerType EMAN;
1788// float scale = get_scale() ;
1789// float inverseScale= 1/scale ;
1790// apply_scale(inverseScale);
1793//}
1794//
1795//
1796//void Transform3D::transpose() // YYY
1797//{
1798// float tempij;
1799// for (int i = 0; i < 3; i++) {
1800// for (int j = 0; j < i; j++) {
1801// tempij= matrix[i][j];
1802// matrix[i][j] = matrix[j][i];
1803// matrix[j][i] = tempij;
1804// }
1805// }
1806//}
1807//
1808//void Transform3D::set_scale(const float& scale) // YYY
1809//{
1810// float OldScale= get_scale();
1811// float Scale2Apply = scale/OldScale;
1812// apply_scale(Scale2Apply);
1813//}
1814//
1815//float Transform3D::get_mag() const //
1816//{
1817// EulerType eulertype= SPIN ;
1818// Dict AA= get_rotation(eulertype);
1819// return AA["omega"];
1820//}
1821//
1822//Vec3f Transform3D::get_finger() const //
1823//{
1824// EulerType eulertype= SPIN ;
1825// Dict AA= get_rotation(eulertype);
1826// return Vec3f(AA["n1"],AA["n2"],AA["n3"]);
1827//}
1828//
1829//Vec3f Transform3D::get_posttrans(int flag) const //
1830//{
1831// if (flag==0){
1832// return Vec3f(matrix[3][0], matrix[3][1], matrix[3][2]);
1833// }
1834// // otherwise as if all the translation was post
1835// return Vec3f(matrix[0][3], matrix[1][3], matrix[2][3]);
1836//}
1837//
1838//Vec3f Transform3D::get_total_posttrans() const {
1839// return get_posttrans(1);
1840//}
1841//
1842//Vec3f Transform3D::get_total_pretrans() const {
1843// return get_pretrans(1);
1844//}
1845//
1846//
1847//Vec3f Transform3D::get_pretrans(int flag) const // Fix Me
1848//{
1850//
1851// Vec3f pretrans;
1852// Vec3f posttrans(matrix[3][0], matrix[3][1], matrix[3][2]);
1853// Vec3f tottrans(matrix[0][3], matrix[1][3], matrix[2][3]);
1854// Vec3f totminuspost;
1855//
1856// totminuspost = tottrans;
1857// if (flag==0) {
1858// totminuspost = tottrans-posttrans;
1859// }
1860//
1861// Transform3D Rinv = inverse();
1862// for (int i=0; i<3; i++) {
1863// float ptnow=0;
1864// for (int j=0; j<3; j++) {
1865// ptnow += Rinv.matrix[i][j]* totminuspost[j] ;
1866// }
1867// pretrans.set_value_at(i,ptnow) ; //
1868// }
1869// return pretrans;
1870//}
1871//
1872//
1873// Vec3f Transform3D::get_center() const // YYY
1874// {
1875// return Vec3f();
1876// }
1877//
1878//
1879//
1880//Vec3f Transform3D::get_matrix3_col(int i) const // YYY
1881//{
1882// return Vec3f(matrix[0][i], matrix[1][i], matrix[2][i]);
1883//}
1884//
1885//
1886//Vec3f Transform3D::get_matrix3_row(int i) const // YYY
1887//{
1888// return Vec3f(matrix[i][0], matrix[i][1], matrix[i][2]);
1889//}
1890//
1891//Vec3f Transform3D::transform(const Vec3f & v3f) const // YYY
1892//{
1894// float x = matrix[0][0] * v3f[0] + matrix[0][1] * v3f[1] + matrix[0][2] * v3f[2] + matrix[0][3] ;
1895// float y = matrix[1][0] * v3f[0] + matrix[1][1] * v3f[1] + matrix[1][2] * v3f[2] + matrix[1][3] ;
1896// float z = matrix[2][0] * v3f[0] + matrix[2][1] * v3f[1] + matrix[2][2] * v3f[2] + matrix[2][3] ;
1897// return Vec3f(x, y, z);
1898//}
1899//
1900//
1901//Vec3f Transform3D::rotate(const Vec3f & v3f) const // YYY
1902//{
1904// float x = matrix[0][0] * v3f[0] + matrix[0][1] * v3f[1] + matrix[0][2] * v3f[2] ;
1905// float y = matrix[1][0] * v3f[0] + matrix[1][1] * v3f[1] + matrix[1][2] * v3f[2] ;
1906// float z = matrix[2][0] * v3f[0] + matrix[2][1] * v3f[1] + matrix[2][2] * v3f[2] ;
1907// return Vec3f(x, y, z);
1908//}
1909//
1910//
1911//Transform3D EMAN::operator*(const Transform3D & M2, const Transform3D & M1) // YYY
1912//{
1915// Transform3D resultant;
1916// for (int i=0; i<3; i++) {
1917// for (int j=0; j<4; j++) {
1918// resultant[i][j] = M2[i][0] * M1[0][j] + M2[i][1] * M1[1][j] + M2[i][2] * M1[2][j];
1919// }
1920// resultant[i][3] += M2[i][3]; // add on the new translation (not included above)
1921// }
1922//
1923// for (int j=0; j<3; j++) {
1924// resultant[3][j] = M2[3][j];
1925// }
1926//
1927// return resultant; // This will have the post_trans of M2
1928//}
1929
1930/* Here starts the pure rotation stuff */
1931
2001//
2002//void Transform3D::set_rotation(const float& az, const float& alt, const float& phi )
2003//{
2004// EulerType euler_type=EMAN;
2005// Dict rot;
2006// rot["az"] = az;
2007// rot["alt"] = alt;
2008// rot["phi"] = phi;
2009// set_rotation(euler_type, rot);
2010//}
2011//
2013//void Transform3D::set_rotation(EulerType euler_type, const float& a1, const float& a2, const float& a3)
2014//{
2015// init();
2016// Dict rot;
2017// switch(euler_type) {
2018// case EMAN:
2019// rot["az"] = a1;
2020// rot["alt"] = a2;
2021// rot["phi"] = a3;
2022// break;
2023// case SPIDER:
2024// rot["phi"] = a1;
2025// rot["theta"] = a2;
2026// rot["psi"] = a3;
2027// break;
2028// case IMAGIC:
2029// rot["alpha"] = a1;
2030// rot["beta"] = a2;
2031// rot["gamma"] = a3;
2032// break;
2033// case MRC:
2034// rot["phi"] = a1;
2035// rot["theta"] = a2;
2036// rot["omega"] = a3;
2037// break;
2038// case XYZ:
2039// rot["xtilt"] = a1;
2040// rot["ytilt"] = a2;
2041// rot["ztilt"] = a3;
2042// break;
2043// default:
2044// throw InvalidValueException(euler_type, "cannot instantiate this Euler Type");
2045// } // ends switch euler_type
2046// set_rotation(euler_type, rot);
2047//}
2048//
2050//void Transform3D::set_rotation(EulerType euler_type, const float& a1, const float& a2, const float& a3, const float& a4)
2051//{
2052// init();
2053// Dict rot;
2054// switch(euler_type) {
2055// case QUATERNION:
2056// rot["e0"] = a1;
2057// rot["e1"] = a2;
2058// rot["e2"] = a3;
2059// rot["e3"] = a4;
2060// break;
2061// case SGIROT:
2062// rot["q"] = a1;
2063// rot["n1"] = a2;
2064// rot["n2"] = a3;
2065// rot["n3"] = a4;
2066// case SPIN:
2067// rot["omega"] = a1;
2068// rot["n1"] = a2;
2069// rot["n2"] = a3;
2070// rot["n3"] = a4;
2071// break;
2072// default:
2073// throw InvalidValueException(euler_type, "cannot instantiate this Euler Type");
2074// } // ends switch euler_type
2075// set_rotation(euler_type, rot);
2076//}
2077//
2078//void Transform3D::set_rotation(EulerType euler_type, const Dict& rotation)
2079//{
2080// double e0 = 0; double e1=0; double e2=0; double e3=0;
2081// double omega=0;
2082// double az = 0;
2083// double alt = 0;
2084// double phi = 0;
2085// double cxtilt = 0;
2086// double sxtilt = 0;
2087// double cytilt = 0;
2088// double sytilt = 0;
2089// bool is_quaternion = 0;
2090// bool is_matrix = 0;
2091//
2092// switch(euler_type) {
2093// case EMAN:
2094// az = (double)rotation["az"] ;
2095// alt = (double)rotation["alt"] ;
2096// phi = (double)rotation["phi"] ;
2097// break;
2098// case IMAGIC:
2099// az = (double)rotation["alpha"] ;
2100// alt = (double)rotation["beta"] ;
2101// phi = (double)rotation["gamma"] ;
2102// break;
2103//
2104// case SPIDER:
2105// az = (double)rotation["phi"] + 90.0;
2106// alt = (double)rotation["theta"] ;
2107// phi = (double)rotation["psi"] - 90.0;
2108// break;
2109//
2110// case XYZ:
2111// cxtilt = cos( (M_PI/180.0f)*(double)rotation["xtilt"]);
2112// sxtilt = sin( (M_PI/180.0f)*(double)rotation["xtilt"]);
2113// cytilt = cos( (M_PI/180.0f)*(double)rotation["ytilt"]);
2114// sytilt = sin( (M_PI/180.0f)*(double)rotation["ytilt"]);
2115// az = (180.0f/M_PI)*atan2(-cytilt*sxtilt,sytilt) + 90.0f ;
2116// alt = (180.0f/M_PI)*acos(cytilt*cxtilt) ;
2117// phi = (double)rotation["ztilt"] +(180.0f/M_PI)*atan2(sxtilt,cxtilt*sytilt) - 90.0f ;
2118// break;
2119//
2120// case MRC:
2121// az = (double)rotation["phi"] + 90.0f ;
2122// alt = (double)rotation["theta"] ;
2123// phi = (double)rotation["omega"] - 90.0f ;
2124// break;
2125//
2126// case QUATERNION:
2127// is_quaternion = 1;
2128// e0 = (double)rotation["e0"];
2129// e1 = (double)rotation["e1"];
2130// e2 = (double)rotation["e2"];
2131// e3 = (double)rotation["e3"];
2132// break;
2133//
2134// case SPIN:
2135// is_quaternion = 1;
2136// omega = (double)rotation["omega"];
2137// e0 = cos(omega*M_PI/360.0f);
2138// e1 = sin(omega*M_PI/360.0f)* (double)rotation["n1"];
2139// e2 = sin(omega*M_PI/360.0f)* (double)rotation["n2"];
2140// e3 = sin(omega*M_PI/360.0f)* (double)rotation["n3"];
2141// break;
2142//
2143// case SGIROT:
2144// is_quaternion = 1;
2145// omega = (double)rotation["q"] ;
2146// e0 = cos(omega*M_PI/360.0f);
2147// e1 = sin(omega*M_PI/360.0f)* (double)rotation["n1"];
2148// e2 = sin(omega*M_PI/360.0f)* (double)rotation["n2"];
2149// e3 = sin(omega*M_PI/360.0f)* (double)rotation["n3"];
2150// break;
2151//
2152// case MATRIX:
2153// is_matrix = 1;
2154// matrix[0][0] = (float)rotation["m11"] ;
2155// matrix[0][1] = (float)rotation["m12"] ;
2156// matrix[0][2] = (float)rotation["m13"] ;
2157// matrix[1][0] = (float)rotation["m21"] ;
2158// matrix[1][1] = (float)rotation["m22"] ;
2159// matrix[1][2] = (float)rotation["m23"] ;
2160// matrix[2][0] = (float)rotation["m31"] ;
2161// matrix[2][1] = (float)rotation["m32"] ;
2162// matrix[2][2] = (float)rotation["m33"] ;
2163// break;
2164//
2165// default:
2166// throw InvalidValueException(euler_type, "unknown Euler Type");
2167// } // ends switch euler_type
2168//
2169//
2170// Vec3f postT = get_posttrans( ) ;
2171// Vec3f preT = get_pretrans( ) ;
2172//
2173//
2174// double azp = fmod(az,360.0)*M_PI/180.0;
2175// double altp = alt*M_PI/180.0;
2176// double phip = fmod(phi,360.0)*M_PI/180.0;
2177//
2178// if (!is_quaternion && !is_matrix) {
2179// matrix[0][0] = (float)(cos(phip)*cos(azp) - cos(altp)*sin(azp)*sin(phip));
2180// matrix[0][1] = (float)(cos(phip)*sin(azp) + cos(altp)*cos(azp)*sin(phip));
2181// matrix[0][2] = (float)(sin(altp)*sin(phip));
2182// matrix[1][0] = (float)(-sin(phip)*cos(azp) - cos(altp)*sin(azp)*cos(phip));
2183// matrix[1][1] = (float)(-sin(phip)*sin(azp) + cos(altp)*cos(azp)*cos(phip));
2184// matrix[1][2] = (float)(sin(altp)*cos(phip));
2185// matrix[2][0] = (float)(sin(altp)*sin(azp));
2186// matrix[2][1] = (float)(-sin(altp)*cos(azp));
2187// matrix[2][2] = (float)cos(altp);
2188// }
2189// if (is_quaternion){
2190// matrix[0][0] = (float)(e0 * e0 + e1 * e1 - e2 * e2 - e3 * e3);
2191// matrix[0][1] = (float)(2.0f * (e1 * e2 + e0 * e3));
2192// matrix[0][2] = (float)(2.0f * (e1 * e3 - e0 * e2));
2193// matrix[1][0] = (float)(2.0f * (e2 * e1 - e0 * e3));
2194// matrix[1][1] = (float)(e0 * e0 - e1 * e1 + e2 * e2 - e3 * e3);
2195// matrix[1][2] = (float)(2.0f * (e2 * e3 + e0 * e1));
2196// matrix[2][0] = (float)(2.0f * (e3 * e1 + e0 * e2));
2197// matrix[2][1] = (float)(2.0f * (e3 * e2 - e0 * e1));
2198// matrix[2][2] = (float)(e0 * e0 - e1 * e1 - e2 * e2 + e3 * e3);
2199// // keep in mind matrix[0][2] is M13 gives an e0 e2 piece, etc
2200// }
2201// // Now do post and pretrans: vfinal = vpost + R vpre;
2202//
2203// matrix[0][3] = postT[0] + matrix[0][0]*preT[0] + matrix[0][1]*preT[1] + matrix[0][2]*preT[2] ;
2204// matrix[1][3] = postT[1] + matrix[1][0]*preT[0] + matrix[1][1]*preT[1] + matrix[1][2]*preT[2] ;
2205// matrix[2][3] = postT[2] + matrix[2][0]*preT[0] + matrix[2][1]*preT[1] + matrix[2][2]*preT[2] ;
2206//}
2207//
2208//
2209//void Transform3D::set_rotation(const float& m11, const float& m12, const float& m13,
2210// const float& m21, const float& m22, const float& m23,
2211// const float& m31, const float& m32, const float& m33)
2212//{
2213// EulerType euler_type = MATRIX;
2214// Dict rot;
2215// rot["m11"] = m11;
2216// rot["m12"] = m12;
2217// rot["m13"] = m13;
2218// rot["m21"] = m21;
2219// rot["m22"] = m22;
2220// rot["m23"] = m23;
2221// rot["m31"] = m31;
2222// rot["m32"] = m32;
2223// rot["m33"] = m33;
2224// set_rotation(euler_type, rot); // Or should it be &rot ?
2225//}
2226//
2227//void Transform3D::set_rotation(const Vec3f & eahat, const Vec3f & ebhat,
2228// const Vec3f & eAhat, const Vec3f & eBhat)
2229//{// this rotation rotates unit vectors a,b into A,B;
2231// Vec3f eahatcp(eahat);
2232// Vec3f ebhatcp(ebhat);
2233// Vec3f eAhatcp(eAhat);
2234// Vec3f eBhatcp(eBhat);
2235//
2236// eahatcp.normalize();
2237// ebhatcp.normalize();
2238// eAhatcp.normalize();
2239// eBhatcp.normalize();
2240//
2241// Vec3f aMinusA(eahatcp);
2242// aMinusA -= eAhatcp;
2243// Vec3f bMinusB(ebhatcp);
2244// bMinusB -= eBhatcp;
2245//
2246//
2247// Vec3f nhat;
2248// float aAlength = aMinusA.length();
2249// float bBlength = bMinusB.length();
2250// if (aAlength==0){
2251// nhat=eahatcp;
2252// }else if (bBlength==0){
2253// nhat=ebhatcp;
2254// }else{
2255// nhat= aMinusA.cross(bMinusB);
2256// nhat.normalize();
2257// }
2258//
2260//
2261// Vec3f neahat = eahatcp.cross(nhat);
2262// Vec3f nebhat = ebhatcp.cross(nhat);
2263// Vec3f neAhat = eAhatcp.cross(nhat);
2264// Vec3f neBhat = eBhatcp.cross(nhat);
2265//
2266// double cosomegaA = (neahat.dot(neAhat)) / (neahat.dot(neahat));
2268// double sinomegaA = (neahat.dot(eAhatcp)) / (neahat.dot(neahat));
2270//
2271// double omegaA = atan2(sinomegaA,cosomegaA);
2273//
2274// EulerType euler_type=SPIN;
2275// Dict rotation;
2276// rotation["n1"]= nhat[0];
2277// rotation["n2"]= nhat[1];
2278// rotation["n3"]= nhat[2];
2279// rotation["omega"] = (float)(omegaA*180.0/M_PI);
2280// set_rotation(euler_type, rotation);
2281//}
2282//
2283//
2284//float Transform3D::get_scale() const // YYY
2285//{
2286// // Assumes uniform scaling, calculation uses Z only.
2287// float scale =0;
2288// for (int i=0; i<3; i++) {
2289// for (int j=0; j<3; j++) {
2290// scale = scale + matrix[i][j]*matrix[i][j];
2291// }
2292// }
2293//
2294// return sqrt(scale/3);
2295//}
2296//
2297//
2298//
2299//Dict Transform3D::get_rotation(EulerType euler_type) const
2300//{
2301// Dict result;
2302//
2303// double max = 1 - ERR_LIMIT;
2304// double sca=get_scale();
2305// double cosalt=matrix[2][2]/sca;
2306//
2307//
2308// double az=0;
2309// double alt = 0;
2310// double phi=0;
2311// double phiS = 0; // like az (but in SPIDER ZXZ)
2312// double psiS =0; // like phi (but in SPIDER ZYZ)
2313//
2314//
2316//
2317// if (cosalt > max) { // that is, alt close to 0
2318// alt = 0;
2319// az=0;
2320// phi = (double)EMConsts::rad2deg * atan2(matrix[0][1], matrix[0][0]);
2321// }
2322// else if (cosalt < -max) { // alt close to pi
2323// alt = 180;
2324// az=0;
2325// phi=360.0f-(double)EMConsts::rad2deg * atan2(matrix[0][1], matrix[0][0]);
2326// }
2327// else {
2328// alt = (double)EMConsts::rad2deg * acos(cosalt);
2329// az = 360.0f+(double)EMConsts::rad2deg * atan2(matrix[2][0], -matrix[2][1]);
2330// phi = 360.0f+(double)EMConsts::rad2deg * atan2(matrix[0][2], matrix[1][2]);
2331// }
2332// az =fmod(az+180.0,360.0)-180.0;
2333// phi=fmod(phi+180.0,360.0)-180.0;
2334//
2336// if (fabs(cosalt) > max) { // that is, alt close to 0
2337// phiS=0;
2338// psiS = az+phi;
2339// }
2340// else {
2341// phiS = az - 90.0;
2342// psiS = phi + 90.0;
2343// }
2344// phiS = fmod((phiS + 360.0 ), 360.0) ;
2345// psiS = fmod((psiS + 360.0 ), 360.0) ;
2346//
2348//
2349// double nphi = (az-phi)/2.0;
2350// // The next is also e0
2351// double cosOover2 = (cos((az+phi)*M_PI/360) * cos(alt*M_PI/360)) ;
2352// double sinOover2 = sqrt(1 -cosOover2*cosOover2);
2353// double cosnTheta = sin((az+phi)*M_PI/360) * cos(alt*M_PI/360) / sqrt(1-cosOover2*cosOover2) ;
2354// double sinnTheta = sqrt(1-cosnTheta*cosnTheta);
2355// double n1 = sinnTheta*cos(nphi*M_PI/180);
2356// double n2 = sinnTheta*sin(nphi*M_PI/180);
2357// double n3 = cosnTheta;
2358// double xtilt = 0;
2359// double ytilt = 0;
2360// double ztilt = 0;
2361//
2362//
2363// if (cosOover2<0) {
2364// cosOover2*=-1; n1 *=-1; n2*=-1; n3*=-1;
2365// }
2366//
2367//
2368// switch (euler_type) {
2369// case EMAN:
2370// result["az"] = az;
2371// result["alt"] = alt;
2372// result["phi"] = phi;
2373// break;
2374//
2375// case IMAGIC:
2376// result["alpha"] = az;
2377// result["beta"] = alt;
2378// result["gamma"] = phi;
2379// break;
2380//
2381// case SPIDER:
2382// result["phi"] = phiS; // The first Euler like az
2383// result["theta"] = alt;
2384// result["psi"] = psiS;
2385// break;
2386//
2387// case MRC:
2388// result["phi"] = phiS;
2389// result["theta"] = alt;
2390// result["omega"] = psiS;
2391// break;
2392//
2393// case XYZ:
2394// xtilt = atan2(-sin((M_PI/180.0f)*phiS)*sin((M_PI/180.0f)*alt),cos((M_PI/180.0f)*alt));
2395// ytilt = asin( cos((M_PI/180.0f)*phiS)*sin((M_PI/180.0f)*alt));
2396// ztilt = psiS*M_PI/180.0f - atan2(sin(xtilt), cos(xtilt) *sin(ytilt));
2397//
2398// xtilt=fmod(xtilt*180/M_PI+540.0,360.0) -180.0;
2399// ztilt=fmod(ztilt*180/M_PI+540.0,360.0) -180.0;
2400//
2401// result["xtilt"] = xtilt;
2402// result["ytilt"] = ytilt*180/M_PI;
2403// result["ztilt"] = ztilt;
2404// break;
2405//
2406// case QUATERNION:
2407// result["e0"] = cosOover2 ;
2408// result["e1"] = sinOover2 * n1 ;
2409// result["e2"] = sinOover2 * n2;
2410// result["e3"] = sinOover2 * n3;
2411// break;
2412//
2413// case SPIN:
2414// result["omega"] =360.0f* acos(cosOover2)/ M_PI ;
2415// result["n1"] = n1;
2416// result["n2"] = n2;
2417// result["n3"] = n3;
2418// break;
2419//
2420// case SGIROT:
2421// result["q"] = 360.0f*acos(cosOover2)/M_PI ;
2422// result["n1"] = n1;
2423// result["n2"] = n2;
2424// result["n3"] = n3;
2425// break;
2426//
2427// case MATRIX:
2428// result["m11"] = matrix[0][0] ;
2429// result["m12"] = matrix[0][1] ;
2430// result["m13"] = matrix[0][2] ;
2431// result["m21"] = matrix[1][0] ;
2432// result["m22"] = matrix[1][1] ;
2433// result["m23"] = matrix[1][2] ;
2434// result["m31"] = matrix[2][0] ;
2435// result["m32"] = matrix[2][1] ;
2436// result["m33"] = matrix[2][2] ;
2437// break;
2438//
2439// default:
2440// throw InvalidValueException(euler_type, "unknown Euler Type");
2441// }
2442//
2443// return result;
2444//}
2445//
2446//Transform3D Transform3D::inverseUsingAngs() const // YYN need to test it for sure
2447//{
2448// // First Find the scale
2449// EulerType eE=EMAN;
2450//
2451//
2452// float scale = get_scale();
2453// Vec3f preT = get_pretrans( ) ;
2454// Vec3f postT = get_posttrans( ) ;
2455// Dict angs = get_rotation(eE);
2456// Dict invAngs ;
2457//
2458// invAngs["phi"] = 180.0f - (float) angs["az"] ;
2459// invAngs["az"] = 180.0f - (float) angs["phi"] ;
2460// invAngs["alt"] = angs["alt"] ;
2461//
2468//
2469// float inverseScale= 1/scale ;
2470//
2471// Transform3D invM;
2472//
2473// invM.set_rotation(EMAN, invAngs);
2474// invM.apply_scale(inverseScale);
2475// invM.set_pretrans(-postT );
2476// invM.set_posttrans(-preT );
2477//
2478//
2479// return invM;
2480//
2481//}
2482//
2483//Transform3D Transform3D::inverse() const // YYN need to test it for sure
2484//{
2485// // This assumes the matrix is 4 by 4 and the last row reads [0 0 0 1]
2486//
2487// double m00 = matrix[0][0]; double m01=matrix[0][1]; double m02=matrix[0][2];
2488// double m10 = matrix[1][0]; double m11=matrix[1][1]; double m12=matrix[1][2];
2489// double m20 = matrix[2][0]; double m21=matrix[2][1]; double m22=matrix[2][2];
2490// double v0 = matrix[0][3]; double v1 =matrix[1][3]; double v2 =matrix[2][3];
2491//
2492// double cof00 = m11*m22-m12*m21;
2493// double cof11 = m22*m00-m20*m02;
2494// double cof22 = m00*m11-m01*m10;
2495// double cof01 = m10*m22-m20*m12;
2496// double cof02 = m10*m21-m20*m11;
2497// double cof12 = m00*m21-m01*m20;
2498// double cof10 = m01*m22-m02*m21;
2499// double cof20 = m01*m12-m02*m11;
2500// double cof21 = m00*m12-m10*m02;
2501//
2502// double Det = m00* cof00 + m02* cof02 -m01*cof01;
2503//
2504// Transform3D invM;
2505//
2506// invM.matrix[0][0] = (float)(cof00/Det);
2507// invM.matrix[0][1] = (float)(-cof10/Det);
2508// invM.matrix[0][2] = (float)(cof20/Det);
2509// invM.matrix[1][0] = (float)(-cof01/Det);
2510// invM.matrix[1][1] = (float)(cof11/Det);
2511// invM.matrix[1][2] = (float)(-cof21/Det);
2512// invM.matrix[2][0] = (float)(cof02/Det);
2513// invM.matrix[2][1] = (float)(-cof12/Det);
2514// invM.matrix[2][2] = (float)(cof22/Det);
2515//
2516// invM.matrix[0][3] = (float)((-cof00*v0 + cof10*v1 - cof20*v2)/Det);
2517// invM.matrix[1][3] = (float)(( cof01*v0 - cof11*v1 + cof21*v2)/Det);
2518// invM.matrix[2][3] = (float)((-cof02*v0 + cof12*v1 - cof22*v2)/Det);
2519//
2520// Vec3f postT = get_posttrans( ) ;
2521// Vec3f invMpre = - postT;
2522// Vec3f invMpost ;
2523// for ( int i = 0; i < 3; i++) {
2524// invMpost[i] = invM.matrix[i][3];
2525// for ( int j = 0; j < 3; j++) {
2526// invMpost[i] += - invM.matrix[i][j]*invMpre[j];
2527// }
2528// invM.matrix[3][i] = invMpost[i];
2529// }
2530//
2531// return invM;
2532//}
2533//
2534//
2535//
2537//
2538//Transform3D Transform3D::get_sym(const string & symname, int n) const
2539//{
2540// int nsym = get_nsym(symname);
2541//
2544//
2545// // see www.math.utah.edu/~alfeld/math/polyhedra/polyhedra.html for pictures
2546// // By default we will put largest symmetry along z-axis.
2547//
2548// // Each Platonic Solid has 2E symmetry elements.
2549//
2550//
2551// // An icosahedron has m=5, n=3, F=20 E=30=nF/2, V=12=nF/m,since vertices shared by 5 triangles;
2552// // It is composed of 20 triangles. E=3*20/2;
2553//
2554//
2555// // An dodecahedron has m=3, n=5 F=12 E=30 V=20
2556// // It is composed of 12 pentagons. E=5*12/2; V= 5*12/3, since vertices shared by 3 pentagons;
2557//
2558//
2559//
2560// // The ICOS symmetry group has the face along z-axis
2561//
2562// float lvl0=0; // there is one pentagon on top; five-fold along z
2563// float lvl1= 63.4349f; // that is atan(2) // there are 5 pentagons with centers at this height (angle)
2564// float lvl2=116.5651f; //that is 180-lvl1 // there are 5 pentagons with centers at this height (angle)
2565// float lvl3=180.f; // there is one pentagon on the bottom
2566// // Notice that 63.439 is the angle between two faces of the dual object
2567//
2568// static double ICOS[180] = { // This is with a pentagon normal to z
2569// 0,lvl0,0, 0,lvl0,288, 0,lvl0,216, 0,lvl0,144, 0,lvl0,72,
2570// 0,lvl1,36, 0,lvl1,324, 0,lvl1,252, 0,lvl1,180, 0,lvl1,108,
2571// 72,lvl1,36, 72,lvl1,324, 72,lvl1,252, 72,lvl1,180, 72,lvl1,108,
2572// 144,lvl1,36, 144,lvl1,324, 144,lvl1,252, 144,lvl1,180, 144,lvl1,108,
2573// 216,lvl1,36, 216,lvl1,324, 216,lvl1,252, 216,lvl1,180, 216,lvl1,108,
2574// 288,lvl1,36, 288,lvl1,324, 288,lvl1,252, 288,lvl1,180, 288,lvl1,108,
2575// 36,lvl2,0, 36,lvl2,288, 36,lvl2,216, 36,lvl2,144, 36,lvl2,72,
2576// 108,lvl2,0, 108,lvl2,288, 108,lvl2,216, 108,lvl2,144, 108,lvl2,72,
2577// 180,lvl2,0, 180,lvl2,288, 180,lvl2,216, 180,lvl2,144, 180,lvl2,72,
2578// 252,lvl2,0, 252,lvl2,288, 252,lvl2,216, 252,lvl2,144, 252,lvl2,72,
2579// 324,lvl2,0, 324,lvl2,288, 324,lvl2,216, 324,lvl2,144, 324,lvl2,72,
2580// 0,lvl3,0, 0,lvl3,288, 0,lvl3,216, 0,lvl3,144, 0,lvl3,72
2581// };
2582//
2583//
2584// // A cube has m=3, n=4, F=6 E=12=nF/2, V=8=nF/m,since vertices shared by 3 squares;
2585// // It is composed of 6 squares.
2586//
2587//
2588// // An octahedron has m=4, n=3, F=8 E=12=nF/2, V=6=nF/m,since vertices shared by 4 triangles;
2589// // It is composed of 8 triangles.
2590//
2591// // We have placed the OCT symmetry group with a face along the z-axis
2592// lvl0=0;
2593// lvl1=90;
2594// lvl2=180;
2595//
2596// static float OCT[72] = {// This is with a face of a cube along z
2597// 0,lvl0,0, 0,lvl0,90, 0,lvl0,180, 0,lvl0,270,
2598// 0,lvl1,0, 0,lvl1,90, 0,lvl1,180, 0,lvl1,270,
2599// 90,lvl1,0, 90,lvl1,90, 90,lvl1,180, 90,lvl1,270,
2600// 180,lvl1,0, 180,lvl1,90, 180,lvl1,180, 180,lvl1,270,
2601// 270,lvl1,0, 270,lvl1,90, 270,lvl1,180, 270,lvl1,270,
2602// 0,lvl2,0, 0,lvl2,90, 0,lvl2,180, 0,lvl2,270
2603// };
2604// // B^4=A^3=1; BABA=1; implies AA=BAB, ABA=B^3 , AB^2A = BBBABBB and
2605// // 20 words with at most a single A
2606// // 1 B BB BBB A BA AB BBA BAB ABB BBBA BBAB BABB ABBB BBBAB BBABB BABBB
2607// // BBBABB BBABBB BBBABBB
2608// // also ABBBA is distinct yields 4 more words
2609// // ABBBA BABBBA BBABBBA BBBABBBA
2610// // for a total of 24 words
2611// // Note A BBB A BBB A reduces to BBABB
2612// // and B A BBB A is the same as A BBB A BBB etc.
2613//
2614// // The TET symmetry group has a face along the z-axis
2615// // It has n=m=3; F=4, E=6=nF/2, V=4=nF/m
2616// lvl0=0; // There is a face along z
2617// lvl1=109.4712f; // that is acos(-1/3) // There are 3 faces at this angle
2618//
2619// static float TET[36] = {// This is with the face along z
2620// 0,lvl0,0, 0,lvl0,120, 0,lvl0,240,
2621// 0,lvl1,60, 0,lvl1,180, 0,lvl1,300,
2622// 120,lvl1,60, 120,lvl1,180, 120,lvl1,300,
2623// 240,lvl1,60, 240,lvl1,180, 240,lvl1,300
2624// };
2625// // B^3=A^3=1; BABA=1; implies A^2=BAB, ABA=B^2 , AB^2A = B^2AB^2 and
2626// // 12 words with at most a single A
2627// // 1 B BB A BA AB BBA BAB ABB BBAB BABB BBABB
2628// // at most one A is necessary
2629//
2630// Transform3D ret;
2631// SymType type = get_sym_type(symname);
2632//
2633// switch (type) {
2634// case CSYM:
2635// ret.set_rotation( n * 360.0f / nsym, 0, 0);
2636// break;
2637// case DSYM:
2638// if (n >= nsym / 2) {
2639// ret.set_rotation((n - nsym/2) * 360.0f / (nsym / 2),180.0f, 0);
2640// }
2641// else {
2642// ret.set_rotation( n * 360.0f / (nsym / 2),0, 0);
2643// }
2644// break;
2645// case ICOS_SYM:
2646// ret.set_rotation((float)ICOS[n * 3 ],
2647// (float)ICOS[n * 3 + 1],
2648// (float)ICOS[n * 3 + 2] );
2649// break;
2650// case OCT_SYM:
2651// ret.set_rotation((float)OCT[n * 3],
2652// (float)OCT[n * 3 + 1],
2653// (float)OCT[n * 3 + 2] );
2654// break;
2655// case TET_SYM:
2656// ret.set_rotation((float)TET[n * 3 ],
2657// (float)TET[n * 3 + 1] ,
2658// (float)TET[n * 3 + 2] );
2659// break;
2660// case ISYM:
2661// ret.set_rotation(0, 0, 0);
2662// break;
2663// default:
2664// throw InvalidValueException(type, symname);
2665// }
2666//
2667// ret = (*this) * ret;
2668//
2669// return ret;
2670//}
2671//
2672//int Transform3D::get_nsym(const string & name)
2673//{
2674// string symname = name;
2675//
2676// for (size_t i = 0; i < name.size(); i++) {
2677// if (isalpha(name[i])) {
2678// symname[i] = (char)tolower(name[i]);
2679// }
2680// }
2681//
2682// SymType type = get_sym_type(symname);
2683// int nsym = 0;
2684//
2685// switch (type) {
2686// case CSYM:
2687// nsym = atoi(symname.c_str() + 1);
2688// break;
2689// case DSYM:
2690// nsym = atoi(symname.c_str() + 1) * 2;
2691// break;
2692// case ICOS_SYM:
2693// nsym = 60;
2694// break;
2695// case OCT_SYM:
2696// nsym = 24;
2697// break;
2698// case TET_SYM:
2699// nsym = 12;
2700// break;
2701// case ISYM:
2702// nsym = 1;
2703// break;
2704// case UNKNOWN_SYM:
2705// default:
2706// throw InvalidValueException(type, name);
2707// }
2708// return nsym;
2709//}
2710//
2711//
2712//
2713//Transform3D::SymType Transform3D::get_sym_type(const string & name)
2714//{
2715// SymType t = UNKNOWN_SYM;
2716//
2717// if (name[0] == 'c') {
2718// t = CSYM;
2719// }
2720// else if (name[0] == 'd') {
2721// t = DSYM;
2722// }
2723// else if (name == "icos") {
2724// t = ICOS_SYM;
2725// }
2726// else if (name == "oct") {
2727// t = OCT_SYM;
2728// }
2729// else if (name == "tet") {
2730// t = TET_SYM;
2731// }
2732// else if (name == "i" || name == "") {
2733// t = ISYM;
2734// }
2735// return t;
2736//}
2737//
2738//vector<Transform3D*>
2739//Transform3D::angles2tfvec(EulerType eulertype, const vector<float> ang) {
2740// int nangles = ang.size() / 3;
2741// vector<Transform3D*> tfvec;
2742// for (int i = 0; i < nangles; i++) {
2743// tfvec.push_back(new Transform3D(eulertype,ang[3*i],ang[3*i+1],ang[3*i+2]));
2744// }
2745// return tfvec;
2746//}
2747//
2748
2749
2750
2751/* Rotation stuff */
2752
V
#define V(i, j)
Definition: analyzer.cpp:697
v0
#define v0(i)
Definition: analyzer.cpp:698
EMAN::Dict::const_iterator
Const iterator support for the Dict object This is just a wrapper, everything is inherited from the m...
Definition: emobject.h:674
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::has_key_ci
bool has_key_ci(const string &key) const
Ask the Dictionary if it as a particular key in a case insensitive way.
Definition: emobject.cpp:1140
EMAN::Dict::end
iterator end()
Definition: emobject.cpp:1061
EMAN::Dict::begin
iterator begin()
Definition: emobject.cpp:1045
EMAN::Dict::get_ci
EMObject get_ci(const string &key) const
Get the EMObject corresponding to the particular key using case insensitivity.
Definition: emobject.cpp:1128
EMAN::EMConsts::rad2deg
static const double rad2deg
Definition: emobject.h:77
EMAN::EMConsts::deg2rad
static const double deg2rad
Definition: emobject.h:78
EMAN::EMObject
EMObject is a wrapper class for types including int, float, double, etc as defined in ObjectType.
Definition: emobject.h:123
EMAN::EMObject::get_type
ObjectType get_type() const
Get the ObjectType This is an enumerated type first declared in the class EMObjectTypes.
Definition: emobject.h:224
EMAN::Factory::get
static T * get(const string &instance_name)
Definition: emobject.h:781
EMAN::Symmetry3D
Symmetry3D - A base class for 3D Symmetry objects.
Definition: symmetry.h:57
EMAN::Symmetry3D::get_nsym
virtual int get_nsym() const =0
The total number of unique symmetry operations that will be return by this object when a calling prog...
EMAN::Symmetry3D::get_sym
virtual Transform get_sym(const int n) const =0
Every Symmetry3D object must provide access to the full set of its symmetry operators via this functi...
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::get_params
Dict get_params(const string &euler_type) const
Get the parameters of the entire transform, using a specific euler convention.
Definition: transform.cpp:479
EMAN::Transform::get_hflip_transform
Transform get_hflip_transform() const
How do I get the transform that will yield the horizontally flipped projection?
Definition: transform.cpp:792
EMAN::Transform::operator!=
bool operator!=(const Transform &rhs) const
Unequality comparision operator.
Definition: transform.cpp:128
EMAN::Transform::get_mirror
bool get_mirror() const
Query whether x_mirroring is occurring.
Definition: transform.cpp:1250
EMAN::Transform::copy_matrix_into_array
void copy_matrix_into_array(float *const) const
Definition: transform.cpp:158
EMAN::Transform::rotate_origin_newBasis
void rotate_origin_newBasis(const Transform &tcs, const float &omega, const float &n1, const float &n2, const float &n3)
Increment the rotation by multipling the rotation bit of the argument transfrom by the rotation part ...
Definition: transform.cpp:738
EMAN::Transform::get_sym
Transform get_sym(const string &sym, int n) const
Apply the symmetry deduced from the function arguments to this Transform and return the result.
Definition: transform.cpp:1389
EMAN::Transform::is_rot_identity
bool is_rot_identity() const
Returns whethers or this matrix rotation is the identity.
Definition: transform.cpp:238
EMAN::Transform::assert_valid_2d
void assert_valid_2d() const
Definition: transform.cpp:1366
EMAN::Transform::matrix
float matrix[3][4]
Definition: transform.h:503
EMAN::Transform::ERR_LIMIT
static const float ERR_LIMIT
Definition: transform.h:78
EMAN::Transform::get_pre_trans_2d
Vec2f get_pre_trans_2d() const
2D version of getting the translation vector as though this object was MSRT_ not MTSR,...
Definition: transform.cpp:1111
EMAN::Transform::get_matrix_string
string get_matrix_string(int precision)
Definition: transform.cpp:169
EMAN::Transform::get_rotation_transform
Transform get_rotation_transform() const
Get the rotation part of the tranformation matrix as a Transform object.
Definition: transform.cpp:688
EMAN::Transform::operator==
bool operator==(const Transform &rhs) const
Equality comparision operator.
Definition: transform.cpp:119
EMAN::Transform::set_mirror
void set_mirror(const bool x_mirror)
Set whether or not x_mirroring is occurring.
Definition: transform.cpp:1238
EMAN::Transform::scale
void scale(const float &scale)
Increment the scale.
Definition: transform.cpp:1159
EMAN::Transform::set_rotation
void set_rotation(const Dict &rotation)
Set a rotation using a specific Euler type and the dictionary interface Works for all Euler types.
Definition: transform.cpp:519
EMAN::Transform::get_params_inverse
Dict get_params_inverse(const string &euler_type) const
Get the parameters of the inverse of the transform as though it were in RSMT order not MTSR.
Definition: transform.cpp:499
EMAN::Transform::get_pre_trans
Vec3f get_pre_trans() const
Get the translation vector as though this object was MSRT_ not MTSR, where T_ is what you want Note M...
Definition: transform.cpp:1100
EMAN::Transform::get_vflip_transform
Transform get_vflip_transform() const
How do I get the transform that will yield the vertically flipped projection?
Definition: transform.cpp:811
EMAN::Transform::permissable_3d_not_rot
static vector< string > permissable_3d_not_rot
Definition: transform.h:509
EMAN::Transform::get_sym_proj
vector< Transform > get_sym_proj(const string &sym) const
Definition: transform.cpp:1407
EMAN::Transform::get_rotation
Dict get_rotation(const string &euler_type="eman") const
Get a rotation in any Euler format.
Definition: transform.cpp:829
EMAN::Transform::get_sym_sparx
Transform get_sym_sparx(const string &sym, int n) const
Definition: transform.cpp:1398
EMAN::Transform::translate
void translate(const float &tx, const float &ty, const float &tz=0)
Increment the current translation by tx, ty, tz.
Definition: transform.cpp:1063
EMAN::Transform::permissable_rot_keys
static map< string, vector< string > > permissable_rot_keys
Definition: transform.h:510
EMAN::Transform::get_trans
Vec3f get_trans() const
Get the post trans as a vec3f.
Definition: transform.cpp:1046
EMAN::Transform::to_identity
void to_identity()
Force the internal matrix to become the identity.
Definition: transform.cpp:207
EMAN::Transform::permissable_2d_not_rot
static vector< string > permissable_2d_not_rot
This map is used to validate keys in the argument maps - e.g. if the type is 2d and the angle is not ...
Definition: transform.h:508
EMAN::Transform::set
void set(int r, int c, float value)
Set the value stored in the internal transformation matrix at at coordinate (r,c) to value.
Definition: transform.h:400
EMAN::Transform::detect_problem_keys
void detect_problem_keys(const Dict &d)
Test to ensure the parametes in the given dictionary are valid Throws if an error is detected Generic...
Definition: transform.cpp:369
EMAN::Transform::get_matrix_4x4
vector< float > get_matrix_4x4() const
Get the 4x4 transformation matrix using a vector.
Definition: transform.cpp:192
EMAN::Transform::negate
Transform negate() const
Negates the Transform - a useful way, for example, for getting an orientation on the opposite side of...
Definition: transform.cpp:781
EMAN::Transform::set_params
void set_params(const Dict &d)
Set the parameters of the entire transform.
Definition: transform.cpp:254
EMAN::Transform::Transform
Transform()
Default constructor Internal matrix is the identity.
Definition: transform.cpp:100
EMAN::Transform::rotate
void rotate(const Transform &by)
Increment the rotation by multipling the rotation bit of the argument transfrom by the current transf...
Definition: transform.cpp:763
EMAN::Transform::get_scale
float get_scale() const
Get the scale that was applied.
Definition: transform.cpp:1145
EMAN::Transform::inverse
Transform inverse() const
Get the inverse of this transformation matrix.
Definition: transform.cpp:1327
EMAN::Transform::get_scale_and_mirror
void get_scale_and_mirror(float &scale, bool &x_mirror) const
Get scale and x_mirror with 1 function call.
Definition: transform.cpp:1260
EMAN::Transform::rotate_origin
void rotate_origin(const Transform &by)
Increment the rotation by multipling the rotation bit of the argument transfrom by the rotation part ...
Definition: transform.cpp:720
EMAN::Transform::set_scale
void set_scale(const float &scale)
Set the scale.
Definition: transform.cpp:1123
EMAN::Transform::operator=
Transform & operator=(const Transform &that)
Assignment operator.
Definition: transform.cpp:110
EMAN::Transform::set_matrix
void set_matrix(const vector< float > &v)
Set the transformation matrix using a vector.
Definition: transform.cpp:147
EMAN::Transform::get_determinant
float get_determinant() const
Get the determinant of the matrix.
Definition: transform.cpp:1279
EMAN::Transform::is_identity
bool is_identity() const
Returns whethers or this matrix is the identity.
Definition: transform.cpp:222
EMAN::Transform::orthogonalize
void orthogonalize()
Reorthogonalize the rotation part of the matrix in place.
Definition: transform.cpp:1179
EMAN::Transform::set_trans
void set_trans(const float &x, const float &y, const float &z=0)
Set the post translation component.
Definition: transform.cpp:1036
EMAN::Transform::init_permissable_keys
void init_permissable_keys()
Called internally to initialize permissable_2d_not_rot, permissable_3d_not_rot, and permissable_rot_k...
Definition: transform.cpp:285
EMAN::Transform::translate_newBasis
void translate_newBasis(const Transform &tcs, const float &tx, const float &ty, const float &tz=0)
Increment the current translation by tx, ty, tz using a non standard basis Actualy what it does is re...
Definition: transform.cpp:1072
EMAN::Transform::invert
void invert()
Get the inverse of this transformation matrix.
Definition: transform.cpp:1293
EMAN::Transform::get_matrix
vector< float > get_matrix() const
Get the transformation matrix using a vector.
Definition: transform.cpp:181
EMAN::Transform::transpose_inplace
void transpose_inplace()
Get the transpose of this transformation matrix.
Definition: transform.cpp:1333
EMAN::Transform::get_trans_2d
Vec2f get_trans_2d() const
Get the degenerant 2D post trans as a vec2f.
Definition: transform.cpp:1088
EMAN::Transform::set_params_inverse
void set_params_inverse(const Dict &d)
Set the parameters of the entire transform as though they there in the inverse format.
Definition: transform.cpp:425
EMAN::Transform::tet_3_to_2
static Transform tet_3_to_2()
Get the transform that moves any tetrahedron generated by eman2 so that it matches the 2-2-2 (MRC,...
Definition: transform.cpp:85
EMAN::Transform::get_nsym
static int get_nsym(const string &sym)
get the number of symmetries associated with the given symmetry name
Definition: transform.cpp:1570
EMAN::Transform::transpose
Transform transpose() const
Get the transpose of this transformation matrix.
Definition: transform.cpp:1346
EMAN::Util::str_to_lower
static string str_to_lower(const string &s)
Return a lower case version of the argument string.
Definition: util.cpp:283
EMAN::Util::hypot3
static float hypot3(int x, int y, int z)
Euclidean distance function in 3D: f(x,y,z) = sqrt(x*x + y*y + z*z);.
Definition: util.h:827
EMAN::Util::apply_precision
static void apply_precision(float &value, const float &precision)
Definition: util.h:1309
EMAN::Vec2
The Vec2 is precisely the same as Vec3 except it works exclusively in 2D Note there are convenient ty...
Definition: vec3.h:710
EMAN::Vec3::normalize
float normalize()
Normalize the vector and return its length before the normalization.
Definition: vec3.h:332
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
InvalidParameterException
#define InvalidParameterException(desc)
Definition: exception.h:361
InvalidValueException
#define InvalidValueException(val, desc)
Definition: exception.h:285
InvalidStringException
#define InvalidStringException(str, desc)
Definition: exception.h:306
UnexpectedBehaviorException
#define UnexpectedBehaviorException(desc)
Definition: exception.h:400
EMAN
E2Exception class.
Definition: aligner.h:40
EMAN::operator*
EMData * operator*(const EMData &em, float n)
Definition: emdata.cpp:3201
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
print_matrix
void print_matrix(gsl_matrix *M, unsigned int r, unsigned int c, const string &message)
Definition: transform.cpp:1167
Generated by doxygen 1.9.2