ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Quaternion.h
Go to the documentation of this file.
1 // @(#)root/mathcore:$Id$
2 // Authors: W. Brown, M. Fischler, L. Moneta 2005
3 
4  /**********************************************************************
5  * *
6  * Copyright (c) 2005 , LCG ROOT FNAL MathLib Team *
7  * *
8  * *
9  **********************************************************************/
10 
11 // Header file for rotation in 3 dimensions, represented by a quaternion
12 // Created by: Mark Fischler Thurs June 9 2005
13 //
14 // Last update: $Id$
15 //
16 #ifndef ROOT_Math_GenVector_Quaternion
17 #define ROOT_Math_GenVector_Quaternion 1
18 
19 
26 
27 #include <algorithm>
28 #include <cassert>
29 
30 
31 namespace ROOT {
32 namespace Math {
33 
34 
35 //__________________________________________________________________________________________
36  /**
37  Rotation class with the (3D) rotation represented by
38  a unit quaternion (u, i, j, k).
39  This is the optimal representation for multiplication of multiple
40  rotations, and for computation of group-manifold-invariant distance
41  between two rotations.
42  See also ROOT::Math::AxisAngle, ROOT::Math::EulerAngles, and ROOT::Math::Rotation3D.
43 
44  @ingroup GenVector
45  */
46 
47 class Quaternion {
48 
49 public:
50 
51  typedef double Scalar;
52 
53  // ========== Constructors and Assignment =====================
54 
55  /**
56  Default constructor (identity rotation)
57  */
59  : fU(1.0)
60  , fI(0.0)
61  , fJ(0.0)
62  , fK(0.0)
63  { }
64 
65  /**
66  Construct given a pair of pointers or iterators defining the
67  beginning and end of an array of four Scalars
68  */
69  template<class IT>
70  Quaternion(IT begin, IT end) { SetComponents(begin,end); }
71 
72  // ======== Construction From other Rotation Forms ==================
73 
74  /**
75  Construct from another supported rotation type (see gv_detail::convert )
76  */
77  template <class OtherRotation>
78  explicit Quaternion(const OtherRotation & r) {gv_detail::convert(r,*this);}
79 
80 
81  /**
82  Construct from four Scalars representing the coefficients of u, i, j, k
83  */
85  fU(u), fI(i), fJ(j), fK(k) { }
86 
87  // The compiler-generated copy ctor, copy assignment, and dtor are OK.
88 
89  /**
90  Re-adjust components to eliminate small deviations from |Q| = 1
91  orthonormality.
92  */
93  void Rectify();
94 
95  /**
96  Assign from another supported rotation type (see gv_detail::convert )
97  */
98  template <class OtherRotation>
99  Quaternion & operator=( OtherRotation const & r ) {
100  gv_detail::convert(r,*this);
101  return *this;
102  }
103 
104  // ======== Components ==============
105 
106  /**
107  Set the four components given an iterator to the start of
108  the desired data, and another to the end (4 past start).
109  */
110  template<class IT>
111 #ifndef NDEBUG
112  void SetComponents(IT begin, IT end) {
113 #else
114  void SetComponents(IT begin, IT ) {
115 #endif
116  fU = *begin++;
117  fI = *begin++;
118  fJ = *begin++;
119  fK = *begin++;
120  assert (end==begin);
121  }
122 
123  /**
124  Get the components into data specified by an iterator begin
125  and another to the end of the desired data (4 past start).
126  */
127  template<class IT>
128 #ifndef NDEBUG
129  void GetComponents(IT begin, IT end) const {
130 #else
131  void GetComponents(IT begin, IT ) const {
132 #endif
133  *begin++ = fU;
134  *begin++ = fI;
135  *begin++ = fJ;
136  *begin++ = fK;
137  assert (end==begin);
138  }
139 
140  /**
141  Get the components into data specified by an iterator begin
142  */
143  template<class IT>
144  void GetComponents(IT begin ) const {
145  *begin++ = fU;
146  *begin++ = fI;
147  *begin++ = fJ;
148  *begin = fK;
149  }
150 
151  /**
152  Set the components based on four Scalars. The sum of the squares of
153  these Scalars should be 1; no checking is done.
154  */
156  fU=u; fI=i; fJ=j; fK=k;
157  }
158 
159  /**
160  Get the components into four Scalars.
161  */
162  void GetComponents(Scalar & u, Scalar & i, Scalar & j, Scalar & k) const {
163  u=fU; i=fI; j=fJ; k=fK;
164  }
165 
166  /**
167  Access to the four quaternion components:
168  U() is the coefficient of the identity Pauli matrix,
169  I(), J() and K() are the coefficients of sigma_x, sigma_y, sigma_z
170  */
171  Scalar U() const { return fU; }
172  Scalar I() const { return fI; }
173  Scalar J() const { return fJ; }
174  Scalar K() const { return fK; }
175 
176  // =========== operations ==============
177 
178  /**
179  Rotation operation on a cartesian vector
180  */
182  XYZVector operator() (const XYZVector & v) const {
183 
184  const Scalar alpha = fU*fU - fI*fI - fJ*fJ - fK*fK;
185  const Scalar twoQv = 2*(fI*v.X() + fJ*v.Y() + fK*v.Z());
186  const Scalar twoU = 2 * fU;
187  return XYZVector ( alpha * v.X() + twoU * (fJ*v.Z() - fK*v.Y()) + twoQv * fI ,
188  alpha * v.Y() + twoU * (fK*v.X() - fI*v.Z()) + twoQv * fJ ,
189  alpha * v.Z() + twoU * (fI*v.Y() - fJ*v.X()) + twoQv * fK );
190  }
191 
192  /**
193  Rotation operation on a displacement vector in any coordinate system
194  */
195  template <class CoordSystem,class Tag>
198  DisplacementVector3D< Cartesian3D<double> > xyz(v.X(), v.Y(), v.Z());
201  vNew.SetXYZ( rxyz.X(), rxyz.Y(), rxyz.Z() );
202  return vNew;
203  }
204 
205  /**
206  Rotation operation on a position vector in any coordinate system
207  */
208  template <class CoordSystem, class Tag>
213  return PositionVector3D<CoordSystem,Tag> ( rxyz );
214  }
215 
216  /**
217  Rotation operation on a Lorentz vector in any 4D coordinate system
218  */
219  template <class CoordSystem>
223  xyz = operator()(xyz);
224  LorentzVector< PxPyPzE4D<double> > xyzt (xyz.X(), xyz.Y(), xyz.Z(), v.E());
225  return LorentzVector<CoordSystem> ( xyzt );
226  }
227 
228  /**
229  Rotation operation on an arbitrary vector v.
230  Preconditions: v must implement methods x(), y(), and z()
231  and the arbitrary vector type must have a constructor taking (x,y,z)
232  */
233  template <class ForeignVector>
234  ForeignVector
235  operator() (const ForeignVector & v) const {
238  return ForeignVector ( rxyz.X(), rxyz.Y(), rxyz.Z() );
239  }
240 
241  /**
242  Overload operator * for rotation on a vector
243  */
244  template <class AVector>
245  inline
246  AVector operator* (const AVector & v) const
247  {
248  return operator()(v);
249  }
250 
251  /**
252  Invert a rotation in place
253  */
254  void Invert() { fI = -fI; fJ = -fJ; fK = -fK; }
255 
256  /**
257  Return inverse of a rotation
258  */
259  Quaternion Inverse() const { return Quaternion(fU, -fI, -fJ, -fK); }
260 
261  // ========= Multi-Rotation Operations ===============
262 
263  /**
264  Multiply (combine) two rotations
265  */
266  /**
267  Multiply (combine) two rotations
268  */
270  return Quaternion ( fU*q.fU - fI*q.fI - fJ*q.fJ - fK*q.fK ,
271  fU*q.fI + fI*q.fU + fJ*q.fK - fK*q.fJ ,
272  fU*q.fJ - fI*q.fK + fJ*q.fU + fK*q.fI ,
273  fU*q.fK + fI*q.fJ - fJ*q.fI + fK*q.fU );
274  }
275 
276  Quaternion operator * (const Rotation3D & r) const;
277  Quaternion operator * (const AxisAngle & a) const;
278  Quaternion operator * (const EulerAngles & e) const;
279  Quaternion operator * (const RotationZYX & r) const;
280  Quaternion operator * (const RotationX & rx) const;
281  Quaternion operator * (const RotationY & ry) const;
282  Quaternion operator * (const RotationZ & rz) const;
283 
284  /**
285  Post-Multiply (on right) by another rotation : T = T*R
286  */
287  template <class R>
288  Quaternion & operator *= (const R & r) { return *this = (*this)*r; }
289 
290 
291  /**
292  Distance between two rotations in Quaternion form
293  Note: The rotation group is isomorphic to a 3-sphere
294  with diametrically opposite points identified.
295  The (rotation group-invariant) is the smaller
296  of the two possible angles between the images of
297  the two totations on that sphere. Thus the distance
298  is never greater than pi/2.
299  */
300 
301  Scalar Distance(const Quaternion & q) const ;
302 
303  /**
304  Equality/inequality operators
305  */
306  bool operator == (const Quaternion & rhs) const {
307  if( fU != rhs.fU ) return false;
308  if( fI != rhs.fI ) return false;
309  if( fJ != rhs.fJ ) return false;
310  if( fK != rhs.fK ) return false;
311  return true;
312  }
313  bool operator != (const Quaternion & rhs) const {
314  return ! operator==(rhs);
315  }
316 
317 private:
318 
323 
324 }; // Quaternion
325 
326 // ============ Class Quaternion ends here ============
327 
328 /**
329  Distance between two rotations
330  */
331 template <class R>
332 inline
333 typename Quaternion::Scalar
334 Distance ( const Quaternion& r1, const R & r2) {return gv_detail::dist(r1,r2);}
335 
336 /**
337  Multiplication of an axial rotation by an AxisAngle
338  */
339 Quaternion operator* (RotationX const & r1, Quaternion const & r2);
340 Quaternion operator* (RotationY const & r1, Quaternion const & r2);
341 Quaternion operator* (RotationZ const & r1, Quaternion const & r2);
342 
343 /**
344  Stream Output and Input
345  */
346  // TODO - I/O should be put in the manipulator form
347 
348 std::ostream & operator<< (std::ostream & os, const Quaternion & q);
349 
350 
351 } // namespace Math
352 } // namespace ROOT
353 
354 #endif // ROOT_Math_GenVector_Quaternion
void GetComponents(IT begin) const
Get the components into data specified by an iterator begin.
Definition: Quaternion.h:144
Quaternion Inverse() const
Return inverse of a rotation.
Definition: Quaternion.h:259
Scalar E() const
return 4-th component (time, or energy for a 4-momentum vector)
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:54
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
DisplacementVector3D< CoordSystem, Tag > & SetXYZ(Scalar a, Scalar b, Scalar c)
set the values of the vector from the cartesian components (x,y,z) (if the vector is held in polar or...
::ROOT::Math::DisplacementVector3D< Cartesian3D< Scalar > > Vect() const
get the spatial components of the Vector in a DisplacementVector based on Cartesian Coordinates ...
#define assert(cond)
Definition: unittest.h:542
Rotation class representing a 3D rotation about the Z axis by the angle of rotation.
Definition: RotationZ.h:43
XYZVector operator()(const XYZVector &v) const
Definition: Quaternion.h:182
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
AVector operator*(const AVector &v) const
Overload operator * for rotation on a vector.
Definition: Quaternion.h:246
Class describing a generic position vector (point) in 3 dimensions.
Rotation class with the (3D) rotation represented by angles describing first a rotation of an angle p...
Definition: RotationZYX.h:71
TArc * a
Definition: textangle.C:12
Quaternion & operator=(OtherRotation const &r)
Assign from another supported rotation type (see gv_detail::convert )
Definition: Quaternion.h:99
bool operator==(const Quaternion &rhs) const
Equality/inequality operators.
Definition: Quaternion.h:306
Rotation class with the (3D) rotation represented by a unit quaternion (u, i, j, k).
Definition: Quaternion.h:47
Scalar J() const
Definition: Quaternion.h:173
Scalar K() const
Definition: Quaternion.h:174
std::ostream & operator<<(std::ostream &os, const AxisAngle &a)
Stream Output and Input.
Definition: AxisAngle.cxx:91
AxisAngle class describing rotation represented with direction axis (3D Vector) and an angle of rotat...
Definition: AxisAngle.h:41
Rotation class representing a 3D rotation about the Y axis by the angle of rotation.
Definition: RotationY.h:43
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
void GetComponents(Scalar &u, Scalar &i, Scalar &j, Scalar &k) const
Get the components into four Scalars.
Definition: Quaternion.h:162
void GetComponents(IT begin, IT end) const
Get the components into data specified by an iterator begin and another to the end of the desired dat...
Definition: Quaternion.h:129
Quaternion(const OtherRotation &r)
Construct from another supported rotation type (see gv_detail::convert )
Definition: Quaternion.h:78
Scalar I() const
Definition: Quaternion.h:172
void SetComponents(Scalar u, Scalar i, Scalar j, Scalar k)
Set the components based on four Scalars.
Definition: Quaternion.h:155
Class describing a generic displacement vector in 3 dimensions.
ROOT::R::TRInterface & r
Definition: Object.C:4
Quaternion(Scalar u, Scalar i, Scalar j, Scalar k)
Construct from four Scalars representing the coefficients of u, i, j, k.
Definition: Quaternion.h:84
SVector< double, 2 > v
Definition: Dict.h:5
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
Rotation class representing a 3D rotation about the X axis by the angle of rotation.
Definition: RotationX.h:43
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
Scalar U() const
Access to the four quaternion components: U() is the coefficient of the identity Pauli matrix...
Definition: Quaternion.h:171
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
Rotation operation on a cartesian vector.
Definition: Quaternion.h:181
void convert(R1 const &, R2 const)
Definition: 3DConversions.h:41
Scalar Z() const
Cartesian Z, converting if necessary from internal coordinate system.
EulerAngles class describing rotation as three angles (Euler Angles).
Definition: EulerAngles.h:43
Scalar Distance(const Quaternion &q) const
Distance between two rotations in Quaternion form Note: The rotation group is isomorphic to a 3-spher...
Definition: Quaternion.cxx:92
bool operator!=(const Quaternion &rhs) const
Definition: Quaternion.h:313
void Rectify()
Re-adjust components to eliminate small deviations from |Q| = 1 orthonormality.
Definition: Quaternion.cxx:35
Quaternion & operator*=(const R &r)
Post-Multiply (on right) by another rotation : T = T*R.
Definition: Quaternion.h:288
void Invert()
Invert a rotation in place.
Definition: Quaternion.h:254
DefaultCoordinateSystemTag Default tag for identifying any coordinate system.
float * q
Definition: THbookFile.cxx:87
Plane3D::Scalar Scalar
Definition: Plane3D.cxx:29
void SetComponents(IT begin, IT end)
Set the four components given an iterator to the start of the desired data, and another to the end (4...
Definition: Quaternion.h:112
Quaternion()
Default constructor (identity rotation)
Definition: Quaternion.h:58
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
Quaternion(IT begin, IT end)
Construct given a pair of pointers or iterators defining the beginning and end of an array of four Sc...
Definition: Quaternion.h:70