HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ImathQuat.h
Go to the documentation of this file.
1 //
2 // SPDX-License-Identifier: BSD-3-Clause
3 // Copyright Contributors to the OpenEXR Project.
4 //
5 
6 //
7 // A quaternion
8 //
9 // "Quaternions came from Hamilton ... and have been an unmixed
10 // evil to those who have touched them in any way. Vector is a
11 // useless survival ... and has never been of the slightest use
12 // to any creature."
13 //
14 // - Lord Kelvin
15 //
16 
17 #ifndef INCLUDED_IMATHQUAT_H
18 #define INCLUDED_IMATHQUAT_H
19 
20 #include "ImathExport.h"
21 #include "ImathNamespace.h"
22 
23 #include "ImathMatrix.h"
24 
25 #include <iostream>
26 
27 IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
28 
29 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
30 // Disable MS VC++ warnings about conversion from double to float
31 # pragma warning(push)
32 # pragma warning(disable : 4244)
33 #endif
34 
35 ///
36 /// The Quat class implements the quaternion numerical type -- you
37 /// will probably want to use this class to represent orientations
38 /// in R3 and to convert between various euler angle reps. You
39 /// should probably use Imath::Euler<> for that.
40 ///
41 
42 template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Quat
43 {
44  public:
45 
46  /// @{
47  /// @name Direct access to elements
48 
49  /// The real part
50  T r;
51 
52  /// The imaginary vector
54 
55  /// @}
56 
57  /// Element access: q[0] is the real part, (q[1],q[2],q[3]) is the
58  /// imaginary part.
59  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T& operator[] (int index) IMATH_NOEXCEPT; // as 4D vector
60 
61  /// Element access: q[0] is the real part, (q[1],q[2],q[3]) is the
62  /// imaginary part.
63  IMATH_HOSTDEVICE constexpr T operator[] (int index) const IMATH_NOEXCEPT;
64 
65  /// @{
66  /// @name Constructors
67 
68  /// Default constructor is the identity quat
70 
71  /// Copy constructor
72  IMATH_HOSTDEVICE constexpr Quat (const Quat& q) IMATH_NOEXCEPT;
73 
74  /// Construct from a quaternion of a another base type
75  template <class S> IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat (const Quat<S>& q) IMATH_NOEXCEPT;
76 
77  /// Initialize with real part `s` and imaginary vector 1(i,j,k)`
78  IMATH_HOSTDEVICE constexpr Quat (T s, T i, T j, T k) IMATH_NOEXCEPT;
79 
80  /// Initialize with real part `s` and imaginary vector `d`
81  IMATH_HOSTDEVICE constexpr Quat (T s, Vec3<T> d) IMATH_NOEXCEPT;
82 
83  /// The identity quaternion
84  IMATH_HOSTDEVICE constexpr static Quat<T> identity() IMATH_NOEXCEPT;
85 
86  /// Assignment
87  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator= (const Quat<T>& q) IMATH_NOEXCEPT;
88 
89  /// Destructor
90  IMATH_HOSTDEVICE ~Quat() IMATH_NOEXCEPT = default;
91 
92  /// @}
93 
94  /// @{
95  /// @name Basic Algebra
96  ///
97  /// Note that the operator return values are *NOT* normalized
98  //
99 
100  /// Quaternion multiplication
101  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator*= (const Quat<T>& q) IMATH_NOEXCEPT;
102 
103  /// Scalar multiplication: multiply both real and imaginary parts
104  /// by the given scalar.
105  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator*= (T t) IMATH_NOEXCEPT;
106 
107  /// Quaterion division, using the inverse()
108  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator/= (const Quat<T>& q) IMATH_NOEXCEPT;
109 
110  /// Scalar division: multiply both real and imaginary parts
111  /// by the given scalar.
112  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator/= (T t) IMATH_NOEXCEPT;
113 
114  /// Quaternion addition
115  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator+= (const Quat<T>& q) IMATH_NOEXCEPT;
116 
117  /// Quaternion subtraction
118  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator-= (const Quat<T>& q) IMATH_NOEXCEPT;
119 
120  /// Equality
121  template <class S> IMATH_HOSTDEVICE constexpr bool operator== (const Quat<S>& q) const IMATH_NOEXCEPT;
122 
123  /// Inequality
124  template <class S> IMATH_HOSTDEVICE constexpr bool operator!= (const Quat<S>& q) const IMATH_NOEXCEPT;
125 
126  /// @}
127 
128 
129  /// @{
130  /// @name Query
131 
132  /// Return the R4 length
133  IMATH_HOSTDEVICE constexpr T length() const IMATH_NOEXCEPT; // in R4
134 
135  /// Return the angle of the axis/angle representation
136  IMATH_HOSTDEVICE constexpr T angle() const IMATH_NOEXCEPT;
137 
138  /// Return the axis of the axis/angle representation
139  IMATH_HOSTDEVICE constexpr Vec3<T> axis() const IMATH_NOEXCEPT;
140 
141  /// Return a 3x3 rotation matrix
142  IMATH_HOSTDEVICE constexpr Matrix33<T> toMatrix33() const IMATH_NOEXCEPT;
143 
144  /// Return a 4x4 rotation matrix
145  IMATH_HOSTDEVICE constexpr Matrix44<T> toMatrix44() const IMATH_NOEXCEPT;
146 
147  /// Return the logarithm of the quaterion
148  IMATH_HOSTDEVICE Quat<T> log() const IMATH_NOEXCEPT;
149 
150  /// Return the exponent of the quaterion
151  IMATH_HOSTDEVICE Quat<T> exp() const IMATH_NOEXCEPT;
152 
153  /// @}
154 
155  /// @{
156  /// @name Utility Methods
157 
158  /// Invert in place: this = 1 / this.
159  /// @return const reference to this.
160  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& invert() IMATH_NOEXCEPT;
161 
162  /// Return 1/this, leaving this unchanged.
163  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> inverse() const IMATH_NOEXCEPT;
164 
165  /// Normalize in place
166  /// @return const reference to this.
167  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& normalize() IMATH_NOEXCEPT;
168 
169  /// Return a normalized quaternion, leaving this unmodified.
170  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> normalized() const IMATH_NOEXCEPT;
171 
172  /// Rotate the given point by the quaterion.
173  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Vec3<T> rotateVector (const Vec3<T>& original) const IMATH_NOEXCEPT;
174 
175  /// Return the Euclidean inner product.
176  IMATH_HOSTDEVICE constexpr T euclideanInnerProduct (const Quat<T>& q) const IMATH_NOEXCEPT;
177 
178  /// Set the quaterion to be a rotation around the given axis by the
179  /// given angle.
180  /// @return const reference to this.
181  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& setAxisAngle (const Vec3<T>& axis, T radians) IMATH_NOEXCEPT;
182 
183  /// Set the quaternion to be a rotation that transforms the
184  /// direction vector `fromDirection` to `toDirection`
185  /// @return const reference to this.
186  IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>&
187  setRotation (const Vec3<T>& fromDirection, const Vec3<T>& toDirection) IMATH_NOEXCEPT;
188 
189  /// @}
190 
191  /// The base type: In templates that accept a parameter `V`, you
192  /// can refer to `T` as `V::BaseType`
193  typedef T BaseType;
194 
195  private:
196  IMATH_HOSTDEVICE void setRotationInternal (const Vec3<T>& f0, const Vec3<T>& t0, Quat<T>& q) IMATH_NOEXCEPT;
197 };
198 
199 template <class T>
200 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> slerp (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT;
201 
202 template <class T>
203 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> slerpShortestArc (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT;
204 
205 template <class T>
206 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>
207 squad (const Quat<T>& q1, const Quat<T>& q2, const Quat<T>& qa, const Quat<T>& qb, T t) IMATH_NOEXCEPT;
208 
209 ///
210 /// From advanced Animation and Rendering Techniques by Watt and Watt,
211 /// Page 366:
212 ///
213 /// computing the inner quadrangle points (qa and qb) to guarantee
214 /// tangent continuity.
215 template <class T>
216 IMATH_HOSTDEVICE void intermediate (const Quat<T>& q0,
217  const Quat<T>& q1,
218  const Quat<T>& q2,
219  const Quat<T>& q3,
220  Quat<T>& qa,
221  Quat<T>& qb) IMATH_NOEXCEPT;
222 
223 template <class T>
224 IMATH_HOSTDEVICE constexpr Matrix33<T> operator* (const Matrix33<T>& M, const Quat<T>& q) IMATH_NOEXCEPT;
225 
226 template <class T>
227 IMATH_HOSTDEVICE constexpr Matrix33<T> operator* (const Quat<T>& q, const Matrix33<T>& M) IMATH_NOEXCEPT;
228 
229 template <class T> std::ostream& operator<< (std::ostream& o, const Quat<T>& q);
230 
231 template <class T>
232 IMATH_HOSTDEVICE constexpr Quat<T> operator* (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
233 
234 template <class T>
235 IMATH_HOSTDEVICE constexpr Quat<T> operator/ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
236 
237 template <class T>
238 IMATH_HOSTDEVICE constexpr Quat<T> operator/ (const Quat<T>& q, T t) IMATH_NOEXCEPT;
239 
240 template <class T>
241 IMATH_HOSTDEVICE constexpr Quat<T> operator* (const Quat<T>& q, T t) IMATH_NOEXCEPT;
242 
243 template <class T>
244 IMATH_HOSTDEVICE constexpr Quat<T> operator* (T t, const Quat<T>& q) IMATH_NOEXCEPT;
245 
246 template <class T>
247 IMATH_HOSTDEVICE constexpr Quat<T> operator+ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
248 
249 template <class T>
250 IMATH_HOSTDEVICE constexpr Quat<T> operator- (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
251 
252 template <class T>
253 IMATH_HOSTDEVICE constexpr Quat<T> operator~ (const Quat<T>& q) IMATH_NOEXCEPT;
254 
255 template <class T>
256 IMATH_HOSTDEVICE constexpr Quat<T> operator- (const Quat<T>& q) IMATH_NOEXCEPT;
257 
258 template <class T>
259 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Vec3<T> operator* (const Vec3<T>& v, const Quat<T>& q) IMATH_NOEXCEPT;
260 
261 /// Quaternion of type float
262 typedef Quat<float> Quatf;
263 
264 /// Quaternion of type double
265 typedef Quat<double> Quatd;
266 
267 //---------------
268 // Implementation
269 //---------------
270 
271 template <class T>
272 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat() IMATH_NOEXCEPT : r (1), v (0, 0, 0)
273 {
274  // empty
275 }
276 
277 template <class T>
278 template <class S>
279 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>::Quat (const Quat<S>& q) IMATH_NOEXCEPT : r (q.r), v (q.v)
280 {
281  // empty
282 }
283 
284 template <class T>
285 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (T s, T i, T j, T k) IMATH_NOEXCEPT : r (s), v (i, j, k)
286 {
287  // empty
288 }
289 
290 template <class T>
291 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (T s, Vec3<T> d) IMATH_NOEXCEPT : r (s), v (d)
292 {
293  // empty
294 }
295 
296 template <class T>
297 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (const Quat<T>& q) IMATH_NOEXCEPT : r (q.r), v (q.v)
298 {
299  // empty
300 }
301 
302 template <class T>
303 IMATH_HOSTDEVICE constexpr inline Quat<T>
305 {
306  return Quat<T>();
307 }
308 
309 template <class T>
310 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
312 {
313  r = q.r;
314  v = q.v;
315  return *this;
316 }
317 
318 template <class T>
319 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
321 {
322  T rtmp = r * q.r - (v ^ q.v);
323  v = r * q.v + v * q.r + v % q.v;
324  r = rtmp;
325  return *this;
326 }
327 
328 template <class T>
329 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
331 {
332  r *= t;
333  v *= t;
334  return *this;
335 }
336 
337 template <class T>
338 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
340 {
341  *this = *this * q.inverse();
342  return *this;
343 }
344 
345 template <class T>
346 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
348 {
349  r /= t;
350  v /= t;
351  return *this;
352 }
353 
354 template <class T>
355 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
357 {
358  r += q.r;
359  v += q.v;
360  return *this;
361 }
362 
363 template <class T>
364 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
366 {
367  r -= q.r;
368  v -= q.v;
369  return *this;
370 }
371 
372 template <class T>
373 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T&
375 {
376  return index ? v[index - 1] : r;
377 }
378 
379 template <class T>
380 IMATH_HOSTDEVICE constexpr inline T
382 {
383  return index ? v[index - 1] : r;
384 }
385 
386 template <class T>
387 template <class S>
388 IMATH_HOSTDEVICE constexpr inline bool
390 {
391  return r == q.r && v == q.v;
392 }
393 
394 template <class T>
395 template <class S>
396 IMATH_HOSTDEVICE constexpr inline bool
398 {
399  return r != q.r || v != q.v;
400 }
401 
402 /// 4D dot product
403 template <class T>
404 IMATH_HOSTDEVICE constexpr inline T
406 {
407  return q1.r * q2.r + (q1.v ^ q2.v);
408 }
409 
410 template <class T>
411 IMATH_HOSTDEVICE constexpr inline T
413 {
414  return std::sqrt (r * r + (v ^ v));
415 }
416 
417 template <class T>
418 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
420 {
421  if (T l = length())
422  {
423  r /= l;
424  v /= l;
425  }
426  else
427  {
428  r = 1;
429  v = Vec3<T> (0);
430  }
431 
432  return *this;
433 }
434 
435 template <class T>
436 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
438 {
439  if (T l = length())
440  return Quat (r / l, v / l);
441 
442  return Quat();
443 }
444 
445 template <class T>
446 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
448 {
449  //
450  // 1 Q*
451  // - = ---- where Q* is conjugate (operator~)
452  // Q Q* Q and (Q* Q) == Q ^ Q (4D dot)
453  //
454 
455  T qdot = *this ^ *this;
456  return Quat (r / qdot, -v / qdot);
457 }
458 
459 template <class T>
460 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
462 {
463  T qdot = (*this) ^ (*this);
464  r /= qdot;
465  v = -v / qdot;
466  return *this;
467 }
468 
469 template <class T>
470 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Vec3<T>
472 {
473  //
474  // Given a vector p and a quaternion q (aka this),
475  // calculate p' = qpq*
476  //
477  // Assumes unit quaternions (because non-unit
478  // quaternions cannot be used to rotate vectors
479  // anyway).
480  //
481 
482  Quat<T> vec (0, original); // temporarily promote grade of original
483  Quat<T> inv (*this);
484  inv.v *= -1; // unit multiplicative inverse
485  Quat<T> result = *this * vec * inv;
486  return result.v;
487 }
488 
489 template <class T>
490 IMATH_HOSTDEVICE constexpr inline T
492 {
493  return r * q.r + v.x * q.v.x + v.y * q.v.y + v.z * q.v.z;
494 }
495 
496 ///
497 /// Compute the angle between two quaternions,
498 /// interpreting the quaternions as 4D vectors.
499 template <class T>
500 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T
501 angle4D (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
502 {
503  Quat<T> d = q1 - q2;
504  T lengthD = std::sqrt (d ^ d);
505 
506  Quat<T> s = q1 + q2;
507  T lengthS = std::sqrt (s ^ s);
508 
509  return 2 * std::atan2 (lengthD, lengthS);
510 }
511 
512 ///
513 /// Spherical linear interpolation.
514 /// Assumes q1 and q2 are normalized and that q1 != -q2.
515 ///
516 /// This method does *not* interpolate along the shortest
517 /// arc between q1 and q2. If you desire interpolation
518 /// along the shortest arc, and q1^q2 is negative, then
519 /// consider calling slerpShortestArc(), below, or flipping
520 /// the second quaternion explicitly.
521 ///
522 /// The implementation of squad() depends on a slerp()
523 /// that interpolates as is, without the automatic
524 /// flipping.
525 ///
526 /// Don Hatch explains the method we use here on his
527 /// web page, The Right Way to Calculate Stuff, at
528 /// http://www.plunk.org/~hatch/rightway.php
529 template <class T>
530 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
531 slerp (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT
532 {
533  T a = angle4D (q1, q2);
534  T s = 1 - t;
535 
536  Quat<T> q = sinx_over_x (s * a) / sinx_over_x (a) * s * q1 +
537  sinx_over_x (t * a) / sinx_over_x (a) * t * q2;
538 
539  return q.normalized();
540 }
541 
542 ///
543 /// Spherical linear interpolation along the shortest
544 /// arc from q1 to either q2 or -q2, whichever is closer.
545 /// Assumes q1 and q2 are unit quaternions.
546 template <class T>
547 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
549 {
550  if ((q1 ^ q2) >= 0)
551  return slerp (q1, q2, t);
552  else
553  return slerp (q1, -q2, t);
554 }
555 
556 ///
557 /// Spherical Cubic Spline Interpolation - from Advanced Animation and
558 /// Rendering Techniques by Watt and Watt, Page 366:
559 ///
560 /// A spherical curve is constructed using three spherical linear
561 /// interpolations of a quadrangle of unit quaternions: q1, qa, qb,
562 /// q2. Given a set of quaternion keys: q0, q1, q2, q3, this routine
563 /// does the interpolation between q1 and q2 by constructing two
564 /// intermediate quaternions: qa and qb. The qa and qb are computed by
565 /// the intermediate function to guarantee the continuity of tangents
566 /// across adjacent cubic segments. The qa represents in-tangent for
567 /// q1 and the qb represents the out-tangent for q2.
568 ///
569 /// The q1 q2 is the cubic segment being interpolated.
570 ///
571 /// The q0 is from the previous adjacent segment and q3 is from the
572 /// next adjacent segment. The q0 and q3 are used in computing qa and
573 /// qb.
574 template <class T>
575 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
576 spline (const Quat<T>& q0, const Quat<T>& q1, const Quat<T>& q2, const Quat<T>& q3, T t) IMATH_NOEXCEPT
577 {
578  Quat<T> qa = intermediate (q0, q1, q2);
579  Quat<T> qb = intermediate (q1, q2, q3);
580  Quat<T> result = squad (q1, qa, qb, q2, t);
581 
582  return result;
583 }
584 
585 ///
586 /// Spherical Quadrangle Interpolation - from Advanced Animation and
587 /// Rendering Techniques by Watt and Watt, Page 366:
588 ///
589 /// It constructs a spherical cubic interpolation as a series of three
590 /// spherical linear interpolations of a quadrangle of unit
591 /// quaternions.
592 template <class T>
593 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
594 squad (const Quat<T>& q1, const Quat<T>& qa, const Quat<T>& qb, const Quat<T>& q2, T t) IMATH_NOEXCEPT
595 {
596  Quat<T> r1 = slerp (q1, q2, t);
597  Quat<T> r2 = slerp (qa, qb, t);
598  Quat<T> result = slerp (r1, r2, 2 * t * (1 - t));
599 
600  return result;
601 }
602 
603 /// Compute the intermediate point between three quaternions `q0`, `q1`,
604 /// and `q2`.
605 template <class T>
606 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
607 intermediate (const Quat<T>& q0, const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
608 {
609  Quat<T> q1inv = q1.inverse();
610  Quat<T> c1 = q1inv * q2;
611  Quat<T> c2 = q1inv * q0;
612  Quat<T> c3 = (T) (-0.25) * (c2.log() + c1.log());
613  Quat<T> qa = q1 * c3.exp();
614  qa.normalize();
615  return qa;
616 }
617 
618 template <class T>
621 {
622  //
623  // For unit quaternion, from Advanced Animation and
624  // Rendering Techniques by Watt and Watt, Page 366:
625  //
626 
627  T theta = std::acos (std::min (r, (T) 1.0));
628 
629  if (theta == 0)
630  return Quat<T> (0, v);
631 
632  T sintheta = std::sin (theta);
633 
634  T k;
635  if (std::abs(sintheta) < 1 && std::abs(theta) >= std::numeric_limits<T>::max() * std::abs(sintheta))
636  k = 1;
637  else
638  k = theta / sintheta;
639 
640  return Quat<T> ((T) 0, v.x * k, v.y * k, v.z * k);
641 }
642 
643 template <class T>
646 {
647  //
648  // For pure quaternion (zero scalar part):
649  // from Advanced Animation and Rendering
650  // Techniques by Watt and Watt, Page 366:
651  //
652 
653  T theta = v.length();
654  T sintheta = std::sin (theta);
655 
656  T k;
657  if (abs (theta) < 1 && abs (sintheta) >= std::numeric_limits<T>::max() * abs (theta))
658  k = 1;
659  else
660  k = sintheta / theta;
661 
662  T costheta = std::cos (theta);
663 
664  return Quat<T> (costheta, v.x * k, v.y * k, v.z * k);
665 }
666 
667 template <class T>
668 IMATH_HOSTDEVICE constexpr inline T
670 {
671  return 2 * std::atan2 (v.length(), r);
672 }
673 
674 template <class T>
675 IMATH_HOSTDEVICE constexpr inline Vec3<T>
677 {
678  return v.normalized();
679 }
680 
681 template <class T>
682 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
684 {
685  r = std::cos (radians / 2);
686  v = axis.normalized() * std::sin (radians / 2);
687  return *this;
688 }
689 
690 template <class T>
691 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
693 {
694  //
695  // Create a quaternion that rotates vector from into vector to,
696  // such that the rotation is around an axis that is the cross
697  // product of from and to.
698  //
699  // This function calls function setRotationInternal(), which is
700  // numerically accurate only for rotation angles that are not much
701  // greater than pi/2. In order to achieve good accuracy for angles
702  // greater than pi/2, we split large angles in half, and rotate in
703  // two steps.
704  //
705 
706  //
707  // Normalize from and to, yielding f0 and t0.
708  //
709 
710  Vec3<T> f0 = from.normalized();
711  Vec3<T> t0 = to.normalized();
712 
713  if ((f0 ^ t0) >= 0)
714  {
715  //
716  // The rotation angle is less than or equal to pi/2.
717  //
718 
719  setRotationInternal (f0, t0, *this);
720  }
721  else
722  {
723  //
724  // The angle is greater than pi/2. After computing h0,
725  // which is halfway between f0 and t0, we rotate first
726  // from f0 to h0, then from h0 to t0.
727  //
728 
729  Vec3<T> h0 = (f0 + t0).normalized();
730 
731  if ((h0 ^ h0) != 0)
732  {
733  setRotationInternal (f0, h0, *this);
734 
735  Quat<T> q;
736  setRotationInternal (h0, t0, q);
737 
738  *this *= q;
739  }
740  else
741  {
742  //
743  // f0 and t0 point in exactly opposite directions.
744  // Pick an arbitrary axis that is orthogonal to f0,
745  // and rotate by pi.
746  //
747 
748  r = T (0);
749 
750  Vec3<T> f02 = f0 * f0;
751 
752  if (f02.x <= f02.y && f02.x <= f02.z)
753  v = (f0 % Vec3<T> (1, 0, 0)).normalized();
754  else if (f02.y <= f02.z)
755  v = (f0 % Vec3<T> (0, 1, 0)).normalized();
756  else
757  v = (f0 % Vec3<T> (0, 0, 1)).normalized();
758  }
759  }
760 
761  return *this;
762 }
763 
764 template <class T>
765 IMATH_HOSTDEVICE inline void
767 {
768  //
769  // The following is equivalent to setAxisAngle(n,2*phi),
770  // where the rotation axis, n, is orthogonal to the f0 and
771  // t0 vectors, and 2*phi is the angle between f0 and t0.
772  //
773  // This function is called by setRotation(), above; it assumes
774  // that f0 and t0 are normalized and that the angle between
775  // them is not much greater than pi/2. This function becomes
776  // numerically inaccurate if f0 and t0 point into nearly
777  // opposite directions.
778  //
779 
780  //
781  // Find a normalized vector, h0, that is halfway between f0 and t0.
782  // The angle between f0 and h0 is phi.
783  //
784 
785  Vec3<T> h0 = (f0 + t0).normalized();
786 
787  //
788  // Store the rotation axis and rotation angle.
789  //
790 
791  q.r = f0 ^ h0; // f0 ^ h0 == cos (phi)
792  q.v = f0 % h0; // (f0 % h0).length() == sin (phi)
793 }
794 
795 template <class T>
796 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
798 {
799  return Matrix33<T> (1 - 2 * (v.y * v.y + v.z * v.z),
800  2 * (v.x * v.y + v.z * r),
801  2 * (v.z * v.x - v.y * r),
802 
803  2 * (v.x * v.y - v.z * r),
804  1 - 2 * (v.z * v.z + v.x * v.x),
805  2 * (v.y * v.z + v.x * r),
806 
807  2 * (v.z * v.x + v.y * r),
808  2 * (v.y * v.z - v.x * r),
809  1 - 2 * (v.y * v.y + v.x * v.x));
810 }
811 
812 template <class T>
813 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
815 {
816  return Matrix44<T> (1 - 2 * (v.y * v.y + v.z * v.z),
817  2 * (v.x * v.y + v.z * r),
818  2 * (v.z * v.x - v.y * r),
819  0,
820  2 * (v.x * v.y - v.z * r),
821  1 - 2 * (v.z * v.z + v.x * v.x),
822  2 * (v.y * v.z + v.x * r),
823  0,
824  2 * (v.z * v.x + v.y * r),
825  2 * (v.y * v.z - v.x * r),
826  1 - 2 * (v.y * v.y + v.x * v.x),
827  0,
828  0,
829  0,
830  0,
831  1);
832 }
833 
834 /// Transform the quaternion by the matrix
835 /// @return M * q
836 template <class T>
837 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
839 {
840  return M * q.toMatrix33();
841 }
842 
843 /// Transform the matrix by the quaterion:
844 /// @return q * M
845 template <class T>
846 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
848 {
849  return q.toMatrix33() * M;
850 }
851 
852 /// Stream output as "(r x y z)"
853 template <class T>
854 std::ostream&
855 operator<< (std::ostream& o, const Quat<T>& q)
856 {
857  return o << "(" << q.r << " " << q.v.x << " " << q.v.y << " " << q.v.z << ")";
858 }
859 
860 /// Quaterion multiplication
861 template <class T>
862 IMATH_HOSTDEVICE constexpr inline Quat<T>
864 {
865  return Quat<T> (q1.r * q2.r - (q1.v ^ q2.v), q1.r * q2.v + q1.v * q2.r + q1.v % q2.v);
866 }
867 
868 /// Quaterion division
869 template <class T>
870 IMATH_HOSTDEVICE constexpr inline Quat<T>
872 {
873  return q1 * q2.inverse();
874 }
875 
876 /// Quaterion division
877 template <class T>
878 IMATH_HOSTDEVICE constexpr inline Quat<T>
880 {
881  return Quat<T> (q.r / t, q.v / t);
882 }
883 
884 /// Quaterion*scalar multiplication
885 /// @return q * t
886 template <class T>
887 IMATH_HOSTDEVICE constexpr inline Quat<T>
889 {
890  return Quat<T> (q.r * t, q.v * t);
891 }
892 
893 /// Quaterion*scalar multiplication
894 /// @return q * t
895 template <class T>
896 IMATH_HOSTDEVICE constexpr inline Quat<T>
898 {
899  return Quat<T> (q.r * t, q.v * t);
900 }
901 
902 /// Quaterion addition
903 template <class T>
904 IMATH_HOSTDEVICE constexpr inline Quat<T>
906 {
907  return Quat<T> (q1.r + q2.r, q1.v + q2.v);
908 }
909 
910 /// Quaterion subtraction
911 template <class T>
912 IMATH_HOSTDEVICE constexpr inline Quat<T>
914 {
915  return Quat<T> (q1.r - q2.r, q1.v - q2.v);
916 }
917 
918 /// Compute the conjugate
919 template <class T>
920 IMATH_HOSTDEVICE constexpr inline Quat<T>
922 {
923  return Quat<T> (q.r, -q.v);
924 }
925 
926 /// Negate the quaterion
927 template <class T>
928 IMATH_HOSTDEVICE constexpr inline Quat<T>
930 {
931  return Quat<T> (-q.r, -q.v);
932 }
933 
934 /// Quaterion*vector multiplcation
935 /// @return v * q
936 template <class T>
937 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Vec3<T>
939 {
940  Vec3<T> a = q.v % v;
941  Vec3<T> b = q.v % a;
942  return v + T (2) * (q.r * a + b);
943 }
944 
945 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
946 # pragma warning(pop)
947 #endif
948 
949 IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
950 
951 #endif // INCLUDED_IMATHQUAT_H
IMATH_HOSTDEVICE constexpr T euclideanInnerProduct(const Quat< T > &q) const IMATH_NOEXCEPT
Return the Euclidean inner product.
Definition: ImathQuat.h:491
SYS_API double cos(double x)
Definition: SYS_FPUMath.h:69
T BaseType
Definition: ImathQuat.h:193
SYS_API double atan2(double y, double x)
Definition: SYS_FPUMath.h:79
T z
Definition: ImathVec.h:310
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat< T > & operator=(const Quat< T > &q) IMATH_NOEXCEPT
Assignment.
Definition: ImathQuat.h:311
#define IMATH_NOEXCEPT
Definition: ImathConfig.h:72
Definition: ImathVec.h:32
*get result *(waiting if necessary)*A common idiom is to fire a bunch of sub tasks at the and then *wait for them to all complete We provide a helper class
Definition: thread.h:623
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > & normalize() IMATH_NOEXCEPT
Definition: ImathQuat.h:419
Definition: ImathQuat.h:42
SIM_API const UT_StringHolder angle
GLboolean invert
Definition: glcorearb.h:549
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat< T > & operator+=(const Quat< T > &q) IMATH_NOEXCEPT
Quaternion addition.
Definition: ImathQuat.h:356
const GLdouble * v
Definition: glcorearb.h:837
IMATH_HOSTDEVICE constexpr Quat< T > operator-(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
Quaterion subtraction.
Definition: ImathQuat.h:913
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7481
GLboolean GLboolean GLboolean GLboolean a
Definition: glcorearb.h:1222
GLdouble s
Definition: glad.h:3009
IMATH_HOSTDEVICE Vec3< T > normalized() const IMATH_NOEXCEPT
Return a normalized vector. Does not modify *this.
Definition: ImathVec.h:1801
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > squad(const Quat< T > &q1, const Quat< T > &q2, const Quat< T > &qa, const Quat< T > &qb, T t) IMATH_NOEXCEPT
Definition: ImathQuat.h:594
IMATH_HOSTDEVICE static constexpr Quat< T > identity() IMATH_NOEXCEPT
The identity quaternion.
Definition: ImathQuat.h:304
Vec3< T > v
The imaginary vector.
Definition: ImathQuat.h:53
ImageBuf OIIO_API min(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > spline(const Quat< T > &q0, const Quat< T > &q1, const Quat< T > &q2, const Quat< T > &q3, T t) IMATH_NOEXCEPT
Definition: ImathQuat.h:576
IMATH_HOSTDEVICE Quat< T > exp() const IMATH_NOEXCEPT
Return the exponent of the quaterion.
Definition: ImathQuat.h:645
**But if you need a result
Definition: thread.h:613
IMATH_HOSTDEVICE constexpr T operator^(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
4D dot product
Definition: ImathQuat.h:405
IMATH_HOSTDEVICE constexpr T angle() const IMATH_NOEXCEPT
Return the angle of the axis/angle representation.
Definition: ImathQuat.h:669
IMATH_HOSTDEVICE constexpr bool operator==(const Quat< S > &q) const IMATH_NOEXCEPT
Equality.
Definition: ImathQuat.h:389
GLdouble GLdouble GLdouble q
Definition: glad.h:2445
#define IMATH_HOSTDEVICE
Definition: ImathConfig.h:102
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T angle4D(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
Definition: ImathQuat.h:501
IMATH_HOSTDEVICE Quat< T > log() const IMATH_NOEXCEPT
Return the logarithm of the quaterion.
Definition: ImathQuat.h:620
T r
The real part.
Definition: ImathQuat.h:50
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat< T > & operator/=(const Quat< T > &q) IMATH_NOEXCEPT
Quaterion division, using the inverse()
Definition: ImathQuat.h:339
IMATH_HOSTDEVICE constexpr bool operator!=(const Quat< S > &q) const IMATH_NOEXCEPT
Inequality.
Definition: ImathQuat.h:397
T x
Definition: ImathVec.h:310
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat< T > & operator-=(const Quat< T > &q) IMATH_NOEXCEPT
Quaternion subtraction.
Definition: ImathQuat.h:365
IMATH_HOSTDEVICE constexpr T length() const IMATH_NOEXCEPT
Return the R4 length.
Definition: ImathQuat.h:412
IMATH_HOSTDEVICE constexpr Quat() IMATH_NOEXCEPT
Default constructor is the identity quat.
Definition: ImathQuat.h:272
IMATH_HOSTDEVICE constexpr Matrix44< T > toMatrix44() const IMATH_NOEXCEPT
Return a 4x4 rotation matrix.
Definition: ImathQuat.h:814
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > & setRotation(const Vec3< T > &fromDirection, const Vec3< T > &toDirection) IMATH_NOEXCEPT
Definition: ImathQuat.h:692
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > & invert() IMATH_NOEXCEPT
Definition: ImathQuat.h:461
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T & operator[](int index) IMATH_NOEXCEPT
Definition: ImathQuat.h:374
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > & setAxisAngle(const Vec3< T > &axis, T radians) IMATH_NOEXCEPT
Definition: ImathQuat.h:683
SYS_API double acos(double x)
Definition: SYS_FPUMath.h:83
GLint GLenum GLboolean normalized
Definition: glcorearb.h:872
GLboolean GLboolean GLboolean b
Definition: glcorearb.h:1222
IMATH_HOSTDEVICE constexpr Quat< T > operator+(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
Quaterion addition.
Definition: ImathQuat.h:905
GLdouble t
Definition: glad.h:2397
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > slerp(const Quat< T > &q1, const Quat< T > &q2, T t) IMATH_NOEXCEPT
Definition: ImathQuat.h:531
GLint j
Definition: glad.h:2733
IMATH_HOSTDEVICE constexpr Quat< T > operator~(const Quat< T > &q) IMATH_NOEXCEPT
Compute the conjugate.
Definition: ImathQuat.h:921
IMATH_HOSTDEVICE constexpr Matrix33< T > operator*(const Matrix33< T > &M, const Quat< T > &q) IMATH_NOEXCEPT
Definition: ImathQuat.h:838
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > inverse() const IMATH_NOEXCEPT
Return 1/this, leaving this unchanged.
Definition: ImathQuat.h:447
GLuint index
Definition: glcorearb.h:786
#define IMATH_EXPORT_TEMPLATE_TYPE
Definition: ImathExport.h:60
ImageBuf OIIO_API max(Image_or_Const A, Image_or_Const B, ROI roi={}, int nthreads=0)
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Vec3< T > rotateVector(const Vec3< T > &original) const IMATH_NOEXCEPT
Rotate the given point by the quaterion.
Definition: ImathQuat.h:471
T y
Definition: ImathVec.h:310
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE constexpr T abs(T a) IMATH_NOEXCEPT
Definition: ImathFun.h:26
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER IMATH_HOSTDEVICE T sinx_over_x(T x)
Definition: ImathMath.h:136
IMATH_HOSTDEVICE constexpr Vec3< T > axis() const IMATH_NOEXCEPT
Return the axis of the axis/angle representation.
Definition: ImathQuat.h:676
GLboolean r
Definition: glcorearb.h:1222
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat< T > & operator*=(const Quat< T > &q) IMATH_NOEXCEPT
Quaternion multiplication.
Definition: ImathQuat.h:320
IMATH_HOSTDEVICE void intermediate(const Quat< T > &q0, const Quat< T > &q1, const Quat< T > &q2, const Quat< T > &q3, Quat< T > &qa, Quat< T > &qb) IMATH_NOEXCEPT
IMATH_HOSTDEVICE constexpr Matrix33< T > toMatrix33() const IMATH_NOEXCEPT
Return a 3x3 rotation matrix.
Definition: ImathQuat.h:797
OIIO_FORCEINLINE T log(const T &v)
Definition: simd.h:7688
IMATH_HOSTDEVICE constexpr Quat< T > operator/(const Quat< T > &q1, const Quat< T > &q2) IMATH_NOEXCEPT
Quaterion division.
Definition: ImathQuat.h:871
OIIO_FORCEINLINE OIIO_HOSTDEVICE T radians(T deg)
Convert degrees to radians.
Definition: fmath.h:669
SYS_API double sin(double x)
Definition: SYS_FPUMath.h:71
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > normalized() const IMATH_NOEXCEPT
Return a normalized quaternion, leaving this unmodified.
Definition: ImathQuat.h:437
constexpr T normalize(UT_FixedVector< T, D > &a) noexcept
IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat< T > slerpShortestArc(const Quat< T > &q1, const Quat< T > &q2, T t) IMATH_NOEXCEPT
Definition: ImathQuat.h:548