Extracting the angular velocity

[Nasma Dasser]

Dear Opensim community,

I am trying to extract a body (in this example the tibia_r) angular velocity using the SimBodyEngine on Python 2.7 without going through the Analyze Tool.
- The model is the Gait2354
- The states file consists of joint coordinates extracted from IK
- I realise the state variable derivatives
- the COM_global is correctly computed
- However the value obtained for the w_local is all zeroes

Has anyone faced this issue or am I missing a crucial step ?

Thank you very much for your help!

Best Regards

This is the Python code that I am running:

##############################################
model = opensim.Model('Gait2354.osim')
bodyset = model.getBodySet()
bodypart = bodyset.get('tibia_r')
states = 'IK_joint_coordinates.mot'

### COM in local
COM_local = opensim.Vec3(0)
COM_local = bodypart.getMassCenter()
ground_body = model.getGround()

MotStorage = opensim.Storage(states)
n_coord = model.getNumCoordinates()

COM1 = []
COM2 = []
COM3 = []
W1 = []
W2 = []
W3 = []

if MotStorage.isInDegrees():
r = radians
else:
r = lambda x: x

## Realize the state
for n_frame in range (0 , MotStorage.getSize()):

## initialize state
state = model.initSystem()

## loop through coordinate set
for i in range(n_coord):
name = model.getCoordinateSet().get(i).getName()
index = MotStorage.getStateIndex(name)
vector = MotStorage.getStateVector(n_frame).getData().get(index)
model.updCoordinateSet().get(i).setValue(state,r(vector))

## generate vector for COM in global
COM_global = opensim.Vec3(0)
w_local = opensim.Vec3(0)

## transform COM to global
model.getSimbodyEngine().transformPosition(state, bodypart, COM_local, ground_body, COM_global)

## realize to the acceleration stage
model.computeStateVariableDerivatives(state)

## get the angular velocity
model.getSimbodyEngine().getAngularVelocityBodyLocal(state, bodypart, w_local);

COM1.append([COM_global.get(0)])
COM2.append([COM_global.get(1)])
COM3.append([COM_global.get(2)])
W1.append([w_local.get(0)])
W2.append([w_local.get(1)])
W3.append([w_local.get(2)])
COM_glob = np.hstack((COM1,COM2,COM3))
W_loc = np.hstack((W1,W2,W3))

[Dimitar Stanev]

Maybe the initializeState call is causing the problem. You should initialize the state only once before the first for loop.

[Nasma Dasser]

Dear Dimitar,

Thank you for taking the time to reply to my questions. Unfortunately even if the state initialisation takes place before the first loop, the results for the angular velocity do not change and remain zero. When trying to compute afterwards:
model.getSimbodyEngine().getAcceleration(state, bodypart, COM_global, acc)

I do manage to extract some values for the acceleration but the angular velocity always remains zero. The angular acceleration also delivers zero values.

[Thomas Uchida]

It looks like you are setting the positions of the coordinates

Code: Select all

model.updCoordinateSet().get(i).setValue(state,r(vector))

but leaving the velocities set to zero. I would expect to see something like

Code: Select all

model.updCoordinateSet().get(i).setSpeedValue( #etc...

https://simtk.org/api_docs/opensim/api_ ... dd7a8febd3

Also, in general, you should enforce the position-level constraints only once all coordinates have been set. Please see the comments for the Coordinate::setValue() method: https://simtk.org/api_docs/opensim/api_ ... f2e88d7c42.

[Nasma Dasser]

Dear Thomas,

Thank you very much for taking the time to reply! Indeed, this was missing. After setting the speed value, I am able to extract the angular velocity. Important thing is to differentiate the coordinates. Here are the few lines of code I added if someone is interested in the same operations:

speed = np.diff(vector,0)
model.updCoordinateSet().get(i).setValue(state,r(vector))
model.updCoordinateSet().get(i).setSpeedValue(state,r(speed))

[Hossein Mokhtarzadeh]

Hello,

Thanks for this post. I am trying to get whole body CoM (center of mass) during a dynamic motion using Maltab Scripts.


I have done the following in Matlab and then I did Body Kinematics Analysis to get CoM position. However, the outcomes for CoM are different and Maltab Script takes a long time to calculate CoM. Here is what I used in Maltab. I wonder if you could help me understand where I am making the mistake. Thank you.

import org.opensim.modeling.*
model_folder = '\OpenSim\4.0\Code\Matlab';
model = Model(fullfile(model_folder, 'subject01_simbody.osim'));
state = model.initSystem();

motion_storage = Storage("W1_ik.mot");MotStorage = motion_storage;
motion_storage.getSize();
n_coord = model.getNumCoordinates();

% ## Realize the state
for n_frame = 0:MotStorage.getSize()-1

% ## initialize state
% state = model.initSystem();

for i=0:n_coord-1

name = model.getCoordinateSet().get(i).getName();
index = MotStorage.getStateIndex(name);
vector = MotStorage.getStateVector(n_frame).getData().get(index);
model.updCoordinateSet().get(i).setValue(state,vector);
model.computeStateVariableDerivatives(state);

end

p = model.calcMassCenterPosition(state);
% COMLoc(i+1,:,n_frame+1) = str2num(char(ArrayDouble.getValuesFromVec3(p)));
COMLoc(n_frame+1, = str2num(char(ArrayDouble.getValuesFromVec3(p)));

end

[Dimitar Stanev]

You can calculate the total CoM directly from the model:

https://simtk.org/api_docs/opensim/api_ ... 8a90aaf4cd

[Hossein Mokhtarzadeh]

Thanks Dimitar! I think I am using the same method you kindly shared. My problem were the following:

1. Why the results are different between this script (that I shared before here) and BodyKinematics analysis in OpenSim GUI?
2. Also, scripting takes a few second to calculate CoM while BodyKinematics analysis is much faster.

I appreciate any ideas.

Cheers,
Hossein

[Dimitar Stanev]

I think the problem is in the setValue, because vector variable is probably in degrees, while Simbody/OpenSim operates in radiance.

The execution time is different because the analyze tool is compiled and optimized for speed, whereas the Matlab script file is possibly not or in some case uses JIT compilation to optimize portions of code.