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

最新推荐文章于 2021-12-21 09:50:10 发布

iTsing203

最新推荐文章于 2021-12-21 09:50:10 发布

阅读量2.2k

点赞数 1

分类专栏：人体运动分析[笔记] 文章标签：人体运动分析关节中心计算经验分享

本文链接：https://blog.csdn.net/qq_44207723/article/details/112064823

版权

人体运动分析[笔记] 专栏收录该内容

7 篇文章 7 订阅

订阅专栏

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

文章目录

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

利用所给的由运动捕捉系统抓取的 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 构建方法可以依据如下数学运算进行

$\bar{Y}_{\mathrm{PELVIS}}=\frac{\text { LASI }-\text { RASI }}{|\overline{\text { LASI }}-\overline{\text { RASI }}|}$

$\bar{Z}_{\mathrm{PELVIS}}=\frac{\overline{M A S I}-\overline{SACR}}{\mid \overline{MASI}-\overline{SACR}\mid}\times \overline{Y}_{\mathrm{PELVIS}}$

$\overline{Y}_{\mathrm{PELVIS}}$ 和 $\overline{Z}_{\mathrm{PELVIS}}$ 分别表示 Pelvis-LCS 的侧向轴和轴向轴，其中 $\overline{MASI}$ 和 $\overline{SACR}$ 分别表示左右髂前上棘的中点和左右髂后上棘的中点。通过 $\overline{Y}_{\mathrm{PELVIS}}$ 和 $\overline{Z}_{\mathrm{PELVIS}}$ 的叉积，即 $\overline{X}_{\mathrm{PELVIS}}=\overline{Z}_{\mathrm{PELVIS}}\times\overline{Y}_{\mathrm{PELVIS}}$ ，即得 Pelvis-LCS 的三个坐标轴， $(\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}})$ 。

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

function [O_PelvisLCS, Rotation_PelvisLCS] = establish_PelvisLCS(RASI,LASI,RPSI,LPSI)
% 函数establish_PelvisLCS：构建骨盆Pelvis-LCS坐标系
% 输出平动向量O_PelvisLCS和旋转矩阵Rotation_PelvisLCS

MASI = ( 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}=-\left[x_{\mathrm{dis}}+r_{\text {marker }}\right] \cos (\beta)+C \cos (\theta) \sin (\beta)$

$Y_{H}=S\left[C \sin (\theta)-\frac{d_{\text {ASIS }}}{2}\right]$

$Z_{H}=-\left[x_{\text {dis }}+r_{\text {marke r}}\right] \sin (\beta)-C \cos (\theta) \cos (\beta)$

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

此时， $X_H, Y_H, Z_H)$ 即为当前 HJC 在 Pelvis-LCS 中的坐标。

2.1.3 HJC在GCS内的坐标

Pelvis-LCS 的原点为 $\overline{MASI}$ ，GCS 坐标系到 Pelvis-LCS 的旋转矩阵 $R_{PelvisLCS}$ 也可以由 $(\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}})$ 并归一化后来表示。由 2.1.2 可以知道，我们已经得到 HJC 在 Pelvis-LCS 中的坐标 $X_H, Y_H, Z_H)$ ，表示为 $P^{'}$ ，由此我们可以求得 HJC 在 GCS 坐标系内的坐标 $P$ 。GCS 坐标系到 Pelvis-LCS 的平动向量

$\overline{O}=\overline{MASI}$

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

$P=R'_{PelvisLCS}P' + \overline{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;
else
    error('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 的向量为 $\overline{Y}$ ，令 Known Joint Center 到 Joint Marker 的向量为 $\overline{K}$ ，通过向量叉积得到前向轴 $\overline{X}$ ，即 $\overline{X} = \overline{Y} \times \overline{K}$ ，轴向轴 $\overline{Z} = \overline{X} \times \overline{Y}$ 。
如2.1.3节，依据 Known Joint Center 可以得到 GCS 到当前所构建 LCS 的平动向量，依据归一化 $(\overline{X}_{\mathrm{PELVIS}},\overline{Y}_{\mathrm{PELVIS}},\overline{Z}_{\mathrm{PELVIS}})$ 可以得到 GCS 到当前所构建 LCS 的旋转矩阵R。

2.2.2 Vitural Marker求解

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

2.2.3 Vitural Marker多解的确定

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

%% 所用代码
% 记得把前面提及的函数代入噢~