FC2ブログ

初心者によるMATLABメモ

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

 

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

[小技] OSの環境変数を使用する

MATLABを使っていて、OSの環境変数を使用できると便利な場面が意外とあると思うのですが、そんな場合はgetenvを使用することで簡単に実現できます。
N = getenv('環境変数名')

例えば、自分はシミュレーションを複数の計算機でかけてるんですが、進捗状況メールがどの計算機から送られてきたのか混乱することも度々…
ということを繰り返した結果、今では進捗状況をメールするときに下記のようにサブジェクトにコンピュータ名を含めています。
subject = strcat(['[', getenv('COMPUTERNAME'), ']', 'テストメール']);


一応書いておくと、Windowsの場合は環境変数が見たければ、コマンドプロンプトで「set」コマンドを引数なしで実行すれば見れます。

参考:MATLAB 関数リファレンス - getenv
スポンサーサイト



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

このエントリに付けられたタグ|MATLAB環境変数getenv小技

[基本] 作業ログの保存などに使う関数diaryの拡張

簡単な作業ログの保存などに使う関数diaryを拡張したdiary2()です。(mファイルのダウンロードはこちら)
基本的にはdiary()と機能、使い方は変わりません。

違う点は、diaryだと追加書き込みに自動的になってしまいますが、diary2だと新規書き込みを行うことができます。
また、開始時と停止時にタイムスタンプを自動的に書き込みます。
あとは現在のdiaryファイルの名前の表示、diaryのon/offの表示の程度の機能が追加されています。

使い方は基本的にはdiary()と同じです。
diary2 … diary('on')と同じ
diary new … 以後の出力をdiaryファイルに書き込みます。(ファイルが存在しても新規書き込み)
diary on … 以後の出力をdiaryファイルに書き込みます。(ファイルが存在すれば追加書き込み)
diary off … 書き込みを中止します。

diary stat … diaryモードのON/OFFを表示します。
diary fname … 現在指定されているdiaryファイル名を表示します。

diary('ファイル名') … diaryファイル名を指定して書き込みを開始します。
diary('ファイル名','on') … diaryファイル名を指定し、追加書き込みを行うことができます。
diary('ファイル名','new') … diaryファイル名を指定し、新規書き込みを行うことができます。

※diary2('file')は、diary2('file','on')と同じです。

on/off/newあたりは普通のdiaryと使い方が同じなので、ここでは省略するとして、一応statとfnameだけ例を載せておきます。
一連の流れとして使用するとこんな感じです。
>> diary2('diary.txt','new')
>> diary2 fname
現在のdiaryのファイル名 ... [ diary.txt ]
>> diary2 stat
現在のdiaryの状態 ... [ on ]
>> diary2 off
>> diary2 stat
現在のdiaryの状態 ... [ off ]

最後にソースを載せておきます。
また、diary2()のヘルプ付きのmファイルのダウンロードはこちらから行えます。
function diary2(varargin) 

%% 引数のチェックと処理
error(nargchk(0,2,nargin));

opt = 0;
diary_path = get(0, 'DiaryFile');
mode = {'on','new','off','stat','fname'};

for ii = 1:nargin 
    in = varargin{ii};

    if ischar(in) 
        hit = find(strcmpi(lower(in), mode));
        if isempty(hit) 
            if ~strcmpi(diary_path,in) 
                diary2 off;
                diary_path = in;
            end;
        else 
            opt = hit(1);
        end;
    else 
        error('引数は文字列である必要があります');
    end;
end;

if opt==0 
    if strcmpi(get(0,'Diary'),'on') 
        opt=3;
    else 
        opt=1;
    end;
end;
%% 


%% 実際の処理
switch opt 
    % 書き込み開始 
    case {1,2} 

        % diaryファイルへの書き込み開始メッセージ書き込み 
        if strcmpi(get(0,'Diary'),'off') 
            if opt==1 
                perm = 'a+';
            else 
                perm = 'w+';
            end;
            fid = fopen(diary_path,perm);
            fprintf(fid, '%s¥n', '%%% =================  diary 書き込み開始 ==================== ');
            fprintf(fid, '%s%s¥n¥n', '+ 現在時間 : ', datestr(now,31));
            fclose(fid);
        end;

        % diary on
        diary(diary_path);
    % 書き込み停止 
    case 3
        % diaryファイルへの書き込み終了メッセージ書き込み 
        if strcmpi(get(0,'Diary'),'on') 
            diary off;

            fid = fopen(diary_path,'r+');

            body = {};
            frewind(fid);
            while 1 
                tline = fgetl(fid);
                if ~ischar(tline), break, end;
                body{length(body)+1} = tline;
            end;

            body{end} = [];

            frewind(fid);
            for ii=1:length(body)
                fprintf(fid, '%s¥n', body{ii});
            end;
            fprintf(fid, '%s%s¥n', '+ 現在時間 : ', datestr(now,31));
            fprintf(fid, '%s¥n', '%% ========== / diary 書き込み停止 ==========');

            fclose(fid);
        end;
    % ステータス 
    case 4 
        disp(sprintf('現在のdiaryの状態 ... [ %s ]', get(0,'Diary')));
    % ファイル名表示 
    case 5 
        disp(sprintf('現在のdiaryのファイル名 ... [ %s ]', get(0,'DiaryFile')));
    otherwise
        help diary2;
end


参考:MATLAB Function Reference - diary
ダウンロード:このエントリのmファイル 左クリックでお願いします。 (リファラを送信しない場合は,アクセス拒否になってしまいます)

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

このエントリに付けられたタグ|MATLABdiary()作業ログ

[数値解析] 非線形方程式を解く(準ニュートン法: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

[小技] 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

[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 - ジャンル:コンピュータ

[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

次のページ