3D rotation algorithm needed

Given two orthogonal unit vectors, A and B, and two different orthogonal unit vectors C and D, I need the 3x3 direction cosine matrix or a quaternion which will rotate A to align with C AND will rotate B to align with D.

The vectors are all 3-vectors (x, y, z).

I have a brute force algorithm, but I am almost certain there is a much simpler solution. And my web searches have not yet revealed one.

I am using C#, but if you have something in C, or FORTRAN, or Basic (etc.) I can convert it. Or, I can use the terms written out mathematically.

The application is determining the orientation needed for a spacecraft so an optical device attached rigidly to it will be properly aligned to take a picture. (Both the bore sight direction and the appropriate rotation of the optic device about bore sight must be attained, thus the need to align TWO vect开发者_运维技巧ors at once.) The computation may be used in a loop fed by real-time sensed data, and the brute force solution is too slow.

I reread your question and the answer below (while right) doesn't give you what you want. This link is about constructing a 3x3 rotation matrix.

since they are both orthogonal unit vectors, you just need to add one more too each to construct a basis (use the cross product for this). So now you have two basis {A, B, AxB} and {C, D, CxD}. A rotation that moves {A, B} onto {C, D} will re-express a vector a1A + a2B + a3(AXB) as b1C + b2D + b3(CxD). because it's linear you only need to know how it behaves on the basis (This uniquely determines the linear transformation). So, taking {A, B, ..} as our basis and letting the transformation be T, we see that T(1, 0, 0) = C, T(0, 1, 0) = D and T(0, 0, 1) = CxD. Remember A = (1, 0, 0) etc. But the columns of this matrix are just M =(C, D, CxD)

To use this matrix as it stands, you have to express every vector in the basis {A, B, CxD} before you left-multiply it by M. You do this in the same way. In fact, it N is the matrix which translates from your normal basis to {A, B, ..} and M above translates that to {C, D...}, then MN (left multiplication here) will translate from your basis to {C, D, ..} and provide the rotation you want.

So now, all of your vectors are expressed in the basis {C, D, ..} :(

The solution is yet another transformation matrix. This one maps from {A, B, ..} to your primary basis and undoes N, also known as the inverse and denoted N^-1. So your final matrix is (N^-1)MN. The good news is that because N is orthogonal, you just need it's transpose.

The trick is to choose your primary basis so that the matrices you deal with most are pretty.

https://alleg.svn.sourceforge.net/svnroot/alleg/allegro_outdated/branches/allegro/src/math3d.c

almost complete code for 3d graphics....check out functions at the bottom starting with..

get_align_matrix_f

and here goes same for quaterian...

https://alleg.svn.sourceforge.net/svnroot/alleg/allegro_outdated/branches/allegro/src/quat.c

also in matrix, it may not give you shortest ( or direct ) path from vactor A to C, so if you are animating the visuals, its better to use quats.

Given that spacecraft cost hundreds of millions of dollars, you might want to find someone who can do this sort of thing in their sleep and ask them to generate a bulletproof and optimized solution, with code, rather than relying on descriptions here. (Unless this is just an exercise.)

Also, the solution you pick should depend very heavily on the thrusters available on the spacecraft; you want the rotation to use as little fuel as possible, and that's going to depend on what they intrinsically do to the spacecraft. I assume that you've set the problem up so that one axis of rotation is around the z axis; can you rotate independently around x and y also, or do you only have one other axis? Are some rotations more expensive than others (e.g. due to a different moment of inertia along some axes)?

Given these caveats, I'm pretty hesitant to give advice like: find (A x C), the axis of rotation needed to move A onto C and rotate (see Wikipedia for rotation along an axis). Then figure out what that does to B (by multiplying the rotation matrix by B) and compute the angle between B and D; finally, rotate along the (B x D) axis--which had better be the same as the C axis at this point--to fix that difference (giving another rotation matrix). Mutliply the two matrices, and voila, you're done.

Somewhat after the original post I know but in case you too are suffering from column/row major and left/right hand confusion issues, if you wanted to encode Mark T's splendid result into WPF you'd do it like this:

static public Matrix3D TwistToNewAxes(Vector3D A, Vector3D B, Vector3D D, Vector3D E)
    {
        Vector3D C = Vector3D.CrossProduct(B, A);
        Vector3D F = Vector3D.CrossProduct(E, D);

        Matrix3D result = Matrix3D.Identity;

        result.M11 = D.X * A.X + E.X * B.X + F.X * C.X;
        result.M21 = D.X * A.Y + E.X * B.Y + F.X * C.Y;
        result.M31 = D.X * A.Z + E.X * B.Z + F.X * C.Z;
        result.M12 = D.Y * A.X + E.Y * B.X + F.Y * C.X;
        result.M22 = D.Y * A.Y + E.Y * B.Y + F.Y * C.Y;
        result.M32 = D.Y * A.Z + E.Y * B.Z + F.Y * C.Z;
        result.M13 = D.Z * A.X + E.Z * B.X + F.Z * C.X;
        result.M23 = D.Z * A.Y + E.Z * B.Y + F.Z * C.Y;
        result.M33 = D.Z * A.Z + E.Z * B.Z + F.Z * C.Z;

        return result;


    }

Many thanks Mark!