Merge pull request #323 from rcnl-org/dev

Merge up v1.2

CHANGELOG.md

@@ -1,5 +1,28 @@
 # Changelog
 
+## v1.2.0 - 2024-05-09
+
+### Added
+- MTP can now use translational coordinates spanned by muscles.
+- Added a plotting function to show the distribution of spring stiffnesses in GCP.
+
+### Fixed
+- Fixed a bug where Treatment Optimization tracking terms were using the wrong time array
+- MTP now properly saves passive model moments when input passive moment data has multiple columns with non-zero data.
+- Fixed a bug with plotting JMP results where the marker names ordering was not respected.
+- Fixed a bug with GCP damping terms not being applied correctly.
+- Fixed a bug where the wrong muscle tendon length file was parsed in Treatment Optimization.
+- Fixed a bug with the incorrect time array being used for the initial guess and dependency finding steps of Treatment Optimization.
+
+### Changed
+- Treatment Optimization now no longer splines results back to the initial time points, instead using the collocation time points as the final results. This reduces numerical inconsistencies between the tracking and verification steps.
+- Some Treatment Optimization plotting functions have been improved.
+- The surrogate model has been sped up.
+- The surrogate model now uses an independent data directory. A script has been added (`src/SurrogateModelCreation/surrogateKinematicsScript.m`) to generate Latin hypercube sample (LHS) kinematics for the surrogate model. This improves accuracy when finding novel motion.
+- .sto file parsing has been sped up.
+- Improved the quality of the Treatment Optimization initial guess.
+- Changed the way Treatment Optimization results are saved to make VO easier.
+
 ## v1.1.0 - 2024-03-10
 
 ### Added

src/GroundContactPersonalization/Analysis/plotGcpStiffnessCoefficients.m

@@ -0,0 +1,82 @@
+% This function is part of the NMSM Pipeline, see file for full license.
+%
+% (struct, struct) -> (None)
+% Plot GCP stiffness values from model and osimx files. 
+
+% ----------------------------------------------------------------------- %
+% The NMSM Pipeline is a toolkit for model personalization and treatment  %
+% optimization of neuromusculoskeletal models through OpenSim. See        %
+% nmsm.rice.edu and the NOTICE file for more information. The             %
+% NMSM Pipeline is developed at Rice University and supported by the US   %
+% National Institutes of Health (R01 EB030520).                           %
+%                                                                         %
+% Copyright (c) 2021 Rice University and the Authors                      %
+% Author(s): Spencer Williams                                             %
+%                                                                         %
+% Licensed under the Apache License, Version 2.0 (the "License");         %
+% you may not use this file except in compliance with the License.        %
+% You may obtain a copy of the License at                                 %
+% http://www.apache.org/licenses/LICENSE-2.0.                             %
+%                                                                         %
+% Unless required by applicable law or agreed to in writing, software     %
+% distributed under the License is distributed on an "AS IS" BASIS,       %
+% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or         %
+% implied. See the License for the specific language governing            %
+% permissions and limitations under the License.                          %
+% ----------------------------------------------------------------------- %
+
+function plotGcpStiffnessCoefficients(modelFileName, osimxFileName, ...
+    surfaceNumber)
+% Parse inputs
+[model, state] = Model(modelFileName);
+osimx = parseOsimxFile(osimxFileName, model);
+if nargin < 3
+    surfaceNumber = 1;
+end
+contactSurface = getFieldByNameOrError(osimx, 'contactSurface');
+assert(iscell(contactSurface), "The given .osimx file does not " + ...
+    "contain contact surfaces.")
+assert(surfaceNumber <= length(contactSurface), "Contact surface " + ...
+    "number is out of range.")
+contactSurface = contactSurface{1};
+
+% Get toes body name
+joints = getBodyJointNames(model, contactSurface.hindfootBodyName);
+assert(length(joints) == 2, "GCP supports two segment foot models only");
+for i = 1 : length(joints)
+    [parent, ~] = getJointBodyNames(model, joints(i));
+    if strcmp(parent, contactSurface.hindfootBodyName)
+        toesJointName = joints(i);
+        break
+    end
+end
+assert(model.getJointSet().get(toesJointName) ...
+    .numCoordinates() == 1, "GCP toes joint must be a pin joint");
+[~, toesBodyName] = getJointBodyNames(model, toesJointName);
+
+% Get spring locations in common reference frame
+calcnToToes = model.getBodySet().get(toesBodyName) ...
+    .findTransformBetween(state, model.getBodySet() ...
+    .get(contactSurface.hindfootBodyName));
+xOffset = calcnToToes.T.get(0);
+zOffset = calcnToToes.T.get(2);
+springX = zeros(1, length(contactSurface.springs));
+springZ = zeros(1, length(contactSurface.springs));
+stiffness = zeros(1, length(contactSurface.springs));
+for i = 1 : length(contactSurface.springs)
+    position = contactSurface.springs{i}.location;
+    offset = strcmp(contactSurface.springs{i}.parentBody, toesBodyName);
+    springX(i) = position(1) + xOffset * offset;
+    springZ(i) = position(3) + zOffset * offset;
+    stiffness(i) = contactSurface.springs{i}.springConstant;
+end
+
+% Plot values
+scatter(springZ, springX, 200, stiffness, "filled")
+set(gca, 'DataAspectRatio', [1, 1, 1])
+title("Stiffness coefficients")
+xlabel("Z location on foot (m)")
+ylabel("X location on foot (m)")
+colormap jet
+colorbar
+end

src/GroundContactPersonalization/Optimizations/ModelCalculation/calcModeledVerticalGroundReactionForce.m

@@ -58,7 +58,7 @@
     constant = -s .* (v .* ymax - c .* log(cosh((ymax + h) ./ c)));
     freglyVerticalGrf = -s .* (v .* height - c .* ...
         log(cosh((height + h) ./ c))) - constant;
-    springForces(2, i) = freglyVerticalGrf * (1 + dampingFactor * ...
+    springForces(2, i) = freglyVerticalGrf * (1 - dampingFactor * ...
         verticalVelocity);
     modeledVerticalGrf = modeledVerticalGrf + springForces(2, i);
 end

src/GroundContactPersonalization/parseGroundContactPersonalizationSettingsTree.m

@@ -164,7 +164,7 @@
 function task = getFootData(tree)
 task.isLeftFoot = strcmpi('true', ...
     getFieldByNameOrError(tree, 'is_left_foot').Text);
-hindfootBodyName = getFieldByName(tree, "hindfoot_body");
+hindfootBodyName = getFieldByName(tree, 'hindfoot_body');
 if ~isstruct(hindfootBodyName)
     throw(MException('', "<toes_coordinate> is replaced by <hindfoot_body> in the GCP settings file."))
 else

src/JointModelPersonalization/Analysis/reportDistanceErrorByMarker.m

@@ -26,7 +26,7 @@ function reportDistanceErrorByMarker(model, markerFileName, ...
 names = ArrayStr();
 names.append('time');
 for i=1:length(markerNames)
-    names.append(markerNames{i});
+    names.append(ikSolver.getMarkerNameForIndex(i-1));
 end
 storage.setColumnLabels(names);
 

src/JointModelPersonalization/parseJointModelPersonalizationSettingsTree.m

@@ -80,9 +80,9 @@
 end
 
 function output = getTask(model, tree)
-output.markerFile = parseElementTextByName(tree, "marker_file_name");
+output.markerFile = parseElementTextByName(tree, 'marker_file_name');
 output.anatomicalMarkers = getBooleanLogicFromField( ...
-    getFieldByName(tree, "anatomical_markers"));
+    getFieldByName(tree, 'anatomical_markers'));
 timeRange = getFieldByName(tree, 'time_range');
 if(isstruct(timeRange))
     timeRange = strsplit(timeRange.Text, ' ');
@@ -91,24 +91,24 @@
 end
 
 try
-    markerNames = parseSpaceSeparatedList(tree, "marker_names");
+    markerNames = parseSpaceSeparatedList(tree, 'marker_names');
     if ~isempty(markerNames)
         output.markerNames = markerNames;
     end
 catch; end
 try
-    freeMarkers = parseSpaceSeparatedList(tree, "free_markers");
+    freeMarkers = parseSpaceSeparatedList(tree, 'free_markers');
 catch
     freeMarkers = [];
 end
 
 output.parameters = {};
-if(isstruct(getFieldByName(tree, "JMPJointSet")))
+if(isstruct(getFieldByName(tree, 'JMPJointSet')))
     output.parameters = getJointParameters(tree.JMPJointSet);
 end
 output.scaling = [];
 output.markers = [];
-if isstruct(getFieldByName(tree, "JMPBodySet")) || ~isempty(freeMarkers)
+if isstruct(getFieldByName(tree, 'JMPBodySet')) || ~isempty(freeMarkers)
     [output.scaling, output.markers] = ...
         getBodyParameters(tree.JMPBodySet, model, freeMarkers);
 end
@@ -134,7 +134,7 @@
 % solve this problem, it's fine
 function inputs = getJointParameters(jointSetTree)
 inputs = {};
-if isfield(jointSetTree, "JMPJoint")
+if isfield(jointSetTree, 'JMPJoint')
     jointTree = jointSetTree.JMPJoint;
     counter = 1; % for index of parameter in output
     for i=1:length(jointTree)
@@ -186,7 +186,7 @@
 
 function [scaling, markers] = getBodyParameters( ...
     bodySetTree, model, freeMarkers)
-if isfield(bodySetTree, "JMPBody")
+if isfield(bodySetTree, 'JMPBody')
     bodyTree = bodySetTree.JMPBody;
     scaling = getScalingBodies(bodyTree);
     markers = getMarkers(bodyTree, model);
@@ -207,7 +207,7 @@
     end
     bodyName = body.Attributes.name;
     scaleBodies = strcmp(getFieldByNameOrError(body, ...
-        "scale_body").Text, "true");
+        'scale_body').Text, "true");
     if(scaleBodies)
         inputs(end + 1) = bodyName;
     end
@@ -224,7 +224,7 @@
         body = bodyTree{i};
     end
     bodyName = body.Attributes.name;
-    axesStrings = parseSpaceSeparatedList(body, "move_markers");
+    axesStrings = parseSpaceSeparatedList(body, 'move_markers');
     axes = zeros(1, 3);
     for j = 1:3
         axes(j) = axesStrings(j) == "true";
@@ -303,15 +303,15 @@
 
 function output = getParams(tree)
 output = struct();
-paramArgs = ["accuracy", "diff_min_change", "optimality_tolerance", ...
-    "function_tolerance", "step_tolerance", "max_function_evaluations"];
+paramArgs = {'accuracy', 'diff_min_change', 'optimality_tolerance', ...
+    'function_tolerance', 'step_tolerance', 'max_function_evaluations'};
 % name in matlab is different, use for output struct arg name
-paramName = ["accuracy", "diffMinChange", "optimalityTolerance", ...
-    "functionTolerance", "stepTolerance", "maxFunctionEvaluations"];
+paramName = {'accuracy', 'diffMinChange', 'optimalityTolerance', ...
+    'functionTolerance', 'stepTolerance', 'maxFunctionEvaluations'};
 for i=1:length(paramArgs)
-    value = getFieldByName(tree, paramArgs(i));
+    value = getFieldByName(tree, paramArgs{i});
     if(isstruct(value))
-        output.(paramName(i)) = str2double(value.Text);
+        output.(paramName{i}) = str2double(value.Text);
     end
 end
 end

src/MuscleTendonLengthInitialization/parseMuscleTendonLengthInitializationSettingsTree.m

@@ -31,19 +31,19 @@
 function inputs = parseMuscleTendonLengthInitializationSettingsTree( ...
     settingsTree)
 if strcmp(getFieldByNameOrError(settingsTree, ...
-        "MuscleTendonLengthInitialization").is_enabled.Text, "true")
+        'MuscleTendonLengthInitialization').is_enabled.Text, "true")
     inputs = getInputs(settingsTree);
     inputs = getMtpModelInputs(inputs);
     inputs = getMuscleVolume(inputs);
-    inputs = rmfield(inputs, "model");
+    inputs = rmfield(inputs, 'model');
 else
     inputs = false;
 end
 end
 
 function inputs = getInputs(tree)
 inputs = parseMtpNcpSharedInputs(tree);
-inputs = getPassiveData(getFieldByNameOrError(tree, "MuscleTendonLengthInitialization"), inputs);
+inputs = getPassiveData(getFieldByNameOrError(tree, 'MuscleTendonLengthInitialization'), inputs);
 inputs = getTask(tree, inputs);
 
 inputs = getNormalizedFiberLengthSettings(tree, inputs);
@@ -146,7 +146,7 @@
 
 function inputs = getNormalizedFiberLengthSettings(tree, inputs)
 normalizedFiberLengthGroupNames = parseSpaceSeparatedList(tree, ...
-    "normalized_fiber_length_muscle_groups");
+    'normalized_fiber_length_muscle_groups');
 inputs.normalizedFiberLengthGroups = groupNamesToGroups( ...
     normalizedFiberLengthGroupNames, inputs.model);
 inputs.numMuscleGroups = numel(inputs.normalizedFiberLengthGroups);

src/MuscleTendonPersonalization/Analysis/extractMtpDataFromSto.m

@@ -44,7 +44,7 @@
     end
     data = cell(1, numel(dataFiles));
     for i = 1:numel(dataFiles)
-        dataStorage = Storage(fullfile(fullfile(resultsDirectory, dataFiles{i})));
+        dataStorage = Storage(fullfile(resultsDirectory, dataFiles{i}));
         columnNames = getStorageColumnNames(dataStorage);
         data{i} = storageToDoubleMatrix(dataStorage)';
     end

src/MuscleTendonPersonalization/Analysis/plotMtpPassiveMomentCurves.m

@@ -29,50 +29,80 @@
 % permissions and limitations under the License.                          %
 % ----------------------------------------------------------------------- %
 
-function plotMtpPassiveMomentCurves(resultsDirectory)
+function plotMtpPassiveMomentCurves(resultsDirectory, figureWidth, ...
+    figureHeight)
 analysisDirectory = fullfile(resultsDirectory, "Analysis");
 [momentNames, passiveMomentsExperimental] = extractMtpDataFromSto( ...
     fullfile(analysisDirectory, "passiveJointMomentsExperimental"));
 [~, passiveMomentsModel] = extractMtpDataFromSto( ...
     fullfile(analysisDirectory, "passiveJointMomentsModeled"));
-momentNames = strrep(momentNames, '_', ' ');
-meanMomentsExperimental = mean(passiveMomentsExperimental, 3);
-stdMomentsExperimental = std(passiveMomentsExperimental, [], 3);
-meanMomentsModel = mean(passiveMomentsModel, 3);
-stdMomentsModel = std(passiveMomentsModel, [], 3);
-maxMoment = max([max(meanMomentsExperimental, [], 'all'), ...
-    max(meanMomentsModel, [], 'all')]);
-minMoment = min([min(meanMomentsExperimental, [], 'all'), ...
-    min(meanMomentsModel, [], 'all')]);
-time = 1:1:size(meanMomentsModel,1);
-
-figureWidth = ceil(sqrt(numel(momentNames)));
-figureHeight = ceil(numel(momentNames)/figureWidth);
+columnsWithAllZeros = all(passiveMomentsExperimental == 0, 1);
+experimentalMomentFilesDirectory = dir(fullfile(analysisDirectory, ...
+    "passiveJointMomentsExperimental"));
+experimentalMomentFilesDirectory = experimentalMomentFilesDirectory(3:end);
+plotLabels = [];
+for i = 1 : size(passiveMomentsExperimental, 3)
+    trialPrefix = strrep(experimentalMomentFilesDirectory(i).name, ...
+        "passiveJointMomentsExperimental.sto", "");
+    plotLabels = [plotLabels, ...
+        strcat(trialPrefix, momentNames(~columnsWithAllZeros(:, :, i)))];
+end
+plotLabels = strrep(plotLabels, "_", " ");
+passiveMomentsModel = passiveMomentsModel(:, ~columnsWithAllZeros(1,:,:));
+passiveMomentsExperimental = ...
+    passiveMomentsExperimental(:, ~columnsWithAllZeros(1,:,:));
 
-figure(Name = strcat(resultsDirectory, " Passive Moment Curves"), ...
+if nargin < 2
+    figureWidth = ceil(sqrt(numel(plotLabels)));
+    figureHeight = ceil(numel(plotLabels)/figureWidth);
+elseif nargin < 3
+    figureHeight = ceil(sqrt(numel(plotLabels)));
+end
+figureSize = figureWidth * figureHeight;
+figure(Name = strcat(resultsDirectory, " Passive Moment Matching"), ...
     Units='normalized', ...
     Position=[0.05 0.05 0.9 0.85])
+subplotNumber = 1;
+figureNumber = 1;
 t = tiledlayout(figureHeight, figureWidth, ...
     TileSpacing='Compact', Padding='Compact');
+minMoment = min([ ...
+    min(passiveMomentsExperimental, [], "all"), ...
+    min(passiveMomentsModel, [], "all")]);
+maxMoment = max([ ...
+    max(passiveMomentsExperimental, [], "all"), ...
+    max(passiveMomentsModel, [], "all")]);
+xlabel(t, "Joint Position")
+ylabel(t, "Joint Moment [Nm]")
 
-for i = 1:numel(momentNames)
-    nexttile(i);
+for i = 1 : numel(plotLabels)
+    if i > figureSize * figureNumber
+        figureNumber = figureNumber + 1;
+        figure(Name="Treatment Optimization Joint Angles", ...
+            Units='normalized', ...
+            Position=[0.05 0.05 0.9 0.85])
+        t = tiledlayout(figureHeight, figureWidth, ...
+            TileSpacing='Compact', Padding='Compact');
+        xlabel(t, "Percent Movement [0-100%]")
+        ylabel(t, "Joint Angle [deg]")
+        subplotNumber = 1;
+    end
+    nexttile(subplotNumber)
     hold on
-    plotMeanAndStd(meanMomentsExperimental(:,i), stdMomentsExperimental(:,i), ...
-        time, 'k-')
-    plotMeanAndStd(meanMomentsModel(:,i), stdMomentsModel(:,i), time, 'r-')
+    plot(passiveMomentsExperimental(:, i), ...
+        LineWidth=3, ...
+        Color = "k")
+    plot(passiveMomentsModel(:, i), ...
+        LineWidth=3, ...
+        Color = "r")
     hold off
-    set(gca, fontsize=11)
-    rmse = rms(passiveMomentsExperimental(:,i) - passiveMomentsModel(:,i));
-    title(sprintf("%s \n RMSE: %.4f", ...
-        momentNames(i), rmse), FontSize=12)
-    axis([time(1) time(end) minMoment maxMoment])
-    if i == 1
-        legend("Experimental", "Model")
-    end
-    if mod(i,figureWidth) == 1
-        ylabel("Moment [Nm]")
+    title(plotLabels(i));
+    if subplotNumber == 1
+        legend("Experimental Moments", "Model Moments")
     end
+    xlim([1 size(passiveMomentsExperimental, 1)])
+    ylim([minMoment, maxMoment])
+    subplotNumber = subplotNumber + 1;
 end
 end