[人体运动分析]关节中心的计算

文章目录

[人体运动分析]关节中心的计算
- 一、数据介绍
- 二、技术方法
- - 2.1 髋关节关节中心估计
  - - 2.1.1 Pelvis-LCS构建方法
    - 2.1.2 Pelvis-LCS内HJC的坐标
    - 2.1.3 HJC在GCS内的坐标
  - 2.2 膝/踝关节计算模板
  - - 2.2.1 以已知关节中心建立LCS
    - 2.2.2 Vitural Marker求解
    - 2.2.3 Vitural Marker多解的确定
- 三、完整代码

利用所给的由运动捕捉系统抓取的 Marker 点坐标数据，使用课程学习知识，通过计算虚拟跟踪标记点 ( Virtual Marker ) 点来估计关节中心。我们一定要画出右腿膝关节中心；我们可以选择性的完成左右踝关节中心的计算，并利用 Matlab 平台绘制动画。

一、数据介绍

变量名	变量含义
RASI/LASI	盆骨髂前上棘
RPSI/LPSI	盆骨髂后上棘
RTHI/LTHI	股骨跟踪标记点
RKNE/LKNE	膝关节跟踪标记点
RTIB/LTIB	胫骨跟踪标记点
RANK/LANK	踝关节跟踪标记点
RHEE/LHEE	足后跟跟踪标记点
RTOE/LTOE	脚趾跟踪标记点
InterAsisDistance	左右髂前上棘之间的距离
AsisTrocDist	髂前上棘到髋关节中心的水平距离
KneeWidth	膝关节直径
AnkleWidth	踝关节直径
Beta	骨盆倾角
Theta	髋关节中心与髂前上棘连线与竖直方向夹角
C	Davis1991算法输入参数

二、技术方法

2.1 髋关节关节中心估计

对髋关节中心的 ( Hip Joint Center，HJC ) 估计我们使用 Davis1991 方法来进行计算。首先，我们要建立盆骨坐标系 ( Pelvis-LCS ) ，在 Pelvis-LCS 下我们才可以使用 Davis1991 所描述的 HJC 估计方法来进行计算。

2.1.1 Pelvis-LCS构建方法

Pelvis-LCS 构建方法可以依据如下数学运算进行

Y ˉ P E L V I S = LASI − RASI ∣ LASI ‾ − RASI ‾ ∣ \bar{Y}_{\mathrm{PELVIS}}=\frac{\text { LASI }-\text { RASI }}{|\overline{\text { LASI }}-\overline{\text { RASI }}|} YˉPELVIS=∣ LASI − RASI ∣ LASI − RASI

Z ˉ P E L V I S = M A S I ‾ − S A C R ‾ ∣ M A S I ‾ − S A C R ‾ ∣ × Y ‾ P E L V I S \bar{Z}_{\mathrm{PELVIS}}=\frac{\overline{M A S I}-\overline{SACR}}{\mid \overline{MASI}-\overline{SACR}\mid}\times \overline{Y}_{\mathrm{PELVIS}} ZˉPELVIS=∣MASI−SACR∣MASI−SACR×YPELVIS

Y ‾ P E L V I S \overline{Y}_{\mathrm{PELVIS}} YPELVIS 和 Z ‾ P E L V I S \overline{Z}_{\mathrm{PELVIS}} ZPELVIS 分别表示 Pelvis-LCS 的侧向轴和轴向轴，其中 M A S I ‾ \overline{MASI} MASI 和 S A C R ‾ \overline{SACR} SACR 分别表示左右髂前上棘的中点和左右髂后上棘的中点。通过 Y ‾ P E L V I S \overline{Y}_{\mathrm{PELVIS}} YPELVIS 和 Z ‾ P E L V I S \overline{Z}_{\mathrm{PELVIS}} ZPELVIS 的叉积，即 X ‾ P E L V I S = Z ‾ P E L V I S × Y ‾ P E L V I S \overline{X}_{\mathrm{PELVIS}}=\overline{Z}_{\mathrm{PELVIS}}\times\overline{Y}_{\mathrm{PELVIS}} XPELVIS=ZPELVIS×YPELVIS，即得 Pelvis-LCS 的三个坐标轴， ( X ‾ P E L V I S , Y ‾ P E L V I S , Z ‾ P E L V I S ) (\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}}) (XPELVIS,YPELVIS,ZPELVIS)。

关于Rotation旋转矩阵的求解，见解释：点击这里跳转

function [O_PelvisLCS, Rotation_PelvisLCS] = establish_PelvisLCS(RASI,LASI,RPSI,LPSI)
% 函数establish_PelvisLCS：构建骨盆Pelvis-LCS坐标系
% 输出平动向量O_PelvisLCS和旋转矩阵Rotation_PelvisLCSMASI = ( RASI+LASI )/2;% 左右髂前上棘中点
SACR = ( RPSI+LPSI )/2;% 左右髂后上棘中点
PelvisLCS_O = MASI;% LCS坐标系原点
PelvisLCS_X = (MASI-SACR)/norm(MASI-SACR);
PelvisLCS_Y = (LASI-RASI)/norm(RASI-LASI);% 侧向轴使用的是左右髂前上棘的连线
PelvisLCS_Z = cross(PelvisLCS_X,PelvisLCS_Y);% 这里z轴是向下的，与教材不同
PelvisLCS_Z = PelvisLCS_Z/norm(PelvisLCS_Z);Rotation_PelvisLCS = [PelvisLCS_X,PelvisLCS_Y,PelvisLCS_Z]';% 关键-获得旋转矩阵;为什么要转置?(见链接)
O_PelvisLCS = PelvisLCS_O;
end

2.1.2 Pelvis-LCS内HJC的坐标

完成 Pelvis-LCS 建立之后，使用 Davis1991 方法，依据公式

X H = − [ x d i s + r marker ] cos ⁡ ( β ) + C cos ⁡ ( θ ) sin ⁡ ( β ) X_{H}=-\left[x_{\mathrm{dis}}+r_{\text {marker }}\right] \cos (\beta)+C \cos (\theta) \sin (\beta) XH=−[xdis+rmarker ]cos(β)+Ccos(θ)sin(β)

Y H = S [ C sin ⁡ ( θ ) − d ASIS 2 ] Y_{H}=S\left[C \sin (\theta)-\frac{d_{\text {ASIS }}}{2}\right] YH=S[Csin(θ)−2dASIS ]

Z H = − [ x dis + r marke r ] sin ⁡ ( β ) − C cos ⁡ ( θ ) cos ⁡ ( β ) Z_{H}=-\left[x_{\text {dis }}+r_{\text {marke r}}\right] \sin (\beta)-C \cos (\theta) \cos (\beta) ZH=−[xdis +rmarke r]sin(β)−Ccos(θ)cos(β)

该公式的由来可以通过如图2.1.2 中的几何关系求得。其中， d A S I S d_{ASIS} dASIS 为左右髂前上棘之间的距离， x d i s x_{dis} xdis 为髂前上棘到髋关节中心的水平距离。公式参数的含义和几何意义分别在 “数据表格” 和 “图2.1.2” 中表示出来，

此时， ( X H , Y H , Z H ) (X_H, Y_H, Z_H) (XH,YH,ZH) 即为当前 HJC 在 Pelvis-LCS 中的坐标。

2.1.3 HJC在GCS内的坐标

Pelvis-LCS 的原点为 M A S I ‾ \overline{MASI} MASI，GCS 坐标系到 Pelvis-LCS 的旋转矩阵 R P e l v i s L C S R_{PelvisLCS} RPelvisLCS 也可以由 ( X ‾ P E L V I S , Y ‾ P E L V I S , Z ‾ P E L V I S ) (\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}}) (XPELVIS,YPELVIS,ZPELVIS) 并归一化后来表示。由 2.1.2 可以知道，我们已经得到 HJC 在 Pelvis-LCS 中的坐标 ( X H , Y H , Z H ) (X_H, Y_H, Z_H) (XH,YH,ZH) ，表示为 P ′ P' P′ ，由此我们可以求得 HJC 在 GCS 坐标系内的坐标 P P P 。GCS 坐标系到 Pelvis-LCS 的平动向量

O ‾ = M A S I ‾ \overline{O}=\overline{MASI} O=MASI

则 HJC 在 GCS 内的坐标依据如下公式可以求得

P = R P e l v i s L C S ′ P ′ + O ‾ P=R'_{PelvisLCS}P' + \overline{O} P=RPelvisLCS′P′+O

function [GCS_HJC_x, GCS_HJC_y, GCS_HJC_z] = Estimation_GCS_HJC(C,strRL,theta,beta,InterAsisDistance,AsisTrocDist,...
O_PelvisLCS, Rotation_PelvisLCS)
% 函数Estimation_GCS_HJC：估计髋关节的关节中心位置
% strRL 用来指示对右/左侧腿的估计r_marker = 7;% 标记点半径
C_R = C;% C_R = 0.115*LegLength_R-0.0153=91.4097;% 判断左右腿的髋关节中心
if strcmp(strRL,'Right')S = 1;
elseif strcmp(strRL,'Left')S = -1;
elseerror('strRL is wrong');
end% 髋关节中心估计公式
% % 所获得的HJC坐标是Pelvis-LCS坐标系下的LCS坐标
LCS_HJC_Y = S*(C_R*sin(theta)-InterAsisDistance/2);
dsis = AsisTrocDist;% AsisTrocDist 髂前上棘到髋关节中心的水平距离
LCS_HJC_X = -( (dsis+r_marker)*cos(beta)-C_R*cos(theta)*sin(beta) );
LCS_HJC_Z = -C_R*cos(theta)*cos(beta)-(dsis+r_marker)*sin(beta);% 获得HJC的GCS坐标(已知LCS求GCS)
GCS = Rotation_PelvisLCS'*[LCS_HJC_X; LCS_HJC_Y; LCS_HJC_Z] + O_PelvisLCS;
GCS_HJC_x = GCS(1); GCS_HJC_y = GCS(2); GCS_HJC_z = GCS(3);end

2.2 膝/踝关节计算模板

对于膝关节和踝关节中心的计算，我们往往要构建虚拟跟踪标记点 Vitural Marker 来帮助预测膝关节和踝关节的中心。我们对 Vitural Marker 的计算有三个约束条件，如 “图2.2” 所示

(1) 该 Vitural Marker 与已知关节中心的连线，以及与当前关节外侧的解剖跟踪标记点的连线应当互相垂直；
      (2) 该 Vitural Marker 必须位于该下肢其他三个跟踪标记点(包含已知关节中心)所确定的平面当中；
      (3) 该 Vitural Marker 到当前关节外侧的解剖跟踪标记点的距离为该关节的直径。
依据上述三个约束条件，可以求得 Vitural Marker 三个位置坐标点。
      我们将通过符号计算的方法来完成 Vitural Marker 点的计算。考虑到通过符号计算的方法，若存在三元方程组，将极大的增加运算时间，因此，我们选择将三维坐标系降阶到二维坐标系。

2.2.1 以已知关节中心建立LCS

如图2所示，我们令到 Plane Definition Marker 的向量为 Y ‾ \overline{Y} Y，令 Known Joint Center 到 Joint Marker 的向量为 K ‾ \overline{K} K，通过向量叉积得到前向轴 X ‾ \overline{X} X，即 X ‾ = Y ‾ × K ‾ \overline{X} = \overline{Y} \times \overline{K} X=Y×K，轴向轴 Z ‾ = X ‾ × Y ‾ \overline{Z} = \overline{X} \times \overline{Y} Z=X×Y。
如2.1.3节，依据 Known Joint Center 可以得到 GCS 到当前所构建 LCS 的平动向量，依据归一化 ( X ‾ P E L V I S , Y ‾ P E L V I S , Z ‾ P E L V I S ) (\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}}) (XPELVIS,YPELVIS,ZPELVIS) 可以得到 GCS 到当前所构建 LCS 的旋转矩阵R。

2.2.2 Vitural Marker求解

根据 2.2 节中约束条件(2)可以知道，Vitural Marker 位于图2所示的蓝色平面中。该蓝色平面由 ( Y ‾ P E L V I S , Z ‾ P E L V I S ) (\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}}) (YPELVIS,ZPELVIS) 确定，也即 Vitural Marker 的 x 轴坐标为 0。此时，我们完成了对原始方程组的降阶。

2.2.3 Vitural Marker多解的确定

在 Matlab 平台内我使用 vpasolve 函数求解我们的方程组。然而 vpasolve 函数无法给出除多项式输入之外的其他方程组的全部解，而 vpasolve 函数的这一局限性，要求我们需要设置解的初始搜寻点。
依据 vpasolve 函数的方法文档，我们引入一个初始搜寻点 ( Initial Point ) 。分析 “图2.2”，我们可以知道，所需要的 Vitural Marker 解靠近 2.2.1 节内所确定 LCS 的 Z ‾ \overline{Z} Z 轴，因此该 Initial Point 确定为在 Z ‾ \overline{Z} Z 轴上的点，即 Y ‾ \overline{Y} Y 轴坐标为 0 。我们设 Initial Point 的 Z ‾ \overline{Z} Z 轴与 Joint Marker 的 Z ‾ \overline{Z} Z 轴坐标相同。

% 符号计算解法
function VirtualMarker = Estimation_GCS_VMarker(KnonwnJointCenter,PlaneDefinitionMarker,JointMarker,Width)
% 函数Estimation_GCS_VMarker：用来求内侧虚拟点VirtualMarker
% KnonwnJointCenter:HJC
% JointMarker:KNE
% PlaneDefinitionMarker:THI
% Width:KneeWidth% % 围绕已知关节中心建立坐标系Known_LCS
Known_LCS_i = PlaneDefinitionMarker-KnonwnJointCenter;
Known_LCS_p = JointMarker-KnonwnJointCenter;
Known_LCS_k = cross(Known_LCS_i,Known_LCS_p);
Known_LCS_j = cross(Known_LCS_k,Known_LCS_i);
Known_LCS_i = Known_LCS_i/norm(Known_LCS_i);
Known_LCS_j = Known_LCS_j/norm(Known_LCS_j);
Known_LCS_k = Known_LCS_k/norm(Known_LCS_k);
% % 求JointMarker在LCS内的坐标
Known_LCS_O = KnonwnJointCenter;% 平动向量
RotationKnownLCS = [Known_LCS_i, Known_LCS_j, Known_LCS_k]';% RotationKnownLCS为GCS到Known_LCS的旋转矩阵
JointMarker_knownlcs = RotationKnownLCS*(JointMarker-Known_LCS_O);% 约束条件1，位于平面内--面法向量相同
syms x y;% 设置待求点的坐标(由于是在xOy平面，所以没有设z值)
p_vk = JointMarker_knownlcs(1:2)-[x,y]';% 构建向量
p_vh = -[x,y]';
% 约束条件2，JointMarker和KnonwnJointCenter与该点的连线是互相垂直
eqn1 = dot(p_vk,p_vh);% 点积为0表示垂直
% 约束条件3，虚拟点与KNE之间的距离为KNEE关节的直径
eqn2 = norm(p_vk)-Width;eqn = [eqn1,eqn2];
var = [x,y];
Pinitial = [0,JointMarker_knownlcs(2)];% 设置初始值
[VM_x,VM_y] = vpasolve(eqn,var,Pinitial);% vpasolve不能返回所有的解
VirtualMarker = [VM_x,VM_y,0]';%ok
VirtualMarker = RotationKnownLCS'*VirtualMarker + Known_LCS_O;% LCS到GCS %wrong
end

三、完整代码

所需的“LowerLimbMarkers.mat”将上传到资源当中，点击这里跳转。

您也可以自己编程，本篇编程的难点在第二部分关于虚拟标记点的计算上(可能会提示“无解析解”)，以及旋转矩阵的构建技巧。
提示“无解析解”一般有三种原因：1). 代码写错；2). 初值点选择有误(就本文代码来说)；3). 的确没有解析解，需要通过数值计算方法来找到数值解。

本篇所附代码的不足在于，没有提供数值解法。

% % 2021.01.15 Jlin.Zheng (Sun Yat-sen University)
clc;clear
load('LowerLimbMarkers.mat');% 每一帧的坐标数据都是横向排列的，为了算法需要，在函数的参数值中选择将其转置为列向量%% 定义GCS坐标系
% 由于缺少旋转矩阵，无法逆推GCS坐标%% 计算膝关节中心
% 先求解RHJC坐标位置 Davis1991
[O_PelvisLCS, Rotation_PelvisLCS] = establish_PelvisLCS(RASI(1,:)',LASI(1,:)',RPSI(1,:)',LPSI(1,:)');
[GCS_HJC_x_Right, GCS_HJC_y_Right, GCS_HJC_z_Right] = Estimation_GCS_HJC(C,'Right',theta,beta,InterAsisDistance,AsisTrocDist_R,...
O_PelvisLCS, Rotation_PelvisLCS);
[GCS_HJC_x_Left, GCS_HJC_y_Left, GCS_HJC_z_Left] = Estimation_GCS_HJC(C,'Left',theta,beta,InterAsisDistance,AsisTrocDist_L,...
O_PelvisLCS, Rotation_PelvisLCS);
% 获得股骨内测虚拟点virtualmarker
GCS_RightHJC = [GCS_HJC_x_Right; GCS_HJC_y_Right; GCS_HJC_z_Right];
VirtualMarker_Right_Knee = Estimation_GCS_VMarker(GCS_RightHJC,RTHI(1,:)',RKNE(1,:)',RightKneeWidth);
VirtualMarker_Right_Knee = double(VirtualMarker_Right_Knee);
GCS_LeftHJC = [GCS_HJC_x_Left; GCS_HJC_y_Left; GCS_HJC_z_Left];
VirtualMarker_Left_Knee = Estimation_GCS_VMarker(GCS_LeftHJC,LTHI(1,:)',LKNE(1,:)',LeftKneeWidth');
VirtualMarker_Left_Knee = double(VirtualMarker_Left_Knee);
% 获得膝关节中心KJC
GCS_RightKJC = ( VirtualMarker_Right_Knee+RKNE(1,:)' )/2;
GCS_LeftKJC = ( VirtualMarker_Left_Knee+LKNE(1,:)' )/2;%% 计算踝关节中心
VirtualMarker_Right_Ankle = Estimation_GCS_VMarker(GCS_RightKJC,RTIB(1,:)',RANK(1,:)',RightAnkleWidth);
VirtualMarker_Right_Ankle = double(VirtualMarker_Right_Ankle);
VirtualMarker_Left_Ankle = Estimation_GCS_VMarker(GCS_LeftKJC,LTIB(1,:)',LANK(1,:)',LeftAnkleWidth);
VirtualMarker_Left_Ankle = double(VirtualMarker_Left_Ankle);
GCS_RightAJC = ( VirtualMarker_Right_Ankle+RANK(1,:)' )/2;
GCS_LeftAJC = ( VirtualMarker_Left_Ankle+LANK(1,:)' )/2;%% 计算全部帧左右膝关节中心和踝关节中心并制作动画
hp = figure;
hold on;
%% 描点
% 跟踪标记点绘制
TracTrocMarker = [RASI(1,:)',LASI(1,:)',RPSI(1,:)',LPSI(1,:)',RKNE(1,:)',LKNE(1,:)',RTHI(1,:)',LTHI(1,:)'...RTIB(1,:)',LTIB(1,:)',RANK(1,:)',LANK(1,:)',...RTOE(1,:)',LTOE(1,:)',RHEE(1,:)',LHEE(1,:)'];
scatter3(TracTrocMarker(1,:)',TracTrocMarker(2,:)',TracTrocMarker(3,:)',...'MarkerEdgeColor','black','MarkerFaceColor',[0.00,0.45,0.74]);
% 虚拟跟踪标记点绘制
TraceVirtualMarker = [GCS_RightHJC,GCS_LeftHJC,GCS_LeftKJC,GCS_RightKJC,GCS_RightAJC,GCS_LeftAJC];
scatter3(TraceVirtualMarker(1,:)',TraceVirtualMarker(2,:)',TraceVirtualMarker(3,:)',...'MarkerEdgeColor','black','MarkerFaceColor',[0.85,0.33,0.10]);
%% 连线
% 骨盆连线
PelvisMarker = [RASI(1,:)',LASI(1,:)',LPSI(1,:)',RPSI(1,:)',RASI(1,:)'];
plot3(PelvisMarker(1,:)',PelvisMarker(2,:)',PelvisMarker(3,:)','Color','black','LineWidth',1.5);
% 骨盆补充连线2
RightPelvisHJC = [RASI(1,:)',GCS_RightHJC,RPSI(1,:)'];
LeftPelvisHJC = [LASI(1,:)',GCS_LeftHJC,LPSI(1,:)'];
plot3(RightPelvisHJC(1,:)',RightPelvisHJC(2,:)',RightPelvisHJC(3,:)','Color','black','LineWidth',1.5);
plot3(LeftPelvisHJC(1,:)',LeftPelvisHJC(2,:)',LeftPelvisHJC(3,:)','Color','black','LineWidth',1.5);
% 下肢连线
RightLegMarker = [GCS_RightHJC, GCS_RightKJC, GCS_RightAJC];
LeftLegMarker = [GCS_LeftHJC, GCS_LeftKJC, GCS_LeftAJC];
plot3(RightLegMarker(1,:)',RightLegMarker(2,:)',RightLegMarker(3,:)','Color','black','LineWidth',1.5);
plot3(LeftLegMarker(1,:)',LeftLegMarker(2,:)',LeftLegMarker(3,:)','Color','black','LineWidth',1.5);
%% 设置图窗大小
xlim([-750,750]);ylim([-750,750]);zlim([0,1500]);
set(hp,'Position',[20,20,400,600]);
view([210,20])%% 所用代码
% 记得把前面提及的函数代入噢~