Quaternion from two 3D Position Vectors

[Jenner]

This came up in the Games forum, I thought I'd post it here since it has a great deal of mathematical relevance.

Given two 3D Position Vectors indicating the positions of two objects on a 3D grid, is it possible to create a quaternion to describe the absolute rotation required for object 1 to be pointing at object 2?

And if so, how?

I've found about 3 or 4 poorly described code samples trying to do this, but they all either error or produce undesirable results, or break in certain conditions (such as when object 1 is at the center point of the grid).

[jemidiah]

Yes, it's possible, presuming a rotation that takes object 1 to object 2 exists. I've never done it myself, but I can explain the Wikipedia page on it, which I presume is correct.

You said that you found code that breaks when object 1 is at the center point of the grid--what would you want it to do there? Rotating anything at the center any amount doesn't do anything.

Using the procedure described on the page I linked:

Mechanically, let p1 and p2 be vectors pointing to objects 1 and 2 from the origin. To rotate p1 onto p2, you'd need to know two things: an axis of rotation, and an amount of rotation. We haven't gotten to quaternions yet.

The axis of rotation must be perpendicular to both p1 and p2--that is, the axis of rotation is the normal vector of the plane containing p1 and p2. A convenient axis of rotation, then, is their cross product, p1 x p2, which by definition yields a vector normal to the plane of the two arguments. We'll want this to be a unit vector later, so define u = axis of rotation = p1 x p2 / |p1 x p2|, where |x| is the magnitude of x.

Now we need an amount of rotation. This is actually pretty easy--we just need to know the angle, in the plane containing both vectors, between them. The dot product handily solves this for us, since p1 dot p2 = |p1||p2|cos(alpha) where alpha is the angle between them. Solve for alpha simply with alpha = arccos(p1 dot p2 / |p1||p2|).

Now comes the quaternion. From that page, in general a quaternion has 4 components,
q = w + xi + yj + zk, where i, j, and k are similar to imaginary components of a complex number. To be honest, for your use, the specifics just don't matter, so I'll give you the algorithmic version.

A rotation quaternion corresponding to a rotation in 3 dimensions about axis u=(ux, uy, uz) [where u is taken to be a unit vector for later convenience], by an angle alpha, is defined to be, using notation that gives the 4 components of that quaternion (which is how you'd store it in an array)

q = <cos(alpha/2), sin(alpha/2)*ux, sin(alpha/2)*uy, sin(alpha/2)*uz>

To rotate a vector p1 by this amount through this angle, compute p2' = q*p1'*(q^-1), where the p1' = <0, p1x, p1y, p1z>, p2' = <0, p2x, p2y, p2z>. The multiplication is quaternion multiplication; pick off the last 3 coordinates of the resulting quaternion p2' to get p2.


To summarize:
1. Find your axis of rotation as a unit vector. Computationally, compute u=unitize(p1 x p2).
2. Find your angle of rotation. Computationally, alpha = arccos(p1 dot p2 / |p1||p2|).
3. Compute your rotation quaternion. Computationally, q = the expression above = cos(alpha/2) + u'*sin(alpha/2), where u' = <0, ux, uy, uz>.
4. Compute the multiplicative inverse of your rotation quaternion. We've made q a unit quaternion to simplify this step slightly. Computationally, if q = <qw, qx, qy, qz> is a rotation quaternion, then q^-1 = <qw, -qx, -qy, -qz>.
5. "Quaternionize" your first object's vector. Computationally, if p1 = <p1x, p1y, p1z>, then p1' = <0, p1x, p1y, p1z>.
6. Compute the conjugate of p1' using the rotation quaternion. Computationally, using quaternion multiplication, this is p2' = q*p1'*(q^-1).
7. "Unquaternionize" your result from the previous step to find the rotated vector. Computationally, at this point p2' = <0, p2x, p2y, p2z>; simply pick off the last 3 components to form p2 = <p2x, p2y, p2z>.

Helper functions you'll need:
a. Cross Product; definition is easy to find
b. Dot Product; even easier
c. Unitize; for any vector v, unitize(v) = v / |v|.
d. Length of a vector, i.e. |v|; for any vector v, |v| = sqrt(v dot v).
e. Quaternionize; described above; turns a vector into a quaternion.
f. Unquaternionize; see (e), inverted.
g. Vector scaling for (c).
h. Quaternion multiplication; defined on the page I linked.
g. Possibly quaternion addition, depending on how you calculate q.

Is this clear? Again I haven't implemented this, but the math isn't that bad (at least to me; I know for people with less experience in abstract algebraic spaces this must be awful).


Edit: I felt like writing a routine for this, assuming the above helper functions are defined. Pseudocode used.

    Rotate3D(point1 p1 as Vector3D, axis u as Vector3D, angle alpha as Float) As Vector3D
    'Rotates the given point p1 through axis u an angle alpha; returns the resulting point.
        Dim q As Quaternion
        u = Unitize(u)
        q = cos(alpha/2) + sin(alpha/2)*Quaternionize(u)
        Rotate3D = Unqaternionize(q*Quaternionize(p1)*q^-1)
    End Rotate3D

    RotateP1toP2(point p1 as Vector3D, point p2 As Vector3D)
        new_p2 = Rotate3D(p1, p1 x p2, arccos(unitize(p1) dot unitize(p2))
        'At this point, new_p2 should be equal to p2, up to numerical
        'accuracy of the computations. Note that I've neglected edge
        'cases--I've assumed all of the above quantities exist and are
        'well-defined.
    End RotateP1ToP2

[Jenner]

You said that you found code that breaks when object 1 is at the center point of the grid--what would you want it to do there? Rotating anything at the center any amount doesn't do anything.

Say you have an 3D coordinate grid with two spaceships on the grid.
The first spaceship is at the center point of the grid, the origin, thus it has a positional vector of (0, 0, 0). Spaceship 2 is at (5, 10, -15) i.e. the lower first quadrant. Some of the routines I found break in this scenario, since the first object's position is at (0, 0, 0), yet, the object can still rotate to be facing the second spaceship.

Your breakdown of the wiki is invaluable. Quite honestly, math notation is something I've always considered myself weak on, and throughout college have constantly struggled until I got to the point where real numbers were involved; which led to many epiphany of "Oh! that's so friggin simple! Why did he waste three weeks before showing one real-world example of it's use!".

I didn't understand any of my first three classes of Calculus beyond the regurgitation I needed for a passing grade until I got to the fourth class: Differential Equations; and since we were finally solving real-number systems and real-world scenarios, everything finally made sense to me! Up until that point, I viewed the subject as useless to me. After that point, Calculus was as easy as multiplication.

I'm going to take your explanation and see if I can't create a new routine that does what I want it to do properly. Thanks much!

[jemidiah]

Some of the routines I found break in this scenario, since the first object's position is at (0, 0, 0), yet, the object can still rotate to be facing the second spaceship.

I hope I haven't misconstrued your question. The algorithm I listed takes a point and rotates it some amount around an arbitrary axis, much like this illustration shows. It uses P1 as I've used it, and its P3 is my P2. The "rotation axis" is my u vector. Alpha is identical.

It has nothing to do with making P1 somehow face P2. I can derive some formulas to tell you which direction to turn to face P2, but that will probably have little to do with quaternions.


Lemme explain my routine another way. Perhaps you're familiar with polar notation for complex numbers. That is, for any complex number a+bi, we can write that number as Re^(i theta) for some numbers R and theta that depend on a and b--it turns out that R and theta correspond exactly to polar coordinates in the complex plane. Anywho, if you take R=1, we're just left with e^(i theta). Ok.... Now let's say I have a vector in a plane, say (1, 1), and I want it rotated 90 degrees [counterclockwise]. How could I do this?

One easy way is to take (1 + 1i) * e^(i*90 degrees) = (1 + i)*(cos(90 deg) + i sin(90 deg)) = (1 + i)*(0+i) = (i - 1) = (-1 + 1i); picking off the coordinates, (1, 1) rotated by 90 degrees is (-1, 1). Drawing it out in my mind--yup, it worked. Magical--but not quite. (1+1i) has some representation in polar coordinates; it so happens that (1+1i) = Sqrt(2)*e^(i*45 degrees). Multiply this by e^(i*90 degrees) to get Sqrt(2)*e^(i*135 degrees)--but this is just our original number with an extra 90 degrees added on to the angle, which is exactly a rotation. Quaternions generalize this operation to three dimensional space (more or less). In 3D, you need to specify an axis of rotation as well as an amount of rotation, so things get significantly more complicated (hence, why quaternions are somewhat obtuse).