FC2ブログ

初心者によるMATLABメモ

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

 

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

[Geometry] 点列から曲率、曲率半径を求める(FastCode)

plotやplot3で表示した曲線の点列から数値的に曲率を求める関数curvature_numeric()です。
先に書いておきますが,接線主法線従法線のときと同様に、データ間の距離があまりにも異なる場合にはこの方法だと誤差が大きくなります。
基本的にある程度等距離に点が並んでいることを前提にした関数です。

まず、使い方です。
kk = curvature_numeric(ppp) … 点列pppから各点における曲率kkを求めます。
kk = curvature_numeric(ppp,dim) … 点ベクトルの形式をdimで指定して,点列pppから各点における曲率kkを求めます。

※点列pppは2次元,もしくは3次元である必要があります。
※第2引数dimは,どの次元でひとつの点を表現するかの整数値です。
※2次元の場合は曲率は正負の符号付で返ります。

ここでは,コードは省いてこの関数における曲率の求め方について簡単に書いておきます。(ダウンロードはこちらから)
基本的には,求めたい曲率を持つ点とその両隣の点の作る三角形の外接円の半径を曲率半径とみなし、その逆数として曲率を計算しています。
曲率 : 外接円
上の図より外接円の直径は下記の式で表されます。ここでRは外接円の半径です。
曲率の式 1 (1)
式(1)より、曲率は下記の式(2)で表すことが出来ます。
曲率の式 2 (2)
まず平面の場合を考えます。
このとき、下記の式(3)が成り立ちます。detは行列式です。
曲率の式 3 (3)
式(2)、(3)より曲率は下記の式(4)と近似できることがわかります。
曲率の式 4 (4)
また、平面でないときには下記の式(5)のように近似しています。
曲率の式 5 (5)
以下、実際の使用例です。
ここでは、この記事の生成例のベジエ曲線の曲率をcurvature_numeric()を使用して、計算してみます。

一応、曲率の理論値は下記です。
こちらの曲率は定義式からちゃんと計算した値です。
曲率のグラフ(理論値)
続いて、curvature_numeric()を使用して計算した曲率をプロットしてみると、下記のようになります。
両者を比較してもらえばわかると思いますが、結構上手く近似されると思います。
曲率のグラフ
一応補足ですが、曲率半径が欲しい場合は下記の関数を使ってください。これは単にcurvature_numeric()の逆数を求めて、絶対値をとっているいるだけです。
rr = curvatureR_numeric(ppp) … 点列pppから各点における曲率半径rrを求めます。
rr = curvatureR_numeric(ppp,dim) … 点ベクトルの形式をdimで指定して,点列pppから各点における曲率半径rrを求めます。

※点列pppは2次元,もしくは3次元である必要があります。
※第2引数dimは,どの次元でひとつの点を表現するかの整数値です。

ヘルプ付きのmファイルは下記よりダウンロードしてください。他の接線などの関数も同封です。

参考:初心者によるMATLABメモ| 点列から接線ベクトルを求める(FastCode)
参考:初心者によるMATLABメモ| 点列から主法線ベクトルを求める(FastCode)
参考:初心者によるMATLABメモ| 点列から従法線ベクトルを求める(FastCode)
参考:初心者によるMATLABメモ| 通過点列を通るC2連続な3次ベジエ曲線の生成
ダウンロード:このエントリのmファイル 左クリックでお願いします。(リファラを送信しない場合は,アクセス拒否になってしまいます)

テーマ:MATLAB - ジャンル:コンピュータ

このエントリに付けられたタグ|MATLABGeometry幾何曲率曲率半径curvaturecurvature_numeric()Fast-code

[Geometry] 点列から主法線ベクトルを求める(FastCode)

この記事の接線、この記事の従法線に続いて,plotやplot3で表示した曲線の点列から数値的に主法線を求める関数です。

まず、使い方です。他の2つと同様です。
nvec = normal_numeric(ppp) … 点列pppから各点における主法線nvecを求めます。
nvec = normal_numeric(ppp,dim) … 点ベクトルの形式をdimで指定して,点列pppから各点における主法線nvecを求めます。

※点列pppは2次元,もしくは3次元である必要があります。
※第2引数dimは,どの次元でひとつの点を表現するかの整数値です。

この関数における従法線の求め方について簡単に書いておきますと、下記のように従法線と接線の外積から求めています。(ダウンロードはこちらから)
法線

一応、接線、主法線、従法線は、下記のような関係になっています。
Frenet-Frame

また、接線についてはtangent_numeric()を利用しているので、下記の式で求めています。
接線

同様に従法線についてもbinormal_numeric()を利用しているので、下記の式で求めています。


以下、使用例です。
数値の例は接線、従法線と同様なのでここでは省略して、実際に表示させてみた例を示します。
螺旋曲線を生成してある点での接線、主法線、従法線を表示しています。
この例の場合もarrow3についてはこの記事中に書いてあるように,標準の関数ではないので、必要に応じてMATLAB CentralのFile Exchangeからダウンロードしてください。(arrow3の使い方はこちらを参照してください。)
% 3次元の例 : 螺旋曲線の接線(赤)、主法線(青)、従法線(緑)
divn = 200;
theta = linspace(0, 5*pi, divn);
theta2 = linspace(0, pi/3, divn);

ppp = 2*[cos(theta).*cos(theta2); sin(theta).*cos(theta2); sin(theta2)];
tvec = tangent_numeric(ppp);
nvec = normal_numeric(ppp);
bvec = binormal_numeric(ppp);

figure( 'InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [2 300 400 300]);
xlabel({'x'}, 'FontName','Times New Roman');
ylabel({'y'}, 'FontName','Times New Roman');
zlabel({'z'}, 'FontName','Times New Roman');
hold on;
nn=140;
plot3(ppp(1,:),ppp(2,:),ppp(3,:));
axis equal;
arrow3(ppp(:,nn)',[ppp(:,nn)+tvec(:,nn)]', 'r-', 0.1, 0.2, 0.1);
arrow3(ppp(:,nn)',[ppp(:,nn)+nvec(:,nn)]', 'b-', 0.1, 0.2, 0.1);
arrow3(ppp(:,nn)',[ppp(:,nn)+bvec(:,nn)]', 'g-', 0.1, 0.2, 0.1);
view(-18,26);



ヘルプ付きのmファイルは下記よりダウンロードしてください。

参考:MATLAB Central File Exchange - Arrow3 version 4.58
参考:初心者によるMATLABメモ| 3次元のグラフィックスで矢印を表示する
参考:初心者によるMATLABメモ| 点列から接線ベクトルを求める(FastCode)
参考:初心者によるMATLABメモ| 点列から従法線ベクトルを求める(FastCode)
ダウンロード:このエントリのmファイル 左クリックでお願いします。(リファラを送信しない場合は,アクセス拒否になってしまいます)

テーマ:MATLAB - ジャンル:コンピュータ

このエントリに付けられたタグ|MATLABGeometry幾何法線主法線主法線ベクトルnormalnormal_numeric()Fast-codeフレネ標構

[Geometry] 点列から従法線ベクトルを求める(FastCode)

この記事の接線に続いて,plotやplot3で表示した曲線の点列から数値的に従法線を求める関数です。
先に書いておきますが,データ間の距離があまりにも異なる場合にはこの方法だと誤差が大きくなります。
基本的にある程度等距離に点が並んでいることを前提にした関数です。

まず、使い方です。
bvec = binormal_numeric(ppp) … 点列pppから各点における従法線bvecを求めます。
bvec = binormal_numeric(ppp,dim) … 点ベクトルの形式をdimで指定して,点列pppから各点における従法線bvecを求めます。

※点列pppは2次元,もしくは3次元である必要があります。
※第2引数dimは,どの次元でひとつの点を表現するかの整数値です。

ここでは,コードは省いてこの関数における従法線の求め方について簡単に書いておきます。(ダウンロードはこちらから)
両端の点以外では,下記のように前後の点を用いて数値的な従法線を求めています。


この式は簡単に言うと前後の点を含めた3点の作る三角形の法線になっています。


また、始点、終点の両端点では下記のように単に隣の従法線を従法線としています。


以下、使用例です。
まず1点を列ベクトル(縦)として表した場合の例です。この例では2次元の点列から従法線を求めています。
>> ppp=[[0;0],[1;1],[3;2]]
ppp =
     0     1     3
     0     1     2

>> binormal_numeric(ppp)
ans =
     0     0     0
     0     0     0
     1     1     1

続いて、1点を行ベクトル(横)として表した場合の例です。この場合にはbinormal_numericに第2引数として次元2を指定してやる必要があります。(一応、指定しないでも上手くいく場合もありますが、明示的に指定してください。)また、この例も同様に2次元の点列から従法線を求めています。
>> ppp=[[0,0];[1,1];[3,2]]
ppp =
     0     0
     1     1
     3     2

>> binormal_numeric(ppp,2)
ans =
     0     0     1
     0     0     1
     0     0     1

この数字の例だけではわかりづらいので,実際に表示させてみた例を示します。
この記事と同様に、螺旋曲線を生成してある点での接線と従法線を表示しています。
この例の場合もarrow3についてはこの記事中に書いてあるように,標準の関数ではないので、必要に応じてMATLAB CentralのFile Exchangeからダウンロードしてください。(arrow3の使い方はこちらを参照してください。)
% 3次元の例 : 螺旋曲線の接線(赤)と従法線(緑)
divn = 200;
theta = linspace(0, 5*pi, divn);
theta2 = linspace(0, pi/3, divn);

ppp = 2*[cos(theta).*cos(theta2); sin(theta).*cos(theta2); sin(theta2)];
tvec = tangent_numeric(ppp);
bvec = binormal_numeric(ppp);

figure( 'InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [2 300 400 300]);
xlabel({'x'}, 'FontName','Times New Roman');
ylabel({'y'}, 'FontName','Times New Roman');
zlabel({'z'}, 'FontName','Times New Roman');
hold on;
nn=140;
plot3(ppp(1,:),ppp(2,:),ppp(3,:));
axis equal;
arrow3(ppp(:,nn)',[ppp(:,nn)+tvec(:,nn)]', 'r-', 0.1, 0.2, 0.1);
arrow3(ppp(:,nn)',[ppp(:,nn)+bvec(:,nn)]', 'g-', 0.1, 0.2, 0.1);
view(-18,26);



ヘルプ付きのmファイルは下記よりダウンロードしてください。

参考:MATLAB Central File Exchange - Arrow3 version 4.58
参考:初心者によるMATLABメモ| 3次元のグラフィックスで矢印を表示する
参考:初心者によるMATLABメモ| 点列から接線ベクトルを求める(FastCode)
ダウンロード:このエントリのmファイル 左クリックでお願いします。(リファラを送信しない場合は,アクセス拒否になってしまいます)

テーマ:MATLAB - ジャンル:コンピュータ

このエントリに付けられたタグ|MATLABGeometry幾何従法線従法線ベクトルbinormalbinormal_numeric()Fast-codeフレネ標構

[Geometry] ベジェ曲面(Beizer Surface)

なんとなく書いたBeizer曲面の関数bezierSurf()です。
ベジエ曲面自体については、このあたりを参考にしてください。
また、コードは綺麗でない、というかかなりわかりづらいので省略します。
興味のある人はダウンロードしたmファイルを覗いてみてください。(ダウンロードはこちらから)
簡単な使い方,使用例の順に示します。

まずは、ベジェ曲面の関数bezierSurf()の簡単な使用方法です。
2つのパラメータと制御メッシュ(Control Points)が引数となっています。
XYZ = bezierSurf(UU, VV, CPS) … CPSを制御メッシュとするN次ベジエ曲面のパラメータ(U,V)における座標を返す.
XYZ = bezierSurf(UU, VV, CPSX, CPSY, CPSZ)  … 制御点のx座標をCPSX,y座標をCPSY,z座標をCPSZとした制御メッシュのN次ベジエ曲面のパラメータ(U,V)における座標を返す.
[X,Y,Z] = bezierSurf(UU, VV, CPS) … CPSを制御メッシュとするN次ベジエ曲面のパラメータ(U,V)における座標x,y,zごとに返す.

パラメータUU, VVは,0以上1以下の実数配列で与える必要がある.
制御メッシュ点列CPSは,1つの3次元の実数行列で与える必要がある.
制御メッシュ点列CPSX, CPSY, CPSZは,3つの同じ次元を持つ2次元の実数行列で与える必要がある.

使用例として、まず単にxy平面となるような双5次ベジェ曲面の制御点を指定して,u=v=0.5としたときを計算した例です。
見て分かるように基本的には3次元行列として座標は出力されます。
(:,:,1)がx座標、(:,:,2)がy座標、(:,:,3)がz座標です。
>> [x,y] = meshgrid([0:5]);
>> z = zeros(size(x));
>> bezierSurf(1/2,1/2,x,y,z)

ans(:,:,1) =
    2.5000

ans(:,:,2) =
    2.5000

ans(:,:,3) =
     0

違う例として、まず単にz座標を乱数とした制御メッシュを持つ双3次ベジェ曲面の座標です。
>> [x,y] = meshgrid([0:3]);
z = rand(size(x));
cps = zeros(4,4,3);
cps(:,:,1) = x;
cps(:,:,2) = y;
cps(:,:,3) = z;
bezierSurf([0, 0.3, 1], [0,0,0], cps)

ans(:,:,1) =
         0    0.9000    3.0000


ans(:,:,2) =
     0     0     0


ans(:,:,3) =
    0.2769    0.4894    0.1869

また、出力を3つにすることでx,y,z座標をばらばらに取得することが可能です。
>> [bx,by,bz] = bezierSurf([0, 0.3, 1], [0,0,0], cps)

bx =
         0    0.9000    3.0000

by =
     0     0     0

bz =
    0.2769    0.4894    0.1869

最後に関数bezierSurf()を用いてベジエ曲面を表示してみるとこんな感じです。
同時に制御メッシュも表示してみました。
% 制御点
nn=6;
pts = [[0,1,2,3,4,5,0,0,0,0,0,0,0.97358,1.6353,0.70145,0.41548,0.45184,0.86041];
[0,1,2,3,4,5,1,1,1,1,1,1,0.87172,1.5897,1.878,0.60249,0.34142,0.36963];
[0,1,2,3,4,5,2,2,2,2,2,2,0.89357,1.2886,1.7519,0.94185,0.45533,1.8098];
[0,1,2,3,4,5,3,3,3,3,3,3,0.6127,0.75722,1.1003,0.46098,0.8714,1.9595];
[0,1,2,3,4,5,4,4,4,4,4,4,1.017,1.6232,1.245,1.6886,0.6222,0.87774];
[0,1,2,3,4,5,5,5,5,5,5,5,1.0215,1.0657,1.1741,0.38953,1.8468,0.22224]];

cps = zeros(nn,nn,3);
cps(:,:,1) = pts(:,1:6);
cps(:,:,2) = pts(:,7:12);
cps(:,:,3) = pts(:,13:18);

% グラフィック表示
figure( 'Name','Bezier Surface', 'Color', [1,1,1],  'Position', [100,200,600,300]);

hold on;
[uu,vv] = meshgrid(linspace(0,1,30));
bez1 = bezierSurf(uu,vv,cps);
ax1 = subplot(1,2,1), surf(bez1(:,:,1),bez1(:,:,2),bez1(:,:,3), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting','phong');
set(ax1, 'Position', [0.05, 0.05, 0.45, 0.9]);
axis equal;
view(-17,40);
title( 'Bezier Surface', 'FontSize', 12, 'FontWeight', 'bold', 'Color', [0,0,0], 'FontName', 'Arial');

ax2 = subplot(1,2,2), mesh(cps(:,:,1), cps(:,:,2), cps(:,:,3), 1.6*ones(nn), 'Marker', 'o', 'FaceColor', 'none');
set(ax2, 'Position', [0.5, 0.05, 0.45, 0.9]);
axis equal;
view(-17,40);
title( 'Control Mesh', 'FontSize', 12, 'FontWeight', 'bold', 'Color', [0,0,0], 'FontName', 'Arial');

caxis([0,2]);

bezierSurf


この関数beizerSurf()のダウンロードは↓から左クリックでお願いします。

参考:Bezier Surface
参考:初心者によるMATLABメモ| ベジェ曲線とバーンスタイン基底関数
ダウンロード:このエントリのmファイル 左クリックでお願いします。 (リファラを送信しない場合は,アクセス拒否になってしまいます)

テーマ:MATLAB - ジャンル:コンピュータ

このエントリに付けられたタグ|MATLABGeometry幾何曲面ベジエ曲面ベジェ曲面BezierSurfaceバーンスタイン関数

[ファイル操作] ディレクトリ(フォルダ)選択ダイアログを使用する

ディレクトリ(フォルダ)選択ダイアログを使用する
GUIを使って読み込むファイルを指定する方法をこちらの記事で、書き込むファイルを指定する方法をこちらの記事で紹介しました。
ここでは今度はディレクトリ(フォルダ)を選択する方法です。ディレクトリに読み込む、書き込むという概念はもちろんないので、どちらの場合でも使えます。

ディレクトリを選択する場合には、uigetdirを使います。基本な使い方は下記です。引数に関してはどちらも省略可能です。
directory_path = uigetdir('表示されたときに選択されているディレクトリ','ダイアログのコメント')

※directory_path … 選択したディレクトリのパス(失敗もしくは未選択、キャンセルなら0が返ります)

以下、実際に使ってみます。
例えば引数を省略して下記とすると,図のようなダイアログが表示され、ディレクトリを選択することが出来ます。返り値として選択したディレクトリのパスが返ります。
ディレクトリ(フォルダ)選択ダイアログ自体に関しては,特に説明は要らないと思うので省略します。
>> dpath = uigetdir
dpath =
C:¥sach1o¥MATLAB

uigetdir002

また,第1引数にパスを与えるとダイアログが表示されたときに、指定したディレクトリが選択された状態になります。
>> dpath = uigetdir('C:¥sach1o¥MATLAB')

また、第2引数は「ダイアログのコメント」と書きましたが、これはメニューバーのタイトルではなくて、下の画像の赤い線で囲んだ部分のコメントになります。これはOSによっても変わるかと思います。
>> dpath = uigetdir('','プロジェクトフォルダの指定')

uigetdir001

実際に使うなら、こんな感じ↓でしょうか?
これは簡単すぎますが、この例では最初にディレクトリを指定し、今回の実行に関するファイルは全て指定したのディレクトリに保存しています。
それとif文の部分でディレクトリがちゃんと指定されたかをチェックしています。
dpath = uigetdir('','プロジェクトフォルダの指定');

if dpath==0 
    error('プロジェクトフォルダが選択されていません。');
else 
    A = rand(3);
    B = rand(3);
    save(fullfile(dpath,'dataA.mat'), 'A');
    save(fullfile(dpath,'dataB.mat'), 'B');
    save(fullfile(dpath,'dataAll.mat'));
end;

実際に指定したディレクトリをdirしてみると下記のようにちゃんと指定したディレクトリにmatファイルが保存されたことが確認できます。
>> dir(dpath)
.            ..           dataA.mat    dataAll.mat  dataB.mat

参考:MATLAB:マニュアル:MATLAB:関数:サイバネットシステム-MATLAB Function Reference - uigetdir
参考:初心者によるMATLABメモ| ファイルダイアログから読み込むファイルを指定する
参考:初心者によるMATLABメモ| ファイルダイアログを使って保存するファイル名を指定する

テーマ:MATLAB - ジャンル:コンピュータ


前のページ 次のページ