motion - 演算法筆記

human motion🚧

human motion

人體運動。你自己全身上下動一動吧。

The content is directly taken from Prof. Tung-Wu Lu. Though I attribute it without permission. I try to hide that fact by rewriting the whole content. I am pretending that I never do anything evil. Keep in mind do not announce it, okay?

使用現成軟體處理人體運動

開源軟體OpenSim、MOKKA。

商業軟體Vicon Nexus、Visual3D。

專案倉庫SimTK。

檔案格式非常多種，維基百科整理了一份詳細列表。此處介紹其中兩種檔案格式。

C3D：人體表面座標、測量儀器訊號。流行於生物力學領域。資源請見生物力學國際協會維護的網站C3D.org。資料庫來自學術期刊scientific data，先找到論文，再到研究資料庫figshare搜尋論文名稱即可下載檔案。資料庫來自大學實驗室CMU、SFU、OSU，但是只含人體表面座標。想要瀏覽檔案可以用開源軟體MOKKA。

二進位檔。檔案包含參數區和資料區。資料區包含許多幀。一幀包含一個時刻人體表面座標（五種數值⨯座標數量）、k個時刻測量儀器訊號（一種數值接連k次⨯訊號種類）。測量儀器訊號取樣頻率是人體表面座標取樣頻率的k倍，時間得以對齊，成為一幀。如果不是整數倍，我不知道會怎樣。

figshare檔案極大卻無法個別下載。如果各位對於下載感到困擾，我這裡直接提供幾個個別檔案，方便大家做實驗。版權宣告請見論文連結。

S1_DP_04_M.c3d (1078KB)
投球與撿球。兩腳分別踩在兩塊測力板。
Three-dimensional motion capture data during repetitive
overarm throwing practice
https://www.nature.com/articles/sdata2018272
https://doi.org/10.6084/m9.figshare.c.4017808.v1

2014001_C1_01.c3d (773KB)
健康成人走過兩塊測力板。
A multimodal dataset of human gait at different walking
speeds established on injury-free adult participants
https://www.nature.com/articles/s41597-019-0124-4
https://doi.org/10.6084/m9.figshare.7734767.v8

tf01_0p6_01.c3d (4076KB)
膝上截肢者在跑步機上行走。
A kinematic and kinetic dataset of 18 above-knee amputees
walking at various speeds
https://www.nature.com/articles/s41597-020-0494-7
https://doi.org/10.6084/m9.figshare.c.4962305.v1

p10_trial8.c3d (3801KB)
健康人在跑步機上行走。
A biomechanics dataset of healthy human walking at various
speeds, step lengths and step widths
https://www.nature.com/articles/s41597-022-01817-1
https://doi.org/10.6084/m9.figshare.c.5897423.v1

BVH：人體表面座標、角色骨架。流行於計算機動畫領域。資料庫來自研究單位Bandai Namco Research、Ubisoft La Forge。

純文字檔。

使用Python處理人體運動

C3D檔案格式的知名函式庫是ezc3d、py-c3d。

ezc3d將資料重新排版，擴充為齊次座標，儲存於多維陣列。py-c3d遵循C3D檔案格式，每次取得一幀。關於陣列格式，ezc3d選擇array[dimension][point][time]適合離線計算，py-c3d選擇array[time][point][dimension]適合即時存取。離線計算、即時存取，其中一者必定經常cache miss，二者不可得兼。

我將C3D檔案轉換成純文字檔案，而且重新排版。我選擇方才沒提到的array[point][time][dimension]適合教學。

import ezc3d
file = '2014001_C1_01.c3d'
c3d = ezc3d.c3d(file, extract_forceplat_data=True);
f = open(file + '.txt', 'w')

point_name    = c3d['parameters']['POINT']['LABELS']['value']
analog_name   = c3d['parameters']['ANALOG']['LABELS']['value']
point_data    = c3d['data']['points']
residual      = c3d['data']['meta_points']['residuals']
camera_mask   = c3d['data']['meta_points']['camera_masks']
analog_data   = c3d['data']['analogs']
platform_data = c3d['data']['platform']

# point data / analog data
print('point name =\n', point_name)
print('point number = ', len(point_name))
print('analog name =\n', analog_name)
print('analog number = ', len(analog_name))
assert(point_data.shape[1] == len(point_name))
assert(analog_data.shape[1] == len(analog_name))
print('point sampling rate = ',
      c3d['parameters']['POINT']['RATE']['value'][0])
print('point sample number = ',
      point_data.shape[2])
print('analog sampling rate = ',
      c3d['parameters']['ANALOG']['RATE']['value'][0])
print('analog sample number = ',
      analog_data.shape[2])
assert(point_data.shape[1] == residual.shape[1])
assert(point_data.shape[2] == residual.shape[2])
assert(point_data.shape[1] == camera_mask.shape[1])
assert(point_data.shape[2] == camera_mask.shape[2])
print('point unit = ',
      c3d['parameters']['POINT']['UNITS']['value'])

# platform data
print('platform number = ',
      len(c3d['data']['platform']))
if (len(c3d['data']['platform']) > 0):
    print('1st platform unit = ',
          c3d['data']['platform'][0]['unit_position'])
    print('1st platform corners = ',
          c3d['data']['platform'][0]['corners'])
    print('1st platform origin = ',
          c3d['data']['platform'][0]['origin'])
    print('2nd platform unit = ',
          c3d['data']['platform'][1]['unit_position'])
    print('2nd platform corners = ',
          c3d['data']['platform'][1]['corners'])
    print('2nd platform origin = ',
          c3d['data']['platform'][1]['origin'])

# force/torque
    print('1st force plate force unit = ',
          c3d['data']['platform'][0]['unit_force'])
    print('1st force plate torque unit = ',
          c3d['data']['platform'][0]['unit_moment'])
#   print('1st force plate force data size = ',
#          c3d['data']['platform'][0]['force'].shape)
#   print('1st force plate torque data size = ',
#          c3d['data']['platform'][0]['moment'].shape)

# point name (e.g., bony landmark)
point_name = c3d['parameters']['POINT']['LABELS']['value']
f.write(f"{len(point_name)}\n")
for line in point_name:
    f.write(f"{line}\n")

# point data
sample_number = point_data.shape[2]
f.write(f"{sample_number}\n")

for n in range(point_data.shape[1]):    # each marker
    for t in range(point_data.shape[2]):# each sample
        x = float(point_data[0][n][t])  # x-coordinate
        y = float(point_data[1][n][t])  # y-coordinate
        z = float(point_data[2][n][t])  # z-coordinate
#       r = float(residual[0][n][t])    # -1.0 when miss
#       c = int8(camera_mask[0...7][n][t]) # -1 when miss
        f.write(f"{x:.2f} {y:.2f} {z:.2f} ")
    f.write("\n") # all samples of one marker form a line

# analog name (e.g., force plates, EMG)
analog_name = c3d['parameters']['ANALOG']['LABELS']['value']
f.write(f"{len(analog_name)}\n")
for line in analog_name:
    f.write(f"{line}\n")

# analog data
sample_number = analog_data.shape[2]
f.write(f"{sample_number}\n")

for n in range(analog_data.shape[1]):    # each device
    for t in range(analog_data.shape[2]):# each sample
        # force/torque/voltage
        value = float(analog_data[0][n][t])
        f.write(f"{value:.2f} ")
    f.write("\n") # all samples of one device form a line

# platform data
platform_number = len(platform_data)
f.write(f"{platform_number}\n")

for n in range(len(platform_data)):      # each plate
    origin = platform_data[n]['origin']
    for d in range(len(origin)):         # each dimension
        f.write(f"{origin[d]:.2f} ")

corners = platform_data[n]['corners']
    for i in range(corners.shape[1]):    # each corner
        for d in range(corners.shape[0]):# each dimension
            f.write(f"{corners[d][i]:.2f} ")
    f.write("\n") # one platform form a line

f.close()

# https://c3d.readthedocs.io/en/master/generated/c3d.Reader.html
import c3d
file = '2014001_C1_01.c3d'
reader = c3d.Reader(open(file, 'rb'))
f = open(file + '.txt', "w")

# information
print('point name =\n', reader.point_labels)
print('analog name =\n', reader.analog_labels)
print('point sampling rate = ', reader.point_rate)
print('analog sampling rate = ', reader.analog_rate)
print('first frame number = ', reader.first_frame)
print('last frame number = ', reader.last_frame)
print('parameter blocks = ', reader.parameter_blocks())

# point name (e.g., bony landmark)
point_name = [s.strip() for s in reader.point_labels]
f.write(f"{len(point_name)}\n")
for line in point_name:
    f.write(f"{line}\n")

# point data
frame_number = reader.last_frame - reader.first_frame + 1
f.write(f"{frame_number}\n")

for n in range(len(point_name)):
    for frame_no, point_data, analog_data in reader.read_frames():
        assert len(point_data) == len(point_name)

x = float(point_data[n][0])    # x-coordinate
        y = float(point_data[n][1])    # y-coordinate
        z = float(point_data[n][2])    # z-coordinate
        r = float(point_data[n][3])    # residual
#       c = int8(point_data[n][4])     # camera mask
        if (r == -1.0): x = y = z = float('nan')
        f.write(f"{x} {y} {z} ")
    f.write("\n") # all samples of one marker form a line

# analog name (e.g., force plates, EMG)
analog_name = [s.strip() for s in reader.analog_labels]
f.write(f"{len(analog_name)}\n")
for line in analog_name:
    f.write(f"{line}\n")

# analog data
multiplier = int(reader.analog_rate / reader.point_rate)
f.write(f"{frame_number * multiplier}\n")

for n in range(len(analog_data)):
    for frame_no, point_data, analog_data in reader.read_frames():
        assert len(analog_data) == len(analog_name)
        assert len(analog_data[n]) == multiplier

for k in range(multiplier):
            # force/torque/voltage
            value = float(analog_data[n][k])
            f.write(f"{value} ")
    f.write("\n") # all samples of one device form a line

f.close()

BVH檔案格式的函式庫仍在開發中。

https://github.com/20tab/bvh-python
https://github.com/Wasserwecken/bvhio
https://github.com/rltonoli/bvhsdk

使用HTML與JavaScript處理人體運動

JavaScript本身沒有處理人體運動的函式庫。目前也沒有現成的函式庫。必須自己造輪子。

至於我是利用Python函式庫將C3D檔案轉換成純文字檔案，再利用JavaScript讀檔。

順帶一提，網路上可以找到將BVH檔案匯入three.js和Babylon.js的程式。

human motion capture🚧

human motion capture

運動捕捉。測量人體運動資訊：位置、方向、……。

human pose estimation與human shape estimation：測量人體肢段方位、測量人體表面座標。已有多種工具，請見產業報告。

一、光學動作捕捉optical motion capture：
　　利用多臺攝影機估計人體表面座標。
　　原理是三角測量triangulation。
　　領域是攝影測量photogrammetry。
二、慣性動作捕捉inertial motion capture：
　　利用穿戴式裝置估計人體肢段方位。
　　例如傳感器（加速規、陀螺儀、磁力計）、機械設備。
　　原理是各種物理定律。
　　領域是微機電系統microelectromechanical system。
三、無標記動作捕捉markerless motion capture：
　　拍攝影片，利用深度學習估計人體肢段方位、人體表面座標。
　　原理是深度神經網路deep neural network。
　　領域是計算機視覺computer vision。

human force estimation：測量人體施力。人體發出的力，受到皮肉緩衝，才是儀器測量結果。已有多種工具。

回、測力板force plate：測量力。
　回、一維測力板1D force plate：
　　　測量壓力。二維陣列，一格一個數值。
　　　腳踩下去，平面上處處有壓力。
　回、三維測力板3D force plate：
　　　測量力、力矩。總共六個數值。主流作法。
　　　測力板四角落有感測器，測得四道力。
　　　測力板內建韌體自動計算測力板中心、合力、合力矩。
　　　得到地面反作用力ground reaction force。
回、動力計dynamometer：測量力矩。
　回、等張肌力儀isotonic dynamometer。
　回、等長肌力儀isometric dynamometer。
　回、等速肌力儀isokinetic dynamometer。
　　　得到肌群施力muscle group force。

human physiological signal measurement：測量人體生理訊號。我只知道肌電圖儀。

回、電圖儀electrography：測量電壓。
　回、肌電圖儀electromyography (EMG)：
　　　得到肌肉激發量muscle excitation。
　　　一連串肌肉激發量訊號，稱作肌電圖electromyogram (EMG)。
回、光體積變化儀photoplethysmography：測量體積。
　回、脈搏血氧儀pulse oximetry：
　　　其中包含了血容積變化量blood volume change。
　　　一連串血容積變化量訊號，
　　　稱作脈搏血氧容積變化圖pulse oximetry plethysmogram (PLETH)。

optical motion capture

光學動作捕捉。需要兩種設備。

一、marker：人體多處貼上標記。

生物力學國際協會已經提出標準規範，告訴大家如何描述人體運動：每個標記的黏貼地點、每個肢段的座標系統、每個關節的旋轉變換（歐拉角三個轉軸）。主要源自這兩篇文章：

ISB recommendation on definitions of joint coordinate system
of various joints for the reporting of human joint motion
--part I: ankle, hip, and spine [Wu et al. 2002]
https://pubmed.ncbi.nlm.nih.gov/11934426/
https://media.isbweb.org/images/documents/standards/isb_jcs_part_i.pdf

--part II: shoulder, elbow, wrist and hand [Wu et al. 2005]
https://pubmed.ncbi.nlm.nih.gov/15844264/
https://media.isbweb.org/images/documents/standards/isb_jcs_part_ii.pdf

下面是黏貼地點圖片、地點名稱清單。圖文截自這篇學術論文。更加詳細的圖文請見Visual3D軟體說明文件。每位研究員習慣自行增減黏貼地點、自行調整地點名稱，未必完全符合標準規範。我知道你想說什麼，這不是人體穴位圖。

name | segment/joint | location
------------------------------------------------------
RFHD | Head          | Apprx. over the Right Temple
RBHD | Head          | Right Back Of The Head
STRN | Trunk         | Suprasternal Notch
XIPH | Trunk         | Xiphoid Process
C7   | Trunk         | 7th Cervical Vertebra
T8   | Trunk         | 8th Thoracic Vertebra
RASI | Pelvis        | Right Ant Superior Iliac Crest
RPSI | Pelvis        | Right Post Superior Iliac Crest
RILC | Pelvis        | Right Iliac Crest
RSHO | Right Humerus | Acromion Process
RSTS | Right Humerus | Trigonum Scapulae
RSAI | Right Humerus | Angulus Inferior
RSAA | Right Humerus | Angulus Acromialis
RELL | Right Elbow   | Lateral Epicondyle
RELM | Right Elbow   | Medial Epicondyle
RWRL | Right Wrist   | Lateral Point of Radial Styloid
RWRM | Right Wrist   | Medial Point of Ulnar Styloid
R3MC | Right Hand    | 3rd Metacarpal Head
R5MC | Right Hand    | 5th Metacarpal Head
RGT  | Right Femur   | Femur Greater Trochanter
RFLE | Right Knee    | Femur Lateral Epicondyle
RFME | Right Knee    | Femur Medial Epicondyle
RFAL | Right Ankle   | Fibula Lateral Malleolus
RTAM | Right Ankle   | Tibia Medial Malleolus
RFM2 | Right Foot    | 2nd Metatarsal Head
RFM5 | Right Foot    | 5th Metatarsal Head
RHEL | Right Foot    | Posterior of Calcaneus
註：省略了左半邊。省略了肢段標記群。

marker分為兩種類型。

一、主動發光：不同顏色的LED燈。利用顏色分辨不同標記。
二、被動反光：額外提供紅外線光源。設計演算法分辨不同標記。

二、camera：使用攝影機拍攝人體運動求得影片，使用軟體讀取影片求得每個時刻每個標記的三維座標。

攝影測量領域已有一套解決方案。設備例如Vicon。演算法請見本站文件camera tracking。

camera calibration需要兩種設備。

一、靜態校正：放置三角板L-frame，於測力板force plate轉角。
　　以便稍後訂立世界座標系。
二、動態校正：揮動T型魔法棒wand。功效如同西洋棋盤。
　　以便確認拍攝空間範圍：視情況重新調整相機位置、方向。
　　（支架高度、鉛直轉角、水平轉角。）
　　以便得到相機外在參數：位置、方向。
　　以便得到相機內在參數：焦距、扭曲係數。

https://www.semanticscholar.org/paper/Photogrammetric-Measurements-for-Rigid-Body-in-Low-Atyam/bdff6e8b19932e35eb5a73eb0b70830fffd59744/figure/4

inertial motion capture

慣性動作補捉。主流設備是慣性測量單元IMU。

IMU的特色是重量輕、體積小。貼在身上、衣服上，不影響人體運動。適合用於室外活動。

IMU內含三種感測器：加速規（加速度）、陀螺儀（角加速度）、磁力計（重力加速度方向）。利用三種感測器的數值，可以估計自身位置和方向。

IMU內含無線傳輸裝置，藉此即時獲取三種感測器的數值。

IMU經常產生誤差。經典問題gyroscope drift。解法是sensor fusion。融合加速規和陀螺儀，套用互補濾波器，計算角度位置。

complementary filter電路圖
https://www.vectornav.com/resources/inertial-navigation-primer/math-fundamentals/math-filtering
https://gunjanpatel.wordpress.com/2016/07/07/complementary-filter-design/
https://pmc.ncbi.nlm.nih.gov/articles/PMC7998881/figure/sensors-21-01937-f003/

IMU容易受到干擾。經典問題accelerator shaking。解法是signal denoising。演算法是Kalman filter。

Kalman filter電路圖

markerless motion capture

無標記動作捕捉。不需要測量設備。大家正在研究當中。

作法是deep learning。事前拍攝大量影片、標記肢段位置、訓練模型參數。事後拍攝任意影片，套用模型參數，計算肢段位置。

例如使用兩台攝影機，固定位置，同時拍攝影片。用deep learning找到二維肢段位置，用triangulation重建三維肢段位置。

1D force plate

一維測力板。測量矩形區域每一處的壓力，得到一個矩陣。

進一步計算壓力中心center of pressure：位置乘以壓力，求平均數。計算公式如同重心。

測量原理主要分為兩種。

一、彈性材料（應變規）：施力大小決定伸縮大小決定電阻大小。
　　缺點是熱脹冷縮、彈性疲乏。
　　適合用於靜態人體。例如醫療診所。
二、壓電材料（壓電感測器）：施力大小決定電壓大小。主流作法。
　　適合用於動態人體。例如戶外運動。

3D force plate

三維測力板。測量力與力矩，得到兩個三維向量。

測力板四角落安裝感測器，得到兩道橫向力、兩道縱向力、四道鉛直力。事先指定參考點（測力板中心），經過計算得到力與力矩。詳情請見測力板說明文件。

進一步計算合力中心（但是大家仍然沿用壓力中心center of pressure這個名稱）。合力中心就是力矩等於零的地方。請見測力板說明文件。

這兩份測力板說明文件，來自生物力學國際協會。

原施力點（視作原點）p₀ = 0、力F、力矩M。
新施力點（壓力中心）p、力F、力矩0、力臂r。
⎰ M = r × F
⎱ r = p - p₀
p是未知數。p₀ M F是已知數。
藉由力矩公式M = p × F，求得壓力中心。
⎰ pₓ = (p₞ Fₓ - M₝) / F₞
⎱ p₝ = (Mₓ + p₞ F₝) / F₞
其中p₞由測力板規格書取得。p₞數值根據感測器材料彈性而定。
p₞數值極小，實務上經常直接視作0。
另外，因為pₓ p₝可能是零，所以不能當作分母。

計算原理是equivalent force couple system：力和力矩可以挪動到物體任意位置，甚至是旋轉中心。

回、力矩公式。旋轉中心p₀、施力點p、力F、力矩M、力臂r。
　　⎰ M = r × F
　　⎱ r = p - p₀
回、力挪至旋轉中心：任意點的力F₁＝旋轉中心的力F₀暨力矩M₀
　　⎧ F₀ = F₁
　　⎨ M₀ = r₁ × F₁
　　⎩ r₁ = p₁ - p₀
回、力暨力矩挪出旋轉中心：假設目的地需要補上力矩M₂。
　　M₀ = r₂ × F₂ + M₂
　　得到公式：
　　⎧ F₂ = F₀
　　⎨ M₂ = M₀ - r₂ × F₂
　　⎩ r₂ = p₂ - p₀
　　教科書將r₂定義成反方向（將目的地p₂想像成參考點）：
　　⎧ F₂ = F₀
　　⎨ M₂ = M₀ + r₂ × F₂
　　⎩ r₂ = p₀ - p₂
　　如果目的地恰是原點：
　　⎰ F₂ = F₀
　　⎱ M₂ = M₀ + p₀ × F₀   since p₂ is zero vector

isokinetic dynamometer

等速肌力儀。我沒有學會。

肌肉運動機制可以分成三類，於是測量儀器隨之分成三類。

http://www.hkpe.net/hkdsepe/human_body/muscular_contraction_types.htm

electromyography

肌電圖儀。測量肌肉激發量muscle excitation。

一連串肌肉激發量，形成訊號，稱作肌電圖electromyogram。

我們需要訊號功率。四個步驟。

https://ouhsc.edu/bserdac/dthompso/web/pk/emg/emg.htm
一、整流rectification：取絕對值。
二、濾波filtering：擷取20-200Hz的訊號。沒有固定說法。
三、積分integration：求平方和，得到能量。再除以訊號長度，得到功率。
四、放大amplification：估計阻抗大小，還原原始能量。

測量手段分為兩種類型。

一、侵入式intramuscular EMG：探針needle、魚鉤fine wire。
二、非侵入式surface EMG：拋棄式電極片disposable electrode。主流作法。

interpreting signal amplitudes in surface electromyography studies in sport and rehabilitation sciences
https://pmc.ncbi.nlm.nih.gov/articles/PMC5758546/

測量步驟請見專著《The ABC of EMG》。

human motion analysis🚧

human motion analysis

骨骼skeleton、關節joint、肌肉muscle、神經nerve。

gait analysis

步態分析。觀察走路的姿勢。

some videoes of cerebrovascular accident gait analysis

想要觀察步態，可以利用地面反作用力。從函數曲線可以初步判斷慢走、奔跑、行動不便，乃至年齡、身高、體重，甚至鑑別身分。

測力板xyz曲線圖，M型圖之類的。
https://www.tekscan.com/blog/medical/how-pressure-mapping-complements-force-measurement-in-gait-analysis
https://www.researchgate.net/figure/Treadmill-and-GRF-plate-measurement-parameters-of-the-gait-cycle-GRF-ground-reaction_fig3_236460965

想要觀察步態，可以利用肌肉活化量。從函數曲線可以進一步判斷行走的每個階段，乃至體能、疲勞、傷害，甚至鑑別疾病。

gait cycle全身肌肉活化圖，就是影印十幾份到處張貼每天背誦的那個圖。
https://cdn.shopify.com/s/files/1/1281/4359/files/E8284_177art.png
https://us.humankinetics.com/blogs/excerpt/see-how-the-muscles-work-to-create-ambulation

想要量化步態，可以利用每個關節的旋轉角度、旋轉角速度、旋轉角加速度。所有關節的旋轉角度，決定了人體姿勢，一一對應。

回、肢段：座標系統。
　回、肢段座標系統。
　　　根據標準規範，取得對應marker。
　　　一個原點（視作位移向量）、三個座標軸（視作旋轉矩陣）。
　　　    ⎡ x ⎤       ⎡ R₁₁ R₁₂ R₁₃ ⎤   ⎡ |  |  | ⎤
　　　t = ⎢ y ⎥   R = ⎢ R₂₁ R₂₂ R₂₃ ⎥ = ⎢ X  Y  Z ⎥
　　　    ⎣ z ⎦       ⎣ R₃₁ R₃₂ R₃₃ ⎦   ⎣ |  |  | ⎦
　回、旋轉矩陣是正規正交矩陣，導致反矩陣等於轉置矩陣。
　　　R⁻¹ = Rᵀ
　回、全域座標系（世界座標系）、局部座標系（肢段座標系）的轉換。
　　　for position vector p_local and p_global,
　　　p_local = Rᵀ (p_global - t)
　　　p_global = R p_local + t
　　　for direction vector v_local and v_global,
　　　v_local = Rᵀ v_global
　　　v_global = R v_local
回、關節：角度。
　回、關節旋轉矩陣。《direction cosine matrix》
　　　關節是相鄰肢段交點。
　　　從身體近端肢段去看身體遠端肢段的相對運動（近端為準）。
　　　遠端三個座標軸，分別從遠端座標系轉換到近端座標系，形成旋轉矩陣。
　　　R_joint = R_proximalᵀ R_distal
　回、關節角度。《angles from rotation matrix》
　　　旋轉矩陣改成歐拉角，三個轉軸分別是三個標準座標軸。
　　　intrinsic rotations of Cardan angles x-y'-z"
　　　https://stackoverflow.com/questions/15022630/
　　　let R_joint = R₞(ψ₞) R₝(ψ₝) Rₓ(ψₓ)
　　　⎧ ψₓ = tan⁻¹(R₂₁/R₁₁)
　　　⎨ ψ₝ = sin⁻¹(-R₃₁)
　　　⎩ ψ₞ = tan⁻¹(R₃₂/R₃₃)
回、關節：角速度、角加速度。
　回、旋轉矩陣一次微分、旋轉矩陣二次微分。《finite difference》
　　　       R[t+1] - R[t]
　　　Ṙ[t] = —————————————
　　　             Δt

　　　       R[t+1] - 2 R[t] + R[t-1]
　　　R̈[t] = ————————————————————————
　　　                (Δt)²
　回、角速度、角加速度。《time derivative of rotation matrix》
　　　https://physics.stackexchange.com/questions/293037/
　　　            ⎡  0  -ω₞  ω₝ ⎤  
　　　Ṙ = ω × R = ⎢  ω₞  0  -ωₓ ⎥ R = ΩR
　　　            ⎣ -ω₝  ωₓ  0  ⎦  
　　　       ⎡  0  -ω₞  ω₝ ⎤
　　　=> Ω = ⎢  ω₞  0  -ωₓ ⎥ = ṘRᵀ
　　　       ⎣ -ω₝  ωₓ  0  ⎦

　　　R̈ = (Ω̇R + ΩṘ) = (Ω̇R + ΩΩR)
　　　=> Ω̇ = (R̈ - ΩΩR)Rᵀ = R̈Rᵀ - ΩΩ

角速度、角加速度還可以更精確【尚待確認】。一、離散數列有限差分改成連續函數微分。為此需要訊號重建。二、旋轉矩陣改成歐拉角／四元數。

回、關節：角速度、角加速度。
　回、旋轉矩陣改寫成歐拉角／四元數。
　　　歐拉角三個轉軸：根據標準規範。
　　　歐拉角三個轉角：維基百科提供了所有12種轉軸順序的公式。
　　　例如膝關節intrinsic rotations of Cardan angles z-x'-y"
　　　let R_joint = R₞(ψ₞) Rₓ(ψₓ) R₝(ψ₝)
　　　⎧ ψ₞ = tan⁻¹(-R₁₂/R₂₂)
　　　⎨ ψₓ = sin⁻¹(R₃₂)
　　　⎩ ψ₝ = tan⁻¹(-R₃₁/R₃₃)
　　　例如肩關節intrinsic rotations of Euler angles y-x'-y"
　　　let R_joint = R₝(α) Rₓ(β) R₝(γ)
　　　⎧ α = tan⁻¹(R₁₂/R₃₂)
　　　⎨ β = cos⁻¹(R₂₂)
　　　⎩ γ = tan⁻¹(-R₂₁/R₂₃)
　回、所有時刻的歐拉角／四元數，建立連續函數，以便微分。三種方法。
　　　《signal reconstruction》
　　回、樣條內插spline interpolation。主流作法。
　　回、局部迴歸local regression。
　　回、傅立葉級數Fourier series。
　回、歐拉角／四元數的一次微分、二次微分。
　回、角速度、角加速度。
　　　例如膝關節intrinsic rotations of Cardan angles z-x'-y"
　　　⎡ ωₓ ⎤   ⎡  cosψ₝   0 -cosψₓsinψ₝ ⎤ ⎡ ψ̇ₓ ⎤
　　　⎢ ω₝ ⎥ = ⎢      0   1      -sinψₓ ⎥ ⎢ ψ̇₝ ⎥
　　　⎣ ω₞ ⎦   ⎣ -sinψ₝   0  cosψₓcosψ₝ ⎦ ⎣ ψ̇₞ ⎦

　　　⎡ αₓ ⎤   ⎡ -ψ̇₝sinψ₝  0  ψ̇ₓsinψₓsinψ₝-ψ̇₝cosψₓcosψ₝ ⎤ ⎡ ψ̇ₓ ⎤
　　　⎢ α₝ ⎥ = ⎢        0  0                    ψ̇ₓcosψₓ ⎥ ⎢ ψ̇₝ ⎥
　　　⎣ α₞ ⎦   ⎣  ψ̇₝cosψ₝  0 -ψ̇ₓsinψₓcosψ₝-ψ̇₝cosψₓsinψ₝ ⎦ ⎣ ψ̇₞ ⎦

　　　         ⎡  cosψ₝   0 -cosψₓsinψ₝ ⎤ ⎡ ψ̈ₓ ⎤
　　　       + ⎢      0   1      -sinψₓ ⎥ ⎢ ψ̈₝ ⎥
　　　         ⎣ -sinψ₝   0  cosψₓcosψ₝ ⎦ ⎣ ψ̈₞ ⎦

生物力學領域習慣使用歐拉角。首倡是Grood–Suntay angles，定義了膝關節歐拉角三個轉軸。

a joint coordinate system for the clinical description of
three-dimensional motions: application to the knee
https://pubmed.ncbi.nlm.nih.gov/6865355/
https://www.researchgate.net/publication/16334377

歐拉角範圍±180°-±90°-±180°。人體每個關節，都至少存在一個轉角無法抵達±90°。妥善規劃轉軸順序，即可符合轉角範圍，避免gimbal lock。肩關節是唯一例外，目前無解。標準規範採用傳統歐拉角（一三轉軸相同），第二轉角範圍超過±90°。有人建議使用其他轉軸順序，以便盡量減少錯誤。

comparison of glenohumeral motion using different
rotation sequences
https://pmc.ncbi.nlm.nih.gov/articles/PMC3042544/

手工插圖
https://www.researchgate.net/figure/The-shoulder-rotations-are-Euler-angles-with-a-YXY-order-The-shoulder-elevation-plane-is_fig2_266104104

歐拉角有嚴重缺陷。一、歐拉角只適合描述旋轉變換，不適合描述旋轉軌跡。三個轉軸依序自旋、三個轉軸同時自旋，旋轉變換不同，旋轉軌跡也不同。歐拉角只能三個轉軸依序自旋。現實世界需要三個轉軸同時自旋。二、即便是等速旋轉，所有時刻的旋轉矩陣9數值／歐拉角3數值／四元數4數值也不會是一次函數／多項式函數，不適合一次迴歸／多項式迴歸。

幸運的是，當今設備足夠先進。取樣頻率極高、取樣時距極短，導致角速度誤差極小。無論哪種方法，計算結果都差不多。

stabilogram analysis

穩定圖分析。所有時刻的center of pressure收集起來，做個principal component analysis，可以知道阿嬤腳麻的情況。

Rambling and Trembling decomposition
https://nbviewer.org/github/BMClab/BMC/blob/master/notebooks/IEP.ipynb
https://www.researchgate.net/publication/11843491

回、測力板：座標系、施力點。
　　讀取C3D檔案，取得每塊測力板的每個頂點座標，建立座標系。
　　原點：測力板中心。測力板四個頂點座標平均值。
　　p = (p₁ + p₂ + p₃ + p₄) / 4
　　座標軸：四個頂點依序對應四個象限。頂點相減可得座標軸。
　　https://www.c3d.org/HTML/Documents/forceplatformcorners.htm
　　⎧ X = normalize(p₁ - p₂)
　　⎨ Y = normalize(p₁ - p₄)
　　⎩ Z = X × Y
　　施力點：新原點到測力板中心的位移量t。
　　https://www.c3d.org/HTML/Documents/forceplatformorigin.htm
　　p' = p - t
回、多塊測力板：力、力矩。《equivalent force couple system》
　　讀取C3D檔案，取得每塊測力板的力、力矩、施力點。
　回、回顧力矩公式。旋轉中心p₀、施力點p、力F、力矩M、力臂r。
　　　⎰ M = r × F
　　　⎱ r = p - p₀
　回、回顧座標系轉換公式。座標軸R₁ R₂、原點p₁ p₂。
　　　⎰ T = R₁ᵀ R₂
　　　⎱ r = p₂ - p₁
　回、多塊測力板，一律轉換為世界座標系。
　回、所有力從旋轉中心挪至原點，進行加總。
　　　⎧ Tᵢ = R₀⁻¹ Rᵢ     (R₀ is identity matrix)
　　　⎨ rᵢ = pᵢ - p₀     (p₀ is zero vector)
　　　⎪ F₀ = sum { Tᵢ Fᵢ }
　　　⎩ M₀ = sum { Tᵢ Mᵢ + rᵢ × (Tᵢ Fᵢ) }
回、多塊測力板：施力點（壓力中心center of pressure）。
　　力從原點挪至零力矩處。
　　⎰ pₓ = (p₞ Fₓ - M₝) / F₞
　　⎱ p₝ = (Mₓ + p₞ F₝) / F₞
　　正確名稱應是合力中心。阿就歷史共業。

joint force estimation

關節力量估計。計算人體關節的力、力矩。

專著《Kinematics of Human Motion》。

專著《Winter's Biomechanics and Motor Control of Human Movement》。

主流作法是inverse kinematics。

已知
一、肢段：施力點、力、力矩
二、肢段遠端關節：施力點、力、力矩
三、肢段近端關節：施力點
求得
三、肢段近端關節：力、力矩

一、肢段：端點。《joint center estimation》
　　關節中心作為肢段端點。估計關節中心，兩種方式。
　回、加權平均數。
　　　大多數關節：肢端末梢附近兩個marker的中點。
　　　髖關節：有許多種方式。請自行查閱論文。
　　　https://pmc.ncbi.nlm.nih.gov/articles/PMC5121588/
　　　肩關節：由兩個關節組成。marker隨皮膚移動，至今無法精算，除非打骨釘。
　回、chord function [Vicon Plug-in Gait®]
　　　https://help.vicon.com/space/Nexus216/11601572/
二、肢段：形狀。《body segment parameter estimation》
　　以少量參數描述形狀。估計肢段參數，三種方式。
　回、大體解剖。
　回、醫學影像：電腦斷層掃描computed tomography scan。
　回、數學模型：圓錐台conical frustum。
三、肢段：質心、移動慣量、轉動慣量。
　　採用Dempster's body segment parameter。
　　http://holmeslab.ca/csb-scb/Software/dempster.pdf
　回、質心：肢段兩個端點的一次內插。肢段長度比例請查表。
　回、移動慣量：體重乘上肢段重量比例。肢段重量比例請查表。
　回、轉動慣量：肢段重量乘上慣性半徑平方。慣性半徑請查表。【尚待確認】
　　採用圓錐台conical frustum。
　　http://www.oandplibrary.org/al/pdf/1964_01_044.pdf
四、肢段：速度、加速度、角速度、角加速度。
　　手法如同先前小節gait analysis。
　回、肢段質心的速度、加速度：
　　　肢段質心，建立連續函數，微分。
　　　剛體處處移動速度皆不同。必須針對質心。
　回、肢段質心的角速度、角加速度：
　　　肢段座標系，視作旋轉矩陣，轉換成歐拉角／四元數，
　　　建立連續函數，微分，轉換成三軸旋轉。
　　　剛體處處旋轉速度皆相同。座標系原點在哪都無所謂。
五、肢段：力、力矩、施力點（質心）。
　　⎧ F = ma + mg
　　⎨ M = Iω
　　⎩ o = mass center
六、測力板：力、力矩、施力點。
　　手法如同先前小節stabilogram analysis。
　回、讀取C3D檔案，取得每塊測力板的每個頂點座標，建立座標系。
　　　讀取C3D檔案，取得每塊測力板的力、力矩、施力點。
　回、一腳踩在多塊測力板。一律轉換為世界座標系，再加總。
七、足底：力、力矩、施力點（壓力中心center of pressure）。
　　手法如同先前小節stabilogram analysis。
　回、正確名稱應是合力中心。阿就歷史共業。
　回、為了區隔左腳和右腳的地面反作用力，
　　　左腳和右腳必須踩在不同測力板。
　　　為了測量自然的步態，
　　　不能告訴受試者要注意地面調整步伐。
　　　實務上是請受試者重複走上百次，等待奇蹟發生，曠日廢時。
　　　等你發明更好的測力板。
　　　其實我覺得可以並排測力板、對準中線走。
八、關節：力、力矩。《inverse kinematics》
　　逆向運動學。從下往上逐步推導每個關節的力、力矩。
　回、肢段相依表。
　　　  | segment | distal   | proximal  | net force
　　　  |         | joint    | joint     |
　　　  | (known) | (known)  | (unknown) | (equation)
　　　--------------------------------------------------
　　　1 | foot f  | ground g | ankle a   | f + g + a = 0
　　　2 | shank s | ankle a  | knee k    | s + a + k = 0
　　　3 | thigh t | knee k   | hip h     | t + k + h = 0
　　　注意到，每當更換肢段，關節的力和力矩必須變號。
　　　例如1和2的ankle a相差一個負號。
　回、力和力矩的詳細公式如下。
　　　範例：已知地面ground g、足肢段foot f，求得踝關節ankle a。
　　　force:
　　　F_f + F_g + F_a = 0
　　　F_a = -(F_f + F_g)
　　　torque:
　　　M_f + r_f × F_f + M_g + r_g × F_g + M_a + r_a × F_a = 0
　　　M_a = -(M_f + r_f × F_f + M_g + r_g × F_g + r_a × F_a)
　　　_f: foot       F: force
　　　_g: ground     M: torque
　　　_a: ankle      r: arm of force
　回、力臂參考點是世界座標系原點、或者人體重心、或者壓力中心。

muscule force estimation

肌肉力量估計。計算人體肌肉的力、力矩。

專著《Biomechanics of Movement: the Science of Sports, Robotics, and Rehabilitation》。

一種作法是constrained optimization。

肌肉數量遠多於肢段數量。不同肌力分配可能導致相同姿勢。未知變數遠多於已知變數。欠定方程組通常導致無限多組解。

欠定方程組可以改寫成最佳化問題，但是我們不知道什麼是好的objective function、好的constraint。這是open problem。

肌肉疲勞將自動改變肌力分配。詳細機制我沒有學會。

人體腳部肌肉圖

constraint optimization example

另一種作法是muscle activation dynamics。

肌肉活化量、肌肉長度、肌肉速度，決定肌肉力量。

激發-活化圖
excitation–activation

長度-速度-肌力圖

利用EMG測量肌肉激發量，進而估計肌肉活化量。但是目前沒有設備可以同時測量肌肉長度和肌肉速度。這是open problem。

EMG-driven models
https://pmc.ncbi.nlm.nih.gov/articles/PMC1357215/

forward / inverse
https://jneuroengrehab.biomedcentral.com/articles/10.1186/s12984-023-01253-1

human motion control🚧

human motion control

QWOP

the man who lost his body

標題應該不是觸覺touch而是本體感覺proprioception。

也許介紹一下外骨骼。

human measurement（anthropometry）

人體測量。測量人體外表。

身高、體重、頭圍、胸圍。

human dissection（anatomy）

人體解剖。測量人體內裡。

頭、胸、腹、上臂、前臂、手、大腿、小腿、足。

human body part
https://www.researchgate.net/figure/Body-scan-post-segmentation-A-total-of-16-body-segements-are-shown_fig18_322919256
https://pmc.ncbi.nlm.nih.gov/articles/PMC8730461/figure/pone.0262296.g003/

解剖學姿勢。全身直立躺平，手心朝向正面。

anatomical plane
coronal plane    冠狀面。人體切成前後兩段的面。
sagittal plane   矢狀面。人體切成左右兩段的面。
transverse plane 橫斷面。人體切成上下兩段的面。

anatomical location
https://zh.wikipedia.org/wiki/解剖學方位

human movement system

人體動作系統。活動的器官的集合。

本文主要討論三種系統：骨骼、肌肉、神經。

骨骼系統　skeletal system
神經系統　nervous system
肌肉系統　muscular system

大家視討論對象重組詞彙。

肌肉骨骼系統　　　musculoskeletal system
神經肌肉系統　　　neuromuscular system
神經肌肉骨骼系統　neuromusculoskeletal system

神經肌肉骨骼系統的構造。

骨頭bone：🦴
骨骼skeleton：所有骨頭，相互銜接。
關節joint：骨頭與骨頭的接口。
肌肉muscle：連接到骨頭。控制骨頭位置。宛如彈簧。
神經nerve：連接到肌肉。控制肌肉伸縮。宛如enter鍵。

軟骨cartilage：聯繫兩塊相鄰骨頭，用來撐開關節。材質是水、有機鹽、骨細胞。
韌帶ligament：聯繫兩塊相鄰骨頭，用來拉緊關節。材質是肌細胞。
肌腱tendon：肌肉末梢。聯繫肌肉與骨頭。材質是肌細胞。
突觸synapse：神經末梢。連繫神經與肌肉。

肌腱、韌帶、軟骨、骨頭，通通都是黏彈性材料。

viscoelastic material
1. elasticity 彈性。正比於應變長度。跟時間無關的恢復力。
2. viscosity  黏性。正比於應變速度。跟時間有關的恢復力。

viscoelastic mechanism
1. creep 潛變。負重、最初迅速變形、持續緩慢變形。坐沙發。
2. stress relaxation 應力鬆弛。持續緩慢恢復原形、應力下降。
                     橡皮筋圈久慢慢變鬆。

human motor system

人體運動系統。神經控制暨肌肉運動的機制。抱歉我沒有學會。

https://en.wikipedia.org/wiki/Motor_system
https://en.wikipedia.org/wiki/Motor_control

gait simulation

ballistic walking model
passive dynamic walking model
wobbling mass model (inverse kinematics)

joint simulation

關節位移運動有三種模式。

joint movement
1. gliding      滑動，關節橫向位移
2. distraction  分散，關節拉開
3. compression  壓縮，關節靠近
https://www.slideshare.net/slideshow/translatoric-joint-playpptx/259499242
https://docs.univr.it/documenti/OccorrenzaIns/matdid/matdid453153.pdf

關節旋轉運動有三種方向。

joint rotation convention
1. flexion-extension  上肢下肢彎曲
2. abduction-adduction 下肢左右搖擺
3. internal-external rotation 下肢自旋

關節自由度測試。

degree of freedom / laxity test
1. stability / mobility   不能移動、能移動
2. unresisted / resisited 移動不受抵抗、受到抵抗（韌帶能用）
3. passive / active       自主意志不能移動、能移動

膝關節：十字韌帶ACL與PCL繫住大腿骨和小腿骨。由於繫住位置不變，而且韌帶長度幾乎不變，形成平面四連桿planar four-bar linkage。顳顎關節也是相同機制。

髕骨：膝關節前方一塊骨頭。膝關節彎曲，產生向前的力。髕骨用來承擔此力，以便提升大小腿肌肉對膝關節的控制靈敏度。髕骨有兩個知名機制。

knee extensor mechanism：
膝關節彎曲，髕骨被大腿骨舉起。
髕骨相對於膝關節重心（大腿骨內部）的力臂變長，更省力。
讓人輕鬆蹲下、輕鬆站起。

screw–home mechanism：
膝關節彎曲，髕骨挪往身體內側。
膝關節伸直，髕骨挪往身體前方。
讓人流暢蹲下、固定站直。

Q-angle：大腿（髂前上棘ASIS與髕骨中心的連線）小腿（髕骨中心與脛骨粗隆的連線）的冠狀面夾角。女生骨盆較寬，Q-angle大，力矩大，十字韌帶較易斷裂。

Lachmann test：膝關節韌帶測試。用手推和拉。

stiffness圖
   force   
     |   | (nonlinear)
     |  /   arctan
----------> x
  /  |
 |   |

肩關節：上臂骨與肩胛骨的關節是球形，用來旋轉。肩胛骨與鎖骨的關節是吊帶，用來聳肩。三塊骨頭的連動機制我沒有學會。

脊椎關節：脊椎關節有如日式木桶，橫向鐵箍圈住縱向木板，中間成分主要是水。脊骨形成波浪狀，使得頭、軀幹、骨盆的重心呈一直線，垂直於地面。

金魚尾巴圖
日式木桶圖
美麗的野獸圖

muscle simulation

Hill's muscle model
https://nbviewer.org/github/demotu/BMC/blob/master/notebooks/MuscleSimulation.ipynb

肌肉-韌帶施力圖
https://nbviewer.org/github/BMClab/BMC/blob/master/notebooks/MusculoskeletaModelingSimulation.ipynb