2

いくつかの参考文献（すなわち、http：//en.wikipedia.org/wiki/Rotation_matrix「軸と角度からの回転行列」、およびFoley etalによる「ComputerGraphics-PrinciplesandPractice」の演習5.15、Cの第2版から） ）、指定されたベクトルを中心に指定された角度で点を回転させる回転行列（以下のオクターブで実装）のこの定義を見てきました。以前は使用していましたが、向きに関連しているように見える回転の問題が発生しています。問題は、次のOctaveコードで再現されます。

  • src（図では緑）とdst（図では赤）の2つの単位ベクトルを取ります。
  • それらの間の角度を計算します：シータ、
  • 両方に垂直なベクトルを計算します：ピボット（図では青）、

  • 最後に、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ベクトルですが、第1象限に回転しています。

ここに画像の説明を入力してください

ここに画像の説明を入力してください

図2に示すように、すべてのベクトルの方向で、srcベクトルが常にdstベクトル上で回転し、赤い円を覆う黒い円で示されることを期待していました。ただし、図1は、srcベクトルがdstベクトル上で回転しない方向を示しています（つまり、黒い円は赤い円の上になく、単位球上にもありません）。

価値のあることとして、回転行列を定義した参考文献は方向の制限について言及していませんでした。私は回転行列の方程式を（数時間と数ページで）導き出し、そこで方向の制限を見つけませんでした。問題が私の側の実装エラーであることを望んでいますが、CとOctaveのどちらの実装でもまだそれを見つけることができませんでした。この回転行列を実装するときに方向の制限を経験しましたか？もしそうなら、どのようにそれらを回避しましたか？必要がなければ、第1象限への余分な変換は避けたいと思います。

ありがとう、
グレッグ

4

1 に答える 1

3

2つのマイナス記号が逃げたようです：

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

これは、はるかにコンパクトで高速で、ロドリゲスの回転式に基づいたバージョンです。

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 に答える