FC2ブログ

初心者によるMATLABメモ

3次元形状やってんのになぜかMATLABを使わされることになった人のメモ。

 

検索エンジンから来た方は読みたい記事が表示されていない可能性があります。
その場合、左メニューのフォームでブログ内を検索すると見たい記事を見れると思います。
もしくは、カテゴリー名をクリックしていただくと各カテゴリーの記事のリストが見れます。

[基本] データの極大値,極小値を取り出す

データの極大値,極小値のインデクスを取り出す方法です。(最小二乗近似曲線ならこちら
datadiff = diff(DATA);
極大値なら
  ind = find(datadiff(2:end).*datadiff(1:end-1)<0 & 0<datadiff(1:end-1)) + 1;
極小値なら
  ind = find(datadiff(2:end).*datadiff(1:end-1)<0 & datadiff(1:end-1))<0 + 1;
として
DATA(ind)

実際に使うならこんな風に両端のデータを加えた方がいいのかも…
ind = [1, ind, length(DATA)];

以下に,実際に例とか。適当なデータを作ってその極大値だけを繋いでプロットしてみる例。
まずこんな感じの適当なデータを作ってみた。
% 適当なデータ
nn = 50;
data = 4*[[1:nn], nn-[1:nn]]/nn;
noise = cos(pi*rand(1, nn*2)*(51/100))/2;
data =data+noise;

プロットするとこんな感じに。
極大値1

で,上に書いた方法で極大値を探索。これで極大値のインデクスが得られます。
% 差分
datadiff = diff(data);

% 極大値探索
ind = find(datadiff(2:end).*datadiff(1:end-1)<0 & 0<datadiff(1:end-1)) + 1;
ind = [1, ind, length(data)];

後はこれをプロットしてみるとこんな感じに。
グラフは上から「データ全体のプロット」,「データ全体のプロットに極大値プロットを重ねてみたもの」,「極大値プロット」です。
まぁ結果は見たまんまな感じで。
極大値2

最後に参考までにこのグラフを書いたコードを載せておきます。
close all;

% 適当なデータ
nn = 50;
data = 4*[[1:nn], nn-[1:nn]]/nn;
noise = cos(pi*rand(1, nn*2)*(51/100))/2;
data =data+noise;

% 差分
datadiff = diff(data);
datadiff(abs(datadiff)<1e-10) = 0;

% 極大値探索
ind = find(datadiff(2:end).*datadiff(1:end-1)<0 & 0<datadiff(1:end-1)) + 1;
ind = [1, ind, length(data)];

% 以下グラフ作成 ============
scrsz = get(0, 'ScreenSize');
figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 scrsz(4)-620 480 480]);

% 軸作成
set(0,'defaultAxesFontSize',8);
set(0,'defaultAxesFontName','Arial');
set(0,'defaultTextFontSize',10);
set(0,'defaultTextFontName','Arial');
axes('Parent', gcf, 'PlotBoxAspectRatio', [1 1 1], 'Box', 'on');

% グラフ -----
subplot(3, 1, 1), plot(1:length(data), data, 'b-');
set(gca, 'Title', ...
    text('String', 'All Data Plot', ...
        'FontSize', 12, 'FontName', 'Times', ...
        'FontWeight', 'Bold', 'FontAngle', 'Italic'), ...
    'XlimMode', 'manual', 'Xlim', [0, length(data)], 'YlimMode', 'manual', 'Ylim', [min(data), max(data)*1.1]);
xlabel('Index');
ylabel('Value');

subplot(3, 1, 2), plot(1:length(data), data, 'b-', ind, data(ind), 'r-');
set(gca, 'Title', ...
    text('String', 'All Data & Local Maximal Values Plot', ...
        'FontSize', 12, 'FontName', 'Times', ...
        'FontWeight', 'Bold', 'FontAngle', 'Italic'), ...
    'XlimMode', 'manual', 'Xlim', [0, length(data)], 'YlimMode', 'manual', 'Ylim', [min(data), max(data)*1.1]);
xlabel('Index');
ylabel('Value');

subplot(3, 1, 3), plot(ind, data(ind), 'r-');
set(gca, 'Title', ...
    text('String', 'Local Maximal Values Plot', ...
        'FontSize', 12, 'FontName', 'Times', ...
        'FontWeight', 'Bold', 'FontAngle', 'Italic'), ...
    'XlimMode', 'manual', 'Xlim', [0, length(data)], 'YlimMode', 'manual', 'Ylim', [min(data), max(data)*1.1]);
xlabel('Index');
ylabel('Value');


diff (MATLAB Function Reference)
スポンサーサイト



このエントリに付けられたタグ|MATLAB極大値極小値グラフ

[その他] 適当はダメすね

「FigureのメタファイルをWord2000で正しく表示させる」の記事に書いてあった関数が良く考えるとパスに半角スペースあると動かないなと気づいたので修正した。
ファイルも差し替えました。

結構適当に書いてるからそういうの多そうだな。

このエントリに付けられたタグ|MATLABファイル操作グラフィックスFigureメタファイル

[Geometry] ベジェ曲線とバーンスタイン基底関数

Geometryってカテゴリに内積しか書いてないことを気づいたので,とりあえずベジェ曲線などを定義してみる。
ベジェ曲線の定義とかはWikipediaでも見てもらうとして,実装だけします。

まず,バーンスタイン基底関数(Bernstein Polynomial)の実装をします。
この関数はこんな感じで定義されます。
二項係数(Binomial Coefficient)に関しては面倒なので,MATLABのFILE EXCHANGEのbinomial.mをDLするのもいいかもしれません。
が,N次ベジェ曲線の場合はどうせ0~Nまで計算するので,コレ↓で十分だと思います。
diag(rot90(pascal(N+1)))

これでなんで二項係数がでてくるかはこれを読んでもらうとして,経過としてはこんな感じ↓。
NN+1のパスカル行列を作ると反対角成分に欲しい二項係数の配列が出るので,これをrot90で90度回転させた後diag()で対角成分を取ってます。

>NN=4;
>A = pascal(NN+1)
A =
     1     1     1     1     1
     1     2     3     4     5
     1     3     6    10    15
     1     4    10    20    35
     1     5    15    35    70
>B = rot90(A)
B =
     1     5    15    35    70
     1     4    10    20    35
     1     3     6    10    15
     1     2     3     4     5
     1     1     1     1     1
>diag(B)
ans =
     1
     4
     6
     4
     1

これを使えばこんな感じでバーンスタイン関数がかけます。
ones(1,length(TT))をかけたりしてるのはTTが配列で与えられたときに対応するものです。
function BE = bernstein(NN, TT)
onest = ones(1, length(TT));
II = [0:NN]'*onest;
BI = diag(rot90(pascal(NN+1)))*onest;
TT = ones(NN+1,1)*TT;

BE = (1 - TT).^(NN - II).*BI.*(TT.^II);

ここまでくれば後は簡単。定義よりベジェ曲線はこれで定義できます。
function BEZ = bezier(TT, CPS)
BEZ = CPS*bernstein(length(CPS(1,:))-1, TT);

以下,実行例です。
とりあえず,平面3次ベジェ曲線です。
> cps=[[0;0],[2;3],[3;4],[4;1]]
cps =
     0     2     3     4
     0     3     4     1
>cds=bezier([0:1/100:1],cps);
plot(cds(1,:),cds(2,:),'b-');
axis equal;
hold on;
plot(cps(1,:),cps(2,:),'r-o');

平面のベジェ曲線

で,続いて3次元空間の5次ベジェ曲線。
> cps=[[0;0;0],[2;3;-1],[3;4;1],[4;4;2],[5.6;2;1],[6;0;0]]
cps =
         0    2.0000    3.0000    4.0000    5.6000    6.0000
         0    3.0000    4.0000    4.0000    2.0000         0
         0   -1.0000    1.0000    2.0000    1.0000         0
> cds=bezier([0:1/100:1],cps);
plot3(cds(1,:),cds(2,:),cds(3,:),'b-');
axis equal;
box on;
hold on;
plot3(cps(1,:),cps(2,:),cps(3,:),'r-o');

空間のベジェ曲線

以上,ベジエ曲線の実装でした。
とりあえず僕が書いたコードでいい人は下記リンクからDLしてください。

pascal (MATLAB Function Reference)
rot90 (MATLAB Function Reference)
diag (MATLAB Function Reference)
このエントリのmファイル(bernstein.m, bezier.m)
このエントリに付けられたタグ|MATLABGeometry幾何ベジエ曲線ベジェ曲線Bezierバーンスタイン関数

[Geometry] 単位ベクトル

2007/12/03 他の形式のベクトルにも対応した修正版がこちらに↓
初心者によるMATLABメモ| 単位ベクトル2 (unitvector)


ベクトル配列を各列ごとに単位ベクトルに変換する方法.
まずこんな感じの5つの3次元ベクトルの配列があるとします.
> A=rand(3,5)
A =
    0.4608    0.4122    0.2974    0.6501    0.4001
    0.4574    0.9016    0.0492    0.9830    0.1988
    0.4507    0.0056    0.6932    0.5527    0.6252

このベクトル配列の各ベクトルの大きさはこれで求まります.
> VS=sqrt(sum(A.*A))
VS =
    0.7903    0.9914    0.7559    1.3017    0.7684

あとは各列をこの配列VSの各要素で割ってやればベクトルは単位ベクトルになります.
で,前処理としてVSをAと同じ次元にしてやります.
> VS = ones(length(A(:,1)),1)*VS
VS =
    0.7903    0.9914    0.7559    1.3017    0.7684
    0.7903    0.9914    0.7559    1.3017    0.7684
    0.7903    0.9914    0.7559    1.3017    0.7684

で,最後に基のベクトル配列AをこのVSで割る。これでもとのベクトルの配列Aの単位ベクトル配列ができます。
> A./VS
ans =
    0.5830    0.4158    0.3935    0.4994    0.5207
    0.5787    0.9094    0.0650    0.7552    0.2587
    0.5703    0.0056    0.9170    0.4246    0.8136

一応,関数化するとこんな感じかと思います。(ソース
function UV = unitvector(V)
% UNITVECTOR ベクトルの単位ベクトル化。
error(nargchk(1, 1, nargin));
if ~isreal(V) error('ベクトルは実数配列である必要があります。'); end;

VS = ones(length(V(:,1)),1)*sqrt(sum(V.*V));
UV = V./VS;

初心者によるMATLABメモ| 同じ行や列を持つ行列/各行(列)ごとに特定の値を四則演算する
初心者によるMATLABメモ| ベクトルの内積(ドット積)
このエントリのmファイルを含むファイル(unitvector.mの最新版)
このエントリに付けられたタグ|MATLABGeometry幾何単位ベクトルunitvector

[ファイル操作] 相対パスから絶対パスに変換する2

前にも書いたんですが,前に書いた方法だとよく考えたら指定したパスのファイルがないときに動かないなと思って,パート2を。
ファイルがないときでも動かないとsave系のコマンドで使うとき困るし。

まぁ基本的には変わんないですけど。
保存先に指定したフォルダ(ディレクトリ)が存在すればこれで指定したファイルが存在しなくても指定したファイル,ディレクトリの絶対パスが取得できます。
ちなみに,この関数は,指定したパスのファイルが存在しない場合でも動きますが,
指定したディレクトリ,もしくは指定したファイルのすぐ上位のディレクトリが存在しない場合はエラーとなります。
function APATH = absolute_path(RPATH) 
% ABSOLUTE_PATH 指定したファイル,ディレクトリの絶対パスを返す。

if ~ischar(RPATH) error('第1引数のパスは文字列である必要があります。'); end; 

[dirname, filename, ext, versn] = fileparts(RPATH);
[num,pathinfo] = fileattrib(dirname);
APATH = fullfile(pathinfo.Name, strcat([filename, ext]));

一応,実行例。
まずカレントディレクトリの絶対パスを取得してみる。
> absolute_path('./')
ans =
C:Program FilesMATLABR2006bwork

続いてファイルがない場合。まずひとつ上をdirしてみるとこんな感じ。test.txtはありません。
> dir ../
.                  demos              ja                 notebook           stateflow          uninstall          
..                 etc                java               patents.txt        sys                work               
MATLAB R2006b.lnk  extern             jhelp              rtw                toolbox            
bin                help               license.txt        simulink           trademarks.txt     

しかし,指定してもちゃんと絶対パスが得られます。
> absolute_path('../test.txt')
ans =
C:Program FilesMATLABR2006btest.txt

適当にエラー処理+ヘルプ付きのファイルが欲しい人は下記リンクにソースがあります。

fullfile (MATLAB Function Reference)
fileattrib (MATLAB Function Reference)
fileparts (MATLAB Function Reference)
初心者によるMATLABメモ| 相対パスから絶対パスに変換する
このエントリのmファイル(absolute_path.m)
このエントリに付けられたタグ|MATLABファイル操作相対パス絶対パスディレクトリフォルダ

[基本] 同じ行や列を持つ行列/各行(列)ごとに特定の値を四則演算する2

前に同じ行や列を持つ行列/各行(列)ごとに特定の値を四則演算するってエントリ書いたんですが,今更衝撃的なことを知りました。
ふとrepmatのドキュメントを見ていたら…

… repmat(A,m,n) は、Aがスカラのとき、Aの値を要素にもつ m行n列の行列を作ります。 これは、m または nが大きいとき、a*ones(m,n)よりも、かなり高速に処理されます。 …
repmat (MATLAB Function Reference)


あ…

a*ones(m,n)よりも、かなり高速に処理されます。

…かな~りの衝撃です。
いや、ちょっとショック。


なわけで,とりあえず前のエントリの例を全てrepmatで書き換えてここに記しておきます。
基本的にはこれです。
A = repmat(A,size(A')) = 1行ベクトルAをn行に増やす.
B = repmat(B,size(B')) = 1行ベクトルBをn列に増やす.
あとはrepmatのドキュメントを見てください。


で、実行例.
行ベクトルAを同じ行を持つn×n行列に.
>> A=[1:5]
>>repmat(A,size(A'))
ans =
     1     2     3     4     5
     1     2     3     4     5
     1     2     3     4     5
     1     2     3     4     5
     1     2     3     4     5

列ベクトルBを同じ列を持つn×n行列に.
>> B=[1:5]'
B =
     1
     2
     3
     4
     5
>> repmat(B,size(B'))
ans=
     1     1     1     1     1
     2     2     2     2     2
     3     3     3     3     3
     4     4     4     4     4
     5     5     5     5     5

j個要素を持つ列ベクトルCを同じ列を持つn×jの行列に.
>> C=rand(1,5)
C =
    0.7266    0.4120    0.7446    0.2679    0.4399
>> repmat(C,3,1)
ans=
    0.7266    0.4120    0.7446    0.2679    0.4399
    0.7266    0.4120    0.7446    0.2679    0.4399
    0.7266    0.4120    0.7446    0.2679    0.4399

Dの各列に同じ値の四則演算をする.(ここでは積)
>> D=[1:5;2:2:10;3:3:15]
D =
     1     2     3     4     5
     2     4     6     8    10
     3     6     9    12    15
>> E=[1:5]
E =
     1     2     3     4     5

そのまま要素ごとの積をするとエラーが出る。
>> D.*E
??? エラー: ==> times
行列の次元は同じである必要があります

で,こうする↓
>> E=repmat(E,length(D(:,1)),1)
E =
     1     2     3     4     5
     1     2     3     4     5
     1     2     3     4     5
>> D.*E
ans=
     1     4     9    16    25
     2     8    18    32    50
     3    12    27    48    75
>> 


しかし,自分が今まで書いた関数とか書き直すのはたるいので基本放置(w

repmat (MATLAB Function Reference)
このエントリに付けられたタグ|MATLAB行列repmat

[Geometry] 同次座標系における平行移動行列

同次座標系の変換行列のまず最初として,平行移動行列を。
ごく簡単に説明します。
↓にように点pがベクトルdによって平行移動され,点になるとします。
平行移動1


このときのベクトルdによる平行移動は行列で表すと下記になります。
平行移動2

これがなんでそうなるかは実際に計算してもらえば分かるかと思います。
で,実装。
ごくごく簡単に実装するならこんな感じ。
d=[1;2;3;0]
eye(4)+[zeros(4,3),d]

まぁしかしこれでは実際に使うとき,例によってMATLAB御法度のforループが出てきますから問題あり。
とりあえず,こんな感じにすれば平行移動ベクトルが複数個与えられても対応できると思います。
TM = repmat(eye(4),1,length(DV(1,:)));
TM(1:3,4:4:end) = DV(1:3,:);

この場合はTMは,平行移動行列を横に連結したものとして得られます。
ここら辺は趣味なのでわかりませんが,3次元配列として欲しいとかセル配列の方が…とかいろいろあるかと思いますが,そこは自分でお願いします。

とりあえず僕が書いた平行移動行列TRANSRATION_MAT(DV)を下記に。
これは同次座標系のDVへの平行移動の行列を返します。
DVがn個指定された場合,平行移動行列を横に連結したものとして返します。
引数DV列ベクトルは3x1でも4x1でも構いません。
function TM = translation_mat(DV)
% TRANSRATION_MAT : 同次座標系における平行移動行列

error(nargchk(1, 1, nargin));
if ~isreal(DV) || all(length(DV(:,1))~=[1, 3, 4]) error('平行移動ベクトルは,実数でスカラもしくは3次元ベクトル,同次座標系のベクトルである必要があります。'); end;
if length(DV(:,1))==1 DV = repmat(DV, 3, 1); end;

TM = repmat(eye(4),1,length(DV(1,:)));
TM(1:3,4:4:end) = DV(1:3,:);

で,例としてはこんな感じ。

X,Y,Z軸正の方向にそれぞれ1平行移動する行列
>> translation_mat(1)
ans =
     1     0     0     1
     0     1     0     1
     0     0     1     1
     0     0     0     1

X軸正の方向に1平行移動する行列
>> translation_mat([1;0;0])
ans =
     1     0     0     1
     0     1     0     0
     0     0     1     0
     0     0     0     1

この例は指定した各列ベクトルに関する平行行列を横に連結したものをAに代入しています。
>> A=translation_mat([[2;3;4],[2;1;3],[10;10;10]])
A =
     1     0     0     2     1     0     0     2     1     0     0    10
     0     1     0     3     0     1     0     1     0     1     0    10
     0     0     1     4     0     0     1     3     0     0     1    10
     0     0     0     1     0     0     0     1     0     0     0     1

実際に使う場合は,こんな感じで次元を変更して使うとか。
>> reshape(A,4,4,length(A(1,:))/4)
ans(:,:,1) =
     1     0     0     2
     0     1     0     3
     0     0     1     4
     0     0     0     1
ans(:,:,2) =
     1     0     0     2
     0     1     0     1
     0     0     1     3
     0     0     0     1
ans(:,:,3) =
     1     0     0    10
     0     1     0    10
     0     0     1    10
     0     0     0     1


こんな感じでブロック対角化するとかでいろいろです。
>> A2 = mat2cell(A,4,4*ones(1,length(A(1,:))/4));
>> A2 = blkdiag(A2{:})
A2 =
     1     0     0     2     0     0     0     0     0     0     0     0
     0     1     0     3     0     0     0     0     0     0     0     0
     0     0     1     4     0     0     0     0     0     0     0     0
     0     0     0     1     0     0     0     0     0     0     0     0
     0     0     0     0     1     0     0     2     0     0     0     0
     0     0     0     0     0     1     0     1     0     0     0     0
     0     0     0     0     0     0     1     3     0     0     0     0
     0     0     0     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     0     0     0     1     0     0    10
     0     0     0     0     0     0     0     0     0     1     0    10
     0     0     0     0     0     0     0     0     0     0     1    10
     0     0     0     0     0     0     0     0     0     0     0     1

一応平行移動する例とか書くと,こんな感じの点群Bがあるときにそれぞれをそれぞれ平行移動させるときには,
>> B=[rand(3);ones(1,3)]
B =
    0.1210    0.8928    0.8656
    0.4508    0.2731    0.2324
    0.7159    0.2548    0.8049
    1.0000    1.0000    1.0000

列ベクトル化して掛け合わせて,reshapeで整形すれば一挙に並行移動できます。
>> A2*[B(:)]
ans =
    2.1210
    3.4508
    4.7159
    1.0000
    2.8928
    1.2731
    3.2548
    1.0000
   10.8656
   10.2324
   10.8049
    1.0000

>> reshape(A2*[B(:)],4,length(A2(1,:))/4)
ans =
    2.1210    2.8928   10.8656
    3.4508    1.2731   10.2324
    4.7159    3.2548   10.8049
    1.0000    1.0000    1.0000

参考文献は下記のやつがそばにあったんで一応書きましたが,別に普通の日本語の本に載ってると思います。

参考文献:J. Foley, A. Van Dam, S. Feiner, J. Hughes, and R. Phillips:Introduction to Computer Graphics, Addison-Wesley, United Staes of America, (1993) 181-182

reshape (MATLAB Function Reference)
blkdiag (MATLAB Function Reference)
このエントリに付けられたタグ|MATLABGeometry幾何同次座標変換平行移動

[Geometry] 同次座標系における拡大、縮小行列(スケーリング)

前のエントリの平行移動行列に続いて,同次座標系の変換行列の拡大、縮小行列。
こちらも,まず簡単に説明しますと,
pがスケーリングベクトルs=[sx,sy,sz]によって拡大、縮小移動され,点になるとすると拡大、縮小行列は下記の式のように表されます。
スケーリング

で,実装していきます。
ごくごく簡単に実装するなら,スケーリングベクトルs=[sx,sy,sz]が単に対角に並んでいればいいんで,これでいい。
>> in=[1;2;3];
>> diag([in;1])
ans =
     1     0     0     0
     0     2     0     0
     0     0     3     0
     0     0     0     1

で、まぁこれでは例のごとく,実際に使うときMATLAB御法度のforループが出てきますから問題あり。
とりあえず,こんな感じにすれば拡大,縮小ベクトルが複数個与えられても対応できる。
function SM = scaling_mat(SV)
% SCALING_MAT : 3次元空間の4×4の拡大,縮小行列

error(nargchk(1, 1, nargin));
if ~isreal(SV) || all(length(SV(:,1))~=[1, 3, 4]) error('スケーリングベクトルは実数の3次元ベクトルである必要があります。'); end;
if length(SV(:,1))==1 SV = repmat(SV, 3, 1); end;

SM = repmat(eye(4),1,length(SV(1,:)));
SM(1,1:4:end) = SV(1,:);
SM(2,2:4:end) = SV(2,:);
SM(3,3:4:end) = SV(3,:);

この場合はSMは,スケーリング行列を横に連結したものとして得られます。
実際にこの関数を使ってみるとこんな感じです。Sで適当に縮小するためのベクトルを作ってます。
>> S=rand(3,5)
S =
    0.9458    0.6673    0.2441    0.3054    0.6317
    0.3404    0.8227    0.0234    0.6011    0.0129
    0.6055    0.3147    0.3797    0.9730    0.2249
>> SM = scaling_mat(S)
ans =
  Columns 1 through 14
    0.9458         0         0         0    0.6673         0         0         0    0.2441         0         0         0    0.3054         0
         0    0.3404         0         0         0    0.8227         0         0         0    0.0234         0         0         0    0.6011
         0         0    0.6055         0         0         0    0.3147         0         0         0    0.3797         0         0         0
         0         0         0    1.0000         0         0         0    1.0000         0         0         0    1.0000         0         0

  Columns 15 through 20
         0         0    0.6317         0         0         0
         0         0         0    0.0129         0         0
    0.9730         0         0         0    0.2249         0
         0    1.0000         0         0         0    1.0000

で、この関数は,スカラを与えた場合x,y,z全方向にそのスカラで拡大,縮小します。
>> scaling_mat(3)
ans =
     3     0     0     0
     0     3     0     0
     0     0     3     0
     0     0     0     1

もうあとは平行移動のときと一緒ですが、
まとめて計算するときは列ベクトル化して掛け合わせて,reshapeで整形すれば一挙に並行移動できます。
下記の例は[4;3;2]のベクトルを縮小してみています。
>> PP=repmat([4;3;2;1],1,4)
PP =
     4     4     4     4
     3     3     3     3
     2     2     2     2
     1     1     1     1
>> SM=mat2cell(SM, 4, 4*ones(1,length(SM(1,:))/4));
>> SM=blkdiag(SM{:});
>> SM*[PP(:)]
ans =
    2.8010
    1.3865
    0.3215
    1.0000
    2.3968
    0.7683
    0.1692
    1.0000
    0.6911
    2.8092
    1.5676
    1.0000
    0.5250
    1.0769
    1.3331
    1.0000

>> reshape(SM*[PP(:)],4,4)
ans =
    2.8010    2.3968    0.6911    0.5250
    1.3865    0.7683    2.8092    1.0769
    0.3215    0.1692    1.5676    1.3331
    1.0000    1.0000    1.0000    1.0000


参考文献:J. Foley, A. Van Dam, S. Feiner, J. Hughes, and R. Phillips:Introduction to Computer Graphics, Addison-Wesley, United Staes of America, (1993) 181-182

diag (MATLAB Function Reference)
reshape (MATLAB Function Reference)
blkdiag (MATLAB Function Reference)
このエントリに付けられたタグ|MATLABGeometry幾何同次座標変換拡大縮小スケーリング

[Geometry] 同次座標系の回転行列(ローテーション)

前のエントリの平行移動行列拡大、縮小行列に続いて,同次座標系の変換行列の回転行列。
まず任意の軸に対する回転行列の基本形は下記。
これは本を見てもらう方がよいと思います。
pが回転軸u=[ux,uy,uz]周りに回転角θで回転され,点になるとするときの回転行列Mは下記の式のように表されます。
回転行列001

で,実装していきます。
単なる1個だけ欲しいならそのまま書いてこれです。
cosR = cos(RAD);
sinR = sin(RAD);
RM = [[cosR + (1 - cosR)*AX(1).^2; AX(1)*AX(2)*(1 - cosR) + AX(3)*sinR; AX(3)*AX(1)*(1 - cosR) - AX(2)*sinR; 0], ...
[AX(1)*AX(2)*(1 - cosR) - AX(3)*sinR; cosR + AX(2)*AX(2).*(1 - cosR); AX(2)*AX(3)*(1 - cosR) + AX(1)*sinR; 0], ...
[AX(3)*AX(1)*(1 - cosR) + AX(2)*sinR; AX(2)*AX(3)*(1 - cosR) - AX(1)*sinR; cosR + AX(3)*AX(3).*(1 - cosR);0], ...
[0;0;0;1]];

で、まぁこれでは例のごとく,MATLAB御法度のforループが出てくる実用に耐えない表現だから問題あり。
実際いろいろな実装方法があるかと思いますが,ここではまず回転行列Mを列表現にした下記M'のような1行の列を横に連結するマトリックを考えます。
回転行列2

基本的な戦略としては,これを横に連結し後から整形します。
なんでこの立場を取るかというと一応回転行列の場合は3×3が欲しい場合もあるので,この方が整形しやすいのが主な理由です。
これを関数化すると下記。
function RM = rotation_mat0(AX, RAD)
% ROTATION_MAT0 : 同次座標系の回転変換行列を1列にした行列を返す


rcols = length(RAD);
[arows, acols] = size(AX);
mcols = max(rcols, acols);

error(nargchk(2, 2, nargin));
if ~isreal(AX) || all(arows~=[3,4]) error('軸ベクトルは3次元の実数配列である必要があります。'); end;
if ~isreal(RAD) error('回転角は実数配列である必要があります。'); end;
if acols~=1 && rcols~=1 && rcols~=acols error('軸が3x1の行列,回転角がスカラ,もしくは軸と回転角の列の数が等しくある必要があります。'); end;

AX = unitvector(AX);

cosR = cos(RAD);
sinR = sin(RAD);
zeroM = zeros(1, mcols);

RM = [cosR + (1 - cosR).*AX(1,:).^2;
AX(1,:).*AX(2,:).*(1 - cosR) + AX(3,:).*sinR;
AX(3,:).*AX(1,:).*(1 - cosR) - AX(2,:).*sinR;
zeroM;
AX(1,:).*AX(2,:).*(1 - cosR) - AX(3,:).*sinR;
cosR + AX(2,:).*AX(2,:).*(1 - cosR)
AX(2,:).*AX(3,:).*(1 - cosR) + AX(1,:).*sinR;
zeroM;
AX(3,:).*AX(1,:).*(1 - cosR) + AX(2,:).*sinR;
AX(2,:).*AX(3,:).*(1 - cosR) - AX(1,:).*sinR;
cosR + AX(3,:).*AX(3,:).*(1 - cosR);
zeroM;
repmat([0;0;0;1], 1, mcols)
];

で,後は整形するだけ。
上の関数を内部で呼び出して実装すると4×4ならこれ
function RM = rotation_mat4(AX, RAD)
% ROTATION_MAT4 : 同座標系の回転変換行列 

error(nargchk(2, 2, nargin));
RM = rotation_mat0(AX, RAD);
RM = reshape(RM, 4, 4*length(RM(1,:)));

3×3ならこれです。
function RM = rotation_mat3(AX, RAD)
% ROTATION_MAT3 : 3次元空間の3×3の回転変換行列(テンソル) 

error(nargchk(2, 2, nargin));
RM = rotation_mat0(AX, RAD);
RM([4:4:end, end-3:end-1],:) =[]; 
RM = reshape(RM, 3, 3*length(RM(1,:)));

この関数がちゃんと動くかは適当に自分で確認してもらえると幸いです。
一応例を示しておくとこんな感じ。
例のRMはz軸周りにπ/3回転,y軸周りにπ/2回転,軸(1,1,1)周りにpi/4回転する回転行列を横に連結したものです。
>> RM=rotation_mat4([[0;0;1],[0;1;0],[1;1;1]], [pi/3,pi/2,pi/4])
RM =
    0.5000   -0.8660         0         0    0.0000         0    1.0000         0    0.8047   -0.3106    0.5059         0
    0.8660    0.5000         0         0         0    1.0000         0         0    0.5059    0.8047   -0.3106         0
         0         0    1.0000         0   -1.0000         0    0.0000         0   -0.3106    0.5059    0.8047         0
         0         0         0    1.0000         0         0         0    1.0000         0         0         0    1.0000

でこれら3つのそれぞれの回転をx軸方向の単位ベクトルに適応したものが下記です。あってるかの確認は自分で…
>> VV = repmat([1;0;0;1],1,3);
>> RMcell = mat2cell(RM,4,4*ones(1,length(RM(1,:))/4));
>> RM = blkdiag(RMcell{:});
>> RV = reshape(RM*[VV(:)], 4, length(VV(1,:)))
RV =
    0.5000    0.0000    0.8047
    0.8660         0    0.5059
         0   -1.0000   -0.3106
    1.0000    1.0000    1.0000


参考文献:J. Foley, A. Van Dam, S. Feiner, J. Hughes, and R. Phillips:Introduction to Computer Graphics, Addison-Wesley, United Staes of America, (1993) p.p.192

reshape (MATLAB Function Reference)
blkdiag (MATLAB Function Reference)
このエントリに付けられたタグ|MATLABGeometry幾何同次座標変換回転変換ローテーション回転行列回転テンソル

[Graphics] 3次元のグラフィックスで矢印を表示する

3次元のグラフィックスのときにベクトルを上手いこと表示するいい関数がないかな。と思って探してたらFILE Exchangeにあるarrow3いい感じですね。
使い方はこんな感じ↓
ARROW3(P1,P2)
ARROW3(P1,P2,S,W,H,IP,ALPHA,BETA)

※P1,P2以外はオプション。
P1:始点座標。
P2:終点座標(こちらに矢印の頭がつく)。
S:PlotのSに似てる。'b-'とかで指定。help arrow3をした方がわかりやすい。
W:矢印の頭の太さ。
H:矢印の頭の長さ。
IP:始点に球を起きたい場合にスカラで半径を指定。
ALPHA:パッチ面の透明度 (FaceAlpha)。
BETA:影の部分と光が当たるところの相対的な明るさ。

適当に使用例。点[0,0,0]に単位フレームを書いてみるとこんな感じ。
AROOW3

で,そのコード。参考までに
pp=[[0,0,0];eye(3)];

WIN=[280,280]; % [横,縦]
scrsz = get(0,'ScreenSize');
figure('InvertHardcopy', 'off', 'Color', [0 0,0], 'Position', [50 scrsz(4)-WIN(2)-80 WIN(1) WIN(2)]);
axes('Parent', gcf, 'PlotBoxAspectRatio', [1 1 1], 'Color', [0, 0, 0], 'Position', [0.100 0.100 0.9 0.9], 'Visible','off');

hold('all');
axis equal;

arrow3(pp(1,:), pp(2,:), 'r-', 0.1, 0.2)
arrow3(pp(1,:), pp(3,:), 'g-', 0.1, 0.2)
arrow3(pp(1,:), pp(4,:), 'b-', 0.1, 0.2)
light('Position',[-10 -10 -10],'Style','local');
light('Position',[60,60,60]), lighting gouraud;

ダウンロードしてパスが通ってるところにおいて置くと便利だと思います。

MATLAB Central File Exchange - Arrow3 version 4.58
このエントリに付けられたタグ|MATLABGraphics矢印ベクトル

[Geometry] 同次座標系の幾何学的な変換

ここ3つくらいのエントリで同次座標系の幾何学的な変換シリーズとして平行移動行列拡大、縮小行列(スケーリング)回転行列(ローテーション)の関数を紹介したので、ここらで暫定的にまとめたファイルとかUPしました。(こっち

とりあえず各変換の説明はこのあたりを見てください。
初心者によるMATLABメモ| 同次座標系における平行移動行列
初心者によるMATLABメモ| 同次座標系における拡大、縮小行列(スケーリング)
初心者によるMATLABメモ| 同次座標系の回転行列(ローテーション)

DLするとこんなファイルが入ってます。
各mファイルには一応helpを書きましたが、かなり適当…
rotation_mat0.m - 同次座標系の回転変換行列を1列にした行列を返す
rotation_mat3.m - 3次元空間の3×3の回転変換行列
rotation_mat4.m - 同座標系の回転変換行列
scaling_mat.m - 同座標系の拡大,縮小行列
transformations_samples.m - サンプルプログラム(要:Arrow3()
translation_mat.m - 同次座標系の平行移動行列
vector_rotation.m - 点、ベクトルの回転変換
vector_scaling.m - 点、ベクトルの拡大、縮小変換
vector_translation.m - 点、ベクトルの平行移動変換

で,下記サンプルプログラム(transformations_samples.m)を実行した例です。
ちなみにこのサンプルプログラムを動かすにはこっちのエントリで紹介してるarrow3が必要なので、MATLAB File Exchangeよりダウンロードしてください。
まず平行移動のデモ。↓が表示されてFigureウインドウが開きます。図中で黄色の矢印方向に平行移動させてます。
>> transformations_samples(1)
初期フレームを原点に平行移動します。
平行移動マトリックス
    1.0000         0         0    0.5000
         0    1.0000         0    1.3000
         0         0    1.0000    0.8000
         0         0         0    1.0000

幾何学的な変換1

続いて拡大、縮小のデモ。↓が表示されてFigureウインドウが開きます。
>> transformations_samples(2)
初期フレームをx軸方向に1.5倍,y方向に2倍, z軸方向に1/2倍に拡大、縮小します。
拡大マトリックス
    1.5000         0         0         0
         0    2.0000         0         0
         0         0    0.5000         0
         0         0         0    1.0000

幾何学的な変換2

最後に回転変換のデモ。↓が表示されてFigureウインドウが開きます。
図中の黄色のベクトルが回転軸です。
>> transformations_samples(3)
初期フレームを軸[1;-1;1;0]周りにpi/4だけ回転します。
回転マトリックス
    0.8047   -0.5059   -0.3106         0
    0.3106    0.8047   -0.5059         0
    0.5059    0.3106    0.8047         0
         0         0         0    1.0000

幾何学的な変換3

なにせ使い方は,各エントリ。もしくはhelp コマンドで見てください。
helpテキトーに書いてますけど…

初心者によるMATLABメモ| 同次座標系における平行移動行列
初心者によるMATLABメモ| 同次座標系における拡大、縮小行列(スケーリング)
初心者によるMATLABメモ| 同次座標系の回転行列(ローテーション)
このエントリの変換行列mファイル
このエントリに付けられたタグ|MATLABGeometry幾何同次座標変換回転変換回転行列平行移動拡大縮小

[雑記,愚痴] Bloglinesで画像が表示されないときが多いので

よくわかんないけどFirefox+BloglinesでFC2ブログをRSSで読んでると,画像がたびたび表示されない。(されたりもする…)
ひょっとしてBODYには,ちょっとだけ書いて,基本的にはEXTENDED BODYの方に本文を書いた方がいいのかしら…
まぁ現状TOPページが見づらい状況なのでその方がいいのかもなぁ。
ただ,RSSで読めなくなる不便なブログと化しますけど…

うーん,過去ログ全部書き換えるのも自動でやってうまくいくのかなぁとか…
まぁ検討で…
急にFC2ブログの画像Bloglinesとか有名なRSSリーダからのリファラは許可しねーかなぁ…
このエントリに付けられたタグ|