Calculating Segment Angular Velocities In the Global Coordinate System

[Steven Hirsch]

I have been trying to write a script that computes the angular velocities of all segments in my Opensim model in the global coordinate system using an existing scaled .osim model file and the joint angle results from the Inverse Kinematics Tool stored in a .mot file. I am currently using Opensim 4.4 in Python 3.8.16.

I have started by building a script using ChatGPT (as I'm not very familiar with the scripting tools in Opensim):

Code: Select all

import opensim as osim

# Load your OpenSim model
model = osim.Model("path/to/your/model.osim")

# Create a Storage object to read in the .mot file
storage = osim.Storage("path/to/your/model_outputs.mot")

# Get the number of bodies in the model
num_bodies = model.getNumBodies()

# Loop through each time step in the .mot file
for i in range(storage.getSize()):
    # Set the model state to the current time step
    time = storage.getTime(i)
    state = model.initSystem()
    osim.TimeSeriesTableUtilities.setAllValues(state, storage.getRowAtIndex(i))

    # Loop through all bodies in the model and get their global angular velocities
    for j in range(num_bodies):
        body = model.getBodySet().get(j)
        ang_vel = body.getAngularVelocityInGround(state)
        print("Time = {}, Body {}: global angular velocity = {}".format(time, j+1, ang_vel))

However, I have run into several errors with this script. I can bypass the error in the line

Code: Select all

time = storage.getTime(i)

by removing that error or replacing it with i/sampling_frequency. But, I then run into two errors on line

Code: Select all

osim.TimeSeriesTableUtilities.setAllValues(state, storage.getRowAtIndex(i))

: 1) AttributeError: module 'opensim' has no attribute 'TimeSeriesTableUtilities', and 2) 'Storage' object has no attribute 'getRowAtIndex'. I am now unsure of which attributes I should be using in order to return each segment's angular velocity in the global coordinate system.

Any feedback on how to modify the existing script, or perhaps use a different script entirely, would be much appreciated!

[Mohammadreza Rezaie]

Hi, Follow these steps:

1. read the model and call initSystem:

Code: Select all

model = osim.Model('out_scaled.osim')
state = model.initSystem()

2. read IK output, get time, filter the columns and convert them to Radians. Personally, I prefer using TimeSeriesTable rather than Storage:

Code: Select all

q = osim.TimeSeriesTable('out_ik.mot')
time = q.getIndependentColumn()
osim.TableUtilities().filterLowpass(q, 20, padData=True)
model.getSimbodyEngine().convertDegreesToRadians(q)
q.trim(time[0], time[-1])

3. Calculate the speed (dq/dt), you can use either numpy.gradient or OpenSim::GCVSpline:

Code: Select all

GCVS = osim.GCVSplineSet(q)
u = q.clone()
for i,column in enumerate(q.getColumnLabels()):
	updu = u.updDependentColumn(column)
	for j,jj in enumerate(time):
		updu[j] = GCVS.evaluate(i,1,jj) # first derivative

4. create an empty output, TimeSeriesTableVec3, and set bodies' name as columns' label:

Code: Select all

angVel = osim.TimeSeriesTableVec3()
angVel.setColumnLabels([i.getName() for i in model.getBodySet()])

5. for each time frame, update the coordinates' value and speed, call assemble followed by realizeVelocity:

Code: Select all

for i,ii in enumerate(time):

	value = osim.RowVector(q.getRowAtIndex(i))
	speed = osim.RowVector(u.getRowAtIndex(i))
	for j,coordinate in enumerate(model.updCoordinateSet()):
		coordinate.setValue(state, value[j], False)
		coordinate.setSpeedValue(state, speed[j])

	model.assemble(state)
	model.realizeVelocity(state)

6. in this loop, create an empty RowVectorVec3 with the size of model's bodies, get the angular velocity of each body in ground frame (in radians) and convert it to degrees(?), then append this row to the output:

Code: Select all

	vector = osim.RowVectorVec3(model.getBodySet().getSize())
	for j,body in enumerate(model.getBodySet()):
		values = body.getAngularVelocityInGround(state)
		vector[j] = values.scalarDivideEq(osim.SimTK_PI/180)
	angVel.appendRow(ii, vector)

7. add meta data to the output and write it to a file:

Code: Select all

angVel.addTableMetaDataString('inDegrees', 'yes')
angVel.addTableMetaDataString('nColumns', str(angVel.getNumColumns()))
angVel.addTableMetaDataString('nRows',    str(angVel.getNumRows()))

osim.STOFileAdapter().write(angVel.flatten(['_x','_y','_z']), 'angVel.sto')

8. finally, compare your values with Analyze>BodyKinematics (GUI) and see if they match.

Hope this helps,
Mohammadreza

[Steven Hirsch]

Thank you very much, Mohammadreza! This seems to work now.