歩行の2次元動作解析(データ処理)【Python/Pandas】

# はじめに
前回Kinoveaを使用して算出した歩行中の各ランドマークをグラフにプロットし、関節角度を求めていきたいと思います。


例えば左側矢状面ではこんな感じです↓
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/5834dfb9-9143-ba7d-3fc5-21348deb2f4f.png)

※元動画の私です
左側下肢の各ランドマークの軌跡がグラフにプロットできていると思います。
![walkSide.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/0bc41a01-34ab-1e20-2f11-12afc6059d7e.gif)
左側踵接地から2回目の左踵接地までを切り出しています。
## 前回の内容

https://qiita.com/naitoyuma/items/1a54e0fbaf749807e167

- [x] Kinoveaを使用した映像内マーカーのトラッキング
- [x] カメラ撮影による誤差範囲の把握 
- [x] 角関節角度の算出に必要なランドマークの検討
- [x] 算出パラメータの検討 
- [x] PythonライブラリOpenPoseによる姿勢定位の紹介 

## 今回の内容
- [x] 前額面、左矢状面、右矢状面の各ランドマークのプロット直線の描画による視覚化
- [x] 体幹、股関節、膝関節、足関節の各運動方向の角度算出
- [x] 両肩峰、ASIS、膝蓋骨中央の位置から各傾斜角度を算出

## 環境構築
歩行時のランドマーク座標のcsvデータを含め、動作解析の過程をGithubで公開しています。

https://github.com/naitoyuma7110/walkAnalysis

Dockerで環境を作りJupyterLabを使って順番に算出していこうかと思います。
```Dockerfile:Dockerfile
FROM python:3.9.7-buster
ARG DEBIAN_FRONTEND=noninteractive

RUN apt-get update && apt-get install -y \
    tzdata \
    libgl1-mesa-glx \
&&  ln -sf /usr/share/zoneinfo/Asia/Tokyo /etc/localtime \
&&  apt-get clean \
&&  rm -rf /var/lib/apt/lists/*

ENV TZ=Asia/Tokyo

RUN python3 -m pip install --upgrade pip \
&&  pip install --no-cache-dir \
    black \
    jupyterlab \
    jupyterlab_code_formatter \
    jupyterlab-git \
    lckr-jupyterlab-variableinspector \
    jupyterlab_widgets \
    ipywidgets \
    import-ipynb

RUN pip install --no-cache-dir \
    numpy \
    pandas \
    scipy \
    matplotlib \
    japanize_matplotlib \
    seaborn \
    requests \
    Pillow
```

## データ分析

### 歩行中のイベント時間
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/1f098269-80e6-6844-c2e3-2be82f86be1c.png)

KinoveaのKey Image機能を使って切り出した各歩行周期のイベント時間を使用
下記の時間をもとに1歩行周期を切り出します。
```python
    key_image = pd.read_csv("/workspace/walkAnalysis/pointData/keyImage.csv")
    # NaNを空白文字列に置き換える
    key_image = key_image.fillna("")
    key_image
```

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>Front</th>
      <th>Time</th>
      <th>Side</th>
      <th>Time.1</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>右踵接地</td>
      <td>5.1</td>
      <td>左踵接地</td>
      <td>3.34</td>
    </tr>
    <tr>
      <th>1</th>
      <td>右全足底接地</td>
      <td>5.2</td>
      <td>左全足底接地</td>
      <td>3.47</td>
    </tr>
    <tr>
      <th>2</th>
      <td>左足趾離地</td>
      <td>5.34</td>
      <td>対側足趾離地</td>
      <td>3.60</td>
    </tr>
    <tr>
      <th>3</th>
      <td>左踵接地</td>
      <td>5.74</td>
      <td>対側踵接地</td>
      <td>4.04</td>
    </tr>
    <tr>
      <th>4</th>
      <td>左全足底接地</td>
      <td>5.87</td>
      <td>対側全足底接地</td>
      <td>4.17</td>
    </tr>
    <tr>
      <th>5</th>
      <td>右足趾離地</td>
      <td>5.9</td>
      <td>左足趾離地</td>
      <td>4.24</td>
    </tr>
    <tr>
      <th>6</th>
      <td>右踵接地(2回目)</td>
      <td>6.34</td>
      <td>左踵接地(2回目)</td>
      <td>4.67</td>
    </tr>
    <tr>
      <th>7</th>
      <td>右全足底接地</td>
      <td>6.47</td>
      <td>右踵接地</td>
      <td>10.34</td>
    </tr>
    <tr>
      <th>8</th>
      <td>左足趾離地(2回目)</td>
      <td>6.54</td>
      <td>右全足底接地</td>
      <td>10.44</td>
    </tr>
    <tr>
      <th>9</th>
      <td>左踵接地(2回目)</td>
      <td>6.97</td>
      <td>対側足趾離地</td>
      <td>10.54</td>
    </tr>
    <tr>
      <th>10</th>
      <td></td>
      <td></td>
      <td>対側踵離地</td>
      <td>10.94</td>
    </tr>
    <tr>
      <th>11</th>
      <td></td>
      <td></td>
      <td>対側全足底接地</td>
      <td>11.08</td>
    </tr>
    <tr>
      <th>12</th>
      <td></td>
      <td></td>
      <td>右足趾離地</td>
      <td>11.11</td>
    </tr>
    <tr>
      <th>13</th>
      <td></td>
      <td></td>
      <td>右踵接地</td>
      <td>11.54</td>
    </tr>
  </tbody>
</table>

### ランドマーク座標の視覚化
 歩行中の下記ランドマークx,y座標
 
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/366e25fe-54fd-0150-9ba9-a6a5da0ef260.png)

### 右矢状面のランドマーク座標
```python
right_side = pd.read_csv("/workspace/walkAnalysis/pointData/sideRightTrim.csv")
right_side = right_side[(right_side['t'] >= 10.34) & (right_side['t'] <= 11.54)]
right_side.head()
```

<table class="scroll-table" border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>t</th>
      <th>R_shoulder_x</th>
      <th>R_shoulder_y</th>
      <th>L_shoulder_x</th>
      <th>L_shoulder_y</th>
      <th>R_asis_x</th>
      <th>R_asis_y</th>
      <th>L_asis_x</th>
      <th>L_asis_y</th>
      <th>R_knee_x</th>
      <th>R_knee_y</th>
      <th>L_knee_x</th>
      <th>L_knee_y</th>
      <th>R_foot_x</th>
      <th>R_foot_y</th>
      <th>L_foot_x</th>
      <th>L_foot_y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.14</td>
      <td>-123.98</td>
      <td>167.38</td>
      <td>18.45</td>
      <td>162.79</td>
      <td>-113.57</td>
      <td>-10.49</td>
      <td>-4.45</td>
      <td>-4.84</td>
      <td>-97.83</td>
      <td>-182.18</td>
      <td>-27.57</td>
      <td>-179.59</td>
      <td>-95.89</td>
      <td>-325.78</td>
      <td>-40.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>1</th>
      <td>5.17</td>
      <td>-124.98</td>
      <td>168.38</td>
      <td>17.45</td>
      <td>161.80</td>
      <td>-114.58</td>
      <td>-10.49</td>
      <td>-5.45</td>
      <td>-6.84</td>
      <td>-97.82</td>
      <td>-184.17</td>
      <td>-27.55</td>
      <td>-181.59</td>
      <td>-94.29</td>
      <td>-331.05</td>
      <td>-40.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>2</th>
      <td>5.20</td>
      <td>-124.98</td>
      <td>170.38</td>
      <td>15.45</td>
      <td>162.80</td>
      <td>-116.58</td>
      <td>-10.49</td>
      <td>-7.44</td>
      <td>-7.84</td>
      <td>-98.81</td>
      <td>-184.19</td>
      <td>-26.57</td>
      <td>-188.61</td>
      <td>-94.66</td>
      <td>-334.34</td>
      <td>-39.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>3</th>
      <td>5.24</td>
      <td>-132.98</td>
      <td>167.38</td>
      <td>14.13</td>
      <td>164.80</td>
      <td>-118.56</td>
      <td>-7.49</td>
      <td>-8.44</td>
      <td>-7.84</td>
      <td>-100.82</td>
      <td>-184.20</td>
      <td>-26.58</td>
      <td>-189.61</td>
      <td>-93.44</td>
      <td>-335.33</td>
      <td>-40.76</td>
      <td>-301.22</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.27</td>
      <td>-133.97</td>
      <td>172.38</td>
      <td>12.51</td>
      <td>168.42</td>
      <td>-120.55</td>
      <td>-3.49</td>
      <td>-9.44</td>
      <td>-6.84</td>
      <td>-101.84</td>
      <td>-182.20</td>
      <td>-26.59</td>
      <td>-189.61</td>
      <td>-93.45</td>
      <td>-335.33</td>
      <td>-37.94</td>
      <td>-295.30</td>
    </tr>
  </tbody>
  <tfoot>
    <tr><th colspan=18>$\hspace{88em}$</th></tr>
    </tfoot>
　</table>

### 左矢状面のランドマーク座標
```python
left_side = pd.read_csv("/workspace/walkAnalysis/pointData/sideLeftTrim.csv")
left_side = left_side[(left_side['t'] >= 3.34) & (left_side['t'] <= 4.67)]
left_side.head()
```

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>t</th>
      <th>L_ancle_x</th>
      <th>L_ancle_y</th>
      <th>L_heel_x</th>
      <th>L_heel_y</th>
      <th>L_toe_x</th>
      <th>L_toe_y</th>
      <th>L_leg_x</th>
      <th>L_leg_y</th>
      <th>L_knee_x</th>
      <th>L_knee_y</th>
      <th>L_hip_x</th>
      <th>L_hip_y</th>
      <th>L_asis_x</th>
      <th>L_asis_y</th>
      <th>L_shoulder_x</th>
      <th>L_shoulder_y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>3.37</td>
      <td>345.83</td>
      <td>-389.40</td>
      <td>348.67</td>
      <td>-422.05</td>
      <td>290.12</td>
      <td>-404.02</td>
      <td>374.13</td>
      <td>-219.73</td>
      <td>373.50</td>
      <td>-208.10</td>
      <td>461.05</td>
      <td>-28.42</td>
      <td>423.92</td>
      <td>24.78</td>
      <td>482.62</td>
      <td>263.05</td>
    </tr>
    <tr>
      <th>1</th>
      <td>3.40</td>
      <td>340.85</td>
      <td>-392.40</td>
      <td>346.68</td>
      <td>-422.04</td>
      <td>286.27</td>
      <td>-414.14</td>
      <td>365.13</td>
      <td>-223.73</td>
      <td>364.51</td>
      <td>-208.11</td>
      <td>448.05</td>
      <td>-31.42</td>
      <td>411.91</td>
      <td>22.79</td>
      <td>473.61</td>
      <td>261.06</td>
    </tr>
    <tr>
      <th>2</th>
      <td>3.44</td>
      <td>335.85</td>
      <td>-393.42</td>
      <td>344.67</td>
      <td>-422.07</td>
      <td>283.43</td>
      <td>-420.21</td>
      <td>355.13</td>
      <td>-224.73</td>
      <td>356.51</td>
      <td>-208.11</td>
      <td>433.05</td>
      <td>-32.42</td>
      <td>399.90</td>
      <td>21.81</td>
      <td>463.61</td>
      <td>259.05</td>
    </tr>
    <tr>
      <th>3</th>
      <td>3.47</td>
      <td>333.84</td>
      <td>-394.43</td>
      <td>346.74</td>
      <td>-423.04</td>
      <td>281.51</td>
      <td>-423.25</td>
      <td>344.52</td>
      <td>-228.60</td>
      <td>346.51</td>
      <td>-208.12</td>
      <td>415.05</td>
      <td>-31.42</td>
      <td>387.91</td>
      <td>22.83</td>
      <td>452.61</td>
      <td>256.05</td>
    </tr>
    <tr>
      <th>4</th>
      <td>3.50</td>
      <td>332.82</td>
      <td>-395.43</td>
      <td>346.74</td>
      <td>-422.06</td>
      <td>279.51</td>
      <td>-424.23</td>
      <td>338.52</td>
      <td>-228.60</td>
      <td>338.59</td>
      <td>-209.75</td>
      <td>409.29</td>
      <td>-30.42</td>
      <td>375.92</td>
      <td>24.85</td>
      <td>441.62</td>
      <td>254.05</td>
    </tr>
  </tbody>
    <tfoot>
    <tr><th colspan=18>$\hspace{88em}$</th></tr>
    </tfoot>
</table>

### 前額面のランドマーク座標
```python
front = pd.read_csv("/workspace/walkAnalysis/pointData/frontTrim.csv")
# 歩行周期の時間に合わせデータフレームを切り出す
front = front[(front['t'] >= 5.10) & (front['t'] <= 6.97)]
front.head()
```

<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>t</th>
      <th>R_shoulder_x</th>
      <th>R_shoulder_y</th>
      <th>L_shoulder_x</th>
      <th>L_shoulder_y</th>
      <th>R_asis_x</th>
      <th>R_asis_y</th>
      <th>L_asis_x</th>
      <th>L_asis_y</th>
      <th>R_knee_x</th>
      <th>R_knee_y</th>
      <th>L_knee_x</th>
      <th>L_knee_y</th>
      <th>R_foot_x</th>
      <th>R_foot_y</th>
      <th>L_foot_x</th>
      <th>L_foot_y</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>5.14</td>
      <td>-123.98</td>
      <td>167.38</td>
      <td>18.45</td>
      <td>162.79</td>
      <td>-113.57</td>
      <td>-10.49</td>
      <td>-4.45</td>
      <td>-4.84</td>
      <td>-97.83</td>
      <td>-182.18</td>
      <td>-27.57</td>
      <td>-179.59</td>
      <td>-95.89</td>
      <td>-325.78</td>
      <td>-40.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>1</th>
      <td>5.17</td>
      <td>-124.98</td>
      <td>168.38</td>
      <td>17.45</td>
      <td>161.80</td>
      <td>-114.58</td>
      <td>-10.49</td>
      <td>-5.45</td>
      <td>-6.84</td>
      <td>-97.82</td>
      <td>-184.17</td>
      <td>-27.55</td>
      <td>-181.59</td>
      <td>-94.29</td>
      <td>-331.05</td>
      <td>-40.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>2</th>
      <td>5.20</td>
      <td>-124.98</td>
      <td>170.38</td>
      <td>15.45</td>
      <td>162.80</td>
      <td>-116.58</td>
      <td>-10.49</td>
      <td>-7.44</td>
      <td>-7.84</td>
      <td>-98.81</td>
      <td>-184.19</td>
      <td>-26.57</td>
      <td>-188.61</td>
      <td>-94.66</td>
      <td>-334.34</td>
      <td>-39.76</td>
      <td>-307.55</td>
    </tr>
    <tr>
      <th>3</th>
      <td>5.24</td>
      <td>-132.98</td>
      <td>167.38</td>
      <td>14.13</td>
      <td>164.80</td>
      <td>-118.56</td>
      <td>-7.49</td>
      <td>-8.44</td>
      <td>-7.84</td>
      <td>-100.82</td>
      <td>-184.20</td>
      <td>-26.58</td>
      <td>-189.61</td>
      <td>-93.44</td>
      <td>-335.33</td>
      <td>-40.76</td>
      <td>-301.22</td>
    </tr>
    <tr>
      <th>4</th>
      <td>5.27</td>
      <td>-133.97</td>
      <td>172.38</td>
      <td>12.51</td>
      <td>168.42</td>
      <td>-120.55</td>
      <td>-3.49</td>
      <td>-9.44</td>
      <td>-6.84</td>
      <td>-101.84</td>
      <td>-182.20</td>
      <td>-26.59</td>
      <td>-189.61</td>
      <td>-93.45</td>
      <td>-335.33</td>
      <td>-37.94</td>
      <td>-295.30</td>
    </tr>
  </tbody>
  <tfoot>
<tr><th colspan=18>$\hspace{88em}$</th></tr>
</tfoot>
</table>

### 各座標間を結ぶ直線を描画する

各身体部位に対応するカラム名を配列にしてループさせmatplotlibを使い直線を描画しました。
下記は前額面ですが、これを左右矢状面でも実施しました。
```python
# グラフの描画
plt.figure(figsize=(10, 6))
plt.title('前額面')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
plt.axis('equal')
plt.grid(True)

landamrks = ['shoulder','asis','knee','foot']

previous_landmark = None

for landmark in landamrks:
  plt.scatter(front[f'R_{landmark}_x'], front[f'R_{landmark}_y'],  marker='.', s=5)
  plt.scatter(front[f'L_{landmark}_x'], front[f'L_{landmark}_y'],  marker='.', s=5)
  
  # 各点をつなぐ線を描画(斜度を算出する直線)
  for i in range(len(front[f'R_{landmark}_x'])):
      if landmark == 'foot':
          break
      plt.plot([front[f'R_{landmark}_x'][i], front[f'L_{landmark}_x'][i]],
              [front[f'R_{landmark}_y'][i], front[f'L_{landmark}_y'][i]],
              color='red', linestyle='-', linewidth=0.1)

  if previous_landmark == None:
    previous_landmark = landmark
  
  else:
    # 各点をつなぐ線を描画(身体部位に対応する直線)
    for i in range(len(front[f'R_{landmark}_x'])):
        plt.plot([front[f'R_{landmark}_x'][i], front[f'R_{previous_landmark}_x'][i]],
                [front[f'R_{landmark}_y'][i], front[f'R_{previous_landmark}_y'][i]],
                color='blue', linestyle='-', linewidth=0.1)
        
        plt.plot([front[f'L_{landmark}_x'][i], front[f'L_{previous_landmark}_x'][i]],
        [front[f'L_{landmark}_y'][i], front[f'L_{previous_landmark}_y'][i]],
        color='blue', linestyle='-', linewidth=0.1)
    
    previous_landmark = landmark

```

※右矢状面
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/3727ee06-23ec-732b-c78d-2d165ae7824a.png)

※左矢状面
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/057a1954-84cb-acb2-3867-eecac23df55b.png)

※前額面
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/267ec8e4-a9fa-39ff-c47f-d31ec3818a74.png)

## 各関節角度の算出
下記のように各角度を求める事ができます。

| 部位   | 運動方向        | 基本軸と移動軸                                       | 
| :----- | :-------------- | :--------------------------------------------------- | 
| 体幹   | 前傾/後傾       | 肩峰-大転子と平面座標の垂直軸のなす角                | 
| 股関節 | 屈曲/伸展   | 大転子-膝関節と腓骨頭-外踝のなす角                   | 
|        | 外転/内転       | 肩峰-上前腸骨棘と上前腸骨棘-膝蓋骨中央のなす角       | 
| 膝関節 | 屈曲/伸展       | 大転子-大腿骨外側顆と腓骨頭-外踝のなす角             | 
|        | 膝内反角度(FTA) | 上前腸骨棘-膝蓋骨中央と膝蓋骨中央-外内踝中央のなす角 | 
| 足関節 | 底屈/背屈       | 腓骨頭‐大腿骨外側顆と踵骨-第五中足骨のなす角        | 

### 角度算出用の関数作成

2点間を結ぶ直線とX軸とのなす角を算出する関数
```python
# 4点の座標から成る2直線のなす角を求める関数
def calculate_angle_between_lines(x1_1, y1_1, x2_1, y2_1, x1_2, y1_2, x2_2, y2_2):
    # 直線1の傾きを計算
    m1 = (y2_1 - y1_1) / (x2_1 - x1_1)
    
    # 直線2の傾きを計算
    m2 = (y2_2 - y1_2) / (x2_2 - x1_2)

    # 2つの直線のなす角を計算
    angle_rad = math.atan(abs((m2 - m1) / (1 + m1 * m2)))

    # ラジアンから度に変換
    angle_deg = math.degrees(angle_rad)

    return angle_deg
```


2点を結ぶ直線同士のなす角を求める関数
```python
# 2点のX,Y座標を与え2点を結ぶ直線とX軸とのなす角を求める関数
def calculate_angle(x1, y1, x2, y2):
    # 2点の座標からベクトルを計算
    dx = x2 - x1
    dy = y2 - y1

    # アークタンジェントを使用して角度を計算（ラジアン）
    angle_rad = math.atan2(dy, dx)

    # ラジアンから度に変換
    angle_deg = math.degrees(angle_rad)

    return angle_deg
```

角度算出をデータフレームの各レコードに使用し動画フレーム毎に関節角度、水平角度を算出する
```python
# 左股関節角度を算出する例
left_side['L_hip_angle'] = left_side.apply(lambda row: calculate_angle_between_lines(row['L_shoulder_x'], row['L_shoulder_y'], row['L_hip_x'], row['L_hip_y'],row['L_hip_x'], row['L_hip_y'], row['L_knee_x'], row['L_knee_y']), axis=1)

# 両肩を結んだ直線の水平角度を算出する例
front['shoulder_horizonal_angle'] = front.apply(lambda row: calculate_angle(row['R_shoulder_x'], row['R_shoulder_y'], row['L_shoulder_x'], row['L_shoulder_y']), axis=1)
```

各データフレームに角度データが追加された
- hip_angle(股関節角度)
- knee_angle(膝関節角度)
- ancle_angle(足関節角度)
- shoulder_horizonal_angle(両肩峰水平角度)
- asis_horizonal_angle(両ASIS水平角度)
- FTA(膝内反角度)


## 関節角度・水平角度の視覚化

各データフレームは踵接地の時間を0として同期しているため、歩行中のイベント時間のずれや左右差といった情報が視覚化できています。

もっと分かりやすく各イベントの時間帯にマーカーを付けてもいいのですが、それは検証や考察の段階で行う事にします。

### 矢状面の各関節角度
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/b9f7eb02-fcd7-0bd8-d34b-ba3fb72f8c7c.png)

### 前額面の各水平角度,FTA
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/601137/d45623ac-b5fe-3598-9d4b-55a069eb559f.png)

左右を比較した時に明らかに関節角度の変化や左右の水平角度に傾きがありそうです。
ここから角加速度やランドマークの移動速度、加速度などのパラメータを算出できます。
データの平滑化の処理も必要になりそうです。

# まとめ
- [x] 前額面、左矢状面、右矢状面の各ランドマークのプロット直線の描画による視覚化
- [x] 体幹、股関節、膝関節、足関節の各運動方向の角度算出
- [x] 両肩峰、ASIS、膝蓋骨中央の位置から各傾斜角度を算出

余談ですが、今回の結果を受けて自分の身体の不均等さにショックを受けています。

引き続き結果を整理していき私の歩行の何がどう悪いのか検証していきたいと思います。