PageRenderTime 60ms CodeModel.GetById 23ms app.highlight 32ms RepoModel.GetById 1ms app.codeStats 0ms

/indra/llmath/llquaternion.h

https://bitbucket.org/lindenlab/viewer-beta/
C++ Header | 588 lines | 381 code | 80 blank | 127 comment | 27 complexity | 2901803c13d86b8f3acf0c3a76f5b305 MD5 | raw file
  1/** 
  2 * @file llquaternion.h
  3 * @brief LLQuaternion class header file.
  4 *
  5 * $LicenseInfo:firstyear=2000&license=viewerlgpl$
  6 * Second Life Viewer Source Code
  7 * Copyright (C) 2010, Linden Research, Inc.
  8 * 
  9 * This library is free software; you can redistribute it and/or
 10 * modify it under the terms of the GNU Lesser General Public
 11 * License as published by the Free Software Foundation;
 12 * version 2.1 of the License only.
 13 * 
 14 * This library is distributed in the hope that it will be useful,
 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 17 * Lesser General Public License for more details.
 18 * 
 19 * You should have received a copy of the GNU Lesser General Public
 20 * License along with this library; if not, write to the Free Software
 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 22 * 
 23 * Linden Research, Inc., 945 Battery Street, San Francisco, CA  94111  USA
 24 * $/LicenseInfo$
 25 */
 26
 27#ifndef LLQUATERNION_H
 28#define LLQUATERNION_H
 29
 30#include <iostream>
 31
 32#ifndef LLMATH_H //enforce specific include order to avoid tangling inline dependencies
 33#error "Please include llmath.h first."
 34#endif
 35
 36class LLVector4;
 37class LLVector3;
 38class LLVector3d;
 39class LLMatrix4;
 40class LLMatrix3;
 41
 42//	NOTA BENE: Quaternion code is written assuming Unit Quaternions!!!!
 43//			   Moreover, it is written assuming that all vectors and matricies
 44//			   passed as arguments are normalized and unitary respectively.
 45//			   VERY VERY VERY VERY BAD THINGS will happen if these assumptions fail.
 46
 47static const U32 LENGTHOFQUAT = 4;
 48
 49class LLQuaternion
 50{
 51public:
 52	F32 mQ[LENGTHOFQUAT];
 53
 54	static const LLQuaternion DEFAULT;
 55
 56	LLQuaternion();									// Initializes Quaternion to (0,0,0,1)
 57	explicit LLQuaternion(const LLMatrix4 &mat);				// Initializes Quaternion from Matrix4
 58	explicit LLQuaternion(const LLMatrix3 &mat);				// Initializes Quaternion from Matrix3
 59	LLQuaternion(F32 x, F32 y, F32 z, F32 w);		// Initializes Quaternion to normalize(x, y, z, w)
 60	LLQuaternion(F32 angle, const LLVector4 &vec);	// Initializes Quaternion to axis_angle2quat(angle, vec)
 61	LLQuaternion(F32 angle, const LLVector3 &vec);	// Initializes Quaternion to axis_angle2quat(angle, vec)
 62	LLQuaternion(const F32 *q);						// Initializes Quaternion to normalize(x, y, z, w)
 63	LLQuaternion(const LLVector3 &x_axis,
 64				 const LLVector3 &y_axis,
 65				 const LLVector3 &z_axis);			// Initializes Quaternion from Matrix3 = [x_axis ; y_axis ; z_axis]
 66
 67	BOOL isIdentity() const;
 68	BOOL isNotIdentity() const;
 69	BOOL isFinite() const;									// checks to see if all values of LLQuaternion are finite
 70	void quantize16(F32 lower, F32 upper);					// changes the vector to reflect quatization
 71	void quantize8(F32 lower, F32 upper);							// changes the vector to reflect quatization
 72	void loadIdentity();											// Loads the quaternion that represents the identity rotation
 73
 74	const LLQuaternion&	set(F32 x, F32 y, F32 z, F32 w);		// Sets Quaternion to normalize(x, y, z, w)
 75	const LLQuaternion&	set(const LLQuaternion &quat);			// Copies Quaternion
 76	const LLQuaternion&	set(const F32 *q);						// Sets Quaternion to normalize(quat[VX], quat[VY], quat[VZ], quat[VW])
 77	const LLQuaternion&	set(const LLMatrix3 &mat);				// Sets Quaternion to mat2quat(mat)
 78	const LLQuaternion&	set(const LLMatrix4 &mat);				// Sets Quaternion to mat2quat(mat)
 79
 80	const LLQuaternion&	setAngleAxis(F32 angle, F32 x, F32 y, F32 z);	// Sets Quaternion to axis_angle2quat(angle, x, y, z)
 81	const LLQuaternion&	setAngleAxis(F32 angle, const LLVector3 &vec);	// Sets Quaternion to axis_angle2quat(angle, vec)
 82	const LLQuaternion&	setAngleAxis(F32 angle, const LLVector4 &vec);	// Sets Quaternion to axis_angle2quat(angle, vec)
 83	const LLQuaternion&	setEulerAngles(F32 roll, F32 pitch, F32 yaw);	// Sets Quaternion to euler2quat(pitch, yaw, roll)
 84
 85	const LLQuaternion&	setQuatInit(F32 x, F32 y, F32 z, F32 w);	// deprecated
 86	const LLQuaternion&	setQuat(const LLQuaternion &quat);			// deprecated
 87	const LLQuaternion&	setQuat(const F32 *q);						// deprecated
 88	const LLQuaternion&	setQuat(const LLMatrix3 &mat);				// deprecated
 89	const LLQuaternion&	setQuat(const LLMatrix4 &mat);				// deprecated
 90	const LLQuaternion&	setQuat(F32 angle, F32 x, F32 y, F32 z);	// deprecated
 91	const LLQuaternion&	setQuat(F32 angle, const LLVector3 &vec);	// deprecated
 92	const LLQuaternion&	setQuat(F32 angle, const LLVector4 &vec);	// deprecated
 93	const LLQuaternion&	setQuat(F32 roll, F32 pitch, F32 yaw);		// deprecated
 94
 95	LLMatrix4	getMatrix4(void) const;							// Returns the Matrix4 equivalent of Quaternion
 96	LLMatrix3	getMatrix3(void) const;							// Returns the Matrix3 equivalent of Quaternion
 97	void		getAngleAxis(F32* angle, F32* x, F32* y, F32* z) const;	// returns rotation in radians about axis x,y,z
 98	void		getAngleAxis(F32* angle, LLVector3 &vec) const;
 99	void		getEulerAngles(F32 *roll, F32* pitch, F32 *yaw) const;
100
101	F32	normalize();	// Normalizes Quaternion and returns magnitude
102	F32	normQuat();		// deprecated
103
104	const LLQuaternion&	conjugate(void);	// Conjugates Quaternion and returns result
105	const LLQuaternion&	conjQuat(void);		// deprecated
106
107	// Other useful methods
108	const LLQuaternion&	transpose();		// transpose (same as conjugate)
109	const LLQuaternion&	transQuat();		// deprecated
110
111	void			shortestArc(const LLVector3 &a, const LLVector3 &b);	// shortest rotation from a to b
112	const LLQuaternion& constrain(F32 radians);						// constrains rotation to a cone angle specified in radians
113
114	// Standard operators
115	friend std::ostream& operator<<(std::ostream &s, const LLQuaternion &a);					// Prints a
116	friend LLQuaternion operator+(const LLQuaternion &a, const LLQuaternion &b);	// Addition
117	friend LLQuaternion operator-(const LLQuaternion &a, const LLQuaternion &b);	// Subtraction
118	friend LLQuaternion operator-(const LLQuaternion &a);							// Negation
119	friend LLQuaternion operator*(F32 a, const LLQuaternion &q);					// Scale
120	friend LLQuaternion operator*(const LLQuaternion &q, F32 b);					// Scale
121	friend LLQuaternion operator*(const LLQuaternion &a, const LLQuaternion &b);	// Returns a * b
122	friend LLQuaternion operator~(const LLQuaternion &a);							// Returns a* (Conjugate of a)
123	bool operator==(const LLQuaternion &b) const;			// Returns a == b
124	bool operator!=(const LLQuaternion &b) const;			// Returns a != b
125
126	friend const LLQuaternion& operator*=(LLQuaternion &a, const LLQuaternion &b);	// Returns a * b
127
128	friend LLVector4 operator*(const LLVector4 &a, const LLQuaternion &rot);		// Rotates a by rot
129	friend LLVector3 operator*(const LLVector3 &a, const LLQuaternion &rot);		// Rotates a by rot
130	friend LLVector3d operator*(const LLVector3d &a, const LLQuaternion &rot);		// Rotates a by rot
131
132	// Non-standard operators
133	friend F32 dot(const LLQuaternion &a, const LLQuaternion &b);
134	friend LLQuaternion lerp(F32 t, const LLQuaternion &p, const LLQuaternion &q);		// linear interpolation (t = 0 to 1) from p to q
135	friend LLQuaternion lerp(F32 t, const LLQuaternion &q);								// linear interpolation (t = 0 to 1) from identity to q
136	friend LLQuaternion slerp(F32 t, const LLQuaternion &p, const LLQuaternion &q); 	// spherical linear interpolation from p to q
137	friend LLQuaternion slerp(F32 t, const LLQuaternion &q);							// spherical linear interpolation from identity to q
138	friend LLQuaternion nlerp(F32 t, const LLQuaternion &p, const LLQuaternion &q); 	// normalized linear interpolation from p to q
139	friend LLQuaternion nlerp(F32 t, const LLQuaternion &q); 							// normalized linear interpolation from p to q
140
141	LLVector3	packToVector3() const;						// Saves space by using the fact that our quaternions are normalized
142	void		unpackFromVector3(const LLVector3& vec);	// Saves space by using the fact that our quaternions are normalized
143
144	enum Order {
145		XYZ = 0,
146		YZX = 1,
147		ZXY = 2,
148		XZY = 3,
149		YXZ = 4,
150		ZYX = 5
151	};
152	// Creates a quaternions from maya's rotation representation,
153	// which is 3 rotations (in DEGREES) in the specified order
154	friend LLQuaternion mayaQ(F32 x, F32 y, F32 z, Order order);
155
156	// Conversions between Order and strings like "xyz" or "ZYX"
157	friend const char *OrderToString( const Order order );
158	friend Order StringToOrder( const char *str );
159
160	static BOOL parseQuat(const std::string& buf, LLQuaternion* value);
161
162	// For debugging, only
163	//static U32 mMultCount;
164};
165
166// checker
167inline BOOL	LLQuaternion::isFinite() const
168{
169	return (llfinite(mQ[VX]) && llfinite(mQ[VY]) && llfinite(mQ[VZ]) && llfinite(mQ[VS]));
170}
171
172inline BOOL LLQuaternion::isIdentity() const
173{
174	return 
175		( mQ[VX] == 0.f ) &&
176		( mQ[VY] == 0.f ) &&
177		( mQ[VZ] == 0.f ) &&
178		( mQ[VS] == 1.f );
179}
180
181inline BOOL LLQuaternion::isNotIdentity() const
182{
183	return 
184		( mQ[VX] != 0.f ) ||
185		( mQ[VY] != 0.f ) ||
186		( mQ[VZ] != 0.f ) ||
187		( mQ[VS] != 1.f );
188}
189
190
191
192inline LLQuaternion::LLQuaternion(void)
193{
194	mQ[VX] = 0.f;
195	mQ[VY] = 0.f;
196	mQ[VZ] = 0.f;
197	mQ[VS] = 1.f;
198}
199
200inline LLQuaternion::LLQuaternion(F32 x, F32 y, F32 z, F32 w)
201{
202	mQ[VX] = x;
203	mQ[VY] = y;
204	mQ[VZ] = z;
205	mQ[VS] = w;
206
207	//RN: don't normalize this case as its used mainly for temporaries during calculations
208	//normalize();
209	/*
210	F32 mag = sqrtf(mQ[VX]*mQ[VX] + mQ[VY]*mQ[VY] + mQ[VZ]*mQ[VZ] + mQ[VS]*mQ[VS]);
211	mag -= 1.f;
212	mag = fabs(mag);
213	llassert(mag < 10.f*FP_MAG_THRESHOLD);
214	*/
215}
216
217inline LLQuaternion::LLQuaternion(const F32 *q)
218{
219	mQ[VX] = q[VX];
220	mQ[VY] = q[VY];
221	mQ[VZ] = q[VZ];
222	mQ[VS] = q[VW];
223
224	normalize();
225	/*
226	F32 mag = sqrtf(mQ[VX]*mQ[VX] + mQ[VY]*mQ[VY] + mQ[VZ]*mQ[VZ] + mQ[VS]*mQ[VS]);
227	mag -= 1.f;
228	mag = fabs(mag);
229	llassert(mag < FP_MAG_THRESHOLD);
230	*/
231}
232
233
234inline void LLQuaternion::loadIdentity()
235{
236	mQ[VX] = 0.0f;
237	mQ[VY] = 0.0f;
238	mQ[VZ] = 0.0f;
239	mQ[VW] = 1.0f;
240}
241
242
243inline const LLQuaternion&	LLQuaternion::set(F32 x, F32 y, F32 z, F32 w)
244{
245	mQ[VX] = x;
246	mQ[VY] = y;
247	mQ[VZ] = z;
248	mQ[VS] = w;
249	normalize();
250	return (*this);
251}
252
253inline const LLQuaternion&	LLQuaternion::set(const LLQuaternion &quat)
254{
255	mQ[VX] = quat.mQ[VX];
256	mQ[VY] = quat.mQ[VY];
257	mQ[VZ] = quat.mQ[VZ];
258	mQ[VW] = quat.mQ[VW];
259	normalize();
260	return (*this);
261}
262
263inline const LLQuaternion&	LLQuaternion::set(const F32 *q)
264{
265	mQ[VX] = q[VX];
266	mQ[VY] = q[VY];
267	mQ[VZ] = q[VZ];
268	mQ[VS] = q[VW];
269	normalize();
270	return (*this);
271}
272
273
274// deprecated
275inline const LLQuaternion&	LLQuaternion::setQuatInit(F32 x, F32 y, F32 z, F32 w)
276{
277	mQ[VX] = x;
278	mQ[VY] = y;
279	mQ[VZ] = z;
280	mQ[VS] = w;
281	normalize();
282	return (*this);
283}
284
285// deprecated
286inline const LLQuaternion&	LLQuaternion::setQuat(const LLQuaternion &quat)
287{
288	mQ[VX] = quat.mQ[VX];
289	mQ[VY] = quat.mQ[VY];
290	mQ[VZ] = quat.mQ[VZ];
291	mQ[VW] = quat.mQ[VW];
292	normalize();
293	return (*this);
294}
295
296// deprecated
297inline const LLQuaternion&	LLQuaternion::setQuat(const F32 *q)
298{
299	mQ[VX] = q[VX];
300	mQ[VY] = q[VY];
301	mQ[VZ] = q[VZ];
302	mQ[VS] = q[VW];
303	normalize();
304	return (*this);
305}
306
307// There may be a cheaper way that avoids the sqrt.
308// Does sin_a = VX*VX + VY*VY + VZ*VZ?
309// Copied from Matrix and Quaternion FAQ 1.12
310inline void LLQuaternion::getAngleAxis(F32* angle, F32* x, F32* y, F32* z) const
311{
312	F32 cos_a = mQ[VW];
313	if (cos_a > 1.0f) cos_a = 1.0f;
314	if (cos_a < -1.0f) cos_a = -1.0f;
315
316    F32 sin_a = (F32) sqrt( 1.0f - cos_a * cos_a );
317
318    if ( fabs( sin_a ) < 0.0005f )
319		sin_a = 1.0f;
320	else
321		sin_a = 1.f/sin_a;
322
323    F32 temp_angle = 2.0f * (F32) acos( cos_a );
324	if (temp_angle > F_PI)
325	{
326		// The (angle,axis) pair should never have angles outside [PI, -PI]
327		// since we want the _shortest_ (angle,axis) solution.
328		// Since acos is defined for [0, PI], and we multiply by 2.0, we
329		// can push the angle outside the acceptible range.
330		// When this happens we set the angle to the other portion of a 
331		// full 2PI rotation, and negate the axis, which reverses the 
332		// direction of the rotation (by the right-hand rule).
333		*angle = 2.f * F_PI - temp_angle;
334    	*x = - mQ[VX] * sin_a;
335    	*y = - mQ[VY] * sin_a;
336    	*z = - mQ[VZ] * sin_a;
337	}
338	else
339	{
340		*angle = temp_angle;
341    	*x = mQ[VX] * sin_a;
342    	*y = mQ[VY] * sin_a;
343    	*z = mQ[VZ] * sin_a;
344	}
345}
346
347inline const LLQuaternion& LLQuaternion::conjugate()
348{
349	mQ[VX] *= -1.f;
350	mQ[VY] *= -1.f;
351	mQ[VZ] *= -1.f;
352	return (*this);
353}
354
355inline const LLQuaternion& LLQuaternion::conjQuat()
356{
357	mQ[VX] *= -1.f;
358	mQ[VY] *= -1.f;
359	mQ[VZ] *= -1.f;
360	return (*this);
361}
362
363// Transpose
364inline const LLQuaternion& LLQuaternion::transpose()
365{
366	mQ[VX] *= -1.f;
367	mQ[VY] *= -1.f;
368	mQ[VZ] *= -1.f;
369	return (*this);
370}
371
372// deprecated
373inline const LLQuaternion& LLQuaternion::transQuat()
374{
375	mQ[VX] *= -1.f;
376	mQ[VY] *= -1.f;
377	mQ[VZ] *= -1.f;
378	return (*this);
379}
380
381
382inline LLQuaternion 	operator+(const LLQuaternion &a, const LLQuaternion &b)
383{
384	return LLQuaternion( 
385		a.mQ[VX] + b.mQ[VX],
386		a.mQ[VY] + b.mQ[VY],
387		a.mQ[VZ] + b.mQ[VZ],
388		a.mQ[VW] + b.mQ[VW] );
389}
390
391
392inline LLQuaternion 	operator-(const LLQuaternion &a, const LLQuaternion &b)
393{
394	return LLQuaternion( 
395		a.mQ[VX] - b.mQ[VX],
396		a.mQ[VY] - b.mQ[VY],
397		a.mQ[VZ] - b.mQ[VZ],
398		a.mQ[VW] - b.mQ[VW] );
399}
400
401
402inline LLQuaternion 	operator-(const LLQuaternion &a)
403{
404	return LLQuaternion(
405		-a.mQ[VX],
406		-a.mQ[VY],
407		-a.mQ[VZ],
408		-a.mQ[VW] );
409}
410
411
412inline LLQuaternion 	operator*(F32 a, const LLQuaternion &q)
413{
414	return LLQuaternion(
415		a * q.mQ[VX],
416		a * q.mQ[VY],
417		a * q.mQ[VZ],
418		a * q.mQ[VW] );
419}
420
421
422inline LLQuaternion 	operator*(const LLQuaternion &q, F32 a)
423{
424	return LLQuaternion(
425		a * q.mQ[VX],
426		a * q.mQ[VY],
427		a * q.mQ[VZ],
428		a * q.mQ[VW] );
429}
430
431inline LLQuaternion	operator~(const LLQuaternion &a)
432{
433	LLQuaternion q(a);
434	q.conjQuat();
435	return q;
436}
437
438inline bool	LLQuaternion::operator==(const LLQuaternion &b) const
439{
440	return (  (mQ[VX] == b.mQ[VX])
441			&&(mQ[VY] == b.mQ[VY])
442			&&(mQ[VZ] == b.mQ[VZ])
443			&&(mQ[VS] == b.mQ[VS]));
444}
445
446inline bool	LLQuaternion::operator!=(const LLQuaternion &b) const
447{
448	return (  (mQ[VX] != b.mQ[VX])
449			||(mQ[VY] != b.mQ[VY])
450			||(mQ[VZ] != b.mQ[VZ])
451			||(mQ[VS] != b.mQ[VS]));
452}
453
454inline const LLQuaternion&	operator*=(LLQuaternion &a, const LLQuaternion &b)
455{
456#if 1
457	LLQuaternion q(
458		b.mQ[3] * a.mQ[0] + b.mQ[0] * a.mQ[3] + b.mQ[1] * a.mQ[2] - b.mQ[2] * a.mQ[1],
459		b.mQ[3] * a.mQ[1] + b.mQ[1] * a.mQ[3] + b.mQ[2] * a.mQ[0] - b.mQ[0] * a.mQ[2],
460		b.mQ[3] * a.mQ[2] + b.mQ[2] * a.mQ[3] + b.mQ[0] * a.mQ[1] - b.mQ[1] * a.mQ[0],
461		b.mQ[3] * a.mQ[3] - b.mQ[0] * a.mQ[0] - b.mQ[1] * a.mQ[1] - b.mQ[2] * a.mQ[2]
462	);
463	a = q;
464#else
465	a = a * b;
466#endif
467	return a;
468}
469
470const F32 ONE_PART_IN_A_MILLION = 0.000001f;
471
472inline F32	LLQuaternion::normalize()
473{
474	F32 mag = sqrtf(mQ[VX]*mQ[VX] + mQ[VY]*mQ[VY] + mQ[VZ]*mQ[VZ] + mQ[VS]*mQ[VS]);
475
476	if (mag > FP_MAG_THRESHOLD)
477	{
478		// Floating point error can prevent some quaternions from achieving
479		// exact unity length.  When trying to renormalize such quaternions we
480		// can oscillate between multiple quantized states.  To prevent such
481		// drifts we only renomalize if the length is far enough from unity.
482		if (fabs(1.f - mag) > ONE_PART_IN_A_MILLION)
483		{
484			F32 oomag = 1.f/mag;
485			mQ[VX] *= oomag;
486			mQ[VY] *= oomag;
487			mQ[VZ] *= oomag;
488			mQ[VS] *= oomag;
489		}
490	}
491	else
492	{
493		// we were given a very bad quaternion so we set it to identity
494		mQ[VX] = 0.f;
495		mQ[VY] = 0.f;
496		mQ[VZ] = 0.f;
497		mQ[VS] = 1.f;
498	}
499
500	return mag;
501}
502
503// deprecated
504inline F32	LLQuaternion::normQuat()
505{
506	F32 mag = sqrtf(mQ[VX]*mQ[VX] + mQ[VY]*mQ[VY] + mQ[VZ]*mQ[VZ] + mQ[VS]*mQ[VS]);
507
508	if (mag > FP_MAG_THRESHOLD)
509	{
510		if (fabs(1.f - mag) > ONE_PART_IN_A_MILLION)
511		{
512			// only renormalize if length not close enough to 1.0 already
513			F32 oomag = 1.f/mag;
514			mQ[VX] *= oomag;
515			mQ[VY] *= oomag;
516			mQ[VZ] *= oomag;
517			mQ[VS] *= oomag;
518		}
519	}
520	else
521	{
522		mQ[VX] = 0.f;
523		mQ[VY] = 0.f;
524		mQ[VZ] = 0.f;
525		mQ[VS] = 1.f;
526	}
527
528	return mag;
529}
530
531LLQuaternion::Order StringToOrder( const char *str );
532
533// Some notes about Quaternions
534
535// What is a Quaternion?
536// ---------------------
537// A quaternion is a point in 4-dimensional complex space.
538// Q = { Qx, Qy, Qz, Qw }
539// 
540//
541// Why Quaternions?
542// ----------------
543// The set of quaternions that make up the the 4-D unit sphere 
544// can be mapped to the set of all rotations in 3-D space.  Sometimes
545// it is easier to describe/manipulate rotations in quaternion space
546// than rotation-matrix space.
547//
548//
549// How Quaternions?
550// ----------------
551// In order to take advantage of quaternions we need to know how to
552// go from rotation-matricies to quaternions and back.  We also have
553// to agree what variety of rotations we're generating.
554// 
555// Consider the equation...   v' = v * R 
556//
557// There are two ways to think about rotations of vectors.
558// 1) v' is the same vector in a different reference frame
559// 2) v' is a new vector in the same reference frame
560//
561// bookmark -- which way are we using?
562// 
563// 
564// Quaternion from Angle-Axis:
565// ---------------------------
566// Suppose we wanted to represent a rotation of some angle (theta) 
567// about some axis ({Ax, Ay, Az})...
568//
569// axis of rotation = {Ax, Ay, Az} 
570// angle_of_rotation = theta
571//
572// s = sin(0.5 * theta)
573// c = cos(0.5 * theta)
574// Q = { s * Ax, s * Ay, s * Az, c }
575//
576//
577// 3x3 Matrix from Quaternion
578// --------------------------
579//
580//     |                                                                    |
581//     | 1 - 2 * (y^2 + z^2)   2 * (x * y + z * w)     2 * (y * w - x * z)  |
582//     |                                                                    |
583// M = | 2 * (x * y - z * w)   1 - 2 * (x^2 + z^2)     2 * (y * z + x * w)  |
584//     |                                                                    |
585//     | 2 * (x * z + y * w)   2 * (y * z - x * w)     1 - 2 * (x^2 + y^2)  |
586//     |                                                                    |
587
588#endif