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:
27 #
28 ###############################################################################
29 */
30 
31 
38 #ifndef EGS_TRANSFORMATIONS_
39 
40 #define EGS_TRANSFORMATIONS_
41 
42 #include "egs_vector.h"
43 #include "egs_libconfig.h"
44 #include "egs_math.h"
45 #include "egs_functions.h"
46 
47 #include <iostream>
48 
49 using namespace std;
50 
51 class EGS_Input;
52 
65 
66 protected:
67 
69  EGS_Float rxx, rxy, rxz,
70  ryx, ryy, ryz,
71  rzx, rzy, rzz;
72 
73 public:
74 
77  rxx=1;
78  rxy=0;
79  rxz=0;
80  ryx=0;
81  ryy=1;
82  ryz=0;
83  rzx=0;
84  rzy=0;
85  rzz=1;
86  };
87 
88  //
89  // hopefully no one will use this constructor
90  // with a non-rotation matrix elements
91  //
93  EGS_RotationMatrix(EGS_Float xx, EGS_Float xy, EGS_Float xz,
94  EGS_Float yx, EGS_Float yy, EGS_Float yz,
95  EGS_Float zx, EGS_Float zy, EGS_Float zz) {
96  rxx=xx;
97  rxy=xy;
98  rxz=xz;
99  ryx=yx;
100  ryy=yy;
101  ryz=yz;
102  rzx=zx;
103  rzy=zy;
104  rzz=zz;
105  };
106 
109  rxx=m.rxx;
110  rxy=m.rxy;
111  rxz=m.rxz;
112  ryx=m.ryx;
113  ryy=m.ryy;
114  ryz=m.ryz;
115  rzx=m.rzx;
116  rzy=m.rzy;
117  rzz=m.rzz;
118  };
119 
127  bool isRotation() const {
128  EGS_Float d = det() - 1;
129  EGS_RotationMatrix t((*this)*inverse());
130  if (fabs(d) > epsilon || !t.isI()) {
131  return false;
132  }
133  return true;
134  };
135 
136  // creates a matrix which, when applied to the vector v,
137  // transforms it into a vector along the z-axis.
138  // why z-axis?
139  // well, it has to be one of the axis and in physics
140  // things usually happen along the z-axis!
149  EGS_Float sinz = v.x*v.x + v.y*v.y;
150  EGS_Float norm = sinz + v.z*v.z;
151  if (norm < epsilon) egsFatal("EGS_RotationMatrix::EGS_RotationMatrix: \n"
152  " no construction from a zero vector possible!\n");
153  norm = sqrt(norm);
154  if (sinz > epsilon) {
155  sinz = sqrt(sinz);
156  EGS_Float cphi = v.x/sinz;
157  register EGS_Float sphi = v.y/sinz;
158  EGS_Float cost = v.z/norm;
159  register EGS_Float sint = -sinz/norm;
160  *this = rotY(cost,sint)*rotZ(cphi,sphi);
161  }
162  else if (v.z < 0.) {
163  // v along negative z is degenerate: a rotation of pi around any axis in the
164  // xy plane satisfies the condition; we pick the x axis
165  *this = rotX(M_PI);
166  }
167  else { // v is along the z-axis => matrix is the unit transformation
168  rxx = 1;
169  rxy = 0;
170  rxz = 0;
171  ryx = 0;
172  ryy = 1;
173  ryz = 0;
174  rzx = 0;
175  rzy = 0;
176  rzz = 1;
177  }
178  };
179 
180  // R_x(alpha) * R_y(beta) * R_z(gamma)
184  EGS_RotationMatrix(EGS_Float alpha, EGS_Float beta, EGS_Float gamma) {
185  *this = rotX(alpha)*rotY(beta)*rotZ(gamma);
186  };
187 
188  // R_z(phi) * R_x(theta)
192  EGS_RotationMatrix(EGS_Float phi, EGS_Float theta) {
193  *this = rotZ(phi)*rotX(theta);
194  };
195 
198  rxx=m.rxx;
199  rxy=m.rxy;
200  rxz=m.rxz;
201  ryx=m.ryx;
202  ryy=m.ryy;
203  ryz=m.ryz;
204  rzx=m.rzx;
205  rzy=m.rzy;
206  rzz=m.rzz;
207  return *this;
208  };
209 
212  return
213  ((rxx == m.rxx) && (rxy == m.rxy) && (rxz == m.rxz) &&
214  (ryx == m.ryx) && (ryy == m.ryy) && (ryz == m.ryz) &&
215  (rzx == m.rxz) && (rzy == m.rzy) && (rzz == m.rzz)) ? true : false;
216  };
217 
220  bool isI() const {
221 
222  // compare with unit matrix components
223  if (fabs(rxy) < epsilon &&
224  fabs(rxz) < epsilon &&
225  fabs(ryx) < epsilon &&
226  fabs(ryz) < epsilon &&
227  fabs(rzx) < epsilon &&
228  fabs(rzy) < epsilon &&
229  fabs(rxx-1) < epsilon &&
230  fabs(ryy-1) < epsilon &&
231  fabs(rzz-1) < epsilon) {
232  return true;
233  }
234  return false;
235  };
236 
238  EGS_Vector operator*(const EGS_Vector &v) const {
239  return EGS_Vector(rxx*v.x + rxy*v.y + rxz*v.z,
240  ryx*v.x + ryy*v.y + ryz*v.z,
241  rzx*v.x + rzy*v.y + rzz*v.z);
242  };
243 
247  return EGS_RotationMatrix(
248  rxx*m.rxx+rxy*m.ryx+rxz*m.rzx, rxx*m.rxy+rxy*m.ryy+rxz*m.rzy,
249  rxx*m.rxz+rxy*m.ryz+rxz*m.rzz,
250  ryx*m.rxx+ryy*m.ryx+ryz*m.rzx, ryx*m.rxy+ryy*m.ryy+ryz*m.rzy,
251  ryx*m.rxz+ryy*m.ryz+ryz*m.rzz,
252  rzx*m.rxx+rzy*m.ryx+rzz*m.rzx, rzx*m.rxy+rzy*m.ryy+rzz*m.rzy,
253  rzx*m.rxz+rzy*m.ryz+rzz*m.rzz);
254  };
255 
256  // r *= m means r = r*m
261  return *this = operator * (m);
262  };
263 
264  // r.multiply(m) means r = m*r;
272  void multiply(const EGS_RotationMatrix &m) {
273  *this = m.operator * (*this);
274  };
275 
282  return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
283  };
284 
287  return EGS_RotationMatrix(rxx,ryx,rzx,rxy,ryy,rzy,rxz,ryz,rzz);
288  };
289 
292  return *this = inverse();
293  };
294 
298  static EGS_RotationMatrix rotX(EGS_Float cphi,EGS_Float sphi) {
299  return EGS_RotationMatrix(
300  (EGS_Float)1,(EGS_Float)0,(EGS_Float)0,
301  (EGS_Float)0, cphi, sphi,
302  (EGS_Float)0, -sphi, cphi);
303  };
304 
308  static EGS_RotationMatrix rotY(EGS_Float cphi,EGS_Float sphi) {
309  return EGS_RotationMatrix(cphi, (EGS_Float)0, -sphi,
310  (EGS_Float)0, (EGS_Float)1, (EGS_Float)0,
311  sphi, (EGS_Float)0, cphi);
312  };
313 
317  static EGS_RotationMatrix rotZ(EGS_Float cphi,EGS_Float sphi) {
318  return EGS_RotationMatrix(cphi, sphi, (EGS_Float)0,
319  -sphi, cphi, (EGS_Float)0,
320  (EGS_Float)0,(EGS_Float)0, (EGS_Float)1);
321  };
322 
324  static EGS_RotationMatrix rotX(EGS_Float phi) {
325  return rotX(cos(phi),sin(phi));
326  };
327 
329  static EGS_RotationMatrix rotY(EGS_Float phi) {
330  return rotY(cos(phi),sin(phi));
331  };
332 
334  static EGS_RotationMatrix rotZ(EGS_Float phi) {
335  return rotZ(cos(phi),sin(phi));
336  };
337 
341  static EGS_RotationMatrix rotV(EGS_Float phi, const EGS_Vector &v) {
342  return rotV(cos(phi),sin(phi),v);
343  };
344 
346  static EGS_RotationMatrix rotV(EGS_Float cphi, EGS_Float sphi,
347  const EGS_Vector &v) {
348  EGS_RotationMatrix m(v);
349  return (m.inverse())*(rotZ(cphi,sphi))*(m);
350  };
351 
353  EGS_Float det() const {
354  return rxx*ryy*rzz + rxy*ryz*rzx + ryx*rzy*rxz -
355  rxz*ryy*rzx - rxy*ryx*rzz - rzy*ryz*rxx;
356  };
357 
361  friend EGS_Vector operator*(const EGS_Vector &v,
362  const EGS_RotationMatrix &m) {
363  return EGS_Vector(v.x*m.rxx+v.y*m.ryx+v.z*m.rzx,
364  v.x*m.rxy+v.y*m.ryy+v.z*m.rzy,
365  v.x*m.rxz+v.y*m.ryz+v.z*m.rzz);
366  };
367 
373  v = v*m;
374  return v;
375  };
376 
377  inline EGS_Float xx() const {
378  return rxx;
379  };
380  inline EGS_Float xy() const {
381  return rxy;
382  };
383  inline EGS_Float xz() const {
384  return rxz;
385  };
386  inline EGS_Float yx() const {
387  return ryx;
388  };
389  inline EGS_Float yy() const {
390  return ryy;
391  };
392  inline EGS_Float yz() const {
393  return ryz;
394  };
395  inline EGS_Float zx() const {
396  return rzx;
397  };
398  inline EGS_Float zy() const {
399  return rzy;
400  };
401  inline EGS_Float zz() const {
402  return rzz;
403  };
404 
405 };
406 
418 
419 protected:
420 
422  EGS_Vector t;
423  bool has_t, has_R;
424 
425 public:
426 
428  EGS_AffineTransform() : R(),t(),has_t(false),has_R(false) {};
429 
432  R(tr.R),t(tr.t),has_t(tr.has_t),has_R(tr.has_R) {};
433 
437  R(m),t(v) {
438  if (t.length2() > 0) {
439  has_t = true;
440  }
441  else {
442  has_t = false;
443  }
444  if (R.isI()) {
445  has_R = false;
446  }
447  else {
448  has_R = true;
449  }
450  };
451 
455  has_t = false;
456  if (R.isI()) {
457  has_R = false;
458  }
459  else {
460  has_R = true;
461  }
462  };
463 
466  EGS_AffineTransform(const EGS_Vector &v) : R(),t(v) {
467  has_R = false;
468  if (t.length2() > 0) {
469  has_t = true;
470  }
471  else {
472  has_t = false;
473  }
474  };
475 
486  return EGS_AffineTransform(R*tr.R,R*tr.t+t);
487  };
488 
493  return EGS_AffineTransform(R*m,t);
494  };
495 
502  return *this = operator * (tr);
503  };
504 
511  return *this = operator * (m);
512  };
513 
514  EGS_AffineTransform operator+(const EGS_Vector &v) const {
515  return EGS_AffineTransform(R,t+v);
516  };
517 
518  EGS_AffineTransform &operator+=(const EGS_Vector &v) {
519  t += v;
520  return *this;
521  };
522 
526  EGS_Vector operator*(const EGS_Vector &v) const {
527  return (R*v + t);
528  };
529 
532  friend EGS_Vector operator*(const EGS_Vector &v,
533  const EGS_AffineTransform &tr) {
534  return ((v-tr.t)*tr.R);
535  };
536 
540  const EGS_AffineTransform &tr) {
541  v = v*tr;
542  return v;
543  };
544 
546  void transform(EGS_Vector &v) const {
547  if (has_R) {
548  v = R*v;
549  }
550  if (has_t) {
551  v += t;
552  }
553  };
554 
556  void inverseTransform(EGS_Vector &v) const {
557  if (has_t) {
558  v -= t;
559  }
560  if (has_R) {
561  v *= R;
562  }
563  //v -= t; v *= R;
564  };
565 
568  EGS_Vector tmp;
569  tmp -= t*R;
570  return EGS_AffineTransform(R.inverse(),tmp);
571  };
572 
574  void rotate(EGS_Vector &v) const {
575  if (has_R) {
576  v = R*v;
577  }
578  };
580  void rotateInverse(EGS_Vector &v) const {
581  if (has_R) {
582  v *= R;
583  }
584  };
586  void translate(EGS_Vector &v) const {
587  v += t;
588  };
589 
592  const EGS_Vector &getTranslation() const {
593  return t;
594  };
598  return R;
599  };
600 
603  bool isI() const {
604  return (!has_R && !has_t);
605  };
606 
609  bool hasTranslation() const {
610  return has_t;
611  };
612 
615  bool hasRotation() const {
616  return has_R;
617  };
618 
653  static EGS_AffineTransform *getTransformation(EGS_Input *inp);
654 
655 };
656 
657 #endif
#define EGS_EXPORT
Export symbols from the egspp library.
Definition: egs_libconfig.h:90
void translate(EGS_Vector &v) const
Applies the translation to 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*=(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 operator*(const EGS_RotationMatrix &m) const
Multiplies the invoking object with m from the right and returns the result.
static EGS_RotationMatrix rotY(EGS_Float phi)
Returns a rotation around the y-axis by the angle phi.
static EGS_RotationMatrix rotV(EGS_Float cphi, EGS_Float sphi, const EGS_Vector &v)
static EGS_RotationMatrix rotX(EGS_Float phi)
Returns a rotation around the x-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 ...
A class providing affine transformations.
EGS_RotationMatrix()
Default constructor, results in a unit matrix object.
A class for vector rotations.
EGS_Vector methods for the manipulation of 3D vectors in cartesian co-ordinates.
EGS_Float y
y-component
Definition: egs_vector.h:61
EGS_RotationMatrix & operator=(const EGS_RotationMatrix &m)
Assignment operator.
bool hasTranslation() const
Returns true if the transformation involves a translation, false otherwise.
A class representing 3D vectors.
Definition: egs_vector.h:56
EGS_AffineTransform(const EGS_RotationMatrix &m)
Constructs an affine transformation object from the rotation m, which has no translation.
Global egspp functions header file.
EGS_Vector operator*(const EGS_Vector &v) const
Returns the rotated a vector .
EGS_AffineTransform inverse() const
Returns the inverse affine transformation.
static EGS_RotationMatrix rotY(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the y-axis by the angle with cphi, sphi = .
bool operator==(const EGS_RotationMatrix &m)
Comparison operator.
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)
Multiplies the invoking object from the right with m. Returns a reference to 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 ...
EGS_RotationMatrix T()
Returns the transposed matrix.
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...
EGS_AffineTransform()
Constructs a unit affine transformation.
EGS_RotationMatrix & invert()
Inverts the matrix and returns a reference to it.
bool isI() const
Returns true, if this object is approximately the unit matrix, false otherwise.
EGS_Float z
z-component
Definition: egs_vector.h:62
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
static EGS_RotationMatrix rotZ(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the z-axis by the angle with cphi, sphi = .
bool hasRotation() const
Returns true if the transformation involves a rotation, false otherwise.
const EGS_Vector & getTranslation() const
Returns the translation vector of the affine transformation object.
EGS_Vector operator*(const EGS_Vector &v) const
Applies the transformation to the vector v from the left and returns the result.
EGS_RotationMatrix inverse() const
Returns the inverse matrix.
void rotate(EGS_Vector &v) const
Applies the rotation to the vector v.
bool isI() const
Returns true if the object is a unity transformation, false otherwise.
bool isRotation() const
Is this object a real rotation matrix?
EGS_AffineTransform & operator*=(const EGS_AffineTransform &tr)
Multiplies the invoking object from the right with tr. Returns a reference to the result...
EGS_RotationMatrix(const EGS_RotationMatrix &m)
Copy constructor.
EGS_AffineTransform operator*(const EGS_AffineTransform &tr) const
Returns the multiplication of the invoking object with tr.
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.
Attempts to fix broken math header files.
void rotateInverse(EGS_Vector &v) const
Applies the inverse rotation to the vector v.
EGS_Float det() const
Calculates and returns the determinant of the matrix.
void transform(EGS_Vector &v) const
Transforms the vector v.
EGS_Float x
x-component
Definition: egs_vector.h:60
void multiply(const EGS_RotationMatrix &m)
Multiplies the invoking object with m from the left and returns the result.
void inverseTransform(EGS_Vector &v) const
Applies the inverse transformation to the vector v.
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.
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:61
EGS_AffineTransform(const EGS_Vector &v)
Constructs an affine transformation object from the translation v, which has no rotation.
A class for storing information in a tree-like structure of key-value pairs. This class is used throu...
Definition: egs_input.h:182
EGS_AffineTransform operator*(const EGS_RotationMatrix &m) const
Returns the affine transformation , where and are the rotation and translation of the invoking obje...
EGS_AffineTransform(const EGS_AffineTransform &tr)
Copy constructor.
Defines the EGS_EXPORT and EGS_LOCAL macros.
EGS_RotationMatrix(const EGS_Vector &v)
Create a rotation matrix from the vector v.
const EGS_RotationMatrix & getRotation() const
Returns the rotation matrix of the affine transformation object.
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...
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...
static EGS_RotationMatrix rotZ(EGS_Float phi)
Returns a rotation around the z-axis by the angle phi.
static EGS_RotationMatrix rotX(EGS_Float cphi, EGS_Float sphi)
Returns a rotation around the x-axis by the angle with cphi, sphi = .