How to set muscle force to calculate moment arms?

[Marc Carmichael]

Hi all,

I am trying to calculate the muscle moment arms using the method described in document "How to compute muscle moment arm using generalized coordinates"
http://wiki.simtk.org/opensim/FrontPage ... entArm.pdf

I am using the Generalized force method as this seems to be the only way to accurately calculate moment arms for complex models without using the finite different method (is this right?).

My problem is that the approach shown in the code on page 9 uses the function

Code: Select all

muscle.setTension()

which appears not to exist, maybe depreciated? To get around this I set the muscle's activation to 1, then save the resulting muscle force. However this is not robust as if the muscle/tendon is slack then the activation produces negligible force and the computation fails.

Is there another way to manually set the muscle tension? I tried using

Code: Select all

muscle.setForce()

but it did not appear to effect the calculated Spatial forces. I am using Holzbaur's upper limb model that uses Schutte1993 muscle models.

Note: I am using the multiplyBySystemJacobianTranspose() function as the calcInternalGradientFromSpatial() does not exist either, are they the same thing?

Any help is greatly appreciated.

Marc

[Michael Sherman]

Hi, Marc.

That paper was used internally as a guide for OpenSim's moment arm calculation that is already available, debugged and tested, in OpenSim 3.0. Is there some reason you want to calculate it yourself?

If so, you can look in the OpenSim source code to find the moment arm calculation -- Ajay Seth can probably point you in the right direction.

And yes, the obscure name calcInternalGradientFromSpatial() was replaced in Simbody 3.0 by the hopefully-less-obscure name multiplyBySystemJacobianTranspose() which is what it actually does!

Regards,
Sherm

[Marc Carmichael]

Hi Sherm,

Sorry for the late reply. I want to explicitly calculate (as quick as possible) for a given pose a matrix containing the moment arms. I use it in my work to relate upper limb muscle strength to the strength at the hand in different directions.

I was using OpenSim's computeMomentArm() function but I found that calculating it using the method in the pdf was actually quicker. I found it very helpful, and have since utilized the coupling matrix "C" to also calculate kinematic Jacobians a lot quicker than I was able to previously.

However I still am unable to set the force output of the individual muscle-tendon actuators to be able to apply it robustly to calculate the moment arms in a similar way. My feeling is that my setting of the muscle force is being over written at some point when the state is realized.

I will have a look at the source code and see exactly how the pros do it

Marc

[Ajay Seth]

You can also use the MomentArmSolver. This is what the GeomtryPath uses to compute its moment-arm. It is not tied to any Force or Actuator and computes the effective moment-arm based on the set of point forces distributed over multiple bodies. You can ask a GeometryPath to supply the set of point force directions, which handles issues of wrapping and conditional path points, or you can create them yourself for non GeometryPath based applications of force, for example from a constraint or gravity.

[Morgan Sangeux]

Hi,

How do you access the MomentArmSolver?

Thanks,

[Marc Carmichael]

Thanks for all the information. I have used a number of different methods in the past. I have used both the MomentArmSolver as well as the muscle.computeMomentArm() function, and both worked well.

However for my application I found that it was quicker for me to implement my own routine which basically reimplements the MomentArmSolver but only for the muscles and coordinates of interest to me. I have just changed the method I detailed in my original post to use the "Point Force Directions" of the muscle geometry path as suggested by Ajay. It works well and is quick.

FYI here are some benchmarks. They are VERY rough (acquired via a MATLAB mex function on an old PC) but they provide an interesting comparison:

• 1) Using manual finite difference method = 8.3 s
• 2) Using muscle.computeMomentArm() = 19.9 s
• 3) Using MomentArmSolver = 15.2 s
• 4) Using method in my first post = 0.74 s
• 5) Using latest method with PFDs = 0.80 s

This was using Holzabur's upper limb model, calculating moment arms for 37 muscles over 4 coordinates (148 moment arms).

I think my method for implementing it is much quicker for me as I only calculate the model's coupling matrix once for each coordinate of interest, then reuse it for each of the muscles. Interestingly method 4) where I set the muscle's activation to 1 then 0 and use the difference in the generalized forces is always slightly quicker than method 5).

The majority of the time for methods 4) and 5) is spent calculating the coupling matrix "C" (around 0.57 seconds in this case). I would be very interested if you guys know of any way of making that calculation quicker?

Marc

[Marc Carmichael]

Hi,

How do you access the MomentArmSolver?

Thanks,

Hi Morgan,

Here is some sample code that implements the MomentArmSolver. The iGenCoords parameter is a SimTK array of ints containing the indicies of the joints I wish to calculate the moment arms for. It returns a [nMuscles x nJoints] matrix of moment arms.

Code: Select all

// Compute muscle Jacobian L using method 2 (OpenSim's MomentArmSolver)
SimTK::Matrix calcL2(OpenSim::Model& model, SimTK::State& si, SimTK::Array_<int> iGenCoords)
{
	int nGenCoords = iGenCoords.size();
	int nMuscles = model.getMuscles().getSize();
	SimTK::State s = si;
	model.getMultibodySystem().realize(s,SimTK::Stage::Position);

	SimTK::Matrix L(nMuscles,nGenCoords);
	
	OpenSim::MomentArmSolver momentArmSolver(model);
    OpenSim::Array<OpenSim::PointForceDirection *> pfds;
	
	// Cycle through muscles
    for (int m = 0 ; m < nMuscles ; m++) { 
        pfds.setSize(0);
        model.getMuscles().get(m).getGeometryPath().getPointForceDirections(s,&pfds);

        for (int q = 0 ; q < nGenCoords ; q++ ) { // Cycle through gen coords
            const OpenSim::Coordinate &coord = model.getCoordinateSet().get(iGenCoords[q]);
            L[m][q] = - momentArmSolver.solve(s,coord,pfds);
        }
    }
    return L;
}

Hope this helps.

Marc

[Michael Sherman]

Hi, Marc. I'm not sure how to calculate C faster, but that doesn't mean there isn't a way to do it. If the velocity-level constraints are given by G(q)*u=0, then C is the vector in the null space of G that has the u-of-interest ui set to 1. It is unique because we assume kinematic coupling among all the relevant u's for a muscle.

We use Simbody's velocity projection method projectU() to solve for the vector of interest. Given a u with ui=1 that doesn't satisfy the constraints, it calculates a least squares change du to u that does satisfy the constraints: G du = G u by calculating du = pinv(G)*G*u where pinv() is the pseudoinverse (calculated in our case with a complete orthogonal factorization of G). Then C = (u-du)/ui.

For a large batch of muscles and lots of constraints it would make sense to calculate the pseudoinverse once and use it repeatedly with different right hand sides. It's hard to say how much that would save though because it isn't clear what fraction of the time is spent in projectU(). Generally G doesn't have many rows so the pinv() calculation doesn't necessarily take very long. If you want to explore this further, please do and let us know what you find. You can get the constraint matrix G from Simbody's SimbodyMatterSubsystem::calcG() method described here.

Regards,
Sherm

[Ajay Seth]

Defintitely computing the kinematic coupling matrix for each muscle is not as efficient as computing it once for all muscles that involve that coordinate. The main use case was getting moment-arms for handful of muscles. It would be a good addition to the MASolver to efficient provide users for the moment arm of all muscles. Once you get the moment-arm matrix, is the first thing you do is multiply it with the column vector of muscle forces to get muscle "moments" (generalized forces)?

If so, a faster approach would be to replace Muscles with PathActuators that have the same GeometryPath as the muscles and provide the tension as a control input. It seems all you want from OpenSim is the muscle geometry and not muscle dynamics. After applying your computed tension to each actuator, you could can ask the multibody system for the generalized forces (getMobilityForces) or get the accelerations (udot). This would be a lot quicker unless you really need to explicitly form the moment arm matrix for some other reason. Please, let us know.