動かして学ぶバイオメカニクス #2 〜多関節構造のプログラミング〜
身体特有の多関節構造はプログラミングにも影響する．そこで身体の力学問題の第１歩は，身体の力学計算とプログラミングとの関係について考える．最初のこの設計が肝心だ．ここでは身体の多関節構造と再帰呼び出しの関係について述べる．
身体の多関節構造とプログラミング
前章
で，バイオメカニクスにおける逆動力学計算（inverse daynamics）の意味について述べた．ここから，早速，身体構造の定義をPythonのプログラムを通して考える．
はじめに
前述のように，身体運動をもっとも力学的に特徴づけているのは身体の多関節構造（articulated structure），すなわち筋骨格系（musculoskeletal system）である．
図１に身体の多関節構造の例を示した．多関節構造の表現方法については後述するが，各部位（例，上腕）の近位側（身体の中心に近い側，根元側）の関節に座標系の原点を固定し，RGBの色で表した座標軸を示している．たとえば，図２の右上腕では右肩関節にRGBで表した座標系の原点が固定されている．
このような身体構造に対する力学計算を行う際に，文献１を参考に簡潔にプログラミングを行う方法について考えていきたい．
身体の分割・構成方法
身体構造はヒューマノイドロボットとも非常に似た構造を持つ．機械工学やロボットではロボットアームの一つ一つをリンク（link）と呼び，リンクとリンクは関節で接続される．このため，ロボットのような構造はリンク機構（link structure，またはリンク系，link system）と呼ばれる．身体の多関節構造もリンク機構と考えて良い．
ヒトの大腿，下腿，上腕，下腿，胸部などの部位が，各リンクに相当する．そこで，まず，図３のようなリンクの構造を定義する．
図３の右のように各リンクを分解し，その各リンクの根本側（近位側）には関節が固定され，各リンクと１つの根元側の関節をセットとする構造を取る．たとえば，右上腕（4 RUARM）の根本には肩関節がセットになって固定されている．しかし，頂点である腰部（1 HIP）には，このルールから除外され，この分割方法において関節は存在しない．ただし，関節を固定する受け口（ソケット）は３つあるが．
図４に示すように，このリンク機構は骨盤を頂点とする親子構造のツリー構造（tree structure）によって，わかりやすく表現できる．ただし，この親子関係のツリー構造は，部位によって子供の数が異なることが，プログラミング上あまり嬉しくない．この場合，３つの場合もあれば，１つしか子供を持たないことがある．
二分木による表現
そこで，次に，親子関係だけ記述するこのツリー構造と似ているが，文献１で述べられている記述方法を参考にし．データ構造とアルゴリズム（data structures and algorithms）で扱われる二分木（binary tree）で記述する（図５）．
ツリー構造では親子関係しか存在しないが，二分木では親子関係と兄弟姉妹関係を与える．ここでは左下に子供，右下に兄弟姉妹の枝分かれで記述している．ツリー構造には存在しない，左右に意味を与えている．
なお，一番上の親をツリー構造の頂点と呼ぶ．ここでは，HIPを頂点としている．
たとえば，図３，４で見比べると，腰部（1）の子供には，胸部（2)，右大腿（10），左大腿（13）の３つの部位が存在し，このときこれらの３つは兄弟姉妹の関係にある．その場合，右下に兄弟姉妹を書く．たとえば，胸部（2）の右下に右大腿（10），右大腿（10）の右下に左大腿（13）を書くことができる．これらの３つの順番は異なっていても良い．そして，終端を０と置くことで，それ以上の計算も終える．
このようなルールに基づいた二分木の構造を表現することで，常に下位の部位が２つになる．そして，全身計算を行う際に，骨盤（HIP）から，後述する再帰呼び出しを利用して，すべての部位を「同じプログラムで」たどることができる．
兄弟姉妹を計算から除外するという（図５の４RUARM以下の右下の７LUARM以下を計算しない）オプション（次章以降にPythonのコードを示すことになるが，sister=0というオプション）を与えば，右上腕から計算を始めることで，右腕（右上腕，右前腕，右手）という部分的な質量を計算することもできる．
図４のツリー構造を記述するのは容易だが，図５の二分木を構成するのは若干面倒であるかもしれない．しかし，身体の構造を決め順番を一度作ってしまえば，二分木の構造を変更する必要はない．仮に構造を変更しても，プログラムに変更は発生しない．
このような構造を利用することで，後ほど，リンク間の接続関係の記述の仕方次第で，プログラミングも効率的に行えることが理解できるだろう．なお，文献１にはMatlabのコードが記載されているので，Matlabを使用される方は参照されると良い．
リンク構造と漸化式
ここで親子関係だけに絞って考える．兄弟姉妹もその延長である．前述の身体全体の関係を配列b_linkに埋め込み表１に示した．この親子関係とたとえば，4 RUARM（右上腕）の子供は5 RFARM（右前腕），その子供は6 RHAND（右手）になる．4 RUARM → 5 RFARM → 6 RHANDとつながる．
これらの部位を4,5,6とIDを割り振り，各部位は質量（mass）を持つので，これをルール化すると
（b_link[4]以下の合計の質量） ＝ ４番目の部位の質量 + （b_link[5]以下の合計の質量）
と書ける．
ここで，４番目の部位（RUARM）の質量をmass[4]とし．事前に
mass[4] = 0.3
mass[5] = 0.2
mass[6] = 0.1などのように値を保持しておけばよい．
これを一般化し，（ｎ以下の合計の質量）をResultant_Mass(ｎ)という関数で置き換えると，
Resultant_Mass(4) = mass[4] + Resultant_Mass(5)
Resultant_Mass(5) = mass[5] + Resultant_Mass(6)
Resultant_Mass(6) = mass[6] + 0
と書ける．最後のResultant_Mass(7)，すなわち7以下の部位はないので0とした．これをさらに数字をnと置いて，一般化すると，
Resultant_Mass(ｎ) = mass[n] + Resultant_Mass(ｎ＋１)
ただし，もしn = 7のとき，Resultant_Mass(n)は０である．
と書ける．Python風プログラムにするなら
def Resultant_Mass(ｎ):
if n == 7:
return 0
else
return mass[n] + Resultant_Mass(ｎ＋１)となる．
ちなみに，Resultant_Mass(4)は
Resultant_Mass(4) = mass[4] + Resultant_Mass(5)
= mass[4] + (mass[5] + Resultant_Mass(6))
= mass[4] + mass[5] +(mass[6]+ Resultant_Mass(7))
= mass[4] + mass[5] + mass[6]+ 0
のように書き換えることができ，部位４以下の合計の質量が計算できていることが式から確認できる．プログラムを実行すれば，
mass[4] + mass[5] + mass[6] = 0.3 + 0.2 + 0.1 = 0.6
となる．n番目の部位の質量に相当するmass[n]を数式の$${m(n)}$$に置き換え，さらに Resultant_Mass(n) を数式の$${a(n)}$$に置き換えれば，親子関係は
$${a(n) = m(n) + a(n+1)}$$
という数式のルールに書き換えることができ，$${a(n)}$$を，$${n}$$と次の$${n+1}$$の式（または関数）だけで記述したこの形式を数学では，漸化式（recurrence relation）と呼ぶ．
このため，リンク機構も子供のリンク機構との関係だけで記述できると，特に加算の演算が簡潔に記述できる．forループを使用しないプログラムの簡潔さは，間違いを減らしてわかりやすくすることに繋がる．ちなみに，データ構造とアルゴリズムを学んだ時，プログラムでforループはできるだけ使用しないのが望ましいと叩き込まれた．
漸化式と再帰
プログラムで親子関係だけで記述し「プログラムの関数の中で，自分自身の関数を呼び出す」記述を，再帰呼び出し（recursive call）と呼ぶ．漸化式そのものである．
以上から，Pythonにおける漸化式の再帰（recursion）表現は，関数Resultant_Mass()を用い，
b_link = [0] * 7
b_link[4] = 0.3
b_link[5] = 0.2
b_link[6] = 0.1
def Resultant_Mass(n):
if n == 7:
return 0
else:
return b_link[n]+Resultant_Mass(n+1)
Resultant_Mass(4)
### 0.6のように定義することで記述できる．
前にも示したように，Resultant_Mass(4)がResultant_Mass(5)を呼び出し，Resultant_Mass(5)がResultant_Mass(6)を呼び出すという連鎖を繰り返して計算する．なお，永遠に計算を繰り返して（無限ループ）は困るので，ここではn=7のときに終了するように０を返して，Resultant_Mass()の計算を終える．なお，ここでは，コードの中で実行結果を###の後に示した．
クラスと再帰呼び出しによるプログラミング
ここでは，各部位の質量を加算する実際のPythonのプログラムを，クラス（class）と組み合わせて再帰呼び出しを記述する例を示す．
親子関係と姉妹関係はchild_valとsister_valに記述するが，とりあえず空欄にし，表２のように与える．
なお，Pythonのクラスの説明をここでも多少行うが，クラスの詳細については文献２などを参考に別途，学習されたい．
class BodyLink:
def __init__(self, l_id:int, name, mass:float):
self.l_id = l_id # linkのID
self.name = name # linkの名前
self.mass = mass # linkの質量分配比
self.child_val = None # 子供の変数
self.sister_val = None # 兄弟姉妹の変数ここで，class BodyLink: によってBodyLinkクラスの定義を行っている．Pythonにおけるクラスの定義がわからなくても，
val_name, l_id, name, sister_val, child_val, massからクラスが，データをしまうこむための表２の構造（入れ物の形）を決めていることが推察できる．これは，Cなどで扱われる構造体（structure）と同じである．このクラスに，先程示した身体の構造を示す二分木構造のchild_val（子）, sister_val（兄弟姉妹．ただし記号はsister）を記述する．
今後，このクラスに，各リンクの慣性モーメントや，関節の原点，姿勢，重心位置などを加えていくことになる．
以下では，表２のデータを格納するコードを以下に示す．b_linkという配列にクラス変数を格納しているが，いったん親子関係（chidl_val），兄弟姉妹（sister_val）欄は空欄（None）にする．
b_link = [0] * 16 #配列の初期化
b_link[0]=BodyLink(0, '', .0) # 0のときストップ
b_link[1]=BodyLink(1, 'Hip', .187)
b_link[2]=BodyLink(2, 'Chest', .302)
b_link[3]=BodyLink(3, 'Head', .069)
b_link[4]=BodyLink(4, 'RUArm', .027)
b_link[5]=BodyLink(5, 'RFArm', .016)
b_link[6]=BodyLink(6, 'RHand', .006)
b_link[7]=BodyLink(7, 'LUArm', .027)
b_link[8]=BodyLink(8, 'LFArm', .016)
b_link[9]=BodyLink(9, 'LHand', .006)
b_link[10]=BodyLink(10, 'RThigh', .110)
b_link[11]=BodyLink(11, 'RShin', .051)
b_link[12]=BodyLink(12, 'RFoot', .011)
b_link[13]=BodyLink(13, 'LThigh', .110)
b_link[14]=BodyLink(14, 'LShin', .051)
b_link[15]=BodyLink(15, 'LFoot', .011)そして，以下の，set_child_sister()を利用して，後から表３のグレー部分のように，親子関係（chidl_val），兄弟姉妹（sister_val）を追加する．
残念ながら，b_linkを与えると同時に，一度に親子関係と兄弟姉妹関係を，b_linkに与えることができないので，後から次のようにchild_valとsister_valに格納している．
# BodyLinkを一旦作成してから，その後に親子，兄弟姉妹関係をBodyLinkに与える．
def set_child_sister(self, child_val, sister_val):
self.child_val= child_val
self.sister_val = sister_val
# 合成の質量（全質量）の計算
# 親子関係と兄弟姉妹関係を下記でb_linkのchild_val, sister_valに追記
b_link[1].set_child_sister(b_link[2], b_link[0])
b_link[2].set_child_sister(b_link[3], b_link[10])
b_link[3].set_child_sister(b_link[0], b_link[4])
b_link[4].set_child_sister(b_link[5], b_link[7])
b_link[5].set_child_sister(b_link[6], b_link[0])
b_link[6].set_child_sister(b_link[0], b_link[0])
b_link[7].set_child_sister(b_link[8], b_link[0])
b_link[8].set_child_sister(b_link[9], b_link[0])
b_link[9].set_child_sister(b_link[0], b_link[0])
b_link[10].set_child_sister(b_link[11], b_link[13])
b_link[11].set_child_sister(b_link[12], b_link[0])
b_link[12].set_child_sister(b_link[0], b_link[0])
b_link[13].set_child_sister(b_link[14], b_link[0])
b_link[14].set_child_sister(b_link[15], b_link[0])
b_link[15].set_child_sister(b_link[0], b_link[0])これで計算の準備が整ったので，以下のメソッドの再帰関数Resultant_Mass()で，全身の質量を計算する．なお，クラス変数massには質量の分配比が格納されている．体重にこの比をかけることで，各部位（例えば右上腕：RURM）の質量を計算する．
# 合成の質量（全質量）の計算
# sister == 0 なら，部分解析（そのリンクの遠位側だけ計算）
# sister != 0 （例えば１）なら，全身解析（姉妹を含める）
# bw: 体重
def Resultant_Mass(self, bw, sister=0):
if self.l_id == 0:
return 0.0
elif sister == 0: # 部分（そこから遠位のみの）解析
return bw*self.mass + (self.child_val).Resultant_Mass(bw, sister)
else: # 全身解析
return (bw*self.mass +(self.child_val).Resultant_Mass(bw, sister) +
(self.sister_val).Resultant_Mass(bw, sister))当たり前だが，各部位のmassの和は１なので，体重（bw = 86. [kg]）を与えると，体重86.0 [kg]の人の全質量和は86.0 [kg]となる．
ここでは，各b_link[n]に，クラスBodyLinkで定義された表３のデータを与えているが，表３からもb_link[1]はリンク機構の頂点の腰（HIP）に相当する．そこで，HIP以下の全リンクに対して，総和の質量を計算する関数Resultant_Mass()で計算をするなら，
b_link[1].Resultant_Mass(86, sister=1)のように
変数 ＋ ドット（ . ）＋ 関数
と書くことで，クラスBodyLinkの変数（ここではb_link[1]）に対して関数Resultant_Mass()を計算することができる．
なお，sisterオプションは，ツリーの子供だけを計算するか（sister = ０の場合），兄弟姉妹も計算するか（sister =1）の違いである．この計算には体重も引数として与える（この例では86.）．
そこで，b_link[1]に対して，それ以下の総和の質量を計算する関数Resultant_Mass()を上記のように適用すると，
###[出力結果] 86.0当たり前だが，86 [kg]を得る（ここでは実行結果を###の後に示した）．一方，右腕だけの質量を，右上腕（RUARM）である b_link[4]を与えて，それ以下の子供だけをsister=0のオプションで指定して，計算すると
b_link[4].Resultant_Mass(86., sister=0)
###[出力結果] 4.214のように，右腕の合計質量 bw*（RUARM+RFARM+RHAND）= 86.*(.027+.016+.006) = 4.214 が計算される．
このように再帰呼び出しを利用することで，全身も部分的な計算も，きれいなコードで記述することができる．
おわりに
ここでは，身体構造を二分木で表現し，それを再帰呼び出しでプログラミングする方法について述べた．この方法に固執する必要はないのだが，forループで計算するよりも利点は大いにある．
たとえば，二分木と再帰を活用することで，上記の質量だけでなく，たとえば
b_link[4].Force(86., sister =0)
b_link[4].Torque(86., sister =0)と書くことで，リンク４（右上腕）に作用する力や関節に作用するトルクを今後定義するForce()やTorque()関数で計算することができ，関数の記述も簡潔になり，誤りが減る．
以後のシリーズでは，このクラスに様々な変数（慣性モーメント，角速度，関節位置）と，様々な逆動力学計算の関数（メソッド）をを加筆して，逆動力学計算を実現することになる．クラスと再帰呼び出しを活用することで，力学計算のコードを簡潔に記述することができるだろう．
コードのまとめ：
クラスBodyLinkの定義（メソッド，set_child_sister，Resultant_Massを含む）：
class BodyLink:
def __init__(self, l_id:int, name, mass:float):
self.l_id = l_id # linkのID
self.name = name # linkの名前
self.mass = mass # linkの質量分配比
self.child_val = None # 子供の変数
self.sister_val = None # 兄弟姉妹の変数
# BodyLinkを一旦作成してから，その後に親子，兄弟姉妹関係をBodyLinkに与える．
def set_child_sister(self, child_val, sister_val):
self.child_val= child_val
self.sister_val = sister_val
# 合成の質量（全質量）の計算
# sister == 0 なら，部分解析（そのリンクの遠位側だけ計算）
# sister != 0 （例えば１）なら，全身解析（姉妹を含める）
# bw: 体重
def Resultant_Mass(self, bw, sister=0):
if self.l_id == 0:
return 0.0
elif sister == 0: # 部分（そこから遠位のみの）解析
return bw*self.mass + (self.child_val).Resultant_Mass(bw, sister)
else: # 全身解析
return (bw*self.mass +(self.child_val).Resultant_Mass(bw, sister) +
(self.sister_val).Resultant_Mass(bw, sister))身体の構造のデータ：
b_link = [0] * 16 #配列の初期化
b_link[0]=BodyLink(0, '', .0) # 0のときストップ
b_link[1]=BodyLink(1, 'Hip', .187)
b_link[2]=BodyLink(2, 'Chest', .302)
b_link[3]=BodyLink(3, 'Head', .069)
b_link[4]=BodyLink(4, 'RUArm', .027)
b_link[5]=BodyLink(5, 'RFArm', .016)
b_link[6]=BodyLink(6, 'RHand', .006)
b_link[7]=BodyLink(7, 'LUArm', .027)
b_link[8]=BodyLink(8, 'LFArm', .016)
b_link[9]=BodyLink(9, 'LHand', .006)
b_link[10]=BodyLink(10, 'RThigh', .110)
b_link[11]=BodyLink(11, 'RShin', .051)
b_link[12]=BodyLink(12, 'RFoot', .011)
b_link[13]=BodyLink(13, 'LThigh', .110)
b_link[14]=BodyLink(14, 'LShin', .051)
b_link[15]=BodyLink(15, 'LFoot', .011)
# 親子関係と兄弟姉妹関係を下記でb_linkのchild_val, sister_valに追記
b_link[1].set_child_sister(b_link[2], b_link[0])
b_link[2].set_child_sister(b_link[3], b_link[10])
b_link[3].set_child_sister(b_link[0], b_link[4])
b_link[4].set_child_sister(b_link[5], b_link[7])
b_link[5].set_child_sister(b_link[6], b_link[0])
b_link[6].set_child_sister(b_link[0], b_link[0])
b_link[7].set_child_sister(b_link[8], b_link[0])
b_link[8].set_child_sister(b_link[9], b_link[0])
b_link[9].set_child_sister(b_link[0], b_link[0])
b_link[10].set_child_sister(b_link[11], b_link[13])
b_link[11].set_child_sister(b_link[12], b_link[0])
b_link[12].set_child_sister(b_link[0], b_link[0])
b_link[13].set_child_sister(b_link[14], b_link[0])
b_link[14].set_child_sister(b_link[15], b_link[0])
b_link[15].set_child_sister(b_link[0], b_link[0])を定義することで，
# b_link[1]（ツリー構造の最上位）に対して計算することで，それ以下のツリーを計算．
# sister = １で兄弟姉妹を含める．
# すなわち全リンク（全身）の質量を以下で計算する．
b_link[1].Resultant_Mass(86, sister=1)
###[出力結果] 86.0
# 4:右上腕以下（右腕全体）の質量を計算する．
# sister = 0で兄弟姉妹を含めない（左腕などを含めない）．
b_link[4].Resultant_Mass(86, sister=0)
###[出力結果] 4.214のような実行結果を得る．
次章
では，具体的に日本人の身体の質量分布モデルを利用し，慣性パラメータを計算する．
参考文献
１）ヒューマノイドロボット（改訂２版），梶田秀司編著，オーム社，2020
２）Python ゼロからはじめるプログラミング，三谷純著，翔泳社，2021
【著作権・転載・免責について】
権利の帰属
本ホームページで提示しているソフトウェアならびにプログラムリストは，スポーツセンシング社の著作物であり，スポーツセンシング社に知的所有権がありますが，自由にご利用いただいて構いません．
本ページに掲載されている記事，ソフトウェア，プログラムなどに関する著作権および工業所有権については，株式会社スポーツセンシングに帰属するものです．非営利目的で行う研究用途に限り，無償での使用を許可します．
転載
本ページの内容の転載については非営利目的に限り，本ページの引用であることを明記したうえで，自由に行えるものとします．
免責
本ページで掲載されている内容は，特定の条件下についての内容である場合があります． ソフトウェアやプログラム等，本ページの内容を参照して研究などを行う場合には，その点を十分に踏まえた上で，自己責任でご利用ください．また，本ページの掲載内容によって生じた一切の損害については，株式会社スポーツセンシングおよび著者はその責を負わないものとします．
【プログラムの内容について】
プログラムや内容に対する質問に対しては，回答できないことのほうが多くなると思いますが，コメントには目は通します．回答は必要最低限にとどめますので，返信はあまり期待しないでいただけると幸いです，
「動かして学ぶ」という大それたタイトルをつけたものの，また，きれいなプログラムに対するこだわりはあるものの，実際のプログラミングのスキルは決して高くありません．最下部の方のコメント欄によるプログラムの間違いのご指摘は歓迎します．できるだけ反映します．
【解析・受託開発について】
スポーツセンシングでは，豊富な知見を持つ，研究者や各種エンジニアが研究・開発のお手伝いをしております．研究・開発でお困りの方は，ぜひスポーツセンシングにご相談ください．
【例】
・データ解析の代行
・受託開発
（ハードウェア、組込みソフトウェア、PC/モバイルアプリ）
・測定システム構築に関するコンサルティング など
その他，幅広い分野をカバーしておりますので，まずはお気軽にお問い合わせください．
株式会社スポーツセンシング
【ホームページ】sports-sensing.com
【Facebook】sports.sensing
【Twitter】Sports_Sensing
【メール】support@sports-sensing.com