motion - 演算法筆記

human motion

plagiarism declaration

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

human motion

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

human measurement（anthropometry）

人體測量。觀察人體外表。例如身高、體重、三圍、足長。

human body segments 人體肢段

head      頭　  ┬ head and neck 頭頸　
neck      頸　  ┘
torso     胸　  ┐ 
abdomen   腹　  ┼ trunk         軀幹
pelvis    骨盆  ┘
upper arm 上臂  ┐
forearm   前臂  ┼ upper limb    上肢
hand      手　  ┘
thigh     大腿  ┐
shank     小腿  ┼ lower limb    下肢
foot      足　  ┘

human body joints 人體關節

shoulder  肩    ┐
elbow     肘    ┼ upper limb 上肢
wrist     腕    ┘
hip       髖    ┐
knee      膝    ┼ lower limb 下肢
ankle     踝    ┘

human dissection（anatomy）

人體解剖。觀察人體內裡。例如骨骼、關節、肌肉、神經。

解剖方位。方便描述人體所有東西的位置。

anatomical position
解剖學姿勢。全身直立站好，手心朝向正面。
（人體躺在解剖台的自然姿勢。）

anatomical planes
解剖學平面。空間三個維度對應三種平面。

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

anatomical directions
解剖學方向。空間三個維度對應六種方向。

superior  上。朝頭。
inferior  下。朝腳。
anterior  前。正面。
posterior 後。背面。
medial    內。中央。
lateral   外。左右。
proximal  近端。靠近軀幹。
distal    遠端。遠離軀幹。

human movement system

人體動作系統。活動的器官的集合。區分許多系統。

心血管系統　　　　cardiovascular system
肺系統　　　　　　pulmonary system
內分泌系統　　　　endocrine system
皮膚系統　　　　　integumentary
神經系統　　　　　nervous system
肌肉骨骼系統　　　musculoskeletal system

本文只涉及肌肉骨骼系統。大家視討論對象重組詞彙。

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

肌肉骨骼系統的構造。

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

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

human motor system

人體運動系統。神經控制暨肌肉運動的機制。

上肢運動系統　upper-limb motor system
下肢運動系統　lower-limb motor system
臉動系統　　　facial motor system
眼動系統　　　oculomotor system

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

開源軟體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]適合教學。

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

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

四、螢光攝影fluorography／放射攝影radiography：
　　利用顯影劑與X光機估計人體內部器官形狀。
　　原理是輻射吸收radiation absorption。
　　領域是放射學radiology。

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

一、測力板force plate：測量力。
　口、足壓儀pedobarograph：
　　　測量壓力。二維陣列，一格一個數值。
　　　得到足壓圖pedobarogram。
　口、三維測力板3D force plate：
　　　測量力、力矩。總共六個數值。主流作法。
　　　得到地面反作用力ground reaction force。
二、動力計dynamometer：測量力矩。
　口、等張肌力儀isotonic dynamometer。
　口、等長肌力儀isometric dynamometer。
　口、等速肌力儀isokinetic dynamometer。
　　　得到肌群力矩muscle group torque。

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

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

順便整理一下各種英文字尾的意義。

-gram   圖表
-graph  製圖的設備
-graphy 製圖的行為
-scope  檢視的設備
-scopy  檢視的行為
-meter  測量的設備
-metry  測量的行為

生物學和醫學當中，-graphy和-scopy經常混用。網路上有一種說法是-graphy又觀察又紀錄，-scopy只觀察不紀錄。網路上還有一種說法是-graphy離散時刻，-scopy連續時間。

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.qualisys.com/accessories/calibration-kits/
https://www.qualisys.com/accessories/calibration-kits/carbon-fibre-calibration-kit/

inertial motion capture

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

IMU的特色是重量輕、體積小。貼在身體上、縫在衣服上，不影響人體運動。IMU的特色是地點不受限，適合用於室外活動。IMU內含無線傳輸裝置，藉此即時獲取感測器數值。

IMU內含三種感測器：加速規（平移加速度）、陀螺儀（旋轉速度）、磁力計（北極方向）。利用三種感測器的數值，得以估計自身位置和方向。演算法Mahony filter或Madgwick filter。

IMU經常產生誤差。經典問題gyroscope drift。解法sensor fusion。演算法 complementary filter。

加速規與陀螺儀，分別利用球座標系與離散積分，計算角位置。然後計算加權平均數，融合成一個角位置。

acc_roll = atan2(acc.y, acc.z);
acc_pitch = atan2(-acc.x, hypot(acc.y, acc.z));
acc_yaw = undefined;
acc_angle = [acc_roll, acc_pitch, acc_yaw];
angle = w * acc_angle + (1-w) * (angle + gyro_vel * Δt);

加速規靜止不動，量到重力加速度，藉此計算角位置。加速規移動越快，角位置越不精準。陀螺儀恰好相對、恰好互補。一種解法是加速規套用低通濾波器（角位置幾乎不變，接近靜止），陀螺儀套用高通濾波器（角位置大幅改變），然後兩者疊加。

時域線性濾波器的情況下，兩者疊加恰是加權平均數、恰是走半步。方便起見，大家都是直接走一步。

low-pass filter
y[n] = w * x[n] + (1-w) * y[n-1];

low-pass filter for accelerometer (a) 
angle = w * acc_angle + (1-w) * angle;

high-pass filter
y[n] = (1-w) * (x[n] - x[n-1]) + (1-w) * y[n-1];

high-pass filter for gyroscope (b)
angle = (1-w) * (gyro_vel * Δt) + (1-w) * angle;

complementary filter ½(a) + ½(b)
angle = w * ½acc_angle + (1-w) * (angle + gyro_vel * ½Δt);

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

markerless motion capture

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

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

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

fluoroscope / radioscope

螢光鏡／放射鏡。使用顯影劑與X光機，估計人體內部器官形狀，得到一段影片／一張圖片。

首先服用或注射染色劑（顯影劑），等待染色劑抵達預定拍攝的人體器官。然後照射輻射（X光），讓染色劑吸收輻射能量（進而產生螢光），導致輻射難以傳遞（穿透）。輻射難以傳遞導致底片不變，輻射容易傳遞導致底片曬黑，以此區分器官形狀。

動作分析主要是拍攝骨頭。骨頭孔隙不勻，輻射難以傳遞。血肉綿密均勻，輻射容易傳遞。

人體照射輻射有副作用，攝影行程需要嚴格管理。因此只有醫療設施才會對人體使用螢光鏡／放射鏡。

順帶一提，X光檢查的螢光鏡（一般稱做螢光透視攝影）、切片檢查的螢光鏡，它們是兩回事。切片檢查當中，切片很薄、顯而易見，得以直接測量螢光，而不是去測量X光。

pedobarograph

足壓儀。測量矩形區域每一處的壓力，得到一個矩陣，稱作足壓圖pedobarogram。

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

測量原理主要分為兩種。

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

3D force plate

三維測力板。測量力與力矩，得到兩個三維向量，稱作地面反作用力ground reaction force。

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

進一步計算合力中心（但是大家仍然沿用壓力中心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₝可能是零，所以不能當作分母。

多臺測力板的壓力中心的計算公式如下。

一、測力板：座標系、施力點。
　　讀取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₞
　　正確名稱應是合力中心。阿就歷史共業。

計算原理是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

等速肌力儀。測量一個關節轉動所產生的力矩，稱作肌群力矩muscle group torque。

肌肉收縮導致關節轉動。一個關節轉動源自一群肌肉收縮。復健醫學當中，肌肉收縮分成兩種：縮短與伸長。伸長也是一種收縮。

concentric contraction 向心收縮　　縮短
eccentric contraction  離心收縮　　伸長

肌肉活化量、肌肉長度、肌肉速度，決定肌肉力量。肌肉收縮機制隨之分成三類。測量儀器隨之分成三類。

isotonic contraction   等張收縮　　肌肉張力不變／阻力不變
isometric contraction  等長收縮　　肌肉長度不變
isokinetic contraction 等速收縮　　肌肉速度不變

isotonic dynamometer   等張肌力儀　固定肢段負重。例如舉放啞鈴。
isometric dynamometer  等長肌力儀　固定關節角度。例如提啞鈴。
isokinetic dynamometer 等速肌力儀　固定關節角速度。以機械設備壓抑速度。

electromyograph

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

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

我們需要訊號功率。四個步驟。細節請見專著《Biomechanics and Motor Control of Human Movement》。

包絡線envelope：
一、整流rectification：取絕對值。
二、平滑smoothing：低通濾波器擷取20-200Hz的訊號。沒有固定說法。
功率power：
三、積分integration：求平方和，得到能量。再除以訊號長度，得到功率。
四、放大amplification：估計阻抗大小，還原原始功率。

商業產品圖片，請自行點擊網址觀看。
https://books.google.com.tw/books?id=_bFHL08IWfwC&pg=PA270

測量手段分為兩種類型。細節請見專著《The ABC of EMG》。

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

human motion analysis

stabilogram analysis

穩定圖分析。三維測力板上站著不動，所有時刻的壓力中心center of pressure收集起來，形成軌跡，稱作穩定圖stabilogram。分析穩定圖，藉此知道阿嬤腳麻的情況。已有多種分析手法。

principal component analysis：壓力中心軌跡的平均數與主成分，建立橢圓範圍，作為晃動範圍。順帶一提，根據統計學信賴區間，真實平均數落在平均數正負兩個標準差構成的橢圓範圍之內的機率是95%。

rambling and trembling decomposition：壓力中心軌跡分離成晃動軌跡rambling trajectory與抖動軌跡trembling trajectory兩者相加。

一、每當壓力任一維度變號（函數曲線穿過時間軸），紀錄此時壓力中心位置。
二、晃動軌跡：上述過零位置做三次樣條內插。
三、抖動軌跡：壓力中心軌跡減去晃動軌跡。

stabilogram diffusion analysis：固定時間間隔，計算兩個壓力中心位置差距長度平方。枚舉所有時刻，計算平均數。各種時間間隔的平均數，繪製成函數曲線。模仿Hurst exponent，在0秒到2.5秒的區段，找到二次導數最小值位置，切成前後兩段函數曲線，分別找到一次迴歸函數的斜率。

gait analysis

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

下圖是步行的一個週期的每個階段名稱。

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

下圖是每個時刻的地面反作用力的XYZ分量，左右前後上下。

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

下圖是每個時刻的下肢肌肉的肌肉激發量。

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

下圖是每個時刻的下肢關節的旋轉角度。

一、肢段：座標系統。
　口、肢段座標系統。
　　　根據標準規範，取得對應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，定義了膝關節歐拉角三個轉軸。

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

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

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

joint force estimation

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

專著《Kinematics of Human Motion》。

主流作法是inverse kinematics。數學模型是wobbling mass model。

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

一、肢段：端點。《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。

最佳化範例請見專著，或者下方課程影片。

另一種作法是利用electromyogram計算而得，數學模型是muscle activation dynamics，請見後面運動控制章節。

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

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

coordination analysis

協調分析。關節們、肌肉們的連動規律。

下圖是每個時刻的兩個關節角度的數對，稱作cyclogram。

human motion control🚧

human motion control

gait

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

大腦感知肌肉狀態、估計肢段位置，稱作本體感覺。閉起眼睛讓兩根食指相碰，仰賴本體感覺。運動、演奏、盲打，仰賴本體感覺。

the man who lost his body是一部紀錄片，描述一位喪失本體感覺的患者如何面對疾病。19歲忽然發病，脖子以下無法如願移動。躺在床上十幾年。經過數年逐漸學會控制身體各個部位。又經過數年學會起身、下床、站起、走路。由於無法感知肌肉狀態，光是拿起蘋果就會導致重心不穩跌倒。做任何動作都必須隨時利用眼睛觀察姿勢，確保肌肉狀態正確。發病之後三十年，終於回到日常生活。

gait simulation

走路模型。

ballistic walking model
passive dynamic walking model

也許介紹一下外骨骼。

joint

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

1. gliding      滑動，關節橫向平移
2. distraction  分散，關節拉開
3. compression  壓縮，關節靠近

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

1. flexion-extension  上肢下肢彎曲伸直
2. abduction-adduction 下肢左右搖擺
3. internal-external rotation 下肢自旋
https://link.springer.com/chapter/10.1007/978-3-031-25322-5_11/figures/23
https://josr-online.biomedcentral.com/articles/10.1186/s13018-020-02047-0/figures/2

關節轉動角度範圍各有不同。

range of motion evaluation
https://www.dshs.wa.gov/sites/default/files/forms/pdf/13-585a.pdf
https://musculoskeletalkey.com/restoring-range-of-motion-and-improving-flexibility/
https://www.amazon.com/Range-Motion-Testing-Set-Charts/dp/0965342522
https://www.semanticscholar.org/paper/7e602a134667ff5f6f5dda8c7608a59f204d662c
https://www.realbodywork.com/joint-range-of-motion/

關節最多有六個自由度degree of freedom。三個平移維度、三個旋轉維度。

degree of freedom of human joint table
https://www.researchgate.net/figure/The-degrees-of-freedom-DoF-and-the-anatomical-segments-of-the-Rugby-Model-Body_fig1_312094367

關節分為六種，根據自由度。

1. ball-and-socket joint 例如肩關節、髖關節
2. hinge joint           例如肘關節、膝關節
3. gliding joint         例如腕關節、踝關節
4. pivot joint           例如第一第二頸椎關節（寰樞關節）
5. saddle joint          例如拇指第一關節（腕掌關節）
6. condyloid joint       例如手指第一關節、腳指第一關節

https://www.pinterest.com/pin/397090892129861559/

關節鬆弛性測試，觀察自由度。

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

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

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

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

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

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

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

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

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

脊椎關節：脊椎中間有管道讓脊髓穿過。椎間盤（軟骨）有如日式木桶，橫向鐵箍圈住縱向木板，中間成分主要是水。脊骨形成波浪狀，使得頭頸、胸腹、骨盆的重心呈一直線，垂直於地面。

金魚尾巴圖
日式木桶圖
美麗的野獸圖
https://bonebro.com/spine-structuretest/

joint simulation

人工關節製作、人工關節置換手術。

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

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

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

muscle

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

神經肌肉構造圖
https://www.youtube.com/watch?v=AMdEjf6IPPk
https://www.youtube.com/watch?v=pxzCSyK7b7g

muscle activation dynamics正向反向圖
https://jneuroengrehab.biomedcentral.com/articles/10.1186/s12984-023-01253-1
https://pmc.ncbi.nlm.nih.gov/articles/PMC5758546/

激發-活化圖
excitation–activation

肌肉纖維長度-力量圖
sarcomere length-force diagram
https://assets.coursehero.com/study-guides/lumen/images/aacc-ap1/2489/Muscle_1611.jpg

長度-速度-肌力圖
https://www.researchgate.net/figure/Normalized-muscle-force-F-as-a-function-of-length-lm-and-velocity-v-A-The-net_fig5_270514950
https://www.researchgate.net/figure/Force-length-and-velocity-relationships-in-a-Hill-type-muscle-model-Left-Normalised_fig2_320894213

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

不同姿勢需要不同肌肉力量。例如握住啞鈴不動，上臂保持水平，前臂向上弓90°、向下弓90°，需要的肌肉力量不同。從結構力學的觀點，兩種姿勢的力矩完全相同。但是從生理結構的觀點，兩種姿勢的肌肉長度略有不同。

muscle simulation

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