Added the option Fischer2019_pelvis incl. example and submodule.

.gitmodules

@@ -0,0 +1,3 @@
+[submodule "STAPLE/external/PelvicLandmarkIdentification"]
+	path = STAPLE/external/PelvicLandmarkIdentification
+	url = https://github.com/RWTHmediTEC/PelvicLandmarkIdentification

Example_use_Fischer2019.m

@@ -0,0 +1,109 @@
+%-------------------------------------------------------------------------%
+% Copyright (c) 2024 MCM Fischer.                                         %
+%    Author:   Maximilian C. M. Fischer,  2024                            %
+%    email:    N/A                                                        %
+% ----------------------------------------------------------------------- %
+% This example demonstrates how to setup a STAPLE workflow using an       %
+% algorithm for the morphological analyses of the pelvis published by     %
+% Fischer et al. https://doi.org/10.1038/s41598-019-49573-4               %
+% ----------------------------------------------------------------------- %
+clear; clc; close all
+addpath(genpath('./STAPLE'));
+
+%----------%
+% SETTINGS %
+%----------%
+output_models_folder = 'opensim_models_examples';
+
+% folder where the various datasets (and their geometries) are located.
+datasets_folder = 'bone_datasets';
+
+% datasets that you would like to process
+dataset_set = {'JIA_MRI'};
+
+% anthropometrics
+% H = ?; %cm
+subj_mass = 76.5; %kg
+
+% Cell array with the bone geometries that you would like to process
+% NOTE: The complete pelvis mesh including both hip bones and the sacrum is
+% required for the algorithm of [Fischer 2019]. Consider the requirements 
+% for the mesh of the pelvis described in the document header of the 
+% fuction pelvicLandmarkID.
+bones_list = {'pelvis','femur_r'};
+
+% visualization geometry format (options: 'stl' or 'obj')
+vis_geom_format = 'obj';
+
+% choose the definition of the joint coordinate systems (see documentation)
+% options: 'Modenese2018' or 'auto2020'
+workflow = 'auto2020';
+%--------------------------------------
+
+tic
+
+% create model folder if required
+if ~isfolder(output_models_folder); mkdir(output_models_folder); end
+
+for n_d = 1:numel(dataset_set)
+
+    % current dataset being processed
+    cur_dataset = dataset_set{n_d};
+
+    % folder from which triangulations will be read
+    tri_folder = fullfile(datasets_folder, cur_dataset,'tri');
+
+    % create geometry set structure for all 3D bone geometries in the dataset
+    triGeom_set = createTriGeomSet(bones_list, tri_folder);
+
+    % get the body side (can also be specified by user as input to funcs)
+    side = inferBodySideFromAnatomicStruct(triGeom_set);
+
+    % model and model file naming
+    model_name = ['auto_',dataset_set{n_d},'_',upper(side)];
+    model_file_name = [model_name, '.osim'];
+
+    % log printout
+    log_file = fullfile(output_models_folder, [model_name, '.log']);
+    logConsolePrintout('on', log_file);
+
+    % create bone geometry folder for visualization
+    geometry_folder_name = [model_name, '_Geometry'];
+    geometry_folder_path = fullfile(output_models_folder,geometry_folder_name);
+
+    % convert geometries in chosen format (30% of faces for faster visualization)
+    writeModelGeometriesFolder(triGeom_set, geometry_folder_path, vis_geom_format,0.3);
+
+    % initialize OpenSim model
+    osimModel = initializeOpenSimModel(model_name);
+
+    % create bodies
+    osimModel = addBodiesFromTriGeomBoneSet(osimModel, triGeom_set, geometry_folder_name, vis_geom_format);
+
+    % process bone geometries (compute joint parameters and identify markers)
+    %                                                         pelvis       femur     tibia 
+    [JCS, BL, CS] = processTriGeomBoneSet(triGeom_set, side, 'Fischer2019','Kai2014','Kai2014');
+
+    % create joints
+    createOpenSimModelJoints(osimModel, JCS, workflow);
+
+    % update mass properties to those estimated using a scale version of
+    % gait2392 with COM based on Winters's book.
+    osimModel = assignMassPropsToSegments(osimModel, JCS, subj_mass);
+
+    % add markers to the bones
+    addBoneLandmarksAsMarkers(osimModel, BL);
+
+    % finalize connections
+    osimModel.finalizeConnections();
+
+    % print
+    osimModel.print(fullfile(output_models_folder, model_file_name));
+
+    % inform the user about time employed to create the model
+    disp('-------------------------')
+    disp(['Model generated in ', num2str(toc)]);
+    disp(['Saved as ', fullfile(output_models_folder, model_file_name),'.']);
+    disp('-------------------------')
+    logConsolePrintout('off');
+end

README.md

@@ -1,6 +1,8 @@
 # STAPLE: Shared Tools for Automatic Personalised Lower Extremity modelling <!-- omit in toc -->
 
-[![License: CC BY-NC 4.0](https://img.shields.io/badge/License-CC%20BY--NC%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc/4.0/)  [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4275285.svg)](https://doi.org/10.5281/zenodo.4275285)   ![visitors](https://visitor-badge.glitch.me/badge?page_id=modenaxe.msk-STAPLE)<!-- omit in toc -->
+[![License: CC BY-NC 4.0](https://img.shields.io/badge/License-CC%20BY--NC%204.0-lightgrey.svg)](https://creativecommons.org/licenses/by-nc/4.0/)  [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.4428103.svg)](https://doi.org/10.5281/zenodo.4428103)   ![visitors](https://visitor-badge.laobi.icu/badge?page_id=modenaxe.msk-STAPLE)<!-- omit in toc -->
+
+### ⚠️ STAPLE is released under a [non-commercial license](#license) and it is free to use for academic purposes only. For any other use please contact the authors. 
 
 # Table of contents <!-- omit in toc -->
 - [What is STAPLE?](#what-is-staple)
@@ -14,7 +16,7 @@
   - [Detailed steps to setup a STAPLE workflow](#detailed-steps-to-setup-a-staple-workflow)
   - [Algorithms for bone morphological analysis](#algorithms-for-bone-morphological-analysis)
   - [Datasets available for testing](#datasets-available-for-testing)
-  - [Further notes on STAPLE](#further-notes-on-staple)
+  - [STAPLE variables and conventions](#STAPLE-variables-and-conventions)
 - [Troubleshooting your workflow](#troubleshooting-your-workflow)
 - [Does STAPLE work only with OpenSim?](#does-staple-work-only-with-opensim)
 - [What are the differences between the STAPLE toolbox, NMSBuilder and the MAP Client for generating OpenSim models?](#what-are-the-differences-between-the-staple-toolbox-nmsbuilder-and-the-map-client-for-generating-opensim-models)
@@ -62,6 +64,10 @@ For example, models of hip, knee and ankle joints can be created as individual m
 
 ![articular_surfaces](./images/artic_surfaces.png)
 
+* **Transformations to local reference systems**: STAPLE performs morphological analyses using the medical images reference system where the bone geometries are segmented, but then it can transform both the bone geometries and model objects to standard segment-based biomechanical reference systems. This is an important feature if you are using STAPLE in OpenSim, as joint reactions are expressed in body reference systems. See the provided [example of this functionality](./Example_BCS_leg_right.m).
+
+![reference_systems](./images/reference_systems.png)
+
 * **Basic identification of bony landmarks**: certain bony landmarks can be easily identified following the morphological analysis of the bone surfaces. These landmarks are intended as first guess for registration with gait analysis data.
 
 * **Merging subject-specific and generic models**: STAPLE includes some basic utilities to merge partial subject-specific skeletal models (obtainable from localized MRI scans) with generic musculoskeletal models. See the provided [advanced examples](./advanced_examples).
@@ -286,11 +292,11 @@ The BCS is a MATLAB structure with the following fields:
 When running a morphological analysis, the joint coordinate system (JCS) associated with the analysed geometry is also returned, together with the BCS. 
 
 The JCS is a MATLAB structure with the following fields:
-1. **<joint name>**: 
+1. **joint name**: the name of the joint if more than one are stored in the variable. 
 2. **Origin** [3x1] vector: the origin of the joint in the bone being analysed
 3. **V** [3x3] matrix: the pose matrix of the joint
 4. **parent_location** [1x3] vector: the origin of the parent coordinate system in OpenSim
-5  **parent_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order.
+5. **parent_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order.
 6. **child_location** [1x3] vector: the origin of the child coordinate system in OpenSim
 7. **child_orientation** [1x3] vector: three angles describing the pose of the parent reference system in OpenSim, obtained by matrix decomposition with XYZ rotation order.
 
@@ -315,6 +321,14 @@ Note that the algorithms used on a single bone normally cannot define a joint co
 
 Bone landmarks are stored in MATLAB structures as [3x1] vectors. Each field name of the structure corresponds to the name of a landmark.
 
+|  Bone  geometry       | Landmarks              	| 
+| ---                   | ---                     	| 
+| pelvis 				| <ul><li>RASI/LASI: anterior superior iliac spine</li><li>RPSI/LPSI: posterior superior iliac spine</li><li>SYMP: pubic symphysis</li></ul> 			| 
+| femur 				| <ul><li>RKNE/LKNE: lateral femoral epicondyle</li><li>RMFC/LMFC: medial femoral epicondyle</li><li>RTRO/LTRO: greater trochanter</li></ul> 			| 
+| tibia 				| <ul><li>RTTB/LTTB: tibial tuberosity</li><li>RHFB/LHFB: head of fibula</li><li>RANK/LANK: lateral malleolus</li><li>RMMA/LMMA: medial malleolus</li><li>RLM/LLM: most distal point of lateral malleolus</li></ul>|
+| patella 				| RLOW/LLOW: most distal point of patella |
+| talus 				| None |
+| foot 					| <ul><li>RHEE/LHEE: heel</li><li>RD5M/LD5M: distal point of 5th metatarsal bone</li><li>RD5MPROX/LD5MPROX: proximal point of 5th metatarsal bone</li><li>RD1M/LD1M: distal point of first metatarsal bone</li><li>R1MGROUND/L1MGROUND: distal point of first metatarsal bone on foot sole</li><li>R5MGROUND/L5MGROUND: distal point of 5th metatarsal bone on foot sole</li><li>RHEEGROUND/LHEEGROUND: calcaneus most distal point on foot sole</li></ul>|
 
 ## Troubleshooting your workflow
 
@@ -323,9 +337,11 @@ Before informing us as suggested in [the contributing guidelines](/CONTRIBUTING.
 - [ ] ensure that the entire `STAPLE` folder is on your MATLAB path
 - [ ] if you are using one of the workflows from the examples, ensure that you are using the correct names for the bone geometries, e.g. `pelvis_no_sacrum`, `femur_r`, `tibia_r`, etc.
 - [ ] ensure that the quality of your bone surface geometries is sufficient for running your selected algorithm. The GIBOC algorithms, in particular, require relatively good quality surface meshes. You can have an idea of what we mean by "good mesh" consulting the table that describes the datasets provided with STAPLE. If necessary, use the filters and tools available on software like [MeshLab](https://www.meshlab.net/) to improve your dataset.
+- [ ] processing the meshes so that they are triangular and contain a reasonable number of vertices (decimation can be applied) also helps with the OpenSim visualization. See [this related tweet](https://twitter.com/JohnRHutchinson/status/1455490276298997766).
 - [ ] if possible, try running alternative algorithms in order to establish if your bone mesh is processable in the first place or if you have encounter a proper bug in the algorithm. You can specify the processing algorithms for each bone as input to the `processTriGeomBoneSet.m` function. For bad quality meshes, we recommend using the `STAPLE` algorithm at the pelvis and `Kai2014` algorithms for femur and tibia. Keep an eye on issue #78 as we will develop an example on how to do this.
 - [ ] verify that processing of your dataset is not failing because of the [current limitations](#current-limitations) of the STAPLE toolbox.
 
+
 ## Does STAPLE work only with OpenSim?
 
 The algorithms collected in the STAPLE toolbox were proposed in publications that did not have modelling focus, and can be applied in broader contexts, such as reference system definition for orthopaedics applications or modelling in other platforms. Similarly, the outputs of a STAPLE workflow, e.g. from `processTriGeomBoneSet.m`, include all the necessary information to create a kinematic or kinetic model in any biomechanical modelling workflow. All our modelling examples, however, rely on OpenSim.

STAPLE/GIBOC-core/SubFunctions/FittingFun/LSSLFitPatellaRidge.m

@@ -1,72 +1,86 @@
-function [ U, Uridge ,LowestPoints_End ] = LSSLFitPatellaRidge( TR,U,nbSlice,StartDist,EndDist, debug_plot)
 %Least Square Straight Line Fit on Patellar Ridge
 %To obtain the direction of the ridge a Least Square Straight Line is
 %fitted on the lowest points on slices of normal U, and the normal U is
 %updated until convergence (of U or number of iterations > 100)
+% U is the vector ~aligned with the ridge
+%-------------------------------------------------------------------------%
+%  Author:   Luca Modenese & Jean-Baptiste Renault. 
+%  Copyright 2020 Luca Modenese & Jean-Baptiste Renault
+%-------------------------------------------------------------------------%
+function [U, Uridge, LowestPoints] = LSSLFitPatellaRidge(   patellaTri,...
+                                                            U,...
+                                                            nbSlice,...
+                                                            startDist,...
+                                                            endDist,...
+                                                            debug_plots)
 
+% initialization
+if nargin<3; nbSlice = 50;                                 end
+if nargin<4; startDist = 0.025*range(patellaTri.Points*U); end %Offset distance from first points
+if nargin<5; endDist = startDist;                          end%Offset distance from last points
+if nargin<6; debug_plots=0;                                 end
 
-%% inputs Tests
-if nargin<3
-    StartDist = 0.025*range(TR.Points*U); %Offset distance from first points
-    EndDist = StartDist; %Offset distance from last points
-    nbSlice = 50;
-end
-
-if nargin<4
-    StartDist = 0.025*range(TR.Points*U); %Offset distance from first points
-    EndDist = StartDist; %Offset distance from last points
-end
-
-if nargin<6
-    debug_plot=0;
-end
-%% Code
+% initialize fitting
 U0 = U;
-reslts = [0;0;0];
-k=0;
+U_old = [0;0;0];
+k=0; 
 test = true;
 while test && k < 100 % Cut at 100 iterations (convergence generally occuring before 10)
     k=k+1;
-    Alt = linspace( min(TR.Points*U)+StartDist , max(TR.Points*U)-EndDist, nbSlice);
+    % points along U where the slices will be performed
+    Alt = linspace( min(patellaTri.Points*U)+startDist , max(patellaTri.Points*U)-endDist, nbSlice);
     LowestPoints = zeros(length(Alt),3);
+    % slice patella along U
     i=0;
     for d = -Alt
         i=i+1;
-        [ Curves , ~ , ~ ] = TriPlanIntersect( TR, U , d );
+        [ Curves , ~ , ~ ] = TriPlanIntersect( patellaTri, U , d );
         EdgePts = vertcat(Curves(:).Pts);
         [~,lowestPointID] = min(EdgePts(:,3));
         LowestPoints(i,:) = EdgePts(lowestPointID(1),:);
+        if debug_plots == 1
+            title('Slicing and LowerPoints on AS')
+            quickPlotTriang(patellaTri);
+            plot3(Curves(:).Pts(:,1), Curves(:).Pts(:,2), Curves(:).Pts(:,3),'k'); hold on
+            plot3(LowestPoints(:,1), LowestPoints(:,2), LowestPoints(:,3),'r'); hold on
+            axis equal
+        end
     end
 
+    % updating direction of U usign the variance of the patellar ridge path
     [Vridge_all,~] = eig(cov(LowestPoints));
-    U = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3);
 
-    U(3) = 0; U = normalizeV( U );
-    reslts(1:3,end+1) = U;
-    test = reslts(1:3,end-1)'*reslts(1:3,end) > (1 - 100*eps);
-end
-
-
-if nargout == 2
-Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3);
-
-elseif nargout == 3
-        Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3);
-        LowestPoints_End = LowestPoints;
+    % vector describing the direction of the ridge
+    Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3);
+
+    % the vector will be in the XY plane where the patella is a "circle"
+    % third component = 0
+%     U = Uridge;
+%     U(3) = 0; 
+    U = normalizeV([Uridge(1:2); 0]);
+    
+    % how much the vector is changing in direction
+    U_old(1:3,end+1) = U;
+    test = U_old(1:3,end-1)'*U_old(1:3,end) > (1 - 100*eps);
 end
 
-%% Plots
-if debug_plot == 1
+% debug plots
+if debug_plots == 1
     figure()
-    trisurf(TR,'facecolor','red','faceAlpha',0.5)
-    hold on
-    grid off
-    axis equal
+    title('Final vector of ridge')
+    trisurf(patellaTri,'facecolor','red','faceAlpha',0.5)
+    hold on;    grid off;     axis equal;
+    LocalCS        = struct();
+    LocalCS.Origin = [0 0 0]';
+    LocalCS.V      = eye(3);
+    quickPlotRefSystem(LocalCS);
+
+    % plot trace of points on the proximal surface
     pl3t(LowestPoints,'gs')
     plotArrow( U, 1.5, mean(LowestPoints), 30, 1, 'k')
-    Uridge = sign(U0'*Vridge_all(:,3))*Vridge_all(:,3);
     plotArrow( Uridge, 1.5, mean(LowestPoints), 50, 1, 'c')
     hold off
 end
+
 end