matlab - 这个旋转矩阵(关于向量的角度)是否仅限于某些方向？

Question

从几个引用文献(即 http://en.wikipedia.org/wiki/Rotation_matrix“轴和角度的旋转矩阵”，以及 Foley 等人的“计算机图形 - 原理与实践”中的练习 5.15，C 语言的第 2 版)，我已经看到了一个旋转矩阵(下面在 Octave 中实现)将点围绕指定向量旋转指定角度。虽然我以前使用过它，但我现在看到似乎与方向有关的旋转问题。在以下 Octave 代码中重新创建了该问题

采用两个单位向量:src(图中绿色)和 dst(图中红色)，

计算它们之间的角度:theta,

计算两者的法向向量:pivot(图中蓝色)，

最后尝试将 src 旋转成 dst，方法是将其绕矢量轴旋转角度 theta。

% This test fails: rotated unit vector is not at expected location and is no longer normalized.
s = [-0.49647; -0.82397; -0.27311]
d = [ 0.43726; -0.85770; -0.27048]
test_rotation(s, d, 1);

% Determine rotation matrix that rotates the source and normal vectors to the x and z axes, respectively.
normal = cross(s, d);
normal /= norm(normal);
R = zeros(3,3);
R(1,:) = s;
R(2,:) = cross(normal, s);
R(3,:) = normal;
R

% After rotation of the source and destination vectors, this test passes.
s2 = R * s
d2 = R * d
test_rotation(s2, d2, 2);

function test_rotation(src, dst, iFig)
    norm_src = norm(src)
    norm_dst = norm(dst)

    % Determine rotation axis (i.e., normal to two vectors) and rotation angle.
    pivot = cross(src, dst);
    theta = asin(norm(pivot))
    theta_degrees = theta * 180 / pi
    pivot /= norm(pivot)

    % Initialize matrix to rotate by an angle theta about pivot vector.
    ct = cos(theta);
    st = sin(theta);
    omct = 1 - ct;

    M(1,1) = ct - pivot(1)*pivot(1)*omct;
    M(1,2) = pivot(1)*pivot(2)*omct - pivot(3)*st;
    M(1,3) = pivot(1)*pivot(3)*omct + pivot(2)*st;
    M(2,1) = pivot(1)*pivot(2)*omct + pivot(3)*st;
    M(2,2) = ct - pivot(2)*pivot(2)*omct; 
    M(2,3) = pivot(2)*pivot(3)*omct - pivot(1)*st;
    M(3,1) = pivot(1)*pivot(3)*omct - pivot(2)*st;
    M(3,2) = pivot(2)*pivot(3)*omct + pivot(1)*st;
    M(3,3) = ct - pivot(3)*pivot(3)*omct;

    % Rotate src about pivot by angle theta ... and check the result.
    dst2 = M * src
    dot_dst_dst2 = dot(dst, dst2)
    if (dot_dst_dst2 >= 0.99999)
        "success"
    else
        "FAIL"
    end

    % Draw the vectors: green is source, red is destination, blue is normal.
    figure(iFig);
    x(1) = y(1) = z(1) = 0;
    ubounds = [-1.25 1.25 -1.25 1.25 -1.25 1.25];
    x(2)=src(1); y(2)=src(2); z(2)=src(3);
    plot3(x,y,z,'g-o');
    hold on
    x(2)=dst(1); y(2)=dst(2); z(2)=dst(3);
    plot3(x,y,z,'r-o');
    x(2)=pivot(1); y(2)=pivot(2); z(2)=pivot(3);
    plot3(x,y,z,'b-o');
    x(2)=dst2(1); y(2)=dst2(2); z(2)=dst2(3);
    plot3(x,y,z,'k.o');
    axis(ubounds, 'square');
    view(45,45);
    xlabel("xd");
    ylabel("yd");
    zlabel("zd");
    hold off
end

这是结果数字。图 1 显示了一个不起作用的方向。图 2 显示了一个有效的方向:相同的 src 和 dst 向量但旋转到第一象限。





对于所有矢量方向，我希望 src 向量始终旋转到 dst 向量上，如图 2 中黑色圆圈覆盖红色圆圈所示。然而，图 1 显示了一个方向，其中 src 向量没有旋转到 dst 向量上(即，黑色圆圈不在红色圆圈的顶部，甚至不在单位球面上)。

值得一提的是，定义旋转矩阵的引用文献没有提到方向限制，而我(在几个小时和几页内)推导出了旋转矩阵方程，并且没有发现任何方向限制。我希望这个问题是我的一个实现错误，但我还没有在我的任何一个实现中找到它:C 和 Octave。在实现此旋转矩阵时，您是否遇到过方向限制？如果是这样，你是如何解决这些问题的？如果没有必要，我宁愿避免额外翻译到第一象限。

谢谢，
格雷格

Answer 1

似乎有两个减号逃脱了:

M(1,1) = ct - P(1)*P(1)*omct;
M(1,2) = P(1)*P(2)*omct - P(3)*st;
M(1,3) = P(1)*P(3)*omct + P(2)*st;

M(2,1) = P(1)*P(2)*omct + P(3)*st;
M(2,2) = ct + P(2)*P(2)*omct;      %% ERR HERE; THIS IS THE CORRECT SIGN
M(2,3) = P(2)*P(3)*omct - P(1)*st;

M(3,1) = P(1)*P(3)*omct - P(2)*st;
M(3,2) = P(2)*P(3)*omct + P(1)*st;
M(3,3) = ct + P(3)*P(3)*omct;      %% ERR HERE; THIS IS THE CORRECT SIGN

这是一个更紧凑、更快的版本，并且同样基于 Rodrigues' rotation formula :

function test

% first test: pass
s  = [-0.49647; -0.82397; -0.27311];
d  = [ 0.43726; -0.85770; -0.27048]
d2 = axis_angle_rotation(s, d)

% Determine rotation matrix that rotates the source and normal vectors to the x and z axes, respectively.
normal = cross(s, d);
normal = normal/norm(normal);

R(1,:) = s;
R(2,:) = cross(normal, s);
R(3,:) = normal;

% Second test: pass
s2 = R * s;
d2 = R * d
d3 = axis_angle_rotation(s2, d2)

end

function vec = axis_angle_rotation(vec, dst)

    % These following commands are just here for the function to act 
    % the same as your original function. Eventually, the function is 
    % probably best defined as 
    %
    %     vec = axis_angle_rotation(vec, axs, angle)
    %
    % or even 
    %
    %     vec = axis_angle_rotation(vec, axs)
    %
    % where the length of axs defines the angle. 
    %     
    axs = cross(vec, dst);
    theta = asin(norm(axs));

    % some preparations
    aa = axs.'*axs;        
    ra = vec.'*axs;

    % location of circle centers
    c = ra.*axs./aa;

    % first coordinate axis on the circle's plane
    u = vec-c;

    % second coordinate axis on the circle's plane
    v = [axs(2)*vec(3)-axs(3)*vec(2)
         axs(3)*vec(1)-axs(1)*vec(3)
         axs(1)*vec(2)-axs(2)*vec(1)]./sqrt(aa);

    % the output vector   
    vec = c + u*cos(theta) + v*sin(theta);        

end