emg and marker tracking

[Mohammad RahimiGoloujeh]

hello i am trying to track emg and markers at the same time i wonder if there is a way to assign different initial and final time for each reference file i.e. the maker set data and emg data as my emg data is processed and start from 0 to 1 s and the marker data i am going to use represent a stride that starts from 33s to 34 s

[Mohammad RahimiGoloujeh]

function exampleEMGTracking_answers
clear; close all; clc;
import org.opensim.modeling.*;
model = Model('JASC02_scaled.osim')
reference_track='Trial06_markers_filteredtest'
solutionname="newCodetest"
emgReference = TimeSeriesTable('walking12ms_jumpexcluded.sto');
w=5;
study = MocoStudy();
study.setName('marker_tracking_emg');
problem = study.updProblem();
modelProcessor = ModelProcessor(model);
modelProcessor.append(ModOpIgnoreTendonCompliance());
modelProcessor.append(ModOpTendonComplianceDynamicsModeDGF('implicit'));
modelProcessor.append(ModOpReplaceMusclesWithDeGrooteFregly2016());
model = modelProcessor.process();
model.initSystem();
problem.setModel(model);
%% marker
markerTrackingCost = MocoMarkerTrackingGoal();
markerTrackingCost.setName('marker_tracking');
markerWeights = SetMarkerWeights();
markerWeights.cloneAndAppend(MarkerWeight("RASI", 20));
markerWeights.cloneAndAppend(MarkerWeight("LASI", 20));
markerWeights.cloneAndAppend(MarkerWeight("RPSI", 20));
markerWeights.cloneAndAppend(MarkerWeight("LPSI", 20));
markerWeights.cloneAndAppend(MarkerWeight("RKNE", 10));
markerWeights.cloneAndAppend(MarkerWeight("RANK", 10));
markerWeights.cloneAndAppend(MarkerWeight("RHEE", 10));
markerWeights.cloneAndAppend(MarkerWeight("RTOE", 2));
markerWeights.cloneAndAppend(MarkerWeight("LKNE", 10));
markerWeights.cloneAndAppend(MarkerWeight("LANK", 10));
markerWeights.cloneAndAppend(MarkerWeight("LHEE", 10));
markerWeights.cloneAndAppend(MarkerWeight("LTOE", 2));
markerWeights.cloneAndAppend(MarkerWeight("C7", 10));
markersRef = MarkersReference(strcat(reference_track,'.trc'), ...
markerWeights);
markerTrackingCost.setMarkersReference(markersRef);
markerTrackingCost.setAllowUnusedReferences(true);
problem.addGoal(markerTrackingCost);
%% EMG
tracking = MocoControlTrackingGoal('emg_tracking');
tracking.setWeight(5);
tracking.setReference(TableProcessor(emgReference));
tracking.setReferenceLabel('/forceset/tib_ant_r', 'tib_ant_r');
tracking.setReferenceLabel('/forceset/lat_gas_r', 'lat_gas_r');
tracking.setReferenceLabel('/forceset/vas_lat_r', 'vas_lat_r');
tracking.setReferenceLabel('/forceset/rect_fem_r', 'rect_fem_r');
tracking.setReferenceLabel('/forceset/bifemlh_r', 'bifemlh_r');
tracking.setReferenceLabel('/forceset/add_mag1_r', 'add_mag1_r');
tracking.setReferenceLabel('/forceset/glut_max1_r', 'glut_max1_r');
tracking.setReferenceLabel('/forceset/tfl_r', 'tfl_r');
tracking.setReferenceLabel('/forceset/glut_med1_r', 'tfl_r');
tracking.setReferenceLabel('/forceset/med_gas_r', 'med_gas_r');
tracking.setReferenceLabel('/forceset/soleus_r', 'soleus_r');
tracking.setReferenceLabel('/forceset/per_long_r', 'per_long_r');

tracking.setReferenceLabel('/forceset/tib_ant_l', 'tib_ant_l');
tracking.setReferenceLabel('/forceset/lat_gas_l', 'lat_gas_l');
tracking.setReferenceLabel('/forceset/vas_lat_l', 'vas_lat_l');
tracking.setReferenceLabel('/forceset/rect_fem_l', 'rect_fem_l');
tracking.setReferenceLabel('/forceset/bifemlh_l', 'bifemlh_l');
tracking.setReferenceLabel('/forceset/add_mag1_l', 'add_mag1_l');
tracking.setReferenceLabel('/forceset/glut_max1_l', 'glut_max1_l');
tracking.setReferenceLabel('/forceset/tfl_l', 'tfl_l');
tracking.setReferenceLabel('/forceset/glut_med1_l', 'glut_med1_l');
tracking.setReferenceLabel('/forceset/med_gas_l', 'med_gas_l');
tracking.setReferenceLabel('/forceset/soleus_l', 'soleus_l');
tracking.setReferenceLabel('/forceset/per_long_l', 'per_long_l');

tracking.addScaleFactor('tib_ant_r_factor', '/forceset/tib_ant_r', [0.01 1.0]);
tracking.addScaleFactor('lat_gas_r_factor', '/forceset/lat_gas_r', [0.01 1.0]);
tracking.addScaleFactor('vas_lat_r_factor', '/forceset/vas_lat_r', [0.01 1.0]);
tracking.addScaleFactor('rect_fem_r_factor', '/forceset/rect_fem_r', [0.01 1.0]);
tracking.addScaleFactor('bifemlh_r_factor', '/forceset/bifemlh_r', [0.01 1.0]);
tracking.addScaleFactor('add_mag1_r_factor', '/forceset/add_mag1_r', [0.01 1.0]);
tracking.addScaleFactor('glut_max1_r_factor', '/forceset/glut_max1_r', [0.01 1.0]);
tracking.addScaleFactor('tfl_r_factor', '/forceset/tfl_r', [0.01 1.0]);
tracking.addScaleFactor('glut_med1_r_factor', '/forceset/glut_med1_r', [0.01 1.0]);
tracking.addScaleFactor('med_gas_r_factor', '/forceset/med_gas_r', [0.01 1.0]);
tracking.addScaleFactor('soleus_r_factor', '/forceset/soleus_r', [0.01 1.0]);
tracking.addScaleFactor('per_long_r_factor', '/forceset/per_long_r', [0.01 1.0]);

tracking.addScaleFactor('tib_ant_l_factor', '/forceset/tib_ant_l', [0.01 1.0]);
tracking.addScaleFactor('lat_gas_l_factor', '/forceset/lat_gas_l', [0.01 1.0]);
tracking.addScaleFactor('vas_lat_l_factor', '/forceset/vas_lat_l', [0.01 1.0]);
tracking.addScaleFactor('rect_fem_l_factor', '/forceset/rect_fem_l', [0.01 1.0]);
tracking.addScaleFactor('bifemlh_l_factor', '/forceset/bifemlh_l', [0.01 1.0]);
tracking.addScaleFactor('add_mag1_l_factor', '/forceset/add_mag1_l', [0.01 1.0]);
tracking.addScaleFactor('glut_max1_l_factor', '/forceset/glut_max1_l', [0.01 1.0]);
tracking.addScaleFactor('tfl_l_factor', '/forceset/tfl_l', [0.01 1.0]);
tracking.addScaleFactor('glut_med1_l_factor', '/forceset/glut_med1_l', [0.01 1.0]);
tracking.addScaleFactor('med_gas_l_factor', '/forceset/med_gas_l', [0.01 1.0]);
tracking.addScaleFactor('soleus_l_factor', '/forceset/soleus_l', [0.01 1.0]);
tracking.addScaleFactor('per_long_l_factor', '/forceset/per_long_l', [0.01 1.0]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

problem.addGoal(tracking)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Control effort term (minimize squared muscle excitations)
effort= MocoControlGoal();
effort.setWeight(w);
effort.setExponent(2);
effort.setDivideByDisplacement(false);
problem.addGoal(effort)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GRF
contactTracking = MocoContactTrackingGoal('contact',6/(75.337*9.81)/(0.1073^2)/29);
contactTracking.setExternalLoadsFile('external_loads.xml');
forceNamesRightFoot = StdVectorString();
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s1_r');
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s2_r');
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s3_r');
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s4_r');
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s5_r');
forceNamesRightFoot.add('/forceset/SmoothSphereHalfSpaceForce_s6_r');
trackRightGRF = MocoContactTrackingGoalGroup(forceNamesRightFoot,'Right_GRF');
trackRightGRF.append_alternative_frame_paths('/bodyset/toes_r');
contactTracking.addContactGroup(trackRightGRF);
forceNamesLeftFoot = StdVectorString();
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s1_l');
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s2_l');
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s3_l');
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s4_l');
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s5_l');
forceNamesLeftFoot.add('/forceset/SmoothSphereHalfSpaceForce_s6_l');
trackLeftGRF = MocoContactTrackingGoalGroup(forceNamesLeftFoot,'Left_GRF');
trackLeftGRF.append_alternative_frame_paths('/bodyset/toes_l');
contactTracking.addContactGroup(trackLeftGRF);
problem.addGoal(contactTracking);

solver = MocoCasADiSolver.safeDownCast(study.updSolver());
solver.resetProblem(problem);

solver.set_parameters_require_initsystem(false);

%% Define the solver and set its options
solver.set_optim_max_iterations(3000);
solver.set_num_mesh_intervals(50);
solver.set_optim_constraint_tolerance(1e-04);
solver.set_optim_convergence_tolerance(1e+01);
solver.set_minimize_implicit_auxiliary_derivatives(true)
solver.set_implicit_auxiliary_derivatives_weight(0.00001)
solver.set_enforce_constraint_derivatives(true)
solver.resetProblem(problem);

solver.set_parameters_require_initsystem(false);

gaitTrackingSolution = study.solve();
gaitTrackingSolution.write(strcat(solutionname,'.sto'));
externalForcesTableFlat = opensimMoco.createExternalLoadsTableForGait(model,gaitTrackingSolution,forceNamesRightFoot,forceNamesLeftFoot);
STOFileAdapter.write(externalForcesTableFlat, strcat(solutionname,'_GRF.sto'))

% Part 3h: Get the values of the optimized scale factors.
trackingSolution = MocoTrajectory('trackingSolution.sto');
gastroc_factor = trackingSolution.getParameter('gastroc_factor');
tibant_factor = trackingSolution.getParameter('tibant_factor');
bifem_factor = trackingSolution.getParameter('bifem_factor');
gluteus_factor = trackingSolution.getParameter('gluteus_factor');

%% Part 4: Plot the EMG-tracking muscle redundancy problem solution.
% Part 4a: Print the scale factor values to the command window.
fprintf('\n')
fprintf('Optimized scale factor values: \n')
fprintf('------------------------------ \n')
fprintf(['gastrocnemius = ' num2str(gastroc_factor) '\n'])
fprintf(['tibialis anterior = ' num2str(tibant_factor) '\n'])
fprintf(['biceps femoris short head = ' num2str(bifem_factor) '\n'])
fprintf(['gluteus = ' num2str(gluteus_factor) '\n'])
fprintf('\n')

% Part 4b: Re-scale the reference data using the optimized scale factors.
gastroc = emgReference.updDependentColumn('gastrocnemius');
tibant = emgReference.updDependentColumn('tibialis_anterior');
bifem = emgReference.updDependentColumn('biceps_femoris');
gluteus = emgReference.updDependentColumn('gluteus');
for t = 0:emgReference.getNumRows() - 1
gastroc.set(t, gastroc_factor * gastroc.get(t));
tibant.set(t, tibant_factor * tibant.get(t));
bifem.set(t, bifem_factor * bifem.get(t));
gluteus.set(t, gluteus_factor * gluteus.get(t));
end

% Part 4c: Generate the plots. Compare results to the effort minimization
% solution.
compareSolutionToEMG(emgReference, 'effortSolution.sto', ...
'trackingSolution.sto');

end

function compareSolutionToEMG(varargin)

import org.opensim.modeling.*;

% Retrieve the inputs.
emgReference = varargin{1};
effortSolution = MocoTrajectory(varargin{2});
if nargin > 2
trackingSolution = MocoTrajectory(varargin{3});
end

% Create a time vector for the EMG data that is consistent with the problem
% time range.
time = effortSolution.getTimeMat();
startIndex = emgReference.getNearestRowIndexForTime(time(1)) + 1;
endIndex = emgReference.getNearestRowIndexForTime(time(end)) + 1;
emgTime = linspace(time(1), time(end), endIndex - startIndex + 1);

% Plot results from the muscles in the left leg.
titles = {'soleus', 'gastroc.', 'tib. ant.', ...
'bi. fem. sh.', 'vastus', 'rec. fem.', 'gluteus', 'per_long'};
emgSignals = {'soleus', 'gastrocnemius', 'tibialis_anterior', ...
'biceps_femoris', 'vastus', ...
'rectus_femoris', 'gluteus', 'perlong'};
muscles = {'soleus_l', 'med_gas_l', 'tib_ant_l', 'bifemlh_l', ...
'vas_lat_l', 'rec_fem_l', 'glut_max1_l', 'per_long_l'};
figure;
for i = 1:length(titles)
subplot(3,3,i)

% If it exists, plot the EMG signal for this muscle.
if ~strcmp(emgSignals{i}, 'none')
col = emgReference.getDependentColumn(emgSignals{i}).getAsMat();
emg = col(startIndex:endIndex)';
patch([emgTime fliplr(emgTime)], [emg fliplr(zeros(size(emg)))], ...
[0.7 0.7 0.7])
alpha(0.5);
end
hold on

% Plot the control from the muscle effort solution.
control = effortSolution.getControlMat(['/forceset/' muscles{i}]);
plot(time, control, '-k', 'linewidth', 2)

% If provided, plot the control from the tracking solution.
if nargin > 2
control = trackingSolution.getControlMat(['/forceset/' muscles{i}]);
plot(time, control, '-r', 'linewidth', 2)
end

% Plot formatting.
title(titles{i})
xlim([0.83 2.0])
ylim([0 1])
if i == 1
if nargin > 2
legend('EMG', 'effort', 'tracking', 'location', 'best')
else
legend('EMG', 'effort', 'location', 'best')
end
end
if i > 6
xlabel('time (s)')
else
xticklabels([])
end
if ~mod(i-1, 3)
ylabel('excitation')
else
yticklabels([])
end
end


end

[Aaron Fox]

Hi Mohammad,

I don't think the offset times will work and there is now way to set different times from what I understand. I think OpenSim will always look for corresponding time stamps across all different data types. The easy solution here is to edit/synchronise the original data files so that the times match (i.e. figure out the offset from marker to EMG data and then subtract that from the time so they sync up).

Aaron

[Mohammad RahimiGoloujeh]

thank you so much,
another scenario I am dealing with is that I i have a reference motion that is corresponding to 1.45m/s with 1s gait cycle duration and I want to use that as a reference to generate 1.1m/s gait with 1.15s cycle duration. I was wondering how I can handle that