EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_transformations.h
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ transformations headers
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Iwan Kawrakow, 2005
25 #
26 # Contributors: Marc Chamberland
27 # Frederic Tessier
28 # Reid Townson
29 # Jan Weidner
30 # Alexandre Demelo
31 #
32 ###############################################################################
33 */
34 
35 
42 #ifndef EGS_TRANSFORMATIONS_
43 
44 #define EGS_TRANSFORMATIONS_
45 
46 #include "egs_vector.h"
47 #include "egs_libconfig.h"
48 #include "egs_math.h"
49 #include "egs_functions.h"
50 
51 #include <iostream>
52 #include <vector>
53 
54 using namespace std;
55 
56 class EGS_Input;
57 
70 
71 protected:
72 
74  EGS_Float rxx, rxy, rxz,
75  ryx, ryy, ryz,
76  rzx, rzy, rzz;
77 
78 public:
79 
82  rxx=1;
83  rxy=0;
84  rxz=0;
85  ryx=0;
86  ryy=1;
87  ryz=0;
88  rzx=0;
89  rzy=0;
90  rzz=1;
91  };
92 
93  //
94  // hopefully no one will use this constructor
95  // with a non-rotation matrix elements
96  //
98  EGS_RotationMatrix(EGS_Float xx, EGS_Float xy, EGS_Float xz,
99  EGS_Float yx, EGS_Float yy, EGS_Float yz,
100  EGS_Float zx, EGS_Float zy, EGS_Float zz) {
101  rxx=xx;
102  rxy=xy;
103  rxz=xz;
104  ryx=yx;
105  ryy=yy;
106  ryz=yz;
107  rzx=zx;
108  rzy=zy;
109  rzz=zz;
110  };
111 
114  rxx=m.rxx;
115  rxy=m.rxy;
116  rxz=m.rxz;
117  ryx=m.ryx;
118  ryy=m.ryy;
119  ryz=m.ryz;
120  rzx=m.rzx;
121  rzy=m.rzy;
122  rzz=m.rzz;
123  };
124 
132  bool isRotation() const {
133  EGS_Float d = det() - 1;
134  EGS_RotationMatrix t((*this)*inverse());
135  if (fabs(d) > epsilon || !t.isI()) {
136  return false;
137  }
138  return true;
139  };
140 
141  // creates a matrix which, when applied to the vector v,
142  // transforms it into a vector along the z-axis.
143  // why z-axis?
144  // well, it has to be one of the axis and in physics
145  // things usually happen along the z-axis!
154  EGS_Float sinz = v.x*v.x + v.y*v.y;
155  EGS_Float norm = sinz + v.z*v.z;
156  if (norm < epsilon) egsFatal("EGS_RotationMatrix::EGS_RotationMatrix: \n"
157  " no construction from a zero vector possible!\n");
158  norm = sqrt(norm);
159  if (sinz > epsilon) {
160  sinz = sqrt(sinz);
161  EGS_Float cphi = v.x/sinz;
162  register EGS_Float sphi = v.y/sinz;
163  EGS_Float cost = v.z/norm;
164  register EGS_Float sint = -sinz/norm;
165  *this = rotY(cost,sint)*rotZ(cphi,sphi);
166  }
167  else if (v.z < 0.) {
168  // v along negative z is degenerate: a rotation of pi around any axis in the
169  // xy plane satisfies the condition; we pick the x axis
170  *this = rotX(M_PI);
171  }
172  else { // v is along the z-axis => matrix is the unit transformation
173  rxx = 1;
174  rxy = 0;
175  rxz = 0;
176  ryx = 0;
177  ryy = 1;
178  ryz = 0;
179  rzx = 0;
180  rzy = 0;
181  rzz = 1;
182  }
183  };
184 
185  // R_x(alpha) * R_y(beta) * R_z(gamma)
189  EGS_RotationMatrix(EGS_Float alpha, EGS_Float beta, EGS_Float gamma) {
190  *this = rotX(alpha)*rotY(beta)*rotZ(gamma);
191  };
192 
193  // R_z(phi) * R_x(theta)
197  EGS_RotationMatrix(EGS_Float phi, EGS_Float theta) {
198  *this = rotZ(phi)*rotX(theta);
199  };
200 
203  rxx=m.rxx;
204  rxy=m.rxy;
205  rxz=m.rxz;
206  ryx=m.ryx;
207  ryy=m.ryy;
208  ryz=m.ryz;
209  rzx=m.rzx;
210  rzy=m.rzy;
211  rzz=m.rzz;
212  return *this;
213  };
214 
217  return
218  ((rxx == m.rxx) && (rxy == m.rxy) && (rxz == m.rxz) &&
219  (ryx == m.ryx) && (ryy == m.ryy) && (ryz == m.ryz) &&
220  (rzx == m.rxz) && (rzy == m.rzy) && (rzz == m.rzz)) ? true : false;
221  };
222 
225  bool isI() const {
226 
227  // compare with unit matrix components
228  if (fabs(rxy) < epsilon &&
229  fabs(rxz) < epsilon &&
230  fabs(ryx) < epsilon &&
231  fabs(ryz) < epsilon &&
232  fabs(rzx) < epsilon &&
233  fabs(rzy) < epsilon &&
234  fabs(rxx-1) < epsilon &&
235  fabs(ryy-1) < epsilon &&
236  fabs(rzz-1) < epsilon) {
237  return true;
238  }
239  return false;
240  };
241 
243  EGS_Vector operator*(const EGS_Vector &v) const {
244  return EGS_Vector(rxx*v.x + rxy*v.y + rxz*v.z,
245  ryx*v.x + ryy*v.y + ryz*v.z,
246  rzx*v.x + rzy*v.y + rzz*v.z);
247  };
248 
252  return EGS_RotationMatrix(
253  rxx*m.rxx+rxy*m.ryx+rxz*m.rzx, rxx*m.rxy+rxy*m.ryy+rxz*m.rzy,
254  rxx*m.rxz+rxy*m.ryz+rxz*m.rzz,
255  ryx*m.rxx+ryy*m.ryx+ryz*m.rzx, ryx*m.rxy+ryy*m.ryy+ryz*m.rzy,
256  ryx*m.rxz+ryy*m.ryz+ryz*m.rzz,
257  rzx*m.rxx+rzy*m.ryx+rzz*m.rzx, rzx*m.rxy+rzy*m.ryy+rzz*m.rzy,
258  rzx*m.rxz+rzy*m.ryz+rzz*m.rzz);
259  };
260 
261  // r *= m means r = r*m
266  return *this = operator * (m);
267  };
268 
269  // r.multiply(m) means r = m*r;
277  void multiply(const EGS_RotationMatrix &m) {
278  *this = m.operator * (*this);
279  };
280 
287  return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
288  };
289 
292  return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
293  };
294 
297  return *this = inverse();
298  };
299 
303  static EGS_RotationMatrix rotX(EGS_Float cphi,EGS_Float sphi) {
304  return EGS_RotationMatrix(
305  (EGS_Float)1,(EGS_Float)0,(EGS_Float)0,
306  (EGS_Float)0, cphi, sphi,
307  (EGS_Float)0, -sphi, cphi);
308  };
309 
313  static EGS_RotationMatrix rotY(EGS_Float cphi,EGS_Float sphi) {
314  return EGS_RotationMatrix(cphi, (EGS_Float)0, -sphi,
315  (EGS_Float)0, (EGS_Float)1, (EGS_Float)0,
316  sphi, (EGS_Float)0, cphi);
317  };
318 
322  static EGS_RotationMatrix rotZ(EGS_Float cphi,EGS_Float sphi) {
323  return EGS_RotationMatrix(cphi, sphi, (EGS_Float)0,
324  -sphi, cphi, (EGS_Float)0,
325  (EGS_Float)0,(EGS_Float)0, (EGS_Float)1);
326  };
327 
329  static EGS_RotationMatrix rotX(EGS_Float phi) {
330  return rotX(cos(phi),sin(phi));
331  };
332 
334  static EGS_RotationMatrix rotY(EGS_Float phi) {
335  return rotY(cos(phi),sin(phi));
336  };
337 
339  static EGS_RotationMatrix rotZ(EGS_Float phi) {
340  return rotZ(cos(phi),sin(phi));
341  };
342 
346  static EGS_RotationMatrix rotV(EGS_Float phi, const EGS_Vector &v) {
347  return rotV(cos(phi),sin(phi),v);
348  };
349 
351  static EGS_RotationMatrix rotV(EGS_Float cphi, EGS_Float sphi,
352  const EGS_Vector &v) {
353  EGS_RotationMatrix m(v);
354  return (m.inverse())*(rotZ(cphi,sphi))*(m);
355  };
356 
358  EGS_Float det() const {
359  return rxx*ryy*rzz + rxy*ryz*rzx + ryx*rzy*rxz -
360  rxz*ryy*rzx - rxy*ryx*rzz - rzy*ryz*rxx;
361  };
362 
366  friend EGS_Vector operator*(const EGS_Vector &v,
367  const EGS_RotationMatrix &m) {
368  return EGS_Vector(v.x*m.rxx+v.y*m.ryx+v.z*m.rzx,
369  v.x*m.rxy+v.y*m.ryy+v.z*m.rzy,
370  v.x*m.rxz+v.y*m.ryz+v.z*m.rzz);
371  };
372 
378  v = v*m;
379  return v;
380  };
381 
382  inline EGS_Float xx() const {
383  return rxx;
384  };
385  inline EGS_Float xy() const {
386  return rxy;
387  };
388  inline EGS_Float xz() const {
389  return rxz;
390  };
391  inline EGS_Float yx() const {
392  return ryx;
393  };
394  inline EGS_Float yy() const {
395  return ryy;
396  };
397  inline EGS_Float yz() const {
398  return ryz;
399  };
400  inline EGS_Float zx() const {
401  return rzx;
402  };
403  inline EGS_Float zy() const {
404  return rzy;
405  };
406  inline EGS_Float zz() const {
407  return rzz;
408  };
409 
410 };
411 
423 
424 protected:
425 
427  EGS_Vector t;
428  bool has_t, has_R;
429 
430 public:
431 
433  EGS_AffineTransform() : R(),t(),has_t(false),has_R(false) {};
434 
437  R(tr.R),t(tr.t),has_t(tr.has_t),has_R(tr.has_R) {};
438 
442  R(m),t(v) {
443  if (t.length2() > 0) {
444  has_t = true;
445  }
446  else {
447  has_t = false;
448  }
449  if (R.isI()) {
450  has_R = false;
451  }
452  else {
453  has_R = true;
454  }
455  };
456 
460  has_t = false;
461  if (R.isI()) {
462  has_R = false;
463  }
464  else {
465  has_R = true;
466  }
467  };
468 
471  EGS_AffineTransform(const EGS_Vector &v) : R(),t(v) {
472  has_R = false;
473  if (t.length2() > 0) {
474  has_t = true;
475  }
476  else {
477  has_t = false;
478  }
479  };
480 
491  return EGS_AffineTransform(R*tr.R,R*tr.t+t);
492  };
493 
498  return EGS_AffineTransform(R*m,t);
499  };
500 
507  return *this = operator * (tr);
508  };
509 
516  return *this = operator * (m);
517  };
518 
519  EGS_AffineTransform operator+(const EGS_Vector &v) const {
520  return EGS_AffineTransform(R,t+v);
521  };
522 
523  EGS_AffineTransform &operator+=(const EGS_Vector &v) {
524  t += v;
525  return *this;
526  };
527 
531  EGS_Vector operator*(const EGS_Vector &v) const {
532  return (R*v + t);
533  };
534 
537  friend EGS_Vector operator*(const EGS_Vector &v,
538  const EGS_AffineTransform &tr) {
539  return ((v-tr.t)*tr.R);
540  };
541 
545  const EGS_AffineTransform &tr) {
546  v = v*tr;
547  return v;
548  };
549 
551  void transform(EGS_Vector &v) const {
552  if (has_R) {
553  v = R*v;
554  }
555  if (has_t) {
556  v += t;
557  }
558  };
559 
561  void inverseTransform(EGS_Vector &v) const {
562  if (has_t) {
563  v -= t;
564  }
565  if (has_R) {
566  v *= R;
567  }
568  //v -= t; v *= R;
569  };
570 
573  EGS_Vector tmp;
574  tmp -= t*R;
575  return EGS_AffineTransform(R.inverse(),tmp);
576  };
577 
579  void rotate(EGS_Vector &v) const {
580  if (has_R) {
581  v = R*v;
582  }
583  };
585  void rotateInverse(EGS_Vector &v) const {
586  if (has_R) {
587  v *= R;
588  }
589  };
591  void translate(EGS_Vector &v) const {
592  v += t;
593  };
594 
597  const EGS_Vector &getTranslation() const {
598  return t;
599  };
603  return R;
604  };
605 
608  bool isI() const {
609  return (!has_R && !has_t);
610  };
611 
614  bool hasTranslation() const {
615  return has_t;
616  };
617 
620  bool hasRotation() const {
621  return has_R;
622  };
623 
658  static EGS_AffineTransform *getTransformation(EGS_Input *inp);
659 
660  static EGS_AffineTransform *getTransformation(vector<EGS_Float> trnsl, vector<EGS_Float> rot);
661 };
662 
663 #endif
A class providing affine transformations.
bool hasRotation() const
Returns true if the transformation involves a rotation, false otherwise.
bool hasTranslation() const
Returns true if the transformation involves a translation, false otherwise.
EGS_AffineTransform(const EGS_Vector &v)
Constructs an affine transformation object from the translation v, which has no rotation.
EGS_AffineTransform(const EGS_RotationMatrix &m, const EGS_Vector &v)
Constructs an affine transformation object from the rotation m and translation v.
EGS_AffineTransform operator*(const EGS_RotationMatrix &m) const
Returns the affine transformation , where and are the rotation and translation of the invoking obje...
void rotate(EGS_Vector &v) const
Applies the rotation to the vector v.
EGS_AffineTransform inverse() const
Returns the inverse affine transformation.
EGS_AffineTransform operator*(const EGS_AffineTransform &tr) const
Returns the multiplication of the invoking object with tr.
EGS_AffineTransform(const EGS_RotationMatrix &m)
Constructs an affine transformation object from the rotation m, which has no translation.
EGS_Vector operator*(const EGS_Vector &v) const
Applies the transformation to the vector v from the left and returns the result.
const EGS_RotationMatrix & getRotation() const
Returns the rotation matrix of the affine transformation object.
EGS_AffineTransform & operator*=(const EGS_RotationMatrix &m)
Multiplies the invoking object from the right with m. Returns a reference to the result.
EGS_AffineTransform(const EGS_AffineTransform &tr)
Copy constructor.
friend EGS_Vector operator*(const EGS_Vector &v, const EGS_AffineTransform &tr)
Applies the transformation tr to the invoking vector from the right and returns the result.
void rotateInverse(EGS_Vector &v) const
Applies the inverse rotation to the vector v.
EGS_AffineTransform & operator*=(const EGS_AffineTransform &tr)
Multiplies the invoking object from the right with tr. Returns a reference to the result.
void inverseTransform(EGS_Vector &v) const
Applies the inverse transformation to the vector v.
const EGS_Vector & getTranslation() const
Returns the translation vector of the affine transformation object.
friend EGS_Vector & operator*=(EGS_Vector &v, const EGS_AffineTransform &tr)
Applies the transformation tr to the invoking vector from the right, assignes the result to v and ret...
bool isI() const
Returns true if the object is a unity transformation, false otherwise.
void translate(EGS_Vector &v) const
Applies the translation to the vector v.
EGS_AffineTransform()
Constructs a unit affine transformation.
void transform(EGS_Vector &v) const
Transforms the vector v.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
A class for vector rotations.
static EGS_RotationMatrix rotV(EGS_Float cphi, EGS_Float sphi, const EGS_Vector &v)
static EGS_RotationMatrix rotZ(EGS_Float phi)
Returns a rotation around the z-axis by the angle phi.
EGS_RotationMatrix(EGS_Float phi, EGS_Float theta)
Constructs a rotation matrix from the angles theta and phi (polar and azimuthal) as .
EGS_RotationMatrix()
Default constructor, results in a unit matrix object.
EGS_RotationMatrix operator*(const EGS_RotationMatrix &m) const
Multiplies the invoking object with m from the right and returns the result.
EGS_RotationMatrix(EGS_Float alpha, EGS_Float beta, EGS_Float gamma)
Constructs a rotation matrix from rotation angles around the x-, y-, and z-axis as .
bool isRotation() const
Is this object a real rotation matrix?
static EGS_RotationMatrix rotY(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the y-axis by the angle with cphi, sphi = .
EGS_RotationMatrix T()
Returns the transposed matrix.
EGS_RotationMatrix(const EGS_Vector &v)
Create a rotation matrix from the vector v.
EGS_Vector operator*(const EGS_Vector &v) const
Returns the rotated a vector .
static EGS_RotationMatrix rotV(EGS_Float phi, const EGS_Vector &v)
Returns a rotation by the angle phi around the axis defined by the vector v.
EGS_RotationMatrix & operator*=(const EGS_RotationMatrix &m)
Multiplies the invoking object with m from the right and assigns the resulting matrix to the invoking...
friend EGS_Vector operator*(const EGS_Vector &v, const EGS_RotationMatrix &m)
Multiplies the invoking vector v from the right with the matrix m and returns the result.
EGS_Float det() const
Calculates and returns the determinant of the matrix.
static EGS_RotationMatrix rotX(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the x-axis by the angle with cphi, sphi = .
EGS_RotationMatrix & operator=(const EGS_RotationMatrix &m)
Assignment operator.
EGS_RotationMatrix(const EGS_RotationMatrix &m)
Copy constructor.
static EGS_RotationMatrix rotY(EGS_Float phi)
Returns a rotation around the y-axis by the angle phi.
friend EGS_Vector & operator*=(EGS_Vector &v, const EGS_RotationMatrix &m)
Multiplies the invoking vector v from the right with the matrix m and assigns the result to the invok...
EGS_RotationMatrix & invert()
Inverts the matrix and returns a reference to it.
static EGS_RotationMatrix rotZ(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the z-axis by the angle with cphi, sphi = .
bool operator==(const EGS_RotationMatrix &m)
Comparison operator.
EGS_RotationMatrix(EGS_Float xx, EGS_Float xy, EGS_Float xz, EGS_Float yx, EGS_Float yy, EGS_Float yz, EGS_Float zx, EGS_Float zy, EGS_Float zz)
Construct a rotation matrix object from 9 floating point numbers.
void multiply(const EGS_RotationMatrix &m)
Multiplies the invoking object with m from the left and returns the result.
bool isI() const
Returns true, if this object is approximately the unit matrix, false otherwise.
EGS_RotationMatrix inverse() const
Returns the inverse matrix.
static EGS_RotationMatrix rotX(EGS_Float phi)
Returns a rotation around the x-axis by the angle phi.
A class representing 3D vectors.
Definition: egs_vector.h:57
EGS_Float y
y-component
Definition: egs_vector.h:62
EGS_Float z
z-component
Definition: egs_vector.h:63
EGS_Float x
x-component
Definition: egs_vector.h:61
Global egspp functions header file.
Defines the EGS_EXPORT and EGS_LOCAL macros.
#define EGS_EXPORT
Export symbols from the egspp library.
Definition: egs_libconfig.h:90
Attempts to fix broken math header files.
EGS_Vector methods for the manipulation of 3D vectors in cartesian co-ordinates.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:62