Rotate a vector about another vector

I am writing a 3d vector class for OpenGL. How do I rotate a vector v1 about another v开发者_Python百科ector v2 by an angle A?

You may find quaternions to be a more elegant and efficient solution.

After seeing this answer bumped recently, I though I'd provide a more robust answer. One that can be used without necessarily understanding the full mathematical implications of quaternions. I'm going to assume (given the C++ tag) that you have something like a Vector3 class with 'obvious' functions like inner, cross, and *= scalar operators, etc...

#include <cfloat>
#include <cmath>

...

void make_quat (float quat[4], const Vector3 & v2, float angle)
{
    // BTW: there's no reason you can't use 'doubles' for angle, etc.
    // there's not much point in applying a rotation outside of [-PI, +PI];
    // as that covers the practical 2.PI range.

    // any time graphics / floating point overlap, we have to think hard
    // about degenerate cases that can arise quite naturally (think of
    // pathological cancellation errors that are *possible* in seemingly
    // benign operations like inner products - and other running sums).

    Vector3 axis (v2);

    float rl = sqrt(inner(axis, axis));
    if (rl < FLT_EPSILON) // we'll handle this as no rotation:
    {
        quat[0] = 0.0, quat[1] = 0.0, quat[2] = 0.0, quat[3] = 1.0;
        return; // the 'identity' unit quaternion.
    }

    float ca = cos(angle);

    // we know a maths library is never going to yield a value outside
    // of [-1.0, +1.0] right? Well, maybe we're using something else -
    // like an approximating polynomial, or a faster hack that's a little
    // rough 'around the edge' cases? let's *ensure* a clamped range:
    ca = (ca < -1.0f) ? -1.0f : ((ca > +1.0f) ? +1.0f : ca);

    // now we find cos / sin of a half-angle. we can use a faster identity
    // for this, secure in the knowledge that 'sqrt' will be valid....

    float cq = sqrt((1.0f + ca) / 2.0f); // cos(acos(ca) / 2.0);
    float sq = sqrt((1.0f - ca) / 2.0f); // sin(acos(ca) / 2.0);

    axis *= sq / rl; // i.e., scaling each element, and finally:

    quat[0] = axis[0], quat[1] = axis[1], quat[2] = axis[2], quat[3] = cq;
}

Thus float quat[4] holds a unit quaternion that represents the axis and angle of rotation, given the original arguments (, v2, A).

Here's a routine for quaternion multiplication. SSE/SIMD can probably speed this up, but complicated transform & lighting are typically GPU-driven in most scenarios. If you remember complex number multiplication as a little weird, quaternion multiplication is more so. Complex number multiplication is a commutative operation: a*b = b*a. Quaternions don't even preserve this property, i.e., q*p != p*q :

static inline void
qmul (float r[4], const float q[4], const float p[4])
{
    // quaternion multiplication: r = q * p

    float w0 = q[3], w1 = p[3];
    float x0 = q[0], x1 = p[0];
    float y0 = q[1], y1 = p[1];
    float z0 = q[2], z1 = p[2];

    r[3] = w0 * w1 - x0 * x1 - y0 * y1 - z0 * z1;
    r[0] = w0 * x1 + x0 * w1 + y0 * z1 - z0 * y1;
    r[1] = w0 * y1 + y0 * w1 + z0 * x1 - x0 * z1;
    r[2] = w0 * z1 + z0 * w1 + x0 * y1 - y0 * x1;
}

Finally, rotating a 3D 'vector' v (or if you prefer, the 'point' v that the question has named v1, represented as a vector), using the quaternion: float q[4] has a somewhat strange formula: v' = q * v * conjugate(q). Quaternions have conjugates, similar to complex numbers. Here's the routine:

static inline void
qrot (float v[3], const float q[4])
{
    // 3D vector rotation: v = q * v * conj(q)

    float r[4], p[4];

    r[0] = + v[0], r[1] = + v[1], r[2] = + v[2], r[3] = +0.0;
    glView__qmul(r, q, r);

    p[0] = - q[0], p[1] = - q[1], p[2] = - q[2], p[3] = q[3];
    glView__qmul(r, r, p);

    v[0] = r[0], v[1] = r[1], v[2] = r[2];
}

Putting it all together. Obviously you can make use of the static keyword where appropriate. Modern optimising compilers may ignore the inline hint depending on their own code generation heuristics. But let's just concentrate on correctness for now:

How do I rotate a vector v1 about another vector v2 by an angle A?

Assuming some sort of Vector3 class, and (A) in radians, we want the quaternion representing the rotation by the angle (A) about the axis v2, and we want to apply that quaternion rotation to v1 for the result:

float q[4]; // we want to find the unit quaternion for `v2` and `A`...

make_quat(q, v2, A);

// what about `v1`? can we access elements with `operator [] (int)` (?)
// if so, let's assume the memory: `v1[0] .. v1[2]` is contiguous.
// you can figure out how you want to store and manage your Vector3 class.

qrot(& v1[0], q);

// `v1` has been rotated by `(A)` radians about the direction vector `v2` ...

Is this the sort of thing that folks would like to see expanded upon in the Beta Documentation site? I'm not altogether clear on its requirements, expected rigour, etc.

This may prove useful:

double c = cos(A);
double s = sin(A);
double C = 1.0 - c;

double Q[3][3];
Q[0][0] = v2[0] * v2[0] * C + c;
Q[0][1] = v2[1] * v2[0] * C + v2[2] * s;
Q[0][2] = v2[2] * v2[0] * C - v2[1] * s;

Q[1][0] = v2[1] * v2[0] * C - v2[2] * s;
Q[1][1] = v2[1] * v2[1] * C + c;
Q[1][2] = v2[2] * v2[1] * C + v2[0] * s;

Q[2][0] = v2[0] * v2[2] * C + v2[1] * s;
Q[2][1] = v2[2] * v2[1] * C - v2[0] * s;
Q[2][2] = v2[2] * v2[2] * C + c;

v1[0] = v1[0] * Q[0][0] + v1[0] * Q[0][1] + v1[0] * Q[0][2];
v1[1] = v1[1] * Q[1][0] + v1[1] * Q[1][1] + v1[1] * Q[1][2];
v1[2] = v1[2] * Q[2][0] + v1[2] * Q[2][1] + v1[2] * Q[2][2];

Use a 3D rotation matrix.

The easiest-to-understand way would be rotating the coordinate axis so that vector v2 aligns with the Z axis, then rotate by A around the Z axis, and rotate back so that the Z axis aligns with v2.

When you have written down the rotation matrices for the three operations, you'll probably notice that you apply three matrices after each other. To reach the same effect, you can multiply the three matrices.

I found this here: http://steve.hollasch.net/cgindex/math/rotvec.html

let
    [v] = [vx, vy, vz]      the vector to be rotated.
    [l] = [lx, ly, lz]      the vector about rotation
          | 1  0  0|
    [i] = | 0  1  0|           the identity matrix        
          | 0  0  1|

          |   0  lz -ly |
    [L] = | -lz   0  lx |
          |  ly -lx   0 |

    d = sqrt(lx*lx + ly*ly + lz*lz)
    a                       the angle of rotation

then

matrix operations gives:

[v] = [v]x{[i] + sin(a)/d*[L] + ((1 - cos(a))/(d*d)*([L]x[L]))} 

I wrote my own Matrix3 class and Vector3Library that implemented this vector rotation. It works absolutely perfectly. I use it to avoid drawing models outside the field of view of the camera.

I suppose this is the "use a 3d rotation matrix" approach. I took a quick look at quaternions, but have never used them, so stuck to something I could wrap my head around.