EDIT: this code snippet only works for MOCO 0.4, if you've updated to 1.0 you need to modify according to the posts by Nick & Ross later on in this thread.
Here is a snippet of Matlab code that extracts all the jointreactions from a GaitPredictionSolution, expresses them in the child body frame, and places them in a structured array, with tidied up column headers as a seperate cell array (for easier plotting/processing). You should be able to place this directly inline with the script that runs a simulation, or you could place this in a function (with inputs: model, study, and a predictionSolution). If you'd like the reactions on parent instead, change every instance of child (or Child) to parent (or Parent).
Hope this helps someone!
Code: Select all
JR_paths=StdVectorString();
JR_paths.add('.*reaction_on_child'); %regexp selecting any output that ends with this expression, so this selects all the joints
JR_outputs_table = opensimMoco.analyzeSpatialVec(model, gaitPredictionSolution, JR_paths).flatten(); %Flatten turns a TimeSeriesTableSpatialVec into a TimeSeriesTable, with 6 columns for each spatialvec
% A SpatialVec in this context provides the following outputs: Mx My Mz Fx Fy Fz
% Get a Matlab struct
JR_outputs = osimTableToStruct (JR_outputs_table);
JRHeaders= fieldnames(JR_outputs);
%%% Place in an array
for i = 1:length(JRHeaders) -1
Joint_reactions (:,i) = JR_outputs.(JRHeaders{i,1}); %These are still expressed in the ground frame
end
statesTraj= gaitPredictionSolution.exportToStatesTrajectory(model); %get a StatesTrajectory, so we can loop through it
[~,n] = size(Joint_reactions);
for i = 1:statesTraj.getSize %loop through each timepoint in the trajectory
state=statesTraj.get(i-1);
model.realizeDynamics(state);
for j = 1:(n/6) %Loop through each joint (joint reactions provide a spatialvec as output, so 6 columns per joint)
Torque_ground = Vec3(Joint_reactions(i,j*6 -5), Joint_reactions(i,j*6 -4), Joint_reactions(i,j*6 -3));
Force_ground = Vec3(Joint_reactions(i,j*6 -2), Joint_reactions(i,j*6 -1), Joint_reactions(i,j*6 ));
ground = model.getGround();
child_frame = model.getJointSet().get(j-1).getChildFrame; %it is worth double checking if these frames
%correspond to the right joints as in JR_outputs_table
% Express the torques and forces in the child frame
Torque_local = ground.expressVectorInAnotherFrame(state, Torque_ground, child_frame);
Force_local = ground.expressVectorInAnotherFrame(state, Force_ground, child_frame);
%place back in an array.
joint_reactions_local(i,j*6 -5) = Torque_local.get(0);
joint_reactions_local(i,j*6 -4) = Torque_local.get(1);
joint_reactions_local(i,j*6 -3) = Torque_local.get(2);
joint_reactions_local(i,j*6 -2) = Force_local.get(0);
joint_reactions_local(i,j*6 -1) = Force_local.get(1);
joint_reactions_local(i,j*6 ) = Force_local.get(2);
end
end
%Improve formatting of the headers
JRHeaders = strrep(JRHeaders,'unlabeled_jointset_','');
JRHeaders = strrep(JRHeaders,'_',' ');
JRHeaders = strrep(JRHeaders,'reaction on child 1','Mx on child');
JRHeaders = strrep(JRHeaders,'reaction on child 2','My on child');
JRHeaders = strrep(JRHeaders,'reaction on child 3','Mz on child');
JRHeaders = strrep(JRHeaders,'reaction on child 4','Fx on child');
JRHeaders = strrep(JRHeaders,'reaction on child 5','Fy on child');
JRHeaders = strrep(JRHeaders,'reaction on child 6','Fz on child');
%place in a struct
data.JRHeaders = JRHeaders(1:end-1,:); %remove the last entry which corresponded to the time vector
data.JR = joint_reactions_local; %joint reactions on the child body, expressed in the child body frame