Logo ROOT   6.16/01
Reference Guide
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
31namespace ROOT {
32namespace 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
48
49public:
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 */
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
317private:
318
323
324}; // Quaternion
325
326// ============ Class Quaternion ends here ============
327
328/**
329 Distance between two rotations
330 */
331template <class R>
332inline
333typename Quaternion::Scalar
334Distance ( const Quaternion& r1, const R & r2) {return gv_detail::dist(r1,r2);}
335
336/**
337 Multiplication of an axial rotation by an AxisAngle
338 */
339Quaternion operator* (RotationX const & r1, Quaternion const & r2);
340Quaternion operator* (RotationY const & r1, Quaternion const & r2);
341Quaternion 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
348std::ostream & operator<< (std::ostream & os, const Quaternion & q);
349
350
351} // namespace Math
352} // namespace ROOT
353
354#endif // ROOT_Math_GenVector_Quaternion
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::R::TRInterface & r
Definition: Object.C:4
#define R(a, b, c, d, e, f, g, h, i)
Definition: RSha256.hxx:110
#define e(i)
Definition: RSha256.hxx:103
float * q
Definition: THbookFile.cxx:87
AxisAngle class describing rotation represented with direction axis (3D Vector) and an angle of rotat...
Definition: AxisAngle.h:41
DefaultCoordinateSystemTag Default tag for identifying any coordinate system.
Class describing a generic displacement vector in 3 dimensions.
Scalar X() const
Cartesian X, converting if necessary from internal coordinate system.
Scalar Y() const
Cartesian Y, converting if necessary from internal coordinate system.
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...
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
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:48
Class describing a generic position vector (point) in 3 dimensions.
Rotation class with the (3D) rotation represented by a unit quaternion (u, i, j, k).
Definition: Quaternion.h:47
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
Scalar I() const
Definition: Quaternion.h:172
Quaternion(const OtherRotation &r)
Construct from another supported rotation type (see gv_detail::convert )
Definition: Quaternion.h:78
bool operator!=(const Quaternion &rhs) const
Definition: Quaternion.h:313
void SetComponents(Scalar u, Scalar i, Scalar j, Scalar k)
Set the components based on four Scalars.
Definition: Quaternion.h:155
Quaternion()
Default constructor (identity rotation)
Definition: Quaternion.h:58
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
DisplacementVector3D< Cartesian3D< double >, DefaultCoordinateSystemTag > XYZVector
Rotation operation on a cartesian vector.
Definition: Quaternion.h:181
void Rectify()
Re-adjust components to eliminate small deviations from |Q| = 1 orthonormality.
Definition: Quaternion.cxx:35
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
Scalar U() const
Access to the four quaternion components: U() is the coefficient of the identity Pauli matrix,...
Definition: Quaternion.h:171
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 & operator*=(const R &r)
Post-Multiply (on right) by another rotation : T = T*R.
Definition: Quaternion.h:288
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
XYZVector operator()(const XYZVector &v) const
Definition: Quaternion.h:182
AVector operator*(const AVector &v) const
Overload operator * for rotation on a vector.
Definition: Quaternion.h:246
Scalar K() const
Definition: Quaternion.h:174
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
void GetComponents(Scalar &u, Scalar &i, Scalar &j, Scalar &k) const
Get the components into four Scalars.
Definition: Quaternion.h:162
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
Scalar J() const
Definition: Quaternion.h:173
void Invert()
Invert a rotation in place.
Definition: Quaternion.h:254
Rotation class with the (3D) rotation represented by a 3x3 orthogonal matrix.
Definition: Rotation3D.h:65
Rotation class representing a 3D rotation about the X axis by the angle of rotation.
Definition: RotationX.h:43
Rotation class representing a 3D rotation about the Y axis by the angle of rotation.
Definition: RotationY.h:43
Rotation class with the (3D) rotation represented by angles describing first a rotation of an angle p...
Definition: RotationZYX.h:61
Rotation class representing a 3D rotation about the Z axis by the angle of rotation.
Definition: RotationZ.h:43
Namespace for new Math classes and functions.
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
void convert(R1 const &, R2 const)
Definition: 3DConversions.h:41
std::ostream & operator<<(std::ostream &os, const AxisAngle &a)
Stream Output and Input.
Definition: AxisAngle.cxx:91
Rotation3D::Scalar Scalar
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
auto * a
Definition: textangle.C:12