Calculate rotation matrix to transform one vector to another?

Question

There are several questions about this on Stack Overflow, but most use Quaternions and I am looking for something that does not. And this question is all over the web, but it is surprisingly hard to find a straightforward example in code. I am looking for a C/C++, or GLSL/Metal solution for the rotation matrix to transform one vector to another. I found this which seems to answer the question: http://www.j3d.org/matrix_faq/matrfaq_latest.html#Q38

However, it is not working for me in practice. With this I am getting the correct transformations at the ninety-degree angles (ie transforming 0,0,1 to 0,1,0), but in between it twists in the wrong direction. Here is the code I am using -- does anyone see what is wrong with it?

In this case I am always transforming from a fixed start normal of (0,0,1).

// matrix to rotate from (0,0,1) to specified direction

float4x4 directionMatrix(float3 direction) {

    float3 normal = float3(0.0,0.0,1.0);
    float3 V = normalize(cross(normal, direction));
    float phi = acos(dot(normal,direction));

    float rcos = cos(phi);
    float rsin = sin(phi);

    float4x4 M;

    M[0].x =        rcos + V.x * V.x * (1.0 - rcos);
    M[1].x =  V.z * rsin + V.y * V.x * (1.0 - rcos);
    M[2].x = -V.y * rsin + V.z * V.x * (1.0 - rcos);
    M[3].x = 0.0;

    M[0].y = -V.z * rsin + V.x * V.y * (1.0 - rcos);
    M[1].y =        rcos + V.y * V.y * (1.0 - rcos);
    M[2].y = -V.x * rsin + V.z * V.y * (1.0 - rcos);
    M[3].y = 0.0;

    M[0].z =  V.y * rsin + V.x * V.z * (1.0 - rcos);
    M[1].z = -V.x * rsin + V.y * V.z * (1.0 - rcos);
    M[2].z =        rcos + V.z * V.z * (1.0 - rcos);
    M[3].z = 0.0;

    M[0].w = 0.0;
    M[1].w = 0.0;
    M[2].w = 0.0;
    M[3].w = 1.0;

    return M;
}

I tried transposing this, in case the original article used the opposite column-row majority of mine, but it did not help.

EDIT: The problem here turned out to be that this code works for a left-handed coordinate system, where mine is right-handed. So to get this to work you just have to multiply phi by -1. I also edited the code to add normalization of the cross vector, as suggested in the comments.

You could simplify to `rcos=dot(normal,direction)` and `rsin=norm(V)`. — Oliv, Sep 05 '18 at 17:09
Then to solve your problem, add this line before computing matrix elements: `V=V/norm(V)`. The axis vector must be normalized ;)! — Oliv, Sep 05 '18 at 17:12
Could this be made for a left-handed coordinate system? I am using Metal and I think I have a right-handed coordinate system. Would that make a difference? — bsabiston, Sep 05 '18 at 21:37
Yes, that was the problem! I read that in order to switch between right and left-handed systems, you can just reverse the angle of rotation. So I multiply phi by -1 and it works, finally. Whew. — bsabiston, Sep 05 '18 at 21:46
An other way to change orientation is to reverse the cross product: `V=cross(direction,normal)`. — Oliv, Sep 06 '18 at 07:03
Also consider the special cases where `norm(V)<1e-6`. And finaly, also take into account error propagation, the computation can be simplified if you consider the fact that `rsin=norm(V)`. — Oliv, Sep 06 '18 at 07:05
what do you mean by norm()? Metal does not have that. I was thinking it was the same thing as length(), except that does not work. I am fine with sin() as long as it works. — bsabiston, Sep 06 '18 at 21:40
Note that there is no *unique* way to rotate from one vector to another: there is always one free continuous parameter, or two if the vectors are opposite. In the one-parameter case, there is, however, a *smallest* rotation, which is usually what is meant. — Davis Herring, Sep 06 '18 at 23:52
FWIW, all the trigonometry can be optimized away: https://iquilezles.org/articles/noacos/ — Tom, Jul 24 '23 at 07:00

Calculate rotation matrix to transform one vector to another?

0 Answers0