2

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

  • 采用两个单位向量：src（图中绿色）和 dst（图中红色），
  • 计算它们之间的角度：theta，
  • 计算两者的法线向量：枢轴（图中为蓝色），

  • 最后尝试通过将 src 围绕矢量枢轴旋转角度 theta 来将 src 旋转到 dst 中。

    % 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。在实现这个旋转矩阵时，您是否遇到过方向限制？如果是这样，您是如何解决这些问题的？如果没有必要，我宁愿避免额外翻译到第一象限。

谢谢，
格雷格

4

1 回答 1

3

似乎两个减号已经逃脱：

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 旋转公式的版本：

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
于 2012-11-26T08:03:00.410 回答