This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/README.md b/README.md
new file mode 100755
index 0000000..ac7f4e1
--- /dev/null
+++ b/README.md
@@ -0,0 +1,114 @@
+# EPFL upper limb musculoskeletal model
+This upper limb musculoskeletal model was developed at the EPFL by David Ingram and Ehsan Sarshari. It is implemented in Matlab. The details are explained in README files within each of the following subfolders.
+
+## Getting Started
+The following instructions will get you a copy of the project locally.
+* Clone the repository locally
+* Start matlab and set path one of the model's directorie (see below).
+
+## Prerequisites
+The model requires Matlab 2015 with the following toolboxes.
+
+## Description of the model
+
+ShoulderModel
+
+This folder contains the musculoskeletal model of the human shoulder and elbow that can replicate healthy subjects and patients with anatomical total shoulder arthroplasty (aTSA).
+
+
+## ShoulderModel_pre_post_processing
+This folder contains scripts and documentations that allow/explain transformation of the patients data obtained from the preoperative Blueprint software and CTs into inputs that are consistent and readable by different versions of the MSM model.
+
+
+## Authors
+* David Ingram (EPFL-LA)
+Development of the 1st version of the model, during his PhD thesis.
+https://infoscience.epfl.ch/record/204692
+
+* Ehsan Sarshari (EPFL-LA/LBO)
+Development of the 2nd version of the model, during his PhD thesis.
+https://infoscience.epfl.ch/record/256396
+
+## Contributors
+* Christoph Engelhardt (EPFL-LBO): Performed anatomical measurements on MRI of healthy volunteer.
+* Yasmine Boulanaache (EPFL-LBO): Performed anatomical measurements on MRI of healthy volunteer.
+* Matteo Mancuso (EPFL-LMAM): Performed kinematics and EMG measurements on healthy volunteer.
+* Fabio Becche (CHUV-RAD): set MRI protocols to build generic model from volunteer
+* Philippe Mullhaupt (EPFL-LA): PhDs' supervision
+
+## Bibliography
+These models are described in the following journal articles.
+
+
+* Sarshari, E., Mancuso, M., Terrier, A., Farron, A., Mullhaupt, P., Pioletti, D., 2020.
+Muscle co-contraction in an upper limb musculoskeletal model: EMG-assisted vs. standard load-sharing. Computer Methods in Biomechanics and Biomedical Engineering 1–14.
+* Sarshari, E., Mancuso, M., Terrier, A., Farron, A., Mullhaupt, P., Pioletti, D., 2020.
+Feasibility of an alternative method to estimate glenohumeral joint center from videogrammetry measurements and CT/MRI of patients. Computer Methods in Biomechanics and Biomedical Engineering 1–10.
+* Sarshari, E., Farron, A., Terrier, A., Pioletti, D., Mullhaupt, P.
+A simulation framework for humeral head translations
+(2017) Medical Engineering and Physics, 49, pp. 140-147.
+DOI: 10.1016/j.medengphy.2017.08.013
+
+* Engelhardt, C., Farron, A., Becce, F., Place, N., Pioletti, D.P., Terrier, A.
+Effects of glenoid inclination and acromion index on humeral head translation and glenoid articular cartilage strain
+(2017) Journal of Shoulder and Elbow Surgery, 26 (1), pp. 157-164.
+DOI: 10.1016/j.jse.2016.05.031
+
+* Engelhardt, C., Ingram, D., Müllhaupt, P., Farron, A., Becce, F., Pioletti, D., Terrier, A.
+Effect of partial-thickness tear on loading capacities of the supraspinatus tendon: a finite element analysis
+(2016) Computer Methods in Biomechanics and Biomedical Engineering, 19 (8), pp. 875-882.
+DOI: 10.1080/10255842.2015.1075012
+
+* Ingram, D., Engelhardt, C., Farron, A., Terrier, A., Müllhaupt, P.
+Modelling of the human shoulder as a parallel mechanism without constraints
+(2016) Mechanism and Machine Theory, 100, pp. 120-137.
+DOI: 10.1016/j.mechmachtheory.2016.02.004
+
+* Ingram, D., Engelhardt, C., Farron, A., Terrier, A., Müllhaupt, P.
+Improving anterior deltoid activity in a musculoskeletal shoulder model – an analysis of the torque-feasible space at the sternoclavicular joint
+(2016) Computer Methods in Biomechanics and Biomedical Engineering, 19 (4), pp. 450-463.
+DOI: 10.1080/10255842.2015.1042465
+
+* Engelhardt, C., Malfroy Camine, V., Ingram, D., Müllhaupt, P., Farron, A., Pioletti, D., Terrier, A.
+Comparison of an EMG-based and a stress-based method to predict shoulder muscle forces
+(2015) Computer Methods in Biomechanics and Biomedical Engineering, 18 (12), pp. 1272-1279.
+DOI: 10.1080/10255842.2014.899587
+
+* Ingram, D., Engelhardt, C., Farron, A., Terrier, A., Müllhaupt, P.
+Muscle moment-arms: a key element in muscle-force estimation
+(2015) Computer Methods in Biomechanics and Biomedical Engineering, 18 (5), pp. 506-513.
+DOI: 10.1080/10255842.2013.818666
+
+* Ingram, D., Engelhardt, C., Farron, A., Terrier, A., Mullhaupt, P.
+A minimal set of coordinates for describing humanoid shoulder motion
+(2013) IEEE International Conference on Intelligent Robots and Systems, art. no. 6697159, pp. 5537-5544.
+DOI: 10.1109/IROS.2013.6697159
+
+* Terrier, A., Aeberhard, M., Michellod, Y., Mullhaupt, P., Gillet, D., Farron, A., Pioletti, D.P.
+A musculoskeletal shoulder model based on pseudo-inverse and null-space optimization
+(2010) Medical Engineering and Physics, 32 (9), pp. 1050-1056.
+DOI: 10.1016/j.medengphy.2010.07.006
+
+* Aeberhard, M., Michellod, Y., Mullhaupt, P., Terrier, A., Pioletti, D.P., Gillet, D.
+Dynamical biomechanical model of the shoulder: Null space based optimization of the overactuated system
+(2008) 2008 IEEE International Conference on Robotics and Biomimetics, ROBIO 2008, art. no. 4912981, pp. 67-73.
+DOI: 10.1109/ROBIO.2009.4912981
+
+## License
+This should be discussed with Technology Transfer office (TTO) of EPFL
+
+## Acknowledgments
+The project was financially supported by the Swiss National Science Foundation (SNSF)
+* K-32K1_122512
+* CR32I2_143704
+* CR32I2_162766
+
+and by the Lausanne Orthopedic Research Foundation (LORF).
+
+We thank Julien Ston, former PhD student, the young healthy volunteer who spent hours in the MRI and in the Laboratory of Movement Analysis and Measurements (EPFL-LMAM).
+% Ribcage Ellipsoid Data for scapulothoracic constraints
+DYDATA.OE = REDATA.Centre;
+DYDATA.AE = [REDATA.TSaxes, REDATA.AIaxes];
+
+% Mass Parameters
+Mc = 0.156;
+Ms = 0.704;
+Mh = 2.052; % The body weight of our male subject is Mh*100/2.4 (=85.5 Kg) according to R. Dumas
+M_u = 0.9012; % not in the paper by Breteler but close to "A multibody biomechanical model of the upper limb including the shoulder girdle"
+Mr = 0.5523; % not in the paper by Breteler "
+
+DYDATA.MHand = 0;% The Radius mass can also contain a weight carried in the hand.
+DYDATA.LHand = 0;% The inertia will be MHand*LHand^2, the position from CP;
+
+% Gravitational Constant
+g = 9.81;
+
+% Inertia Parameters
+Ict = 0.001; Icl = 0.003; % from before, but I do not like them as they are not consistent with the other segments values or the portugease paper
+Ist = 0.007; Isl = 0.007; % the same as above
+Iht = 0.0165; Ihl = 0.0033; % these values are quite different that what it was before.
+Iut = 0.0060; Iul = 9.6437e-04; % defined using the SUBJECT_TOOL_Update_Weight by giving 85.5kg and 1.86 m, and the values are comparable to those of Portugease
+Irt = 0.0037; Irl = 5.9107e-04; % the same as above
+
+% ATTENTION!! One needs to apply the parallel axis theorem to be correct.
+% initialize glenoid fossa orientation for the subject specific toolbox [version inclination GC_SA]
+BLDATA.Glenoid_Orientations_Generic = [-4.28 9.77 25.0523 -6.4111 -8.4591 19.7584];% for the visualization we keep this to have the initial values unchanged to visualize
+% Ribcage Ellipsoid Data for scapulothoracic constraints
+DYDATA.OE = REDATA.Centre;
+DYDATA.AE = [REDATA.TSaxes, REDATA.AIaxes];
+
+% the mass and moment of inertia for the generic model are from
+% M.D. Klein Breteler et al., Journal of Biomechanics 32 (1999) 1191.
+% The body weight is not mentioned in this paper so I approximated it using
+% the R. Dumas et al., Journal of Biomechanics 40 (2007) 543?553.
+
+% Mass Parameters [Kg]
+Mc = 0.156;
+Ms = 0.704;
+Mh = 2.052; % The body weight of our male subject is Mh*100/2.4 (=85.5 Kg) according to R. Dumas
+M_u = 0.9012; % not in the paper by Breteler but close to "A multibody biomechanical model of the upper limb including the shoulder girdle"
+Mr = 0.5523; % not in the paper by Breteler "
+
+
+% Inertia Parameters [Kg M^2]
+Ict = 0.001; Icl = 0.003; % from before, but I do not like them as they are not consistent with the other segments values or the portugease paper
+Ist = 0.007; Isl = 0.007; % the same as above
+Iht = 0.0165; Ihl = 0.0033; % these values are quite different that what it was before.
+Iut = 0.0060; Iul = 9.6437e-04; % defined using the SUBJECT_TOOL_Update_Weight by giving 85.5kg and 1.86 m, and the values are comparable to those of Portugease
+Irt = 0.0037; Irl = 5.9107e-04; % the same as above
+
+DYDATA.MHand = 0;% The Radius mass can also contain a weight carried in the hand.
+DYDATA.LHand = 270e-3;% The inertia will be MHand*LHand^2, the position from CP;
+
+% Gravitational Constant
+g = 9.81;
+
+% ATTENTION!! One needs to apply the parallel axis theorem to be correct.
+ % The muscle has wrapping and the object must be rotated into
+ % the correct local bone frame.
+ if isempty(MWDATA{Mident,1}.MSCInfo.ObjectCentre) == 0
+ % rotate object centre
+ for idx = 1:size(MWDATA{Mident,1}.MSCInfo.ObjectRef,2)%we use the idx for the case of double cylinder when we have to objects with obviously two different centers
+if isequal(task, 'Build Rotation Matrices From Euler Angles')
+
+ % Get the Joint Angles
+ JEA = varargin{1,1};% varargin is a 1-by-N cell array, where N is the number of inputs that the function receives after the explicitly declared inputs
+
+ BLDATA = varargin{1,2}; % lately added to feed in the BLDATA for the rotation matrices of the ulna and radius regarding their proximal segments
+
+ % Clavicula Rotation Matrix
+ Rcx = [1.00, 0.0000, 0.0000;
+ 0.00, cos(JEA(1)), -sin(JEA(1));
+ 0.00, sin(JEA(1)), cos(JEA(1))];
+
+ Rcy = [cos(JEA(2)), 0.00, -sin(JEA(2));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(2)), 0.00, cos(JEA(2))];
+
+ Rcz = [cos(JEA(3)), -sin(JEA(3)), 0.00;
+ sin(JEA(3)), cos(JEA(3)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rc = Rcz*Rcy'*Rcx;
+
+ % Scapula Rotation Matrix
+ Rsx = [1.00, 0.0000, 0.0000;
+ 0.00, cos(JEA(4)), -sin(JEA(4));
+ 0.00, sin(JEA(4)), cos(JEA(4))];
+
+ Rsy = [cos(JEA(5)), 0.00, -sin(JEA(5));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(5)), 0.00, cos(JEA(5))];
+
+ Rsz = [cos(JEA(6)), -sin(JEA(6)), 0.00;
+ sin(JEA(6)), cos(JEA(6)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rs = Rsz*Rsy'*Rsx;
+
+ % Humerus Rotation Matrix
+ Rhz = [cos(JEA(7)), -sin(JEA(7)), 0.00;
+ sin(JEA(7)), cos(JEA(7)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rhy = [cos(JEA(8)), 0.00, -sin(JEA(8));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(8)), 0.00, cos(JEA(8))];
+
+ Rhzz = [cos(JEA(9)), -sin(JEA(9)), 0.00;
+ sin(JEA(9)), cos(JEA(9)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rh = Rhzz*Rhy'*Rhz;
+
+ %Ulna Rotation Matrix
+ Rux = [1.00, 0.00, 0.00;
+ 0.00, cos(JEA(10)), -sin(JEA(10));
+ 0.00, sin(JEA(10)), cos(JEA(10))];% ulna to the proximal segment part 1
+
+ Ru_h = BLDATA.Original_Matrices_L2A.Rh'*BLDATA.Original_Matrices_L2A.Ru; % ulna to the proximal segment part 1
+
+ Ru = Rh*Ru_h*Rux; % ulna to the inertial frame (thorax)
+
+ %Radius Rotation Matrix
+ Rrz = [cos(JEA(11)), -sin(JEA(11)), 0.00;
+ sin(JEA(11)), cos(JEA(11)), 0.00;
+ 0.00, 0.00, 1.00]; % radius to the proximal segment part 1
+
+ Rr_u = BLDATA.Original_Matrices_L2A.Ru'*BLDATA.Original_Matrices_L2A.Rr; % radius to the proximal segment part 2
+
+ Rr = Ru*Rr_u*Rrz; % radius to the inertial frame (thorax)
+elseif isequal(task, 'Update Initial Bony Landmark Data from Joint Rotation Matrices')
+ % Get the SC joint rotation matrix
+ Rc = varargin{1,1};
+
+ % Get the AC joint rotation matrix
+ Rs = varargin{1,2};
+
+ % Get the GH joint rotation matrix
+ Rh = varargin{1,3};
+
+ % Get the HU joint rotation matrix
+ Ru = varargin{1,4};
+
+ % Get the RU joint rotation matrix
+ Rr = varargin{1,5};
+
+ % Get the BLDATA input
+ BLDATA = varargin{1,6};
+
+ % Get Initial Rotation Matrices
+ Rci = BLDATA.Original_Matrices_L2A.Rc;
+ Rsi = BLDATA.Original_Matrices_L2A.Rs;
+ Rhi = BLDATA.Original_Matrices_L2A.Rh;
+ Rui = BLDATA.Original_Matrices_L2A.Ru;
+ Rri = BLDATA.Original_Matrices_L2A.Rr;
+
+ % Ge all the points (these points are in the thorax frame i.e. matlab frame)
+ IJ = BLDATA.Original_Points.IJ;
+ PX = BLDATA.Original_Points.PX;
+ T8 = BLDATA.Original_Points.T8;
+ C7 = BLDATA.Original_Points.C7;
+ SC = BLDATA.Original_Points.SC;
+ AC = BLDATA.Original_Points.AC;
+ AA = BLDATA.Original_Points.AA;
+ TS = BLDATA.Original_Points.TS;
+ AI = BLDATA.Original_Points.AI;
+ GH = BLDATA.Original_Points.GH;
+ HU = BLDATA.Original_Points.HU;
+ EL = BLDATA.Original_Points.EL;
+ EM = BLDATA.Original_Points.EM;
+ CP = BLDATA.Original_Points.CP;
+ US = BLDATA.Original_Points.US;
+ RS = BLDATA.Original_Points.RS;
+
+ % Set the Points
+ BLDATA.Initial_Points.IJ = IJ;
+ BLDATA.Initial_Points.PX = PX;
+ BLDATA.Initial_Points.T8 = T8;
+ BLDATA.Initial_Points.C7 = C7;
+ BLDATA.Initial_Points.SC = SC;
+ BLDATA.Initial_Points.AC = Rc*Rci'*(AC - SC) + BLDATA.Initial_Points.SC;% Rci'*(AC - SC) is AC in Clavicle coordinate, so with the rotation matrix Rc we can define its new position in the inertial frame as such.
+ 'string', '<html> <b>- SUBJECT SPECIFIC TOOL -----------------------------</b><br><p style="font-size: 90%;"> <i> Tool for adapting the generic model to a specific subject</i></p>',...
+ 'fontsize', 14,...
+ 'callback', MAIN_TOOL_script_generator('Open Subject Specific Tool'));
+
+% Push Button for opening the muscle wrapping toolbox
+MAINGUIHandle.Muscle_Tool_Button = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.85, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- MUSCLE WRAPPING TOOL ---------------------------------</b><br><p style="font-size: 90%;"> <i> Tool for defining the muscle segement structure</i></p>',...
+ 'string', '<html> <b>- JOINT INITIALIZATION TOOL -----------------------------</b><br><p style="font-size: 90%;"> <i> Tool for setting the initial configuration of the joints</i></p>',...
+ 'string', '<html> <b>- CHECK INDIVIDUAL MOMENT ARMS ---------------------</b><br><p style="font-size: 90%;"> <i>Tool for checking that the moment arms are continuous</i></p>',...
+ 'fontsize', 14,...
+ 'callback', MAIN_TOOL_script_generator('Check Moment Arms'));
+ 'string', '<html><b>- SAVE VISUALISATION ----------------------------------------</b><br><p style="font-size: 90%;"> <i> Saves the visualization to an EPS file</i></p>',...
+ 'string', '<html><b>- SAVE DATA ---------------------------------------------------------</b><br><p style="font-size: 90%;"> <i> Saves the current data structures</i></p>',...
+ 'fontsize', 14,...
+ 'callback', MAIN_TOOL_script_generator('Save All Data Structures'));
+
+% Push Button for Closing the GUI
+MAINGUIHandle.Close_GUI_Button = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.0, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html><b>- CLOSE GUI ---------------------------------------------------------</b><br><p style="font-size: 90%;"><i>Closes and saves the current data structures</i></p>',...
+This readme details how to use the EPFL musculoskeletal model through its GUI.
+This version of the model allows replicating patients/subjects with anatomical glenohumeral joint (no implant) as well as patients with anatomical total shoulder arthroplasty (aTSA).
+
+Author: ehsan.sarshari@epfl.ch
+Revision 1, 26 Nov 2018
+——————————————————————————————
+Some notes before getting started:
+
+- For running the model you need to have Matlab installed.
+- The model is developed such that it is consistent with different platforms (Unix and Windows).
+- It is strongly recommended to run the model on Matlab 2014 or 2015 to have the best performance.
+- Depending on your screen size/resolution you may have to resize the GUI windows to full screen view so that the menus/buttons are readable.
+
+——————————————————————————————
+By following the steps bellow, you would be able to:
+ - adapt/scale the model to a specific patient
+ - simulate a measured motion
+ - estimate the muscle and joint reaction forces
+ - understand the overall structure of the model
+
+Now, let’s get started.
+
+1- Open your Matlab.
+
+2- Set the directory to /ShoulderModel.
+
+3- Open MODEL_MAIN_FILE.m and run it.
+
+The MasterHead window appears first including information about the contributors in developing the model. You can close it right-away.
+
+Meantime the model is constructed and two messages are printed one after another in the Command Window as the following.
+
+Building the model data, please wait ...
+The model has been constructed, ...
+
+Following the two messages, the main GUI window (EPFL-LBO Upper Extremity Model MAIN GUI) now appears.
+
+In this window, access to six model toolbox is granted. The toolboxes are:
+
+ - SUBJECT SPECIFIC TOOL
+ To adapt/scale the generic model to a specific patient/subject
+
+ - MUSCLE WRAPPING TOOL
+ To visually check the path of the muscles in different joint configurations and set the number of strings by which a muscle should be replicated.
+
+ - JOINT INITIALIZATION TOOL
+ To change the initial configuration of the model in terms of its joint angles.
+
+ - KINEMATICS TOOL
+ To reconstruct a motion either when measured data are available from videogrammetry systems or when no measured data is available.
+
+ - CHECK INDIVIDUAL MOMENT ARMS
+ To calculate the moment arms of all the muscles around each joint for the reconstructed motion.
+
+ - ESTIMATE FORCES
+ To transform the measured EMG signals to muscle forces and estimate the muscle and joint reaction forces using inverse dynamics and static optimization.
+
+The main GUI window also provides interactive options for changing/saving the model visualizations.
+
+You can close the model and clean up the Workspace by hitting the CLOSE GUI button.
+
+4- Click on the “SUBJECT SPECIFIC TOOL” button.
+
+SUBJECT SPECIFIC TOOLBOX window appears that allows specifying numbers of subject-specific parameters so that the generic model can be adapted to a specific patient/subject. These parameters are as the following:
+
+ - Gender
+ - Weight
+ - Height
+ - PCSA of muscles
+ - Glenoid fossa orientation defined using glenoid version angle, glenoid inclination angle, glenoid center point, and implant/humeral head diameter (a definition of these parameters is available in “Terrier, A., et al. Measurements of three-dimensional glenoid erosion when planning the prosthetic replacement of osteoarthritic shoulders, The bone & joint journal 96.4 (2014): 513-518).
+
+A protocol named MSM_input_output_preprocessing.pdf/docx is available in smb://lbovenus/shoulder/methods/matlab/MSM/ShoulderModel_pre_post_processing in order to allow you feed the subject-specific tool with patients’/subjects’ data obtainable from CT or blueprint data. For patients with anatomical glenohumeral joint refer to section 1-1 and for patients with anatomical total shoulder arthroplasty refer to section 1-2.
+
+5- Once you input the subject-specific parameters, you can load videogrametry-based measured kinematics by hitting the gray button below “Load Measured Kinematics”.
+
+A set of measured motions in their raw format (without filtration) are provided together with the model that can be obtained in /ShoulderModel/Generic_Measurements/kinematics. It includes the following activities in a desired format/structure that the model can read.
+
+ act1: abduction frontal plane with 2 kg
+ act2: elevation sagittal plane with 2 kg
+ act3: abduction scapula plane with 2 kg
+ act4: fast abduction scapula plane
+ act5: slow abduction scapula plane
+ act6: put 2 kg in a shelf at head height
+ act7: hand behind the head
+ act8: touch the other shoulder
+ act9: french canes
+ act10: counter external rotation (static pose, no motion)
+ act11: counter internal rotation (static pose, no motion)
+ act12: french cranes a second try
+
+Once the kinematic data is loaded, a message box will pop up to show if the provided data was consistent or something was missing or went wrong. You can simply close it in case it is SUCCESSFUL.
+
+6- Click on the “SCALE MEASURED KINEMATICS” to scale the generic measurements to the specific subject that you input his/her data.
+
+“The measured kinematics were scaled.” will be printed once this is done.
+
+7- Click “SCALE USING GENDER, BW, BH” to scale the model dynamics and kinematics based on gender, weight, and hight and according to the anthropometric studies that the model uses in the background.
+
+“The model skeletal morphology and muscle architecture and propertis were scaled” will be printed in the Command Window once this is done.
+
+8- Click “SCALE RIBCAGE ELLIPSOID” to scale the ellipsoids approximating the ribcage.
+
+These ellipsoids are used to constrain the scapula to gliding over the ribcage (for more detail see “Sarshari, Ehsan, A Closed-Loop EMG-Assisted Shoulder Model, No. THESIS. EPFL, 2018” and “Ingram, David, Musculoskeletal Model of the Human Shoulder for Joint Force Estimation, No. THESIS_LIB. EPFL, 2015”.
+
+This step takes longer than the other last two steps and when it is done the following message is printed in the Command Window.
+“Ribcage ellipsoides and dynamic model were scaled…”.
+
+9- Click “VISUALIZE THE SCALED DATA” to have a visual comparison between the scaled-generic model (will be shown in blue) and the generic model (shown in red).
+
+10- Click the “CLOSE TOOL” button to permanently save the subject-specific changes in the model. Otherwise, you can leave this window open so that in case you want to reset the model to its generic form, you can hit “RESET TO GENERIC MODEL” button. Evidently, you can always start from step 3 to have the generic model.
+
+11- Hit the “MUSCLE WRAPPING TOOL” in the MAIN GUI window.
+
+The “MUSCLE WRAPPING TOOLBOX” appears.
+
+By default all the muscles are replicated with one single segment shown in red. The origin of the muscles are shown on the bones in green and the insertions are in blue.
+
+On the right hand-side of this window you can vary 11 joint angles to visually check the path taken by all the muscles.
+
+You can reset the joint angles to the model’s initial configuration by hitting the “RESET ANGLES” button.
+
+By default all the muscles, their origins, and their insertions are shown. You can change the visualization of each muscle or all of the muscles through “Muscle Properties” menu and “Set All Properties” menu in the menu bar, respectively. You need to hit the “UPDATE VISUALIZATION” button to put the changes in effect.
+
+You can also visualize the wrapping object(s) of each muscle through Muscle Properties menu. First select the muscle, and tick the “Object Visible?” option. Then, click on the “UPDATE VISUALIZATION” button.
+
+The obstacle set method of Garner and Pandy with some modifications in implementation was used to construct the muscle paths, while the insertions, origins, and the wrapping objects were defined based on our generic subject’s MRI and CTs with the help of a professional radiologist. For more details see “Garner, B. A., & Pandy, M. G. (2000). The obstacle-set method for representing muscle paths in musculoskeletal models. Computer methods in biomechanics and biomedical engineering, 3(1), 1-30.” and “Ingram, D. (2015). Musculoskeletal Model of the Human Shoulder for Joint Force Estimation (No. THESIS_LIB). EPFL.”.
+
+12- On the “MUSCLE WRAPPING TOOLBOX” window, click on the “Set All Properties” menu and set the “Number of Segments” option to 3. Then click on the “UPDATE VISUALIZATION” button. This would increase the number of segments used for replicating each one of the muscles from 1 to 3.
+
+13- Click on the “CLOSE TOOL” button to close the “MUSCLE WRAPPING TOOLBOX”.
+
+14- Click the “KINEMATICS TOOL” in the Main GUI window.
+
+The “KINEMATICS TOOLBOX” window appears.
+
+On the top right-hand corner, you have two options to reconstruct a motion.
+
+ - Option 1: Use this option in case you do not have videogrammetry-based measured motions. This option constructs the joint space using 9 minimal coordinates that inherently satisfy the kinematic constraints regarding scapula gliding over the ribcage ellipsoids.
+
+The minimal coordinates (M1 to M9) are detailed in the bottom of the top right-hand corner of the “KINEMATICS TOOLBOX”.
+
+To construct a motion using this option you should first set the initial and final values of the minimal coordinates through the boxes available on the left hand-side corner of the toolbox. Then set the number of frames that you want to have between the initial and final values through “Number Of Points”, i.e. in how many steps (frames) you want the model to go from the configuration defined by the initial values to the configuration defined using the final values of the minimal coordinates. Once you set these, hit the “BUILD MOTION (Option 1)” button to construct a motion.
+
+The constructed motion would appear in terms of 11 joint angles on the lower half of the window. You can click on each graph to have it open in a separate window.
+
+In order to replicate a specific motion using the minimal coordinates you need to set proper initial and final values for them. To this end, you require a more thorough understanding of these coordinates and also some tweaking afterwards. To find out more about the minimal coordinates check the followings “Ingram, D. (2015). Musculoskeletal Model of the Human Shoulder for Joint Force Estimation (No. THESIS_LIB). EPFL.” and “Ingram, D., Engelhardt, C., Farron, A., Terrier, A., & Müllhaupt, P. (2016). Modelling of the human shoulder as a parallel mechanism without constraints. Mechanism and Machine Theory, 100, 120-137.”.
+
+ - Option 2: Once you loaded a videogrammetry-based measured kinematics in the “SUBJECT-SPECIFIC TOOLBOX” you can use this option to reconstruct a motion. Motion reconstruction here refers to deriving the joint angles from the videogrammetry-based measured trajectories of skin-fixed markers. The joint angles are defined such that the sum of distance between the model landmarks and the measured markers is minimized for each frame of measured data while the kinematic constraint of scapula gliding over the ribcage is satisfied.
+
+15- Click on the “BUILD MOTION (Option 2)” button.
+
+“BUILD MOTION (Option 2), reconstruction of a measured motion” window appears.
+
+This window allows you to
+ - analyze the measured data to design a proper filter
+ - design a filter
+ - filter the data
+ - estimate missing markers (see step 18)
+ - perform inverse-kinematics
+ - see the reconstructed motion in terms of an animated motion as well as the joint angles graphs.
+
+16- Click on either “residual Analysis” or “Harmonic Analysis” buttons to analyze the measured kinematic data in terms of the amount of noise they have. For each case a graph showing the results of the analysis will appear. Based on these graphs you can define the cutoff frequency of the low-pass filter. If you are not familiar with these two standard analyses, you can find more details in “Winter, D. A. (2009). Biomechanics and motor control of human movement. John Wiley & Sons.”.
+
+17- Set the “Cutoff Frequency” and “Filter Order” in the associated boxes and then hit the “Filter Measured Data” button.
+
+A success message would pop up, that you can close, and a graph showing the filtered data together with the raw data. You can use this graph to have an idea of how long the measurements are and set a proper time for your simulation in the “Motion time span” box. Otherwise, you can keep the default value. It is recommended to consider a proper value to avoid simulation of the measured motion while the subject was in the transition pause.
+
+18- Click on the “Estimate Missing Landmarks”.
+
+In the measurements provided with the model, only the trajectories of 11 markers were measured. However, the model consists of 14 landmarks. The 3 missing markers were associated with the center of the glenohumeral joint (GH) and two points on the border of the scapula (TS and AI). These three markers cannot effectively be palpated. On the other hand, in order to be able to evaluate the distance between the model landmarks and the measured markers a one-to-one association between the model landmarks and the measured markers is required. To this end, the missing markers are estimated. Our approach in estimating these markers can be found here “Sarshari, E. (2018). A Closed-Loop EMG-Assisted Shoulder Model (No. THESIS). EPFL.”. Alternatively, if a set of kinematic data including these missing landmarks is provided for instance by using a scapula kinematic measurement-device, the estimation step can be simply skipped.
+
+
+While the estimation is on course, the progress is printed in the Command Window. Once the estimation is finished, the following message is printed in the Command Window “Show the estimated/measured landmarks in GREEN for the first time Stamp” and the configuration of the subject for the first frame of the measured data is displayed in green together with the model configuration (based on the posture of the generic subject in the MRI machine in blue).
+
+19- Set a proper/meaningful time span for the simulation in the “Motion time Span” that excludes portions of the motion that are not necessary/interesting. Practically speaking, there is always a few seconds at the end of our measurements where the subject is not performing any motion and instead waiting for the indication to perform the next task.
+
+20- Hit the “Perform Inverse Kinematics” button to reconstruct the motion.
+
+It takes some time to complete and meanwhile the progress can be observed through the Command Window.
+
+Once it is completed, the following messages are printed in the Command Window “Motion reconstruction was finished and the results smoothened.
+Plot the resulted generalized coordinates (joint angels)” and the 11 joint angles are plotted in the lower half of the window.
+
+Each one of the plots will be plotted separately upon clicking on them. Furthermore, a smooth fitted polynomial and its first and second order derivatives will be also shown in the same figure.
+
+Click on the “View Motion” button if you would like to see the reconstructed motion. A video with AVI format will be save in the main directory of the model at the same time.
+
+21- Click on the “CLOSE TOOL” and proceed to the MAIN GUI window.
+
+
+22- Click on the “ESTIMATE FORCES” button on the MAIN GUI window.
+
+The “MUSCEL FORCE ESTIMATION TOOLBOX” window appears.
+
+This window allows you to estimate the muscle and joint reaction forces. There are three options possible for performing this estimation as the following that can be activated through the “Estimation Options/Options” menu in the menu bar.
+
+ - Glenohumeral Stability Constraint: forces the glenohumeral joint reaction force to always point toward the glenoid fossa, more details here “Ingram, D. (2015). Musculoskeletal Model of the Human Shoulder for Joint Force Estimation (No. THESIS_LIB). EPFL.” and “Sarshari, E. (2018). A Closed-Loop EMG-Assisted Shoulder Model (No. THESIS). EPFL.”.
+
+ - EMG-Assisted: introduces the muscle forces of 15 muscles calculated based on the measured EMG excitations and a Hill-type muscle model as upper and lower bounds in the load-sharing problem, more details in Step 24 and here “Sarshari, E. (2018). A Closed-Loop EMG-Assisted Shoulder Model (No. THESIS). EPFL.”.
+
+ - Measured Kinematics: considers the reconstructed motion from the measured kinematics in the inverse dynamics.
+
+
+23- Check the tick for the three options in “Estimation Options/Options” menu in the menu bar.
+
+24- Hit the gray button below the “Load EMG data”.
+
+This allows you to load the EMG excitations of 15 superficial muscles that were measured on the same subject during the same motion. The measured EMG data can be found in ShoulderModel/Generic_Measurements/EMG.
+
+The idea is to use the measured EMG excitations to estimate the muscle forces of 15 muscles using a Hill type muscle model for the reconstructed motion. These EMG-based muscle forces will be used as upper and lower bounds in the load-sharing to shrink its feasible set. We have shown that this can improve the model’s muscle force estimations.
+
+25- Select the measured EMG data for the associating motion (activity).
+
+A success message would appear that you can dismiss.
+
+Note that the EMG data have been already processed, more details here “Sarshari, E. (2018). A Closed-Loop EMG-Assisted Shoulder Model (No. THESIS). EPFL.”.
+
+26- Set the time of the simulation through the box under “Motion Time Span” according to the time that you set in Step 19.
+
+27- If the subject carried a weight during the measurements you should indicate it through the box under the “Weight in Hand”.
+
+28- Click on the “EMG to Force” button.
+
+This will calculate the muscle forces for each one of the 15 muscles for which you provided EMG signals by integrating a Hill type muscle-tendon model. More details here “Sarshari, E. (2018). A Closed-Loop EMG-Assisted Shoulder Model (No. THESIS). EPFL.”.
+
+You can follow the progress in the Command Window given that it might take quite some time depending on the motion time span.
+
+Once it is completed, a graph showing the EMG based muscle forces will be plotted.
+
+You can dismiss this graph and proceed to the next step.
+
+29- Click the “ESTIMATE FORCES” button to solve the load-sharing problem given the 3 options that you selected.
+
+The progress can be observed through the Command Window.
+
+Once it is completed, a message box would pop up that can assist you to check if the load-sharing was performed successfully.
+
+The muscle and joint reaction forces will be plotted on the “MUSCEL FORCE ESTIMATION TOOLBOX” window. By clicking on each one of these plots you can visualize them separately.
+
+The save button saves the muscle and joint reaction force estimations results (JRDATA and EFDATA respectively) in the folder “Data_Structures_and_Documentation”.
+
+30- The following details the simulations outcome that you find in the Workspace.
+
+ - SSDATA: includes all subject-specific data, all data associating with the measurements (EMG/kinematics), and the model kinematics (joint angles) if inverse kinematics (Option 2) is used to reconstruct a motion.
+ - KEDATA: includes all the data associating with the model kinematics if minimal coordinates (Option 1) is used to construct a motion.
+ - JRDATA: is a 16 column vertical matrix that includes the 3 components of the GH joint reaction force in the thorax coordinate system, the magnitude of the GH joint reaction force, the three coordinates of the intersection of the GH joint reaction force with the glenoid fossa in the fossa coordinate system, the stability ratio, the 3 components of the HU joint reaction force, the magnitude of the HU joint reaction force, the 3 components of the RU joint reaction force and the magnitude of the RU joint reaction force, respectively.
+ - EFDATA: includes estimations of muscle forces and their moment arms.
+ - DYDATA: includes all the data associating with the model inertial properties.
+ - BLDATA: includes all the data associating with the model morphologies. The coordinate systems of different bone segment can be found here.
+ - MWDATA: includes all the data associating with the muscle paths and their physiological constants.
+ - REDATA: includes all the data associating with the ribcage ellipsoids.
+ 'FLCU';'INFS']; % here is the full list of the landmarks that indeed here 'IJ';'PX';'T8';'C7';'SC';'AC';'GH';'HU';'TS';'AI';'AA';'EL';'EM';'CP';'US';'RS'
+ % NOTE: if the imported list has more landmarks than what it's specified
+ % it's not a problem. If, once I wanted to differentiate with that to
+ % force the accepting only the files with exact number and name of landmarks
+ % we have to first use fieldnames(filename) to have them in a cell array and
+ % then draw the parallel within the for loop as it is now. In otehr words,
+ % the eval only allows a whole to an invidual comparision not a one to one
+ % comparision.
+
+% make sure that a right file is chosen by the user
+[~, filename, ext] = fileparts(fullfileName); % returns the path, file name, and file name extension for the specified "filename"
+
+ if strcmp(ext,'.mat') % if the extension is mat then load it otherwise do nothing
+ load(fullfileName);
+
+ % Check whether the chosen file contains the required fields.
+ for muscle_check_id = 1:length(required_muscle_list)
+
+ if ~isfield(eval(filename),required_muscle_list(muscle_check_id,:))
+ fprintf('muscle %s was not found in the loaded data. \n Please try a different EMG data file and make sure it \n has the requirements mentioned in the help menu. \n', required_muscle_list(muscle_check_id,:));
+ %disp(['Bony landmark ''' required_landmarks_list(landmarks_check_id,:) ''' was not found in the loaded data',...
+ % '.\n Please try a different kinematic file after reading the file requirement from the help meneu']);
+% Get initial and current scapula rotation matrices
+Rs0 = BLDATA.Initial_Matrices_L2A.Rs;
+Rs = BLDATA.Current_Matrices_L2A.Rs;
+
+
+% Build the normal vectors
+for i = 1:size(Nout,1)
+ % first make the normal vectors for the generic model
+ % points on the fossa for the generic model in the cone frame
+ C1 = C(1);
+ y = P(2,i);
+ z = P(3,i);
+
+ % old method (the same but more calculations are required on paper)
+ %dxdy = -(C1^2*y)/(DYDATA.ConeDimensions(1)^2*((C1^2*y^2)/DYDATA.ConeDimensions(1)^2 + (C1^2*z^2)/DYDATA.ConeDimensions(2)^2)^(1/2)); % the minus is due to the position of C_x that is defined as ((C1^2*y^2)/DYDATA.ConeDimensions(1)^2 + (C1^2*z^2)/DYDATA.ConeDimensions(2)^2)^(1/2) however we could have just put C1 here and it's fine. For the details I have the derivations and one can see them.
+ title('intersection locus of the GHJ reaction force and the glenoid fossa (green o = start, red o = end)', 'fontsize', 12, 'fontname', 'sansserif', 'fontweight', 'bold');
+ title('intersection locus of the GHJ reaction force and the glenoid fossa (green o = start, red o = end)', 'fontsize', 12, 'fontname', 'sansserif', 'fontweight', 'bold');
+ 'string', '<html> <b>- CLOSE GUI -------------------------------------------------</b><br><p style="font-size: 90%;"> <i>Closes this GUI and updates BLDATA</i></p>',...
+ 'string', '<html> <b>- SAVE VISUALISATION ------------------------------------</b><br><p style="font-size: 90%;"> <i>Saves the current axis content</i></p>',...
+ 'string', '<html><b>- SAVE VISUALISATION ----------------------------</b><br><p style="font-size: 90%;"><i>Save''s the GUI''s current visualisation.</i></p>',...
+ 'string', '<html><b>- RESET CONE DATA -------------------------------</b><br><p style="font-size: 90%;"><i>Restes the data to original values.</i></p>',...
+ 'string', '<html> <b>- CLOSE TOOL ---------------------------------------</b><br><p style="font-size: 90%;"> <i>Closes the tool and deletes all the local variables</i></p>',...
+ 'string', '<html> <b>- VIEW MOTION -------------------------------------</b><br><p style="font-size: 90%;"> <i>Launches a small GUI where the motion is visualised</i></p>',...
+% Text for captioning the reconstructed generalized coordinates
+KRGUIHandle.GCTxt = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.29, 0.45, 0.5, 0.03],...
+ 'style', 'text',...
+ 'string', 'Reconstructed generalized coordinates (click on each of the plots below for further plots including the associated generalized velocities and accelerations)',...
+ disp(['Motion reconstruction, Time Stamp : ', num2str(TimeId), '/', num2str(length(SSDATA.Measured_Kinematics.IJ)-1)]);
+
+ % initial guess for fmincon
+% if TimeId == 1
+% q0 = zeros(11,1);%KEDATA.Initial_Minimal_Coordinate % the optimization performance highly depends on the initial guess, subject calibration therefore plays a crucial role in the accuracy of the inverse kinematics method used here to reconstruct a measured motion.
+% else
+% q0 = q_r(:,TimeId-1);
+% end
+
+ if TimeId == 1
+ p0 = zeros(88,1);%KEDATA.Initial_Minimal_Coordinate % the optimization performance highly depends on the initial guess, subject calibration therefore plays a crucial role in the accuracy of the inverse kinematics method used here to reconstruct a measured motion.
+% Text for captioning the reconstructed generalized coordinates
+KRGUIHandle.GCTxt = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.29, 0.45, 0.5, 0.03],...
+ 'style', 'text',...
+ 'string', 'Reconstructed generalized coordinates (click on each of the plots below for further plots including the associated generalized velocities and accelerations)',...
+% Push Button for filtering the measured kinematics
+KRGUIHandle.Filter_Button = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.78, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- Filter Measured Data -----------------------------</b><br><p style="font-size: 90%;"> <i> Filters and transforms the measured data provided the necessary settings</i></p>',...
+ 'string', '<html> <b>- Perform Inverse Kinematics --------------------</b><br><p style="font-size: 90%;"> <i> Defines the generalized coordinates such that the model tracks the measured data</i></p>',...
+% Push Button for viewing the reconstructed motion together with the
+% measured motion
+KRGUIHandle.View_Button = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.62, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- View Motions ---------------------------------------</b><br><p style="font-size: 90%;"> <i> Shows the reconstructed motion resulted from IK together with the meaured+estimated motion</i></p>',...
+% this time is not really important during the IK, except when we want to
+% plot the resulted generalized coordinates along with time or when we want
+% to define/plot their differntiation i.e. generalized velocities and
+% accelerations. Otherwise, this value is of importance for force
+% prediction. As it defines the dynamics of motion (accelerations and
+% velocities).
+
+% Text for imposing the final time (we know it from measurements)
+KRGUIHandle.TimeFrameTxt = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.57, 0.68, 0.13, 0.04],...
+ 'style', 'text',...
+ 'string', 'Motion time span [S]:',...
+ 'backgroundcolor', 'white',...
+ 'fontsize', 12,...
+ 'HorizontalAlignment','left',...
+ 'fontweight', 'bold');
+
+% Edit for imposing the final time (we know it from measurements)
+KRGUIHandle.TimeFrameEdit = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.57, 0.68, 0.12, 0.02],...
+ 'style', 'edit',...
+ 'string', num2str(length(SSDATA.Measured_Kinematics.IJ)/100),... % the number of sample devided by the fs (sampling frequency) gives the total time span of the motion
+ 'string', '<html> <b>- CLOSE TOOL -------------------------------------</b><br><p style="font-size: 90%;"> <i> Closes the ellipsoid tool and returns to main GUI</i></p>',...
+ 'string', '<html> <b>- SAVE VISUALISATION ---------------------------</b><br><p style="font-size: 90%;"> <i> Saves the current visualisation to en .eps file</i></p>',...
+ 'AA';'EL';'EM';'US';'RS']; % here is the full list of the landmarks that indeed here 'IJ';'PX';'T8';'C7';'SC';'AC';'GH';'HU';'TS';'AI';'AA';'EL';'EM';'CP';'US';'RS'
+ % NOTE: if the imported least has more landmarks than what it's specified
+ % it's not a problem. If, once I wanted to differentiate with that to
+ % force the accepting only the files with exact number and name of landmarks
+ % we have to first use fieldnames(filename) to have them in a cell array and
+ % then draw the parallel within the for loop as it is now. In otehr words,
+ % the eval only allows a whole to an invidual comparision not a one to one
+ % comparision.
+
+% make sure that a right file is chosen by the user
+[~, filename, ext] = fileparts(fullfileName); % returns the path, file name, and file name extension for the specified "filename"
+
+ if strcmp(ext,'.mat') % if the extension is mat then load it otherwise do nothing
+ load(fullfileName);
+
+ % Check whether the chosen file contains the required fields.
+ for landmarks_check_id = 1:length(required_landmarks_list)
+
+ if ~isfield(eval(filename),required_landmarks_list(landmarks_check_id,:))
+ fprintf('Bony landmark %s was not found in the loaded data. \n Please try a different kinematic file and make sure it \n has the requirements mentioned in the help menu. \n', required_landmarks_list(landmarks_check_id,:));
+ %disp(['Bony landmark ''' required_landmarks_list(landmarks_check_id,:) ''' was not found in the loaded data',...
+ % '.\n Please try a different kinematic file after reading the file requirement from the help meneu']);
+ DYDATA.Inertia(1,4) = 0.62*1.7e-2*Subject_Weight; % Mu: This one is a bit messy, because the forearm is mentioned in the BSIP studies and I deceided to consider 62% of the forearm mass as the ulna
+ DYDATA.Inertia(1,5) = 0.38*1.7e-2*Subject_Weight; % Mr: The same as above.
+% DYDATA.Inertia(2,1) = % Ict: Clavicle and scapula keep their generic inertia values
+% DYDATA.Inertia(2,2) = % Icl
+% DYDATA.Inertia(3,1) = % Ist
+% DYDATA.Inertia(3,2) = % Isl
+
+ DYDATA.Segments_Lenghts.Humerus = (27.1/177)*Subject_Height; % scale the humerus lenght according to the body height
+ DYDATA.Segments_Lenghts.Ulna = (28.3/177)*Subject_Height; % scale the ulna lenght according to the body height
+ DYDATA.Segments_Lenghts.Radius = (28.3/177)*Subject_Height; % scale the radius lenght according to the body height
+ %P_inT = DYDATA.Cone_Rb*R_opt*P(:,i)*1e3 + BLDATA.Initial_Points.GH; % first rotate the fossa points by the R_opt and then transfer them to the thorax frame to be drawn
+S1 = Height_Scaling_Factor*[-4.437187576293945e+001; 1.044776678085327e+000; 1.171667938232422e+002]; % the most lateral point on the supraspinatus fossa
+S2 = Height_Scaling_Factor*[-2.524987030029297e+001; 3.014913558959961e+001; 1.212621765136719e+002]; % one of the other 4 points on the supraspinatus fossa
+S3 = Height_Scaling_Factor*[-4.073601150512695e+001; 9.907344818115234e+000; 1.174581298828125e+002]; % one of the other 4 points on the supraspinatus fossa
+S4 = Height_Scaling_Factor*[-3.028522109985352e+001; 2.512048149108887e+001; 1.201423034667969e+002]; % one of the other 4 points on the supraspinatus fossa
+S5 = Height_Scaling_Factor*[-3.626570129394531e+001; 1.977687454223633e+001; 1.190045471191406e+002]; % one of the other 4 points on the supraspinatus fossa
+% Ribcage Ellipsoid Data for scapulothoracic constraints
+DYDATA.OE = REDATA.Centre;
+DYDATA.AE = [REDATA.TSaxes, REDATA.AIaxes];
+
+% the mass and moment of inertia for the generic model are from
+% M.D. Klein Breteler et al., Journal of Biomechanics 32 (1999) 1191.
+% The body weight is not mentioned in this paper so I approximated it using
+% the R. Dumas et al., Journal of Biomechanics 40 (2007) 543?553.
+%{
+% Mass Parameters [Kg]
+Mc = 0.156;
+Ms = 0.704;
+Mh = 2.052; % The body weight of our male subject is Mh*100/2.4 (=85.5 Kg) according to R. Dumas
+M_u = 0.9012; % not in the paper by Breteler but close to "A multibody biomechanical model of the upper limb including the shoulder girdle"
+Mr = 0.5523; % not in the paper by Breteler "
+
+
+% Inertia Parameters [Kg M^2]
+Ict = 0.001; Icl = 0.003; % from before, but I do not like them as they are not consistent with the other segments values or the portugease paper
+Ist = 0.007; Isl = 0.007; % the same as above
+Iht = 0.0165; Ihl = 0.0033; % these values are quite different that what it was before.
+Iut = 0.0060; Iul = 9.6437e-04; % defined using the SUBJECT_TOOL_Update_Weight by giving 85.5kg and 1.86 m, and the values are comparable to those of Portugease
+Irt = 0.0037; Irl = 5.9107e-04; % the same as above
+
+DYDATA.MHand = 0;% The Radius mass can also contain a weight carried in the hand.
+DYDATA.LHand = 270e-3;% The inertia will be MHand*LHand^2, the position from CP;
+
+% Gravitational Constant
+g = 9.81;
+
+% ATTENTION!! One needs to apply the parallel axis theorem to be correct.
+ % The muscle has wrapping and the object must be rotated into
+ % the correct local bone frame.
+ if isempty(MWDATA{Mident,1}.MSCInfo.ObjectCentre) == 0
+ % rotate object centre
+ for idx = 1:size(MWDATA{Mident,1}.MSCInfo.ObjectRef,2)%we use the idx for the case of double cylinder when we have to objects with obviously two different centers
+ 'string', '<html> <b>- CLOSE TOOL -------------------------------------</b><br><p style="font-size: 90%;"> <i> Closes the subject specific tool and returns to main GUI</i></p>',...
+% Push Button for resetting the values enteretd by the user to the generic
+% model
+SGUIHandle.Reset_Button = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.3, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- RESET TO GENERIC MODEL -------------------</b><br><p style="font-size: 90%;"> <i>Resets the anthropometry data entered to those of the generic model</i></p>',...
+ 'fontsize', 14,...
+ 'fontweight', 'bold',...
+ 'callback', SUBJECT_TOOL_script_generator('RESET TO GENERIC MODEL'));
+
+% Push Button for visualization of the scaled data
+SGUIHandle.Visulization = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.38, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- VISULIZE THE SCALED DATA ---------------</b><br><p style="font-size: 90%;"> <i>Visulizes the scaled data together with the generic model</i></p>',...
+% Push Button for updating everything based on the gender, BW, and BH
+SGUIHandle.Update = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.48, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- SCALE USING GENDER, BW, BH ------------</b><br><p style="font-size: 90%;"> <i>Scales the model based on the patient specific data provided above</i></p>',...
+% Scale based on Measured Kinematics data (Pushbutton)
+% SGUIHandle.Visulization = uicontrol(...
+% 'units', 'normalized',...
+% 'position', [0.7, 0.55, 0.3, 0.05],...
+% 'style', 'pushbutton',...
+% 'string', '<html> <b>- SCALE USING MEASURED KINEMATICS-----</b><br><p style="font-size: 90%;"> <i>Scale using kienamtics data uploaded by user. It overwites scaling based on subject height.</i></p>',...
+% Scale based on Measured Kinematics data (Pushbutton)
+SGUIHandle.Kinematic_scaling = uicontrol(...
+ 'units', 'normalized',...
+ 'position', [0.7, 0.55, 0.3, 0.05],...
+ 'style', 'pushbutton',...
+ 'string', '<html> <b>- SCALE MEASURED KINEMATICS-----</b><br><p style="font-size: 90%;"> <i>Scale measured kienamtics data recorded on the generic subject.</i></p>',...
+if isequal(task, 'Build Rotation Matrices From Euler Angles')
+
+ % Get the Joint Angles
+ JEA = varargin{1,1};% varargin is a 1-by-N cell array, where N is the number of inputs that the function receives after the explicitly declared inputs
+
+ BLDATA = varargin{1,2}; % lately added to feed in the BLDATA for the rotation matrices of the ulna and radius regarding their proximal segments
+
+ % Clavicula Rotation Matrix
+ Rcx = [1.00, 0.0000, 0.0000;
+ 0.00, cos(JEA(1)), -sin(JEA(1));
+ 0.00, sin(JEA(1)), cos(JEA(1))];
+
+ Rcy = [cos(JEA(2)), 0.00, -sin(JEA(2));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(2)), 0.00, cos(JEA(2))];
+
+ Rcz = [cos(JEA(3)), -sin(JEA(3)), 0.00;
+ sin(JEA(3)), cos(JEA(3)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rc = Rcz*Rcy'*Rcx;
+
+ % Scapula Rotation Matrix
+ Rsx = [1.00, 0.0000, 0.0000;
+ 0.00, cos(JEA(4)), -sin(JEA(4));
+ 0.00, sin(JEA(4)), cos(JEA(4))];
+
+ Rsy = [cos(JEA(5)), 0.00, -sin(JEA(5));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(5)), 0.00, cos(JEA(5))];
+
+ Rsz = [cos(JEA(6)), -sin(JEA(6)), 0.00;
+ sin(JEA(6)), cos(JEA(6)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rs = Rsz*Rsy'*Rsx;
+
+ % Humerus Rotation Matrix
+ Rhz = [cos(JEA(7)), -sin(JEA(7)), 0.00;
+ sin(JEA(7)), cos(JEA(7)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rhy = [cos(JEA(8)), 0.00, -sin(JEA(8));
+ 0.0000, 1.00, 0.0000;
+ sin(JEA(8)), 0.00, cos(JEA(8))];
+
+ Rhzz = [cos(JEA(9)), -sin(JEA(9)), 0.00;
+ sin(JEA(9)), cos(JEA(9)), 0.00;
+ 0.0000, 0.0000, 1.00];
+
+ Rh = Rhzz*Rhy'*Rhz;
+
+ %Ulna Rotation Matrix
+ Rux = [1.00, 0.00, 0.00;
+ 0.00, cos(JEA(10)), -sin(JEA(10));
+ 0.00, sin(JEA(10)), cos(JEA(10))];% ulna to the proximal segment part 1
+
+ Ru_h = BLDATA.Original_Matrices_L2A.Rh'*BLDATA.Original_Matrices_L2A.Ru; % ulna to the proximal segment part 1
+
+ Ru = Rh*Ru_h*Rux; % ulna to the inertial frame (thorax)
+
+ %Radius Rotation Matrix
+ Rrz = [cos(JEA(11)), -sin(JEA(11)), 0.00;
+ sin(JEA(11)), cos(JEA(11)), 0.00;
+ 0.00, 0.00, 1.00]; % radius to the proximal segment part 1
+
+ Rr_u = BLDATA.Original_Matrices_L2A.Ru'*BLDATA.Original_Matrices_L2A.Rr; % radius to the proximal segment part 2
+
+ Rr = Ru*Rr_u*Rrz; % radius to the inertial frame (thorax)
+elseif isequal(task, 'Update Initial Bony Landmark Data from Joint Rotation Matrices')
+ % Get the SC joint rotation matrix
+ Rc = varargin{1,1};
+
+ % Get the AC joint rotation matrix
+ Rs = varargin{1,2};
+
+ % Get the GH joint rotation matrix
+ Rh = varargin{1,3};
+
+ % Get the HU joint rotation matrix
+ Ru = varargin{1,4};
+
+ % Get the RU joint rotation matrix
+ Rr = varargin{1,5};
+
+ % Get the BLDATA input
+ BLDATA = varargin{1,6};
+
+ % Get Initial Rotation Matrices
+ Rci = BLDATA.Original_Matrices_L2A.Rc;
+ Rsi = BLDATA.Original_Matrices_L2A.Rs;
+ Rhi = BLDATA.Original_Matrices_L2A.Rh;
+ Rui = BLDATA.Original_Matrices_L2A.Ru;
+ Rri = BLDATA.Original_Matrices_L2A.Rr;
+
+ % Ge all the points (these points are in the thorax frame i.e. matlab frame)
+ IJ = BLDATA.Original_Points.IJ;
+ PX = BLDATA.Original_Points.PX;
+ T8 = BLDATA.Original_Points.T8;
+ C7 = BLDATA.Original_Points.C7;
+ SC = BLDATA.Original_Points.SC;
+ AC = BLDATA.Original_Points.AC;
+ AA = BLDATA.Original_Points.AA;
+ TS = BLDATA.Original_Points.TS;
+ AI = BLDATA.Original_Points.AI;
+ GH = BLDATA.Original_Points.GH;
+ HU = BLDATA.Original_Points.HU;
+ EL = BLDATA.Original_Points.EL;
+ EM = BLDATA.Original_Points.EM;
+ CP = BLDATA.Original_Points.CP;
+ US = BLDATA.Original_Points.US;
+ RS = BLDATA.Original_Points.RS;
+
+ % Set the Points
+ BLDATA.Initial_Points.IJ = IJ;
+ BLDATA.Initial_Points.PX = PX;
+ BLDATA.Initial_Points.T8 = T8;
+ BLDATA.Initial_Points.C7 = C7;
+ BLDATA.Initial_Points.SC = SC;
+ BLDATA.Initial_Points.AC = Rc*Rci'*(AC - SC) + BLDATA.Initial_Points.SC;% Rci'*(AC - SC) is AC in Clavicle coordinate, so with the rotation matrix Rc we can define its new position in the inertial frame as such.
+% plot the points to be sure they hold correspondence with their anatomical
+% landmarks
+Current_Folder = pwd; % get the current folder
+if isunix % load the scapula mesh that is in scapula coordinate system
+ load([Current_Folder, '/Scapula_Mesh0.mat']);
+else
+ load([Current_Folder, '\Scapula_Mesh0.mat']);
+end
+
+tri = Scapula_Mesh0.tri; % save the surfaces and the points separately
+points = 1.e3*Scapula_Mesh0.points; % scale the points to mm
+
+BLDATA = MAIN_INITIALISATION_build_data_bony_landmark(); % get the BLDATA in which we have R_S2T
+
+points_T = MAIN_TOOL_geometry_functions('Rotate Points From Local To Global Frame (Current)', points', BLDATA, 2); % rotate the points of the mesh to the thorax coordinate system