Scripting - Inconsistent muscle and ACL length compared to GUI

[Sven Paldauf]

Hi all,

I am currently trying to compute the ACL strain of the right knee ('Knee model with ligaments') via Matlab.
The procedure is to update the coordinates from a .mot file computed by the IK Tool and further to calculate the ACL length for each frame.
Since the GUI does not provide a function to calculate ACL strain, I also extract the muscle length of the biceps femoris computed by Matlab to compare with the results from the GUI.

See the code:

Code: Select all

import org.opensim.modeling.*
% Import the model
model = Model('MODELDIR\ResultsScale\orth_no_simbody_tester.osim');
state = model.initSystem();

&Load the .mot file generated by IK in the GUI
motSto = Storage('FILEDIR\IKResults\orth_no_gen_02_pointed_ik.mot');
nmotSto = motSto.getSize();
timeStep=motSto.getMinTimeStep();

coords = model.getCoordinateSet();
nCoords = coords.getSize();

%Update each model coordinate for each frame
for j=1:nmotSto

        for i=0:nCoords-16					% only update coordinates from right knee
            
            coordvalue = ArrayDouble();
            currentcoord = coords.get(i);
            motSto.getDataColumn(currentcoord.getName(),coordvalue);
            q = coordvalue.getitem(j-1);
            
            model.updCoordinateSet().get(i).setValue(state, q);
        
        end
        
    model.realizePosition(state);
    
%%%% aACL_r %%%%
% access ligament in models forceset; get ligament length after safedowncast
aACL_r = model.getForceSet().get(0);
aACL_r_as_lig = Ligament.safeDownCast(aACL_r);
aACL_r_length(j) = aACL_r_as_lig.getLength(state);

%%%% bicfemlh_r %%%%
bifemlh_r = model.getMuscles().get(3);
bifemlh_r_length(j) = bifemlh_r.getLength(state);

end

%% Plot results
%%%% bicfemlh_r %%%%
%  plot absolute length 

figure
plot(bifemlh_r_length);

%%%% aACL_r %%%%
% calculate strain, plot strain and absolute length 
aACL_r_L0 = aACL_r_as_lig.getRestingLength();
aACL_r_strain = (aACL_r_length-aACL_r_L0)/aACL_r_L0;

figure
plot(aACL_r_strain);

In Matlab, the following error appears, which is not present in the GUI:

AssemblySolver::track() attempt Failed: SimTK Exception thrown at assembler.cpp:973:
Method Assembler::track() failed because:
Unable to achieve required assembly error tolerance.
Assembly error tolerance achieved: 2.0932349631408442e-08 required: 1e-10.
Model unable to assemble: AssemblySolver::assemble() Failed: SimTK Exception thrown at assembler.cpp:897:
Method Assembler::assemble() failed because:
Unable to achieve required assembly error tolerance.
Assembly error tolerance achieved: 2.0939234790517958e-08 required: 1e-10.
Model relaxing constraints and trying again.

Also, ACL strains are way to high (>50%) and the length of the biceps femoris does not correspond to the GUI result (see ATTACHMENT).

It seems that the coordinates are updated incorrectly, but I don't know what I'm doing wrong here. I think it has also something to do with the Assembly error.

Can anyone help me with this issue?

Thanks to all and greetings from Germany,

Sven


Model:
https://simtk.org/projects/kneeligament

[jimmy d]

When you use the IKTool to generate a motion, it will report Rotational coordinates in degrees, but the model itself uses Radians. So when you are setting the coordinate values, you will have to convert back to Radians.

Code: Select all

...
for i=0:nCoords-16 % only update coordinates from right knee
            coordvalue = ArrayDouble();
            currentcoord = coords.get(i);
            motSto.getDataColumn(currentcoord.getName(),coordvalue);
            q = coordvalue.getitem(j-1);
            if strcmp('Rotational', char(currentcoord.getMotionType()))
                q = deg2rad(q);
            end
            model.updCoordinateSet().get(i).setValue(state, q);
        end
 ...

-James

[Sven Paldauf]

Exactly what I was looking for.
This also solved the assembly problem.

James, you're the best!

Thank you very much,

Sven

[Sven Paldauf]

One more question, sorry:)

After I implemented the lines from James,the code works.
Unfortunately, I observed that computation time is more than 2x higher now.
Since number of frames is around 3000, it takes quite a long time.

I assume that the conversion for each 'q' from deg to rad increases computation time.


Are there any ideas how to improve computation time in my code?

I thought about converting the Storage before entering the loop from deg to rad, but couldn't find a function in the Storage class. And not sure if this will improve time....

Thanks again, I really appreciate your help.

Sven

Code: Select all

import org.opensim.modeling.*
% Import the model
model = Model('MODELDIR\ResultsScale\orth_no_simbody_tester.osim');
state = model.initSystem();

%Load the .mot file generated by IK in the GUI
motSto = Storage('FILEDIR\IKResults\orth_no_gen_02_pointed_ik.mot');
nmotSto = motSto.getSize();
timeStep=motSto.getMinTimeStep();

coords = model.getCoordinateSet();
nCoords = coords.getSize();

%Update each model coordinate for each frame
for j=1:nmotSto

        for i=0:nCoords-16					% only update coordinates from right knee
            
            coordvalue = ArrayDouble();
            currentcoord = coords.get(i);
            motSto.getDataColumn(currentcoord.getName(),coordvalue);
            q = coordvalue.getitem(j-1);
            if strcmp('Rotational', char(currentcoord.getMotionType()))
                q = deg2rad(q);
            end
            model.updCoordinateSet().get(i).setValue(state, q);
        
        end
        
    model.realizePosition(state);
    
%%%% aACL_r %%%%
% access ligament in models forceset; get ligament length after safedowncast
aACL_r = model.getForceSet().get(0);
aACL_r_as_lig = Ligament.safeDownCast(aACL_r);
aACL_r_length(j) = aACL_r_as_lig.getLength(state);

%%%% bicfemlh_r %%%%
bifemlh_r = model.getMuscles().get(3);
bifemlh_r_length(j) = bifemlh_r.getLength(state);

end

[Thomas Uchida]

A potential issue is that you're projecting the solution onto the constraint manifold each time you're setting one of the Coordinates (line 26). You should only be enforcing the constraints when setting the last Coordinate (see the documentation for Coordinate::setValue() here: https://simtk.org/api_docs/opensim/api_ ... f2e88d7c42).

I assume that the conversion for each 'q' from deg to rad increases computation time.

You could try identifying the slow line(s) by either commenting them out one-by-one (does the script speed up?) or calling them twice (does the script slow down?). I wouldn't expect deg2rad() to be the bottleneck.

[Sven Paldauf]

Hey Thomas,

much better now, computation time significantly decreased!

That is what I was looking for.

Great how you guys manage to reply in only a few hours!

Thanks!