Path: utzoo!attcan!uunet!lll-winken!arisia!sgi!gavin@krypton.SGI.COM
From: gavin@krypton.SGI.COM (Gavin Bell)
Newsgroups: comp.graphics
Subject: Re: 3D Rotations/Instancing
Summary: Using Euler paramaters for rotations
Message-ID: <25633@sgi.SGI.COM>
Date: 27 Jan 89 02:11:01 GMT
References: <65@sdcc10.ucsd.EDU> <3390@cloud9.Stratus.COM> <3391@cloud9.Stratus.COM>
Sender: daemon@sgi.SGI.COM
Organization: Silicon Graphics, Inc., Mountain View, CA
Lines: 123


	Here's another way of keeping error out of your rotation matrices:
don't store rotations as a matrix; instead, convert to a matrix only when
you have to, and store rotations as 'Euler paramaters'.
	Euler paramaters take up 1/4 of the room of a rotation matrix (4 floats
instead of 16), and are a lot easier to normalize.
	So, here's the code.  I know that it works; I've used it to spin a cube
at 30 frames/second for over a week, with no deformation of the cube.  The
only thing missing is the vector routines, which are boring anyway, and
are left as an exercise for the reader.  Note that this code is far from
optimized; in my case, optimization wasn't necessary, since most of the
program's time was spent either interacting with the user or drawing the
objects to be displayed.
	If you really want me to, I could dig up the reference to Euler
paramaters.  And, if you are really lazy, I could even be convinced to
send you the source to the vector routines.

/*
 * Given an axis and angle, compute euler paramaters
 * a is the axis (3 floats), phi is angle in radians, and
 * e is where to put the resulting Euler paramaters (4 floats)
 */
void
axis_to_euler(float *a, float phi, float *e)
{
	vnormal(a);	/* Normalize axis */
	vcopy(a, e);
	vscale(e, fsin(phi/2.0));
	e[3] = fcos(phi/2.0);
}
/*
 *	Given two rotations, e1 and e2, expressed as Euler paramaters,
 * figure out the equivalent single rotation and stuff it into dest.
 * 
 * This routine also normalizes the result every COUNT times it is
 * called, to keep error from creeping in.
 */
#define COUNT 100
void
add_eulers(float *e1, float *e2, float *dest)
{
	static int count=0;
	register int i;
	float t1[3], t2[3], t3[3];
	float tf[4];

	vcopy(e1, t1); vscale(t1, e2[3]);
	vcopy(e2, t2); vscale(t2, e1[3]);
	vcross(e2, e1, t3);	/* t3 = e2 cross e1 */
	vadd(t1, t2, tf);	/* tf = t1 + t2, if this was C++ */
	vadd(t3, tf, tf);	/* tf += t3 */
	tf[3] = e1[3] * e2[3] - vdot(e1, e2);

	/* I use tf[] so that e1 or e2 can be the same as dest */
	for (i=0; i < 4; i++) dest[i] = tf[i] ;
	if (++count > COUNT)
	{
		count = 0;
		normalize_euler(dest);
	}
}
/*
 * Euler paramaters always obey:  a^2 + b^2 + c^2 + d^2 = 1.0
 * We'll normalize based on this formula.  Also, normalize greatest
 * component, to avoid problems that occur when the component we're
 * normalizing gets close to zero (and the other components may add up
 * to more than 1.0 because of rounding error).
 */
void
normalize_euler(float *e)
{	/* Normalize result */
	int which, i;
	float gr;

	which = 0;
	gr = e[which];
	for (i = 1 ; i < 4 ; i++)
	{
		if (fabs(e[i]) > fabs(gr))
		{
			gr = e[i];
			which = i;
		}
	}
	e[which] = 0.0;	/* So it doesn't affect next operation */

	e[which] = fsqrt(1.0 - (e[0]*e[0] + e[1]*e[1] +
		e[2]*e[2] + e[3]*e[3]));

	/* Check to see if we need negative square root */
	if (gr < 0.0)
		e[which] = -e[which];
}

/*
 * Build a rotation matrix, given Euler paramaters.
 */
void
build_rotmatrix(Matrix m, float *e)
{
	m[0][0] = 1 - 2.0 * (e[1] * e[1] + e[2] * e[2]);
	m[0][1] = 2.0 * (e[0] * e[1] - e[2] * e[3]);
	m[0][2] = 2.0 * (e[2] * e[0] + e[1] * e[3]);
	m[0][3] = 0.0;

	m[1][0] = 2.0 * (e[0] * e[1] + e[2] * e[3]);
	m[1][1] = 1 - 2.0 * (e[2] * e[2] + e[0] * e[0]);
	m[1][2] = 2.0 * (e[1] * e[2] - e[0] * e[3]);
	m[1][3] = 0.0;

	m[2][0] = 2.0 * (e[2] * e[0] - e[1] * e[3]);
	m[2][1] = 2.0 * (e[1] * e[2] + e[0] * e[3]);
	m[2][2] = 1 - 2.0 * (e[1] * e[1] + e[0] * e[0]);
	m[2][3] = 0.0;

	m[3][0] = 0.0;
	m[3][1] = 0.0;
	m[3][2] = 0.0;
	m[3][3] = 1.0;
}
-------------
Insert Standard Disclaimer here
--gavin    (gavin@sgi.com)