FC2ブログ

初心者によるMATLABメモ

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

 

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

[その他] MATLABライセンスサーバのバージョンアップ

MATLABライセンスサーバのバージョンアップ作業のまとめです。
以下は自分の環境の場合で、ライセンスサーバがWindows XPのときの場合です。
作業はAdministratorでログインして行う。
下記が作業の流れ。
  1. 前のバージョンのMATLABを削除。(クライアントは残したままでも上手くいくかもしれないけど一応。FLEXlmは削除しない方がいいと思う。)
  2. PLPコードとライセンスファイルの取得(ここから
  3. 新しいバージョンのMATLABのインストール
  4. 起動するか確認
  5. Windowsの再起動
  6. 「Matlab Licence Server」サービスが起動しているか確認
  7. サービスの依存関係の変更(必要な場合)
  8. クライアントマシン用のライセンスファイルを作る。
  9. クライアントマシンにMATLABをインストールし、起動するか確認。

以下、もう少し詳細に書いてみたものです。

1. 前のバージョンのMATLABを削除。
Windowsの通常通り「プログラムの追加と削除」から削除

2. PLPコードとライセンスファイルの取得
The MathWorksのサイトにログイン。My MathWorks Accountの「Get Passcode」をクリック。
するとずらずらテキストが表示されるので,このような部分を探す。この数字文字列がPLPコード。
これをテキストファイルに保存しておく。
2) Install the software for Network Concurrent User installations

Use this Personal License Password and the License File below for
Windows Network Concurrent User installations:

18-99999-99999-99999-99999-99999-99999-99999-99999-99999-99999-99999-99999

ライセンスファイルはもう少し下の下記の「# BEGIN------cut here------CUT HERE-----BEGIN」の行から「# END------cut here------CUT HERE-----END」までの部分。
この部分をコピー&ペーストしてlicence.datとして保存する。
Use this License file for Windows, UNIX, Linux, and Mac Network
Concurrent User installations:

# BEGIN------cut here------CUT HERE-----BEGIN


省略


# END--------cut here-----CUT HERE-------END

3. 新しいバージョンのMATLABのインストール
これについては下記参照のこと。
参考MATLAB:技術サポート:インストールウィザード:フローティングライセンス:サーバーインストール (Windows):サイバネットシステム

4. 起動するか確認
インストールが完了したらこの段階で,一応MATLABを立ち上げ起動するかどうか確認。
起動しない場合,
  • FLEXlmの「Matlab Licence Server」のライセンスファイルに対するパスなどが合っているか、チェック。(前のバージョンのままになっていないか、どうか?)

5. Windowsの再起動
(略)。普通に再起動してください。

6. 「Matlab Licence Server」サービスが起動しているか確認
「プログラム→管理ツール→サービス」より「MATLAB License Server」が開始状態になっているか、確認。
起動しないなら,ファイヤーウォールの設定あたりが怪しい。
起動するなら7.の作業をやる。

7. サービスの依存関係の変更(必要な場合)
Windows起動時に他のサービスとMATLABのライセンスサーバ管理サービスが衝突してMATLABのライセンスサーバが正常に起動していない場合,衝突しているサービスをMATLABのライセンスサーバ管理サービスに依存するように変更する。
まず,regeditでレジストリエディタを起動し,下記のキーを一度削除。(MATLAB License Server以外と関連付けられている場合もあるので、その場合もそのサービスについても処理する必要がある。)
HKEY_LOCAL_MACHINE¥SYSTEM¥CurrentControlSet¥Services¥衝突しているサービス¥DependOnService

再び新規のキー「DependOnService」を複数行文字列として作成し、値を「MATLAB License Server」とする。
HKEY_LOCAL_MACHINE¥SYSTEM¥CurrentControlSet¥Services¥衝突しているサービス¥DependOnService=MATLAB License Server

Windowsを再起動して,再び「6.」をチェック。

8. クライアントマシン用のライセンスファイルを作る。
クライアントマシン用のライセンスファイル(license.dat)を作成する。
ここでのlicense.datはライセンスサーバのインストールの時に用いたものではなく,下記の2行だけのテキストファイルのこと。
1行目については、ライセンスサーバのライセンスファイルの1行目を参照すればよいかと思います。ちなみに「(」や「)」はいりません。
SERVER (ライセンスサーバのプライベートIP) (ライセンスサーバのPhysical Address) (ポート番号)
USE_SERVER

9. クライアントマシンにMATLABをインストールし、起動するか確認。
クライアントマシンにMATLABをインストールし,起動するか確認する。
参考MATLAB:技術サポート:インストールウィザード:フローティングライセンス:クライアントインストール (Windows):サイバネットシステム

だいたい、以上のような流れです。
スポンサーサイト



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

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

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

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

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


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

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

>> tangent_numeric(ppp)
ans =
    0.7071    0.8321    0.8944
    0.7071    0.5547    0.4472

続いて、1点を行ベクトル(横)として表した場合の例です。この場合にはtangent_numericに第2引数として2を指定してやる必要があります。また、この例も同様に2次元の点列から接線を求めています。
>> ppp=[[0,0];[1,1];[3,2]]
ppp =
     0     0
     1     1
     3     2


>> tangent_numeric(ppp,2)
ans =
    0.7071    0.7071
    0.8321    0.5547
    0.8944    0.4472

この数字の例だけではわかりづらいので,実際に表示させてみたいと思います。
まず、平面(2次元)の場合の例。ここでは円弧を書いてさらに途中の点での接線を表示させています。
また、ここで使用されている関数arrow3は標準の関数ではないので、必要に応じてMATLAB CentralのFile Exchangeからダウンロードしてください。(arrow3の使い方はこちらを参照してください。)
% 2次元の例 : 円の接線
theta = linspace(0,pi,30);
ppp = 2*[cos(theta);sin(theta)];
tvec = tangent_numeric(ppp);

figure( 'InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [2 300 400 240]);
xlabel({'x'}, 'FontName','Times New Roman');
ylabel({'y'}, 'FontName','Times New Roman');
hold on;
tvec(:,10)
plot(ppp(1,:),ppp(2,:));
axis equal;
arrow3(ppp(:,10)',[ppp(:,10)+tvec(:,10)]', 'r-', 0.1, 0.2, 0.1);



続いて、3次元の例。螺旋曲線を生成してある点での接線を表示しています。
この例の場合も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);

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);
view(-18,26);



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

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

[ファイル操作] ファイルダイアログから読み込むファイルを指定する

ファイルを読み込むときにファイルダイアログを表示させてGUIからファイルを指定したい場合もあるかと思います。
そんなときに使う関数uigetfileを紹介しておきます。
まずは、基本的な使い方。
[filename, dirpath, filterindex] = uigetfile( fext, dtitle, defaultname);

[引数]
fext … ファイルの拡張子のセル配列
dtitle … ダイアログのタイトル
defaultname … ファイル名ならデフォルトのファイル名が補完されます。また、絶対パスで指定した場合は最初に指定したディレクトリが表示されます。その場合単にディレクトリを指定する場合には最後にWindowsなら「¥」やLinuxなら「/」を付けてください。

[出力]
filename … ファイルの名前
pathname … ディレクトリのパス
filterindex … ファイルのパスの取得に成功したかどうか(失敗もしくは未選択なら0)

以下、実際に使ってみます。
例えば引数を与えず下記とすると…
[fname, dpath] = uigetfile;

こんな感じのダイアログが表示されます。


uigetfileを使う場合は、当然その後ファイルを読み込むと思いますが、そのためにはこの場合下記のどちらかの方法で、ファイルのパスを完全なものにしてやる必要があります。
この例の場合にはfpathに指定したファイルの正しいパスが入ります。
意味合い的にはfullfileを使った方がよいとは思いますが。
fpath = fullfile(dpath, fname)
もしくは
fpath = strcat([dpath, fname])

別に直接下記でロードしても構いません。
load(fullfile(dpath, fname))

ただ、実際にはファイルが指定されたかをチェックするべきです。
そんなときは、下記のような方法でチェックできます。
[fname, dpath] = uigetfile('*.mat');
if ~isequal(fname,0) & ~isequal(dpath,0)
    fpath = fullfile(dpath, fname);
    load(fpath);
else 
    disp('ファイルがありません。');
end;

出力を3つにすれば第3出力を使って下記の方法でチェックすることも可能です。
[fname, dpath, findex] = uigetfile('*.mat');
if findex~=0 
    fpath = fullfile(dpath, fname);
    load(fpath);
else 
    disp('ファイルがありません。');
end;

また、ファイルの種類を限定したい場合には、下記のようにワイルドカードを使った表現とその説明を組にしたセル配列で指定することが出来ます。
[fname, dpath, findex] = uigetfile( ...
    {'*.m;*.fig;*.mat;*.mdl', 'すべてのMATLABファイル(*.m, *.fig, *.mat, *.mdl)';
    '*.m',  'M-ファイル(*.m)'; ...
    '*.fig','Figures (*.fig)'; ...
    '*.mat','MAT-ファイル (*.mat)'; ...
    '*.mdl','モデルファイル(*.mdl)'; ...
    '*.*',  'すべてのファイル(*.*)'} ...
    , 'MATLABファイルの選択');

ファイルの説明を省略することも可能ですが、その場合には拡張子がMATLABファイルでないとダイアログ中で説明が省略されてしまいます。(画像参照)
[fname, dpath] = uigetfile({'*.m';'*.txt';'*.xls';'*.*'}, 'ファイルの選択');



また、第3の引数にファイル名を与えることで、デフォルトの名前を指定することができます。
[fname, dpath, findex] = uigetfile('*.*', 'ファイルの選択','aaaa.m');



また、絶対パスを第3の引数に与えると、ダイアログが表示されたときに最初に表示されるディレクトリを指定することができます。
気をつけるべきことは、ディレクトリを指定する場合には最後にWindowsなら「¥」やLinuxなら「/」を付ける必要があることです。
[fname, dpath, findex] = uigetfile('*.*', 'ファイルの選択','C:¥icons¥');



関数uigetfileを使って、ファイルダイアログ内で複数のファイル選択したい場合にはこちらの記事を参考にしてください。
また、ファイルを保存するときにファイルダイアログを表示させたいとき場合はこちらの記事を見てください。

参考:MATLAB Function Reference - uigetfile
参考:初心者によるMATLABメモ| ファイルダイアログを使って保存するファイル名を指定する
参考:初心者によるMATLABメモ| ディレクトリ(フォルダ)選択ダイアログを使用する

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

このエントリに付けられたタグ|MATLABファイル操作ファイルダイアログGUI読み込みuigetfile()

[ファイル操作] ファイルダイアログでの複数ファイルの選択

そういえば先程の記事で書かなかったのですが、関数uigetfileを使ったときにファイルダイアログ内で複数のファイル選択することもできます。
使い方は簡単で、uigetfileの引数に'MultiSelect', 'on'を与えるだけです。(デフォルト時は'off'になっています。)
[fname, dpath, findex] = uigetfile( ...
{'*.mat','MAT-files (*.mat)'; ...
'*.mdl','Models (*.mdl)'; ...
'*.*',  'All Files (*.*)'}, ...
'ファイルの選択', ...
'MultiSelect', 'on');

実際に'MultiSelect', 'on'とすると下記のように選択したファイル名がリスト化され、セル配列として取得できます。
>> [fname, dpath, findex] = uigetfile( '*.*', 'ファイルの選択', 'MultiSelect', 'on')

fname = 
  Columns 1 through 3
    'uigetfile001.png'    'uigetfile002.png'    'uigetfile003.png'

  Column 4
    'uigetfile004.png'

dpath =
C:¥sach1o¥MATLAB¥uidemo¥

findex =
     1



実際に処理するときには、下記のような感じでforループで回してあげればよいかと思います。
if findex==0 
    disp( 'ファイルが選択されていません。' );
else 
    for ind = 1:length(fname) 
        fpath = fullfile(dpath, fname{ind});
        disp([ '指定したファイルの絶対パスは', fpath, 'です。' ]);
    end;
end;

一応上記を実行した結果はこれです。
指定したファイルの絶対パスはC:¥sach1o¥MATLAB¥uidemo¥uigetfile001.pngです。
指定したファイルの絶対パスはC:¥sach1o¥MATLAB¥uidemo¥uigetfile002.pngです。
指定したファイルの絶対パスはC:¥sach1o¥MATLAB¥uidemo¥uigetfile003.pngです。
指定したファイルの絶対パスはC:¥sach1o¥MATLAB¥uidemo¥uigetfile004.pngです。


参考:MATLAB Function Reference - uigetfile
参考:初心者によるMATLABメモ| ファイルダイアログから読み込むファイルを指定する

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

[ファイル操作] ファイルダイアログを使って保存するファイル名を指定する

ファイルを読み込むときにファイルダイアログを表示させてGUIからファイルを指定したい場合に使う関数uigetfileはこちらの記事で紹介しました。
続いて今度はファイルを保存するときにファイルダイアログを表示させたいときに使う関数uiputfileについてです。

基本的な使い方はuigetfileと同じです。違うところは複数選択できないところでしょうか。
[filename, dirpath, filterindex] = uiputfile(fext, dtitle, defaultname);

[引数]
fext … ファイルの拡張子のセル配列
dtitle … ダイアログのタイトル
defaultname … ファイル名ならデフォルトのファイル名が補完されます。また、パスで指定した場合は最初に指定したディレクトリが表示されます。その場合単にディレクトリを指定する場合には最後にWindowsなら「\」やLinuxなら「/」を付けてください。

[出力]
filename … ファイルの名前(失敗もしくは未選択なら0)
pathname … ディレクトリのパス(失敗もしくは未選択なら0)
filterindex … ファイルのパスの取得に成功したかどうか(失敗もしくは未選択なら0)

以下、簡単に例を載せますが、使い方についてはuigetfileについての説明の記事も参考にしてください。
シンプルな一例として下記のコードを載せておきます。
この場合はファイルダイアログを表示させて入力した名前のMATファイルにワークスペースのデータを保存します。ダイアログがキャンセルされた場合は保存されません。
[fname, pname] = uiputfile('*.mat');
if ~isequal(fname,0) & ~isequal(dpath,0)
    fpath = fullfile(dpath, fname);
    save(fpath);
else 
    disp('保存しませんでした。');
end;

uigetfileのときと同様に、出力を3つにすれば第3出力を使って下記の方法でチェックすることも可能です。
[fname, dpath, findex] = uiputfile('*.mat');
if findex~=0 
        fpath = fullfile(dpath, fname);
    save(fpath);
else 
    disp('保存しませんでした。');
end;

さらにファイルダイアログのタイトルを自分で決めたい場合や
[fname, pname, findex] = uiputfile('', '名前をつけて保存');

デフォルトの名前を決めておきたい場合もuigetfileのときと同様です。
[fname, pname, findex] = uiputfile('','名前をつけて保存','default.m')

で、結局uigetfileとuiputfileは、別にそれ自体が保存や読み込みをしてくれるわけではなくパスを取得(もしくは生成)するだけです。
では、uigetfileとuiputfileのどっちを使ってもいいのでは?と思う方もいるかもしれませんが、なにが違うのかというと、まずすぐわかることは、下記のような見た目の違いが挙げられます。

・ファイルダイアログのタイトルのデフォルトがuigetfileでは「Select File to Open」、uiputfileでは「Select File to Write」。
・ファイルダイアログのOKボタンがuigetfileでは「開く」、uiputfileでは「保存」。

で、まぁそんな見た目だけの違いならやっぱりどっちでもいいんでは?という話になりますが、
ファイル読み込み用のuigetfileでは、存在しないファイル(この例ではhoge.m)を指定すると下記のようなダイアログがでます。
これでは新規ファイルを作って保存するときには使えません。




というわけで、ファイル書き込み用にはuiputfileを使うわけです。
逆にuiputfileで存在するファイル(この例ではhoge.m)を指定すると上書きするかどうかをとちゃんと聞いてくれます。




まぁここからへんは一般的なWindowsアプリケーションのファイルダイアログはこのような挙動をしますから、深く考える必要はないかと思いますけど。

参考:MATLAB Function Reference - uigetfile
参考:MATLAB Function Reference - uiputfile
参考:初心者によるMATLABメモ| ファイルダイアログから読み込むファイルを指定する
参考:初心者によるMATLABメモ| ディレクトリ(フォルダ)選択ダイアログを使用する

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

このエントリに付けられたタグ|MATLABファイル操作ファイルダイアログGUI書き込みuigetfile()uiputfile()

[Geometry] 円錐のZ座標をX,Y座標から求める

これはちょっと人に聞かれて作ったもので、実用性があるのかないのかよくわかりませんが、Z軸を軸とする円錐のZ座標をX,Y座標から求める関数cornz()です。
簡単な使い方,コード、使用例の順に示します。(ダウンロードはこちらから)

まずは、簡単な使用方法です。
Z = cornz(X, Y, R, H) … 原点を中心とし、最大半径Rで高さHの円錐のZ座標をX,Yから計算し、Zとして取得します。(円錐の外側にあるX,Yが指定された場合にはNaN)
Z = cornz(X, Y, R, H, M) … この場合,円錐の外側にあるX,Yが指定された場合には0が返ります。

コードはこんな感じです。エラー処理とヘルプ部分は省いてあります。
あくまでも最大半径の円がx,y平面上にあり、Z軸を軸とした円錐のZ座標を求めるための関数です。
平行移動、回転などをさせたい場合は、このあたりを参考にしてください。
function z = cornz(x, y, rm, hh, varargin) 
    rr = sqrt(x.^2 + y.^2);
    z = (hh/rm)*(rm - rr);

    if nargin==5
        z(find(rr>rm)) = 0;
    else
        z(find(rr>rm)) = NaN;
    end;

使用例として、実際に-15<x,y<15の範囲で,半径12、高さ15の円錐を書いてみました。
% プロットするxyの範囲
range = 15;

% 分割数
nn=301;

% 円錐の最大半径
rmax =12;

% 円錐の高さ
hmax =15;

[x, y] = meshgrid(linspace(-range,range,nn));
z = cornz(x, y, rmax, hmax);

figure('Position',[100,100,380,260]);
mesh(x, y, z, 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting','phong');
axis equal;

円錐(cornz,NaNモード)

上記の例の場合円錐の外側の場合にはZ=NaNとなり表示されませんが、下記のように第5引数を与えると円錐の外側の場合にはZ=0として表示します。
z = cornz(x, y, rmax, hmax, 1);
figure('Position',[100,100,380,260]);
mesh(x, y, z, 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting','phong');
axis equal;

円錐(cornz,0モード)

余談ですが、私に尋ねた方は下記のようなイメージを書きたかったようです。用途を聞いたのですがスペクトルがうんぬん言われて僕には意味不明だったので聞くのはやめましたw
figure('Position',[100,100,380,260]);
imagesc(z);
axis equal;
colorbar;

円錐(cornz、image表示)

参考:MATLAB Function Reference - atan2
参考:初心者によるMATLABメモ| 同次座標系の幾何学的な変換
ダウンロード:このエントリのmファイル(cornz()) 左クリックでお願いします。 (リファラを送信しない場合は,アクセス拒否になってしまいます)

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

このエントリに付けられたタグ|MATLABGeometry幾何円錐cornz

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

ディレクトリ(フォルダ)選択ダイアログを使用する
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 - ジャンル:コンピュータ

[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バーンスタイン関数

[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] 点列から主法線ベクトルを求める(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で表示した曲線の点列から数値的に曲率を求める関数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

[小技] MATLABからメールを送信する

意外と知らない人もいますが、MATLABからメールを送信することができます。
これが意外と便利なので紹介しておきます。
例えば、どれくらいかかるかわからない長い計算をかけている場合、「終わったかな?」とわざわざチェックするのは面倒です。
そこで、僕は計算が終わったらメールが来るようにしています。

以下、MATLABからメールを送信する方法です。
まず、基本的なSMTPサーバ(送信用のメールサーバ)の設定を行います。
ここらへんはstartup.mに書いておくと便利でしょう。ただし、パスワードも記入する場合には取り扱いに注意してください。
% SMTPサーバの設定
setpref( 'Internet', 'SMTP_Server', 'SMTPサーバ');
% SMTPサーバでのユーザ名(必要なら)
setpref( 'Internet', 'SMTP_Username', 'ユーザ名');
% SMTPサーバでのパスワード(必要なら)
setpref( 'Internet', 'SMTP_Password', 'パスワード');
% 自分のアドレス(Fromになる)
setpref('Internet', 'E_mail', '自分のアドレス');

あとは下記のようにsendmail関数を使用すればメールを送ることが出来ます。
sendmail('相手のメールアドレス','題名','メッセージ')

また、添付ファイルをつけたい場合は以下です。
sendmail('相手のメールアドレス','題名','メッセージ','添付ファイルのパス')

その他、複数の人に送る場合、複数の添付ファイルを付けたい場合などはこちらを参考にしてください。

参考:MATLAB:マニュアル:サイバネットシステム - MATLAB Function Reference - sendmail
参考:MATLAB:マニュアル:サイバネットシステム - MATLAB Programming - E メールの送信
参考:初心者によるMATLABメモ| 起動時にデフォルトカラーなどを設定する (起動時にスクリプトを実行)

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

このエントリに付けられたタグ|MATLAB小技メール送信sendmailsetprefsmtp

[数値解析] 非線形方程式を解く(準ニュートン法:csolve)

たまたまここでで見たんですが、↓だそうです。
2月最後の授業で言ったように、Matlabのoptimization toolboxを持っている場合、この中の関数fsolveを用いることによって非線形の連立方程式の解を求めることができます。しかし、その後調べたところ、数値計算に関する教科書を著しているKenn Juddは彼のホームページでこの関数はイマイチ(”poor”)であると述べ(リンク先参照)、それよりもChristopher Simsの書いた関数コードであるcsolve.mの利用を推奨しています。…

参考:コメント:Matlab関数fsolveについて

で、実際にリンクをたどってみると、確かに書いてありますね。そんなものを売るな!という話もありますが、まぁそれはそれ。無料の方が強力なら、それもよしです。
A Matlab Nonlinear Equation Solver The nonlinear equation solver in Matlab, fsolve, is a poor one based on minimizing the sum of squares of the functions. Chris Sims offers an improvement; see his ftp page.…

参考:Chapter 5 public code

とりあえずcsolve()が欲しいなら、ここからダウンロードできます。

以下で、csolve()の使い方を書いておきます。
非線形の方程式ソルバーcolve()は,準ニュートン法ベースで、通常の解法で極値,振動などで解が求まらない場合,ランダム方向に解を検索することで解を改善しようとします.
実際の使い方は下記です。
[x,rc] = csolve(FUN,x0,gradfun,crit,itmax)
[x,rc] = csolve(FUN,x0,gradfun,crit,itmax,varargin)

% [引数]
FUN        :    解かれる非線形システム方程式。
            FUNはベクトルxを引数として計算した非線形方程式のベクトルFを返します。
            ベクトルFはベクトルxと常に同じサイズである必要があります。
            関数FUNは、関数ハンドルを使って指定することができます。
x0        :    解の探索のための初期値.
gradfun    :    ヤコビ行列の関数.
crit    :    FUNが返す絶対値の合計がこの値より小さいならば、解が収束したと判断されます。
itmax    :    許可する繰り返しの最大回数.もし繰り返しの回数がこの値に達した場合、第2出力rcに4が返ります。
varargin:    関数FUNとgradfunに渡す、ベクトルx以外のパラメータを指定します。

この関数は実際に使った例を示します。
ここでは、下記の非線形方程式を解いてみます。

この式はfsolveのときと同じ方程式です。
方程式の関数は下記です。
% 方程式の関数
function F = funcs(x) 
    F = zeros(numel(x),1);
    F(1) = x(1)^2+x(2)^2-4;
    F(2) = x(1)*x(2)-x(3);
    F(3) = x(1)*x(3)-1;
    F = F(:);

% 関数funcsのヤコビ行列
function J = funcsj(x) 
    J = zeros(numel(x));
    J(1,1) = 2*x(1);
    J(2,1) = x(2);
    J(3,1) = x(3);

    J(1,2) = 2*x(2);
    J(2,2) = x(1);
    J(3,2) = 0;

    J(1,3) = 0;
    J(2,3) = -1;
    J(3,3) = x(1);

実行例は、下記です。fsolveのときと同じ解になっていることがわかります。
>> [x,exitflag] = mcsolve(@funcs, ones(3,1), @funcsj, 1e-12, 100, 4, 1)
itct 1, af 1.0208, lambda 0.36, rc 0
   x        0.64       1.72       1.36 
   f      -0.632    -0.2592    -0.1296 
itct 2, af 0.0452359, lambda 1, rc 0
   x    0.735407    1.86822    1.35976 
   f   0.0310718  0.0141413 -2.28221e-005 
itct 3, af 8.75519e-005, lambda 1, rc 0
   x    0.733074    1.86082     1.3641 
   f  6.0167e-005 1.72542e-005 -1.01306e-005 
itct 4, af 3.64242e-010, lambda 1, rc 0
   x    0.733077    1.86081    1.36411 
   f  2.9954e-010 -4.20057e-011 2.26965e-011 
itct 5, af 8.88178e-016, lambda 1, rc 0
   x    0.733077    1.86081    1.36411 
   f  8.88178e-016          0          0 
x =
    0.7331
    1.8608
    1.3641


exitflag =
     0

ただし、この関数csolve()は少々使いづらい部分もありますね。そのあたりは次のエントリに書きたいと思います。

参考:コメント:Matlab関数fsolveについて
参考:csolve()のダウンロード先
参考:【文献調査】BFGSの基礎
参考:初心者によるMATLABメモ| 非線形方程式を解く(fsolve)

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

このエントリに付けられたタグ|MATLABOptimizationToolbox方程式非線形方程式fsolvecsolve