Deriving axis/angle rotation from two pairs of three points (or, two pairs of two vectors)

Question

Please note that although it sounds similar, this is not the common "how to rotate one vector onto another" question.

I would like to derive an affine transform (either in matrix or a quaternion+vector form) from two sets of 3 points. Those can be regarded as "marker points" on rigid bodies, or as the endpoints of "forward and up" vectors. Translation and rotation is necessary, scaling is not necessary. Also, a quaternion+vector solution would be a plus, since it would allow me to cram 1/3 more instances into a drawing batch (8 uniforms instead of 12). The intent is to have a system for determining the pose of (articulated or not) ridid bodies in an intuitive way, without needing to maintain and walk a complicated hierarchy.

The first obvious simplification is to eleminate the translational part by picking one of the points and subtracting the "destination" from the respective "start" point. Now we only need to deal with a rotation.

There is a well-known, computionally inexpensive solution of constructing a quaternion that rotates one vector onto another, namely q(cross(v1,v2) ; sqrt(v1.len_sq * v2.len_sq) + dot(v1,v2)) or q(cross(v1,v2) ; 1 + dot(v1,v2)) for unit-length vectors. Unluckily, this method has no notion of an "up direction", and therefore always roates on the shortest arc (which will misalign objects). The naive thing to do would be to simply use this method for both vectors and multiply the quaternions together, but it obviously won't work that easily. What one would need to do is pick one of the two vectors (let's call that one "forward"), and create a quaternion for this one, then rotate the other ("up") vector using this quaternion, then construct a second quaternion for the rotated "up" vector (and the target "up" vector), and finally multiply the second to the first quaternion. This will be correct as far as I can tell, but it is also horribly complicated.

Now... as for rotation matrices, I am aware of the "triad method" which I understand as follows: - Orthonormalize the vector pairs (both start and end) - This results in two orthonormal bases which are the respective rotational matrices for start and end from a "common reference frame". It does not matter what reference frame exactly this is, all that matters is that it is the same for both. - S is the transform from the "common frame" to the start frame, and D is the transform to the end frame, respectively. - Therefore, S^-1 * D * v transforms any point from the start to the end coordinate system (going via the common reference frame). - S^-1 == S^T since it is an orthonormal matrix, and S^T * x = x * S - Therefore: S^T * D * v = D * S * v

This should work, but it still seems quite complicated for something that should actually be really, really simple.

Is there an easier, more straightforward solution?

Answer 1

Your "horribly complicated" quaternion solution would not generally work. You would have to project the second pair of vectors on the plane orthogonal to the axis of the first rotation to make sure that the second rotation is orthogonal to the first one.

I have described the principle in my blog if you are interested: http://robokitchen.tumblr.com/post/67060392720/finding-a-rotation-quaternion-from-two-pairs-of-vectors

Before rotation: u0, v0. After rotation: u2, v2.

Quaternion q2 = Quaternion::fromTwoVectors(u0, u2);
Vector v1 = v2.rotate(q2.conjugate());
Vector v0_proj = v0.projectPlane(u0);
Vector v1_proj = v1.projectPlane(u0);
Quaternion q1 = Quaternion::fromTwoVectors(v0_proj, v1_proj);
return (q2 * q1).normalized();

I'm not sure this is the solution you want but the code runs surprisingly fast.

Answer 2

To deal with the rotation part only, your second method will work and I suspect it will work well. Alternatively, you could use a hybrid of the two methods which might be a bit easier. Assume the two pairs of two vectors you constructed above, each pair in its own vector space. Compute the orthormal basis of each pair and call them X₀ and X₁ in one vector space, and their corresponding vectors Y₀ and Y₁ in the other vector space. You now have to compute two quaternion rotations:

1) q₀ rotates X₀, and X₁ to X'₀ and X'₁ respectively such that X'₀ = Y₀. X'₁ and Y₁ should now be coplanar with plane normal X'₀ = Y₀.

2) q₁ rotates X'₁ to X''₁ = Y₁. All you need to do is compute the angle between the vectors since you already know the rotation vector will just be X'₁ x Y₁ = X'₀ = Y₀

You can calculate q = q₁ * q₀ to perform the rotation in a single step.