Thread: How to make a spaceship match rotation and translation with world view in OpenGL

Thank you Nominal Animal! That was awesome! The world view was really cool as well.

It didn't work on Mac OS X because the CLOCK_REALTIME isn't recognized on it but it worked fine on Linux. In addition to my Desktop iMac, I have a laptop with Ubuntu on it to test OpenGL stuff to make sure it will work on both environments.

The main part of the problem with what I was doing was tying the Euler angles to my Quaternion stuff and I was basically still limping on a fundamental problem with Euler angles. Instead of limping with a cedar limp, I was limping on a Titanium one. In the end, I was still limping. It was kind of a slap in the face to Quaternions what I was trying to do with my code.

You have no idea how much you made my day with your example code. I was almost going to give up on this. I'm going to be learning from your example and I will slowly implement it into my existing program so that I get the full psychological benefit of learning how it works. I owe you one big time man.

Thank you again.

Originally Posted by Progger

It didn't work on Mac OS X because the CLOCK_REALTIME isn't recognized on it

Ah, I thought Apple had fixed it already.

Could you try this one? This uses a workaround, but I don't have Mac myself to test with. It would help me make better examples in the future. If it fails, the exact error would help pinpoint the reason (unless you can fix it to work on your Mac).

Code:

#define  _POSIX_C_SOURCE 200809L
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

#if defined(__APPLE__) && defined(__MACH__)
#include <OpenGL/gl.h>
#include <GLUT/glut.h>
#include <mach/clock.h>
#include <mach/mach.h>
#elif defined(__linux__)
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#else
#error Unsupported operating system.
#endif

#define UNUSED __attribute__((unused))

typedef struct {
    GLdouble    x;  /* i component */
    GLdouble    y;  /* j component */
    GLdouble    z;  /* k component */
    GLdouble    w;  /* Scalar component */
} rotationq;

typedef struct {
    GLdouble    x;  /* Rotation axis x component */
    GLdouble    y;  /* Rotation axis y component */
    GLdouble    z;  /* Rotation axis z component */
    GLdouble    r;  /* Radians per half a second */
} vernier;

typedef enum {
    CONTROL_WORLD = 0,
    CONTROL_PITCH_POS,
    CONTROL_PITCH_NEG,
    CONTROL_ROLL_POS,
    CONTROL_ROLL_NEG,
    CONTROL_YAW_POS,
    CONTROL_YAW_NEG,
    CONTROL_THRUST,
    CONTROL_COUNT
} control;

/* Keyboard controls */
static unsigned char  control_state[CONTROL_COUNT] = { 0U };
static unsigned char  control_keyboard[CONTROL_COUNT] = { 0U };
static int            control_special[CONTROL_COUNT] = { 0U };

/* Ship properties */
static GLdouble       ship_location[3] = { 100.0, 100.0, 0.0 };
static GLdouble       ship_matrix[16] = { 0 };
static GLdouble       ship_inverse[16] = { 0 };
static rotationq      ship_orientation = { 0.0, 0.0, 0.0, 1.0 };
static vernier        ship_vernier[CONTROL_COUNT] = {{0}};
static GLdouble       ship_thrust = 100.0; /* Coordinate units per second */

#if defined(__APPLE__) && defined(__MACH__)
static clock_serv_t   system_clock;
#endif

static void define_matrix(GLdouble *const matrix, const rotationq *const rq, const GLdouble *const center)
{
    const GLdouble origin[3] = { center[0], center[1], center[2] };
    /* Column-major order:
     *      matrix[0] matrix[4] matrix[8]  matrix[12]
     *      matrix[1] matrix[5] matrix[9]  matrix[13]
     *      matrix[2] matrix[6] matrix[10] matrix[14]
     *      matrix[3] matrix[7] matrix[11] matrix[15]
    */
    matrix[0]  = 1.0 - 2.0 * (rq->y * rq->y + rq->z * rq->z);
    matrix[1]  =       2.0 * (rq->x * rq->y + rq->z * rq->w);
    matrix[2]  =       2.0 * (rq->x * rq->z - rq->y * rq->w);
    matrix[3]  = 0.0;
    matrix[4]  =       2.0 * (rq->x * rq->y - rq->z * rq->w);
    matrix[5]  = 1.0 - 2.0 * (rq->x * rq->x + rq->z * rq->z);
    matrix[6]  =       2.0 * (rq->y * rq->z + rq->x * rq->w);
    matrix[7]  = 0.0;
    matrix[8]  =       2.0 * (rq->x * rq->z + rq->y * rq->w);
    matrix[9]  =       2.0 * (rq->y * rq->z - rq->x * rq->w);
    matrix[10] = 1.0 - 2.0 * (rq->x * rq->x + rq->y * rq->y);
    matrix[11] = 0.0;
    matrix[12] = -(matrix[0] * origin[0] + matrix[4] * origin[1] + matrix[8] * origin[2]);
    matrix[13] = -(matrix[1] * origin[0] + matrix[5] * origin[1] + matrix[9] * origin[2]);
    matrix[14] = -(matrix[2] * origin[0] + matrix[6] * origin[1] + matrix[10]* origin[2]);
    matrix[15] = 1.0;
}

static void define_inverse(GLdouble *const matrix, const rotationq *const rq, const GLdouble *const center)
{
    /* Inverting the non-scalar components in the rotation quaternion inverts the rotation.
     * This also moves the origin to the center. Thus, it is the opposite of define_matrix().
    */
    matrix[0]  = 1.0 - 2.0 * (rq->y * rq->y + rq->z * rq->z);
    matrix[1]  =       2.0 * (rq->x * rq->y - rq->z * rq->w);
    matrix[2]  =       2.0 * (rq->x * rq->z + rq->y * rq->w);
    matrix[3]  = 0.0;
    matrix[4]  =       2.0 * (rq->x * rq->y + rq->z * rq->w);
    matrix[5]  = 1.0 - 2.0 * (rq->x * rq->x + rq->z * rq->z);
    matrix[6]  =       2.0 * (rq->y * rq->z - rq->x * rq->w);
    matrix[7]  = 0.0;
    matrix[8]  =       2.0 * (rq->x * rq->z - rq->y * rq->w);
    matrix[9]  =       2.0 * (rq->y * rq->z + rq->x * rq->w);
    matrix[10] = 1.0 - 2.0 * (rq->x * rq->x + rq->y * rq->y);
    matrix[11] = 0.0;
    matrix[12] = center[0];
    matrix[13] = center[1];
    matrix[14] = center[2];
    matrix[15] = 1.0;
}

static void apply_vernier(rotationq *const rq, const vernier *const v, const GLdouble duration)
{
    const double vn = sqrt(v->x * v->x + v->y*v->y + v->z*v->z);
    if (vn > 0.0) {
        const GLdouble  s = sin(v->r * duration) / vn;
        const GLdouble  c = cos(v->r * duration);
        const rotationq r1 = { s * v->x, s * v->y, s * v->z, c };
        const rotationq r2 = *rq;
        const rotationq r = { r1.w*r2.x + r1.x*r2.w + r1.y*r2.z - r1.z*r2.y,
                              r1.w*r2.y - r1.x*r2.z + r1.y*r2.w + r1.z*r2.x,
                              r1.w*r2.z + r1.x*r2.y - r1.y*r2.x + r1.z*r2.w,
                              r1.w*r2.w - r1.x*r2.x - r1.y*r2.y - r1.z*r2.z };
        const double rn = sqrt(r.x*r.x + r.y*r.y + r.z*r.z + r.w*r.w);
        if (rn > 0.0) {
            rq->x = r.x / rn;
            rq->y = r.y / rn;
            rq->z = r.z / rn;
            rq->w = r.w / rn;
        } else {
            rq->x = 0.0;
            rq->y = 0.0;
            rq->z = 0.0;
            rq->w = 1.0;
        }
    }
}

static void display_window(void)
{
    const int world_mode = control_state[CONTROL_WORLD];
    double                 duration; /* Frame duration in seconds */
#if defined(__APPLE__) && defined(__MACH__)
    static mach_timespec_t prev_update = { 0 };
    static mach_timespec_t curr_update = { 0 };
    prev_update = curr_update;
    clock_get_time(system_clock, &curr_update);
    duration = (double)(curr_update.tv_sec - prev_update.tv_sec)
             + (double)(curr_update.tv_nsec - prev_update.tv_nsec) / 1000000000.0;
#else
    static struct timespec prev_update = { 0 };
    static struct timespec curr_update = { 0 };
    prev_update = curr_update;
    clock_gettime(CLOCK_MONOTONIC, &curr_update);
    duration = (double)(curr_update.tv_sec - prev_update.tv_sec)
             + (double)(curr_update.tv_nsec - prev_update.tv_nsec) / 1000000000.0;
#endif
    if (duration < 0.0 || duration > 1.0)
        duration = 0.0;

    /* Inertialess drive. Eh. */
    if (duration > 0.0) {
        size_t i = CONTROL_COUNT;
        while (i-->0)
            if (control_state[i])
                apply_vernier(&ship_orientation, ship_vernier + i, duration);
        if (control_state[CONTROL_THRUST]) {
            const double effect = duration * ship_thrust;
            ship_location[0] -= effect * ship_matrix[2];
            ship_location[1] -= effect * ship_matrix[6];
            ship_location[2] -= effect * ship_matrix[10];
        }
    }

    define_matrix(ship_matrix, &ship_orientation, ship_location);
    define_inverse(ship_inverse, &ship_orientation, ship_location);

    if (world_mode) {
        glClearColor(0.0, 0.0, 0.125, 0.0);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        glMatrixMode(GL_MODELVIEW);
        glLoadIdentity();
    } else {
        glClearColor(0.0, 0.0, 0.0, 0.0);
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
        glMatrixMode(GL_MODELVIEW);
        glLoadMatrixd(ship_matrix);
    }

    glPushMatrix();
    glTranslatef(-100.0f, -200.0f, -1000.0f);
    glColor3f(0.5f, 0.5f, 0.5f);
    glutWireSphere(400.0, 16, 16);
    glPopMatrix();

    if (world_mode) {
        glLoadIdentity();
        glLoadMatrixd(ship_inverse);
        glBegin(GL_LINES);
        glColor3f(   1.0f,   0.0f,  0.0f);
        glVertex3f(  0.0f,   0.0f,  0.0f);
        glVertex3f(  0.0f,   0.0f, 50.0f);
        glColor3f(   0.0f,   1.0f,  0.0f);
        glVertex3f(  0.0f, -25.0f,  0.0f);
        glVertex3f(  0.0f,  50.0f,  0.0f);
        glColor3f(   0.0f,   0.0f,  1.0f);
        glVertex3f(-25.0f,   0.0f,  0.0f);
        glVertex3f( 50.0f,   0.0f,  0.0f);
        glColor3f(   1.0f,   1.0f,  1.0f);
        glVertex3f(  0.0f,   0.0f,  0.0f);
        glVertex3f(  0.0f,   0.0f,-50.0f);
        glEnd();
    }        

    glFinish();
    glutSwapBuffers();
    glutPostRedisplay();
}

static void reshape_window(int width, int height)
{
    glViewport(0, 0, width, height);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    gluPerspective(60.0f, (GLdouble)width / (GLdouble)height, 1.0, 1048576.0);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glEnable(GL_DEPTH_TEST);
    glDepthFunc(GL_LESS);
    glShadeModel(GL_FLAT);
}

static void special_down(int key, int x UNUSED, int y UNUSED)
{
    if (key) {
        size_t i = CONTROL_COUNT;
        while (i-->0)
            if (key == control_special[i])
                control_state[i] = 1;
    }            
}

static void special_up(int key, int x UNUSED, int y UNUSED)
{
    if (key) {
        size_t i = CONTROL_COUNT;
        while (i-->0)
            if (key == control_special[i])
                control_state[i] = 0;
    }
}

static void keyboard_down(unsigned char key, int x UNUSED, int y UNUSED)
{
    if (key) {
        size_t i = CONTROL_COUNT;
        while (i-->0)
            if (key == control_keyboard[i])
                control_state[i] = 1;
    }
}

static void keyboard_up(unsigned char key, int x UNUSED, int y UNUSED)
{
    if (key) {
        size_t i = CONTROL_COUNT;
        while (i-->0)
            if (key == control_keyboard[i])
                control_state[i] = 0;
    }
}

static void define_vernier(vernier *const v,
                           const double x, const double y, const double z,
                           const double degrees_per_second)
{
    v->x = x;
    v->y = y;
    v->z = z;
    v->r = degrees_per_second * 3.14159265358979323846 / 360.0;
}


int main(int argc, char *argv[])
{
    glutInit(&argc, argv);
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
    glutInitWindowSize(768, 384);
    glutInitWindowPosition(128, 128);
    glutCreateWindow("Example");

#if defined(__APPLE__) && defined(__MACH__)
    host_get_clock_service(mach_host_self(), SYSTEM_CLOCK, &system_clock);
#endif

    glutDisplayFunc(display_window);
    glutReshapeFunc(reshape_window);

    glutSpecialFunc(special_down);
    glutSpecialUpFunc(special_up);
    glutKeyboardFunc(keyboard_down);
    glutKeyboardUpFunc(keyboard_up);
    glutIgnoreKeyRepeat(1);

    control_special[CONTROL_PITCH_POS] = GLUT_KEY_UP;
    control_special[CONTROL_PITCH_NEG] = GLUT_KEY_DOWN;
    control_special[CONTROL_YAW_POS] = GLUT_KEY_RIGHT;
    control_special[CONTROL_YAW_NEG] = GLUT_KEY_LEFT;
    control_keyboard[CONTROL_ROLL_POS] = 'z';
    control_keyboard[CONTROL_ROLL_NEG] = 'x';
    control_keyboard[CONTROL_THRUST] = ' ';
    control_keyboard[CONTROL_WORLD] = '\r';

    fprintf(stdout, "Controls:\n");
    fprintf(stdout, "\tUp, Down    - Pitch\n");
    fprintf(stdout, "\tLeft, Right - Yaw\n");
    fprintf(stdout, "\tZ, X        - Roll\n");
    fprintf(stdout, "\tSpace       - Thrust\n");
    fprintf(stdout, "\tEnter       - World view (blue)\n");
    fflush(stdout);

    define_vernier(ship_vernier + CONTROL_PITCH_POS, 1.0, 0.0, 0.0,  +90.0);
    define_vernier(ship_vernier + CONTROL_PITCH_NEG, 1.0, 0.0, 0.0,  -90.0);
    define_vernier(ship_vernier + CONTROL_YAW_POS,   0.0, 1.0, 0.0,  +90.0);
    define_vernier(ship_vernier + CONTROL_YAW_NEG,   0.0, 1.0, 0.0,  -90.0);
    define_vernier(ship_vernier + CONTROL_ROLL_POS,  0.0, 0.0, 1.0,  -90.0);
    define_vernier(ship_vernier + CONTROL_ROLL_NEG,  0.0, 0.0, 1.0,  +90.0);
    define_vernier(ship_vernier + CONTROL_THRUST,    0.0, 0.0, 0.0,    0.0);
    define_vernier(ship_vernier + CONTROL_WORLD,     0.0, 0.0, 0.0,    0.0);

    glutMainLoop();

#if defined(__APPLE__) && defined(__MACH__)
    mach_port_deallocate(mach_task_self(), system_clock );
#endif

    return EXIT_SUCCESS;
}

Originally Posted by Progger

a fundamental problem with Euler angles.

Yeah, I remember my astonishment years ago; I had tried hours and hours on end to get matrix normalization working (the best did it iteratively in a small fixed number of steps), and then having seen how simply it works with versors aka unit quaternions aka rotation quaternions.. well, it blew my mind.

As to the quaternion-to-matrix calculation in define_matrix(): It is the same as the latter form in section Quaternion-derived rotation matrix in the Quaternions and spatial rotation Wikipedia article, but with column-major order, as the comment shows. The translation vector is (matrix[12], matrix[13], matrix[14]), but because the translation occurs before the rotation, the center vector is multiplied by the rotation matrix and negated, to get the translation vector. (The result is the same if you were to subtract the center vector from all points before rotation.)

The inverse one relies on sin() being odd and cos() even. Negating the rotation angle, which inverts the rotation part of the matrix (but is more precise than calculating the matrix inverse, numerically -- mathematically both give the exact same answer), has the same effect as negating the vector components of the rotation quaternion. That in turn is just a few sign changes in how the matrix is calculated, simple!
Of course, in the reverse case the translation to center occurs after the rotation, so in this case the center vector is not multiplied by the rotation matrix.

Originally Posted by Progger

I'm going to be learning from your example and I will slowly implement it into my existing program so that I get the full psychological benefit of learning how it works.

Excellent: that is the best I could ask for. (I am here, after all, to help others learn, and learn myself too.) If you get stuck on some detail at some point, do let me know, and I'll try to clarify.

Originally Posted by Nominal Animal

Ah, I thought Apple had fixed it already.

Could you try this one? This uses a workaround, but I don't have Mac myself to test with. It would help me make better examples in the future. If it fails, the exact error would help pinpoint the reason (unless you can fix it to work on your Mac).

This one compiled and ran perfectly. Just to be on the safe side, I went through every rotation forwards and backwards in both first person view and world view and it's perfect.

I came across a similar situation when I saw example code for a custom DrawText function using glutBitmapString (a freeglut function and not the original glut) and my compiler kept saying it was undefined. I ended up having to make a modified version using glutBitmapCharacter with a char array in order for it to work on my Mac. My Linux uses the freeglut library but my Mac uses the original stuff. It's kind of nice to have both environments to test my programs on so that I know it will be able to run on almost any machine with Linux and Mac OS X. Once I add music and sound to whatever game I eventually make, that's when they will have to part ways and one version will be Objective-C++ so that I can use Mac's Cocoa library for that stuff. And I will have to figure out the API's on Linux for the Linux version.

Excellent: that is the best I could ask for. (I am here, after all, to help others learn, and learn myself too.) If you get stuck on some detail at some point, do let me know, and I'll try to clarify.

I definitely want to learn. C/C++/OpenGL is my hobby now. I watched lectures on YouTube about Quaternions and Axis-Angle rotations. I was on Khan Academy several times learning about vector dot product and vector cross product, matrix multiplication and trig. I know there are guys on here that are apprehensive about sharing code with others because they don't want someone to just copy their program verbatim and improve on it and pretend they understand how it works. In all candor, I'm still studying the printout and slowly absorbing what is going on with the example program you gave me. It will take me a few days to understand it but at least I have a working Quaternion example program to learn from.

Thank you.

Originally Posted by Progger

I know there are guys on here that are apprehensive about sharing code with others because they don't want someone to just copy their program verbatim and improve on it and pretend they understand how it works.

Well, yes and no. Yes, showing complete code is not usually done. No, it's not apprehension about sharing, it's just that showing complete code is usually counterproductive: it does more harm than good.

You did the right thing from the get go by showing your own efforts first. Your questions were to the point, and showed that you understand what the code does (and thus is not just something you copied verbatim from somewhere -- even if you had copied it from somewhere else, you certainly had spent the effort to understand how it works).

The changes I would suggest were too extensive to make it useful, and that's why I wrote a counter-example from scratch. There are three key differences to yours. First is the obvious, the way the OpenGL transformations are done; it's not that much code, but understanding it will take a bit of time. The second one is how my code retains a state of keys being pressed at any point of time. Your code acts upon the key whenever it gets pressed, mine applies the change the key is in charge of, whenever a new frame is drawn. The third one is that my code measures the duration of the last frame, and uses that as the time step (for both movement and rotation).

Originally Posted by Progger

I'm still studying the printout and slowly absorbing what is going on with the example program you gave me.

Play with it. Change stuff, take it apart.

I now realize that I should have included the link to the Wikipedia articles on Quaternions and spatial rotation and Rotation matrix in the comments.

It is also important (not just for the quaternion stuff in the example, but for OpenGL stuff in general) to understand that the 4x4 transformation matrices use the following logic,

Xx	Xy	Xz	Tx		px		Vx
Yx	Yy	Yz	Ty	.	py	=	Vy
Zx	Zy	Zz	Tz		pz		Vz
0	0	0	1		1		1

or

V_x = X_x p_x + X_y p_y + X_z p_z + T_x
V_y = Y_x p_x + Y_y p_y + Y_z p_z + T_y
V_z = Z_x p_x + Z_y p_y + Z_z p_z + T_z

Using the original coordinates, the X-axis unit vector in the rotated coordinate system is (X_x, X_y, X_z); similarly for the Y and Z-axis unit vectors. The vector (T_x, T_y, T_z) specifies the translation after rotation, where the origin in the rotated coordinates will be using the old coordinate system.

The order OpenGL expects to see these in memory (in an 1D array) is X_x Y_x Z_x 0 X_y Y_y Z_y 0 X_z Y_z Z_z 0 T_x T_y T_z 1.

If you need translation B before the rotation, it is calculated by applying the rotation part:

T_x = X_x B_x + Y_x B_y + Z_x B_z
T_y = X_y B_x + Y_y B_y + Z_y B_z
T_z = X_z B_x + Y_z B_y + Z_z B_z

Various translations are simply added together. (So, if you want to rotate around a specific point B, without moving the origin, just negate the above, and add B. That way B is moved to origin prior to rotation, and back to where it was relative to origin after the rotation.)

Since the ship view is towards the negative Z axis, using the original coordinates, the ship points towards (Z_x, Z_y, Z_z). (In the program, these are ship_matrix[2], ship_matrix[6], ship_matrix[10].)

Finally, the above reflects my current understanding. If something sounds wrong, verify it: my explanation above might just be transposed. (I do not think these things in words, but directly as operations; I find rotating and transforming 3D views and objects visually in my mind easier than describing them in words. Thus, my words can easily be wrong, even if the thing I am trying to describe is correct. Although I'm pretty often wrong, too.) So, don't take any of this as fact; verify everything until your expectations based on your understanding matches the effects you get.

Originally Posted by Nominal Animal

Well, yes and no. Yes, showing complete code is not usually done. No, it's not apprehension about sharing, it's just that showing complete code is usually counterproductive: it does more harm than good.

You're right. I can appreciate that. I was really desperate in this particular case because unless I went and enrolled in a 3D game programming course or 3D graphics course at a college or something, I never would have been able to figure this out without the help of fine people like yourself. What is explained in mathematics circles on the websites talking about Quaternions, and what has to be coded into a GL to make a simulation work is tricky to figure out without professional help.

You did the right thing from the get go by showing your own efforts first. Your questions were to the point, and showed that you understand what the code does (and thus is not just something you copied verbatim from somewhere -- even if you had copied it from somewhere else, you certainly had spent the effort to understand how it works).

I kind of had an itch to learn this stuff for a long time. I started learning computer programming in college back in 1991 to 1994 and back then, they only taught C, COBOL, REXX, JCL and mainframe assembly language. I learned Java a few years later at another college night school course and concentrated on that since networking, music, sound and animation are easy to do with it, but I always had a soft-spot for C and learning C++ and OpenGL later on.

The changes I would suggest were too extensive to make it useful, and that's why I wrote a counter-example from scratch. There are three key differences to yours. First is the obvious, the way the OpenGL transformations are done; it's not that much code, but understanding it will take a bit of time. The second one is how my code retains a state of keys being pressed at any point of time. Your code acts upon the key whenever it gets pressed, mine applies the change the key is in charge of, whenever a new frame is drawn. The third one is that my code measures the duration of the last frame, and uses that as the time step (for both movement and rotation).

Yeah, I kept my interface very crude and rudimentary for the time-being because I wanted to focus on getting a working sim of rotating on every axis and upside down in a simulated 3D space because if that works, I can take my time figuring out everything else. It will take me a bit more time to understand the other keyboard functions you implemented but I was impressed with how it worked in your example.

Play with it. Change stuff, take it apart.

I basically focused on the functions and variables that were required to help me understand how and in what order the quaternions had to be multiplied together and I changed my program to adapt to what I learned from your example with regard to quaternions and the formulas for creating the final matrix. The one I was using before wasn't accurate. I got it working now thanks to you. I added a chase cam view and a world view to my sim as well to help me to visualize better how objects are moving and appear from alternate viewpoints. I was having trouble getting the spaceship to match the viewpoint with the translation information mixed in with the rotation information so I kept two separate matrices for the rotation and translation so that I could invert each and apply them in opposite order for the object that my view is attached to. I'm sure there was probably another way to do it but this was the easiest way that I knew how to figure it out for now.

No offence to memory referencing arithmetic and pointers but I kept those out of my program (other than for passing arrays to functions which are always passed by reference rather than passing a second copy of all the values), because it improves the readability of the code for me and keeps it easier to follow. And it will be easier for me to port to Java or some other programming language in the future that may not allow arithmetic on memory addresses or pointers. I know that from a performance standpoint, pointers are more efficient but for me, that's a bit of a last resort.

Let me know if you want me to post the code for my updated program. It's obviously not great as far as the input method is concerned for the controls but I will eventually want to use specialized mouse and keyboard events so that it's more professionally controlled.

Hi Nominal Animal.

I don't know if you are still checking this forum out but I seemed to have hit another snag. Any help you could give me would be great.

The good news is that I was able to get a drone ship to face any target I wanted to (even me if I want it to track me). I cheated a bit because I do a gluLookAt from the ship to the target vector and then after copying the matrix, inverting it and using it to re-orient the object. This is fine for accuracy, but bad for slowly rotating it rather than just instantly snapping to its target.

When I researched converting a rotation matrix to a quaternion, I created the function to do it and it matches two sites I checked out. Unfortunately when I do the quaternion multiplication, it doesn't match up with what the inverted gluLookAt rotation matrix does. I want to get this to work because then I can use the formula to interpolate between two quaternions (which seems easy to do) which is better than trying to do it between two rotation matrices.

Did you ever get a rotation matrix to quaternion conversion to work perfectly? Otherwise, did you ever have to make another object track you while moving around in 3D?

If I can get this rotation matrix to quaternion conversion to work perfectly, it would be a major hurdle to put behind me once and for all.

Here is the code for the conversion that I have been using:

Code:

void Ship::rotationMatrixToQuaternion(GLfloat inp_mat[16])
{
    GLfloat t = 1 + inp_mat[0] + inp_mat[5] + inp_mat[10];
    GLfloat s = 0.0;
    
    if (t > 0.0)
    {
        s = sqrt(t) * 2;
        rmq_x = (inp_mat[9] - inp_mat[6]) / s;
        rmq_y = (inp_mat[2] - inp_mat[8]) / s;
        rmq_z = (inp_mat[4] - inp_mat[1]) / s;
        rmq_w = 0.25 * s;
    }
    else if ((inp_mat[0] > inp_mat[5]) && (inp_mat[0] > inp_mat[10]))
    {
        s = sqrt(1.0 + inp_mat[0] - inp_mat[5] - inp_mat[10]) * 2;
        rmq_x = 0.25 * s;
        rmq_y = (inp_mat[4] + inp_mat[1] ) / s;
        rmq_z = (inp_mat[2] + inp_mat[8] ) / s;
        rmq_w = (inp_mat[9] - inp_mat[6] ) / s;
    }
    else if (inp_mat[5] > inp_mat[10])
    {
        s = sqrt(1.0 + inp_mat[5] - inp_mat[0] - inp_mat[10]) * 2;
        rmq_x = (inp_mat[4] + inp_mat[1]) / s;
        rmq_y = 0.25 * s;
        rmq_z = (inp_mat[9] + inp_mat[6]) / s;
        rmq_w = (inp_mat[2] - inp_mat[8]) / s;
    }
    else
    {
        s = sqrt(1.0 + inp_mat[10] - inp_mat[0] - inp_mat[5]) * 2;
        rmq_x = (inp_mat[2] + inp_mat[8]) / s;
        rmq_y = (inp_mat[9] + inp_mat[6]) / s;
        rmq_z = 0.25 * s;
        rmq_w = (inp_mat[4] - inp_mat[1]) / s;
    }
}

rmq stands for "rotation matrix quaternion" and the 4 variables are global variables. This is OpenGL I'm using so the matrix is in column major order. [0] to [4] is the first column on the left from top to bottom, followed by [5] to [8] for the second column from top to bottom, so on and so forth.

I would appreciate any help you could give me Nominal Animal. This is driving me crazy.

Originally Posted by Progger

I don't know if you are still checking

I took a little hiatus. Probably will take a longer one.

Originally Posted by Progger

When I researched converting a rotation matrix to a quaternion, I created the function to do it and it matches two sites I checked out.

Don't trust the web. Unless the round-trip conversions produce the original quaternion -- say, dot product between the two quaternions within 0.9999999 and 1.0000001 --, one of the conversions is wrong.

Consider the following test bench I whipped up:

Code:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define   EPSILON       0.000000001
#define   SQR_EPSILON  (EPSILON*EPSILON)

#if defined(__APPLE__) && defined(__MACH__)
#include <OpenGL/gl.h>
#include <GLUT/glut.h>
#include <mach/clock.h>
#include <mach/mach.h>
#elif defined(__linux__)
#include <GL/gl.h>
#include <GL/glu.h>
#include <GL/glut.h>
#else
#error Unsupported operating system.
#endif

typedef  double  vec4d __attribute__((vector_size (4 * sizeof (double))));

typedef struct {
    union {
        vec4d           c;
        struct {
            double      x;
            double      y;
            double      z;
            double      w;
        };
    };
} rotation;

/* For OpenGL, you can supply transform.m as the 16-element array.
*/
typedef struct {
    union {
        vec4d           c[4];
        double          m[16];
        struct {
            vec4d       x;
            vec4d       y;
            vec4d       z;
            vec4d       t;
        };
    };
} transform;

/* Normalize a rotation.
 * Use before rotation_to_transform().
*/
static void rotation_normalize(rotation *const r)
{
    if (r) {
        const vec4d rr = r->c * r->c;
        const double rsqr = rr[0] + rr[1] + rr[2] + rr[3];
        if (rsqr <= SQR_EPSILON) {
            r->c[0] = 0.0;
            r->c[1] = 0.0;
            r->c[2] = 0.0;
            r->c[3] = 1.0;
        } else
        if (rsqr != 1.0) {
            const double scale = (r->c[3] >= 0.0) ? 1.0 / sqrt(rsqr) : -1.0 / sqrt(rsqr);
            const vec4d rs = { scale, scale, scale, scale };
            r->c *= rs;
        }
    }
}

/* Populate rotation part of transform t with rotation r.
*/
static void rotation_to_transform(transform *const t, const rotation *const r)
{
    if (t) {
        if (r) {
            const vec4d rx = { r->c[0], r->c[0], r->c[0], r->c[0] };
            const vec4d ry = { r->c[1], r->c[1], r->c[1], r->c[1] };
            const vec4d rz = { r->c[2], r->c[2], r->c[2], r->c[2] };
            const vec4d rc = r->c + r->c;
            const vec4d x  = rc * rx;
            const vec4d y  = rc * ry;
            const vec4d z  = rc * rz;
            t->c[0][0] = 1.0 - y[1] - z[2]; /* 1 - 2 * r->y * r->y - 2 * r->z * r->z */
            t->c[0][1] =       x[1] - z[3]; /*     2 * r->x * r->y - 2 * r->z * r->w */
            t->c[0][2] =       x[2] + y[3]; /*     2 * r->x * r->z + 2 * r->y * r->w */
            t->c[1][0] =       x[1] + z[3]; /*     2 * r->x * r->y + 2 * r->y * r->w */
            t->c[1][1] = 1.0 - x[0] - z[2]; /* 1 - 2 * r->x * r->x - 2 * r->z * r->z */
            t->c[1][2] =       y[2] - x[3]; /*     2 * r->y * r->z - 2 * r->x * r->w */
            t->c[2][0] =       x[2] - y[3]; /*     2 * r->x * r->z - 2 * r->y * r->w */
            t->c[2][1] =       y[2] + x[3]; /*     2 * r->y * r->z + 2 * r->x * r->w */
            t->c[2][2] = 1.0 - x[0] - y[1]; /* 1 - 2 * r->x * r->x - 2 * r->y * r->y */
        } else {
            t->c[0][0] = 1.0;
            t->c[0][1] = 0.0;
            t->c[0][2] = 0.0;
            t->c[1][0] = 0.0;
            t->c[1][1] = 1.0;
            t->c[1][2] = 0.0;
            t->c[2][0] = 0.0;
            t->c[2][1] = 0.0;
            t->c[2][2] = 0.0;
        }
    }
}

/* Calculate rotation r from transformation t.
*/
static void transform_to_rotation(rotation *const r, const transform *const t)
{
    if (r) {
        if (t) {
            const double det = t->c[0][0] * ( t->c[1][1] * t->c[2][2] - t->c[1][2] * t->c[2][1] )
                             + t->c[0][1] * ( t->c[1][2] * t->c[2][0] - t->c[1][0] * t->c[2][2] )
                             + t->c[0][2] * ( t->c[1][0] * t->c[2][1] - t->c[1][1] * t->c[2][0] );
            const double one = cbrt(det);
            const double xsqr = one + t->c[0][0] - t->c[1][1] - t->c[2][2];
            const double ysqr = one - t->c[0][0] + t->c[1][1] - t->c[2][2];
            const double zsqr = one - t->c[0][0] - t->c[1][1] + t->c[2][2];
            const double wsqr = one + t->c[0][0] + t->c[1][1] + t->c[2][2];

            if (t->c[2][1] > t->c[1][2])
                r->c[0] = (xsqr > 0.0) ?  0.5 * sqrt(xsqr) :  0.0;
            else
            if (t->c[2][1] < t->c[1][2])
                r->c[0] = (xsqr > 0.0) ? -0.5 * sqrt(xsqr) : -0.0;
            else
                r->c[0] = 0.0;

            if (t->c[0][2] > t->c[2][0])
                r->c[1] = (ysqr > 0.0) ?  0.5 * sqrt(ysqr) :  0.0;
            else
            if (t->c[0][2] < t->c[2][0])
                r->c[1] = (ysqr > 0.0) ? -0.5 * sqrt(ysqr) : -0.0;
            else
                r->c[1] = 0.0;

            if (t->c[1][0] > t->c[0][1])
                r->c[2] = (zsqr > 0.0) ?  0.5 * sqrt(zsqr) :  0.0;
            else
            if (t->c[1][0] < t->c[0][1])
                r->c[2] = (zsqr > 0.0) ? -0.5 * sqrt(zsqr) : -0.0;
            else
                r->c[2] = 0.0;

            r->c[3] = (wsqr > 0.0) ? 0.5 * sqrt(wsqr) : 0.0;

        } else {
            r->c[0] = 0.0;
            r->c[1] = 0.0;
            r->c[2] = 0.0;
            r->c[3] = 1.0;
        }
    }
}

static double rand_one(void)
{
    return (rand() / (double)RAND_MAX) * 2.0 - 1.0;
}

#define MAX_ERROR 1.0e-9
#define TESTS     1000000

int main(void)
{
    rotation from, to;
    transform t;
    int i;

    for (i = 1; i <= TESTS; i++) {

        /* Random unit quaternion: */
        from.x = rand_one();
        from.y = rand_one();
        from.z = rand_one();
        from.w = rand_one();
        rotation_normalize(&from);
        if (fabs(from.c[0]*from.c[0] + from.c[1]*from.c[1] + from.c[2]*from.c[2] + from.c[3]*from.c[3] - 1.0) > MAX_ERROR) {
            fprintf(stderr, "Error in rotation_normalize(): rotation quaternion is not unit length.\n");
            return EXIT_FAILURE;
        }

        /* Construct the transformation matrix. */
        rotation_to_transform(&t, &from);

        /* Verify transformation matrix. */
        if (fabs(t.c[0][0] * t.c[0][0] + t.c[0][1] * t.c[0][1] + t.c[0][2] * t.c[0][2] - 1.0) > MAX_ERROR ||
            fabs(t.c[1][0] * t.c[1][0] + t.c[1][1] * t.c[1][1] + t.c[1][2] * t.c[1][2] - 1.0) > MAX_ERROR ||
            fabs(t.c[2][0] * t.c[2][0] + t.c[2][1] * t.c[2][1] + t.c[2][2] * t.c[2][2] - 1.0) > MAX_ERROR ||
            fabs(t.c[0][0] * t.c[0][0] + t.c[1][0] * t.c[1][0] + t.c[2][0] * t.c[2][0] - 1.0) > MAX_ERROR ||
            fabs(t.c[0][1] * t.c[0][1] + t.c[1][1] * t.c[1][1] + t.c[2][1] * t.c[2][1] - 1.0) > MAX_ERROR ||
            fabs(t.c[0][2] * t.c[0][2] + t.c[1][2] * t.c[1][2] + t.c[2][2] * t.c[2][2] - 1.0) > MAX_ERROR ||
            fabs(t.c[0][0] * t.c[1][0] + t.c[0][1] * t.c[1][1] + t.c[0][2] * t.c[1][2])       > MAX_ERROR ||
            fabs(t.c[0][0] * t.c[2][0] + t.c[0][1] * t.c[2][1] + t.c[0][2] * t.c[2][2])       > MAX_ERROR ||
            fabs(t.c[1][0] * t.c[2][0] + t.c[1][1] * t.c[2][1] + t.c[1][2] * t.c[2][2])       > MAX_ERROR) {
            fprintf(stderr, "Error in rotation_to_transform(): matrix is not orthonormal.\n");
            return EXIT_FAILURE;
        }

        /* Reconstruct the rotation quaternion. */
        transform_to_rotation(&to, &t);

        /* Verify rotation quaternion length. */
        if (fabs(to.c[0]*to.c[0] + to.c[1]*to.c[1] + to.c[2]*to.c[2] + to.c[3]*to.c[3] - 1.0) > MAX_ERROR) {
            fprintf(stderr, "Error in transform_to_rotation(): rotation quaternion is not unit length.\n");
            return EXIT_FAILURE;
        }

        /* Verify quaternions match (inner product is 1). */
        if (fabs(from.c[0]*to.c[0] + from.c[1]*to.c[1] + from.c[2]*to.c[2] + from.c[3]*to.c[3] - 1.0) > MAX_ERROR) {
            fprintf(stderr, "Error in transform_to_rotation(): Quaternion differs from original.\n");
            return EXIT_FAILURE;
        }

        fprintf(stdout, "\r%d of %d OK     \r", i, TESTS);
        fflush(stdout);
    }

    return EXIT_SUCCESS;
}

If you run the program, it generates a million test unit quaternions, verifies that each quaternion can be converted to matrix notation (i.e. that the generated matrix is orthonormal) using rotation_to_transform(), and that transform_to_rotation() generates the original quaternion (with nine significant digits or so).
I'd prefer to use a better pseudo-random number generator than rand(), but I thought that alone should suffice for quick testing.

I'd like to explain the exact math behind this, but if I do, someone like Hodor will eventually complain how that is waffling beyond the purpose of this discussion board, and that I'm surely not seriously expecting that anyone cares. If you want to discuss that further, my e-mail address is on my home page.

The code uses the GCC vector notation, in the hopes that whatever architecture it is compiled on, GCC can vectorize some of the operations. The transform type (use the m member) is compatible with OpenGL transformation matrices.

The above functions only touch the rotation sub-matrix, so for OpenGL, you need to ensure the fourth components are 0, and set the translation (and fourth component of translation to 1). For example, you can set them:

Code:

    transform.m[3] = 0.0;
    transform.m[7] = 0.0;
    transform.m[11] = 0.0;
    transform.m[15] = 1.0;
    transform.c[3][0] = 0.0; /* X translation */
    transform.c[3][1] = 0.0; /* Y translation */
    transform.c[3][2] = 0.0; /* Z translation */

The transform structure has the unit vectors in .c[0], .c[1], and .c[2], and should have their fourth components zeroed. The translation vector (after rotation) is .c[3], and its fourth component should be 1.

Originally Posted by Progger

I want to get this to work because then I can use the formula to interpolate between two quaternions (which seems easy to do) which is better than trying to do it between two rotation matrices.

Yes.

Above, I've taken the liberty of normalizing the quaternions with the scalar component nonnegative. This essentially means that rotation_normalize() ensures the rotation is at most 180 degrees. Larger rotations simply have the axis and angle negated. (Remember, rotation quaternions can be negated without it affecting the actual rotation (or matrix) at all.)
I don't think that induces any issues in practice, and it should ensure large turns use the shorter isocircle segment, but if I am wrong, the issues should be easily visible when blending over a large difference.

The blend itself really is trivial; you just interpolate between the two components, and normalize the resulting quaternion:

Code:

static void rotation_blend(rotation *const r, const rotation *const r0, const rotation *const r1, const double phase)
{
    if (r && r0 && r1) {
        if (phase <= 0.0)
            r->c = r0->c;
        else
        if (phase >= 1.0)
            r->c = r1->c;
        else {
            const vec4d p0 = { 1.0 - phase, 1.0 - phase, 1.0 - phase, 1.0 - phase };
            const vec4d p1 = {       phase,       phase,       phase,       phase };
            r->c = r0->c * p0 + r1->c * p1;
        }
    } else
    if (r && r0)
        r->c = r0->c;
    else
    if (r && r1)
        r->c = r1->c;
    else
    if (r) {
        r->c[0] = 0.0;
        r->c[1] = 0.0;
        r->c[2] = 0.0;
        r->c[3] = 1.0;
    }

    /* Finish with rotation_normalize(r): */
    if (r) {
        const vec4d rr = r->c * r->c;
        const double rsqr = rr[0] + rr[1] + rr[2] + rr[3];
        if (rsqr <= SQR_EPSILON) {
            r->c[0] = 0.0;
            r->c[1] = 0.0;
            r->c[2] = 0.0;
            r->c[3] = 1.0;
        } else
        if (rsqr != 1.0) {
            const double scale = (r->c[3] >= 0.0) ? 1.0 / sqrt(rsqr) : -1.0 / sqrt(rsqr);
            const vec4d rs = { scale, scale, scale, scale };
            r->c *= rs;
        }
    }
}

Above, phase is 0 for the first rotation, 1 for the latter. If you want a smooth transition (so that the "camera" movest fastest in the middle, and starts and stops smoothly), use

3*phase*phase - 2*phase*phase

as the fourth parameter to the function instead.

It's over-carefully written, too: if you know you'll always supply valid parameters to it, you can shave it quite a bit smaller. (Applies to typically all my code, really. It seems I care about error cases much more than the average programmer.)

Originally Posted by Progger

Did you ever get a rotation matrix to quaternion conversion to work perfectly?

Yes, see above. The test program even lets you verify it works.

Inverting an orthogonal matrix is trivial: it's the same as transpose. You only need to swap twelve elements in the transformation matrix. For pure rotation, as in the functions above, the transformation matrix is orthonormal.

In my earlier example, I do believe I used the inverse Z vector (that is, the Z vector of the transposed transformation matrix) to show the orientation of the ship in the world view.

Originally Posted by Progger

did you ever have to make another object track you while moving around in 3D?

Not really.

I don't see any real issues with that, really. If it is an "enemy", I'd write a simple AI to control it. Just having it smoothly follow something looks robotic, and might require non-Newtonian mechanics; using a bit of logic to control that like the user controls theirs makes things much more interesting. (For example, humans often prefer a "horizon", with their enemies either to the side, or to "above" or "below". This may dictate attitude thruster positions, and might mean that forcing the AI pilot to first rotate before turning makes it more human-like.)

Originally Posted by Nominal Animal

I took a little hiatus. Probably will take a longer one.


Don't trust the web. Unless the round-trip conversions produce the original quaternion -- say, dot product between the two quaternions within 0.9999999 and 1.0000001 --, one of the conversions is wrong.

Consider the following test bench I whipped up:

Wow, thank you! I wasn't expecting all that code.

My intuition was telling me that something wasn't adding up. I didn't know how I was going to solve the problem but judging from how the quaternion information was plugged into a rotation matrix, I knew that doing the reverse was going to go beyond basic algebra. It seemed kind of like one of those scenarios where it's probably easier to bang three protons out of a lead atom to make gold than it is to insert three protons into a gold atom to make lead. If the inclination was there of course.

I won't be able to experiment with the code for a couple of days because business started getting a little busy again after the holidays but I do appreciate it. Thank you very much.

I don't see any real issues with that, really. If it is an "enemy", I'd write a simple AI to control it. Just having it smoothly follow something looks robotic, and might require non-Newtonian mechanics; using a bit of logic to control that like the user controls theirs makes things much more interesting. (For example, humans often prefer a "horizon", with their enemies either to the side, or to "above" or "below". This may dictate attitude thruster positions, and might mean that forcing the AI pilot to first rotate before turning makes it more human-like.)

Oh yeah, the AI part I know exactly what you mean. I used to think that I never took human reasoning and intelligence for granted until I had to try and write code out of it.

I was just wondering if you ever had to use code to rotate to track a user moving around in 3D. All this AI stuff is going to be more artsy fartsy logic than purely an absolute results-based problem. In the back of my mind, I've been breaking down all sorts of ideas involving formation flying, interpreting user movements to adapt to a different strategy, etc.

My brother was telling me that someone tried to write code that was smart enough to know when a commercial was on the screen, rather than a movie you were watching and it's impossible to teach a computer to be smart enough to recognize a commercial it hasn't seen before rather than a movie you were watching. A human can figure it out no problem as we all know. I supposed the big commercial entities won't ever allow mandatory coding in the signal so that software can tell the difference either. It's hard to fund movies and television shows that way by ticking off your sponsors so we can only dream.

Anyway, thank you again Nominal Animal. I appreciate it.