Retrieve muscle parameters using API

[Mohammadreza Rezaie]

Dear Experts,

I hope all is well.

I'm trying to calculate a set of muscle properties such as muscle moment arm and strength using OpenSim 4.4 API in MATLAB (for custom static optimization).

Here is the script:

Code: Select all

% q2 ==> Kinematics_q.sto file (2D array [time,coordinate])
% u2 ==> Kinematics_u.sto file (2D array [time,coordinate])

model = Model('subject01_simbody.osim');
state = model.initSystem();
muscles = model.updMuscles();
nMuscles = muscles.getSize();
coordinates = model.getCoordinateSet();
nCoordinates = coordinates.getSize();

for i=1:frameTime
	Q = Vector(); Q.resize(nCoordinates);
	U = Vector(); U.resize(nCoordinates);
    for j=1:nCoordinates
		Q.set(j-1, q2(i,j));
		U.set(j-1, u2(i,j));
    end
	state.setQ(Q);
	state.setU(U);
    
%     model.realizePosition(state);
%     model.realizeVelocity(state);
%     model.equilibrateMuscles(state);
    model.realizeDynamics(state);

    for j=1:nMuscles
        muscle = Millard2012EquilibriumMuscle.safeDownCast(muscles.get(j-1));
        muscle.set_ignore_activation_dynamics(true); % disable activation dynamics
        muscle.set_ignore_tendon_compliance(true); % disable tendon compliance

%         muscle.computeFiberEquilibrium(state);
%         muscle.computeEquilibrium(state);
%         muscle.computeActuation(state);
        model.equilibrateMuscles(state);

        L(i,j)  = muscle.getLength(state);
        FL(i,j) = muscle.getFiberLength(state);
        TL(i,j) = muscle.getTendonLength(state);
        OFL(j)  = muscle.getOptimalFiberLength();
        TSL(j)  = muscle.getTendonSlackLength();
        TL(i,j) = muscle.getTendonLength(state);
        MIF(j)  = muscle.getMaxIsometricForce();
        FF(i,j) = muscle.getFiberForce(state);
        AFF(i,j)= muscle.getActiveFiberForce(state);
        PFF(i,j)= muscle.getPassiveFiberForce(state);
        TF(i,j) = muscle.getTendonForce(state);

        muscle.setActivation(state, 1);
        S(i,j) = muscle.getFiberForce(state); % muscle strength

%         for k=1:nCoordinates
% 			coordinate = coordinates.get(k-1);
% 			coordinate.setValue(state, q2(i,k))
% 			MA(i,k,j) = muscle.computeMomentArm(state, coordinate);
%         end
    end
end

This is the result of soleus_r strength for a stance phase of the right foot in MATLAB which seems incorrect:

MatLab_soleus.jpg (35.72 KiB) Viewed 505 times

But in Python, same script gives me this plot:

Python_soleus.jpg (20.61 KiB) Viewed 505 times

I'm a little confused here, because the snippets are the same. I was wondering if anyone could guide me whether the snippet is correct and why the results are not identical.

Here are the complete files:

test.zip

(109.02 KiB) Downloaded 13 times

Any help is much appreciated.

Regards,
Mohammadreza

[Ayman Habib]

Hello,

While it's hard to debug specific code snippets extracted from a code block, generally you should be aware of the following:

Methods that operates on properties (e.g. set_ignore_activation_dynamics) typically invalidate the system under the model as they can lead to a change in the state layout or change in the number of state variables. Accordingly, trying to reuse a state across these calls is not supported and will lead to unpredictable results.

I'd be surprised if this behavior is different between Matlab and Python but I'd address this issue first.

Best regards,
-Ayman

[Varun Joshi]

Hi,

Your python code converts input coordinates and velocities to radians from degrees, but the Matlab code doesn't.

If you change lines 35 - 38 to -

Code: Select all

for i=1:nCoordinates
    q2(:,i) = q.(nameCoordinates{i})(timeBool)*2*pi/180;
    u2(:,i) = u.(nameCoordinates{i})(timeBool)*2*pi/180;
end

then the Matlab code produces the same figure.

[Mohammadreza Rezaie]

Dear Ayman and Varun,

I deeply appreciate your responses. My bad, I should have paid attention to the units. All solved.

Also, I put the set_ignore_tendon_compliance() and set_ignore_activation_dynamics() at the beginning of the code and recalled the state (model.initSystem(), there was no difference.
But the results of using setIgnoreActivationDynamics() and setIgnoreTendonCompliance() were different. Which one of these must be used for this purpose (customStaticOptimization practice)?