PageRenderTime 77ms CodeModel.GetById 2ms app.highlight 67ms RepoModel.GetById 2ms app.codeStats 0ms

/src/FreeImage/Source/OpenEXR/Imath/ImathQuat.h

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