python - OpenGL Rotation from Velocity Vector

Question

This should be easy, but I've been all over trying to find a simple explanation that I can grasp. I have an object that I'd like to represent in OpenGL as a cone. The object has x, y, z coordinates and a velocity vector vx, vy, and vz. The cone should point in the direction of the velocity vector.

So, I think my PyOpenGL code should look something like this:

glPushMatrix()
glTranslate(x, y, z)
glPushMatrix()

# do some sort of rotation here #

glutSolidCone(base, height, slices, stacks)
glPopMatrix()
glPopMatrix()

So, is that correct (so far)? What do I put in place of the "# do some sort of rotation here #" ?

In my world, the Z-axis points up (0, 0, 1) and, without any rotations, so does my cone.

---

Okay, Reto Koradi's answer seems to be the approach that I should take, but I'm not sure of some of the implementation details and my code is not working.

If I understand correctly, the rotation matrix should be a 4x4. Reto shows me how to get a 3x3, so I'm assuming that the 3x3 should be the upper-left corner of a 4x4 identity matrix. Here's my code:

import numpy as np

def normalize(v):
    norm = np.linalg.norm(v)
    if norm > 1.0e-8:  # arbitrarily small
        return v/norm
    else:
        return v

def transform(v):
    bz = normalize(v)
    if (abs(v[2]) < abs(v[0])) and (abs(v[2]) < abs(v[1])):
        by = normalize(np.array([v[1], -v[0], 0]))
    else:
        by = normalize(np.array([v[2], 0, -v[0]]))
        #~ by = normalize(np.array([0, v[2], -v[1]]))

    bx = np.cross(by, bz)
    R =  np.array([[bx[0], by[0], bz[0], 0],
                   [bx[1], by[1], bz[1], 0],
                   [bx[2], by[2], bz[2], 0],
                   [0,     0,     0,     1]], dtype=np.float32)

    return R

and here is the way it gets inserted into the rendering code:

glPushMatrix()
glTranslate(x, y, z)
glPushMatrix()

v = np.array([vx, vy, vz])
glMultMatrixf(transform(v))

glutSolidCone(base, height, slices, stacks)
glPopMatrix()
glPopMatrix()

Unfortunately, this isn't working. My test case cones just do not point correctly and I can't identify the failure mode. Without the "glutMultMatrixf(transform(v)" line, the cones align along the z-axis, as expected.

---

It's working. Reto Koradi correctly identified that the rotation matrix needed to be transposed in order to match the column order of OpenGL. The code should look like this (before optimization):

def transform(v):
    bz = normalize(v)
    if (abs(v[2]) < abs(v[0])) and (abs(v[2]) < abs(v[1])):
        by = normalize(np.array([v[1], -v[0], 0]))
    else:
        by = normalize(np.array([v[2], 0, -v[0]]))
        #~ by = normalize(np.array([0, v[2], -v[1]]))

    bx = np.cross(by, bz)
    R =  np.array([[bx[0], by[0], bz[0], 0],
                   [bx[1], by[1], bz[1], 0],
                   [bx[2], by[2], bz[2], 0],
                   [0,     0,     0,     1]], dtype=np.float32)

    return R.T

Answer 1

A helpful concept to remember here is that a linear transformation can also be interpreted as a change of coordinate systems. In other words, instead of picturing points being transformed within a coordinate system, you can just as well picture the points staying in place, but their coordinates being expressed in a new coordinate system. When looking at the matrix expressing the transformation, the base vectors of this new coordinate system are the column vectors of the matrix.

In the following, the base vectors of the new coordinate system are named bx, by and bz. Since the columns of a rotation matrix need to be orthonormal, bx, by and bz need to form an orthonormal set of vectors.

In this case, the original cone is oriented along the z-axis. Since you want the cone to be oriented along (vx, vy, vz) instead, we use this vector as the z-axis of our new coordinate system. Since we want an orthonormal coordinate system, the only thing left to do to obtain bz is to normalize this vector:

               [vx]
bz = normalize([vy])
               [vz]

Since the cone is rotationally symmetrical, it does not really matter how the remaining two base vectors are chosen, just as long as they are both orthogonal to bz, and orthogonal to each other. A simple way to find an arbitrary orthogonal vector to a given vector is to keep one coordinate 0, swap the other two coordinates, and change the sign of one of those two coordinates. Again, the vector needs to be normalized. Vectors we could choose with this approach include:

               [ vy]                  [ vz]                  [ 0 ]
by = normalize([-vx])  by = normalize([ 0 ])  by = normalize([ vz])
               [ 0 ]                  [-vx]                  [-vy]

The dot product of each of these vectors with (vx, vy, vz) is zero, which means that the vectors are orthogonal.

While the choice between these (or other variations) is mostly arbitrary, care must be taken to not end up with a degenerate vector. For example if vx and vy are both zero, using the first of these vector would be bad. To avoid choosing a (near) degenerate vector, a simple strategy is to use the first of these three vectors if vz is smaller than both vx and vy, and one of the other two otherwise.

With two new base vectors in place, the third is the cross product of the other two:

bx = by x bz

All that's left is to populate the rotation matrix with column vectors bx, by and bz, and the rotation matrix is complete:

    [ bx.x by.x bz.x ]
R = [ bx.y by.y bz.y ]
    [ bx.z by.z bz.z ]

If you need a 4x4 matrix, e.g. because you are using the legacy fixed function OpenGL pipeline, you can extend this to a 4x4 matrix:

    [ bx.x by.x bz.x 0 ]
R = [ bx.y by.y bz.y 0 ]
    [ bx.z by.z bz.z 0 ]
    [  0    0    0   1 ]