|
October
31, 2005
Implementing
the Quaternion Method of Rotation
Wow, three
spatial dimensions are at hand! I’m
still not real handy with OpenGL, but I’ve got some working code, so let me
see if I can’t lay out here in a way that at least I can understand…oh, go
ahead, follow along if you dare!
So here’s
some of the global stuff, which is similar to the 2D stuff, but with and
extra D. Extraordinary, no? Note that the global variable matRot is type-defined as a
4x4 matrix of floats…why not GLfloats?
|
#include <GL/glut.h>
#include <cstdlib>
#include <cmath>
#include <iostream>
using namespace std;
/* Set initial display-window size. */
GLsizei winWidth = 600,
winHeight = 600;
/* Set range for world coordinates. */
GLfloat xwcMin = -200.0,
xwcMax = 200.; // 225.0;
GLfloat ywcMin = -200.0,
ywcMax = 200.; //225.0;
class wcPt3D {
public:
GLfloat x, y, z;
};
typedef float Matrix4x4 [4][4];
// The
global mother matrix for rotation
Matrix4x4 matRot;
|
Now we need
some functionality. This is set up in
H/B very sensibly as a more-or-less identity matrix constructor which is then
used as a basis for premultiplying (the second function below) by various
translation and rotation matrices to build up matRot.
|
/* Construct
the 4 by 4 identity matrix. */
void matrix4x4SetIdentity
(Matrix4x4 matIdent4x4) {
GLint row, col;
for
(row = 0; row < 4; row++)
for
(col = 0; col < 4 ; col++)
matIdent4x4 [row][col] = (row == col);
}
/*
Premultiply matrix m1 times matrix m2, store result in m2. */
void matrix4x4PreMultiply
(Matrix4x4 m1, Matrix4x4 m2) {
GLint row, col;
Matrix4x4 matTemp;
for (row = 0; row < 4;
row++)
for (col = 0; col < 4 ; col++)
matTemp
[row][col] = m1 [row][0] * m2 [0][col]
+ m1
[row][1] * m2 [1][col]
+ m1 [row][2]
* m2 [2][col]
+ m1
[row][3] * m2 [3][col];
for
(row = 0; row < 4; row++)
for (col = 0; col < 4; col++)
m2 [row][col] = matTemp [row][col];
}
|
…and here
are the translation and rotation functions.
The heavy computation in the rotate3D function is done by one square
root, one sine and one cosine evaluation, which are named and then used to
build the factors for the 3x3 rotation part of matRot. Each of the diagonal terms of the matrix
require two multiplications and the others require 3.
I’m getting
about a 10 second delay between when I request a run (after compiling) and
when the picture starts coming out. Is
it doing these computations then?
|
void translate3D (GLfloat
tx, GLfloat ty, GLfloat tz) {
Matrix4x4 matTransl3D;
/* Initialize translation matrix to
identity. */
matrix4x4SetIdentity (matTransl3D);
// put translation data in 4th
column of 4x4 identity
matTransl3D [0][3] = tx;
matTransl3D [1][3] = ty;
matTransl3D [2][3] = tz;
// Concatenate translation matrix
with matRot.
matrix4x4PreMultiply (matTransl3D, matRot);
}
//////////
rotate3D ////////
void rotate3D (wcPt3D p1,
wcPt3D p2, GLfloat radianAngle) {
Matrix4x4 matQuaternionRot;
GLfloat axisVectLength =
sqrt ((p2.x - p1.x) * (p2.x - p1.x)
+ (p2.y
- p1.y) * (p2.y - p1.y)
+ (p2.z
- p1.z) * (p2.z - p1.z));
GLfloat cosA = cos(radianAngle);
GLfloat oneC = 1 - cosA;
GLfloat sinA = sin(radianAngle);
GLfloat ux = (p2.x - p1.x) / axisVectLength;
GLfloat uy = (p2.y - p1.y) /
axisVectLength;
GLfloat uz = (p2.z - p1.z) /
axisVectLength;
/* Set up translation matrix for moving p1
to origin. */
translate3D (-p1.x, -p1.y, -p1.z);
/* Initialize matQuaternionRot to identity
matrix. */
matrix4x4SetIdentity
(matQuaternionRot);
matQuaternionRot [0][0] = ux*ux*oneC +
cosA;
matQuaternionRot [0][1] = ux*uy*oneC -
uz*sinA;
matQuaternionRot [0][2] = ux*uz*oneC +
uy*sinA;
matQuaternionRot [1][0] = uy*ux*oneC +
uz*sinA;
matQuaternionRot [1][1] = uy*uy*oneC +
cosA;
matQuaternionRot [1][2] = uy*uz*oneC -
ux*sinA;
matQuaternionRot [2][0] = uz*ux*oneC -
uy*sinA;
matQuaternionRot [2][1] = uz*uy*oneC +
ux*sinA;
matQuaternionRot [2][2] = uz*uz*oneC +
cosA;
/* Combine matQuaternionRot with translation
matrix. */
matrix4x4PreMultiply (matQuaternionRot,
matRot);
/* Set up inverse matTransl3D and
concatenate with
*
product of previous two matrices.
*/
translate3D (p1.x, p1.y, p1.z);
}
|
Next I
copied the same simple init()
from a the 2D program and modified triangle()
and transformVerts3D() as shown below:
|
///////
init()
void init (void) {
// Set color of display window to white.
glClearColor (1.0, 1.0, 1.0, 0.0);
}
///////
triangle draws a triangle
void triangle (wcPt3D
*verts)
{
GLint k;
glBegin (GL_TRIANGLES);
for
(k = 0; k < 3; k++)
glVertex3f (verts [k].x, verts [k].y, verts [k].z);
glEnd ( );
}
/* Using the
composite matrix, calculate transformed coordinates. */
void transformVerts3D (GLint
nVerts, wcPt3D * verts)
{
GLint k;
GLfloat temp1, temp2;
for (k
= 0; k < nVerts; k++) {
temp1 = matRot [0][0] * verts [k].x
+ matRot [0][1] * verts [k].y
+ matRot [0][2] * verts
[k].z;
temp2 = matRot [1][0] * verts [k].x
+ matRot [1][1] * verts [k].y
+ matRot [1][2] * verts
[k].z;
verts [k].z =
matRot [2][0] * verts [k].x
+ matRot [2][1] * verts
[k].y
+ matRot [2][2] * verts
[k].z;
verts [k].x = temp1;
verts [k].y = temp2;
}
}
|
Note that
only the 3x3 submatrix of matRot is needed for the rotation.
Now for the
big millstone: displayFunc() . The H/B
text does a good thing by just barely outlining what you might want here, but
enough so that you feel compelled to get in there and write something to make
a pretty picture demonstrating what’s possible. My first attempt where I feel like I’ve got
some control over iterating rotations is shown below.
I’ve set a
constant angle of rotation theta
= .25 , set up a color scheme so I can update the color on each
iteration and created a triangle. Since
p1 is at the origin, it’s evident the axis of rotation is a unit vector
directed along the positive z
axis. Thus it’s kind of like viewing a
2D rotation, which is a simplification I found I needed to get some initial
perspective. The triangle tri is a 3D
array of wcPt3D, initiated with convenient “=” overloading constructor deal,
if I’m not completely whacky? That is
pretty neat. It’s drawn in red.
Then we
initialize matRot with that constructor-like construct of making it the
identity matrix. It’s the repeated
calls to rotate3D that multiply this by the appropriate quaternion rotation
matrix. transformVerts3D then updates the triangle
vertices by multiplying it by the new rotation matrix.
|
void displayFcn (void)
{
// Input
rotation parameters.
GLfloat theta = .25;
GLfloat red = 1.0, redInc
= -0.05,
green = 0.0, greenInc = 0.03,
blue = 0.0, blueInc = 0.04;
// Input two points to determine
the axis of rotation
wcPt3D p1 = {0.0,0.0,0.0};
wcPt3D p2 = {0.0,0.0,1.0};
GLint nVerts = 3;
// Input the 3D vertices of a
triangle
wcPt3D tri[3] = { {10.0, 10.0, 1.0},
{150.0, 0.0,
0.0},
{75.0, 200.0, 2.0}
};
// Draw the initial triangle red
glColor3f(1.0,0.0,0.0);
triangle(tri);
// Initialize matRot to identity matrix:
matrix4x4SetIdentity (matRot);
for(float i = 0; i < 6.3; i += theta) {
blue += blueInc ;
green += greenInc; red +=
redInc;
glColor3f(red,green,blue);
transformVerts3D (nVerts, tri);
// Pass rotation parameters
to rotate3D.
rotate3D(p1,p2,theta);
triangle(tri);
/*cout <<
"\nrotation matrix = \n";
printMatrix(3,3,matRot);
cout << "\n\ntriangle
= \n";
cout<<tri[0].x<<"\t"<<tri[0].y<<"\t"<<tri[0].z<<endl;
cout<<tri[1].x<<"\t"<<tri[1].y<<"\t"<<tri[1].z<<endl;
cout<<tri[2].x<<"\t"<<tri[2].y<<"\t"<<tri[2].z<<endl;
cin >> ch;*/
/* Display rotated object. */
//glFlush();
}
/* Display rotated object. */
glFlush();
}
|
A simple
main routine
|
void winReshapeFcn (GLint
newWidth, GLint newHeight) {
glMatrixMode (GL_PROJECTION);
glLoadIdentity ( );
gluOrtho2D (xwcMin, xwcMax, ywcMin,
ywcMax);
glClear (GL_COLOR_BUFFER_BIT);
}
void main (int argc, char **
argv)
{
glutInit (&argc, argv);
glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
glutInitWindowPosition (50, 50);
glutInitWindowSize (winWidth,
winHeight);
glutCreateWindow ("Geometric
Transformation Sequence");
cout << "\nHowdy, pardner.";
init();
glutDisplayFunc (displayFcn);
glutReshapeFunc (winReshapeFcn);
glutMainLoop ( );
}
|
Ah, I
see. This wasn’t what I had intended,
and it’s because I’m both computing the rotation matrix repeatedly and
multiplying the new triangle by the increased rotation—which results in a
sort of angular acceleration, capiche?
If I make a slight adjustment, simply taking the rotation function out
of the loop and so computing it only once, voila:
Now
consider what happens when I give it a slight wobble out of the implicit z-range by changing the line
wcPt3D p2 = {0.0,0.0,1.0};
in displayFcn() to
wcPt3D p2 = {0.0,0.05,1.0};.
Produces
this picture:
This
picture seems to show that, as triangles twist out of the z layer the initial triangle lies in
(0<z<2) it’s extremities are
clipped off as well as changing apparent shape due to perspective. Somehow the rotations in the southern
hemisphere are completely occluded…don’t quite get why, but they seem to
reappear as the rotation comes back
into the first octant—the lightest blue triangle has its ends shorn to make a
blue pentagon.
With just a
slight wobble
wcPt3D p2 = {0.0,0.002,1.00};
and working
the loop like so
rotate3D(p1,p2,0.25);
for(float i = 0; i
< 28.3; i += theta)
we get a
picture without any clipped triangles:
If I simply
tilt the windmill’s axis slightly back instead of forward, with
wcPt3D p2 = {0.0,-0.002,1.00}
I get this:
How about
axis vector
p2 = {0.001,-0.002,1.00};
with
rotate3D(p1,p2,0.6);
for(float i = 0; i
< 12.3; i += theta)
for which
we get
|