forked from modenaxe/3d-muscles
-
Notifications
You must be signed in to change notification settings - Fork 0
/
b2_plot_momArms_Fig5_as_Blemker2005.m
211 lines (188 loc) · 8.37 KB
/
b2_plot_momArms_Fig5_as_Blemker2005.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
%-------------------------------------------------------------------------%
% Copyright (c) 2020 Modenese L., Kohout J. %
% %
% Licensed under the Apache License, Version 2.0 (the "License"); %
% you may not use this file except in compliance with the License. %
% You may obtain a copy of the License at %
% http://www.apache.org/licenses/LICENSE-2.0. %
% %
% Unless required by applicable law or agreed to in writing, software %
% distributed under the License is distributed on an "AS IS" BASIS, %
% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or %
% implied. See the License for the specific language governing %
% permissions and limitations under the License. %
% %
% Author: Luca Modenese %
% email: [email protected] %
% ----------------------------------------------------------------------- %
% This script plots the moment arms of Glut Max and Glut Med from the highly
% discretised musculoskeletal models in a fashion similar to Blemker and
% Delp, Ann of Biom Eng (2005).
% Data from previous studies is also plotted.
% The input data are generated by the MATLAB script
% a_compute_biomech_moving_viapoints.m
% NB: all coordinates are assumed in degrees
% ----------------------------------------------------------------------- %
clear; clc; close all
% OpenSim libraries
import org.opensim.modeling.*
addpath('./_support_functions');
%-------------- SETTINGS -------------------------
% model
model = 'LHDL';
method = 'L_Method';
% task to use in creating the models
task_set = {'hip_flexion','hip_adduction', 'hip_rotation'};
indip_coord_set = {'hip_flexion_r','hip_adduction_r', 'hip_rotation_r'};
% folder where the biomechanical results are stored
biomech_res_folder = ['c_fibre_biomechanics/', method];
% folder where the MSK model results are stored
MSK_folder = './a_MSK_model';
% folder where the figures will be saved
temp_figure_folder = 'd_paper_figures/temp';
% external data folder (from previous publications)
literature_data_folder = '_literature_data';
%-------------------------------------------------
% name of muscles (depend on model)
muscles_list = {'r_gluteus_maximus', 'r_gluteus_medius'};
title_list = {'Gluteus Maximus', 'Gluteus Medius'};
% in_deg = 'yes';
% check folder existance
checkFolder(temp_figure_folder)
n_plot = 1;
h = figure; set(h, 'Position', [410 89 536 768]);
for n_task = 1:3
% setting variables for computations
cur_task = task_set{n_task};
sim_name = [model,'_',cur_task];
coord_name = indip_coord_set{n_task};
kinematicsMotFile = [MSK_folder, filesep, 'b_kinematics/',cur_task,'.mot'];
MSK_results_file = [MSK_folder, filesep, 'c_MSK_results/',model, '_', cur_task,'.mot'];
switch cur_task
case 'hip_flexion'
xlabel_str = 'hip flexion angle [deg]';
ylabel_str = {'hip extension'; 'moment arm [cm]'};
x_lim = [-10.0 60.0];
y_lim = [-3.0 10.0];
coeff = 1.0;% SURE
coeff_os = -1.0;
case 'hip_adduction'
xlabel_str = 'hip adduction angle [deg]';
ylabel_str = {'hip adduction'; 'moment arm [cm]'};
x_lim = [-40.0 40.0];
y_lim = [-8.0 8.5];
coeff = -1; % SURE!
coeff_os = 1.0;
case 'hip_rotation'
xlabel_str = 'hip internal rotation angle [deg]';
ylabel_str = {'hip internal rotation'; 'moment arm [cm]'};
x_lim = [-30.0 30.0];
y_lim = [-4.25 3.5];
coeff = -1; %SURE
coeff_os = 1.0;
end
% load results for simulated task
load(fullfile(biomech_res_folder, ['Results_',sim_name]));
% load mom arms from OpenSim model
osim_MA = sto2MatStruct(MSK_results_file);
% angles
q = getValueColumnForHeader(osim_MA, coord_name);
% glut max
os_glutmax(:,1) = getValueColumnForHeader(osim_MA, 'glut_max1_r');
os_glutmax(:,2) = getValueColumnForHeader(osim_MA, 'glut_max2_r');
os_glutmax(:,3) = getValueColumnForHeader(osim_MA, 'glut_max3_r');
% glut med
os_glutmed(:,1) = getValueColumnForHeader(osim_MA, 'glut_med1_r');
os_glutmed(:,2) = getValueColumnForHeader(osim_MA, 'glut_med2_r');
os_glutmed(:,3) = getValueColumnForHeader(osim_MA, 'glut_med3_r');
% creating structure with kinematics
kinMat = sto2MatStruct(kinematicsMotFile);
% kinematics
angles = getValueColumnForHeader(kinMat, coord_name);
N_mus = numel(muscles_list);
for n_m = 1:N_mus
title_str = title_list{n_m};
cur_mus = muscles_list{n_m};
% create title
tit1 = strrep(cur_mus(3:end),'_',' ');
cur_tit = [upper(tit1(1)),tit1(2:end)];
disp(['Processing ', cur_mus]);
cur_mus_res = ResultsSummary.(cur_mus);
% get Blemker color for same muscle
switch cur_mus
case 'r_gluteus_medius'
col = [0.8 0 0];
title_str = 'Gluteus Medius';
os_ = os_glutmed;
case 'r_gluteus_maximus'
col = [0 0.8 0.8];
title_str = 'Gluteus Maximus';
os_ = os_glutmax;
end
% plot
subplot(3,N_mus,n_plot);
plot( cur_mus_res.ind_coord(1:end-1)*180/pi, coeff*(cur_mus_res.mom_arm_mat*100),'Color',col); hold on
plot(q, coeff_os*os_*100,'k', 'Linewidth', 1.0)
if n_task == 1
title(title_str);
end
if mod(n_plot,2)~=0
ylabel(ylabel_str);
end
xlabel(xlabel_str);
xlim(x_lim);
ylim(y_lim)
box off
n_plot = n_plot+1;
end
end
% additional data
load([literature_data_folder, filesep,'Nemeth1985.mat']);
load([literature_data_folder, filesep,'Dostal1986.mat']);
%----------- GLUTEUS MAXIMUS ------------
% add line from Nemeth and Ohlsen 1985 (flex)
ax = subplot(3,2,1);
% line was digitised
plot(ax, Nemeth1985.data(:,1), Nemeth1985.data(:,2)/10, '-k');
% points were from paper: added first point
x = [Nemeth1985.data(1,1), 5:5:90];
y = [Nemeth1985.data(1,2), 75, 74, 72, 70, 69, 66, 64, 62, 59, 56, 54, 51, 48, 44, 41, 38, 34, 31]/10;
plot(ax, x, y, 'sk', 'MarkerFaceColor',[1 1 1],'MarkerSize',6)
% add Dostal data (flexion)
subplot(3,2,1)
plot(0.0, 4.6, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
% add Dostal data (adduction)
subplot(3,2,3)
plot(0.0, 0.7, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
% add Dostal data (int/ext rot)
subplot(3,2,5)
plot(0.0, -2.1, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
%----------- GLUT MED -----------------
% flex-ext
ax = subplot(3,2,2);
plot(ax, Dostal1986.data(:,1), -Dostal1986.data(:,2), '-k', 'Linewidth', 1.0)
plot(ax, Dostal1986.data(1:5:end,1), -Dostal1986.data(1:5:end,2),'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3);
% from Table1 of Dostal (add/abd)
subplot(3,2,4)
plot(0.0, -6.7, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
plot(0.0, -6.0, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
plot(0.0, -4.3, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
%int/ext rot
subplot(3,2,6)
plot(0.0, 2.3, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
plot(0.0, 0.1, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
plot(0.0, -2.4, 'ok', 'Linewidth', 1.0,'MarkerFaceColor',[0 0 0],'MarkerSize',3)
% from Nemeth and Ohlsen 1989
subplot(3,2,3) % glut max adductor mom arm
plot(0.0, 1.6, 'sk', 'Linewidth', 1.0,'MarkerFaceColor',[1 1 1],'MarkerSize',6)
% glut medius abductor
subplot(3,2,4) % glut med abductor mom arm
plot(0.0, -5.4, 'sk', 'Linewidth', 1.0,'MarkerFaceColor',[1 1 1],'MarkerSize',6)
%----------------------------------------------------------
% setting appropriately the figure for saving the file
set(h,'PaperPositionMode','Auto');
saveas(h, fullfile(temp_figure_folder,[model,'_Fig4_temp.fig']));
close all
delete(h);
% remove dir from path
rmpath('./_support_functions');