diff --git a/ShoulderCase/@Humerus/plot.m b/ShoulderCase/@Humerus/plot.m
new file mode 100644
index 0000000..fba7f96
--- /dev/null
+++ b/ShoulderCase/@Humerus/plot.m
@@ -0,0 +1,12 @@
+function output = plot(obj)
+  [x, y, z] = sphere;
+  x = obj.radius*x + obj.center(1);
+  y = obj.radius*y + obj.center(2);
+  z = obj.radius*z + obj.center(3);
+  plotHandle = surf(x, y, z,...
+    "FaceColor", "magenta",...
+    "FaceAlpha", 0.3,...
+    "EdgeColor", "magenta",...
+    "EdgeAlpha", 1);
+  output = plotHandle;
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/Muscle.m b/ShoulderCase/@Muscle/Muscle.m
index ef058b4..4132f85 100644
--- a/ShoulderCase/@Muscle/Muscle.m
+++ b/ShoulderCase/@Muscle/Muscle.m
@@ -1,85 +1,87 @@
 classdef Muscle < handle
 % The Muscle class is linked to a segmented muscle file (a mask)
 % and to the slice it has been segmented out of.
 %
 % Then, this class can measured values linked to the PCSA and the
 % muscle's degeneration.
 
   properties
-    name = '';
-    segmentationSet = '';
-    sliceName = '';
+    name = "";
+    segmentationSet = "";
+    sliceName = "";
     
-    PCSA = [];
-    atrophy = [];
-    fat = [];
-    osteochondroma = [];
-    degeneration = [];
+    PCSA = nan;
+    atrophy = nan;
+    fat = nan;
+    osteochondroma = nan;
+    degeneration = nan;
+
+    container = nan;
+    insertions = nan;
+    forceVector = nan;
   end
 
-
-
-  properties (Hidden = true)
-    dataPath = '';
-    maskDataPath = '';
-    contourDataPath = '';
-    container = [];
-  end
-
-
   methods
-
-    function obj = Muscle(musclesContainer,muscleName,slicesWithMuscleName)
+    function obj = Muscle(musclesContainer,muscleName)
       obj.container = musclesContainer;
       obj.name = muscleName;
-      obj.sliceName = slicesWithMuscleName;
-    end
 
-    function propagateDataPath(obj)
-      % Update current dataPath
-      % Propagate the path to objects in properties
+      [~, ~] = mkdir(obj.dataPath);
+      [~, ~] = mkdir(obj.dataMaskPath);
+      [~, ~] = mkdir(obj.dataContourPath);
+    end
 
-      obj.dataPath = fullfile(obj.container.dataPath,obj.name);
-      obj.maskDataPath = fullfile(obj.dataPath,'mask');
-      obj.contourDataPath = fullfile(obj.dataPath,'contour');
+    function output = dataPath(obj)
+      output = fullfile(obj.container.dataPath, obj.name);
+    end
 
-      if not(isfolder(obj.dataPath))
-        mkdir(obj.dataPath);
-      end
-      if not(isfolder(obj.maskDataPath))
-        mkdir(obj.maskDataPath);
-      end
-      if not(isfolder(obj.contourDataPath))
-        mkdir(obj.contourDataPath);
-      end
+    function output = dataMaskPath(obj)
+      output = fullfile(obj.dataPath, "mask");
     end
 
+    function output = dataContourPath(obj)
+      output = fullfile(obj.dataPath, "contour");
+    end
 
     function output = getValues(obj)
       output = cellfun(@(property) obj.(property), properties(obj),...
-                       'UniformOutput',false);
+        'UniformOutput',false);
     end
 
-    function measure(obj,segmentationSet)
+    function output = summary(obj)
+      summary.name = string(obj.name);
+      summary.segmentationSet = string(obj.segmentationSet);
+      summary.sliceName = string(obj.sliceName);
+      summary.PCSA = obj.PCSA;
+      summary.atrophy = obj.atrophy;
+      summary.fat = obj.fat;
+      summary.osteochondroma = obj.osteochondroma;
+      summary.degeneration = obj.degeneration;
+      output = struct2table(summary);
+    end
+
+    function measureDegeneration(obj,sliceName, segmentationSet)
+      obj.sliceName = sliceName;
       obj.segmentationSet = segmentationSet;
-      measurer = MuscleMeasurer(obj);
 
-      obj.PCSA = measurer.getPCSA;
-      obj.atrophy = measurer.getRatioAtrophy;
-      obj.fat = measurer.getRatioFat;
-      obj.osteochondroma = measurer.getRatioOsteochondroma;
-      obj.degeneration = measurer.getRatioDegeneration;
+      try
+        measurer = MuscleMeasurer(obj);
+        obj.PCSA = measurer.getPCSA;
+        obj.atrophy = measurer.getRatioAtrophy;
+        obj.fat = measurer.getRatioFat;
+        obj.osteochondroma = measurer.getRatioOsteochondroma;
+        obj.degeneration = measurer.getRatioDegeneration;
+      end
     end
 
     function output = loadMask(obj)
-      output = imread(fullfile(obj.maskDataPath,obj.segmentationSet + "_Mask.png")) > 0;
+      output = imread(fullfile(obj.dataMaskPath,obj.segmentationSet + "_Mask.png")) > 0;
     end
 
     function output = loadPixelCoordinates(obj)
-      output = load(fullfile(obj.container.slicesDataPath,...
+      output = load(fullfile(obj.container.dataSlicesPath,...
         obj.sliceName + "_PixelCoordinates.mat")).imagesPixelCoordinates;
     end
-
   end
 
 end
diff --git a/ShoulderCase/@Muscle/measureForceVector.m b/ShoulderCase/@Muscle/measureForceVector.m
new file mode 100644
index 0000000..434ab2b
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureForceVector.m
@@ -0,0 +1,5 @@
+function measureForceVector(obj)
+  force = Vector(obj.getCentroid().coordinates, obj.insertions);
+  force = force*obj.PCSA / force.norm;
+  obj.forceVector = force;
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureInsertionTangentToSphere.m b/ShoulderCase/@Muscle/measureInsertionTangentToSphere.m
new file mode 100644
index 0000000..e664d3e
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureInsertionTangentToSphere.m
@@ -0,0 +1,6 @@
+function measureInsertionTangentToSphere(obj, sphere, sliceNormal, polarPointOrientation)
+  obj.insertions = sphere.polarTangentPoint(...
+    obj.getCentroid().coordinates,...
+    sliceNormal,...
+    polarPointOrientation);
+end
\ No newline at end of file
diff --git a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
index 7efa616..762cab7 100644
--- a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
+++ b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
@@ -1,119 +1,119 @@
 classdef MuscleMeasurer < handle
 % Used to perfom PCSA and degeneration measurements of a Muscle object.
 % Slices and mask images are expected to be found at some specific places
 % and have specific names.
 %
 % This class is the results of extracting measurements methods from the
 % project of Nathan Donini.  
 
   properties
     muscle
     segmentedArea
     segmentedImage
 
     rangeHU = double([-1000 1000]);
     muscleThreshold = 0;
     fatThreshold = 30;
     osteochondromaThreshold = 166;
   end
 
 
 
   methods
 
     function obj = MuscleMeasurer(muscle)
       obj.muscle = muscle;
       obj.segmentedArea = obj.getSegmentedArea;
       obj.segmentedImage = obj.getSegmentedImage;
 
       obj.normalizeThresholds;
     end
 
     function output = getSegmentedArea(obj)
-      output = logical(imread(fullfile(obj.muscle.maskDataPath,[obj.muscle.segmentationSet '_Mask.png'])));
+      output = logical(imread(fullfile(obj.muscle.dataMaskPath,obj.muscle.segmentationSet + "_Mask.png")));
     end
 
     function output = getSegmentedImage(obj)
-      load(fullfile(obj.muscle.container.slicesDataPath,[obj.muscle.sliceName '_ForMeasurements.mat']));
+      load(fullfile(obj.muscle.container.dataSlicesPath,obj.muscle.sliceName + "_ForMeasurements.mat"));
       normalizedImage = obj.getNormalizedImage(imageForMeasurements);
 
       segmentedImage = normalizedImage.*obj.segmentedArea;
 
       output = segmentedImage;
     end
 
     function output = getNormalizedImage(obj,image)
       image(image < obj.rangeHU(1)) = obj.rangeHU(1);
       image(image > obj.rangeHU(2)) = obj.rangeHU(2);
       output = (double(image)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
     end
 
     function normalizeThresholds(obj)
       obj.muscleThreshold = ((obj.muscleThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
       obj.fatThreshold = ((obj.fatThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
       obj.osteochondromaThreshold = ((obj.osteochondromaThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
     end
 
     function output = getPCSA(obj)
-      load(fullfile(obj.muscle.container.slicesDataPath,[obj.muscle.sliceName '_PixelSpacings.mat']));
+      load(fullfile(obj.muscle.container.dataSlicesPath,obj.muscle.sliceName + "_PixelSpacings.mat"));
       pixelSurface = (imagesPixelSpacings(1) * imagesPixelSpacings(2)) / 100; % cm^2
       output = pixelSurface * bwarea(obj.segmentedArea);
     end
 
     function output = getRatioAtrophy(obj)
       output = obj.getRatioAreas(obj.getAreaAtrophy,obj.segmentedArea);
     end
 
     function output = getRatioFat(obj)
       output = obj.getRatioAreas(obj.getAreaFat,obj.segmentedArea);
     end
 
     function output = getRatioOsteochondroma(obj)
       output = obj.getRatioAreas(obj.getAreaOsteochondroma,obj.segmentedArea);
     end
 
     function output = getRatioDegeneration(obj)
       output = obj.getRatioAtrophy + obj.getRatioFat + obj.getRatioOsteochondroma;
     end
 
     function output = getRatioAreas(obj,partialArea,totalArea)
       output = bwarea(partialArea)/bwarea(totalArea);
     end
 
     function output = getAreaMuscle(obj)
       % The area "Muscle" is the area of the segmented image that is not atrophied.
       % Looking for a better name.
       areaMuscle = imbinarize(obj.segmentedImage,obj.muscleThreshold);
       areaMuscle = imfill(areaMuscle,'holes');
       % Keep biggest island only
       islands = bwconncomp(areaMuscle);
       islandsSizes = cellfun(@numel,islands.PixelIdxList);
       [~,biggestIslandIndex] = max(islandsSizes);
       areaMuscle = false(size(areaMuscle));
       areaMuscle(islands.PixelIdxList{biggestIslandIndex}) = true;
 
       output = areaMuscle;
     end
 
     function output = getAreaAtrophy(obj)
       output = obj.segmentedArea & not(obj.getAreaMuscle);
     end
 
     function output = getAreaFat(obj)
       muscleImage = obj.segmentedImage.*obj.getAreaMuscle;
       areaFat = obj.getAreaMuscle & not(imbinarize(muscleImage,obj.fatThreshold));
 
       output = areaFat;
     end
 
     function output = getAreaOsteochondroma(obj)
       muscleImage = obj.segmentedImage.*obj.getAreaMuscle;
       areaOsteochondroma = imbinarize(muscleImage,obj.osteochondromaThreshold);
       areaOsteochondroma = imfill(areaOsteochondroma,'holes');
 
       output = areaOsteochondroma;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@MusclesContainer/MusclesContainer.m b/ShoulderCase/@MusclesContainer/MusclesContainer.m
index 3cbdc73..130e173 100644
--- a/ShoulderCase/@MusclesContainer/MusclesContainer.m
+++ b/ShoulderCase/@MusclesContainer/MusclesContainer.m
@@ -1,122 +1,103 @@
 classdef MusclesContainer < handle
 % Replaces the deprecated Muscles class.
 %
 % One key features of MusclesContainer is to create the dicom
 % slices and to call the automatic rotator cuff segmentation
 % on them. These methods might be extracted at some point.
 %
 % The plot() methods create a figure with the rotator cuff slice
 % with the segmented muscles contour superposed.
 %
 % Also handle a list of Muscle instances a display a summary of
 % these muscles measurements.
 %
 % New muscles can easily be added and measured through this class,
 % however their slicing/segmentation will remain to be implemented.
 
   properties
     summary = [];
     list = {};
-  end
-
-
-
-  properties (Hidden = true)
-    dataPath = '';
-    slicesDataPath
     shoulder
   end
 
-
-
   methods
 
     function obj = MusclesContainer(shoulder)
       obj.shoulder = shoulder;
       obj.clearSummary();
-    end
-
-    function propagateDataPath(obj)
-      % Update current dataPath
-      % Propagate the path to objects in properties
 
-      obj.dataPath = fullfile(obj.shoulder.dataPath,'muscles');
-      obj.slicesDataPath = fullfile(obj.dataPath,'oblique_slices');
+      [~, ~] = mkdir(obj.dataPath);
+      [~, ~] = mkdir(obj.dataSlicesPath);
+    end
 
-      if not(isfolder(obj.dataPath))
-        mkdir(obj.dataPath);
-      end
-      if not(isfolder(obj.slicesDataPath))
-        mkdir(obj.slicesDataPath);
-      end
+    function output = dataPath(obj)
+      output = fullfile(obj.shoulder.dataPath, "muscles");
+    end
 
-      if obj.isEmpty
-        return 
-      end
-      
-      structfun(@(muscle) muscle.propagateDataPath, obj.list);
+    function output = dataSlicesPath(obj)
+      output = fullfile(obj.dataPath, "oblique_slices");
     end
 
     function addMuscle(obj,muscleName,sliceWithMuscleFilename)
       newMuscle = Muscle(obj,muscleName,sliceWithMuscleFilename);
       newMuscle.propagateDataPath;
       obj.list.(muscleName) = newMuscle;
 
       obj.updateSummary;
     end
 
     function updateSummary(obj)
       muscleValues = obj.getAllMusclesValues();
       if isempty(muscleValues)
         obj.clearSummary();
         return
       end
 
       muscleProperties = properties('Muscle');
       obj.summary = cell2table(muscleValues(:,2:end),...
                                'VariableNames',muscleProperties(2:end),...
                                'RowNames',muscleValues(:,1));
     end
     
     function clearSummary(obj)
       muscleProperties = properties('Muscle');
       obj.summary = array2table(zeros(1,length(muscleProperties)-1),...
                                 'VariableNames',muscleProperties(2:end));
     end
 
     function output = getAllMusclesValues(obj)
       output = {};
       if obj.isEmpty
         return
       end
 
       muscleNames = fields(obj.list);
       for i = 1:length(muscleNames)
         output(end+1,:) = obj.list.(muscleNames{i}).getValues';
       end
     end
 
     function measureAllMuscles(obj,sliceName,maskName)
       if obj.isEmpty
         return 
       end
 
       muscleNames = fields(obj.list);
       for i = 1:length(fields(obj.list))
         try
           obj.addMuscle(muscleNames{i},sliceName); % update muscle with requested slice
           obj.list.(muscleNames{i}).measure(maskName);
         catch ME
           warning(ME.message);
         end
       end
       obj.updateSummary;
     end
 
     function output = isEmpty(obj)
       output = isempty(obj.list);
     end 
 
   end
 
 end
diff --git a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m
index 98945d5..1efce3f 100644
--- a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m
+++ b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m
@@ -1,97 +1,97 @@
 function createRotatorCuffSlicesMatthieu(obj,sliceName)
   % This function and its results might be reffered as 
   % Matthieu algorithm and Matthieu slices in some studies.
   %
   % This method has been written to simplify, improve, and make reusable 
   % the Nathan's slicing algorithm written in createRotatorCuffSlicesNathan(). 
   % This means:
   %     - taking advantage of the MATLAB obliqueslice() method wich wasn't available
   %       when Nathan's algorithm has been written.
   %     - encapsulating the slicing data and methods in a DicomVolumeSlicer object.
   %     - being able to orient a slice according to dicom points coordinates.
   %     - being able to crop the slice according to dicom points coordinates
   %       and cropping dimensions given in the dicom set units. 
   %
   % Images and pixel spacings are saved at:
   % shoulder/dataDev/*/*/*/SCase-IPP/CT-SCase-IPP/matlab/shoulder/muscles/oblique_slice/
 
   assert(not(obj.shoulder.scapula.coordSys.isEmpty),'Scapula coordinate system not measured.');
   [imageForSegmentation,imageForMeasurements,imagesPixelSpacings,imagesPixelCoordinates] = getSlices(obj);
   saveImages(obj,imageForSegmentation,imageForMeasurements,sliceName);
   saveImagesPixelSpacings(obj,imagesPixelSpacings,sliceName);
   saveImagesPixelCoordinates(obj,imagesPixelCoordinates,sliceName);
 end
 
 
 
 function [slice8Bit,slice16Bit,slicePixelSpacings,slicePixelCoordinates] = getSlices(obj)
   SCase = obj.shoulder.SCase;
   scapula = obj.shoulder.scapula;
   path = SCase.dataDicomPath;
   try
     rotatorCuffSlicer = DicomVolumeSlicer(SCase.getSmoothDicomPath);
   catch
     rotatorCuffSlicer = DicomVolumeSlicer(SCase.dataDicomPath);
   end
 
   % The 3D volume must be normalised for the slices not to be distorted
   rotatorCuffSlicer.setMinimalResolutionDividedBy(2);
   rotatorCuffSlicer.normaliseVolume;
   rotatorCuffSlicer.slice(scapula.coordSys.origin,scapula.coordSys.ML);
 
   % The following lines set the slice to display a 300x300 mm area around the
   % scapula with the "Y" shape of the scapula oriented up
   rotatorCuffSlicer.orientSliceUpwardVector(scapula.coordSys.origin,scapula.coordSys.IS);
   rotatorCuffSlicer.crop(scapula.coordSys.origin-100*scapula.coordSys.IS,300,300);
 
   % Rotator cuff segmentation code requires 1024x1024 pixels images
   if size(rotatorCuffSlicer.sliced,1) > size(rotatorCuffSlicer.sliced,2)
     rotatorCuffSlicer.resize([1024 nan]);
   else
     rotatorCuffSlicer.resize([nan 1024]);
   end
   rotatorCuffSlicer.addEmptyBackgroundToSlice(1024,1024);
 
 
   rawSlice = rotatorCuffSlicer.sliced;
 
   rotatorCuffSlicer.rescaleSliceToUint8;
   slice8Bit = rotatorCuffSlicer.sliced;
   rotatorCuffSlicer.sliced = rawSlice;
 
   rotatorCuffSlicer.rescaleSliceToInt16;
   slice16Bit = rotatorCuffSlicer.sliced;
 
   slicePixelCoordinates = rotatorCuffSlicer.getSlicedPixelCoordinates();
 
   % Flip the slices for right and left shoulders to give similar slices
   if obj.shoulder.side=='L'
     slice8Bit=rot90(flipud(slice8Bit),2);
     slice16Bit=rot90(flipud(slice16Bit),2);
     slicePixelCoordinates=rot90(flipud(slicePixelCoordinates),2);
   end
 
   slicePixelSpacings = rotatorCuffSlicer.slicedPixelSpacings;
 end
 
 
 
 function saveImages(obj,imageForSegmentation,imageForMeasurements,sliceName)
   imwrite(imageForSegmentation,...
-    fullfile(obj.slicesDataPath,sliceName + "_ForSegmentation.png"));
-  save(fullfile(obj.slicesDataPath,sliceName + "_ForMeasurements.mat"),...
+    fullfile(obj.dataSlicesPath,sliceName + "_ForSegmentation.png"));
+  save(fullfile(obj.dataSlicesPath,sliceName + "_ForMeasurements.mat"),...
     'imageForMeasurements');
 end
 
 
 
 function saveImagesPixelSpacings(obj,imagesPixelSpacings,sliceName)
-  save(fullfile(obj.slicesDataPath,sliceName + "_PixelSpacings.mat"),...
+  save(fullfile(obj.dataSlicesPath,sliceName + "_PixelSpacings.mat"),...
     'imagesPixelSpacings');
 end
 
 
 function saveImagesPixelCoordinates(obj,imagesPixelCoordinates,sliceName)
-  save(fullfile(obj.slicesDataPath,sliceName + "_PixelCoordinates.mat"),...
+  save(fullfile(obj.dataSlicesPath,sliceName + "_PixelCoordinates.mat"),...
     'imagesPixelCoordinates');
 end
\ No newline at end of file
diff --git a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
index 2fd2d24..94052f2 100644
--- a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
+++ b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
@@ -1,368 +1,368 @@
 function createRotatorCuffSlicesNathan(obj,sliceName)
   % This function and its results might be reffered as 
   % Nathan algorithm and Nathan slices in some studies.
   %
   %One needs to create two empty files in the current directory :
   %data_degenerescence.mat (a table) and Didnotwork.mat (a
   %string)
   %
   % This script calculates the degeneration of the 4 muscles
   % of the rotator cuff : SS, SC, IS, TM. It is composed of four
   % different steps :
   % - Load CT images (DICOM files) in order to create a 3D matrix
   % consituted with HU (Hounsfield unit) values
   % - Create a slice of this 3D matrix along a patient-specific
   % plane. Do image processing and set saving format usable by
   % automatic segmentation code (Python's code from ARTORG team
   % (UNIBE))
   % - Execute automatic segmentation code on the new created
   % images
   % - Apply degeneration code (from Laboratoire de biom�canique
   % orthop�dique (EPFL)) on new images with automatic
   % segmentation and on image with manual segmentation (from CHUV
   % radiologist). Save degeneration values under the "muscles"
   % properties of each case in the LBO database
   %
   % Some "try-catch" are present in order to catch the errors
   % while keeping running the script on the rest of the database.
   % They are summarized underneath :
   % (DicomReadVolume) : do not manage to load the dicoms
   % Nothing : Automatic segmentation does not work (generally for TM)
   % * : slice function does not work
   % ** : does not want to crop image
   % *** : black image, don't want to save png
 
   % Main code starts here !
 
   %% Load dicom
   SCase = obj.shoulder.SCase;
   dataCTPath = SCase.dataCTPath;
 
   % Create "muscles" folder (and its 4 subfolders) if it does not exist
   folder_muscle=[dataCTPath '/muscles'];
 
   if isfolder(folder_muscle)
   else
   mkdir(dataCTPath,'muscles');
   mkdir(folder_muscle,'IS');
   mkdir(folder_muscle,'SS');
   mkdir(folder_muscle,'SC');
   mkdir(folder_muscle,'TM');
   end
 
   % Use folder "CT-2" (smoothed kernel) if it exists
   dataCTPath2=dataCTPath;
   dataCTPath2(1,end)='2';
 
   if isfolder(dataCTPath2)
       dataCTPath0=dataCTPath2;
   else
       dataCTPath0=dataCTPath;
   end
 
   dicomFolder = [dataCTPath0 '/dicom'];
 
   %Extract the images as a 4D matrix where spatial is the
   %positions of the Pixels in the CT coordinate system
 
   % Return the case's number if it does not manage to load the DICOM files
   try
       [V,spatial,dim] = dicomreadVolume(dicomFolder);
   catch
       outputArg=1;
       return;
   end
 
   %V needs to be 3D matrix containing the value in HU of the
   %corresponding voxel
   V = squeeze(V) ;
 
   %Read the information from the middle dicom image (all image
   %informations are the same for all dicom files) to avoid a problem
   %if other files are located at the beginning or end of the
   %folder
   dicomList = dir(dicomFolder) ;
   size_list = size(dicomList) ;
   dicom_information = dicominfo([dicomFolder '/' dicomList(round(size_list(1)/2)).name]);
 
   %Take into account that the data will need to be linearly transformed
   %from stored space to memory to come back to true HU units afterwards:
   Rescale_slope = dicom_information.RescaleSlope;
   Rescale_intercept = dicom_information.RescaleIntercept;
 
   %Convert to double type to be able to multiply it with a int64
   %matrix afterwards
   V = double(V);
 
   % To filter hard kernel images
 
   %             if isfolder(dataCTPath2)
   %             else
   %             V = imgaussfilt3(V,1) ;
   %             end
 
   size_V = size(V);
 
   %Extract position of each first pixel of each image (upper left corner)
   PatientPositions = spatial.PatientPositions ;
 
   %Get the distance between each pixel in neighboring rows and
   %columns
   PixelSpacings = spatial.PixelSpacings;
 
 
   % Get local coordinate system from parent scapula object(glenoid center)
   origin = obj.shoulder.scapula.coordSys.origin; % do not mix up with origin_image, "origin" is the origin of the Scapula Coordinate System
   xAxis = obj.shoulder.scapula.coordSys.PA;
   yAxis = obj.shoulder.scapula.coordSys.IS;
   zAxis = obj.shoulder.scapula.coordSys.ML;
 
 
   %% Slice from the CT data set
 
   %Get the position of the upper left pixel of the first image
   %and calculate the position of each pixel in the x,y
   %and z direction starting from the "origin" pixel
   origin_image = PatientPositions(1,:); % coordinates of the first pixel
   x_max = origin_image(1)+PixelSpacings(1,1)*(size_V(1)-1);
   y_max = origin_image(2)+PixelSpacings(1,2)*(size_V(2)-1);
   z_max = PatientPositions(end,3);
 
   %Calculate the coefficient of the linear transformation (CT-scan
   %coordinate system (x,y,z) to images coordinates system (i,j,k))
   coefficients_i = polyfit([1 size_V(1)],[origin_image(1) , x_max],1);
   a_i = coefficients_i(1);
   b_i = coefficients_i(2);
 
   coefficients_j = polyfit([1 size_V(2)],[origin_image(2) , y_max],1);
   a_j = coefficients_j(1);
   b_j = coefficients_j(2);
 
   coefficients_k = polyfit([1 size_V(3)],[origin_image(3) , z_max],1);
   a_k = coefficients_k(1);
   b_k = coefficients_k(2);
 
   Pixel_pos_x=(1:size_V(1))*a_i+b_i; % Create the position vector along x-axis
   Pixel_pos_y=(1:size_V(2))*a_j+b_j;
   Pixel_pos_z=(1:size_V(3))*a_k+b_k;
 
   PixSize_z=(z_max-origin_image(3))/(size_V(3)-1); % Space between 2 slices
 
   % The oblique slice to be created is perpendicular to the z
   % axis of the scapula CS and passes through the origin point.
   % We will apply the matlab function "slice", to create the
   % oblique slice. One must first create a slice parallel to 2
   % system's axis. Here we take a plane parallel to the (x,y)
   % plane, passing through the point "origin". Then, we need to
   % rotate this plane so that it becomes normal to the zAxis. V1
   % is the vector normal to the horizontal slice and v2 the
   % vector normal to the final oblique slice, parallel to the
   % zAxis. We need to rotate around a vector v3 and from an angle
   % given by "vrrotvec(v1,v2)".
   v1=[0,0,1];
   v2=zAxis/norm(zAxis);
   rotation_info=vrrotvec(v1,v2);
   v3=rotation_info(1:3);
   theta=rotation_info(4)/pi*180; % Convert it in degrees
 
   % Create the plane to be rotated
   [x,y,z] = meshgrid(Pixel_pos_x,Pixel_pos_y,Pixel_pos_z);
   s=1000; % Extension factor of the slice in order to cut the whole matrix.Set arbitrarily but sufficiently big so that it works for every patient
   Pixel_pos_x_enlarged=[origin_image(1)+(-PixelSpacings(1,1)*s:PixelSpacings(1,1):-PixelSpacings(1,1)) Pixel_pos_x x_max+(PixelSpacings(1,1):PixelSpacings(1,1):PixelSpacings(1,1)*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged=[origin_image(2)+(-PixelSpacings(1,2)*s:PixelSpacings(1,2):-PixelSpacings(1,2)) Pixel_pos_y y_max+(PixelSpacings(1,2):PixelSpacings(1,2):PixelSpacings(1,2)*s)];
 
   figure('visible','off');
   hsp = surf(Pixel_pos_x_enlarged,Pixel_pos_y_enlarged,origin(3)*ones(length(Pixel_pos_y_enlarged),length((Pixel_pos_x_enlarged))));
 
   rotate(hsp,v3,theta,origin);
 
   xd_0 = hsp.XData; % x-coordinates of the plane's points
   yd_0 = hsp.YData; % y-coordinates of the plane's points
   zd_0 = hsp.ZData; % z-coordinates of the plane's points
 
   % Slice the 3D matrix of CT-scan
   try % Return case's number in case slice does not work
   figure('visible','off');
   hss=slice(x,y,z,V,xd_0,yd_0,zd_0,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
   catch
   outputArg=2;
   return
   end
 
   HU_0 = get(hss,'CData'); % HU values at each point of the plane
 
   NaNplaces_0=find(isnan(HU_0));
 
   HU_1=HU_0*Rescale_slope+Rescale_intercept; % getting true HU values
   p_0=10e6; % value outside 3D matrix, set to a high distinguishable number
   HU_1(NaNplaces_0)=p_0;
 
   % If it's a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_1=flipud(HU_1');
   else
   HU_1=HU_1';
   end
 
   % Find the new pixel's sizes necessary to obtain a 1024*1024
   % shoulder slice. The largest diagonal's length is set to 1024
   % pixels
   [left_1,left_2]=find(HU_1~=p_0,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_1,2)~=p_0,1);
   right_r=right_1;
   right_c=size(HU_1,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_1)~=p_0,1);
   up_c=size(HU_1)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_1,3)~=p_0,1);
   down_r=size(HU_1,1)-down_2+1;
   down_c=down_1;
 
   % In case the slice does not cut the whole matrix (which
   % normally never happens), ask to double "s" parameter
 
   if or(or(right_c==size(HU_1,2),left_c==1),or(up_r==1,down_r==size(HU_1,1)))
     outputArg = 3;
     return
   end
 
   % Find which diagonal of the matrix' effective slice is the largest
   w=right_c-left_c; % width
   h=down_r-up_r; % height
   if h>w % Set the largest diagonal at 1024 pixels
       l=h;
   else
       l=w;
   end
 
   % Define new pixel's size
   NewPixSize_x=PixelSpacings(1,1)*l/1023;
   NewPixSize_y=PixelSpacings(1,2)*l/1023;
 
   % New plane with the new pixel's sizes
   Pixel_pos_x_enlarged_2=[origin_image(1)+(-NewPixSize_x*s:NewPixSize_x:-NewPixSize_x) origin_image(1):NewPixSize_x:x_max x_max+(NewPixSize_x:NewPixSize_x:NewPixSize_x*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged_2=[origin_image(2)+(-NewPixSize_y*s:NewPixSize_y:-NewPixSize_y) origin_image(2):NewPixSize_y:y_max y_max+(NewPixSize_y:NewPixSize_y:NewPixSize_y*s)];
 
   figure('visible','off');
   hsp2 = surf(Pixel_pos_x_enlarged_2,Pixel_pos_y_enlarged_2,origin(3)*ones(length(Pixel_pos_y_enlarged_2),length((Pixel_pos_x_enlarged_2))));
   rotate(hsp2,v3,theta,origin);
 
   xd = hsp2.XData; % x-coordinates of the rotated plane's points.
   yd = hsp2.YData; % y-coordinates of the rotated plane's points.
   zd = hsp2.ZData; % z-coordinates of the rotated plane's points.
 
   % New slice with the good resolution (1024*1024)
   figure('visible','off');
   hss=slice(x,y,z,V,xd,yd,zd,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
 
   HU_2 = get(hss,'CData'); % HU values at each point of the plane
   NaNplaces=find(isnan(HU_2));
   HU_slice=HU_2*Rescale_slope+Rescale_intercept;
   HU_slice_8_bit=HU_slice;
   p_tif=0; % value outside 3D matrix, for the tif image
   p_png=85; % for the png image
   HU_slice(NaNplaces)=p_tif;
 
   % Scaling to fit the training set (for the png only, grayscale)
   HU_interval_png=[-100 160];
   lin_coef=[HU_interval_png(1) 1 ; HU_interval_png(2) 1]^(-1)*[0 255]';
 
   HU_slice_8_bit=lin_coef(1)*HU_slice_8_bit+lin_coef(2);
   HU_slice_8_bit(NaNplaces)=p_png;
 
   % If it a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_slice=flipud(HU_slice');
   HU_slice_8_bit=flipud(HU_slice_8_bit');
   else
   HU_slice=HU_slice';
   HU_slice_8_bit=HU_slice_8_bit';
   end
 
   % Determine the location of the 4 corners of the slice
   [left_1,left_2]=find(HU_slice~=p_tif,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_slice,2)~=p_tif,1);
   right_r=right_1;
   right_c=size(HU_slice,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_slice)~=p_tif,1);
   up_c=size(HU_slice,2)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_slice,3)~=p_tif,1);
   down_r=size(HU_slice,1)-down_2+1;
   down_c=down_1;
 
   Folder=[dataCTPath '/matlab'];
 
   % Crop the image
   try % If cropping does not work, return case's number
   HU_slice=HU_slice(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice_8_bit=HU_slice_8_bit(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice=int16(HU_slice); % convert into 16-bit
   catch
   outputArg=4;
   return
   end
 
   % Create '.tif' file
   file_name_tif=[SCase.id '.tif'];
 
   t = Tiff(fullfile(Folder,file_name_tif),'w');
   setTag(t,'Photometric',Tiff.Photometric.MinIsBlack);
   setTag(t,'Compression',Tiff.Compression.None);
   setTag(t,'BitsPerSample',16);
   setTag(t,'SamplesPerPixel',1);
   setTag(t,'ImageLength',1024);
   setTag(t,'ImageWidth',1024);
   setTag(t,'PlanarConfiguration',Tiff.PlanarConfiguration.Chunky);
   setTag(t,'SampleFormat',Tiff.SampleFormat.Int);
   write(t,HU_slice);
   close(t);
 
   % Create 8 bit '.png' file and save it at the right places
   % (matlab folder and python's input folder)
   file_name_png=[SCase.id '.png'];
 
   HU_slice_8_bit=uint8(HU_slice_8_bit);
   figure('visible','off');
   image(HU_slice_8_bit);
   colormap(gray(256));
 
   try % return case's number if does not manage to write png in the "matlab" folder because it is a black image
   imwrite(HU_slice_8_bit,[Folder '/' file_name_png]);
   catch
   outputArg=5;
   return
   end
 
   pythonDir = ConfigFileExtractor.getVariable('pythonDir');
   Folder_input_py= [pythonDir '/input'];
   imwrite(HU_slice_8_bit,[Folder_input_py '/' file_name_png]); % write png file in the python input folder
   % imwrite(HU_slice_8_bit,['./rcseg/' file_name_png]); % write png file in local /rceg folder
 
   % Create .txt file with pixel dimensions
   file_name_info=[SCase.id '.info'];
   text1='# Amira Stacked Slices' ;
   text2=['pixelsize ', num2str(NewPixSize_x),' ',num2str(NewPixSize_y)];
   fileID = fopen([Folder '/' file_name_info],'w');
   fprintf(fileID,'%s\r\n\r\n',text1);
   fprintf(fileID,'%s\n',text2);
   fclose(fileID);
   close all;
 
 
-  imwrite(HU_slice_8_bit,fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']));
+  imwrite(HU_slice_8_bit,fullfile(obj.dataSlicesPath,[sliceName '_ForSegmentation.png']));
   imageForMeasurements = HU_slice;
-  save(fullfile(obj.slicesDataPath,[sliceName '_ForMeasurements.mat']),'imageForMeasurements');
+  save(fullfile(obj.dataSlicesPath,[sliceName '_ForMeasurements.mat']),'imageForMeasurements');
   imagesPixelSpacings = [NewPixSize_x NewPixSize_y];
-  save(fullfile(obj.slicesDataPath,[sliceName '_PixelSpacings.mat']),'imagesPixelSpacings');
+  save(fullfile(obj.dataSlicesPath,[sliceName '_PixelSpacings.mat']),'imagesPixelSpacings');
 end
diff --git a/ShoulderCase/@MusclesContainer/plot.m b/ShoulderCase/@MusclesContainer/plot.m
index 01485e5..a0fb6c1 100644
--- a/ShoulderCase/@MusclesContainer/plot.m
+++ b/ShoulderCase/@MusclesContainer/plot.m
@@ -1,39 +1,39 @@
 function output = plot(obj,varargin)
   % For now this function is specific to plotting the rotator cuff segmentation
   % results
 
   if nargin == 2
     error('Provide the slice name and the mask name or provide nothing (auto slice and segmentations will be plotted).')
   elseif nargin == 3
     sliceName = [varargin{1} '_ForSegmentation.png'];
     maskName = [varargin{2} '_Mask.png'];
   else
     sliceName = 'rotatorCuffMatthieu_ForSegmentation.png';
     maskName = 'autoMatthieu_Mask.png';
   end
 
   SCase = obj.shoulder.SCase;
   muscleNames = fields(obj.list);
 
-  rotatorCuffCrossSection = imread(fullfile(obj.slicesDataPath,sliceName));
+  rotatorCuffCrossSection = imread(fullfile(obj.dataSlicesPath,sliceName));
 
   leg = {};
   plotHandle(1) = imshow(rotatorCuffCrossSection);
   hold on
   colors = {'g','r','b','y'};
 
   for i = 1:length(muscleNames)
     try
-      segmentedImage = imread(fullfile(obj.list.(muscleNames{i}).maskDataPath,maskName));
+      segmentedImage = imread(fullfile(obj.list.(muscleNames{i}).dataMaskPath,maskName));
     catch
       continue
     end
 
     if (max(segmentedImage,[],'all') > 0)   % if segmentation result is not empty
       leg{end+1} = muscleNames{i};
       [~,plotHandle(i+1)] = contour(segmentedImage,'color',colors{i});
     end
   end
   plotHandle(i+2) = legend(leg);
   output = plotHandle;
 end
diff --git a/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m b/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m
index 587d386..dcfddca 100644
--- a/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m
+++ b/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m
@@ -1,72 +1,72 @@
 function segmentRotatorCuffMuscles(obj,sliceName,maskName)
   % sliceName is the part before the '_ForSegmentation.png' part of the name of
   % the file that is sent to rcseg.
   % maskName is the part before the '_Mask.png' part of the name of
   % the file that is saved in the SCase folder.
   cleanRotatorCuffSegmentationWorkspace(obj);
   sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName);
   callRotatorCuffSegmentation(obj);
 
   % (Re)initialize the muscle objects to refer to the new segmentation results
   obj.addMuscle('IS',sliceName);
   obj.addMuscle('SS',sliceName);
   obj.addMuscle('SC',sliceName);
   obj.addMuscle('TM',sliceName);
 
   saveSegmentationResults(obj,maskName);
   cleanRotatorCuffSegmentationWorkspace(obj);
 end
 
 
 
 function cleanRotatorCuffSegmentationWorkspace(obj)
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
   delete(fullfile(rotatorCuffSegmentationPath,'input','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'IS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SC','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'TM','*'));
 end
 
 
 
 function sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName)
   SCase = obj.shoulder.SCase;
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
-  imageForSegmentationPath = fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']);
+  imageForSegmentationPath = fullfile(obj.dataSlicesPath,[sliceName '_ForSegmentation.png']);
   copyfile(imageForSegmentationPath,fullfile(rotatorCuffSegmentationPath,'input',[SCase.id '.png']));
 end
 
 
 
 function callRotatorCuffSegmentation(obj)
   %% Apply UNIBE method to automatically segment normal muscle
 
   pythonDir = ConfigFileExtractor.getVariable('pythonDir');
 
   pythonCommandActivateEnvironment= ['source ' fullfile(pythonDir,'venv','bin','activate') ';'];
   pythonCommandMoveToSegmentationWorkspace = ['cd ' pythonDir ';'];
   pythonCommandExecuteSegmentation = [fullfile(pythonDir,'rcseg.py') ' segment input' ';'];
 
   [~, pythonCommandOutput] = system([pythonCommandActivateEnvironment,...
                                      pythonCommandMoveToSegmentationWorkspace,...
                                      pythonCommandExecuteSegmentation]);
   disp(pythonCommandOutput);
 end
 
 
 
 function saveSegmentationResults(obj,maskName)
   SCaseID = obj.shoulder.SCase.id;
   muscleNames = fields(obj.list);
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
 
   for i = 1:length(muscleNames)
     try
       copyfile(fullfile(rotatorCuffSegmentationPath,muscleNames{i},[SCaseID '.png']),...
-               fullfile(obj.list.(muscleNames{i}).maskDataPath,[maskName '_Mask.png']) );
+               fullfile(obj.list.(muscleNames{i}).dataMaskPath,[maskName '_Mask.png']) );
     catch ME
       warning(ME.message)
     end
   end
 end
diff --git a/ShoulderCase/@RotatorCuff/RotatorCuff.m b/ShoulderCase/@RotatorCuff/RotatorCuff.m
new file mode 100644
index 0000000..dbbb6ce
--- /dev/null
+++ b/ShoulderCase/@RotatorCuff/RotatorCuff.m
@@ -0,0 +1,100 @@
+classdef RotatorCuff < handle
+
+  properties
+    SC
+    SS
+    IS
+    TM
+    anteroPosteriorImbalance = nan;
+
+    shoulder;
+  end
+
+  methods
+    function obj = RotatorCuff(shoulder)
+      obj.shoulder = shoulder;
+
+      obj.SC = Muscle(obj, "SC");
+      obj.SS = Muscle(obj, "SS");
+      obj.IS = Muscle(obj, "IS");
+      obj.TM = Muscle(obj, "TM");
+
+      [~, ~] = mkdir(obj.dataPath);
+      [~, ~] = mkdir(obj.dataSlicesPath);
+    end
+
+    function output = dataPath(obj)
+      output = fullfile(obj.shoulder.dataPath, "muscles");
+    end
+
+    function output = dataSlicesPath(obj)
+      output = fullfile(obj.dataPath, "oblique_slices");
+    end
+
+    function output = summary(obj)
+      summary = [];
+      summary = [summary; obj.SC.summary];
+      summary = [summary; obj.SS.summary];
+      summary = [summary; obj.IS.summary];
+      summary = [summary; obj.TM.summary];
+      output = summary;
+    end
+
+    function output = measure(obj, sliceName, segmentationSet)
+      try
+        obj.measureDegenerations(sliceName, segmentationSet);
+        obj.measureInsertions();
+        obj.measureAnteroPosteriorImbalance();
+      end
+      output = obj.summary;
+    end
+
+    function measureDegenerations(obj, sliceName, segmentationSet)
+      for muscleName = ["SC", "SS", "IS", "TM"]
+        obj.(muscleName).measureDegeneration(sliceName, segmentationSet);
+      end
+    end
+
+    function measureInsertions(obj)
+      humerus = obj.shoulder.humerus;
+      humeralHead = Sphere(humerus.center, humerus.radius);
+
+      % IS humeral insertion estimation
+      obj.IS.measureInsertionTangentToSphere(...
+        humeralHead,...
+        obj.shoulder.scapula.coordSys.IS,...
+        -obj.shoulder.scapula.coordSys.PA);
+      obj.IS.measureForceVector();
+
+      % SC humeral insertion estimation
+      obj.SC.measureInsertionTangentToSphere(...
+        humeralHead,...
+        obj.shoulder.scapula.coordSys.IS,...
+        obj.shoulder.scapula.coordSys.PA);
+      obj.SC.measureForceVector();
+    end
+
+    function measureAnteroPosteriorImbalance(obj)
+      forceSC = obj.SC.forceVector;
+      forceIS = obj.IS.forceVector;
+      scapula = obj.shoulder.scapula;
+      inferoSuperior = Vector(scapula.coordSys.origin,... 
+        scapula.coordSys.origin + scapula.coordSys.IS);
+
+      % Take infero-superior orthogonal component <=> project on scapular axial plane
+      forceSCAxial = inferoSuperior.orthogonalComplementToPoint(...
+        inferoSuperior.origin + forceSC.vector);
+      forceISAxial = inferoSuperior.orthogonalComplementToPoint(...
+        inferoSuperior.origin + forceIS.vector);
+
+      medioLateral = Vector(scapula.coordSys.origin,...
+        scapula.coordSys.origin + scapula.coordSys.ML);
+      posteroAnterior = Vector(scapula.coordSys.origin,...
+        scapula.coordSys.origin + scapula.coordSys.PA);
+      imbalanceVector = forceSCAxial + forceISAxial;
+      obj.anteroPosteriorImbalance = ...
+        sign(dot(imbalanceVector, -posteroAnterior)) * rad2deg(medioLateral.angle(imbalanceVector));
+    end
+  end
+
+end
\ No newline at end of file
diff --git a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m b/ShoulderCase/@RotatorCuff/createSliceMatthieu.m
similarity index 86%
copy from ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m
copy to ShoulderCase/@RotatorCuff/createSliceMatthieu.m
index 98945d5..201c716 100644
--- a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesMatthieu.m
+++ b/ShoulderCase/@RotatorCuff/createSliceMatthieu.m
@@ -1,97 +1,98 @@
-function createRotatorCuffSlicesMatthieu(obj,sliceName)
+function createSliceMatthieu(obj)
   % This function and its results might be reffered as 
   % Matthieu algorithm and Matthieu slices in some studies.
   %
   % This method has been written to simplify, improve, and make reusable 
-  % the Nathan's slicing algorithm written in createRotatorCuffSlicesNathan(). 
+  % the Nathan's slicing algorithm written in createSlicesNathan(). 
   % This means:
   %     - taking advantage of the MATLAB obliqueslice() method wich wasn't available
   %       when Nathan's algorithm has been written.
   %     - encapsulating the slicing data and methods in a DicomVolumeSlicer object.
   %     - being able to orient a slice according to dicom points coordinates.
   %     - being able to crop the slice according to dicom points coordinates
   %       and cropping dimensions given in the dicom set units. 
   %
   % Images and pixel spacings are saved at:
   % shoulder/dataDev/*/*/*/SCase-IPP/CT-SCase-IPP/matlab/shoulder/muscles/oblique_slice/
 
   assert(not(obj.shoulder.scapula.coordSys.isEmpty),'Scapula coordinate system not measured.');
   [imageForSegmentation,imageForMeasurements,imagesPixelSpacings,imagesPixelCoordinates] = getSlices(obj);
+
+  sliceName = "rotatorCuffMatthieu";
   saveImages(obj,imageForSegmentation,imageForMeasurements,sliceName);
   saveImagesPixelSpacings(obj,imagesPixelSpacings,sliceName);
   saveImagesPixelCoordinates(obj,imagesPixelCoordinates,sliceName);
 end
 
 
 
 function [slice8Bit,slice16Bit,slicePixelSpacings,slicePixelCoordinates] = getSlices(obj)
   SCase = obj.shoulder.SCase;
   scapula = obj.shoulder.scapula;
-  path = SCase.dataDicomPath;
   try
     rotatorCuffSlicer = DicomVolumeSlicer(SCase.getSmoothDicomPath);
   catch
     rotatorCuffSlicer = DicomVolumeSlicer(SCase.dataDicomPath);
   end
 
   % The 3D volume must be normalised for the slices not to be distorted
   rotatorCuffSlicer.setMinimalResolutionDividedBy(2);
   rotatorCuffSlicer.normaliseVolume;
   rotatorCuffSlicer.slice(scapula.coordSys.origin,scapula.coordSys.ML);
 
   % The following lines set the slice to display a 300x300 mm area around the
   % scapula with the "Y" shape of the scapula oriented up
   rotatorCuffSlicer.orientSliceUpwardVector(scapula.coordSys.origin,scapula.coordSys.IS);
   rotatorCuffSlicer.crop(scapula.coordSys.origin-100*scapula.coordSys.IS,300,300);
 
   % Rotator cuff segmentation code requires 1024x1024 pixels images
   if size(rotatorCuffSlicer.sliced,1) > size(rotatorCuffSlicer.sliced,2)
     rotatorCuffSlicer.resize([1024 nan]);
   else
     rotatorCuffSlicer.resize([nan 1024]);
   end
   rotatorCuffSlicer.addEmptyBackgroundToSlice(1024,1024);
 
 
   rawSlice = rotatorCuffSlicer.sliced;
 
   rotatorCuffSlicer.rescaleSliceToUint8;
   slice8Bit = rotatorCuffSlicer.sliced;
   rotatorCuffSlicer.sliced = rawSlice;
 
   rotatorCuffSlicer.rescaleSliceToInt16;
   slice16Bit = rotatorCuffSlicer.sliced;
 
   slicePixelCoordinates = rotatorCuffSlicer.getSlicedPixelCoordinates();
 
   % Flip the slices for right and left shoulders to give similar slices
   if obj.shoulder.side=='L'
     slice8Bit=rot90(flipud(slice8Bit),2);
     slice16Bit=rot90(flipud(slice16Bit),2);
     slicePixelCoordinates=rot90(flipud(slicePixelCoordinates),2);
   end
 
   slicePixelSpacings = rotatorCuffSlicer.slicedPixelSpacings;
 end
 
 
 
 function saveImages(obj,imageForSegmentation,imageForMeasurements,sliceName)
   imwrite(imageForSegmentation,...
-    fullfile(obj.slicesDataPath,sliceName + "_ForSegmentation.png"));
-  save(fullfile(obj.slicesDataPath,sliceName + "_ForMeasurements.mat"),...
+    fullfile(obj.dataSlicesPath,sliceName + "_ForSegmentation.png"));
+  save(fullfile(obj.dataSlicesPath,sliceName + "_ForMeasurements.mat"),...
     'imageForMeasurements');
 end
 
 
 
 function saveImagesPixelSpacings(obj,imagesPixelSpacings,sliceName)
-  save(fullfile(obj.slicesDataPath,sliceName + "_PixelSpacings.mat"),...
+  save(fullfile(obj.dataSlicesPath,sliceName + "_PixelSpacings.mat"),...
     'imagesPixelSpacings');
 end
 
 
 function saveImagesPixelCoordinates(obj,imagesPixelCoordinates,sliceName)
-  save(fullfile(obj.slicesDataPath,sliceName + "_PixelCoordinates.mat"),...
+  save(fullfile(obj.dataSlicesPath,sliceName + "_PixelCoordinates.mat"),...
     'imagesPixelCoordinates');
 end
\ No newline at end of file
diff --git a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m b/ShoulderCase/@RotatorCuff/createSliceNathan.m
similarity index 94%
copy from ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
copy to ShoulderCase/@RotatorCuff/createSliceNathan.m
index 2fd2d24..d8635d9 100644
--- a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
+++ b/ShoulderCase/@RotatorCuff/createSliceNathan.m
@@ -1,368 +1,368 @@
-function createRotatorCuffSlicesNathan(obj,sliceName)
+function createSliceNathan(obj)
   % This function and its results might be reffered as 
   % Nathan algorithm and Nathan slices in some studies.
   %
   %One needs to create two empty files in the current directory :
   %data_degenerescence.mat (a table) and Didnotwork.mat (a
   %string)
   %
   % This script calculates the degeneration of the 4 muscles
   % of the rotator cuff : SS, SC, IS, TM. It is composed of four
   % different steps :
   % - Load CT images (DICOM files) in order to create a 3D matrix
   % consituted with HU (Hounsfield unit) values
   % - Create a slice of this 3D matrix along a patient-specific
   % plane. Do image processing and set saving format usable by
   % automatic segmentation code (Python's code from ARTORG team
   % (UNIBE))
   % - Execute automatic segmentation code on the new created
   % images
   % - Apply degeneration code (from Laboratoire de biom�canique
   % orthop�dique (EPFL)) on new images with automatic
   % segmentation and on image with manual segmentation (from CHUV
   % radiologist). Save degeneration values under the "muscles"
   % properties of each case in the LBO database
   %
   % Some "try-catch" are present in order to catch the errors
   % while keeping running the script on the rest of the database.
   % They are summarized underneath :
   % (DicomReadVolume) : do not manage to load the dicoms
   % Nothing : Automatic segmentation does not work (generally for TM)
   % * : slice function does not work
   % ** : does not want to crop image
   % *** : black image, don't want to save png
 
   % Main code starts here !
 
   %% Load dicom
   SCase = obj.shoulder.SCase;
   dataCTPath = SCase.dataCTPath;
 
   % Create "muscles" folder (and its 4 subfolders) if it does not exist
   folder_muscle=[dataCTPath '/muscles'];
 
   if isfolder(folder_muscle)
   else
   mkdir(dataCTPath,'muscles');
   mkdir(folder_muscle,'IS');
   mkdir(folder_muscle,'SS');
   mkdir(folder_muscle,'SC');
   mkdir(folder_muscle,'TM');
   end
 
   % Use folder "CT-2" (smoothed kernel) if it exists
   dataCTPath2=dataCTPath;
   dataCTPath2(1,end)='2';
 
   if isfolder(dataCTPath2)
       dataCTPath0=dataCTPath2;
   else
       dataCTPath0=dataCTPath;
   end
 
   dicomFolder = [dataCTPath0 '/dicom'];
 
   %Extract the images as a 4D matrix where spatial is the
   %positions of the Pixels in the CT coordinate system
 
   % Return the case's number if it does not manage to load the DICOM files
   try
       [V,spatial,dim] = dicomreadVolume(dicomFolder);
   catch
       outputArg=1;
       return;
   end
 
   %V needs to be 3D matrix containing the value in HU of the
   %corresponding voxel
   V = squeeze(V) ;
 
   %Read the information from the middle dicom image (all image
   %informations are the same for all dicom files) to avoid a problem
   %if other files are located at the beginning or end of the
   %folder
   dicomList = dir(dicomFolder) ;
   size_list = size(dicomList) ;
   dicom_information = dicominfo([dicomFolder '/' dicomList(round(size_list(1)/2)).name]);
 
   %Take into account that the data will need to be linearly transformed
   %from stored space to memory to come back to true HU units afterwards:
   Rescale_slope = dicom_information.RescaleSlope;
   Rescale_intercept = dicom_information.RescaleIntercept;
 
   %Convert to double type to be able to multiply it with a int64
   %matrix afterwards
   V = double(V);
 
   % To filter hard kernel images
 
   %             if isfolder(dataCTPath2)
   %             else
   %             V = imgaussfilt3(V,1) ;
   %             end
 
   size_V = size(V);
 
   %Extract position of each first pixel of each image (upper left corner)
   PatientPositions = spatial.PatientPositions ;
 
   %Get the distance between each pixel in neighboring rows and
   %columns
   PixelSpacings = spatial.PixelSpacings;
 
 
   % Get local coordinate system from parent scapula object(glenoid center)
   origin = obj.shoulder.scapula.coordSys.origin; % do not mix up with origin_image, "origin" is the origin of the Scapula Coordinate System
   xAxis = obj.shoulder.scapula.coordSys.PA;
   yAxis = obj.shoulder.scapula.coordSys.IS;
   zAxis = obj.shoulder.scapula.coordSys.ML;
 
 
   %% Slice from the CT data set
 
   %Get the position of the upper left pixel of the first image
   %and calculate the position of each pixel in the x,y
   %and z direction starting from the "origin" pixel
   origin_image = PatientPositions(1,:); % coordinates of the first pixel
   x_max = origin_image(1)+PixelSpacings(1,1)*(size_V(1)-1);
   y_max = origin_image(2)+PixelSpacings(1,2)*(size_V(2)-1);
   z_max = PatientPositions(end,3);
 
   %Calculate the coefficient of the linear transformation (CT-scan
   %coordinate system (x,y,z) to images coordinates system (i,j,k))
   coefficients_i = polyfit([1 size_V(1)],[origin_image(1) , x_max],1);
   a_i = coefficients_i(1);
   b_i = coefficients_i(2);
 
   coefficients_j = polyfit([1 size_V(2)],[origin_image(2) , y_max],1);
   a_j = coefficients_j(1);
   b_j = coefficients_j(2);
 
   coefficients_k = polyfit([1 size_V(3)],[origin_image(3) , z_max],1);
   a_k = coefficients_k(1);
   b_k = coefficients_k(2);
 
   Pixel_pos_x=(1:size_V(1))*a_i+b_i; % Create the position vector along x-axis
   Pixel_pos_y=(1:size_V(2))*a_j+b_j;
   Pixel_pos_z=(1:size_V(3))*a_k+b_k;
 
   PixSize_z=(z_max-origin_image(3))/(size_V(3)-1); % Space between 2 slices
 
   % The oblique slice to be created is perpendicular to the z
   % axis of the scapula CS and passes through the origin point.
   % We will apply the matlab function "slice", to create the
   % oblique slice. One must first create a slice parallel to 2
   % system's axis. Here we take a plane parallel to the (x,y)
   % plane, passing through the point "origin". Then, we need to
   % rotate this plane so that it becomes normal to the zAxis. V1
   % is the vector normal to the horizontal slice and v2 the
   % vector normal to the final oblique slice, parallel to the
   % zAxis. We need to rotate around a vector v3 and from an angle
   % given by "vrrotvec(v1,v2)".
   v1=[0,0,1];
   v2=zAxis/norm(zAxis);
   rotation_info=vrrotvec(v1,v2);
   v3=rotation_info(1:3);
   theta=rotation_info(4)/pi*180; % Convert it in degrees
 
   % Create the plane to be rotated
   [x,y,z] = meshgrid(Pixel_pos_x,Pixel_pos_y,Pixel_pos_z);
   s=1000; % Extension factor of the slice in order to cut the whole matrix.Set arbitrarily but sufficiently big so that it works for every patient
   Pixel_pos_x_enlarged=[origin_image(1)+(-PixelSpacings(1,1)*s:PixelSpacings(1,1):-PixelSpacings(1,1)) Pixel_pos_x x_max+(PixelSpacings(1,1):PixelSpacings(1,1):PixelSpacings(1,1)*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged=[origin_image(2)+(-PixelSpacings(1,2)*s:PixelSpacings(1,2):-PixelSpacings(1,2)) Pixel_pos_y y_max+(PixelSpacings(1,2):PixelSpacings(1,2):PixelSpacings(1,2)*s)];
 
   figure('visible','off');
   hsp = surf(Pixel_pos_x_enlarged,Pixel_pos_y_enlarged,origin(3)*ones(length(Pixel_pos_y_enlarged),length((Pixel_pos_x_enlarged))));
 
   rotate(hsp,v3,theta,origin);
 
   xd_0 = hsp.XData; % x-coordinates of the plane's points
   yd_0 = hsp.YData; % y-coordinates of the plane's points
   zd_0 = hsp.ZData; % z-coordinates of the plane's points
 
   % Slice the 3D matrix of CT-scan
   try % Return case's number in case slice does not work
   figure('visible','off');
   hss=slice(x,y,z,V,xd_0,yd_0,zd_0,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
   catch
   outputArg=2;
   return
   end
 
   HU_0 = get(hss,'CData'); % HU values at each point of the plane
 
   NaNplaces_0=find(isnan(HU_0));
 
   HU_1=HU_0*Rescale_slope+Rescale_intercept; % getting true HU values
   p_0=10e6; % value outside 3D matrix, set to a high distinguishable number
   HU_1(NaNplaces_0)=p_0;
 
   % If it's a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_1=flipud(HU_1');
   else
   HU_1=HU_1';
   end
 
   % Find the new pixel's sizes necessary to obtain a 1024*1024
   % shoulder slice. The largest diagonal's length is set to 1024
   % pixels
   [left_1,left_2]=find(HU_1~=p_0,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_1,2)~=p_0,1);
   right_r=right_1;
   right_c=size(HU_1,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_1)~=p_0,1);
   up_c=size(HU_1)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_1,3)~=p_0,1);
   down_r=size(HU_1,1)-down_2+1;
   down_c=down_1;
 
   % In case the slice does not cut the whole matrix (which
   % normally never happens), ask to double "s" parameter
 
   if or(or(right_c==size(HU_1,2),left_c==1),or(up_r==1,down_r==size(HU_1,1)))
     outputArg = 3;
     return
   end
 
   % Find which diagonal of the matrix' effective slice is the largest
   w=right_c-left_c; % width
   h=down_r-up_r; % height
   if h>w % Set the largest diagonal at 1024 pixels
       l=h;
   else
       l=w;
   end
 
   % Define new pixel's size
   NewPixSize_x=PixelSpacings(1,1)*l/1023;
   NewPixSize_y=PixelSpacings(1,2)*l/1023;
 
   % New plane with the new pixel's sizes
   Pixel_pos_x_enlarged_2=[origin_image(1)+(-NewPixSize_x*s:NewPixSize_x:-NewPixSize_x) origin_image(1):NewPixSize_x:x_max x_max+(NewPixSize_x:NewPixSize_x:NewPixSize_x*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged_2=[origin_image(2)+(-NewPixSize_y*s:NewPixSize_y:-NewPixSize_y) origin_image(2):NewPixSize_y:y_max y_max+(NewPixSize_y:NewPixSize_y:NewPixSize_y*s)];
 
   figure('visible','off');
   hsp2 = surf(Pixel_pos_x_enlarged_2,Pixel_pos_y_enlarged_2,origin(3)*ones(length(Pixel_pos_y_enlarged_2),length((Pixel_pos_x_enlarged_2))));
   rotate(hsp2,v3,theta,origin);
 
   xd = hsp2.XData; % x-coordinates of the rotated plane's points.
   yd = hsp2.YData; % y-coordinates of the rotated plane's points.
   zd = hsp2.ZData; % z-coordinates of the rotated plane's points.
 
   % New slice with the good resolution (1024*1024)
   figure('visible','off');
   hss=slice(x,y,z,V,xd,yd,zd,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
 
   HU_2 = get(hss,'CData'); % HU values at each point of the plane
   NaNplaces=find(isnan(HU_2));
   HU_slice=HU_2*Rescale_slope+Rescale_intercept;
   HU_slice_8_bit=HU_slice;
   p_tif=0; % value outside 3D matrix, for the tif image
   p_png=85; % for the png image
   HU_slice(NaNplaces)=p_tif;
 
   % Scaling to fit the training set (for the png only, grayscale)
   HU_interval_png=[-100 160];
   lin_coef=[HU_interval_png(1) 1 ; HU_interval_png(2) 1]^(-1)*[0 255]';
 
   HU_slice_8_bit=lin_coef(1)*HU_slice_8_bit+lin_coef(2);
   HU_slice_8_bit(NaNplaces)=p_png;
 
   % If it a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_slice=flipud(HU_slice');
   HU_slice_8_bit=flipud(HU_slice_8_bit');
   else
   HU_slice=HU_slice';
   HU_slice_8_bit=HU_slice_8_bit';
   end
 
   % Determine the location of the 4 corners of the slice
   [left_1,left_2]=find(HU_slice~=p_tif,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_slice,2)~=p_tif,1);
   right_r=right_1;
   right_c=size(HU_slice,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_slice)~=p_tif,1);
   up_c=size(HU_slice,2)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_slice,3)~=p_tif,1);
   down_r=size(HU_slice,1)-down_2+1;
   down_c=down_1;
 
   Folder=[dataCTPath '/matlab'];
 
   % Crop the image
   try % If cropping does not work, return case's number
   HU_slice=HU_slice(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice_8_bit=HU_slice_8_bit(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice=int16(HU_slice); % convert into 16-bit
   catch
   outputArg=4;
   return
   end
 
   % Create '.tif' file
   file_name_tif=[SCase.id '.tif'];
 
   t = Tiff(fullfile(Folder,file_name_tif),'w');
   setTag(t,'Photometric',Tiff.Photometric.MinIsBlack);
   setTag(t,'Compression',Tiff.Compression.None);
   setTag(t,'BitsPerSample',16);
   setTag(t,'SamplesPerPixel',1);
   setTag(t,'ImageLength',1024);
   setTag(t,'ImageWidth',1024);
   setTag(t,'PlanarConfiguration',Tiff.PlanarConfiguration.Chunky);
   setTag(t,'SampleFormat',Tiff.SampleFormat.Int);
   write(t,HU_slice);
   close(t);
 
   % Create 8 bit '.png' file and save it at the right places
   % (matlab folder and python's input folder)
   file_name_png=[SCase.id '.png'];
 
   HU_slice_8_bit=uint8(HU_slice_8_bit);
   figure('visible','off');
   image(HU_slice_8_bit);
   colormap(gray(256));
 
   try % return case's number if does not manage to write png in the "matlab" folder because it is a black image
   imwrite(HU_slice_8_bit,[Folder '/' file_name_png]);
   catch
   outputArg=5;
   return
   end
 
   pythonDir = ConfigFileExtractor.getVariable('pythonDir');
   Folder_input_py= [pythonDir '/input'];
   imwrite(HU_slice_8_bit,[Folder_input_py '/' file_name_png]); % write png file in the python input folder
   % imwrite(HU_slice_8_bit,['./rcseg/' file_name_png]); % write png file in local /rceg folder
 
   % Create .txt file with pixel dimensions
   file_name_info=[SCase.id '.info'];
   text1='# Amira Stacked Slices' ;
   text2=['pixelsize ', num2str(NewPixSize_x),' ',num2str(NewPixSize_y)];
   fileID = fopen([Folder '/' file_name_info],'w');
   fprintf(fileID,'%s\r\n\r\n',text1);
   fprintf(fileID,'%s\n',text2);
   fclose(fileID);
   close all;
 
-
-  imwrite(HU_slice_8_bit,fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']));
+  sliceName = "rotatorCuffNathan";
+  imwrite(HU_slice_8_bit,fullfile(obj.dataSlicesPath,sliceName + "_ForSegmentation.png"));
   imageForMeasurements = HU_slice;
-  save(fullfile(obj.slicesDataPath,[sliceName '_ForMeasurements.mat']),'imageForMeasurements');
+  save(fullfile(obj.dataSlicesPath,sliceName + "_ForMeasurements.mat"),'imageForMeasurements');
   imagesPixelSpacings = [NewPixSize_x NewPixSize_y];
-  save(fullfile(obj.slicesDataPath,[sliceName '_PixelSpacings.mat']),'imagesPixelSpacings');
+  save(fullfile(obj.dataSlicesPath,sliceName + "_PixelSpacings.mat"),'imagesPixelSpacings');
 end
diff --git a/ShoulderCase/@RotatorCuff/plot.m b/ShoulderCase/@RotatorCuff/plot.m
new file mode 100644
index 0000000..e6644e1
--- /dev/null
+++ b/ShoulderCase/@RotatorCuff/plot.m
@@ -0,0 +1,38 @@
+function output = plot(obj,varargin)
+
+  if nargin == 2
+    error('Provide the slice name and the mask name or provide nothing (auto slice and segmentations will be plotted).')
+  elseif nargin == 3
+    sliceName = varargin{1} + "_ForSegmentation.png";
+    maskName = varargin{2} + "_Mask.png";
+  else
+    sliceName = "rotatorCuffMatthieu_ForSegmentation.png";
+    maskName = "autoMatthieu_Mask.png";
+  end
+
+  rotatorCuffCrossSection = imread(fullfile(obj.dataSlicesPath,sliceName));
+
+  leg = {};
+  plotHandle(1) = imshow(rotatorCuffCrossSection);
+  hold on
+  muscleColor = containers.Map;
+  muscleColor("SC") = "b";
+  muscleColor("SS") = "r";
+  muscleColor("IS") = "g";
+  muscleColor("TM") = "y";
+
+  for muscleName = ["SC", "SS", "IS", "TM"]
+    try
+      segmentedImage = imread(fullfile(obj.(muscleName).dataMaskPath,maskName));
+    catch
+      continue
+    end
+
+    if (max(segmentedImage,[],'all') > 0)   % if segmentation result is not empty
+      leg{end+1} = muscleName;
+      [~,plotHandle(end+1)] = contour(segmentedImage,'color',muscleColor(muscleName));
+    end
+  end
+  plotHandle(end+1) = legend(leg);
+  output = plotHandle;
+end
diff --git a/ShoulderCase/@RotatorCuff/plot3.m b/ShoulderCase/@RotatorCuff/plot3.m
new file mode 100644
index 0000000..fc824ae
--- /dev/null
+++ b/ShoulderCase/@RotatorCuff/plot3.m
@@ -0,0 +1,38 @@
+function output = plot3(obj)
+  muscleColor = containers.Map;
+  muscleColor("SC") = "b";
+  muscleColor("SS") = "r";
+  muscleColor("IS") = "g";
+  muscleColor("TM") = "y";
+  plotHandle = [];
+
+  for muscleName = ["SC", "SS", "IS", "TM"]
+    try
+      % Plot segmentations
+      centroidPoint = obj.(muscleName).getCentroid().coordinates;
+      contourPoints = obj.(muscleName).getContour().coordinates;
+      plotHandle(end+1) = scatter3(centroidPoint(:,1), centroidPoint(:,2), centroidPoint(:,3),...
+        "MarkerEdgeColor", muscleColor(muscleName),...
+        "MarkerFaceColor", muscleColor(muscleName));
+      hold on
+      plotHandle(end+1) = scatter3(contourPoints(:,1), contourPoints(:,2), contourPoints(:,3),...
+        "MarkerEdgeColor", muscleColor(muscleName));
+
+        
+      % Plot muscle insertion
+      humeralInsertion = obj.(muscleName).insertions;
+      plotHandle(end+1) = plot(...
+        Vector(obj.(muscleName).forceVector.origin, humeralInsertion),...
+        "Color", muscleColor(muscleName),...
+        "LineWidth", 2);
+      plotHandle(end+1) = scatter3(...
+        humeralInsertion(1),...
+        humeralInsertion(2),...
+        humeralInsertion(3),...
+        "MarkerFaceColor", muscleColor(muscleName),...
+        "MarkerEdgeColor", "magenta");
+    end
+  end
+  
+  output = plotHandle;
+end
\ No newline at end of file
diff --git a/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m b/ShoulderCase/@RotatorCuff/segmentMuscles.m
similarity index 69%
copy from ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m
copy to ShoulderCase/@RotatorCuff/segmentMuscles.m
index 587d386..302f293 100644
--- a/ShoulderCase/@MusclesContainer/segmentRotatorCuffMuscles.m
+++ b/ShoulderCase/@RotatorCuff/segmentMuscles.m
@@ -1,72 +1,64 @@
-function segmentRotatorCuffMuscles(obj,sliceName,maskName)
+function segmentMuscles(obj, sliceName, maskName)
   % sliceName is the part before the '_ForSegmentation.png' part of the name of
   % the file that is sent to rcseg.
   % maskName is the part before the '_Mask.png' part of the name of
   % the file that is saved in the SCase folder.
   cleanRotatorCuffSegmentationWorkspace(obj);
   sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName);
   callRotatorCuffSegmentation(obj);
-
-  % (Re)initialize the muscle objects to refer to the new segmentation results
-  obj.addMuscle('IS',sliceName);
-  obj.addMuscle('SS',sliceName);
-  obj.addMuscle('SC',sliceName);
-  obj.addMuscle('TM',sliceName);
-
   saveSegmentationResults(obj,maskName);
   cleanRotatorCuffSegmentationWorkspace(obj);
 end
 
 
 
 function cleanRotatorCuffSegmentationWorkspace(obj)
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
   delete(fullfile(rotatorCuffSegmentationPath,'input','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'IS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SC','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'TM','*'));
 end
 
 
 
 function sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName)
   SCase = obj.shoulder.SCase;
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
-  imageForSegmentationPath = fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']);
+  imageForSegmentationPath = fullfile(obj.dataSlicesPath,sliceName + "_ForSegmentation.png");
   copyfile(imageForSegmentationPath,fullfile(rotatorCuffSegmentationPath,'input',[SCase.id '.png']));
 end
 
 
 
 function callRotatorCuffSegmentation(obj)
   %% Apply UNIBE method to automatically segment normal muscle
 
   pythonDir = ConfigFileExtractor.getVariable('pythonDir');
 
   pythonCommandActivateEnvironment= ['source ' fullfile(pythonDir,'venv','bin','activate') ';'];
   pythonCommandMoveToSegmentationWorkspace = ['cd ' pythonDir ';'];
   pythonCommandExecuteSegmentation = [fullfile(pythonDir,'rcseg.py') ' segment input' ';'];
 
   [~, pythonCommandOutput] = system([pythonCommandActivateEnvironment,...
-                                     pythonCommandMoveToSegmentationWorkspace,...
-                                     pythonCommandExecuteSegmentation]);
+    pythonCommandMoveToSegmentationWorkspace,...
+    pythonCommandExecuteSegmentation]);
   disp(pythonCommandOutput);
 end
 
 
 
 function saveSegmentationResults(obj,maskName)
   SCaseID = obj.shoulder.SCase.id;
-  muscleNames = fields(obj.list);
   rotatorCuffSegmentationPath = ConfigFileExtractor.getVariable('pythonDir');
-
-  for i = 1:length(muscleNames)
+  
+  for muscleName = ["SC", "SS", "IS", "TM"]
     try
-      copyfile(fullfile(rotatorCuffSegmentationPath,muscleNames{i},[SCaseID '.png']),...
-               fullfile(obj.list.(muscleNames{i}).maskDataPath,[maskName '_Mask.png']) );
+      copyfile(fullfile(rotatorCuffSegmentationPath,muscleName,SCaseID + ".png"),...
+        fullfile(obj.(muscleName).dataMaskPath,maskName + "_Mask.png") );
     catch ME
       warning(ME.message)
     end
   end
 end
diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m
index b30b44b..e4363ff 100644
--- a/ShoulderCase/@Shoulder/Shoulder.m
+++ b/ShoulderCase/@Shoulder/Shoulder.m
@@ -1,24 +1,25 @@
 classdef (Abstract) Shoulder < handle
 % Contains all the shoulder parts (bones) and their measurements.
 %
 % Can be used to plot an overview of the case.
 
   properties
     side = ""; % R or L
     scapula
     humerus
-    muscles
+    rotatorCuff
     CTscan = "";
     comment = "";
     SCase
   end
 
   methods
     function obj = Shoulder(SCase)
       obj.SCase   = SCase;
       obj.humerus = Humerus(obj);
       obj.muscles = MusclesContainer(obj);
+      obj.rotatorCuff = RotatorCuff(obj);
     end
   end
 
 end
diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m
index b31d9a0..deb1397 100755
--- a/ShoulderCase/@ShoulderCase/ShoulderCase.m
+++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m
@@ -1,837 +1,839 @@
 classdef ShoulderCase < handle
     % Properties and methods associated to the SoulderCase object.
 
     % Author: Alexandre Terrier, EPFL-LBO
     %         Matthieu Boubat, EPFL-LBO
     % Creation date: 2018-07-01
     % Revision date: 2020-07-23
 
     % TODO:
     % Need method to load properties (from amira, matlab, excel, SQL)
     % Need method to save properties (in matlab, excel, SQL)
     % Remove datapath from constants
     %
     properties
         id = [];    % id of the shoulder case, as it appears in the database (Pnnn)
         id4C = []; % SCaseId with P/N followed by 3 digits --> 4 char
         diagnosis = [];
         treatment = [];
         outcome = [];
         patient = [];
         shoulderManual = [];
         shoulderAuto = [];
         study = [];
         comment = [];
 
         dataCTPath = [];
     end
 
 
     methods (Access = ?ShoulderCaseLoader)
 
         function obj = ShoulderCase(SCaseID,dataCTPath)
             % SCaseID validation
             rawSCaseID = SCaseIDParser(SCaseID);
             assert(rawSCaseID.isValidID,'The input argument is not a valid SCaseID.')
             obj.id = SCaseID;
             obj.id4C = rawSCaseID.getIDWithNumberOfDigits(3);
 
             % path attribution
             obj.dataCTPath = dataCTPath;
 
             % Folder creation
             if not(isfolder(obj.dataMatlabPath))
                 mkdir(obj.dataMatlabPath);
             end
 
             % Initialsation
             obj.patient   = Patient(obj);
             obj.shoulderManual  = ShoulderManual(obj);
             obj.shoulderAuto  = ShoulderAuto(obj);
         end
 
       end
 
 
 
       methods
 
         function output = dataPath(obj)
             output = fileparts(obj.dataCTPath);
         end
 
         function output = dataAmiraPath(obj)
             output = fullfile(obj.dataCTPath, "amira");
         end
 
         function output = dataMatlabPath(obj)
             output = fullfile(obj.dataCTPath, "matlab");
         end
 
         function output = dataDicomPath(obj)
             output = fullfile(obj.dataCTPath, "dicom");
         end
 
         function outputArg = output(obj, varargin)
             % Below is copy/past from previous version
             % This function is used to output the variable Results, used
             % by scapula_measurePKG the modified funcion scapula_measure.
             % inputArg could be as in scapula_calculation
             % 'density', 'References', 'obliqueSlice', 'display'
 
             Result.SCase_id               = obj.id;                                      %  1
             Result.glenoid_Radius         = obj.shoulderManual.scapula.glenoid.radius;         %  2
             Result.glenoid_radiusRMSE     = obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE;     %  3                              %  3
             %Result.glenoidSphericity     = ' ';
             %Result.glenoid_biconcave     = ' ';
             Result.glenoid_depth          = obj.shoulderManual.scapula.glenoid.depth;          %  4
             Result.glenoid_width          = obj.shoulderManual.scapula.glenoid.width;          %  5
             Result.glenoid_height         = obj.shoulderManual.scapula.glenoid.height;         %  6
             Result.glenoid_center_PA      = obj.shoulderManual.scapula.glenoid.centerLocal(1); %  7
             Result.glenoid_center_IS      = obj.shoulderManual.scapula.glenoid.centerLocal(2); %  8
             Result.glenoid_center_ML      = obj.shoulderManual.scapula.glenoid.centerLocal(3); %  9
             Result.glenoid_version_ampl   = obj.shoulderManual.scapula.glenoid.versionAmplitude;    % 10
             Result.glenoid_version_orient = obj.shoulderManual.scapula.glenoid.versionOrientation;  % 11
             Result.glenoid_Version        = obj.shoulderManual.scapula.glenoid.version;        % 12
             Result.glenoid_Inclination    = obj.shoulderManual.scapula.glenoid.inclination;    % 13
             Result.humerus_joint_radius   = ' ';                                           % 14
             Result.humeral_head_radius    = obj.shoulderManual.humerus.radius;                 % 15
             % Result.humerus_GHsublux_2D = ' ';
             % Result.humerus_SHsublux_2D = ' ';
             Result.humerus_GHsubluxation_ampl   = obj.shoulderManual.humerus.GHSAmpl;          % 16
             Result.humerus_GHsubluxation_orient = obj.shoulderManual.humerus.GHSOrient;        % 17
             Result.humerus_SHsubluxation_ampl   = obj.shoulderManual.humerus.SHSAmpl;          % 18
             Result.humerus_SHsubluxation_orient = obj.shoulderManual.humerus.SHSOrient;        % 19
             Result.scapula_CSA            = obj.shoulderManual.scapula.acromion.criticalShoulderAngle; % radCSA;
             % 5 Lines below should be updated
             Result.scapula_CTangle = 0; %CTorientation; %20
             Result.scapula_CTangleVersion = 0; % WearPlaneAngle; %21
             Result.scapula_CTangleSHS = 0; % SHSPlaneAngle; % 22
             Result.scapula_CTangleGHS = 0; % GHSPlaneAngle; % 23
             Result.scapula_PlaneRMSE = 0; %PlaneRMSE;%24
 
             Result.scapula_AI             = obj.shoulderManual.scapula.acromion.acromionIndex;
 
             outputArg = Result;
         end
 
         function outputArg = saveMatlab(obj)
             % Save SCase to matlab file
             dir = obj.dataMatlabPath;
             % Create dir if not exist
             try
                 if ~exist(dir, 'dir') % create directory if it does not exist
                     mkdir(dir);
                 end
             catch
                 fprintf('Error creating the matlab directory \n'); % Should be in log
             end
 
             % Save SCase in matlab directoty, in a file named SCaseCNNN.m
             filename = 'SCase';
-            filename = [dir '/' filename '.mat'];
+            filename = fullfile(obj.dataMatlabPath, filename + ".mat");
+            dataCTPath = obj.dataCTPath;
             try
                 SCase = obj;
                 % Delete CT path wich must be set when loading/creating the 
                 % SCase to avoid messing with paths on different systems.
                 SCase.dataCTPath = [];
                 save(filename, 'SCase');
                 outputArg = 1;
             catch
                 fprintf('Error creating SCase matlab file \n'); % Should be in log
                 outputArg = -1;
             end
+            obj.dataCTPath = dataCTPath;
         end
 
         function outputArg = appendToCSV(obj,filename)
             % Save SCase to csv file
             logFid  = fopen('log/measureSCase.log', 'a');
 
             dataDir = ConfigFileExtractor.getVariable('dataDir');
             xlsDir  = [dataDir '/Excel/xlsFromMatlab'];
 
             fid = fopen([xlsDir '/' filename],'a');
             fprintf(fid,[...
                 obj.id ','...                                                          % SCase_id
                 obj.shoulderManual.side ','...                                           % shoulder_side
                 num2str(obj.shoulderManual.scapula.glenoid.radius) ','...                % glenoid_radius
                 num2str(obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE) ','...     % glenoid_sphereRMSE
                 num2str(obj.shoulderManual.scapula.glenoid.depth) ','...                 % glenoid_depth
                 num2str(obj.shoulderManual.scapula.glenoid.width) ','...                 % glenoid_width
                 num2str(obj.shoulderManual.scapula.glenoid.height) ','...                % glenoid_height
                 num2str(obj.shoulderManual.scapula.glenoid.centerLocal.x) ','...         % glenoid_centerPA
                 num2str(obj.shoulderManual.scapula.glenoid.centerLocal.y) ','...         % glenoid_centerIS
                 num2str(obj.shoulderManual.scapula.glenoid.centerLocal.z) ','...         % glenoid_centerML
                 num2str(obj.shoulderManual.scapula.glenoid.versionAmplitude) ','...           % glenoid_versionAmpl
                 num2str(obj.shoulderManual.scapula.glenoid.versionOrientation) ','...         % glenoid_versionOrientation
                 num2str(obj.shoulderManual.scapula.glenoid.version) ','...               % glenoid_version
                 num2str(obj.shoulderManual.scapula.glenoid.inclination) ','...           % glenoid_inclination
                 num2str(obj.shoulderManual.humerus.jointRadius) ','...                   % humerus_jointRadius
                 num2str(obj.shoulderManual.humerus.radius) ','...                        % humerus_headRadius
                 num2str(obj.shoulderManual.humerus.GHSAmpl) ','...                       % humerus_GHSAmpl
                 num2str(obj.shoulderManual.humerus.GHSOrient) ','...                     % humerus_GHSOrient
                 num2str(obj.shoulderManual.humerus.SHSAmpl) ','...                       % humerus_SHSAmpl
                 num2str(obj.shoulderManual.humerus.SHSOrient) ','...                     % humerus_SHSOrient
                 num2str(obj.shoulderManual.humerus.SHSAngle) ','...                      % humerus_SHSAgle
                 num2str(obj.shoulderManual.humerus.SHSPA) ','...                         % humerus_SHSPA
                 num2str(obj.shoulderManual.humerus.SHSIS) ','...                         % humerus_SHSIS
                 num2str(obj.shoulderManual.scapula.acromion.AI) ','...                   % acromion_AI
                 num2str(obj.shoulderManual.scapula.acromion.CSA) ','...                  % acromion_CSA
                 num2str(obj.shoulderManual.scapula.acromion.PSA) ','...                  % acromion_PSA
                 num2str(obj.shoulderManual.scapula.acromion.AAA) '\n'...                 % acromion_AAA
                 ]);
             fclose(fid);
             fclose(logFid);
 
             outputArg = 1;
 
         end
 
         function outputArg = saveExcel(~)
             % Save SCase to Excel file
 
 
             outputArg = 1;
 
         end
 
         function outputArg = saveSQL(~)
             % Save SCase to MySQL database
 
 
             outputArg = 1;
 
         end
         
         function outputArg = SCaseDataToTable(obj, varargin)
             % This method export all of the measurements and data of a
             % SCase to a table, wgich then can be used to export all of the
             % measurements as csv or excel file (see exportSCaseData function)
             
             %char(n), nan(a, b) and "" are where the calculations have not
             %been done. I used these based on the field of data in other
             %SCases. For examnple where we want char, I used char(nan).
           
             Result.SCase_id                             = obj.id;     
             if isnan(obj.diagnosis)
                 Result.SCase_diagnosis                      = char(nan);
             else
                 Result.SCase_diagnosis                      = obj.diagnosis;
             end
             if isempty(obj.treatment)
                 Result.SCase_treatment                      = char(nan);
             elseif isnan(obj.treatment) 
                 Result.SCase_treatment                      = char(nan);
             else
                 Result.SCase_treatment                      = obj.treatment;
             end    
             if isempty(obj.outcome)
                 Result.SCase_outcome                        = char(nan);
             else
                 Result.SCase_outcome                        = obj.outcome;
             end
             %Result.patient_id                           = obj.patient.id;
             %Result.patient_idMed                        = obj.patient.idMed;
             Result.patient_gender                       = obj.patient.gender;
             Result.patient_age                          = obj.patient.age;
             Result.patient_ethnicity                    = obj.patient.ethnicity;
             Result.patient_weight                       = obj.patient.weight;
             Result.patient_height                       = obj.patient.height;
             
             if isempty(obj.patient.BMI)
                 Result.patient_BMI                      = nan;
             elseif obj.patient.BMI(1) == '#'
                 Result.patient_BMI                      = nan;
             else
                 Result.patient_BMI                      = obj.patient.BMI;
             end    
             
             Result.patient_comment                      = obj.patient.comment;
             
             %shoulderManual
             Result.shoulderManual_side                  = obj.shoulderManual.side;
             if isempty(obj.shoulderManual.scapula.angulusInferior)
                 Result.shoulderManual_scapula_angulusInferior   = nan(1,3);
                 Result.shoulderManual_scapula_angulusInferiorLocal = nan(1,3);
                 Result.shoulderManual_scapula_trigonumSpinae    = nan(1,3);
                 Result.shoulderManual_scapula_processusCoracoideus = nan(1,3);
                 Result.shoulderManual_scapula_acromioClavicular = nan(1,3);
                 Result.shoulderManual_scapula_angulusAcromialis = nan(1,3);
                 Result.shoulderManual_scapula_spinoGlenoidNotch = nan(1,3);
             
             else
                 Result.shoulderManual_scapula_angulusInferior   = obj.shoulderManual().scapula.angulusInferior;
                 Result.shoulderManual_scapula_angulusInferiorLocal = obj.shoulderManual().scapula.angulusInferiorLocal;
                 Result.shoulderManual_scapula_trigonumSpinae    = obj.shoulderManual().scapula.trigonumSpinae;
                 Result.shoulderManual_scapula_processusCoracoideus = obj.shoulderManual().scapula.processusCoracoideus;
                 Result.shoulderManual_scapula_acromioClavicular = obj.shoulderManual().scapula.acromioClavicular;
                 Result.shoulderManual_scapula_angulusAcromialis = obj.shoulderManual().scapula.angulusAcromialis;
                 Result.shoulderManual_scapula_spinoGlenoidNotch = obj.shoulderManual().scapula.spinoGlenoidNotch;
  
             end
             if isempty(obj.shoulderManual.scapula.coordSys.ML)
                 Result.shoulderManual_scapula_coordSys_ML = nan(1,3);
                 Result.shoulderManual_scapula_coordSys_PA = nan(1,3);
                 Result.shoulderManual_scapula_coordSys_IS = nan(1,3);
                 Result.shoulderManual_scapula_coordSys_origin = nan(1,3);
             else
                 Result.shoulderManual_scapula_coordSys_ML = obj.shoulderManual.scapula.coordSys.ML;
                 Result.shoulderManual_scapula_coordSys_PA = obj.shoulderManual.scapula.coordSys.PA;
                 Result.shoulderManual_scapula_coordSys_IS = obj.shoulderManual.scapula.coordSys.IS;
                 Result.shoulderManual_scapula_coordSys_origin = obj.shoulderManual.scapula.coordSys.origin;
             end 
             if isempty(obj.shoulderManual.scapula.glenoid.center)
                 Result.shoulderManual_scapula_glenoid_center = nan(1,3);
                 Result.shoulderManual_scapula_glenoid_centerLocal_x = nan;
                 Result.shoulderManual_scapula_glenoid_centerLocal_y = nan;
                 Result.shoulderManual_scapula_glenoid_centerLocal_z = nan;
                 Result.shoulderManual_scapula_glenoid_radius = nan;  
                 Result.shoulderManual_scapula_glenoid_centerLine = nan(1,3);
             else
                 Result.shoulderManual_scapula_glenoid_center = obj.shoulderManual.scapula.glenoid.center;
                 Result.shoulderManual_scapula_glenoid_centerLocal_x = obj.shoulderManual.scapula.glenoid.centerLocal.x;
                 Result.shoulderManual_scapula_glenoid_centerLocal_y = obj.shoulderManual.scapula.glenoid.centerLocal.y;
                 Result.shoulderManual_scapula_glenoid_centerLocal_z = obj.shoulderManual.scapula.glenoid.centerLocal.z;
                 Result.shoulderManual_scapula_glenoid_radius = obj.shoulderManual.scapula.glenoid.radius;   
                 Result.shoulderManual_scapula_glenoid_centerLine = obj.shoulderManual.scapula.glenoid.centerLine;
             end    
             if isempty(obj.shoulderManual.scapula.glenoid.depth)
                 Result.shoulderManual_scapula_glenoid_depth = nan;
             else 
                 Result.shoulderManual_scapula_glenoid_depth = obj.shoulderManual.scapula.glenoid.depth; 
             end   
             if isempty(obj.shoulderManual.scapula.glenoid.width)
                 Result.shoulderManual_scapula_glenoid_width = nan;
             else 
                 Result.shoulderManual_scapula_glenoid_width = obj.shoulderManual.scapula.glenoid.width; 
             end 
             if isempty(obj.shoulderManual.scapula.glenoid.height)
                 Result.shoulderManual_scapula_glenoid_height = nan;
             else 
                 Result.shoulderManual_scapula_glenoid_height = obj.shoulderManual.scapula.glenoid.height; 
             end 
             if isempty(obj.shoulderManual.scapula.glenoid.anteroSuperiorAngle)
                 Result.shoulderManual_scapula_glenoid_anteroSuperiorAngle = nan;
             else 
                 Result.shoulderManual_scapula_glenoid_anteroSuperiorAngle = obj.shoulderManual.scapula.glenoid.anteroSuperiorAngle; 
             end 
             if isempty(obj.shoulderManual.scapula.glenoid.versionAmplitude)
                 Result.shoulderManual_scapula_glenoid_versionAmplitude = nan;
                 Result.shoulderManual_scapula_glenoid_versionOrientation = nan;
                 Result.shoulderManual_scapula_glenoid_version = nan;
                 Result.shoulderManual_scapula_glenoid_inclination = nan;
             else 
                 Result.shoulderManual_scapula_glenoid_versionAmplitude = obj.shoulderManual.scapula.glenoid.versionAmplitude; 
                 Result.shoulderManual_scapula_glenoid_versionOrientation = obj.shoulderManual.scapula.glenoid.versionOrientation;
                 Result.shoulderManual_scapula_glenoid_version = obj.shoulderManual.scapula.glenoid.version;
                 Result.shoulderManual_scapula_glenoid_inclination = obj.shoulderManual.scapula.glenoid.inclination;
             end
             if isempty(obj.shoulderManual.scapula.glenoid.density)
                 Result.shoulderManual_scapula_glenoid_density = nan(1, 6);
             else
                 Result.shoulderManual_scapula_glenoid_density = obj.shoulderManual.scapula.glenoid.density;
             end    
             if isempty(obj.shoulderManual.scapula.glenoid.comment)
                 Result.shoulderManual_scapula_glenoid_comment     = nan;
             else
                 Result.shoulderManual_scapula_glenoid_comment     = obj.shoulderManual.scapula.glenoid.comment;
             end
             if isempty(obj.shoulderManual.scapula.glenoid.fittedSphere.center)
                 Result.shoulderManual_scapula_glenoid_fittedSphere_center = nan(1,3);
                 Result.shoulderManual_scapula_glenoid_fittedSphere_radius = nan;
                 Result.shoulderManual_scapula_glenoid_fittedSphere_R2 = nan;
                 Result.shoulderManual_scapula_glenoid_fittedSphere_RMSE = nan;
             else
                 Result.shoulderManual_scapula_glenoid_fittedSphere_center = ...
                     obj.shoulderManual.scapula.glenoid.fittedSphere.center;
                 Result.shoulderManual_scapula_glenoid_fittedSphere_radius = ...
                     obj.shoulderManual.scapula.glenoid.fittedSphere.radius;
                 Result.shoulderManual_scapula_glenoid_fittedSphere_R2 = ...
                     obj.shoulderManual.scapula.glenoid.fittedSphere.R2;
                 Result.shoulderManual_scapula_glenoid_fittedSphere_RMSE = ...
                     obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE;
             end    
             if isempty(obj.shoulderManual.scapula.glenoid.walch)
                 Result.shoulderManual_scapula_glenoid_walch  = char(nan);
             elseif isnan(obj.shoulderManual.scapula.glenoid.walch) 
                 Result.shoulderManual_scapula_glenoid_walch  = char(nan);
             else
                 Result.shoulderManual_scapula_glenoid_walch  = obj.shoulderManual.scapula.glenoid.walch;
             end
             if isempty(obj.shoulderManual.scapula.acromion.AI)
                 Result.shoulderManual_scapula_acromion_AI = nan;
                 Result.shoulderManual_scapula_acromion_CSA = nan;
                 Result.shoulderManual_scapula_acromion_PSA = nan;
                 Result.shoulderManual_scapula_acromion_PSL = nan;
                 Result.shoulderManual_scapula_acromion_AAA = nan;
                 Result.shoulderManual_scapula_acromion_AAL = nan;
                 Result.shoulderManual_scapula_acromion_AAx = nan;
                 Result.shoulderManual_scapula_acromion_AAy = nan;
                 Result.shoulderManual_scapula_acromion_AAz = nan;
                 Result.shoulderManual_scapula_acromion_ACx = nan;
                 Result.shoulderManual_scapula_acromion_ACy = nan;
                 Result.shoulderManual_scapula_acromion_ACz = nan;
             else    
                 Result.shoulderManual_scapula_acromion_AI = obj.shoulderManual.scapula.acromion.AI;
                 Result.shoulderManual_scapula_acromion_CSA = obj.shoulderManual.scapula.acromion.CSA;
                 Result.shoulderManual_scapula_acromion_PSA = obj.shoulderManual.scapula.acromion.PSA;
                 Result.shoulderManual_scapula_acromion_PSL = obj.shoulderManual.scapula.acromion.PSL;
                 Result.shoulderManual_scapula_acromion_AAA = obj.shoulderManual.scapula.acromion.AAA;
                 Result.shoulderManual_scapula_acromion_AAL = obj.shoulderManual.scapula.acromion.AAL;
                 Result.shoulderManual_scapula_acromion_AAx = obj.shoulderManual.scapula.acromion.AAx;
                 Result.shoulderManual_scapula_acromion_AAy = obj.shoulderManual.scapula.acromion.AAy;
                 Result.shoulderManual_scapula_acromion_AAz = obj.shoulderManual.scapula.acromion.AAz;
                 Result.shoulderManual_scapula_acromion_ACx = obj.shoulderManual.scapula.acromion.ACx;
                 Result.shoulderManual_scapula_acromion_ACy = obj.shoulderManual.scapula.acromion.ACy;
                 Result.shoulderManual_scapula_acromion_ACz = obj.shoulderManual.scapula.acromion.ACz;
             end   
             if isempty(obj.shoulderManual.humerus.center)
                 Result.shoulderManual_humerus_center = nan(1,3);
                 Result.shoulderManual_humerus_radius = nan;
             else
                 Result.shoulderManual_humerus_center = obj.shoulderManual.humerus.center;
                 Result.shoulderManual_humerus_radius = obj.shoulderManual.humerus.radius;
             end    
 
             if isempty(obj.shoulderManual.humerus.jointRadius)
                 Result.shoulderManual_humerus_jointRadius     = nan;
             else
                 Result.shoulderManual_humerus_jointRadius     = obj.shoulderManual.humerus.jointRadius;
             end
             if isempty(obj.shoulderManual.humerus.SHSAngle)
                 Result.shoulderManual_humerus_SHSAngle = nan;
                 Result.shoulderManual_humerus_SHSPA = nan;
                 Result.shoulderManual_humerus_SHSIS = nan;
                 Result.shoulderManual_humerus_SHSAmpl = nan;
                 Result.shoulderManual_humerus_SHSOrient = nan;
                 Result.shoulderManual_humerus_GHSAmpl = nan;
                 Result.shoulderManual_humerus_GHSOrient = nan;
             else    
                 Result.shoulderManual_humerus_SHSAngle = obj.shoulderManual.humerus.SHSAngle;
                 Result.shoulderManual_humerus_SHSPA = obj.shoulderManual.humerus.SHSPA;
                 Result.shoulderManual_humerus_SHSIS = obj.shoulderManual.humerus.SHSIS;
                 Result.shoulderManual_humerus_SHSAmpl = obj.shoulderManual.humerus.SHSAmpl;
                 Result.shoulderManual_humerus_SHSOrient = obj.shoulderManual.humerus.SHSOrient;
                 Result.shoulderManual_humerus_GHSAmpl = obj.shoulderManual.humerus.GHSAmpl;
                 Result.shoulderManual_humerus_GHSOrient = obj.shoulderManual.humerus.GHSOrient;
             end
-            if  isempty(obj.shoulderManual.muscles.list)
+            if  isempty(obj.shoulderManual.rotatorCuff)
                 Result.shoulderManual_muscles_IS_segmentationSet = char(nan);
                 Result.shoulderManual_muscles_IS_sliceName = char(nan);
                 Result.shoulderManual_muscles_IS_PCSA = "";
                 Result.shoulderManual_muscles_IS_atrophy = "";
                 Result.shoulderManual_muscles_IS_fat = "";
                 Result.shoulderManual_muscles_IS_osteochondroma = "";
                 Result.shoulderManual_muscles_IS_degeneration = "";
                 Result.shoulderManual_muscles_SS_segmentationSet = char(nan);
                 Result.shoulderManual_muscles_SS_sliceName = char(nan);
                 Result.shoulderManual_muscles_SS_PCSA = "";
                 Result.shoulderManual_muscles_SS_atrophy = "";
                 Result.shoulderManual_muscles_SS_fat = "";
                 Result.shoulderManual_muscles_SS_osteochondroma = "";
                 Result.shoulderManual_muscles_SS_degeneration = "";
                 Result.shoulderManual_muscles_SC_segmentationSet = char(nan);
                 Result.shoulderManual_muscles_SC_sliceName = char(nan);
                 Result.shoulderManual_muscles_SC_PCSA = "";
                 Result.shoulderManual_muscles_SC_atrophy = "";
                 Result.shoulderManual_muscles_SC_fat = "";
                 Result.shoulderManual_muscles_SC_osteochondroma = "";
                 Result.shoulderManual_muscles_SC_degeneration = "";
                 Result.shoulderManual_muscles_TM_segmentationSet = char(nan);
                 Result.shoulderManual_muscles_TM_sliceName = char(nan);
                 Result.shoulderManual_muscles_TM_PCSA = "";
                 Result.shoulderManual_muscles_TM_atrophy = "";
                 Result.shoulderManual_muscles_TM_fat = "";
                 Result.shoulderManual_muscles_TM_osteochondroma = "";
                 Result.shoulderManual_muscles_TM_degeneration = "";   
             else
                 try
-                    Result.shoulderManual_muscles_IS_segmentationSet = obj.shoulderManual.muscles.list.IS.segmentationSet;
-                    Result.shoulderManual_muscles_IS_sliceName = obj.shoulderManual.muscles.list.IS.sliceName;
-                    Result.shoulderManual_muscles_IS_PCSA = obj.shoulderManual.muscles.list.IS.PCSA;
-                    Result.shoulderManual_muscles_IS_atrophy = obj.shoulderManual.muscles.list.IS.atrophy;
-                    Result.shoulderManual_muscles_IS_fat = obj.shoulderManual.muscles.list.IS.fat;
-                    Result.shoulderManual_muscles_IS_osteochondroma = obj.shoulderManual.muscles.list.IS.osteochondroma;
-                    Result.shoulderManual_muscles_IS_degeneration = obj.shoulderManual.muscles.list.IS.degeneration;
-                    Result.shoulderManual_muscles_SS_segmentationSet = obj.shoulderManual.muscles.list.SS.segmentationSet;
-                    Result.shoulderManual_muscles_SS_sliceName = obj.shoulderManual.muscles.list.SS.sliceName;
-                    Result.shoulderManual_muscles_SS_PCSA = obj.shoulderManual.muscles.list.SS.PCSA;
-                    if isempty(obj.shoulderManual.muscles.list.SS.atrophy)
+                    Result.shoulderManual_muscles_IS_segmentationSet = obj.shoulderManual.rotatorCuff.IS.segmentationSet;
+                    Result.shoulderManual_muscles_IS_sliceName = obj.shoulderManual.rotatorCuff.IS.sliceName;
+                    Result.shoulderManual_muscles_IS_PCSA = obj.shoulderManual.rotatorCuff.IS.PCSA;
+                    Result.shoulderManual_muscles_IS_atrophy = obj.shoulderManual.rotatorCuff.IS.atrophy;
+                    Result.shoulderManual_muscles_IS_fat = obj.shoulderManual.rotatorCuff.IS.fat;
+                    Result.shoulderManual_muscles_IS_osteochondroma = obj.shoulderManual.rotatorCuff.IS.osteochondroma;
+                    Result.shoulderManual_muscles_IS_degeneration = obj.shoulderManual.rotatorCuff.IS.degeneration;
+                    Result.shoulderManual_muscles_SS_segmentationSet = obj.shoulderManual.rotatorCuff.SS.segmentationSet;
+                    Result.shoulderManual_muscles_SS_sliceName = obj.shoulderManual.rotatorCuff.SS.sliceName;
+                    Result.shoulderManual_muscles_SS_PCSA = obj.shoulderManual.rotatorCuff.SS.PCSA;
+                    if isempty(obj.shoulderManual.rotatorCuff.SS.atrophy)
                         Result.shoulderManual_muscles_SS_atrophy = nan;
                     else
-                        Result.shoulderManual_muscles_SS_atrophy = obj.shoulderManual.muscles.list.SS.atrophy;
+                        Result.shoulderManual_muscles_SS_atrophy = obj.shoulderManual.rotatorCuff.SS.atrophy;
                     end  
-                    if isempty(obj.shoulderManual.muscles.list.SS.fat)
+                    if isempty(obj.shoulderManual.rotatorCuff.SS.fat)
                         Result.shoulderManual_muscles_SS_fat = nan;
                     else
-                        Result.shoulderManual_muscles_SS_fat = obj.shoulderManual.muscles.list.SS.fat;
+                        Result.shoulderManual_muscles_SS_fat = obj.shoulderManual.rotatorCuff.SS.fat;
                     end
-                    if isempty(obj.shoulderManual.muscles.list.SS.osteochondroma)
+                    if isempty(obj.shoulderManual.rotatorCuff.SS.osteochondroma)
                         Result.shoulderManual_muscles_SS_osteochondroma = nan;
                     else
-                        Result.shoulderManual_muscles_SS_osteochondroma = obj.shoulderManual.muscles.list.SS.osteochondroma;
+                        Result.shoulderManual_muscles_SS_osteochondroma = obj.shoulderManual.rotatorCuff.SS.osteochondroma;
                     end
-                    if isempty(obj.shoulderManual.muscles.list.SS.degeneration)
+                    if isempty(obj.shoulderManual.rotatorCuff.SS.degeneration)
                         Result.shoulderManual_muscles_SS_degeneration = nan;
                     else
-                        Result.shoulderManual_muscles_SS_degeneration = obj.shoulderManual.muscles.list.SS.degeneration;
+                        Result.shoulderManual_muscles_SS_degeneration = obj.shoulderManual.rotatorCuff.SS.degeneration;
                     end
-                    Result.shoulderManual_muscles_SC_segmentationSet = obj.shoulderManual.muscles.list.SC.segmentationSet;
-                    Result.shoulderManual_muscles_SC_sliceName = obj.shoulderManual.muscles.list.SC.sliceName;
-                    Result.shoulderManual_muscles_SC_PCSA = obj.shoulderManual.muscles.list.SC.PCSA;
-                    Result.shoulderManual_muscles_SC_atrophy = obj.shoulderManual.muscles.list.SC.atrophy;
-                    Result.shoulderManual_muscles_SC_fat = obj.shoulderManual.muscles.list.SC.fat;
-                    Result.shoulderManual_muscles_SC_osteochondroma = obj.shoulderManual.muscles.list.SC.osteochondroma;
-                    Result.shoulderManual_muscles_SC_degeneration = obj.shoulderManual.muscles.list.SC.degeneration;
-                    Result.shoulderManual_muscles_TM_segmentationSet = obj.shoulderManual.muscles.list.TM.segmentationSet;
-                    Result.shoulderManual_muscles_TM_sliceName = obj.shoulderManual.muscles.list.TM.sliceName;
-                    Result.shoulderManual_muscles_TM_PCSA = obj.shoulderManual.muscles.list.TM.PCSA;
-                    if isempty(obj.shoulderManual.muscles.list.TM.atrophy)
+                    Result.shoulderManual_muscles_SC_segmentationSet = obj.shoulderManual.rotatorCuff.SC.segmentationSet;
+                    Result.shoulderManual_muscles_SC_sliceName = obj.shoulderManual.rotatorCuff.SC.sliceName;
+                    Result.shoulderManual_muscles_SC_PCSA = obj.shoulderManual.rotatorCuff.SC.PCSA;
+                    Result.shoulderManual_muscles_SC_atrophy = obj.shoulderManual.rotatorCuff.SC.atrophy;
+                    Result.shoulderManual_muscles_SC_fat = obj.shoulderManual.rotatorCuff.SC.fat;
+                    Result.shoulderManual_muscles_SC_osteochondroma = obj.shoulderManual.rotatorCuff.SC.osteochondroma;
+                    Result.shoulderManual_muscles_SC_degeneration = obj.shoulderManual.rotatorCuff.SC.degeneration;
+                    Result.shoulderManual_muscles_TM_segmentationSet = obj.shoulderManual.rotatorCuff.TM.segmentationSet;
+                    Result.shoulderManual_muscles_TM_sliceName = obj.shoulderManual.rotatorCuff.TM.sliceName;
+                    Result.shoulderManual_muscles_TM_PCSA = obj.shoulderManual.rotatorCuff.TM.PCSA;
+                    if isempty(obj.shoulderManual.rotatorCuff.TM.atrophy)
                         Result.shoulderManual_muscles_TM_atrophy     = nan;
                     else
-                        Result.shoulderManual_muscles_TM_atrophy     = obj.shoulderManual.muscles.list.TM.atrophy;
+                        Result.shoulderManual_muscles_TM_atrophy     = obj.shoulderManual.rotatorCuff.TM.atrophy;
                     end
-                    if isempty(obj.shoulderManual.muscles.list.TM.fat)
+                    if isempty(obj.shoulderManual.rotatorCuff.TM.fat)
                         Result.shoulderManual_muscles_TM_fat     = nan;
                     else
-                        Result.shoulderManual_muscles_TM_fat     = obj.shoulderManual.muscles.list.TM.fat;
+                        Result.shoulderManual_muscles_TM_fat     = obj.shoulderManual.rotatorCuff.TM.fat;
                     end
-                    if isempty(obj.shoulderManual.muscles.list.TM.osteochondroma)
+                    if isempty(obj.shoulderManual.rotatorCuff.TM.osteochondroma)
                         Result.shoulderManual_muscles_TM_osteochondroma     = nan;
                     else
-                        Result.shoulderManual_muscles_TM_osteochondroma     = obj.shoulderManual.muscles.list.TM.osteochondroma;
+                        Result.shoulderManual_muscles_TM_osteochondroma     = obj.shoulderManual.rotatorCuff.TM.osteochondroma;
                     end
-                    if isempty(obj.shoulderManual.muscles.list.TM.degeneration)
+                    if isempty(obj.shoulderManual.rotatorCuff.TM.degeneration)
                         Result.shoulderManual_muscles_TM_degeneration     = nan;
                     else
-                        Result.shoulderManual_muscles_TM_degeneration     = obj.shoulderManual.muscles.list.TM.degeneration;
+                        Result.shoulderManual_muscles_TM_degeneration     = obj.shoulderManual.rotatorCuff.TM.degeneration;
                     end
                 catch 
                     Result.shoulderManual_muscles_IS_segmentationSet = 0;
                     Result.shoulderManual_muscles_IS_sliceName = 0;
                     Result.shoulderManual_muscles_IS_PCSA = 0;
                     Result.shoulderManual_muscles_IS_atrophy = 0;
                     Result.shoulderManual_muscles_IS_fat = 0;
                     Result.shoulderManual_muscles_IS_osteochondroma = 0;
                     Result.shoulderManual_muscles_IS_degeneration = 0;
                     Result.shoulderManual_muscles_SS_segmentationSet = 0;
                     Result.shoulderManual_muscles_SS_sliceName = 0;
                     Result.shoulderManual_muscles_SS_PCSA = 0;
                     Result.shoulderManual_muscles_SS_atrophy = 0;
                     Result.shoulderManual_muscles_SS_fat = 0;
                     Result.shoulderManual_muscles_SS_osteochondroma = 0;
                     Result.shoulderManual_muscles_SS_degeneration = 0;
                     Result.shoulderManual_muscles_SC_segmentationSet = 0;
                     Result.shoulderManual_muscles_SC_sliceName = 0;
                     Result.shoulderManual_muscles_SC_PCSA = 0;
                     Result.shoulderManual_muscles_SC_atrophy = 0;
                     Result.shoulderManual_muscles_SC_fat = 0;
                     Result.shoulderManual_muscles_SC_osteochondroma = 0;
                     Result.shoulderManual_muscles_SC_degeneration = 0;
                     Result.shoulderManual_muscles_TM_segmentationSet = 0;
                     Result.shoulderManual_muscles_TM_sliceName = 0;
                     Result.shoulderManual_muscles_TM_PCSA = 0;
                     Result.shoulderManual_muscles_TM_atrophy = 0;
                     Result.shoulderManual_muscles_TM_fat = 0;
                     Result.shoulderManual_muscles_TM_osteochondroma = 0;
                     Result.shoulderManual_muscles_TM_degeneration = 0;
                 end 
             end
             Result.shoulderAuto_side                  = obj.shoulderAuto().side;
             if isempty(obj.shoulderAuto().scapula.angulusInferior)
                  Result.shoulderAuto_scapula_angulusInferior    = nan(1,3);
             else
                  Result.shoulderAuto_scapula_angulusInferior     = obj.shoulderAuto().scapula.angulusInferior;
             end
             if isempty(obj.shoulderAuto().scapula.angulusInferiorLocal)
                  Result.shoulderAuto_scapula_angulusInferiorLocal = nan(1,3);
             else
                  Result.shoulderAuto_scapula_angulusInferiorLocal = obj.shoulderAuto().scapula.angulusInferiorLocal;
             end
             if isempty(obj.shoulderAuto().scapula.trigonumSpinae)
                  Result.shoulderAuto_scapula_trigonumSpinae = nan(1,3);
             else
                  Result.shoulderAuto_scapula_trigonumSpinae = obj.shoulderAuto().scapula.trigonumSpinae;
             end
             if isempty(obj.shoulderAuto().scapula.processusCoracoideus)
                  Result.shoulderAuto_scapula_processusCoracoideus = nan(1,3);
                  Result.shoulderAuto_scapula_acromioClavicular = nan(1,3);
                  Result.shoulderAuto_scapula_angulusAcromialis = nan(1,3);
                  Result.shoulderAuto_scapula_spinoGlenoidNotch = nan(1,3);
             else
                  Result.shoulderAuto_scapula_processusCoracoideus = obj.shoulderAuto.scapula.processusCoracoideus;
                  Result.shoulderAuto_scapula_acromioClavicular = obj.shoulderAuto.scapula.acromioClavicular;
                  Result.shoulderAuto_scapula_angulusAcromialis = obj.shoulderAuto.scapula.angulusAcromialis;
                  Result.shoulderAuto_scapula_spinoGlenoidNotch = obj.shoulderAuto.scapula.spinoGlenoidNotch;
             end
             if isempty(obj.shoulderAuto.scapula.coordSys.ML)
                  Result.shoulderAuto_scapula_coordSys_ML = nan(1,3);
             else
                  Result.shoulderAuto_scapula_coordSys_ML = obj.shoulderAuto.scapula.coordSys.ML;
             end
             if isempty(obj.shoulderAuto.scapula.coordSys.PA)
                  Result.shoulderAuto_scapula_coordSys_PA = nan(1,3);
             else
                  Result.shoulderAuto_scapula_coordSys_PA = obj.shoulderAuto.scapula.coordSys.PA;
             end
             if isempty(obj.shoulderAuto.scapula.coordSys.IS)
                  Result.shoulderAuto_scapula_coordSys_IS = nan(1,3);
             else
                  Result.shoulderAuto_scapula_coordSys_IS= obj.shoulderAuto.scapula.coordSys.IS;
             end
             if isempty(obj.shoulderAuto.scapula.coordSys.origin)
                  Result.shoulderAuto_scapula_coordSys_origin = nan(1,3);
             else
                  Result.shoulderAuto_scapula_coordSys_origin = obj.shoulderAuto.scapula.coordSys.origin;
             end
             if isempty(obj.shoulderAuto.scapula.glenoid.center)
                  Result.shoulderAuto_scapula_glenoid_center = nan(1,3);
                  Result.shoulderAuto_scapula_glenoid_centerLocal_x = nan;
                  Result.shoulderAuto_scapula_glenoid_centerLocal_y = nan;
                  Result.shoulderAuto_scapula_glenoid_centerLocal_z = nan;
                  Result.shoulderAuto_scapula_glenoid_radius = nan;
                  Result.shoulderAuto_scapula_glenoid_centerLine = nan(1,3);
                  Result.shoulderAuto_scapula_glenoid_depth = nan;
                  Result.shoulderAuto_scapula_glenoid_width = nan;
                  Result.shoulderAuto_scapula_glenoid_height = nan;
                  Result.shoulderAuto_scapula_glenoid_anteroSuperiorAngle = nan;
                  Result.shoulderAuto_scapula_glenoid_versionAmplitude = nan;
                  Result.shoulderAuto_scapula_glenoid_versionOrientation = nan;
                  Result.shoulderAuto_scapula_glenoid_version = nan;
                  Result.shoulderAuto_scapula_glenoid_inclination = nan;
      
             else
                  Result.shoulderAuto_scapula_glenoid_center = obj.shoulderAuto.scapula.glenoid.center;
                  Result.shoulderAuto_scapula_glenoid_centerLocal_x = obj.shoulderAuto.scapula.glenoid.centerLocal.x;
                  Result.shoulderAuto_scapula_glenoid_centerLocal_y = obj.shoulderAuto.scapula.glenoid.centerLocal.y;
                  Result.shoulderAuto_scapula_glenoid_centerLocal_z = obj.shoulderAuto.scapula.glenoid.centerLocal.z;
                  Result.shoulderAuto_scapula_glenoid_radius = obj.shoulderAuto.scapula.glenoid.radius;
                  Result.shoulderAuto_scapula_glenoid_centerLine = obj.shoulderAuto.scapula.glenoid.centerLine;
                  Result.shoulderAuto_scapula_glenoid_depth = obj.shoulderAuto.scapula.glenoid.depth;
                  Result.shoulderAuto_scapula_glenoid_width = obj.shoulderAuto.scapula.glenoid.width;
                  Result.shoulderAuto_scapula_glenoid_height = obj.shoulderAuto.scapula.glenoid.height;
                  Result.shoulderAuto_scapula_glenoid_anteroSuperiorAngle = obj.shoulderAuto.scapula.glenoid.anteroSuperiorAngle;
                  Result.shoulderAuto_scapula_glenoid_versionAmplitude = obj.shoulderAuto.scapula.glenoid.versionAmplitude;
                  Result.shoulderAuto_scapula_glenoid_versionOrientation = obj.shoulderAuto.scapula.glenoid.versionOrientation;
                  Result.shoulderAuto_scapula_glenoid_version = obj.shoulderAuto.scapula.glenoid.version;
                  Result.shoulderAuto_scapula_glenoid_inclination = obj.shoulderAuto.scapula.glenoid.inclination;
               
             end
             if isempty(obj.shoulderAuto.scapula.glenoid.density)
                 Result.shoulderAuto_scapula_glenoid_density = nan(1, 6);
             else
                 Result.shoulderAuto_scapula_glenoid_density = obj.shoulderAuto.scapula.glenoid.density;
             end   
             if isempty(obj.shoulderAuto.scapula.glenoid.comment)
                  Result.shoulderAuto_scapula_glenoid_comment     = nan;
             else
                  Result.shoulderAuto_scapula_glenoid_comment     = obj.shoulderAuto.scapula.glenoid.comment;
             end
             if isempty(obj.shoulderAuto.scapula.glenoid.fittedSphere.center)
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_center = nan(1,3);
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_radius = nan;
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_R2 = nan;
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_RMSE = nan;
             else
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_center = ...
                     obj.shoulderAuto.scapula.glenoid.fittedSphere.center;
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_radius = ...
                     obj.shoulderAuto.scapula.glenoid.fittedSphere.radius;
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_R2 = ...
                     obj.shoulderAuto.scapula.glenoid.fittedSphere.R2;
                 Result.shoulderAuto_scapula_glenoid_fittedSphere_RMSE = ...
                     obj.shoulderAuto.scapula.glenoid.fittedSphere.RMSE;
             end    
 
             if isempty(obj.shoulderAuto.scapula.glenoid.walch)
                 Result.shoulderAuto_scapula_glenoid_walch  = char(nan);
             elseif isnan(obj.shoulderAuto.scapula.glenoid.walch) 
                 Result.shoulderAuto_scapula_glenoid_walch  = char(nan);
             else
                 Result.shoulderAuto_scapula_glenoid_walch  = obj.shoulderAuto.scapula.glenoid.walch;
             end
             if isempty(obj.shoulderAuto.scapula.acromion.AI)
                 Result.shoulderAuto_scapula_acromion_AI = nan;
                 Result.shoulderAuto_scapula_acromion_CSA = nan;
                 Result.shoulderAuto_scapula_acromion_PSA = nan;
                 Result.shoulderAuto_scapula_acromion_PSL = nan;
                 Result.shoulderAuto_scapula_acromion_AAA = nan;
                 Result.shoulderAuto_scapula_acromion_AAL = nan;
                 Result.shoulderAuto_scapula_acromion_AAx = nan;
                 Result.shoulderAuto_scapula_acromion_AAy = nan;
                 Result.shoulderAuto_scapula_acromion_AAz = nan;
                 Result.shoulderAuto_scapula_acromion_ACx = nan;
                 Result.shoulderAuto_scapula_acromion_ACy = nan;
                 Result.shoulderAuto_scapula_acromion_ACz = nan;
             else
                 Result.shoulderAuto_scapula_acromion_AI = obj.shoulderAuto.scapula.acromion.AI;
                 Result.shoulderAuto_scapula_acromion_CSA = obj.shoulderAuto.scapula.acromion.CSA;
                 Result.shoulderAuto_scapula_acromion_PSA = obj.shoulderAuto.scapula.acromion.PSA;
                 Result.shoulderAuto_scapula_acromion_PSL = obj.shoulderAuto.scapula.acromion.PSL;
                 Result.shoulderAuto_scapula_acromion_AAA = obj.shoulderAuto.scapula.acromion.AAA;
                 Result.shoulderAuto_scapula_acromion_AAL = obj.shoulderAuto.scapula.acromion.AAL;
                 Result.shoulderAuto_scapula_acromion_AAx = obj.shoulderAuto.scapula.acromion.AAx;
                 Result.shoulderAuto_scapula_acromion_AAy = obj.shoulderAuto.scapula.acromion.AAy;
                 Result.shoulderAuto_scapula_acromion_AAz = obj.shoulderAuto.scapula.acromion.AAz;
                 Result.shoulderAuto_scapula_acromion_ACx = obj.shoulderAuto.scapula.acromion.ACx;
                 Result.shoulderAuto_scapula_acromion_ACy = obj.shoulderAuto.scapula.acromion.ACy;
                 Result.shoulderAuto_scapula_acromion_ACz = obj.shoulderAuto.scapula.acromion.ACz;   
             end    
             if isempty(obj.shoulderAuto.scapula.glenoid.center)
                 Result.shoulderAuto_humerus_center = nan(1,3);
                 Result.shoulderAuto_humerus_radius = nan;
             else
                 Result.shoulderAuto_humerus_center = obj.shoulderAuto.scapula.glenoid.center;
                 Result.shoulderAuto_humerus_radius = obj.shoulderAuto.scapula.glenoid.radius;
             end    
             if isempty(obj.shoulderAuto.scapula.glenoid.comment)
                  Result.shoulderAuto_humerus_jointRadius     = nan;
             else
                  Result.shoulderAuto_humerus_jointRadius     = obj.shoulderAuto.scapula.glenoid.comment;
             end
             if isempty(obj.shoulderAuto.humerus.SHSAngle)
                 Result.shoulderAuto_humerus_SHSAngle = nan;
                 Result.shoulderAuto_humerus_SHSPA = nan;
                 Result.shoulderAuto_humerus_SHSIS = nan;
                 Result.shoulderAuto_humerus_SHSAmpl = nan;
                 Result.shoulderAuto_humerus_SHSOrient = nan;
                 Result.shoulderAuto_humerus_GHSAmpl = nan;
                 Result.shoulderAuto_humerus_GHSOrient = nan;
                 Result.shoulderAuto_humerus_GHSAmpl = nan;     
             else
                 Result.shoulderAuto_humerus_SHSAngle = obj.shoulderAuto.humerus.SHSAngle;
                 Result.shoulderAuto_humerus_SHSPA = obj.shoulderAuto.humerus.SHSPA;
                 Result.shoulderAuto_humerus_SHSIS = obj.shoulderAuto.humerus.SHSIS;
                 Result.shoulderAuto_humerus_SHSAmpl = obj.shoulderAuto.humerus.SHSAmpl;
                 Result.shoulderAuto_humerus_SHSOrient = obj.shoulderAuto.humerus.SHSOrient;
                 Result.shoulderAuto_humerus_GHSAmpl = obj.shoulderAuto.humerus.GHSAmpl;
                 Result.shoulderAuto_humerus_GHSOrient = obj.shoulderAuto.humerus.GHSOrient;
                 Result.shoulderAuto_humerus_GHSAmpl = obj.shoulderAuto.humerus.GHSAmpl;
             end         
             try 
-                Result.shoulderAuto_muscles_IS_segmentationSet = obj.shoulderAuto.muscles.list.IS.segmentationSet;
-                Result.shoulderAuto_muscles_IS_sliceName = obj.shoulderAuto.muscles.list.IS.sliceName;
-                Result.shoulderAuto_muscles_IS_PCSA = obj.shoulderAuto.muscles.list.IS.PCSA;
-                if isempty(obj.shoulderAuto.muscles.list.IS.atrophy)
+                Result.shoulderAuto_muscles_IS_segmentationSet = obj.shoulderAuto.rotatorCuff.IS.segmentationSet;
+                Result.shoulderAuto_muscles_IS_sliceName = obj.shoulderAuto.rotatorCuff.IS.sliceName;
+                Result.shoulderAuto_muscles_IS_PCSA = obj.shoulderAuto.rotatorCuff.IS.PCSA;
+                if isempty(obj.shoulderAuto.rotatorCuff.IS.atrophy)
                     Result.shoulderAuto_muscles_IS_atrophy = nan;
                 else
-                    Result.shoulderAuto_muscles_IS_atrophy = obj.shoulderAuto.muscles.list.IS.atrophy;
+                    Result.shoulderAuto_muscles_IS_atrophy = obj.shoulderAuto.rotatorCuff.IS.atrophy;
                 end  
-                if isempty(obj.shoulderAuto.muscles.list.IS.fat)
+                if isempty(obj.shoulderAuto.rotatorCuff.IS.fat)
                     Result.shoulderAuto_muscles_IS_fat = nan;
                 else
-                    Result.shoulderAuto_muscles_IS_fat = obj.shoulderAuto.muscles.list.IS.fat;
+                    Result.shoulderAuto_muscles_IS_fat = obj.shoulderAuto.rotatorCuff.IS.fat;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.IS.osteochondroma)
+                if isempty(obj.shoulderAuto.rotatorCuff.IS.osteochondroma)
                     Result.shoulderAuto_muscles_IS_osteochondroma = nan;
                 else
-                    Result.shoulderAuto_muscles_IS_osteochondroma = obj.shoulderAuto.muscles.list.IS.osteochondroma;
+                    Result.shoulderAuto_muscles_IS_osteochondroma = obj.shoulderAuto.rotatorCuff.IS.osteochondroma;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.IS.degeneration)
+                if isempty(obj.shoulderAuto.rotatorCuff.IS.degeneration)
                     Result.shoulderAuto_muscles_IS_degeneration = nan;
                 else
-                    Result.shoulderAuto_muscles_IS_degeneration = obj.shoulderAuto.muscles.list.IS.degeneration;
+                    Result.shoulderAuto_muscles_IS_degeneration = obj.shoulderAuto.rotatorCuff.IS.degeneration;
                 end
-                Result.shoulderAuto_muscles_SS_segmentationSet = obj.shoulderAuto.muscles.list.SS.segmentationSet;
-                Result.shoulderAuto_muscles_SS_sliceName = obj.shoulderAuto.muscles.list.SS.sliceName;
-                Result.shoulderAuto_muscles_SS_PCSA = obj.shoulderAuto.muscles.list.SS.PCSA;
-                if isempty(obj.shoulderAuto.muscles.list.SS.atrophy)
+                Result.shoulderAuto_muscles_SS_segmentationSet = obj.shoulderAuto.rotatorCuff.SS.segmentationSet;
+                Result.shoulderAuto_muscles_SS_sliceName = obj.shoulderAuto.rotatorCuff.SS.sliceName;
+                Result.shoulderAuto_muscles_SS_PCSA = obj.shoulderAuto.rotatorCuff.SS.PCSA;
+                if isempty(obj.shoulderAuto.rotatorCuff.SS.atrophy)
                     Result.shoulderAuto_muscles_SS_atrophy = nan;
                 else
-                    Result.shoulderAuto_muscles_SS_atrophy = obj.shoulderAuto.muscles.list.SS.atrophy;
+                    Result.shoulderAuto_muscles_SS_atrophy = obj.shoulderAuto.rotatorCuff.SS.atrophy;
                 end  
-                if isempty(obj.shoulderAuto.muscles.list.SS.fat)
+                if isempty(obj.shoulderAuto.rotatorCuff.SS.fat)
                     Result.shoulderAuto_muscles_SS_fat = nan;
                 else
-                    Result.shoulderAuto_muscles_SS_fat = obj.shoulderAuto.muscles.list.SS.fat;
+                    Result.shoulderAuto_muscles_SS_fat = obj.shoulderAuto.rotatorCuff.SS.fat;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.SS.osteochondroma)
+                if isempty(obj.shoulderAuto.rotatorCuff.SS.osteochondroma)
                     Result.shoulderAuto_muscles_SS_osteochondroma = nan;
                 else
-                    Result.shoulderAuto_muscles_SS_osteochondroma = obj.shoulderAuto.muscles.list.SS.osteochondroma;
+                    Result.shoulderAuto_muscles_SS_osteochondroma = obj.shoulderAuto.rotatorCuff.SS.osteochondroma;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.SS.degeneration)
+                if isempty(obj.shoulderAuto.rotatorCuff.SS.degeneration)
                     Result.shoulderAuto_muscles_SS_degeneration = nan;
                 else
-                    Result.shoulderAuto_muscles_SS_degeneration = obj.shoulderAuto.muscles.list.SS.degeneration;
+                    Result.shoulderAuto_muscles_SS_degeneration = obj.shoulderAuto.rotatorCuff.SS.degeneration;
                 end
-                Result.shoulderAuto_muscles_SC_segmentationSet = obj.shoulderAuto.muscles.list.SC.segmentationSet;
-                Result.shoulderAuto_muscles_SC_sliceName = obj.shoulderAuto.muscles.list.SC.sliceName;
-                Result.shoulderAuto_muscles_SC_PCSA = obj.shoulderAuto.muscles.list.SC.PCSA;
-                if isempty(obj.shoulderAuto.muscles.list.SC.atrophy)
+                Result.shoulderAuto_muscles_SC_segmentationSet = obj.shoulderAuto.rotatorCuff.SC.segmentationSet;
+                Result.shoulderAuto_muscles_SC_sliceName = obj.shoulderAuto.rotatorCuff.SC.sliceName;
+                Result.shoulderAuto_muscles_SC_PCSA = obj.shoulderAuto.rotatorCuff.SC.PCSA;
+                if isempty(obj.shoulderAuto.rotatorCuff.SC.atrophy)
                     Result.shoulderAuto_muscles_SC_atrophy = nan;
                 else
-                    Result.shoulderAuto_muscles_SC_atrophy = obj.shoulderAuto.muscles.list.SC.atrophy;
+                    Result.shoulderAuto_muscles_SC_atrophy = obj.shoulderAuto.rotatorCuff.SC.atrophy;
                 end  
-                if isempty(obj.shoulderAuto.muscles.list.SC.fat)
+                if isempty(obj.shoulderAuto.rotatorCuff.SC.fat)
                     Result.shoulderAuto_muscles_SC_fat = nan;
                 else
-                    Result.shoulderAuto_muscles_SC_fat = obj.shoulderAuto.muscles.list.SC.fat;
+                    Result.shoulderAuto_muscles_SC_fat = obj.shoulderAuto.rotatorCuff.SC.fat;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.SC.osteochondroma)
+                if isempty(obj.shoulderAuto.rotatorCuff.SC.osteochondroma)
                     Result.shoulderAuto_muscles_SC_osteochondroma = nan;
                 else
-                    Result.shoulderAuto_muscles_SC_osteochondroma = obj.shoulderAuto.muscles.list.SC.osteochondroma;
+                    Result.shoulderAuto_muscles_SC_osteochondroma = obj.shoulderAuto.rotatorCuff.SC.osteochondroma;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.SC.degeneration)
+                if isempty(obj.shoulderAuto.rotatorCuff.SC.degeneration)
                     Result.shoulderAuto_muscles_SC_degeneration = nan;
                 else
-                    Result.shoulderAuto_muscles_SC_degeneration = obj.shoulderAuto.muscles.list.SC.degeneration;
+                    Result.shoulderAuto_muscles_SC_degeneration = obj.shoulderAuto.rotatorCuff.SC.degeneration;
                 end
-                Result.shoulderAuto_muscles_TM_segmentationSet = obj.shoulderAuto.muscles.list.TM.segmentationSet;
-                Result.shoulderAuto_muscles_TM_sliceName = obj.shoulderAuto.muscles.list.TM.sliceName;
-                Result.shoulderAuto_muscles_TM_PCSA = obj.shoulderAuto.muscles.list.TM.PCSA;
-                if isempty(obj.shoulderAuto.muscles.list.TM.atrophy)
+                Result.shoulderAuto_muscles_TM_segmentationSet = obj.shoulderAuto.rotatorCuff.TM.segmentationSet;
+                Result.shoulderAuto_muscles_TM_sliceName = obj.shoulderAuto.rotatorCuff.TM.sliceName;
+                Result.shoulderAuto_muscles_TM_PCSA = obj.shoulderAuto.rotatorCuff.TM.PCSA;
+                if isempty(obj.shoulderAuto.rotatorCuff.TM.atrophy)
                     Result.shoulderAuto_muscles_TM_atrophy     = nan;
                 else
-                    Result.shoulderAuto_muscles_TM_atrophy     = obj.shoulderAuto.muscles.list.TM.atrophy;
+                    Result.shoulderAuto_muscles_TM_atrophy     = obj.shoulderAuto.rotatorCuff.TM.atrophy;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.TM.fat)
+                if isempty(obj.shoulderAuto.rotatorCuff.TM.fat)
                     Result.shoulderAuto_muscles_TM_fat     = nan;
                 else
-                    Result.shoulderAuto_muscles_TM_fat     = obj.shoulderAuto.muscles.list.TM.fat;
+                    Result.shoulderAuto_muscles_TM_fat     = obj.shoulderAuto.rotatorCuff.TM.fat;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.TM.osteochondroma)
+                if isempty(obj.shoulderAuto.rotatorCuff.TM.osteochondroma)
                     Result.shoulderAuto_muscles_TM_osteochondroma     = nan;
                 else
-                    Result.shoulderAuto_muscles_TM_osteochondroma     = obj.shoulderAuto.muscles.list.TM.osteochondroma;
+                    Result.shoulderAuto_muscles_TM_osteochondroma     = obj.shoulderAuto.rotatorCuff.TM.osteochondroma;
                 end
-                if isempty(obj.shoulderAuto.muscles.list.TM.degeneration)
+                if isempty(obj.shoulderAuto.rotatorCuff.TM.degeneration)
                     Result.shoulderAuto_muscles_TM_degeneration     = nan;
                 else
-                    Result.shoulderAuto_muscles_TM_degeneration     = obj.shoulderAuto.muscles.list.TM.degeneration;
+                    Result.shoulderAuto_muscles_TM_degeneration     = obj.shoulderAuto.rotatorCuff.TM.degeneration;
                 end
             catch ME
                 Result.shoulderAuto_muscles_IS_segmentationSet = char(nan);
                 Result.shoulderAuto_muscles_IS_sliceName = char(nan);
                 Result.shoulderAuto_muscles_IS_PCSA = "";
                 Result.shoulderAuto_muscles_IS_atrophy = "";
                 Result.shoulderAuto_muscles_IS_fat = "";
                 Result.shoulderAuto_muscles_IS_osteochondroma = "";
                 Result.shoulderAuto_muscles_IS_degeneration = "";
                 Result.shoulderAuto_muscles_SS_segmentationSet = char(nan);
                 Result.shoulderAuto_muscles_SS_sliceName = char(nan);
                 Result.shoulderAuto_muscles_SS_PCSA = "";
                 Result.shoulderAuto_muscles_SS_atrophy = "";
                 Result.shoulderAuto_muscles_SS_fat = "";
                 Result.shoulderAuto_muscles_SS_osteochondroma = "";
                 Result.shoulderAuto_muscles_SS_degeneration = "";
                 Result.shoulderAuto_muscles_SC_segmentationSet = char(nan);
                 Result.shoulderAuto_muscles_SC_sliceName = char(nan);
                 Result.shoulderAuto_muscles_SC_PCSA = "";
                 Result.shoulderAuto_muscles_SC_atrophy = "";
                 Result.shoulderAuto_muscles_SC_fat = "";
                 Result.shoulderAuto_muscles_SC_osteochondroma = "";
                 Result.shoulderAuto_muscles_SC_degeneration = "";
                 Result.shoulderAuto_muscles_TM_segmentationSet = char(nan);
                 Result.shoulderAuto_muscles_TM_sliceName = char(nan);
                 Result.shoulderAuto_muscles_TM_PCSA = "";
                 Result.shoulderAuto_muscles_TM_atrophy = "";
                 Result.shoulderAuto_muscles_TM_fat = "";
                 Result.shoulderAuto_muscles_TM_osteochondroma = "";
                 Result.shoulderAuto_muscles_TM_degeneration = "";                
             end                     
             outputArg = struct2table(Result, 'AsArray', 1);
         end
     end
 end
\ No newline at end of file
diff --git a/ShoulderCase/@ShoulderCasePlotter/plot.m b/ShoulderCase/@ShoulderCasePlotter/plot.m
index fb26f4e..5720b7d 100644
--- a/ShoulderCase/@ShoulderCasePlotter/plot.m
+++ b/ShoulderCase/@ShoulderCasePlotter/plot.m
@@ -1,94 +1,100 @@
 function plot(obj)
   % set method colors
   autoColors = {'b','cyan'};
   manualColors = {'#D95319','#ffcc00'};
 
 
 
   % plot scapula data
   axes(obj.axesHandle('scapula'));  
   hold on;
   grid on
   axis vis3d
 
   auto = containers.Map;
   try
     auto('landmarks') = obj.SCase.shoulderAuto.scapula.plotLandmarks(autoColors{2},autoColors{2});
   end
   try
     auto('scapula surface') = obj.SCase.shoulderAuto.scapula.plotSurface(autoColors{1},autoColors{1});
   end
   try
     auto('glenoid surface') = obj.SCase.shoulderAuto.scapula.glenoid.plot();
   end
   try
     auto('coordinate system') = obj.SCase.shoulderAuto.scapula.coordSys.plot(autoColors{1},false);
   end
   
   manual = containers.Map;
   try
     manual('landmarks') = obj.SCase.shoulderManual.scapula.plotLandmarks(manualColors{2},manualColors{2});
   end
   try
     manual('scapula surface') = obj.SCase.shoulderManual.scapula.plotSurface(manualColors{1},manualColors{1});
   end
   try
     manual('glenoid surface') = obj.SCase.shoulderManual.scapula.glenoid.plot();
   end
   try
     manual('coordinate system') = obj.SCase.shoulderManual.scapula.coordSys.plot(manualColors{1},false);
   end
   
   outliersNormLimit = 10;
   difference = containers.Map;
   try
     difference('difference') = obj.SCase.plotManualAutoDifferences(outliersNormLimit);
   end
 
 
   
 
   % plot centered coordinate system
   axes(obj.axesHandle('centered coordinate system'));    
   hold on;
   grid on
   axis vis3d
   
   try
     auto('centered coordinate system') = obj.SCase.shoulderAuto.scapula.coordSys.plot(autoColors{1},true);
   end
   try
     manual('centered coordinate system') = obj.SCase.shoulderManual.scapula.coordSys.plot(manualColors{1},true);
   end
   
   
   
   % plot muscles segmentation
   axes(obj.axesHandle('auto muscles'));
   title('Auto segmentation with auto landmarks');
   hold on;
   try
-    auto('rotator cuff') = obj.SCase.shoulderAuto.muscles.plot();
+    auto('rotator cuff') = obj.SCase.shoulderAuto.rotatorCuff.plot();
   end
   axes(obj.axesHandle('scapula'))
   try
-    auto("rotator cuff") = [auto("rotator cuff") obj.SCase.shoulderAuto.muscles.plot3()];
+    auto("rotator cuff") = [auto("rotator cuff") obj.SCase.shoulderAuto.rotatorCuff.plot3()];
+  end
+  try
+    auto("rotator cuff") = [auto("rotator cuff") obj.SCase.shoulderAuto.humerus.plot()];
   end
   
   axes(obj.axesHandle('manual muscles'));
   title('Auto segmentation with manual landmarks');
   hold on;
   try
-    manual('rotator cuff') = obj.SCase.shoulderManual.muscles.plot();
+    manual('rotator cuff') = obj.SCase.shoulderManual.rotatorCuff.plot();
   end
   axes(obj.axesHandle('scapula'))
   try
-    manual("rotator cuff") = [manual("rotator cuff") obj.SCase.shoulderManual.muscles.plot3()];
+    manual("rotator cuff") = [manual("rotator cuff") obj.SCase.shoulderManual.rotatorCuff.plot3()];
+  end
+  try
+    manual("rotator cuff") = [manual("rotator cuff") obj.SCase.shoulderManual.humerus.plot()];
   end
 
 
   % store the data handles
   obj.plotHandle.auto = auto;
   obj.plotHandle.manual = manual;
   obj.plotHandle.difference = difference;
 end
\ No newline at end of file
diff --git a/ShoulderCase/@Sphere/Sphere.m b/ShoulderCase/@Sphere/Sphere.m
index e65cf5f..6ab1ec7 100644
--- a/ShoulderCase/@Sphere/Sphere.m
+++ b/ShoulderCase/@Sphere/Sphere.m
@@ -1,54 +1,54 @@
 classdef Sphere < handle
 % Simple sphere object.
 % Used to be fitted to glenoid points.
 % 
 % This class could be used elsewhere and the constructor
 % access might be widen.
 
   properties
     center = [];
     radius = [];
     residuals = [];
     R2 = [];
     RMSE = [];
   end
 
   methods
     function obj = Sphere(varargin)
       if nargin == 2
         obj.center = varargin{1};
         obj.radius = varargin{2};
       end
     end
 
     function fitTo(obj,points)
       [center,radius,residuals,R2] = fitSphere(points);
       obj.center            = center';
       obj.radius            = radius;
       obj.residuals         = residuals;
       obj.R2                = R2;
       obj.RMSE              = norm(residuals)/sqrt(length(residuals));
     end
 
     function output = isempty(obj)
       output = any(cellfun(@(field)isempty(obj.(field)),fields(obj)));
     end
 
     function output = polarTangentPoint(obj, polePoint, slicePlaneNormal,...
       polarPointOrientation)
       poleToCenter = Vector(polePoint, obj.center);
       assert(poleToCenter.norm > obj.radius, "Pole point is inside sphere and doesn't have any tangent.");
 
       poleToCenterPerpendicular = poleToCenter.cross(slicePlaneNormal).orientToward(polarPointOrientation).direction;
       tangentPointAngle = acos(obj.radius/poleToCenter.norm);
 
       polarTangentPoint = obj.center + ...
-        sin(tangentPointAngle)*poleToCenterPerpendicular + ...
-        cos(tangentPointAngle)*(-poleToCenter.direction);
+        obj.radius*sin(tangentPointAngle)*poleToCenterPerpendicular + ...
+        obj.radius*cos(tangentPointAngle)*(-poleToCenter.direction);
       
       output = polarTangentPoint;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@Vector/Vector.m b/ShoulderCase/@Vector/Vector.m
index f77d156..a9387a2 100644
--- a/ShoulderCase/@Vector/Vector.m
+++ b/ShoulderCase/@Vector/Vector.m
@@ -1,150 +1,162 @@
 classdef Vector < handle 
 
   properties 
     origin
     target
   end 
 
   methods
 
     function obj = Vector(varargin) 
       if nargin == 1
         origin = [0 0 0];
         target = varargin{1};
       elseif nargin == 2
         origin = varargin{1};
         target = varargin{2};
       else
         error("Vector constructor only accepts up to 2 arguments");
       end
       obj.origin = origin;
       obj.target = target;
     end
 
     function output = copy(obj)
       output = Vector(obj.origin, obj.target);
     end
 
     function set.origin(obj, value)
       assert(isequal(size(value),[1,3]),"Origin must be a 1x3 double array")
       obj.origin = value;
     end
 
     function set.target(obj, value)
       assert(isequal(size(value),[1,3]),"Target must be a 1x3 double array")
       obj.target = value;
     end
 
     function output = points(obj)
       output = [obj.origin; obj.target];
     end
 
     function output = vector(obj)
       output = obj.target - obj.origin;
     end
 
     function output = norm(obj)
       output = norm(obj.vector);
     end
 
     function output = direction(obj)
       output = obj.vector / obj.norm;
     end
 
     function output = orientToward(obj, orientationElement)
       if isequal(class(orientationElement),"Vector")
         orientation = orientationElement.vector;
       else
         orientation = orientationElement;
       end
 
       if dot(obj.vector, orientation) >= 0
         output = obj.copy();
       else
         output = -obj.copy();
       end
     end
 
     function output = project(obj, projectedElement)
       % projectedElement must be either a point in 3D space (1x3 double array)
       % or another Vector.
       if isequal(class(projectedElement),"Vector")
         projectedVector = projectedElement;
       elseif isequal(size(projectedElement),[1,3])
         projectedVector = Vector(obj.origin, projectedElement);
       end
       output = Vector(obj.origin, obj.origin + obj.direction*dot(projectedVector.vector,obj.direction));
     end
     
     function output = orthogonalComplementToPoint(obj, point)
       output = -obj.project(point) + Vector(obj.origin, point);
     end
 
-    function plot(obj)
+    function output = plot(obj, varargin)
       points = obj.points;
-      plot3(points(:,1), points(:,2), points(:,3));
+      plotHandle = plot3(points(:,1), points(:,2), points(:,3), varargin{:});
+      output = plotHandle;
     end
 
     function output = uminus(obj)
       output = Vector(obj.target, obj.origin);
     end
 
     function output = plus(obj, addedVector)
       if isequal(class(addedVector),"Vector")
         output = Vector(obj.origin, obj.target + addedVector.vector);
       elseif isequal(size(addedVector),[1,3])
         output = Vector(obj.origin, obj.target + addedVector);
       end
     end
     
     function output = minus(obj, substractedVector)
       if isequal(class(substractedVector),"Vector")
         output = obj + -substractedVector;
       elseif isequal(size(substractedVector),[1,3])
         output = Vector(obj.origin, obj.target - substractedVector);
       end
     end
 
-    function output = mtimes(left,right)
+    function output = mtimes(left, right)
       if isequal(class(left),"Vector")
         vectorObject = left;
         coefficient = right;
       else
         vectorObject = right;
         coefficient = left;
       end
       output = Vector(vectorObject.origin, vectorObject.origin + coefficient*vectorObject.vector);
     end
 
+    function output = mrdivide(left, right)
+      if isequal(class(left),"Vector")
+        vectorObject = left;
+        coefficient = right;
+      else
+        vectorObject = right;
+        coefficient = left;
+      end
+      output = Vector(vectorObject.origin, vectorObject.origin + vectorObject.vector/coefficient);
+    end
+
     function output = dot(obj, inputVector)
       if isequal(class(inputVector),"Vector")
         array = inputVector.vector;
       else
         array = inputVector;
       end
       output = dot(obj.vector, array);
     end
 
     function output = cross(obj, inputVector)
       if isequal(class(inputVector),"Vector")
         array = inputVector.vector;
       else
         array = inputVector;
       end
       output = Vector(obj.origin, obj.origin + cross(obj.vector, array));
     end
 
     function output = angle(obj, inputVector)
       if isequal(class(inputVector),"Vector")
         array = inputVector.vector;
       else
         array = inputVector;
       end
       rotation = vrrotvec(obj.vector, array);
       output = rotation(4);
     end
 
 
   end
 
 end
\ No newline at end of file
diff --git a/measureSCase.m b/measureSCase.m
index ca4171a..947895e 100755
--- a/measureSCase.m
+++ b/measureSCase.m
@@ -1,689 +1,689 @@
 function [status,message] = measureSCase(varargin)
 % MEASURESCASE Update the entire SCase DB with anatomical measurements.
 
 % The function should be run from the server (lbovenus), from the user account containing the git repository,
 % with the following shell command
 % matlab -nodisplay -nosplash -r "measureSCase;quit"
 % For long calculations, it can be run in the background using tmux, followed
 % by Ctrl+b d to detach the session and keep the process running.
 
 % Progress can be checked through log file with following shell command
 % tail -f log/measureSCase.log
 
 % It take ~30 min for 700 cases when run from server.
 % After run from server, the output file
 % /shoulder/data/Excel/xlsFromMatlab/anatomy.csv
 % must be open from Excel and saved
 % as xls at the same location. Excel file
 % /shoulder/data/Excel/ShoulderDataBase.xlsx should also be open, updated, and saved.
 
 % The script measureSCase was replacing 3 previous functions
 % 1) rebuildDatabaseScapulaMeasurements.m --> list of all SCases --> execute scapula_measure.m
 % 2) scapula_measure.m --> execute scapula_calculation and save results in Excel and SQL
 % 3) scapula_calculation.m --> perform anatomical analysis --> results variable
 
 
 % Input arguments can be one or several among:
 %     {Empty,'glenoidDensity', 'update', 'musclesDegeneration' '[NP][1-999]', '[NP]'}
 %
 %     If there is no argument the function will measure every new case.
 %     The 'N' or 'P' given as an argument will load all the Normal or Pathological cases.
 %
 %     The measurements are divided into three measurements: "morphology", "glenoidDensity",
 %     and "musclesDegeneration". The default measurement is the "morphology" one, the
 %     two others has to be given in arguments in order to be done.%
 %
 %     Examples:
 %              measureSCase()
 %                 is a valid instruction where on every new case is done the
 %                 "morphology" measurement.
 %              measureSCase('glenoidDensity','musclesDegeneration')
 %                 is a valid instruction where on every cases are done the
 %                 "morphology","glenoidDensity", and "musclesDegeneration" measurements
 %                 if one of these measurements is missing.
 %              measureSCase('P400','N520')
 %                 is a valid instruction where on the cases 'P400' and 'N520' is
 %                 done the "morphology" measurement  only if they have not already
 %                 been measured yet (new cases).
 %              measureSCase('P400','N520','update')
 %                 is a valid instruction where on the cases 'P400' and 'N520' is
 %                 done the "morphology" measurement whether they have already been
 %                 measured yet or not.
 %              measureSCase('glenoidDensity','P400','N520','update','P')
 %                 is a valid instruction where on all the 'P' cases and the case N520 
 %                 are done the "morphology" and "glenoidDensity" measurements whether 
 %                 they have been  measured yet or not.
 
 % Output
 %    File /shoulder/data/Excel/xlsFromMatlab/anatomy.csv
 %    File /shoulder/data/Excel/xlsFromMatlab/anatomy.xls (only windows)
 %    File /home/shoulder/data/matlab/SCaseDB.mat
 %    File {SCaseId}.mat in each {SCaseId}/{CT}/matlab directory
 %    logs in /log/measureSCase.log
 
 % Authors: Alexandre Terrier, EPFL-LBO
 %          Matthieu Boubat, EPFL-LBO
 % Creation date: 2018-07-01
 % Revision date: 2020-06-09
 
 % TO DO:
 %  Read first last anatomy.csv and/or SCaseDB.mat
 %  Fill with patient and clinical data, from excel or MySQL
 %  use argument to redo or update
 %  Save csv within the main loop (not sure it's a good idea)
 %  Add more message in log
 %  Add more test to assess validity of data
 %  Update MySQL
 %  Rename all objects with a capital
 %  Make it faster by not loading scapula surface if already in SCaseDB.mat
 %  Add argument 'load' to force reload files even if they are alerady in DB
 %  Add argument 'measure' to force re-measuring event if measurement is already in DB
 %  Add "scapula_addmeasure.m" --> should be coded in ShoulderCase class
 
 startTime = datetime('now');
 
 %% Add environment to path
 addpath(genpath(pwd));
 
 %% Open log file
 logFileID = openLogFile('measureSCase.log');
 
 %% Set the data directory from the configuration file config.txt
 dataDir = ConfigFileExtractor.getVariable('dataDir');
 
 %% Location of the XLS ShoulderDatabase
 xlsDir    = [dataDir '/Excel/xlsFromMatlab'];
 matlabDir = [dataDir '/matlab'];
 
 %% Instanciate case loader
 database = ShoulderCaseLoader();
 
 %% Get list of all SCase
 % Get list of SCases from varargin
 
 % Delault value for varargin
 calcGlenoidDensity = false;
 calcMusclesDegeneration = false;
 updateCalculations = false;
 specificCases = true;
 
 fprintf(logFileID, '\n\nList of SCase');
 
 allCases = database.getAllCasesID();
 assert(not(isempty(allCases)),'The database is empty, check that the provided data path is valid.');
 requestedCases = {};
 
 for argument = 1:nargin
   switch lower(varargin{argument})                                              % Test every argument in varargin.
     case 'glenoiddensity'
       calcGlenoidDensity = true;
     case 'musclesdegeneration'
       calcMusclesDegeneration = true;
     case 'update'
       updateCalculations = true;
     case regexp(lower(varargin{argument}),'[p]','match')                        % Detect if argument is 'P'.
       requestedCases = [requestedCases database.getAllPathologicalCasesIDs()];
     case regexp(lower(varargin{argument}),'[n]','match')                        % Detect if argument is 'N'.
       requestedCases = [requestedCases database.getAllNormalCasesIDs()];
     case regexp(lower(varargin{argument}),'[np][1-9][0-9]{0,2}','match')        % Test for a correct argument [NP][1-999]
       if database.containsCase(varargin{argument})
         requestedCases = [requestedCases varargin{argument}];
       end
     otherwise
       warning(['\n "%s" is not a valid argument and has not been added to the'...
                ' list of cases.\n measureSCase arguments format must be "[npNP]",'...
                ' "[npNP][1-999]", "GlenoidDensity", or "update"'],varargin{argument})
   end
 end
 
 SCaseList = allCases(find(ismember(allCases,requestedCases)));    % Remove cases that appears multiple times 
 if isempty(SCaseList)
   disp('No valid ShoulderCase have been given in argument. All the cases found in the database will be measured.');
   SCaseList = allCases;
 end
 
 %% Load current xls database (for clinical data)
 fprintf(logFileID, '\nLoad xls database\n\n');
 addpath('XLS_MySQL_DB');
 filename = [dataDir '/Excel/ShoulderDataBase.xlsx'];
 excelRaw = rawFromExcel(filename); % cell array
 
 excelSCaseID = excelRaw(2:end, 1);
 excelDiagnosis = excelRaw(2:end, 4);
 excelPatientHeight = excelRaw(2:end, 23);
 excelGlenoidWalch = excelRaw(2:end, 55);
 
 
 % Lines below adapted from XLS_MySQL_DB/MainExcel2XLS
 % Excel.diagnosisList = diagnosisListFromExcel(Excel.Raw);
 % Excel.treatmentList = treatmentListFromExcel(Excel.Raw);
 % ExcelData.Patient = patientFromExcel(excelRaw); % Structure with patient data
 % Excel.shoulder = shoulderFromExcel(Excel.patient, Excel.Raw);
 % [Excel.SCase, Excel.diagnosis, Excel.treatment, Excel.outcome, Excel.study] = SCaseFromExcel(...
 %     Excel.shoulder, Excel.patient, Excel.diagnosisList, Excel.treatmentList, Excel.Raw);
 % fprintf(logFileID, ': OK');
 
 %% Add path to ShoulderCase class
 addpath('ShoulderCase');
 
 %% Instance of a ShoulderCase object
 if (exist('SCase','var') > 0)
     clear SCase; % Check for delete method
 end
 SCaseDB = {};
 SCase = [];
 
 
 % %%%% definition of analysis variables
 % %
 % %    % manual
 % %
 % ok = {};
 % amira = {};
 % scapulaLoad = {};
 % scapulaCoord = {};
 % degeneration = {};
 % glenoidLoad = {};
 % glenoidMorphology = {};
 % glenoidDensity = {};
 % humerusLoad = {};
 % humerusSubluxation = {};
 % acromionMorphology = {};
 % %
 % %    % auto
 % %
 % autoOk = {};
 % autoScapulaLoad = {};
 % autoScapulaCoord = {};
 % %autoDegeneration = {};
 % autoGlenoidLoad = {};
 % autoGlenoidMorphology = {};
 % autoGlenoidDensity = {};
 % %autoHumerusLoad = {};
 % %autoHumerusSubluxation = {};
 % autoAcromionMorphology = {};
 % %
 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
 %% Start loop over SCases in SCaseList
 nSCase = length(SCaseList); % Number of SCases
 for iSCaseID = 1:nSCase
   SCaseID = SCaseList{iSCaseID};
 
   if not(updateCalculations)
     try
       SCase = database.loadCase(SCaseID);
     catch ME
       SCase = database.createEmptyCase(SCaseID);
       warning(ME.message);
     end
   else
     SCase = database.createEmptyCase(SCaseID);
   end
 
   % newSCase = ShoulderCase(SCaseID);
   % SCase = newSCase;
   % SCase.dataPath = dataDir; % Set dataDir for SCase
   manualMeasurementAvailable = not(isempty(dir(SCase.dataAmiraPath))); % Set data path of this SCase
   autoMeasurementsAvailable = not(isempty(dir([SCase.dataMatlabPath '/scapulaLandmarksAuto*.mat'])));
 
   percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f');
   fprintf(logFileID, ['\n___________________________________\n|\n| SCaseID: ' SCaseID ' (' percentProgress '%%)']);
 
 
 
   %% The following section is a temporary solution to the arguments handling.
   %% It splits the measurements process into six parts (shoulderMorphology,
   %% glenoidDensity, musclesDegeneration for both manual and auto cases).
   %% Each part are independently 'ordered' based on some conditions.
 
   orderMorphologyMeasurementManual = false;
   orderDensityMeasurementManual = false;
   orderDegenerationMeasurementManual = false;
 
   orderMorphologyMeasurementAuto = false;
   orderDensityMeasurementAuto = false;
   orderDegenerationMeasurementAuto = false;
 
   if or(not(exist([SCase.dataMatlabPath '/SCase.mat'],'file')), updateCalculations)
     orderMorphologyMeasurementManual = true;
     orderMorphologyMeasurementAuto = true;
     if (calcGlenoidDensity)
       orderDensityMeasurementManual = true;
       orderDensityMeasurementAuto = true;
     end
     if (calcMusclesDegeneration)
       orderDegenerationMeasurementManual = true;
       orderDegenerationMeasurementAuto = true;
     end
   else
     loadedSCase = database.loadCase(SCase.id);
-    musclesManual = loadedSCase.shoulderManual.muscles;
-    musclesAuto = loadedSCase.shoulderAuto.muscles;
+    rotatorCuffManual = loadedSCase.shoulderManual.rotatorCuff;
+    rotatorCuffAuto = loadedSCase.shoulderAuto.rotatorCuff;
     if and(calcGlenoidDensity, isempty(loadedSCase.shoulderManual.scapula.glenoid.density))
       orderMorphologyMeasurementManual = true;
       orderDensityMeasurementManual = true;
       SCase = loadedSCase;
     end
     if and(calcGlenoidDensity, isempty(loadedSCase.shoulderAuto.scapula.glenoid.density))
       orderMorphologyMeasurementAuto = true;
       orderDensityMeasurementAuto = true;
       SCase = loadedSCase;
     end
-    if and(calcMusclesDegeneration, not(all(table2array(musclesManual.summary(:,{'PCSA'})))) )
+    if and(calcMusclesDegeneration, any(isnan(rotatorCuffManual.summary.PCSA)) )
       orderMorphologyMeasurementManual = true;
       orderDegenerationMeasurementManual = true;
       SCase = loadedSCase;
     end
-    if and(calcMusclesDegeneration, not(all(table2array(musclesAuto.summary(:,{'PCSA'})))) )
+    if and(calcMusclesDegeneration, any(isnan(rotatorCuffAuto.summary.PCSA)) )
       orderMorphologyMeasurementAuto = true;
       orderDegenerationMeasurementAuto = true;
       SCase = loadedSCase;
     end
   end
 
   %%
   %% The section above (and basically all the measureSCase function) will be removed
   %% in the upcoming implementation of ShoulderCase measurements.
   %% Matthieu Boubat, 06.03.2020
   %%
 
 
 
   % There are 3 parts within this SCase loop:
   % 1) Load & analyses manual data from amira
   % 2) Load & analyses auto   data from matlab (Statistical Shape Model)
   % 3) Load clinical data from Excel database, and set them to SCase
   % 4) Save SCase in file SCaseID/matlab/SCase.mat
 
 
 
   % 1) Load & analyses manual data from amira
 
   if manualMeasurementAvailable
     tic;
     fprintf(logFileID, '\n|  Segmentation manual ');
 
     try
       if orderMorphologyMeasurementManual
         measureShoulderMorphology(SCase.shoulderManual,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     try
       if orderDensityMeasurementManual
         measureGlenoidDensity(SCase.shoulderManual,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     try
       if orderDegenerationMeasurementManual
         measureMusclesDegeneration(SCase.shoulderManual,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     fprintf(logFileID, '\n|  Elapsed time (s): %s\n|', num2str(toc));
   end
 
 
 
 
   % 2) Load & analyses auto data from matlab
 
   if autoMeasurementsAvailable
     tic;
     fprintf(logFileID, '\n|  Segmentation auto ');
 
     try
       if orderMorphologyMeasurementAuto
         measureShoulderMorphology(SCase.shoulderAuto,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     try
       if orderDensityMeasurementAuto
         measureGlenoidDensity(SCase.shoulderAuto,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     try
       if orderDegenerationMeasurementAuto
         measureMusclesDegeneration(SCase.shoulderAuto,logFileID);
       end
     catch ME
       fprintf(logFileID,ME.message);
     end
 
     fprintf(logFileID, '\n|  Elapsed time (s): %s\n|', num2str(toc));
   end
 
 
 
   % 3) Set clinical data from Excel database to SCase
   fprintf(logFileID, '\n|\n|  Set clinical data from Excel');
   % Get idx of excelRaw for SCaseID
   % Use cell2struct when dots are removed from excel headers
   try
     idx = find(strcmp(excelRaw, SCaseID)); % Index of SCaseID in excelRaw
     SCase.diagnosis = excelRaw{idx,4};
     SCase.treatment = excelRaw{idx,9};
     SCase.patient.gender = excelRaw{idx,19};
     SCase.patient.age = excelRaw{idx,21};
     SCase.patient.height = excelRaw{idx,23};
     SCase.patient.weight = excelRaw{idx,24};
     SCase.patient.BMI = excelRaw{idx,25};
     SCase.shoulderAuto.scapula.glenoid.walch = excelRaw{idx,55};
     SCase.shoulderManual.scapula.glenoid.walch = excelRaw{idx,55};
     fprintf(logFileID, ': OK');
   catch
     fprintf(logFileID, ': ERROR');
   end
 
   % 4) Save SCase
   try
     fprintf(logFileID, '\n| Save SCase.mat');
     SCase.saveMatlab;
     fprintf(logFileID, ': OK');
   catch
     fprintf(logFileID, ': ERROR');
   end
   fprintf(logFileID, '\n|___________________________________\n\n');
   % 4) Load & analyses auto data from matlab (Statistical Shape Model)
   % Load scapula surface and landmarks
   % save('measureSCase_analysis.mat',...
   %      'ok','amira','scapulaLoad','scapulaCoord',...
   %      'degeneration','glenoidLoad','glenoidMorphology','glenoidDensity',...
   %      'humerusLoad','humerusSubluxation','acromionMorphology',...
   %      'autoOk','autoScapulaLoad','autoScapulaCoord','autoGlenoidLoad',...
   %      'autoGlenoidMorphology','autoGlenoidDensity','autoAcromionMorphology')
 
 end % End of loop on SCaseList
 
 %% Write csv file
 % This might be a function (SCase.csvSave()). The input would be a filename and a structure
 % data
 
 
 % Replace header and data by structure. Currently not working
 %{
 txtFilename = [xlsDir, '/anatomy.txt']; % Name of the txt file
 
 DataStruc = struc(...
     'SCase_id',              SCase.id, ...
     'shoulder_side',         SCase.shoulder.side,...
     'glenoid_radius',        SCase.shoulder.scapula.glenoid.radius,...
     'glenoid_sphereRMSE',    SCase.shoulder.scapula.glenoid.sphereRMSE,...
     'glenoid_depth',         SCase.shoulder.scapula.glenoid.depth,...
     'glenoid_width',         SCase.shoulder.scapula.glenoid.width,...
     'glenoid_height',        SCase.shoulder.scapula.glenoid.height,...
     'glenoid_centerPA',      SCase.shoulder.scapula.glenoid.centerLocal(1),...
     'glenoid_centerIS',      SCase.shoulder.scapula.glenoid.centerLocal(2),...
     'glenoid_centerML',      SCase.shoulder.scapula.glenoid.centerLocal(3),...
     'glenoid_versionAmpl',   SCase.shoulder.scapula.glenoid.versionAmpl,...
     'glenoid_versionOrient', SCase.shoulder.scapula.glenoid.versionOrient,...
     'glenoid_version',       SCase.shoulder.scapula.glenoid.version,...
     'glenoid_inclination',   SCase.shoulder.scapula.glenoid.inclination,...
     'humerus_jointRadius',   SCase.shoulder.humerus.jointRadius,...
     'humerus_headRadius',    SCase.shoulder.humerus.radius,...
     'humerus_GHSAmpl',       SCase.shoulder.humerus.GHSAmpl,...
     'humerus_GHSOrient',     SCase.shoulder.humerus.GHSOrient,...
     'humerus_SHSAmpl',       SCase.shoulder.humerus.SHSAmpl,...
     'humerus_SHSOrient',     SCase.shoulder.humerus.SHSOrient,...
     'humerus_SHSAngle',      SCase.shoulder.humerus.SHSAngle,...
     'humerus_SHSPA',         SCase.shoulder.humerus.SHSPA,...
     'humerus_SHSIS',         SCase.shoulder.humerus.SHSIS,...
     'acromion_AI',           SCase.shoulder.scapula.acromion.AI,...
     'acromion_CSA',          SCase.shoulder.scapula.acromion.CSA,...
     'acromion_PS',           SCase.shoulder.scapula.acromion.PS...
     );
 DataTable = strct2table(Data);
 writetable(DataTable,filename);
 %}
 
 % Header of the csv file
 
 dataHeader = [...
     'SCase_id,' ...
     'shoulder_side,' ...
     'glenoid_radius,' ...
     'glenoid_sphereRMSE,' ...
     'glenoid_depth,' ...
     'glenoid_width,' ...
     'glenoid_height,' ...
     'glenoid_centerPA,' ...
     'glenoid_centerIS,' ...
     'glenoid_centerML,' ...
     'glenoid_versionAmpl,' ...
     'glenoid_versionOrient,' ...
     'glenoid_version,' ...
     'glenoid_inclination,' ...
     'humerus_jointRadius,' ...
     'humerus_headRadius,' ...
     'humerus_GHSAmpl,' ...
     'humerus_GHSOrient,' ...
     'humerus_SHSAmpl,' ...
     'humerus_SHSOrient,' ...
     'humerus_SHSAngle,' ...
     'humerus_SHSPA,' ...
     'humerus_SHSIS,' ...
     'acromion_AI,' ...
     'acromion_CSA,' ...
     'acromion_PSA,'...
     'acromion_AAA\n'...
     ];
 
 
 fprintf(logFileID, '\n\nSave anatomy.csv file');
 fid = fopen([xlsDir, '/anatomy.csv'],'w');      %Last csv is deleted and a new one, updated is being created
 fprintf(fid,dataHeader);
 fclose(fid);
 
 % The following lines re-create the SCaseDB.mat and contruct the anatomy.csv file
 SCaseDB = database.loadAllCases();
 cellfun(@(SCase)SCase.appendToCSV('anatomy.csv'),SCaseDB);     % Call method appendToCSV() on all element of the SCaseDB cell array
 
 fprintf(logFileID, ': OK');
 
 %% Save the entire SCaseDB array as a matlab file
 fprintf(logFileID, '\n\nSave SCase database');
 filename = 'SCaseDB';
 filename = [matlabDir '/' filename '.mat'];
 try
   save(filename, 'SCaseDB');
   fprintf(logFileID, ': OK');
 catch
   fprintf(logFileID, ': Error');
 end
 
 %fprintf(logFileID, [': ' csvFilename]);
 
 %% Write xls file  (Windows only)
 % [~,message] = csv2xlsSCase(csvFilename);
 % fprintf(logFileID, message);
 
 % If run from non-windows system, only csv will be produced. The xls can be
 % produced by opening the csv from Excel and save as xls.
 % Could be a function with 1 input (cvsFilenames
 
 
 
 %{
 xlsFilename = strrep(csvFilename, '.csv', '.xls'); % Replace .csv by .xls
 if ispc % Check if run from windows!
     % Read cvs file as table and write the table as xls
     cvsTable = readtable(csvFilename); % Get table from csv (only way to read non-numeric data)
     cvsCell = table2cell(cvsTable); % Tranform table cell array
     dataHeader = cvsTable.Properties.VariableNames; % Get header
     cvsCell = [dataHeader; cvsCell]; % Add header
 
     sheet = 'anatomy'; % Sheet name of the xls file
 
     [status,message] = xlswrite(xlsFilename, cvsCell, sheet); % Save xls
     if status
         fprintf(logFileId, ['\nSaved ' xlsFilename]);
     else
         fprintf(logFileId, ['\nCould not save ' xlsFilename]);
         fprintf(logFileId, ['\n' message.identifier]);
         fprintf(logFileId, ['\n' message.message]);
     end
 else
     fprintf(logFileId, ['\n' xlsFilename ' not saved. Open and save as xls from Excel']);
 end
 %}
 
 
 %% Close log file
 fprintf(logFileID, '\n\nSCase measured: %s', num2str(length(SCaseList))); % Number of SCases measured
 
 elapsedTime = char(datetime('now')-startTime);
 elapsedTime = [elapsedTime(1:2) ' hours ' elapsedTime(4:5) ' minutes ' elapsedTime(7:8) ' seconds'];
 
 fprintf(logFileID, ['\nTotal Elapsed time: ' elapsedTime ]);
 fprintf(logFileID, '\n');
 fclose(logFileID); % Close log file
 
 %% Output of the SCase if required by argument
 % SCase.output % inputArg = abaqus/density/references/
 
 % update mySQL
 % SCase.sqlSave
 
 
 status  = 1;
 message = 'OK';
 
 end
 
 
 
 function measureShoulderMorphology(shoulder,logFileID)
 
   fprintf(logFileID, '\n|   Load scapula surface and landmarks');
   output = shoulder.scapula.load;
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % scapulaLoad{end+1} = SCaseID;
           return
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Set coord. syst.');
   output = shoulder.scapula.coordSysSet;
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % scapulaCoord{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Load glenoid surface');
   output = shoulder.scapula.glenoid.load;
 
         if ~output
           fprintf(logFileID, '\n|   glenoid surface loading error');
           % glenoidLoad{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Measure glenoid anatomy');
   output = shoulder.scapula.glenoid.morphology;
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % glenoidMorphology{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Load humerus data');
   output = shoulder.humerus.load;
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % humerusLoad{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Measure humerus subluxation');
   output = shoulder.humerus.subluxation(shoulder.scapula);
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % humerusSubluxation{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
   fprintf(logFileID, '\n|   Measure acromion anatomy');
   % Should be run after SCase.shoulderManual.humerus (for AI)
   output = shoulder.scapula.acromion.morphology;
 
         if ~output
           fprintf(logFileID, ': ERROR');
           % acromionMorphology{end+1} = SCaseID;
           return;
         end
 
   fprintf(logFileID, ': OK');
 end
 
 
 
 
 function measureGlenoidDensity(shoulder,logFileID)
   fprintf(logFileID, '\n|   Measure glenoid density');
   output = shoulder.scapula.glenoid.calcDensity;
 
   if ~output
     fprintf(logFileID, ': Density calculus error. Could not read dicom');
     return;
   end
 
   fprintf(logFileID, ': OK');
 end
 
 
 function measureMusclesDegeneration(shoulder,logFileID)
   fprintf(logFileID, '\n|   Measure muscle''s degeneration');
 
   % auto segmentation with most recent slicing methods
   try
-    shoulder.muscles.createRotatorCuffSlicesMatthieu('rotatorCuffMatthieu');
-    shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuffMatthieu','autoMatthieu');
+    shoulder.rotatorCuff.createSliceMatthieu();
+    shoulder.rotatorCuff.segmentMuscles('rotatorCuffMatthieu','autoMatthieu');
   catch ME
     fprintf(logFileID,ME.message);
   end
 
   % auto segmentation with former slicing methods
   try
-    shoulder.muscles.createRotatorCuffSlicesNathan('rotatorCuffNathan');
-    shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuffNathan','autoNathan');
+    shoulder.rotatorCuff.createSliceNathan();
+    shoulder.rotatorCuff.segmentMuscles('rotatorCuffNathan','autoNathan');
   catch ME
     fprintf(logFileID,ME.message);
   end
 
 
   % muscles measurements
   try
-    shoulder.muscles.measureAllMuscles('rotatorCuffMatthieu','autoMatthieu');
+    shoulder.rotatorCuff.measure('rotatorCuffMatthieu','autoMatthieu');
   catch ME
     fprintf(logFileID,ME.message);
   end
   
   fprintf(logFileID, ': OK');
 end