動かして学ぶバイオメカニクス#8 〜身体の慣性モーメント〜
回転の運動方程式を解く上で，個人ごとの慣性モーメントが必要となる．ここでは，第３章で紹介した同じモデルを使用して，身体各部位の慣性モーメントを計算する方法について述べる．
おさらい
第３章で身体重心と質量分布のモデルを示した．このモデルは日本人のアスリートの体形を計測し，体形を均一な密度の楕円の厚みのある板で近似し，身体各部の重心と質量を計算し，体重と身体各部位の長さで記述された統計モデルが構築されていた．
第３章では，各部位の重心位置と質量を算出したが，ここでも同じモデルで計算された慣性モーメントの算出モデルについて述べていく（文献１参照）．
日本人アスリートを対象とした慣性モーメントの統計モデル
文献１より，身体各部位の慣性モーメントの計算に，体重と身体各部位の長さの情報が必要となる．
各軸の各部位の慣性モーメントは
$$
I_i = a_i w + b_i l + c_i~~~~~i=x,y,z
$$
の線形な統計モデルから計算される．ここで，$${I_i}$$は軸$${i}$$の慣性モーメント，$${W}$$は体重，$${l}$$は部位の長さ（図２，３に示す赤色の点の近位端と遠位端間の長さ），$${a_i, b_i, c_i}$$はパラメータで，表１，２に示した．文献１に示された値は小数点以下の桁が一定ではないので，ここでは小数点第二位まで示した．
また，軸$${i=x,y,z}$$の各軸は，$${x}軸を前額軸 (左右軸)方向，$${y}$$軸を矢状軸（前後軸）方向，$${z}$$軸を長軸とした．ただし足部の軸方向は，文献１から改変している．
各部位を構成する，近位端と遠位端の位置を示した第３章の図を図２，３に再度掲載する．図中の赤色の点が，近位端と遠位端を構成する解剖学的特徴点となり，これらの点間の距離を$${l}$$として使用する．モーションキャプチャであれば，ポーズを計測し，その際の距離の平均を使用するか，運動中の平均を使用することになるだろう．なお，マーカの貼付の仕方やモーションキャプチャのシステムにも依存するが，運動中の$${l}}$の変化は小さくないかもしれない．
表１：慣性モーメント推定係数例：男子 （文献１）
表２：慣性モーメント推定係数例：女子 （文献１）
動かして学ぶ慣性モーメントの計算
第６章
のコードに，３つの慣性主軸の慣性モーメントを計算するコードを加筆する．
なお，コードは，これまで同様に，以下のGoogle Colaboratoryのリンクから，ブラウザでPythonコードを実行できるようにした．なお，ログインするためには，Googleアカウントが必要となるので注意をされたい．
また，後半で使用するモーションキャプチャのサンプルデータ（sample_optitrack_221027.csv）とフォースプレートのサンプルデータ（sample_fp_221027.csv）は，このリンクからダウンロードして，Google Driveをマウントし，マウントしたドライブにそれらをコピーしてご使用いただきたい．これは，第６章で使用したデータと同じであるので，すでにダウンロ度済みの方は，サイドダウンロードする必要はない．
パッケージの読み込み
Pandas（データの読み込み），Numpy（行列演算），matplotlib（グラフ）のパッケージの読み込みを行う．
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
#from scipy.interpolate import interp1d
import copy
from scipy.interpolate import UnivariateSplineデータをダウンロードして，ユーザのファイルの置き場所にあわせて以下を変更する．
from google.colab import drive
drive.mount('/content/drive')以下は各自の環境にあわせて設定していただきたい．前述のように，すでにダウンロード済みの方は，同じファイルを使用されたい．
file_path_op_8 = '/content/drive/MyDrive/.../sample_optitrack_221027.csv'motion capture の csv file からマーカなどの情報を抽出する関数群
以下は，これまでと同様に，motion captureのデータからマーカの位置ベクトルを抽出する関数．
# CSVデータからヘッダだけ抽出
def extract_df_header(file_path):
return pd.read_csv(file_path, header=None, nrows=5, index_col=1).T
# skeletonの名前を抽出
def extract_skeleton_name(file_path):
df = extract_df_header(file_path)
marker_name = df['Name'][2]
position_col = marker_name.find(':')
return marker_name[:position_col]
# ファイル名とマーカの名前を引数に与えて各Linkの重心の加速度データを抽出する関数
def extract_marker(file_path, marker_name):
sk_name = extract_skeleton_name(file_path) # skeleton
df_main = pd.read_csv(file_path, header=[4]) # データだけ抽出
df = extract_df_header(file_path) # headerだけ抽出
df_selected = df[df['Name']== sk_name+':'+marker_name] # ヘッダがNameと一致する部分だけ抽出
selected_rows = df_selected.index # 一致する列
marker_data = np.array(df_main.iloc[:, selected_rows[-3:]].values) # 後ろから３個データを取得す
return marker_data[:, [2,0,1]] # OptiTrackの座標がYup（Y軸が上）のため，座標を入れ替えることで，Zupに変更するさらに，以下は以前と同様に，モーションキャプチャのデータから位置ベクトルを抽出し，関節位置（遠位端を含む）を計算する関数である．
# (new6)
# リンクの原点の位置ベクトル
def joint_vec(b_link, l_id, path):
if l_id == 0:
return .0
else:
joint_name = b_link[l_id].name
return extract_marker(path, joint_name)class BodyLinkの定義（逆動力学計算の本体）
class BodyLinkの定義
新しいインスタンスとして，inertia_moment（慣性モーメントの配列）とp（各部位の原点＝関節の位置ベクトル）を加えた．
# new8：追記・変更箇所
class BodyLink:
def __init__(self, l_id:int, name, mass_ratio:float, cg_ratio:float):
self.l_id = l_id # linkのID
self.name = name # linkの名前
self.mass_ratio = mass_ratio # linkの質量分配比
self.cg_ratio = cg_ratio # linkの重心位置の分配比 (new6)
self.inertia_moment = None # linkの慣性モーメント (new8)
self.child_val = None # 子供の変数
self.sister_val = None # 兄弟姉妹の変数
self.distal_val = None # linkの遠位端の変数 (new6)
self.p = None # linkの関節の原点（近位端） (new8)
self.cg = None # linkの重心位置
self.acc = None # linkの重心の加速度
self.force = None # linkの関節に作用する力
# BodyLinkを一旦作成してから，その後に親子，兄弟姉妹関係をBodyLinkに与える．
def set_child_sister(self, child_val, sister_val):
self.child_val= child_val
self.sister_val = sister_val
# (new6)
# リンク（部位）の重心を計算するためには，リンクの原点（近位端）と遠位端の位置ベクトルが必要
# ただし，末端のリンク（部位）には子供のリンク（部位）がないので，リンクの遠位端（関節）の情報を得られない
# BodyLinkを一旦作成してから，その後に，各リンクの遠位端の位置をBodyLinkに与える．親子関係とは区別して利用
def set_distal_end(self, distal_val): # (new6)
self.distal_val= distal_val
# m * (acc - g)
def CgForce(self, bw):
gravity = [0., 0., -9.8] # 重力加速度ベクトル
return bw * self.mass_ratio * (self.acc - gravity)
# Inverse Dynamics Newton equation（改訂）
def ID_Newton(self, bw, sister=False, fp=True):
if self.name == 'RToe': # (new6)右足の場合，床反力2を返す
return self.force
elif self.name == 'LToe': # (new6)左足の場合，床反力1を返す
return self.force
elif self.l_id == 0:
return 0.0
elif sister == False: # 部分（そこから遠位：子供のみ）の解析
f = self.CgForce(bw)+\
self.child_val.ID_Newton(bw, sister)
self.force = f
return f
elif sister == True: # 全身解析（兄弟姉妹を含む）
f = self.CgForce(bw)+\
self.child_val.ID_Newton(bw, sister)+\
self.sister_val.ID_Newton(bw, sister)
self.force = f
return f
else:
print('Put True or False as sister (2nd) option')
print('If sister == True, this calculate also sister tree')
print('If sister == False (default), this calculate only child tree')
# 合成の質量（全質量）の計算
# sister == False なら，部分解析（そのリンクの遠位側だけ計算）
# sister == True なら，全身解析（兄弟を含める）
# bw: 体重
def Resultant_Mass(self, bw, sister=False):
if self.l_id == 0:
return 0.0
elif sister == False: # 部分（そこから遠位のみの）解析
return bw * self.mass_ratio + (self.child_val).Resultant_Mass(bw, sister)
elif sister == True: # 全身解析
return (bw * self.mass_ratio +(self.child_val).Resultant_Mass(bw, sister) +
(self.sister_val).Resultant_Mass(bw, sister))
else:
print('Put True or False as sister (2nd) option')
print('If sister == True, this calculate also sister tree')
print('If sister == False (default), this calculate only child tree')
# m * cg
def WeightedJoint(self, bw):
return bw * self.mass_ratio * self.cg
# 合成の（全）WeightedJoint
# Σ m_i * cg_i
# sister == False なら，部分解析（そのリンクの遠位側だけ計算）
# sister == True なら，全身解析（兄弟を含める）
def Resultant_WeightedJoint(self, bw, sister=False):
if self.l_id == 0:
return 0.0
elif sister == False: # 部分（そこから遠位のみの）解析
return self.WeightedJoint(bw)+\
(self.child_val).Resultant_WeightedJoint(bw, sister) #, sister)
elif sister == True: # 全身解析
return self.WeightedJoint(bw)+\
self.child_val.Resultant_WeightedJoint(bw, sister)+\
self.sister_val.Resultant_WeightedJoint(bw, sister)
else:
print('Put True or False as sister (2nd) option')
print('If sister == True, this calculate also sister tree')
print('If sister == False (default), this calculate only child tree')
# 合成重心
# Resultant_WeightedJoint() / Resultant_Mass()
# = Σ(m_i * cg_i) / Σ m_i
# sister == False なら，部分解析（そのリンクの遠位側だけ計算）
# sister == True なら，全身解析（兄弟を含める）
def Resultant_Cg(self, bw, sister=False):
return self.Resultant_WeightedJoint(bw, sister)/self.Resultant_Mass(bw, sister)
b_link = [0] * 21 # 配列の初期化
b_link[0] = BodyLink(0, '', .0, .0) # 0のときストップ
b_link[1] = BodyLink(1, 'Hip', .187, (1.-.609))
b_link[2] = BodyLink(2, 'Chest', .302, (1.-.428))
b_link[3] = BodyLink(3, 'Neck', .069, (1.-.821))
b_link[4] = BodyLink(4, 'RUArm', .027, .529)
b_link[5] = BodyLink(5, 'RFArm', .016, .415)
b_link[6] = BodyLink(6, 'RHand', .006, .891)
b_link[7] = BodyLink(7, 'LUArm', .027, .529)
b_link[8] = BodyLink(8, 'LFArm', .016, .415)
b_link[9] = BodyLink(9, 'LHand', .006, .891)
b_link[10] = BodyLink(10, 'RThigh', .110, .475)
b_link[11] = BodyLink(11, 'RShin', .051, .406)
b_link[12] = BodyLink(12, 'RFoot', .011, (1.-.595))
b_link[13] = BodyLink(13, 'LThigh', .110, .475)
b_link[14] = BodyLink(14, 'LShin', .051, .406)
b_link[15] = BodyLink(15, 'LFoot', .011, (1.-.595))
# 各リンクの遠位端も定義（遠位端の位置ベクトルを計算するためだけに必要な変数 (new6)
b_link[16] = BodyLink(16, 'Head', .0, .0) # 3 Neckの遠位端
b_link[17] = BodyLink(17, 'RFIN', .0, .0) # 6 RHandの遠位端
b_link[18] = BodyLink(18, 'LFIN', .0, .0) # 9 LHandの遠位端
b_link[19] = BodyLink(19, 'RToe', .0, .0) # 12 RFootの遠位端
b_link[20] = BodyLink(20, 'LToe', .0, .0) # 15 LFootの遠位端前回同様，親子関係，兄弟姉妹関係，遠位端の情報をクラスb_linkに格納する．
# 親子関係と兄弟姉妹関係を下記でb_linkのchild_val, sister_valに追記
child_sister_list = [
[1,2,0], [2,3,10], [3,0,4], [4,5,7], [5,6,0], [6,0,0], [7,8,0],
[8,9,0], [9,0,0], [10,11,13], [11,12,0], [12,0,0], [13,14,0],
[14,15,0], [15,0,0]]
[b_link[i[0]].set_child_sister(b_link[i[1]], b_link[i[2]]) for i in child_sister_list];
# 各リンクの遠位端の情報を格納．b_linkのcdistal_valに追記 (new6)
distal_end_list = [
[1,2], [2,3], [3,16], [4,5], [5,6], [6,17], [7,8],
[8,9], [9,18], [10,11], [11,12], [12,19], [13,14],
[14,15], [15,20]]
[b_link[i[0]].set_distal_end(b_link[i[1]]) for i in distal_end_list];各部位の長さを計算
ここから，慣性モーメントを計算する準備を行う．以下はベクトルのノルムを計算する関数と，インスタンス変数pに，関節の位置ベクトルを代入する関数を準備した．
# (new8)
# ノルム（大きさ）
def norm(vec_array):
return np.linalg.norm(vec_array, axis=1)
# (new8)
# 関節位置をBodyLinkクラスに格納
# 部位を指定し，関節位置ベクトルを計算し，BodyLinkのpに格納する関数
def set_link_p_list(b_link, l_id_list, path):
for n in l_id_list:
b_link[n].p = joint_vec(b_link, n, path)これを利用して，以下に関節の位置ベクトルを代入する．
# 関節位置をBodyLinkのpに格納
set_link_p_list(b_link, range(21), file_path_op_8)以下は，クラスにて定義されている，各リンク（部位）の近位端と遠位端の位置ベクトルから，各部位の長さを計算する関数．
# (new8)
# linkの長さ
def link_length(b_link, link_id):
return norm(b_link[b_link[link_id].distal_val.l_id].p - b_link[link_id].p)以下で，実際に運動中の全身の各部位の長さをプロットした．
# 運動中の全linkの長さ変化をプロット
[plt.plot(link_length(b_link, i)) for i in range(1,16)];
以下で，このグラフのデータの平均を計算し，b_lengthに格納した．
# 各部位の長さの平均
b_length = [np.average(link_length(b_link, i)) for i in range(1,16)]以下に，表１，２で示した慣性モーメントを計算するパラメータを辞書として与えて，
# 慣性モーメントを計算するためのパラメータ： para[0]*bw + para[1]*link_length + para[2]
# これらで計算した結果の慣性モーメントの単位は kg cm^2なので，10^(-4)をかけてkg m^2に変換する必要がある
# 0：左右軸，1：前後軸，2：長軸
# 男子用
i_m_para0 = [[.0, .0, .0], [.0, .0, .0], [.0, .0, .0]]
i_m_para1 = [[22.62, 5588.38, -1687.06], [27.7 , 6516.01, -1982.55],[29.08, 2246.6 , -1376.85]] # HIP（下胴）
i_m_para2 = [[58.01, 15247.8 , -6157.42], [71.52, 15063., -6423.4 ], [48.9 , 1516.61, -2016.55]] # Chest（上胴）
i_m_para3 = [[2.71, 2843.24, -367.9], [2.49, 2681.7 , -354.1 ], [1.25, 1307.37, -138.96]] # NECK（頭部）
i_m_para4_7 = [[1.85, 1007.85, -317.68], [1.74, 999.69, -312.14], [0.71, -44.88, -11.1]] # UpperArm（上腕）
i_m_para5_8 = [[0.86, 562.22, -145.87], [0.8, 576.66, -146.45], [0.24, 26.38, -13.48]] # Fore Arm（前腕）
i_m_para6_9 = [[ 0.11, 80.36, -6.37], [0.14, 82.07, -7.31], [0.05, 9.08, -1.67]] # Hand（手）
i_m_para10_13 = [[11.83, 5684.2 , -2127.91], [10.65, 5547.75, -2043.38], [6.63, 418.34, -350.31]] # Thigh（大腿）
i_m_para11_14 = [[5.27, 3093.33, -1190.24], [5.19, 3048.1, -1174.66], [1.11, 104.75, -62.79]] # Sin（下腿）
# 0：左右（関節）軸，1：上下軸，2：長軸
i_m_para12_15 = [[ 0.01, 214.58, -38.93], [0.01, 228.14, -40.98], [0.01, 37.67, -6.3 ]] # Foot（足）
inertia_para_dict_m = {1:i_m_para1, 2:i_m_para2, 3:i_m_para3,
4:i_m_para4_7, 5:i_m_para5_8, 6:i_m_para6_9, 7:i_m_para4_7, 8:i_m_para5_8, 9:i_m_para6_9,
10:i_m_para10_13, 11:i_m_para11_14, 12:i_m_para12_15, 13:i_m_para10_13, 14:i_m_para11_14, 15:i_m_para12_15}
# 女子用
i_w_para0 = [[.0, .0, .0], [.0, .0, .0], [.0, .0, .0]]
i_w_para1 = [[19.08, 5368.76, -1407.99], [25.86, 7013.95, -1884.32], [27.44, 3909.55, -1472.05]] # HIP（下胴）
i_w_para2 = [[36.67, 10495.1 , -3537.95], [44.82, 10796.7 , -3751.31], [33.43, 1489.2 , -1214.28]] # Chest（上胴）
i_w_para3 = [[2.03, 2109.44, -234.04], [1.81, 1973.49, -224.25], [0.97, 1085.02, -98.47]] # NECK（頭部）
i_w_para4_7 = [[1.28, 888.26, -243.45], [1.13, 890.14, -239.38], [0.46, 16.56, -15.05]]# UpperArm（上腕）
i_w_para5_8 = [[0.53, 387.15, -85.15], [0.49, 401.2, -86.13], [0.14, 22.58, -6.2 ]] # Fore Arm（前腕）
i_w_para6_9 = [[0.1 , 68.32, -6.56], [0.14, 69.02, -6.67], [0.04, 8.99, -1.07]] # Hand（手）
i_w_para10_13 = [[19.53, 4292.37, -1895.67], [17.66, 4347.35, -1851.78], [9.01, -14.12, -262.08]] # Thigh（大腿）
i_w_para11_14 = [[4.82, 2373.89, -840.81], [4.76, 2342.82, -830.82], [1.17, 67.14, -46.44]] # Sin（下腿）
# 0：左右（関節）軸，1：上下軸，2：長軸
i_w_para12_15 = [[0.14, 148.98, -30.67], [0.16, 153.66, -32.22], [0.04, 31.2, -6.5]] # Foot（足）
inertia_para_dict_w = {1:i_w_para1, 2:i_w_para2, 3:i_w_para3,
4:i_w_para4_7, 5:i_w_para5_8, 6:i_w_para6_9, 7:i_w_para4_7, 8:i_w_para5_8, 9:i_w_para6_9,
10:i_w_para10_13, 11:i_w_para11_14, 12:i_w_para12_15, 13:i_w_para10_13, 14:i_w_para11_14, 15:i_w_para12_15}
# (new8)
# 部位に対応するパラメータ，その部位の長さ，体重を与えて，慣性モーメントを計算
def inertia_moment(i_para, body_weight, ave_link_length):
return np.array([i_para[i][0]*body_weight + i_para[i][1]*ave_link_length + i_para[i][2] for i in range(3) ])*(10.**(-4))以上で，慣性モーメントを計算する準備が整った．
例題：
体重86 kgの被験者の慣性モーメントを実際に計算し，クラスb_linkの各インスタンス変数inertia_momentに代入した．
# 体重86.0kgの被験者の慣性モーメントをb_linkのクラスに格納
bw = 86.
for l_id in range(1,16):
b_link[l_id].inertia_moment = inertia_moment(inertia_para_dict_m[l_id], bw, b_length[l_id-1])以下で，その全身の慣性モーメントを表示した．
# 全慣性モーメント
print(np.array([b_link[n].inertia_moment for n in range(1,16)]))array([[ 2.09384436e-01, 2.53992787e-01, 1.86195831e-01],
[ 2.94663261e-01, 3.79263738e-01, 2.59816428e-01],
[ 2.67309421e-02, 2.39341115e-02, 1.53455128e-02],
[ 1.33606613e-02, 1.27320941e-02, 3.69488027e-03],
[ 6.60482241e-03, 6.38515288e-03, 1.36331563e-03],
[ 1.12814837e-03, 1.30957923e-03, 3.55556835e-04],
[ 1.33606612e-02, 1.27320940e-02, 3.69488027e-03],
[ 8.18850666e-03, 8.00951231e-03, 1.43762389e-03],
[ 1.05176442e-03, 1.23156988e-03, 3.46926094e-04],
[ 1.29989395e-01, 1.22508140e-01, 3.97269943e-02],
[ 5.43066745e-02, 5.33049594e-02, 7.60178118e-03],
[-5.18645896e-04, -5.15844229e-04, 3.32779341e-05],
[ 1.29989388e-01, 1.22508133e-01, 3.97269938e-02],
[ 5.43066894e-02, 5.33049741e-02, 7.60178168e-03],
[-5.18646060e-04, -5.15844404e-04, 3.32779053e-05]])次章以降について
以上でオイラーの運動方程式を解く上で，慣性モーメントの計算を終えた．次回は各部位の角速度を計算することを目標として，まずは各部位の座標系の定義方法について述べて行く予定である．
参考文献
１）阿江, 湯, 横井，日本人アスリートの身体部分慣性特性の推定，バイオメカニズム，Vol. 11, pp. 23--33, 1992
【著作権・転載・免責について】
権利の帰属
本ホームページで提示しているソフトウェアならびにプログラムリストは，スポーツセンシング社の著作物であり，スポーツセンシング社に知的所有権がありますが，自由にご利用いただいて構いません．
本ページに掲載されている記事，ソフトウェア，プログラムなどに関する著作権および工業所有権については，株式会社スポーツセンシングに帰属するものです．非営利目的で行う研究用途に限り，無償での使用を許可します．
転載
本ページの内容の転載については非営利目的に限り，本ページの引用であることを明記したうえで，自由に行えるものとします．
免責
本ページで掲載されている内容は，特定の条件下についての内容である場合があります． ソフトウェアやプログラム等，本ページの内容を参照して研究などを行う場合には，その点を十分に踏まえた上で，自己責任でご利用ください．また，本ページの掲載内容によって生じた一切の損害については，株式会社スポーツセンシングおよび著者はその責を負わないものとします．
【プログラムの内容について】
プログラムや内容に対する質問に対しては，回答できないことのほうが多くなると思いますが，コメントには目は通します．回答は必要最低限にとどめますので，返信はあまり期待しないでいただけると幸いです，
「動かして学ぶ」という大それたタイトルをつけたものの，また，きれいなプログラムに対するこだわりはあるものの，実際のプログラミングのスキルは決して高くありません．最下部の方のコメント欄によるプログラムの間違いのご指摘は歓迎します．できるだけ反映します．
【解析・受託開発について】
スポーツセンシングでは，豊富な知見を持つ，研究者や各種エンジニアが研究・開発のお手伝いをしております．研究・開発でお困りの方は，ぜひスポーツセンシングにご相談ください．
【例】
・データ解析の代行
・受託開発
（ハードウェア、組込みソフトウェア、PC/モバイルアプリ）
・測定システム構築に関するコンサルティング など
その他，幅広い分野をカバーしておりますので，まずはお気軽にお問い合わせください．
株式会社スポーツセンシング
【ホームページ】sports-sensing.com
【Facebook】sports.sensing
【Twitter】Sports_Sensing
【メール】support@sports-sensing.com