Upd os4.x

MATLAB_tool/MuscleParOptTool/getBodyJoint.m

@@ -0,0 +1,71 @@
+%-------------------------------------------------------------------------%
+%    Copyright (c) 2021 Modenese L.                                       %
+%    Author:   Luca Modenese,  2021                                       %
+%    email:    [email protected]                                  %
+% ----------------------------------------------------------------------- %
+% Function replacing the old getJoint() method that we all loved in OpenSim
+% 3.3 but currently non available in OpenSim 4.x.
+% Returns the joint for which the specified body would be the child.
+% ----------------------------------------------------------------------- %
+function bodyJoint = getBodyJoint(osimModel, aBodyName, debug_printout)
+
+import org.opensim.modeling.*;
+
+if nargin<3; debug_printout=0; end
+
+% default
+bodyJoint = [];
+
+if strcmp(aBodyName, 'ground')
+    return
+end
+
+% check if body is included in the model
+if osimModel.getBodySet().getIndex(aBodyName)<0
+    error(['getBodyJoint.m The specified body ', aBodyName,' is not included in the OpenSim model'])
+end
+
+% get jointset
+jointSet = osimModel.getJointSet();
+
+% loop through jointset
+nj = 1;
+
+for n_joint = 0:jointSet.getSize()-1
+
+    % get cur joint
+    cur_joint = jointSet.get(n_joint);
+
+    % child frame from joint
+    child_frame = cur_joint.getChildFrame();
+
+    % link back to base frame: this could be a body
+    body_of_frame = child_frame.findBaseFrame();
+
+    % get base frame name
+    possible_body_name = char(body_of_frame.getName());
+
+    % if body with that name exist than the joint is that body's joint
+    if osimModel.getBodySet.getIndex(possible_body_name)>=0 && strcmp(aBodyName, possible_body_name)
+
+        % save the joints with the specified body as Child 
+        jointName(nj) = {char(cur_joint.getName())};
+
+        if debug_printout
+            disp([aBodyName, ' is parent frame on joint: ',jointName{nj}]);
+        end
+        % update counter
+        nj = nj + 1;
+        continue
+    end
+end
+
+% return a string if there is only one body
+if numel(jointName)==1
+    jointName = jointName{1};
+    bodyJoint = osimModel.getJointSet.get(jointName);
+else
+    error(['getBodyJoint.m More than one joint connected to body of interest', aBodyName,'. This function is design to work as getJoint() in OpenSim 3.3.'])
+end
+
+end

MATLAB_tool/MuscleParOptTool/getOpenSimVersion.m

@@ -0,0 +1,36 @@
+%-------------------------------------------------------------------------%
+%    Copyright (c) 2021 Modenese L.                                       %
+%    Author:   Luca Modenese,  2021                                       %
+%    email:    [email protected]                                  %
+% ----------------------------------------------------------------------- %
+% This script returns the OpenSim version
+% ----------------------------------------------------------------------- %
+function [osim_version_float, osim_version_string] = getOpenSimVersion()
+
+import org.opensim.modeling.*
+
+% read env variable
+try
+    % method available only from 4.x
+    os_version = char(opensimCommonJNI.GetVersion());
+
+    % get the field separators, e.g. '-
+    sep_set = strfind(os_version, '-');
+
+    % use the last file separator to get the OpenSim installation folder name
+    version = os_version(1:sep_set(1));
+
+    % get the string for the opensim version
+    osim_version_string = strtrim(strrep(lower(version),'opensim', ''));
+
+    % transform in float
+    osim_version_float = str2double(osim_version_string);
+
+catch
+
+    % GetVersion is not available in earlier OpenSim versions
+    osim_version_float = 3.3;
+    osim_version_string = '3.3';
+    return
+end
+end

MATLAB_tool/MuscleParOptTool/private/getJointsSpannedByMuscle.m

@@ -42,10 +42,6 @@
 BodySet = osimModel.getBodySet();
 muscle  = osimModel.getMuscles.get(OSMuscleName);
 
-% additions BAK
-% load a jointStrucute detailing bone and joint configurations
-jointStructure = getModelJointDefinitions(osimModel);
-
 % Extracting the PathPointSet via GeometryPath
 musclePath = muscle.getGeometryPath();
 musclePathPointSet = musclePath.getPathPointSet();
@@ -88,38 +84,49 @@
 n_spanJointNoDof = 1;
 NoDofjointNameSet = {};
 jointNameSet = {};
+disp(['   spanned joints: '])
 while ~strcmp(bodyName,ProximalBodyName)
 
-    %spannedJoint = body.getJoint();
-    %spannedJointName = char(spannedJoint.getName());
+    if getOpenSimVersion()<4.0
+        spannedJoint = body.getJoint();
+    else
+        spannedJoint = getBodyJoint(osimModel, char(body.getName()), 0);
+    end
+
+    % check if spannedJoint is empty (meaning it is "ground")
+    if isempty(spannedJoint)
+        warndlg(['The path of muscle ',OSMuscleName,' has reached ground without finding its origin, meaning that its path is not going proximal only. This case is not supported. Please check the logs and assess if the computations are acceptable for your case.'])
+        return
+    end
+
+    spannedJointName = char(spannedJoint.getName());
+    disp(['    * ', spannedJointName]);
 
-    % BAK implementation
-    spannedJointName = getChildBodyJoint(jointStructure, body.getName());
-    spannedJoint = osimModel.getJointSet().get(spannedJointName);
-
     if strcmp(spannedJointName, spannedJointNameOld)
-         %body =  spannedJoint.getParentBody();
-         % BAK implementation
-         body = jointStructure.spannedJoint.parentBody;
-
+         body =  spannedJoint.getParentBody();
          spannedJointNameOld = spannedJointName;
     else
-            %if spannedJoint.getCoordinateSet().getSize()~=0
-            if spannedJoint.numCoordinates()~=0
+        if getOpenSimVersion()<4.0
+            nr_coords = spannedJoint.getCoordinateSet().getSize();
+        else
+            nr_coords = spannedJoint.numCoordinates();
+        end
+            if nr_coords~=0
+
         jointNameSet{n_spanJoint} =  spannedJointName;
         n_spanJoint = n_spanJoint+1;
             else
                 NoDofjointNameSet{n_spanJointNoDof} =  spannedJointName;
                 n_spanJointNoDof = n_spanJointNoDof+1;
             end
         spannedJointNameOld = spannedJointName;
-        %body =  spannedJoint.getParentBody();
-        bodyName = jointStructure.(char(spannedJointName)).parentBody;
-        body = osimModel.getBodySet().get(bodyName);
-
+        if getOpenSimVersion()<4.0
+            body =  spannedJoint.getParentBody();
+        else
+            body =  spannedJoint.getParentFrame().findBaseFrame();
+        end
     end
-    bodyName = body.getName();
-
+    bodyName = char(body.getName());
 end
 
 if isempty(jointNameSet)

MATLAB_tool/MuscleParOptTool/private/sampleMuscleQuantities.m

@@ -56,10 +56,7 @@
 % limit (1) or not (0) the discretization of the joint space sampling
 limit_discr = 0; 
 % minimum angular discretization
-%min_increm_in_deg = 2.5;
-min_increm_in_deg = 1;
-
-
+min_increm_in_deg = 2.5;
 %=======================
 
 % initialize the model
@@ -68,7 +65,7 @@
 % getting the joint crossed by a muscle
 muscleCrossedJointSet = getJointsSpannedByMuscle(osimModel, OSMuscle);
 
-% index f+or effective dofs
+% index for effective dofs
 n_dof = 1;
 DOF_Index = [];
 for n_joint = 1:size(muscleCrossedJointSet,2)
@@ -77,12 +74,15 @@
     curr_joint = muscleCrossedJointSet{n_joint};
 
     % get CoordinateSet for the crossed joint
-    %curr_joint_CoordinateSet = osimModel.getJointSet().get(curr_joint).getCoordinateSet();
+    if getOpenSimVersion()<4.0
+        curr_joint_CoordinateSet = osimModel.getJointSet().get(curr_joint).getCoordinateSet();
+        % Initial estimation of the nr of Dof of the CoordinateSet for that
+        % joint before checking for locked and constraint dofs.
+        nDOF = osimModel.getJointSet().get(curr_joint).getCoordinateSet().getSize();
+    else
+        nDOF = osimModel.getJointSet().get(curr_joint).numCoordinates();
+    end
 
-    % Initial estimation of the nr of Dof of the CoordinateSet for that
-    % joint before checking for locked and constraint dofs.
-    %nDOF = osimModel.getJointSet().get(curr_joint).getCoordinateSet().getSize();
-    nDOF = osimModel.getJointSet().get(curr_joint).numCoordinates();
     % skip welded joint and removes welded joint from muscleCrossedJointSet
     if nDOF == 0;
         continue;
@@ -93,9 +93,11 @@
     for n_coord = 0:nDOF-1
 
         % get coordinate
-        %curr_coord = curr_joint_CoordinateSet.get(n_coord);
-        curr_coord = osimModel.getJointSet().get(curr_joint).get_coordinates(n_coord);
-
+        if getOpenSimVersion()<4.0
+            curr_coord = curr_joint_CoordinateSet.get(n_coord);
+        else
+            curr_coord = osimModel.getJointSet().get(curr_joint).get_coordinates(n_coord);
+        end
         curr_coord_name = char(curr_coord.getName());
 
         % skip dof if locked

README.md

@@ -1,3 +1,17 @@
+# Table of Contents <!-- omit in toc -->
+- [Introduction](#introduction)
+- [Requirements](#requirements)
+- [Contents from the paper](#contents-from-the-paper)
+  - [Algorithm description](#algorithm-description)
+  - [Considered Cases](#considered-cases)
+- [Versions of the tool](#versions-of-the-tool)
+  - [MATLAB version](#matlab-version)
+  - [Python tool (currently under revision)](#python-tool-currently-under-revision)
+  - [OpenSim C++ plugin and User Interface Menu](#opensim-c-plugin-and-user-interface-menu)
+- [Notes on upgrading to OpenSim 4.0 and backward compatibility](#notes-on-upgrading-to-opensim-40-and-backward-compatibility)
+- [Contributors](#contributors)
+- [References](#references)
+
 # Introduction
 
 This repository contains a MATLAB package implementing an algorithm for optimizing the parameters of Hill-type muscle models defined by adimensional force-lenght-velocity curves, as described by _Zajac (1989)_.
@@ -56,7 +70,7 @@ __Please note__ that reproducing the sensitivity study can be time-consuming, de
 
 ## Python tool (currently under revision)
 
-A Python version of the tool has been written and kindly shared by @eravera. It is currently **under assessment** although tested and working.
+A Python version of the tool has been written and kindly shared by @eravera.
 
 
 ## OpenSim C++ plugin and User Interface Menu
@@ -66,12 +80,22 @@ A generic tool to optimize musculotendon parameters in musculoskeletal models is
 * as menu extension of the OpenSim GUI (graphical user interface).
 Please refer directly to the repository and to the nice documentation available at [this website](http://muscleoptimizer.github.io/MuscleOptimizer/).
 
+# Notes on upgrading to OpenSim 4.0 and backward compatibility
+
+The tool can be used with both OpenSim 3.3 and 4.0 without modifications.
+
+Please be aware that the results provided by the two versions are slightly different, as noticed when running our tests.
+
+To investigate further please consult the tests folder and scripts.
+
+
 # Contributors
 
 Special thanks to:
 
 * **Bryce Killen** from KU Leuven, Belgium, for updating the code for OpenSim 4.1!
 * **Emiliano Ravera** from Instituto de Investigación y Desarrollo en Bioingenieria y Bioinformática, IBB (CONICET-UNER), Argentina.
+* **Dario Sciacca** from KU Leuven, Belgium, for his feedback on version 4.x.
 
 # References
 * Zajac, F.E. Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control. Critical Reviews in Biomedical Engineering. 17: 359-411, 1989. [LINK](https://www.ncbi.nlm.nih.gov/pubmed/2676342)

tests/TESTS_MusOptTool.m

@@ -36,8 +36,15 @@
 % optimize
 [osimModel_opt, SimsInfo] = optimMuscleParams(reference_model, target_model, N_eval, log_folder);
 
-% test
-expected_model  = './models/test2_expected.osim';
+% test: note that from OpenSim 3.3 to 4.x results are slightly different
+% if compared to 'test2_expected.osim', a model optimized with OpenSim 4
+% would have a loss function ~ 0.001.
+if getOpenSimVersion()<4
+    expected_model  = './models/test2_expected.osim';
+else
+    expected_model  = './models/test2_expected_v4.osim';
+end
+
 osimExp = Model(expected_model);
 % load muscles
 exp_mus = osimExp.getMuscles();