スポーツバイオメカニクス MatlabとOpenSIM

スポーツバイオメカニクスの分析によく用いられるMatlabとOpenSIMの使用に関する備忘録

Matlabが起動しない?OS X Yosemiteに更新する前に確認を。

アップルユーザーの皆さん,こんにちは.
OS X Yosemiteがリリースされましたね.
私もさっそくアップデートしました.
そして,今回のテーマである,Matlabが起動しないという問題にぶち当たりました(笑)
起動しようとすると,こんなエラーメッセージが表示されます.
Matlab error
ちなみに,アイコンはこんな感じになっています.
Matlab icon
結論から言うと,現段階では2013以前のバージョンを起動するようにする方法はわかりません.
素直に2014にバージョンアップするのが賢明です.
ただし,2014バージョンにも2014aと2014bがあります.
公式にサポートが表明されているのは2014bのみです.
したがって,2014bを購入しなければなりません.
しかし,ここでひとつの問題が浮上します.
それは・・・
2014bにstudentが存在しないということです.
通常の個人用を買うとおよそ15000円です.
15000円にして,toolboxはついていません.
一方のstudentは,およそ10000円にしてtoolboxがいくつか付属しています.
さらに,toolboxを必要としない人であれば,5000円足らずで購入できます.
toolboxも同時購入なら,1つ1000円ほどで購入できますので,必要とするtoolboxが少ないのなら,この方法で買った方がお得です.
2014bを買える人は概ね問題なく使用でき,2014aなら学生も購入できます.
解決したかのように見えますが,実はもうひとつ問題が潜んでいます.
なんと...2014aはインストールしただけでは使えないのです.
まず,パッチを適用しないといけません.
patch
ダウンロードしたら,指示に従ってインストールすればOKです.
次に,2014aを起動する際には一工夫が必要で,ターミナルから起動しないといけないのです.
パッケージの中にあるStartMATLABでも起動できます.

f:id:simulate:20170131195830p:plain

f:id:simulate:20170131195851p:plain

f:id:simulate:20170131195924p:plain

最小二乗法

今回は最小二乗法による回帰式の求め方を紹介します。

最小二乗法は一般的には、疑似逆行列を用いる方法ではなく、偏微分を用いる方が多いですが、Matlabの処理では疑似逆行列を用いた方が簡単に行えます。

Matlabは行列計算ソフトなので、当然と言えば当然です。

 

y=ax+b

という形の回帰式を求めたい場合、aとbが求まれば回帰式が作れるわけです。

なので、最小二乗法でaとbを求めます。

そのためには下のような行列式を作ります。

f:id:simulate:20170131200300p:plain

この式のxとyにそれぞれの値を代入していきます。

行列式は次のように置き換えて考えることができます。

f:id:simulate:20170131200322p:plain

 

行列式で表記すると複雑な気がしますが、AX=Bとおくと、とても簡単な式に見えます。

これをx=の形にしてやると、答えが出てきます。

なので、

f:id:simulate:20170131200338p:plain

と言うように解いてしまいたいのですが、Aは正則行列ではないので逆行列を作ることができません。

そこで、疑似逆行列を作ります。

Matlabには疑似逆行列を計算する関数がありますので

pinv(A);

とすれば良いだけです。

なので、疑似逆行列をA+として

f:id:simulate:20170131200232p:plain

 

という計算ができます。

これはMatlabでの表記は、さきほどの疑似逆行列の部分と合わせて書くと

ans=pinv(A)*B;

となります。

すると、出てきた答えは2行1列の行列になっているはずです。

念のために説明しておくと、1行目がa、2行目がbの値になっています。

同様の方法で最小二乗法でいろいろなものを求めることができます。

 

疑似逆行列を用いた最小二乗法の説明はあまり見かけないので、簡単に説明してみました。

私の理解も不十分でわかりにくい点もあったかもしれませんが、参考にしてみてください。

角度の算出 —atan2—

前回はacosを使った余弦定理を用いた角度算出方法を紹介しました。

今回はatan2を使った正接の角度算出方法を紹介します。

前回は2つのベクトルのなす角の求め方でしたが、今回は1つのベクトルのある平面内のひとつの軸に対する角度を求めます。例えば、XY平面の中でX軸となす角度などです。

f:id:simulate:20170131200514p:plain

関数atan2の使い方は上の図の通りです。

では、これをどのように使うとベクトルの角度の算出方法になるのか説明していきます。

f:id:simulate:20170131200556p:plain

仮に、上の図のように斜辺がベクトルaだとすると、そのx成分とy成分を使ってその角度を求めることができます。

ベクトルaの成分がx,yの順番で引数に入っているとすると、上の図のような書き方で角度を求めることができます。
なお、これで算出される角度はradianですので、degreeに変換する場合は
「*180/pi」
と最後に付け足してやる必要があります。

この図では「θ」と記号を使っていますが、Matlabの引数に記号は使えないので「theta」などとすればそのまま使用可能です。

これを使うと、セグメント角度の算出が可能になります。

例えば、大腿セグメントの角度を算出する場合、まずは大腿セグメントベクトルを作らなければなりません。
近位から遠位に向かうベクトルとして、大腿セグメントベクトルを作る場合は、膝関節の座標値から股関節の座標値を引きます。

ベクトルの引き算は高校の数学Bの範囲ですので、わからない方は調べてみてください。
これはそんなに難しい概念ではありません。


次に、さきほど引き算で作成したベクトルの角度を算出します。
さっき説明したように、大腿セグメントベクトルのx成分とy成分をatan2関数に入れてやるだけです。

ちなみに、atan2では-180度から180度の間で角度が求められます。
したがって、atan2(y,x)とそのまま入れた場合は

f:id:simulate:20170131200622p:plain

というような角度の求め方になります。

x成分とy成分の順番を入れ替えたり、-を付けたりすることでどこを基準に角度を求めるかが変わります。
求めたい角度が概ね何度くらいかわかるものから試してみて、どこで求めるのが望ましいか考えなら行うことが良いと思います。

2つのベクトルのなす角

バイオメカニクスでは2つのベクトルのなす角を求めることが求められます。

例えば、関節角度です。

今回はその方法の1つである余弦定理を用いた方法を紹介したいと思います。

余弦定理とは、高校生でも知っているあの余弦定理です。

余弦定理の詳細についてはどこかのサイトを探してもらう、もしくは数学の教科書に載っています。

なので、実際の運用方法について書いていきたいと思います。

余弦定理で2つのベクトルがなす角を求めるために必要なのは、2つのベクトルの内積とそれぞれの大きさ(ノルム)です。

では、作業の方法に入ります。

まずは、座標データから2本のベクトルを作ります。

仮に、膝関節角度を求める場合は膝関節から股関節に向かうベクトルと膝関節から足関節に向かうベクトルの2本を作成します。

なので・・・

vector1=股関節-膝関節
vector2=足関節-膝関節

となります。

rowという引数の1行目から3行目に股関節のX、Y、Z座標、4行目から6行目に膝関節のX、Y、Z座標が入っているとした場合は

vector1=row(1:3,:)-row(4:6,:);

となります。

よりわかりやすくするために、

hip=row(1:3,:);
knee=row(4:6,:);

というように、わかりやすい引数にそれぞれの座標値を入れておくと次のようにより簡単に作れます。

vector1=hip-knee;

次に内積を求めます。

内積は1コマごとに求める必要があります。

したがって、1コマ目の内積

InnerProduct=dot(vector1(:,1),vector2(:,1));

で求められます。

dotは内積の関数です。

これをforを使って全コマ分求めることができます。

コマ数をFrameとしたとき、

for time=1:Frame
InnerProduct(time)=dot(vector1(:,time),vector2(:,time));
end

「for」の使い方は他のサイトを参照してもらうか、後々書きたいと思います(笑

内積が求められたら、次はそれぞれのベクトルのノルムです。

for time=1:Frame
Norm1=norm(vector1(:,time));
Norm2=norm(vector2(:,time));
end

これでノルムが全コマ分求められます。

ノルムはコマによって変わらないような気がしますが、ディジタイズ誤差によって1コマごとに変動していますので、全コマ分別に求める必要があります。

最後に余弦定理で角度を求めます。

余弦定理の式を変形すると、cosの値が得られます。

つまり、ここまでに用いてきた引数で書くと

cos=InnerProduct/(Norm1*Norm2)

になります。

これをforで全コマ分のcosの値を求めると

for time=1:Frame
cos(time)=InnerProduct(time)/(Norm1(time)*Norm2(time));
end

となります。

そして、求められたcosの値を「acos」に入れると、角度になります。

ただし、ここで算出された角度はradianなので、degree(度)に直すためには180/πをかけてやる必要があります。

したがって、膝関節の角度は・・・

angle=acos(cos)*180/pi;

で求められることになります。

いかがでしたでしょうか?

少し難しかったかもしれませんが挑戦してみてください。

なお、この方法では180度を超える角度は求められません。

180度を超える角度を求める方法はまた今度紹介したいと思います。

WindowsとMac

Matlabを使うにあたって、OSによる悩みがあるのではないかと思います。

 

ほとんどの方はWindowsを使っているでしょうからあまり困ることはないでしょう。

しかし、研究者にはMacが多いということも耳にします。

また、Macが使いたいという方も多いのではないでしょうか。

 

MatlabWindowsでもMacでも動きます。

唯一の違いはMacでは英語しか使えないということです。

Windowsでは日本語に対応しているのでエラーも日本語で返されます。ただし、よくわからない日本語になっていることも・・・。

一方、Macでは日本語に対応していないので、メニューもエラーもすべて英語です。

そして、使用するフォルダ、ディレクトリもすべて半角英数字でなければなりません。

したがって、現在のフォルダに指定するフォルダはすべて英語のディレクトリでなければならないのです。

ユーザー名が日本語の場合、書類フォルダの中に英語の名前のフォルダを作ってもそのフォルダを使うことはできません。

 

もし、MacMatlabが使いたくて買ったものの、うまく動かないという方はこれが問題になっている可能性が考えられます。

 

 

MacMatlabを使う場合、もうひとつ問題が生じることがあります。

それは、ディレクトリ名に「¥」が使えないことです。

Windowsでは「¥」によってディレクトリが記されますが、Macでは「¥」は「\」として表示されます。

どちらにせよ、円マークでもバックスラッシュでもMacではディレクトリを正しく読み取ることができません。

そのような場合は、「\(バックスラッシュ)」を「/(スラッシュ)」に置き換えてみてください。

そもそもディレクトリを書き間違えていない限り、これで正しく動くはずです。

なお、Windowsでも「/」は正常に動作するので、バックスラッシュよりもスラッシュを使うことをお勧めします。

 

データの取り扱い

スポーツバイオメカニクスの分析のためにMatlabを使う場合、まずはデータの扱い方から勉強しなければなりません。

バイオメカニクスでは主に座標値のデータを扱います。

X、Y、Zそれぞれの座標データがポイント数分あります。

基本的にディジタイズで得られた座標値は23点か25点を用いると思いますので、69もしくは75個の数字×コマ数分あることになります。

これらのデータは行にポイント、列に時間として扱うのが良いと考えています。

f:id:simulate:20170131200937p:plain

なので、この例では1行目は1ポイント目のX座標、2行目がY座標、3行目がZ座標となっており、4行目は2ポイント目のX座標、5行目がY座標、6行目がZ座標・・・というようになっています。

具体例を挙げると、上から頭頂のX・Y・Z座標、耳珠点のX・Y・Z座標という具合です。

ちなみに、Frame-DIASからエクスポートしたデータをそのまま読み込むと行は時間、列はポイントになっています。

なので、読み込んだものを転置するとこの形になります。

したがって、読み込んでこの形にするには次のような方法を取ります。

Frame-DIASでエクスポートしたデータが「data.3d」というファイルで現在のフォルダーに入っているします。

これをrowという引数に入れるときは・・・

FDdata=dlmread('data.3d');

[m,n]=size(FDdata);

row=FDdata(2:m,:)';

の3行でできます。

1行目でデータをFDdataに読み込みます。

2行目でデータのサイズを取得します。

そして、3行目で不要な部分を切り取り、転置することで冒頭で説明した形になります。

さらに・・・

[Point,Frame]=size(row)

と次の行に入れておけば、ポイント数とコマ数をそれぞれPointとFrameという引数に取得することができます。

しかし、実はこんなことをしなくてもFDdataの1行目にこれらのデータに加えて刻み時間も入っているので、そちらを取得した方が賢いかもしれません。