/////////////////////////////////////////////////////////////////////////// // // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas // Digital Ltd. LLC // // All rights reserved. // // Redistribution and use in source and binary forms, with or without // modification, are permitted provided that the following conditions are // met: // * Redistributions of source code must retain the above copyright // notice, this list of conditions and the following disclaimer. // * Redistributions in binary form must reproduce the above // copyright notice, this list of conditions and the following disclaimer // in the documentation and/or other materials provided with the // distribution. // * Neither the name of Industrial Light & Magic nor the names of // its contributors may be used to endorse or promote products derived // from this software without specific prior written permission. // // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. // /////////////////////////////////////////////////////////////////////////// #ifndef INCLUDED_IMATHQUAT_H #define INCLUDED_IMATHQUAT_H //---------------------------------------------------------------------- // // template class Quat // // "Quaternions came from Hamilton ... and have been an unmixed // evil to those who have touched them in any way. Vector is a // useless survival ... and has never been of the slightest use // to any creature." // // - Lord Kelvin // // This class implements the quaternion numerical type -- you // will probably want to use this class to represent orientations // in R3 and to convert between various euler angle reps. You // should probably use Imath::Euler<> for that. // //---------------------------------------------------------------------- #include "ImathExc.h" #include "ImathMatrix.h" #include namespace Imath { #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER // Disable MS VC++ warnings about conversion from double to float #pragma warning(disable:4244) #endif template class Quat; template Quat slerp (const Quat &q1,const Quat &q2, T t); template Quat squad (const Quat &q1,const Quat &q2, const Quat &qa,const Quat &qb, T t); template void intermediate (const Quat &q0, const Quat &q1, const Quat &q2, const Quat &q3, Quat &qa, Quat &qb); template class Quat { public: T r; // real part Vec3 v; // imaginary vector //----------------------------------------------------- // Constructors - default constructor is identity quat //----------------------------------------------------- Quat() : r(1), v(0,0,0) {} template Quat( const Quat& q) : r(q.r), v(q.v) {} Quat( T s, T i, T j, T k ) : r(s), v(i,j,k) {} Quat( T s, Vec3 d ) : r(s), v(d) {} static Quat identity() { return Quat(); } //------------------------------------------------ // Basic Algebra - Operators and Methods // The operator return values are *NOT* normalized // // operator^ is 4D dot product // operator/ uses the inverse() quaternion // operator~ is conjugate -- if (S+V) is quat then // the conjugate (S+V)* == (S-V) // // some operators (*,/,*=,/=) treat the quat as // a 4D vector when one of the operands is scalar //------------------------------------------------ const Quat& operator= (const Quat&); const Quat& operator*= (const Quat&); const Quat& operator*= (T); const Quat& operator/= (const Quat&); const Quat& operator/= (T); const Quat& operator+= (const Quat&); const Quat& operator-= (const Quat&); T& operator[] (int index); // as 4D vector T operator[] (int index) const; template bool operator == (const Quat &q) const; template bool operator != (const Quat &q) const; Quat& invert(); // this -> 1 / this Quat inverse() const; Quat& normalize(); // returns this Quat normalized() const; T length() const; // in R4 //----------------------- // Rotation conversion //----------------------- Quat& setAxisAngle(const Vec3& axis, T radians); Quat& setRotation(const Vec3& fromDirection, const Vec3& toDirection); T angle() const; Vec3 axis() const; Matrix33 toMatrix33() const; Matrix44 toMatrix44() const; Quat log() const; Quat exp() const; }; //-------------------- // Convenient typedefs //-------------------- typedef Quat Quatf; typedef Quat Quatd; //--------------- // Implementation //--------------- template inline const Quat& Quat::operator= (const Quat& q) { r = q.r; v = q.v; return *this; } template inline const Quat& Quat::operator*= (const Quat& q) { T rtmp = r * q.r - (v ^ q.v); v = r * q.v + v * q.r + v % q.v; r = rtmp; return *this; } template inline const Quat& Quat::operator*= (T t) { r *= t; v *= t; return *this; } template inline const Quat& Quat::operator/= (const Quat& q) { *this = *this * q.inverse(); return *this; } template inline const Quat& Quat::operator/= (T t) { r /= t; v /= t; return *this; } template inline const Quat& Quat::operator+= (const Quat& q) { r += q.r; v += q.v; return *this; } template inline const Quat& Quat::operator-= (const Quat& q) { r -= q.r; v -= q.v; return *this; } template inline T& Quat::operator[] (int index) { return index ? v[index-1] : r; } template inline T Quat::operator[] (int index) const { return index ? v[index-1] : r; } template template inline bool Quat::operator == (const Quat &q) const { return r == q.r && v == q.v; } template template inline bool Quat::operator != (const Quat &q) const { return r != q.r || v != q.v; } template inline T operator^ (const Quat& q1,const Quat& q2) { return q1.r * q2.r + (q1.v ^ q2.v); } template inline T Quat::length() const { return Math::sqrt( r * r + (v ^ v) ); } template inline Quat& Quat::normalize() { if ( T l = length() ) { r /= l; v /= l; } else { r = 1; v = Vec3(0); } return *this; } template inline Quat Quat::normalized() const { if ( T l = length() ) return Quat( r / l, v / l ); return Quat(); } template inline Quat Quat::inverse() const { // 1 Q* // - = ---- where Q* is conjugate (operator~) // Q Q* Q and (Q* Q) == Q ^ Q (4D dot) T qdot = *this ^ *this; return Quat( r / qdot, -v / qdot ); } template inline Quat& Quat::invert() { T qdot = (*this) ^ (*this); r /= qdot; v = -v / qdot; return *this; } template Quat slerp(const Quat &q1,const Quat &q2, T t) { // // Spherical linear interpolation. // // NOTE: Assumes q1 and q2 are normalized and that 0 <= t <= 1. // // This method does *not* interpolate along the shortest arc // between q1 and q2. If you desire interpolation along the // shortest arc, then consider flipping the second quaternion // explicitly before calling slerp. The implementation of squad() // depends on a slerp() that interpolates as is, without the // automatic flipping. // T cosomega = q1 ^ q2; if (cosomega >= (T) 1.0) { // // Special case: q1 and q2 are the same, so just return one of them. // This also catches the case where cosomega is very slightly > 1.0 // return q1; } T sinomega = Math::sqrt (1 - cosomega * cosomega); Quat result; if (sinomega * limits::max() > 1) { T omega = Math::acos (cosomega); T s1 = Math::sin ((1.0 - t) * omega) / sinomega; T s2 = Math::sin (t * omega) / sinomega; result = s1 * q1 + s2 * q2; } else if (cosomega > 0) { // // omega == 0 // T s1 = 1.0 - t; T s2 = t; result = s1 * q1 + s2 * q2; } else { // // omega == -pi // result.v.x = - q1.v.y; result.v.y = q1.v.x; result.v.z = - q1.r; result.r = q1.v.z; T s1 = Math::sin ((0.5 - t) * M_PI); T s2 = Math::sin (t * M_PI); result = s1 * q1 + s2 * result; } return result; } template Quat spline(const Quat &q0, const Quat &q1, const Quat &q2, const Quat &q3, T t) { // Spherical Cubic Spline Interpolation - // from Advanced Animation and Rendering // Techniques by Watt and Watt, Page 366: // A spherical curve is constructed using three // spherical linear interpolations of a quadrangle // of unit quaternions: q1, qa, qb, q2. // Given a set of quaternion keys: q0, q1, q2, q3, // this routine does the interpolation between // q1 and q2 by constructing two intermediate // quaternions: qa and qb. The qa and qb are // computed by the intermediate function to // guarantee the continuity of tangents across // adjacent cubic segments. The qa represents in-tangent // for q1 and the qb represents the out-tangent for q2. // // The q1 q2 is the cubic segment being interpolated. // The q0 is from the previous adjacent segment and q3 is // from the next adjacent segment. The q0 and q3 are used // in computing qa and qb. // Quat qa = intermediate (q0, q1, q2); Quat qb = intermediate (q1, q2, q3); Quat result = squad(q1, qa, qb, q2, t); return result; } template Quat squad(const Quat &q1, const Quat &qa, const Quat &qb, const Quat &q2, T t) { // Spherical Quadrangle Interpolation - // from Advanced Animation and Rendering // Techniques by Watt and Watt, Page 366: // It constructs a spherical cubic interpolation as // a series of three spherical linear interpolations // of a quadrangle of unit quaternions. // Quat r1 = slerp(q1, q2, t); Quat r2 = slerp(qa, qb, t); Quat result = slerp(r1, r2, 2*t*(1-t)); return result; } template Quat intermediate(const Quat &q0, const Quat &q1, const Quat &q2) { // From advanced Animation and Rendering // Techniques by Watt and Watt, Page 366: // computing the inner quadrangle // points (qa and qb) to guarantee tangent // continuity. // Quat q1inv = q1.inverse(); Quat c1 = q1inv*q2; Quat c2 = q1inv*q0; Quat c3 = (T) (-0.25) * (c2.log() + c1.log()); Quat qa = q1 * c3.exp(); qa.normalize(); return qa; } template inline Quat Quat::log() const { // // For unit quaternion, from Advanced Animation and // Rendering Techniques by Watt and Watt, Page 366: // T theta = Math::acos (std::min (r, (T) 1.0)); if (theta == 0) return Quat (0, v); T sintheta = Math::sin (theta); T k; if (abs (sintheta) < 1 && abs (theta) >= limits::max() * abs (sintheta)) k = 0; else k = theta / sintheta; return Quat ((T) 0, v.x * k, v.y * k, v.z * k); } template inline Quat Quat::exp() const { // // For pure quaternion (zero scalar part): // from Advanced Animation and Rendering // Techniques by Watt and Watt, Page 366: // T theta = v.length(); T sintheta = Math::sin (theta); T k; if (abs (theta) < 1 && abs (sintheta) >= limits::max() * abs (theta)) k = 0; else k = sintheta / theta; T costheta = Math::cos (theta); return Quat (costheta, v.x * k, v.y * k, v.z * k); } template inline T Quat::angle() const { return 2.0*Math::acos(r); } template inline Vec3 Quat::axis() const { return v.normalized(); } template inline Quat& Quat::setAxisAngle(const Vec3& axis, T radians) { r = Math::cos(radians/2); v = axis.normalized() * Math::sin(radians/2); return *this; } template Quat& Quat::setRotation(const Vec3& from, const Vec3& to) { // // Ported from SbRotation // T cost = from.dot(to) / Math::sqrt(from.dot(from) * to.dot(to)); // check for degeneracies if (cost > 0.99999) { // // Vectors are parallel. // r = 1.0; v = Vec3(0); } else if (cost < -0.99999) { // // Vectors are opposite. Find an axis to rotate around, // which should be perpendicular to the original axis. // Vec3 frm = from.normalized(); v = frm.cross(Vec3(1, 0, 0)); if (v.length() < 0.00001) v = frm.cross(Vec3(0, 1, 0)); r = 0; v.normalize(); } else { // // Use half-angle formulae: // cos^2 t = ( 1 + cos (2t) ) / 2 // w part is cosine of half the rotation angle // r = Math::sqrt(0.5 * (1.0 + cost)); // // sin^2 t = ( 1 - cos (2t) ) / 2 // Do the normalization of the axis vector at the same time so // we only call sqrt once. // v = from.cross(to); v *= Math::sqrt((0.5 * (1.0 - cost))/(v.dot(v))); } return *this; } template Matrix33 Quat::toMatrix33() const { return Matrix33(1. - 2.0 * (v.y * v.y + v.z * v.z), 2.0 * (v.x * v.y + v.z * r), 2.0 * (v.z * v.x - v.y * r), 2.0 * (v.x * v.y - v.z * r), 1. - 2.0 * (v.z * v.z + v.x * v.x), 2.0 * (v.y * v.z + v.x * r), 2.0 * (v.z * v.x + v.y * r), 2.0 * (v.y * v.z - v.x * r), 1. - 2.0 * (v.y * v.y + v.x * v.x)); } template Matrix44 Quat::toMatrix44() const { return Matrix44(1. - 2.0 * (v.y * v.y + v.z * v.z), 2.0 * (v.x * v.y + v.z * r), 2.0 * (v.z * v.x - v.y * r), 0., 2.0 * (v.x * v.y - v.z * r), 1. - 2.0 * (v.z * v.z + v.x * v.x), 2.0 * (v.y * v.z + v.x * r), 0., 2.0 * (v.z * v.x + v.y * r), 2.0 * (v.y * v.z - v.x * r), 1. - 2.0 * (v.y * v.y + v.x * v.x), 0., 0., 0., 0., 1.0 ); } template inline Matrix33 operator* (const Matrix33 &M, const Quat &q) { return M * q.toMatrix33(); } template inline Matrix33 operator* (const Quat &q, const Matrix33 &M) { return q.toMatrix33() * M; } template std::ostream& operator<< (std::ostream &o, const Quat &q) { return o << "(" << q.r << " " << q.v.x << " " << q.v.y << " " << q.v.z << ")"; } template inline Quat operator* (const Quat& q1, const Quat& q2) { // (S1+V1) (S2+V2) = S1 S2 - V1.V2 + S1 V2 + V1 S2 + V1 x V2 return Quat( q1.r * q2.r - (q1.v ^ q2.v), q1.r * q2.v + q1.v * q2.r + q1.v % q2.v ); } template inline Quat operator/ (const Quat& q1, const Quat& q2) { return q1 * q2.inverse(); } template inline Quat operator/ (const Quat& q,T t) { return Quat(q.r/t,q.v/t); } template inline Quat operator* (const Quat& q,T t) { return Quat(q.r*t,q.v*t); } template inline Quat operator* (T t, const Quat& q) { return Quat(q.r*t,q.v*t); } template inline Quat operator+ (const Quat& q1, const Quat& q2) { return Quat( q1.r + q2.r, q1.v + q2.v ); } template inline Quat operator- (const Quat& q1, const Quat& q2) { return Quat( q1.r - q2.r, q1.v - q2.v ); } template inline Quat operator~ (const Quat& q) { return Quat( q.r, -q.v ); // conjugate: (S+V)* = S-V } template inline Quat operator- (const Quat& q) { return Quat( -q.r, -q.v ); } #if (defined _WIN32 || defined _WIN64) && defined _MSC_VER #pragma warning(default:4244) #endif } // namespace Imath #endif