2 // SPDX-License-Identifier: BSD-3-Clause
3 // Copyright Contributors to the OpenEXR Project.
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
17 #ifndef INCLUDED_IMATHQUAT_H
18 #define INCLUDED_IMATHQUAT_H
20 #include "ImathExport.h"
21 #include "ImathNamespace.h"
23 #include "ImathMatrix.h"
27 IMATH_INTERNAL_NAMESPACE_HEADER_ENTER
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)
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.
42 template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Quat
47 /// @name Direct access to elements
52 /// The imaginary vector
57 /// Element access: q[0] is the real part, (q[1],q[2],q[3]) is the
59 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T& operator[] (int index) IMATH_NOEXCEPT; // as 4D vector
61 /// Element access: q[0] is the real part, (q[1],q[2],q[3]) is the
63 IMATH_HOSTDEVICE constexpr T operator[] (int index) const IMATH_NOEXCEPT;
66 /// @name Constructors
68 /// Default constructor is the identity quat
69 IMATH_HOSTDEVICE constexpr Quat() IMATH_NOEXCEPT;
72 IMATH_HOSTDEVICE constexpr Quat (const Quat& q) IMATH_NOEXCEPT;
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;
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;
80 /// Initialize with real part `s` and imaginary vector `d`
81 IMATH_HOSTDEVICE constexpr Quat (T s, Vec3<T> d) IMATH_NOEXCEPT;
83 /// The identity quaternion
84 IMATH_HOSTDEVICE constexpr static Quat<T> identity() IMATH_NOEXCEPT;
87 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator= (const Quat<T>& q) IMATH_NOEXCEPT;
90 ~Quat () IMATH_NOEXCEPT = default;
95 /// @name Basic Algebra
97 /// Note that the operator return values are *NOT* normalized
100 /// Quaternion multiplication
101 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator*= (const Quat<T>& q) IMATH_NOEXCEPT;
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;
107 /// Quaterion division, using the inverse()
108 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator/= (const Quat<T>& q) IMATH_NOEXCEPT;
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;
114 /// Quaternion addition
115 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator+= (const Quat<T>& q) IMATH_NOEXCEPT;
117 /// Quaternion subtraction
118 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 const Quat<T>& operator-= (const Quat<T>& q) IMATH_NOEXCEPT;
121 template <class S> IMATH_HOSTDEVICE constexpr bool operator== (const Quat<S>& q) const IMATH_NOEXCEPT;
124 template <class S> IMATH_HOSTDEVICE constexpr bool operator!= (const Quat<S>& q) const IMATH_NOEXCEPT;
132 /// Return the R4 length
133 IMATH_HOSTDEVICE constexpr T length() const IMATH_NOEXCEPT; // in R4
135 /// Return the angle of the axis/angle representation
136 IMATH_HOSTDEVICE constexpr T angle() const IMATH_NOEXCEPT;
138 /// Return the axis of the axis/angle representation
139 IMATH_HOSTDEVICE constexpr Vec3<T> axis() const IMATH_NOEXCEPT;
141 /// Return a 3x3 rotation matrix
142 IMATH_HOSTDEVICE constexpr Matrix33<T> toMatrix33() const IMATH_NOEXCEPT;
144 /// Return a 4x4 rotation matrix
145 IMATH_HOSTDEVICE constexpr Matrix44<T> toMatrix44() const IMATH_NOEXCEPT;
147 /// Return the logarithm of the quaterion
148 IMATH_HOSTDEVICE Quat<T> log() const IMATH_NOEXCEPT;
150 /// Return the exponent of the quaterion
151 IMATH_HOSTDEVICE Quat<T> exp() const IMATH_NOEXCEPT;
156 /// @name Utility Methods
158 /// Invert in place: this = 1 / this.
159 /// @return const reference to this.
160 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& invert() IMATH_NOEXCEPT;
162 /// Return 1/this, leaving this unchanged.
163 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> inverse() const IMATH_NOEXCEPT;
165 /// Normalize in place
166 /// @return const reference to this.
167 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& normalize() IMATH_NOEXCEPT;
169 /// Return a normalized quaternion, leaving this unmodified.
170 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> normalized() const IMATH_NOEXCEPT;
172 /// Rotate the given point by the quaterion.
173 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Vec3<T> rotateVector (const Vec3<T>& original) const IMATH_NOEXCEPT;
175 /// Return the Euclidean inner product.
176 IMATH_HOSTDEVICE constexpr T euclideanInnerProduct (const Quat<T>& q) const IMATH_NOEXCEPT;
178 /// Set the quaterion to be a rotation around the given axis by the
180 /// @return const reference to this.
181 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T>& setAxisAngle (const Vec3<T>& axis, T radians) IMATH_NOEXCEPT;
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;
191 /// The base type: In templates that accept a parameter `V`, you
192 /// can refer to `T` as `V::BaseType`
196 IMATH_HOSTDEVICE void setRotationInternal (const Vec3<T>& f0, const Vec3<T>& t0, Quat<T>& q) IMATH_NOEXCEPT;
200 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> slerp (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT;
203 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Quat<T> slerpShortestArc (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT;
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;
210 /// From advanced Animation and Rendering Techniques by Watt and Watt,
213 /// computing the inner quadrangle points (qa and qb) to guarantee
214 /// tangent continuity.
216 IMATH_HOSTDEVICE void intermediate (const Quat<T>& q0,
221 Quat<T>& qb) IMATH_NOEXCEPT;
224 IMATH_HOSTDEVICE constexpr Matrix33<T> operator* (const Matrix33<T>& M, const Quat<T>& q) IMATH_NOEXCEPT;
227 IMATH_HOSTDEVICE constexpr Matrix33<T> operator* (const Quat<T>& q, const Matrix33<T>& M) IMATH_NOEXCEPT;
229 template <class T> std::ostream& operator<< (std::ostream& o, const Quat<T>& q);
232 IMATH_HOSTDEVICE constexpr Quat<T> operator* (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
235 IMATH_HOSTDEVICE constexpr Quat<T> operator/ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
238 IMATH_HOSTDEVICE constexpr Quat<T> operator/ (const Quat<T>& q, T t) IMATH_NOEXCEPT;
241 IMATH_HOSTDEVICE constexpr Quat<T> operator* (const Quat<T>& q, T t) IMATH_NOEXCEPT;
244 IMATH_HOSTDEVICE constexpr Quat<T> operator* (T t, const Quat<T>& q) IMATH_NOEXCEPT;
247 IMATH_HOSTDEVICE constexpr Quat<T> operator+ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
250 IMATH_HOSTDEVICE constexpr Quat<T> operator- (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT;
253 IMATH_HOSTDEVICE constexpr Quat<T> operator~ (const Quat<T>& q) IMATH_NOEXCEPT;
256 IMATH_HOSTDEVICE constexpr Quat<T> operator- (const Quat<T>& q) IMATH_NOEXCEPT;
259 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 Vec3<T> operator* (const Vec3<T>& v, const Quat<T>& q) IMATH_NOEXCEPT;
261 /// Quaternion of type float
262 typedef Quat<float> Quatf;
264 /// Quaternion of type double
265 typedef Quat<double> Quatd;
272 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat() IMATH_NOEXCEPT : r (1), v (0, 0, 0)
279 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>::Quat (const Quat<S>& q) IMATH_NOEXCEPT : r (q.r), v (q.v)
285 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (T s, T i, T j, T k) IMATH_NOEXCEPT : r (s), v (i, j, k)
291 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (T s, Vec3<T> d) IMATH_NOEXCEPT : r (s), v (d)
297 IMATH_HOSTDEVICE constexpr inline Quat<T>::Quat (const Quat<T>& q) IMATH_NOEXCEPT : r (q.r), v (q.v)
303 IMATH_HOSTDEVICE constexpr inline Quat<T>
304 Quat<T>::identity() IMATH_NOEXCEPT
310 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
311 Quat<T>::operator= (const Quat<T>& q) IMATH_NOEXCEPT
319 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
320 Quat<T>::operator*= (const Quat<T>& q) IMATH_NOEXCEPT
322 T rtmp = r * q.r - (v ^ q.v);
323 v = r * q.v + v * q.r + v % q.v;
329 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
330 Quat<T>::operator*= (T t) IMATH_NOEXCEPT
338 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
339 Quat<T>::operator/= (const Quat<T>& q) IMATH_NOEXCEPT
341 *this = *this * q.inverse();
346 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
347 Quat<T>::operator/= (T t) IMATH_NOEXCEPT
355 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
356 Quat<T>::operator+= (const Quat<T>& q) IMATH_NOEXCEPT
364 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline const Quat<T>&
365 Quat<T>::operator-= (const Quat<T>& q) IMATH_NOEXCEPT
373 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T&
374 Quat<T>::operator[] (int index) IMATH_NOEXCEPT
376 return index ? v[index - 1] : r;
380 IMATH_HOSTDEVICE constexpr inline T
381 Quat<T>::operator[] (int index) const IMATH_NOEXCEPT
383 return index ? v[index - 1] : r;
388 IMATH_HOSTDEVICE constexpr inline bool
389 Quat<T>::operator== (const Quat<S>& q) const IMATH_NOEXCEPT
391 return r == q.r && v == q.v;
396 IMATH_HOSTDEVICE constexpr inline bool
397 Quat<T>::operator!= (const Quat<S>& q) const IMATH_NOEXCEPT
399 return r != q.r || v != q.v;
404 IMATH_HOSTDEVICE constexpr inline T
405 operator^ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
407 return q1.r * q2.r + (q1.v ^ q2.v);
411 IMATH_HOSTDEVICE constexpr inline T
412 Quat<T>::length() const IMATH_NOEXCEPT
414 return std::sqrt (r * r + (v ^ v));
418 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
419 Quat<T>::normalize() IMATH_NOEXCEPT
436 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
437 Quat<T>::normalized() const IMATH_NOEXCEPT
440 return Quat (r / l, v / l);
446 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
447 Quat<T>::inverse() const IMATH_NOEXCEPT
451 // - = ---- where Q* is conjugate (operator~)
452 // Q Q* Q and (Q* Q) == Q ^ Q (4D dot)
455 T qdot = *this ^ *this;
456 return Quat (r / qdot, -v / qdot);
460 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
461 Quat<T>::invert() IMATH_NOEXCEPT
463 T qdot = (*this) ^ (*this);
470 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Vec3<T>
471 Quat<T>::rotateVector (const Vec3<T>& original) const IMATH_NOEXCEPT
474 // Given a vector p and a quaternion q (aka this),
475 // calculate p' = qpq*
477 // Assumes unit quaternions (because non-unit
478 // quaternions cannot be used to rotate vectors
482 Quat<T> vec (0, original); // temporarily promote grade of original
484 inv.v *= -1; // unit multiplicative inverse
485 Quat<T> result = *this * vec * inv;
490 IMATH_HOSTDEVICE constexpr inline T
491 Quat<T>::euclideanInnerProduct (const Quat<T>& q) const IMATH_NOEXCEPT
493 return r * q.r + v.x * q.v.x + v.y * q.v.y + v.z * q.v.z;
497 /// Compute the angle between two quaternions,
498 /// interpreting the quaternions as 4D vectors.
500 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline T
501 angle4D (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
504 T lengthD = std::sqrt (d ^ d);
507 T lengthS = std::sqrt (s ^ s);
509 return 2 * std::atan2 (lengthD, lengthS);
513 /// Spherical linear interpolation.
514 /// Assumes q1 and q2 are normalized and that q1 != -q2.
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.
522 /// The implementation of squad() depends on a slerp()
523 /// that interpolates as is, without the automatic
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
530 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
531 slerp (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT
533 T a = angle4D (q1, q2);
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;
539 return q.normalized();
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.
547 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
548 slerpShortestArc (const Quat<T>& q1, const Quat<T>& q2, T t) IMATH_NOEXCEPT
551 return slerp (q1, q2, t);
553 return slerp (q1, -q2, t);
557 /// Spherical Cubic Spline Interpolation - from Advanced Animation and
558 /// Rendering Techniques by Watt and Watt, Page 366:
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.
569 /// The q1 q2 is the cubic segment being interpolated.
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
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
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);
586 /// Spherical Quadrangle Interpolation - from Advanced Animation and
587 /// Rendering Techniques by Watt and Watt, Page 366:
589 /// It constructs a spherical cubic interpolation as a series of three
590 /// spherical linear interpolations of a quadrangle of unit
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
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));
603 /// Compute the intermediate point between three quaternions `q0`, `q1`,
606 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>
607 intermediate (const Quat<T>& q0, const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
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();
619 IMATH_HOSTDEVICE inline Quat<T>
620 Quat<T>::log() const IMATH_NOEXCEPT
623 // For unit quaternion, from Advanced Animation and
624 // Rendering Techniques by Watt and Watt, Page 366:
627 T theta = std::acos (std::min (r, (T) 1.0));
630 return Quat<T> (0, v);
632 T sintheta = std::sin (theta);
635 if (std::abs(sintheta) < 1 && std::abs(theta) >= std::numeric_limits<T>::max() * std::abs(sintheta))
638 k = theta / sintheta;
640 return Quat<T> ((T) 0, v.x * k, v.y * k, v.z * k);
644 IMATH_HOSTDEVICE inline Quat<T>
645 Quat<T>::exp() const IMATH_NOEXCEPT
648 // For pure quaternion (zero scalar part):
649 // from Advanced Animation and Rendering
650 // Techniques by Watt and Watt, Page 366:
653 T theta = v.length();
654 T sintheta = std::sin (theta);
657 if (abs (theta) < 1 && abs (sintheta) >= std::numeric_limits<T>::max() * abs (theta))
660 k = sintheta / theta;
662 T costheta = std::cos (theta);
664 return Quat<T> (costheta, v.x * k, v.y * k, v.z * k);
668 IMATH_HOSTDEVICE constexpr inline T
669 Quat<T>::angle() const IMATH_NOEXCEPT
671 return 2 * std::atan2 (v.length(), r);
675 IMATH_HOSTDEVICE constexpr inline Vec3<T>
676 Quat<T>::axis() const IMATH_NOEXCEPT
678 return v.normalized();
682 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
683 Quat<T>::setAxisAngle (const Vec3<T>& axis, T radians) IMATH_NOEXCEPT
685 r = std::cos (radians / 2);
686 v = axis.normalized() * std::sin (radians / 2);
691 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Quat<T>&
692 Quat<T>::setRotation (const Vec3<T>& from, const Vec3<T>& to) IMATH_NOEXCEPT
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.
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
707 // Normalize from and to, yielding f0 and t0.
710 Vec3<T> f0 = from.normalized();
711 Vec3<T> t0 = to.normalized();
716 // The rotation angle is less than or equal to pi/2.
719 setRotationInternal (f0, t0, *this);
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.
729 Vec3<T> h0 = (f0 + t0).normalized();
733 setRotationInternal (f0, h0, *this);
736 setRotationInternal (h0, t0, q);
743 // f0 and t0 point in exactly opposite directions.
744 // Pick an arbitrary axis that is orthogonal to f0,
750 Vec3<T> f02 = f0 * f0;
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();
757 v = (f0 % Vec3<T> (0, 0, 1)).normalized();
765 IMATH_HOSTDEVICE inline void
766 Quat<T>::setRotationInternal (const Vec3<T>& f0, const Vec3<T>& t0, Quat<T>& q) IMATH_NOEXCEPT
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.
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.
781 // Find a normalized vector, h0, that is halfway between f0 and t0.
782 // The angle between f0 and h0 is phi.
785 Vec3<T> h0 = (f0 + t0).normalized();
788 // Store the rotation axis and rotation angle.
791 q.r = f0 ^ h0; // f0 ^ h0 == cos (phi)
792 q.v = f0 % h0; // (f0 % h0).length() == sin (phi)
796 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
797 Quat<T>::toMatrix33() const IMATH_NOEXCEPT
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),
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),
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));
813 IMATH_HOSTDEVICE constexpr inline Matrix44<T>
814 Quat<T>::toMatrix44() const IMATH_NOEXCEPT
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),
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),
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),
834 /// Transform the quaternion by the matrix
837 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
838 operator* (const Matrix33<T>& M, const Quat<T>& q) IMATH_NOEXCEPT
840 return M * q.toMatrix33();
843 /// Transform the matrix by the quaterion:
846 IMATH_HOSTDEVICE constexpr inline Matrix33<T>
847 operator* (const Quat<T>& q, const Matrix33<T>& M) IMATH_NOEXCEPT
849 return q.toMatrix33() * M;
852 /// Stream output as "(r x y z)"
855 operator<< (std::ostream& o, const Quat<T>& q)
857 return o << "(" << q.r << " " << q.v.x << " " << q.v.y << " " << q.v.z << ")";
860 /// Quaterion multiplication
862 IMATH_HOSTDEVICE constexpr inline Quat<T>
863 operator* (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
865 return Quat<T> (q1.r * q2.r - (q1.v ^ q2.v), q1.r * q2.v + q1.v * q2.r + q1.v % q2.v);
868 /// Quaterion division
870 IMATH_HOSTDEVICE constexpr inline Quat<T>
871 operator/ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
873 return q1 * q2.inverse();
876 /// Quaterion division
878 IMATH_HOSTDEVICE constexpr inline Quat<T>
879 operator/ (const Quat<T>& q, T t) IMATH_NOEXCEPT
881 return Quat<T> (q.r / t, q.v / t);
884 /// Quaterion*scalar multiplication
887 IMATH_HOSTDEVICE constexpr inline Quat<T>
888 operator* (const Quat<T>& q, T t) IMATH_NOEXCEPT
890 return Quat<T> (q.r * t, q.v * t);
893 /// Quaterion*scalar multiplication
896 IMATH_HOSTDEVICE constexpr inline Quat<T>
897 operator* (T t, const Quat<T>& q) IMATH_NOEXCEPT
899 return Quat<T> (q.r * t, q.v * t);
902 /// Quaterion addition
904 IMATH_HOSTDEVICE constexpr inline Quat<T>
905 operator+ (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
907 return Quat<T> (q1.r + q2.r, q1.v + q2.v);
910 /// Quaterion subtraction
912 IMATH_HOSTDEVICE constexpr inline Quat<T>
913 operator- (const Quat<T>& q1, const Quat<T>& q2) IMATH_NOEXCEPT
915 return Quat<T> (q1.r - q2.r, q1.v - q2.v);
918 /// Compute the conjugate
920 IMATH_HOSTDEVICE constexpr inline Quat<T>
921 operator~ (const Quat<T>& q) IMATH_NOEXCEPT
923 return Quat<T> (q.r, -q.v);
926 /// Negate the quaterion
928 IMATH_HOSTDEVICE constexpr inline Quat<T>
929 operator- (const Quat<T>& q) IMATH_NOEXCEPT
931 return Quat<T> (-q.r, -q.v);
934 /// Quaterion*vector multiplcation
937 IMATH_HOSTDEVICE IMATH_CONSTEXPR14 inline Vec3<T>
938 operator* (const Vec3<T>& v, const Quat<T>& q) IMATH_NOEXCEPT
942 return v + T (2) * (q.r * a + b);
945 #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER
946 # pragma warning(pop)
949 IMATH_INTERNAL_NAMESPACE_HEADER_EXIT
951 #endif // INCLUDED_IMATHQUAT_H