Q:

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

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 == ST since it is an orthonormal matrix, and ST * x = x * S
- Therefore: ST * 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?

A:

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 X0 and X1 in one vector space, and their corresponding vectors Y0 and Y1 in the other vector space.  You now have to compute two quaternion rotations:
1) q0 rotates X0, and X1 to X'0 and X'1 respectively such that X'0 = Y0.  X'1 and Y1 should now be coplanar with plane normal X'0 = Y0.
2) q1 rotates X'1 to X''1 = Y1.  All you need to do is compute the angle between the vectors since you already know the rotation vector will just be X'1 x Y1 = X'0 = Y0
You can calculate q = q1 * q0 to perform the rotation in a single step.

