MATLABで楕円体を並進・回転させて描画する

MATLABで楕円体を同次変換するときに少々手こずったのでメモしておきます．

楕円体とは

楕円体は，2次曲面の一種です．2次元において，次の方程式：

で表現される図形を楕円と呼びますが，これの次元へ拡張したものと捉えて問題ありません． より厳密な呼び分けとしては，のときのみ楕円体と呼び，のとき超楕円体と呼ぶ場合もあるようです． 分野によって呼称の差異があるようですのでご注意ください．

（3次元の）楕円体

今回取り扱うのはの場合です．残念ながら，普通の静的な描画では，3次元が限度ですので悪しからず． 楕円体の方程式は次のように表されます．

ここでは基準座標系の各軸に対応する方向の径です．

いずれか2つの径が等しいとき，つまりのいずれかの条件が満たされているとき，図形は特に回転楕円体と呼ばれます． 地球は回転楕円体としてよく近似されます．

また全ての径が等しいとき，つまりが満たされているとき，その楕円体は球となります．

以下にの場合の楕円体を描画したものを載せます．

MATLABにおける楕円体の生成

任意の原点における楕円体

楕円体を描画するためには面の情報が必要です． そのため，plot3よりsurfの方が適していると言えるでしょう．

surfに与える面の情報は，ellipsoid関数を使うことで生成できます． この関数は，楕円体の原点の座標と，先程で記述した各座標に対応する径を引数として，面の情報であるメッシュを生成します．

原点の座標をそれぞれx0, y0, z0とし，対応する径をa, b, cとすれば，次の2行で簡単に楕円体が描画できてしまうわけですね．

[x,y,z] = ellipsoid(x0, y0, z0, a, b, c);
surf(x, y, z)

サーフェスの回転

さて，ellipsoid関数は，の各軸方向を主軸とする楕円体を生成しました．

しかし，一般に楕円体は，主軸方向がある座標系と常に一致するとは限りません． 今回は3次元なので，楕円体に3次の回転行列()をかける操作によって，任意の主軸方向での楕円体を作り出すことができます．

注意するべき点は，

1. 今回扱うデータが表面であること
2. 回転の原点が座標系の原点と一致しているとは限らないこと

です．

これらは，rotate関数を使うことで解決できます． rotate関数はグラフィックスオブジェクトに3次元の回転を作用させることができ，また回転の原点も引数として設定が可能です．

ただし，この関数は，オイラー角のような，直接的なパラメトリック指定ができません． rotate関数に与えられるのは回転軸のベクトルと，その軸回りの回転量です．

もし回転行列が既知であり，かつ単位行列（回転量がゼロ）でなければ，Rodriguesの式（回転公式）の対数を取ることでこれらを導けます．

R: 3d rotational matrix

theta = acos((trace(R) - 1) / 2);
a = [R(3,2) - R(2,3);
     R(1,3) - R(3,1);
     R(2,1) - R(1,2)] / 2 / sin(theta);

ここでaが正規化された回転軸のベクトル，thetaがその軸回りの回転量となります．

先程の楕円体表面をs = surf(x, y, z);のように変数として保持しておけば，

rotate(s, a, theta*180/pi, [x0 y0 z0])

とすることで，楕円体の中央を回転中心として，回転行列Rだけ楕円体を回転させることができます． 回転量が，引数では degreeであることに注意が必要です．

サンプル

楕円体中心をとして，その場で回転させます．

コード

clear;
close all;

% center position
x0 = 2;
y0 = 3;
z0 = 0;

% length of axes
a = 1;
b = 2;
c = 3;

% creates an ellipsoid
[x,y,z] = ellipsoid(x0,y0,z0,a,b,c);

% draws the surface
figure;
s = surf(x,y,z);

% rotational parameters
R = Rx(pi/6) * Ry(pi/6) * Rz(pi/4);

theta = acos((trace(R) - 1) / 2);

a = [R(3,2) - R(2,3);
     R(1,3) - R(3,1);
     R(2,1) - R(1,2)] / 2 / sin(theta);

% rotates
rotate(s,a,theta*180/pi,[x0 y0 z0])

% figure settings
set(groot,'DefaultTextInterpreter', 'Latex');

grid on;
xlabel('$x$','FontSize',16)
ylabel('$y$','FontSize',16)
zlabel('$z$','FontSize',16)
axis equal
view([10;10;10])


%% functions
function R = Rx(a)
R = [1     0        0;
     0 cos(a) -sin(a);
     0 sin(a)  cos(a)];
end

function R = Ry(a)
R = [ cos(a) 0 sin(a);
          0  1      0;
     -sin(a) 0 cos(a)];
end

function R = Rz(a)
R = [cos(a) -sin(a) 0;
     sin(a)  cos(a) 0;
          0       0 1];
end

出力

参考文献

jp.mathworks.com

jp.mathworks.com