diff --git a/ShoulderCase/@CasesFinder/CasesFinder.m b/ShoulderCase/@CasesFinder/CasesFinder.m
index b6fb4a8..621ee01 100644
--- a/ShoulderCase/@CasesFinder/CasesFinder.m
+++ b/ShoulderCase/@CasesFinder/CasesFinder.m
@@ -1,151 +1,151 @@
 classdef (Abstract) CasesFinder < handle
-  % The CaseFinder class is looking for data related to a ShoulderCase in a folder tree.
+  % The CasesFinder class is looking for data related to a ShoulderCase in a folder tree.
   % It has been designed to work with the datastructure found in lbovenus.epfl.ch/shoulder/.
   %
   % A concrete subclass of CasesFinder must implement the expression related to the
   % 'branch' of the datastructure it is looking at.
   % Then, the CasesFinder will list all the cases that can be found following a 
   % fullfile(dataRootPath,casesSpecificExpression) path.
   % The found cases are then filtered according to the dicom files found following their
   % path.
 
 
   properties (Access = protected)
     dataRootPath
     casesSpecificExpression
     casesDirectories
   end
 
 
 
   methods (Abstract)
     setCasesSpecificExpression(obj);
   end
 
 
 
   methods
 
     function obj = CasesFinder(dataRootPath)
       obj.dataRootPath = dataRootPath;
       obj.casesSpecificExpression = [];
       obj.casesDirectories = containers.Map;
 
       obj.setCasesSpecificExpression;
       obj.mapCasesDirectories;
       obj.removeCasesWithNoCTDirectory;
       obj.reachBestCTDirectory;
     end
 
   end
 
 
 
   methods (Access = protected)
 
     function mapCasesDirectories(obj)
       cases = dir(fullfile(obj.dataRootPath,obj.casesSpecificExpression));
       cases = obj.addSCaseIDToDirOutput(cases);
       obj.updateMapWithCases(cases);
     end
 
     function dirOutput = addSCaseIDToDirOutput(obj,dirOutput)
       for i = 1:length(dirOutput)
         dirOutput(i).SCaseID = regexp(dirOutput(i).name,'[NP]\d{1,3}','match','once');
       end
     end
 
 
     function updateMapWithCases(obj,cases)
       for i = 1:length(cases)
         if not(isempty(cases(i).SCaseID))
           obj.casesDirectories(cases(i).SCaseID) = fullfile(cases(i).folder,cases(i).name);
         end
       end
     end
 
     function removeCasesWithNoCTDirectory(obj)
       SCases = obj.casesDirectories.keys;
       for i = 1:length(SCases)
         SCase = SCases{i};
         if not(obj.containsCTDirectories(obj.casesDirectories(SCase)));
           obj.casesDirectories.remove(SCase);
         end
       end
     end
 
     function output = containsCTDirectories(obj,directory)
       if isempty(dir(fullfile(directory,'CT*')))
         output = false;
       else
         output = true;
       end
     end
 
     function reachBestCTDirectory(obj)
       SCases = obj.casesDirectories.keys;
       for i = 1:length(SCases)
         SCase = SCases{i};
         caseWithBestDirectory = obj.getBestCTDirectoryForSCase(SCase);
         caseWithBestDirectory = obj.addSCaseIDToDirOutput(caseWithBestDirectory);
         obj.updateMapWithCases(caseWithBestDirectory)
       end
     end
 
     function output = getBestCTDirectoryForSCase(obj,SCase)
       %  Each SCase directory should at least contain a
       %  directory with a name staring with 'CT' and ending with either
       %  '-1' : preop/shoulder/sharp/noPhantom
       %  '-2' : preop/shoulder/smooth/noPhantom
       %  '-3' : preop/shoulder/sharp/phantom
       %  '-4' : preop/shoulder/smooth/phantom
       %  and containing a 'dicom' subdirectory. The priority is '1', '2', '3', '4',
       %  unless there is a amira folder already present in one of these diretories.
       observedCTDirectories = dir(fullfile(obj.casesDirectories(SCase),'CT*'));
       for i = 1:length(observedCTDirectories)
         if obj.containsAmiraDirectory(observedCTDirectories(i))
           output = observedCTDirectories(i);
           return
         end
       end
       output = observedCTDirectories(1);
     end
 
     function output = containsAmiraDirectory(obj,directory)
       if isempty(dir(fullfile(directory.folder,directory.name,'amira')))
         output = false;
       else
         output = true;
       end
     end
 
   end
 
 
 
   methods
 
     function output = getPathTo(obj,SCase)
       if obj.containsCase(SCase)
         output = obj.casesDirectories(SCase);
       else
         output = [];
       end
     end
 
     function output = getAllCasesIDs(obj)
       output = obj.casesDirectories.keys;
     end
 
     function output = getNumberOfCases(obj)
       output = obj.casesDirectories.Count;
     end
 
     function output = containsCase(obj,SCase)
       output = obj.casesDirectories.isKey(SCase);
     end
   end
 
 
 
 end
diff --git a/ShoulderCase/@Circle/Circle.m b/ShoulderCase/@Circle/Circle.m
index b89b8c6..69145f7 100644
--- a/ShoulderCase/@Circle/Circle.m
+++ b/ShoulderCase/@Circle/Circle.m
@@ -1,40 +1,45 @@
 classdef Circle < handle
+% Create an array of points on a circle in the 3D space.
+% 
+% Designed to plot planes of shoulder cases.
+
+
 
   properties
     points
     numberOfPoints
     radius
     center
     normal
   end
 
 
 
   methods
     function obj = Circle(center,radius,normal,numberOfPoints)
       obj.center = center;
       obj.radius = radius;
       obj.numberOfPoints = numberOfPoints;
       obj.normal = normal;
 
       obj.createPoints;
     end
 
     function createPoints(obj)
       obj.points = obj.radius*[cos(2*pi*(1:obj.numberOfPoints)/obj.numberOfPoints)',...
                                sin(2*pi*(1:obj.numberOfPoints)/obj.numberOfPoints)',...
                                zeros(obj.numberOfPoints,1)];
        obj.rotate(vrrotvec([0 0 1],obj.normal));
        obj.translate(obj.center);
     end
 
     function rotate(obj,rotationVector)
       obj.points = (vrrotvec2mat(rotationVector)*obj.points')';
     end
 
     function translate(obj,vector)
       obj.points = obj.points + vector;
     end
   end
 
 end
diff --git a/ShoulderCase/@ConfigFileExtractor/ConfigFileExtractor.m b/ShoulderCase/@ConfigFileExtractor/ConfigFileExtractor.m
index f5cf806..66fab0b 100644
--- a/ShoulderCase/@ConfigFileExtractor/ConfigFileExtractor.m
+++ b/ShoulderCase/@ConfigFileExtractor/ConfigFileExtractor.m
@@ -1,18 +1,21 @@
 classdef ConfigFileExtractor < FileExtractor
+% Extract variables from a 'config.txt' file in the current directory.
+
+
 
   properties (Constant)
     filename = 'config.txt';
   end
 
 
   methods (Static)
 
     function output = getVariable(variableName)
       output = FileExtractor.getVariableFromTextFile(variableName,ConfigFileExtractor.filename);
       assert(not(isempty(output)),'%s is not defined in the config file',variableName);
     end
 
   end
 
 
 end
diff --git a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m
index 4b3c16e..042a191 100644
--- a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m
+++ b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m
@@ -1,141 +1,149 @@
 classdef CoordinateSystem < handle
+% Useful to define, use, and test a coordinate system.
+%
+% Can be tested on orthogonality, right or left handedness.
+%
+% Can be used to express global points in local coordinates
+% and vice versa.
+% Can be used to project points on planes spanned by the
+% CoordinateSystem axes.
   
   properties
     origin = [];
   end
 
 
   properties (Hidden = true)
     xAxis = [];
     yAxis = [];
     zAxis = [];
 
     rotationMatrix = [];
   end
 
 
 
   methods
 
     function obj = CoordinateSystem()
     end
 
     function output = isEmpty(obj)
       output = isempty(obj.origin) ||...
                isempty(obj.xAxis) ||...
                isempty(obj.yAxis) ||...
                isempty(obj.zAxis);
     end
 
     function set.xAxis(obj,value)
       obj.xAxis = value;
       obj.rotationMatrix(:,1) = value'/norm(value);
     end
 
     function set.yAxis(obj,value)
       obj.yAxis = value;
       obj.rotationMatrix(:,2) = value'/norm(value);
     end
 
     function set.zAxis(obj,value)
       obj.zAxis = value;
       obj.rotationMatrix(:,3) = value'/norm(value);
     end
 
     function output = getRotationMatrix(obj)
       assert(obj.isOrthogonal,"CoordinateSystem is not orthogonal. Can't get its rotation matrix");
       output = obj.rotationMatrix;
     end
 
     function localPoints = express(obj,globalPoints)
       % Give the coordinates of global cartesian points in present coordinate system
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       assert(size(globalPoints,2) == 3,'Points have to be a Nx3 double array.')
       localPoints = (globalPoints-obj.origin) * obj.rotationMatrix;
     end
 
     function globalPoints = getGlobalCoordinatesOfLocalPoints(obj,localPoints)
         % Give the coordinates of local points in global cartesian coordinate system
         assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
         assert(size(localPoints,2) == 3,'Points have to be a Nx3 double array.')
         globalPoints = (localPoints * obj.rotationMatrix') + obj.origin;
     end
 
     function projected = projectOnXYPlane(obj,points)
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.xAxis);(obj.origin + obj.yAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function projected = projectOnYZPlane(obj,points)
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.yAxis);(obj.origin + obj.zAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function projected = projectOnZXPlane(obj,points)
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.zAxis);(obj.origin + obj.xAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function output = isOrthogonal(obj)
       % The axes comparison can't be done looking for strict equality due to
       % rounded values. Therefore the values must be evaluated with a given
       % tolerance
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of
                                  % a millionth of a radian
       if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
       if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle)
         output = false;
         return
       end
       if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isRightHanded(obj)
       % The comparison between a theoretical right-handed axis and the actual
       % z Axis can't be done looking for strict vector equality due to rounded
       % values. Therefore the norm of the sum of the two vectors is evaluated.
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       rightHandedAxis = cross(obj.xAxis,obj.yAxis);
       if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isLeftHanded(obj)
       assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.');
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       if obj.isRightHanded
         output = false;
         return
       end
       output = true;
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@CoordinateSystemScapula/CoordinateSystemScapula.m b/ShoulderCase/@CoordinateSystemScapula/CoordinateSystemScapula.m
index f60a0e4..1e65cc5 100644
--- a/ShoulderCase/@CoordinateSystemScapula/CoordinateSystemScapula.m
+++ b/ShoulderCase/@CoordinateSystemScapula/CoordinateSystemScapula.m
@@ -1,35 +1,35 @@
 classdef CoordinateSystemScapula < CoordinateSystem
-
+% Add the anatomical axes as properties.
 
 
   properties
     ML = [];
     PA = [];
     IS = [];
   end
 
 
   
   methods
 
     function obj = CoordinateSystemScapula()
     end
 
     function set.PA(obj,value)
       obj.xAxis = value;
       obj.PA = value;
     end
 
     function set.IS(obj,value)
       obj.yAxis = value;
       obj.IS = value;
     end
 
     function set.ML(obj,value)
       obj.zAxis = value;
       obj.ML = value;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@DicomVolume/DicomVolume.m b/ShoulderCase/@DicomVolume/DicomVolume.m
index d93d3a3..200d23a 100644
--- a/ShoulderCase/@DicomVolume/DicomVolume.m
+++ b/ShoulderCase/@DicomVolume/DicomVolume.m
@@ -1,24 +1,28 @@
 classdef (Abstract) DicomVolume < handle
+% Used to manipulate the data given by dicomreadVolume() and dicominfo().
+% 
+% Can be used to retrieve a point coordinates given its indices, and vice
+% versa.
 
   properties
     volume = []
     spatial = []
     dicomInfo = []
     dicomFolderPath = []
   end
 
 
 
   methods
     function set.volume(obj,volume)
       assert(length(size(volume))>=3,'Volume must be a 3D array')
       obj.volume = double(squeeze(volume));
     end
 
     function set.spatial(obj,spatial)
       requiredStructureField = {'PatientPositions','PixelSpacings','PatientOrientations'};
       assert(all(contains(fields(spatial),requiredStructureField)),'Invalid spatial structure');
       obj.spatial = spatial;
     end
   end
 end
diff --git a/ShoulderCase/@DicomVolume/getPointIndexInVolume.m b/ShoulderCase/@DicomVolume/getPointIndexInVolume.m
index 147a384..1ec0573 100644
--- a/ShoulderCase/@DicomVolume/getPointIndexInVolume.m
+++ b/ShoulderCase/@DicomVolume/getPointIndexInVolume.m
@@ -1,9 +1,8 @@
 function output = getPointIndexInVolume(obj,point)
   % the output are the three indices of the given point in the given volume.
   % i,j,k: output indices
   k = find(abs(obj.spatial.PatientPositions(:,3)-point(3)) == min(abs(obj.spatial.PatientPositions(:,3)-point(3))));
   i = round( (point(1) - obj.spatial.PatientPositions(k,1)) / obj.spatial.PixelSpacings(k,1) );
   j = round( (point(2) - obj.spatial.PatientPositions(k,2)) / obj.spatial.PixelSpacings(k,2) );
-  % assert(all([i j k] < size(obj.volume)),'Indices found are outside the volume boundaries');
   output = [i,j,k];
 end
diff --git a/ShoulderCase/@DicomVolume/setVolumeToSubvolume.m b/ShoulderCase/@DicomVolume/setVolumeToSubvolume.m
index 1d92741..9ee5dd9 100644
--- a/ShoulderCase/@DicomVolume/setVolumeToSubvolume.m
+++ b/ShoulderCase/@DicomVolume/setVolumeToSubvolume.m
@@ -1,26 +1,29 @@
 function setVolumeToSubvolume(obj,center,x,y,z)
+  % NOT SAFE TO USE -
+  % NEED FURTHER TEST AND FIX
+  % 
   % Set the volume around the center point with x, y, and z given in milimeters.
-
+  
   % find boundaries
   volumeSize = size(obj.volume);
   minXYZ = obj.getPointIndexInVolume(center - [x/2 y/2 z/2]);
   maxXYZ = obj.getPointIndexInVolume(center + [x/2 y/2 z/2]);
   left = max(1,minXYZ(1));
   right = min(volumeSize(1),maxXYZ(1));
   front = max(1,minXYZ(2));
   rear = min(volumeSize(2),maxXYZ(2));
   bottom = max(1,minXYZ(3));
   top = min(volumeSize(3),maxXYZ(3));
 
   % set subvolume
   obj.volume = obj.volume(left:right,front:rear,bottom:top);
 
   % update spatial
   obj.spatial.PixelSpacings = obj.spatial.PixelSpacings(bottom:top,:);
 
   obj.spatial.PatientPositions = obj.spatial.PatientPositions(bottom:top,:);
   obj.spatial.PatientPositions(:,1) = obj.spatial.PatientPositions(:,1) + left*obj.spatial.PixelSpacings(:,1);
   obj.spatial.PatientPositions(:,2) = obj.spatial.PatientPositions(:,2) + front*obj.spatial.PixelSpacings(:,2);
 
   obj.spatial.PatientOrientations = obj.spatial.PatientOrientations(:,:,bottom:top);
 end
diff --git a/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
index 81a2684..da5f425 100644
--- a/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
+++ b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m
@@ -1,17 +1,34 @@
 classdef DicomVolumeNormaliser < handle & DicomVolume
+% Used to normalise dicom volumes spacings.
+%
+% Dicom volume are not necessarily isotrope, often slices spacing
+% is not equal to the slices' pixel spacing.
+% This class normalise the spacings in the three directions
+% based on the smallest spacing found divided by a chosen
+% factor.
+% The higher the factor, the more time will take the computation
+% of the normalisation.
+%
+% Once volumic data are loaded with DicomeVolume methods, the
+% volume normalisation is achieved in two steps:
+%
+%   divisionFactor = 2;   % 2 has been a good enough tradeoff
+%                         % while testing this class
+%   loadedDicomVolume.setMinimalResolutionDividedBy(divisionFactor);
+%   loadedDicomVolume.normaliseVolume();
 
   properties
     resolution = [];
   end
 
 
 
   methods
     function obj = DicomVolumeNormaliser(varargin)
       if (nargin == 1)
         obj.loadDataFromFolder(varargin{1});
       end
     end
   end
 
 end
diff --git a/ShoulderCase/@DicomVolumeNormaliser/normaliseVolume.m b/ShoulderCase/@DicomVolumeNormaliser/normaliseVolume.m
index 62e2708..c69826c 100644
--- a/ShoulderCase/@DicomVolumeNormaliser/normaliseVolume.m
+++ b/ShoulderCase/@DicomVolumeNormaliser/normaliseVolume.m
@@ -1,14 +1,14 @@
 function normaliseVolume(obj)
   normalisedSizeZ = round(abs(obj.spatial.PatientPositions(1,3) - obj.spatial.PatientPositions(end,3))/obj.resolution);
   normalisedSizeX = round(max(obj.spatial.PixelSpacings(:,1),[],'all')*size(obj.volume,1)/obj.resolution);
   normalisedSizeY = round(max(obj.spatial.PixelSpacings(:,2),[],'all')*size(obj.volume,2)/obj.resolution);
   try
     obj.volume = imresize3(obj.volume,[normalisedSizeX normalisedSizeY normalisedSizeZ]);
   catch
-    % Retry in case of "out of memory" error
+    % Retry with casting volume data type from 'double' to 'single' in case of "out of memory" error.
     obj.volume = imresize3(single(obj.volume),[normalisedSizeX normalisedSizeY normalisedSizeZ]);
   end
   newZPositions = (obj.spatial.PatientPositions(1,3):obj.resolution:obj.spatial.PatientPositions(end,3))';
   obj.spatial.PatientPositions = [repmat(obj.spatial.PatientPositions(1,1:2),length(newZPositions),1) newZPositions];
   obj.spatial.PixelSpacings = repmat(obj.resolution,size(obj.spatial.PatientPositions,1),2);
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
index 46bb45b..7a400bb 100644
--- a/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
+++ b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m
@@ -1,28 +1,34 @@
 classdef DicomVolumeSlicer < handle & DicomVolume & DicomVolumeNormaliser
+% Use a normalised dicom volume to create slices.
+% 
+% The rescale methods have been extracted form the project
+% of Nathan Donini to produce formatted slices data that
+% can be used with the rotator cuff auto segmentation 
+% system and with the MuscleMeasurer class.
 
   properties
     sliced
   end
 
 
 
   properties (Access = public)
     slicedPixelSpacings
     slicedPlane
     slicedX
     slicedY
     slicedZ
     pointsInVolume
   end
 
 
 
   methods
     function obj = DicomVolumeSlicer(varargin)
       if (nargin == 1)
         obj.loadDataFromFolder(varargin{1});
       end
     end
   end
 
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/addEmptyBackgroundToSlice.m b/ShoulderCase/@DicomVolumeSlicer/addEmptyBackgroundToSlice.m
index c4870d5..e3447e0 100644
--- a/ShoulderCase/@DicomVolumeSlicer/addEmptyBackgroundToSlice.m
+++ b/ShoulderCase/@DicomVolumeSlicer/addEmptyBackgroundToSlice.m
@@ -1,32 +1,34 @@
 function addEmptyBackgroundToSlice(obj,height,width)
+  % The background is evenly added around the slice.
+  % height and width must be given in pixels.
+
   sliceHeight = size(obj.sliced,1);
   sliceWidth = size(obj.sliced,2);
   height = max(height,sliceHeight);
   width = max(width,sliceWidth);
 
   top = 1+round(height/2-sliceHeight/2);
   left = 1+round(width/2-sliceWidth/2);
   bottom = top + sliceHeight-1;
   right = left + sliceWidth-1;
 
-  % background = nan(height,width);
   newSlice = repmat(min(obj.sliced,[],'all'),[height width]);
   newSlice(top:bottom,left:right) = obj.sliced;
   obj.sliced = newSlice;
 
   newSliceX = -ones(height,width);
   newSliceX(top:bottom,left:right) = obj.slicedX;
   obj.slicedX = newSliceX;
 
   newSliceY = -ones(height,width);
   newSliceY(top:bottom,left:right) = obj.slicedY;
   obj.slicedY = newSliceY;
 
   newSliceZ = -ones(height,width);
   newSliceZ(top:bottom,left:right) = obj.slicedZ;
   obj.slicedZ = newSliceZ;
 
   newPointsInVolume = false(height,width);
   newPointsInVolume(top:bottom,left:right) = obj.pointsInVolume;
   obj.pointsInVolume = newPointsInVolume;
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/getPointIndexInSlice.m b/ShoulderCase/@DicomVolumeSlicer/getPointIndexInSlice.m
index b98eb2b..d779c45 100644
--- a/ShoulderCase/@DicomVolumeSlicer/getPointIndexInSlice.m
+++ b/ShoulderCase/@DicomVolumeSlicer/getPointIndexInSlice.m
@@ -1,13 +1,25 @@
 function output = getPointIndexInSlice(obj,point)
+  % The current implementation actually returns the point in
+  % the slice that is the closest to the given point, even if
+  % the latter is clearly not in the slice (!).
+  %
+  % The former implementation has been kept but commented out.
+  % It was close to the solution I (Matthieu Boubat) proposed
+  % here https://ch.mathworks.com/matlabcentral/answers/524181-how-can-i-retrieve-x-y-coordinates-of-a-sliced-plane-through-a-3d-volume-using-the-obliqueslice-func.
+  % This solution has been discarded because it failed sometime
+  % (did not find point indices).
+  %
+  % This method might be refactored at some point.
+  
   pointIndices = obj.getPointIndexInVolume(point);
 
   coordinatesError = (abs(obj.slicedX-pointIndices(1)) + abs(obj.slicedY-pointIndices(2)) + abs(obj.slicedZ-pointIndices(3)));
   [i,j] = ind2sub(size(obj.sliced),...
                   find( coordinatesError == min(coordinatesError,[],'all') ));
   % [i,j] = ind2sub(size(obj.sliced),...
   %                 find( (fix(obj.slicedX) == pointIndices(1) | ceil(obj.slicedX) == pointIndices(1)) &...
   %                       (fix(obj.slicedY) == pointIndices(2) | ceil(obj.slicedY) == pointIndices(2)) &...
   %                       (fix(obj.slicedZ) == pointIndices(3) | ceil(obj.slicedZ) == pointIndices(3)) ));
   assert(not(isempty([i,j])),'Point has not been found in slice');
   output = [i(1),j(1)];
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/measureSlicedPixelSpacings.m b/ShoulderCase/@DicomVolumeSlicer/measureSlicedPixelSpacings.m
index ebe7373..1a19a8d 100644
--- a/ShoulderCase/@DicomVolumeSlicer/measureSlicedPixelSpacings.m
+++ b/ShoulderCase/@DicomVolumeSlicer/measureSlicedPixelSpacings.m
@@ -1,90 +1,69 @@
 function measureSlicedPixelSpacings(obj)
+  % Find top, bottom, leftest, and rightest points in the slice for which
+  % coordinates in the volume can be found. Then measure the top to bottom
+  % and the leftest to rightest vectors' norms.
+
+  % Find the row and the column in the slice that contains the most points
+  % for which coordinates in the volume can be found.
   longestRowIndex = find(sum(obj.pointsInVolume,2) == max(sum(obj.pointsInVolume,2)));
   longestColumnIndex = find(sum(obj.pointsInVolume) == max(sum(obj.pointsInVolume)));
 
   longestRowElements = find(obj.pointsInVolume(longestRowIndex(1),:));
   longestColumnElements = find(obj.pointsInVolume(:,longestColumnIndex(1)));
 
 
 
 
-
-
+  % Find the leftest point
   leftSliceIndex = [longestRowIndex(1),longestRowElements(1)];
   leftVolumeIndex = [round(obj.slicedX(leftSliceIndex(1),leftSliceIndex(2))),...
-                     round(obj.slicedY(leftSliceIndex(1),leftSliceIndex(2))),...
-                     round(obj.slicedZ(leftSliceIndex(1),leftSliceIndex(2)))];
+  round(obj.slicedY(leftSliceIndex(1),leftSliceIndex(2))),...
+  round(obj.slicedZ(leftSliceIndex(1),leftSliceIndex(2)))];
   % Volume boundaries can be exceeded at this point but must be asserted.
   leftVolumeIndex = max(leftVolumeIndex,[1 1 1]);
   leftVolumeIndex = min(leftVolumeIndex,size(obj.volume));
   leftPoint = obj.getPointCoordinates(leftVolumeIndex);
-
-
-
-
-
+  
+  
+  
+  % Find the rightest point
   rightSliceIndex = [longestRowIndex(1),longestRowElements(end)];
   rightVolumeIndex = [round(obj.slicedX(rightSliceIndex(1),rightSliceIndex(2))),...
-                      round(obj.slicedY(rightSliceIndex(1),rightSliceIndex(2))),...
-                      round(obj.slicedZ(rightSliceIndex(1),rightSliceIndex(2)))];
+  round(obj.slicedY(rightSliceIndex(1),rightSliceIndex(2))),...
+  round(obj.slicedZ(rightSliceIndex(1),rightSliceIndex(2)))];
   % Volume boundaries can be exceeded at this point but must be asserted.
   rightVolumeIndex = max(rightVolumeIndex,[1 1 1]);
   rightVolumeIndex = min(rightVolumeIndex,size(obj.volume));
   rightPoint = obj.getPointCoordinates(rightVolumeIndex);
-
-
-
-
-
+  
+  
+  
+  % Find the top point
   topSliceIndex = [longestColumnElements(1),longestColumnIndex(1)];
   topVolumeIndex = [round(obj.slicedX(topSliceIndex(1),topSliceIndex(2))),...
-                    round(obj.slicedY(topSliceIndex(1),topSliceIndex(2))),...
-                    round(obj.slicedZ(topSliceIndex(1),topSliceIndex(2)))];
+  round(obj.slicedY(topSliceIndex(1),topSliceIndex(2))),...
+  round(obj.slicedZ(topSliceIndex(1),topSliceIndex(2)))];
   % Volume boundaries can be exceeded at this point but must be asserted.
   topVolumeIndex = max(topVolumeIndex,[1 1 1]);
   topVolumeIndex = min(topVolumeIndex,size(obj.volume));
   topPoint = obj.getPointCoordinates(topVolumeIndex);
-
-
-
-
-
+  
+  
+  
+  % Find the bottom point
   bottomSliceIndex = [longestColumnElements(end),longestColumnIndex(1)];
   bottomVolumeIndex = [round(obj.slicedX(bottomSliceIndex(1),bottomSliceIndex(2))),...
                        round(obj.slicedY(bottomSliceIndex(1),bottomSliceIndex(2))),...
                        round(obj.slicedZ(bottomSliceIndex(1),bottomSliceIndex(2)))];
   % Volume boundaries can be exceeded at this point but must be asserted.
   bottomVolumeIndex = max(bottomVolumeIndex,[1 1 1]);
   bottomVolumeIndex = min(bottomVolumeIndex,size(obj.volume));
   bottomPoint = obj.getPointCoordinates(bottomVolumeIndex);
 
 
-  % window = true(100,100);
-  % pointsInVolumeSubarray = findSubarray(obj.pointsInVolume,window);
-  % assert(not(isempty(pointsInVolumeSubarray)),'No subarray found');
-  % row = pointsInVolumeSubarray(1,1);
-  % column = pointsInVolumeSubarray(1,2);
-  %
-  % originIndex = [round(obj.slicedX(row,column)),...
-  %                round(obj.slicedY(row,column)),...
-  %                round(obj.slicedZ(row,column))];
-  % originPoint = obj.getPointCoordinates(originIndex);
-  %
-  % bottomLeftIndex = [round(obj.slicedX(row+size(window,1)-1,column)),...
-  %                    round(obj.slicedY(row+size(window,1)-1,column)),...
-  %                    round(obj.slicedZ(row+size(window,1)-1,column))];
-  % bottomLeftPoint = obj.getPointCoordinates(bottomLeftIndex);
-  %
-  % topRightIndex = [round(obj.slicedX(row,column+size(window,2)-1)),...
-  %                  round(obj.slicedY(row,column+size(window,2)-1)),...
-  %                  round(obj.slicedZ(row,column+size(window,2)-1))];
-  % topRightPoint = obj.getPointCoordinates(topRightIndex);
-  %
-  % obj.slicedPixelSpacings = [norm(topRightPoint-originPoint)/(size(window,2)-1),...
-  %                            norm(bottomLeftPoint-originPoint)/(size(window,1)-1)];
 
 
 
   obj.slicedPixelSpacings = [norm(rightPoint - leftPoint)/(length(longestRowElements)-1),...
                              norm(bottomPoint - topPoint)/(length(longestColumnElements)-1)];
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/orientSliceUpwardVector.m b/ShoulderCase/@DicomVolumeSlicer/orientSliceUpwardVector.m
index 540e886..1d20146 100644
--- a/ShoulderCase/@DicomVolumeSlicer/orientSliceUpwardVector.m
+++ b/ShoulderCase/@DicomVolumeSlicer/orientSliceUpwardVector.m
@@ -1,11 +1,15 @@
 function orientSliceUpwardVector(obj,point,vector)
+  % Give a vector in the volume and rotate the slice such that
+  % the projection of this vector onto the slice is pointing 
+  % upward.
+  
   pointIndices = obj.getPointIndexInSlice(point);
   vectorIndices = obj.getPointIndexInSlice(obj.slicedPlane.projectOnPlane(point + 10*vector/norm(vector))) - pointIndices;
   rotation = vrrotvec([vectorIndices 0],[-1 0 0]);
   obj.sliced = imrotate(obj.sliced,rotation(3)*180*rotation(4)/pi);
   obj.slicedX = imrotate(obj.slicedX,rotation(3)*180*rotation(4)/pi);
   obj.slicedY = imrotate(obj.slicedY,rotation(3)*180*rotation(4)/pi);
   obj.slicedZ = imrotate(obj.slicedZ,rotation(3)*180*rotation(4)/pi);
   obj.pointsInVolume = imrotate(obj.pointsInVolume,rotation(3)*180*rotation(4)/pi);
   obj.measureSlicedPixelSpacings;
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToInt16.m b/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToInt16.m
index 2f77384..5ab188c 100644
--- a/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToInt16.m
+++ b/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToInt16.m
@@ -1,4 +1,6 @@
 function rescaleSliceToInt16(obj)
+  % Format used to run muscle measurements
+
   obj.sliced(not(obj.pointsInVolume)) = 0;
   obj.sliced = int16(obj.sliced);
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToUint8.m b/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToUint8.m
index b902d8c..70b0024 100644
--- a/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToUint8.m
+++ b/ShoulderCase/@DicomVolumeSlicer/rescaleSliceToUint8.m
@@ -1,6 +1,8 @@
 function rescaleSliceToUint8(obj)
+  % Format used by the rotator cuff segmentation system
+
   HU_interval_png = [-100 160];
   lin_coef = [HU_interval_png(1) 1 ; HU_interval_png(2) 1]^(-1)*[0 255]';
   obj.sliced = lin_coef(1)*obj.sliced+lin_coef(2);
   obj.sliced = uint8(obj.sliced);
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/resize.m b/ShoulderCase/@DicomVolumeSlicer/resize.m
index baf09ea..80c66d5 100644
--- a/ShoulderCase/@DicomVolumeSlicer/resize.m
+++ b/ShoulderCase/@DicomVolumeSlicer/resize.m
@@ -1,8 +1,12 @@
 function resize(obj,rowsAndColumns)
+  % rowAndColumns should be a 1x2 array specifying the number of pixels wanted
+  % in rows an columns. For further information check the documentation
+  % of imresize().
+
   obj.sliced = imresize(obj.sliced,rowsAndColumns);
   obj.slicedX = imresize(obj.slicedX,rowsAndColumns);
   obj.slicedY = imresize(obj.slicedY,rowsAndColumns);
   obj.slicedZ = imresize(obj.slicedZ,rowsAndColumns);
   obj.pointsInVolume = imresize(obj.pointsInVolume,rowsAndColumns);
   obj.measureSlicedPixelSpacings;
 end
diff --git a/ShoulderCase/@DicomVolumeSlicer/slice.m b/ShoulderCase/@DicomVolumeSlicer/slice.m
index ae7a5d5..2ec22fb 100644
--- a/ShoulderCase/@DicomVolumeSlicer/slice.m
+++ b/ShoulderCase/@DicomVolumeSlicer/slice.m
@@ -1,15 +1,24 @@
 function slice(obj,point,normalVector)
+  % Use the obliqueslice() MATLAB method
+
+  % The Plane object is used to project points on the slice
   obj.slicedPlane = Plane();
   obj.slicedPlane.point = point;
   obj.slicedPlane.normal = normalVector;
+
+  % obliqueslice() uses volumic indices and not absolute coordinates
   pointIndices = obj.getPointIndexInVolume(point);
   vectorIndices = obj.getPointIndexInVolume(point+50*normalVector/norm(normalVector)) - pointIndices;
   [obj.sliced,x,y,z] = obliqueslice(obj.volume,pointIndices,vectorIndices,...
                                     'FillValues',nan);
+
+  % The following properties are used to manipulate the slice and 
+  % must be modified the same was the slice is modified.
   obj.pointsInVolume = not(isnan(obj.sliced));
   obj.sliced = (obj.sliced * obj.dicomInfo.RescaleSlope) + obj.dicomInfo.RescaleIntercept;
   obj.slicedX = x;
   obj.slicedY = y;
   obj.slicedZ = z;
+  
   obj.measureSlicedPixelSpacings
 end
diff --git a/ShoulderCase/@FileExtractor/FileExtractor.m b/ShoulderCase/@FileExtractor/FileExtractor.m
index a5689d8..3ba63d2 100644
--- a/ShoulderCase/@FileExtractor/FileExtractor.m
+++ b/ShoulderCase/@FileExtractor/FileExtractor.m
@@ -1,31 +1,46 @@
 classdef (Abstract) FileExtractor < handle
+% Extract the variables defined in a .txt file.
+%
+% Example:
+%   given an example.txt file wich contains only one line:
+%
+%     foo = 42
+%
+%   running:
+%
+%     test = FileExtractor.getVariableFromTextFile('foo','example.txt');
+%
+%   will attribute the value 42 to the variable 'test'. The same result
+%   that would have been given by running:
+% 
+%     test = 42; 
 
   properties (Abstract, Constant)
     filename
   end
 
 
 
   methods (Static)
 
     function output = getVariableFromTextFile(variableName,filename)
       output = [];
       assert(ischar(filename),'filename must be a char array')
       assert(ischar(variableName),'variableName must be a char array')
       assert( (exist(filename,'file') == 2),'File %s doesn''t exist',filename);
       fid = fopen(filename,'r');
       while and(isempty(output), not(feof(fid)))
         line = fgetl(fid);
         variableFound = extractAfter(line,[variableName ' = ']);
         if not(isempty(variableFound))
           eval(['output = ' variableFound ';']);
         end
       end
       assert(not(isempty(output)),'%s has not been found in %s',variableName,filename);
       fclose(fid);
     end
 
   end
 
 
 end
diff --git a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
index c2e5058..aedc83c 100644
--- a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
+++ b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
@@ -1,36 +1,38 @@
 classdef GlenoidAuto < Glenoid
+% To be used with ShoulderAuto data.
+% The load() method requires a specific implementation.
 
   methods
 
     function outputArg = load(obj)
       % LOADAUTO Load genoid surface from auto segmentation
       %   Load the points of the glenoid surface from matlab dir
 
       SCase = obj.scapula.shoulder.SCase;
       SCaseId4C = SCase.id4C;
 
       matlabDir = SCase.dataMatlabPath;
       side = SCase.shoulderAuto.side;
 
       % Import glenoid surface points
       fileName = ['glenoidSurfaceAuto' side '.ply'];
       fileName = [matlabDir '/' fileName];
       if exist(fileName,'file') == 2
         [points,faces]=loadPly(fileName,1);
         obj.surface.points = points;
         obj.surface.meanPoint = mean(points);
         obj.surface.faces = faces;
       else
         outputArg = 0;
         return
         % error(['Could not find file: ' fileName])
       end
 
       outputArg = 1; % Should report on loading result
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@GlenoidManual/GlenoidManual.m b/ShoulderCase/@GlenoidManual/GlenoidManual.m
index 8032847..cc14619 100644
--- a/ShoulderCase/@GlenoidManual/GlenoidManual.m
+++ b/ShoulderCase/@GlenoidManual/GlenoidManual.m
@@ -1,37 +1,39 @@
 classdef GlenoidManual < Glenoid
+% To be used with ShoulderManual data.
+% The load() method requires a specific implementation.
 
   methods
 
     function outputArg = load(obj)
       % LOAD Summary of this method goes here
       %   Load the points of the glenoid surface from the amira
       %   directoty.
 
       SCase = obj.scapula.shoulder.SCase;
       SCaseId4C = SCase.id4C;
 
       amiraDir = SCase.dataAmiraPath;
 
       % Import glenoid surface points
       fileName = ['ExtractedSurface' SCaseId4C '.stl'];
       fileName = [amiraDir '/' fileName];
       if exist(fileName,'file') == 2
         [points,faces,~]=loadStl(fileName,1);
         obj.surface.points = points;
         obj.surface.meanPoint = mean(points);
         obj.surface.faces = faces;
       else
         warning(['Could not find file: ' fileName]);
         outputArg = 0;
         return;
        %  error(['Could not find file: ' fileName])
       end
 
       outputArg = 1; % Should report on loading result
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@LinkedList/LinkedList.m b/ShoulderCase/@LinkedList/LinkedList.m
index 7245b00..e795819 100644
--- a/ShoulderCase/@LinkedList/LinkedList.m
+++ b/ShoulderCase/@LinkedList/LinkedList.m
@@ -1,187 +1,188 @@
 classdef LinkedList < handle
+% Simple linked list implementation
 
   properties (Access = private, Hidden = true)
     ID
     value
     previous
     next
   end
 
 
   methods
 
     function obj = LinkedList
       obj.ID = 'head';
       obj.value = obj;
       obj.previous = [];
       obj.next = [];
     end
 
     function output = append(obj,newValue,varargin)
       if nargin == 3
         newID = varargin{1};
       else
         newID = num2str(obj.length+1);
       end
       newElement = LinkedList;
       obj.setLastElement(newElement);
       newElement.setID(newID);
       newElement.setValue(newValue);
       output = newID;
     end
 
     function remove(obj,IndexOrID)
       element = obj.getElement(IndexOrID);
       obj.removeElement(element);
     end
 
     function output = pop(obj,varargin)
       element = obj.getLastElement;
       if (nargin == 2)
         element = obj.getElement(varargin{1});
       end
       output = element.value;
       obj.removeElement(element);
     end
 
     function output = find(obj,ID)
       element = obj.getFirstElement;
       index = 1;
       foundLocations = [];
       if isequal(ID,element.ID)
         foundLocations = [foundLocations index];
       end
       while not(isempty(element.next))
         element = element.next;
         index = index + 1;
         if isequal(ID,element.ID)
           foundLocations = [foundLocations index];
         end
       end
       output = foundLocations;
     end
 
     function output = getValue(obj,IndexOrID)
       element = obj.getElement(IndexOrID);
       output = element.value;
     end
 
     function output = getID(obj,IndexOrID)
       element = obj.getElement(IndexOrID);
       output = element.ID;
     end
 
     function output = getAllValues(obj)
       element = obj.getFirstElement;
       output = {};
       while not(isempty(element.next))
         element = element.next;
         output{end+1} = element.value;
       end
     end
 
     function output = getAllIDs(obj)
       element = obj.getFirstElement;
       output = {};
       while not(isempty(element.next))
         element = element.next;
         output{end+1} = element.ID;
       end
     end
 
     function output = isEmpty(obj)
       output = (obj.length == 0);
     end
 
     function output = length(obj)
       element = obj.getFirstElement;
       count = 0;
       while not(isempty(element.next))
         element = element.next;
         count = count+1;
       end
       output = count;
     end
 
     function clear(obj)
       element = obj.getFirstElement;
       element.next = [];
     end
 
   end
 
 
 
   methods (Access = private, Hidden = true)
 
     function setID(obj,ID)
       assert(ischar(ID),'ID should be a character array.')
       obj.ID = ID;
     end
 
     function setValue(obj,value)
       obj.value = value;
     end
 
     function setLastElement(obj,element)
       lastElement = obj.getLastElement;
       lastElement.next = element;
       element.previous = lastElement;
     end
 
     function output = getElement(obj,IndexOrID)
       if ischar(IndexOrID)
         output = obj.getElementWithID(IndexOrID);
       elseif isnumeric(IndexOrID)
         output = obj.getElementWithIndex(IndexOrID);
       end
     end
 
     function output = getElementWithID(obj,ID)
       index = obj.find(ID);
       assert(not(isempty(index)),'There is no element with this ID.')
       index = index(1);
       output = obj.getElementWithIndex(index);
     end
 
     function output = getElementWithIndex(obj,index)
       assert((index >= 0),'The given element index should be greater or equal to 0.')
       element = obj.getFirstElement;
       for i = 2:index
         assert(not(isempty(element.next)),'Index exceeds the number of list elements.')
         element = element.next;
       end
       output = element;
     end
 
     function removeElement(obj,element)
       if isequal(obj,element)
         return
       end
       if not(isempty(element.next))
         element.next.previous = element.previous;
       end
       element.previous.next = element.next;
     end
 
     function output = getLastElement(obj)
       element = obj;
       while not(isempty(element.next))
         element = element.next;
       end
       output = element;
     end
 
     function output = getFirstElement(obj)
       element = obj;
       while not(isempty(element.previous))
         element = element.previous;
       end
       output = element;
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@Logger/Logger.m b/ShoulderCase/@Logger/Logger.m
index 412dd7d..5b96dc4 100644
--- a/ShoulderCase/@Logger/Logger.m
+++ b/ShoulderCase/@Logger/Logger.m
@@ -1,164 +1,174 @@
 classdef Logger < handle
+% Class to create and manage log files.
+%
+% Log files are created automatically and include starting and ending timstamps.
+% Methods are available to create section and delimitations in the 
+% handled log file.
+% 
+% The log file name is chosen based on the logFileCategory argument
+% given to the contructor. This argument is also used to delete
+% the oldest log files according to the maxNumberOfLogFiles (default = 5)
+% used with a Logger(logFileCategory,maxNumberOfLogFiles).
 
   properties (Access = private)
     startTime
     logFileCategory
     maxNumberOfLogFiles
     logFid
     prefix
     suffix
     horizontalDelimiter
   end
 
 
   methods
 
     function obj = Logger(logFileCategory,varargin)
       obj.startTime = datetime('now');
       obj.logFileCategory = logFileCategory;
       if nargin == 2
         obj.maxNumberOfLogFiles = varargin{1};
       else
         obj.maxNumberOfLogFiles = 5;
       end
       obj.prefix = '';
       obj.suffix = '';
       obj.horizontalDelimiter = '_______________________________________';
       obj.createNewLog;
     end
 
     function log(obj,textToLog,varargin)
       % to be used like sprintf().
       % write arguments in current log file.
       textToLog = obj.replaceSpecialCharacter(textToLog);
       textToLog = sprintf(textToLog,varargin{:});
       fprintf(obj.logFid,[textToLog obj.suffix]);
     end
 
     function logn(obj,textToLog,varargin)
       obj.suffix = [obj.suffix '\n'];
       try
         obj.log(textToLog,varargin{:});
       catch ME
         warning(ME.message);
       end
       obj.suffix = obj.suffix(1:end-2);
       fprintf(obj.logFid,obj.prefix);
     end
 
     function newSection(obj,title)
       if isempty(obj.prefix)
         obj.newBlock;
       end
       obj.logn('');
       obj.prefix = [obj.prefix '  '];
       obj.logn(title);
     end
 
 
     function newDelimitedSection(obj,title)
       if not(isempty(obj.prefix))
         prefixSave = obj.prefix;
         obj.prefix = obj.prefix(1);
         obj.logn('');
         obj.logn(obj.horizontalDelimiter);
         obj.prefix = prefixSave;
         obj.logn('');
         obj.prefix = [obj.prefix '  '];
         obj.log(title);
       else
         obj.newSection(title);
       end
     end
 
     function closeSection(obj)
       if (length(obj.prefix) > 3)
         obj.prefix = obj.prefix(1:end-2);
         obj.logn('');
       elseif (length(obj.prefix) > 0)
         obj.closeBlock;
       end
     end
 
     function closeBlock(obj)
       if (length(obj.prefix) > 0)
         obj.prefix = obj.prefix(1);
         obj.logn('');
         obj.log(obj.horizontalDelimiter);
         obj.prefix = '';
         obj.logn(obj.horizontalDelimiter);
         obj.logn('');
         obj.logn('');
       end
     end
 
     function closeLog(obj)
       obj.delete;
     end
 
     function deleteExtraLogFiles(obj)
       existingLogFiles = dir(['log\' obj.logFileCategory '_*.log']);
       numberOfExtraFiles = (length(existingLogFiles) - obj.maxNumberOfLogFiles);
       if (numberOfExtraFiles > 0)
         for i = 1:numberOfExtraFiles
           delete(fullfile(existingLogFiles(i).folder,existingLogFiles(i).name));
         end
       end
     end
 
   end
 
 
 
   methods (Access = private, Hidden = true)
 
     function createNewLog(obj)
       if isempty(dir('log'))
         mkdir('log');
       end
       obj.createLogFile;
       obj.deleteExtraLogFiles;
       obj.logHeader;
     end
 
     function createLogFile(obj)
       obj.startTime.Format = '_yyyy_MM_dd_HHmmss';
       obj.logFid = fopen(fullfile('log',[obj.logFileCategory char(obj.startTime) '.log']),'w');
     end
 
     function logHeader(obj)
       obj.logn([obj.logFileCategory '.log']);
       obj.startTime.Format = 'default';
       obj.logn('Date: %s',char(obj.startTime));
       obj.logn('');
     end
 
     function str = replaceSpecialCharacter(obj,str)
       str = strrep(str,'\','/');
     end
 
     function newBlock(obj)
       obj.logn('');
       obj.logn('');
       obj.log(obj.horizontalDelimiter);
       obj.prefix = ['|  '];
       obj.logn(obj.horizontalDelimiter);
       obj.logn('');
     end
 
     function delete(obj)
       obj.closeBlock;
       obj.logn('');
       obj.logn('');
       obj.startTime.Format = 'default';
       obj.logn('Date: %s',char(obj.startTime));
       obj.log('End of log file.');
       fclose(obj.logFid);
       obj.logFid = -1;
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@MethodsMonitor/MethodsMonitor.m b/ShoulderCase/@MethodsMonitor/MethodsMonitor.m
index bcec1ae..9563b07 100644
--- a/ShoulderCase/@MethodsMonitor/MethodsMonitor.m
+++ b/ShoulderCase/@MethodsMonitor/MethodsMonitor.m
@@ -1,116 +1,125 @@
 classdef (Abstract) MethodsMonitor < handle
+% NOT SAFE TO USE -
+% IMPLEMENTATION ONGOING
+%
+% This class purpose is to become a superclass for all other classes we want to monitor.
+% Here, monitoring consist in analysing the subclass' methods and their execution.
+% Monitored data would be the status of the methods ('done','failed','available'...),
+% the error they raised, their running time.
+% It will also be possible to 'order' methods (meaning calling them through MethodsMonitor)
+% and set ordering requirements (e.g. the status of other methods). 
 
   properties (SetObservable)
     methodsMonitor
   end
 
 
   methods
     function startMethodsMonitoring(obj)
       objectMethods = obj.getObjectMethods;
 
       obj.methodsMonitor = table();
       for i = 1:length(objectMethods)
         obj.methodsMonitor = [obj.methodsMonitor, {'';'';''}];
       end
       obj.methodsMonitor.Properties.VariableNames = objectMethods;
       obj.methodsMonitor.Properties.RowNames = {'state';'runtime';'comment'};
 
       addlistener(obj,'methodsMonitor','PostSet',@obj.launchOrderedMethods);
       obj.setAllMethodsUnavailable;
     end
 
     function orderMethod(obj,method)
       obj.setMethodState(method,'ordered')
     end
 
     function setAllMethodsAvailable(obj)
       obj.setAllMethodsState('available');
     end
 
     function setMethodAvailable(obj,method)
       obj.setMethodState(method,'available')
     end
 
     function setAllMethodsUnavailable(obj)
       obj.setAllMethodsState('unavailable');
     end
 
     function setMethodUnavailable(obj,method)
       obj.setMethodState(method,'unavailable')
     end
 
     function setMethodComment(obj,method,comment)
       obj.setMethodVariableValue(method,'comment',comment);
     end
 
     function output = isMethodAvailable(obj,method)
       output = isequal(obj.(method)('state'),{'available'});
     end
   end
 
 
 
   methods (Access = protected)
     function output = getObjectMethods(obj)
       objectMethods = methods(obj);
       handleMethods = methods('handle');
       MethodsMonitorMethods = methods('MethodsMonitor');
       objectMethods = cellArraysDifference(objectMethods,handleMethods,MethodsMonitorMethods);
 
       output = objectMethods;
     end
 
     function launchOrderedMethods(obj,~,~)
       method = obj.getNextOrderedMethod;
       while not(isempty(method))
         obj.runMethod(method);
         method = obj.getNextOrderedMethod;
       end
     end
 
     function setAllMethodsState(obj,state)
       objectMethods = obj.getObjectMethods;
       for i = 1:length(objectMethods)
         obj.setMethodState(objectMethods{i},state);
       end
     end
 
     function setMethodState(obj,method,state)
       obj.setMethodVariableValue(method,'state',state);
     end
 
     function setMethodRuntime(obj,method,runtime)
       obj.setMethodVariableValue(method,'runtime',runtime);
     end
 
     function setMethodVariableValue(obj,method,variable,value)
       obj.methodsMonitor.(method)(variable) = {value};
     end
 
     function runMethod(obj,method)
       Timer.start;
       try
         eval(['obj.' method ';']);
         obj.setMethodState(method,'done');
         obj.setMethodComment(method,'');
       catch ME
         obj.setMethodState(method,'failed');
         obj.setMethodComment(method,ME.message)
       end
       obj.setMethodRuntime(method,Timer.stop);
     end
 
     function output = getNextOrderedMethod(obj)
       output = [];
       objectMethods = obj.methodsMonitor.Properties.VariableNames;
       for i = 1:length(objectMethods)
         if isequal(obj.methodsMonitor.(objectMethods{i})('state'),{'ordered'})
           output = objectMethods{i};
           return;
         end
       end
     end
   end
 
 end
diff --git a/ShoulderCase/@Muscle/Muscle.m b/ShoulderCase/@Muscle/Muscle.m
index b7af6ef..cc8a1b6 100644
--- a/ShoulderCase/@Muscle/Muscle.m
+++ b/ShoulderCase/@Muscle/Muscle.m
@@ -1,71 +1,76 @@
 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 = '';
     
     PCSA = [];
     atrophy = [];
     fat = [];
     osteochondroma = [];
     degeneration = [];
   end
 
 
 
   properties (Hidden = true)
     dataPath = '';
     maskDataPath = '';
     contourDataPath = '';
     container = [];
   end
 
 
   methods
 
     function obj = Muscle(musclesContainer,muscleName,slicesWithMuscleName)
       obj.container = musclesContainer;
       obj.name = muscleName;
       obj.sliceName = slicesWithMuscleName;
     end
 
     function propagateDataPath(obj)
       % Update current dataPath
       % Propagate the path to objects in properties
 
       obj.dataPath = fullfile(obj.container.dataPath,obj.name);
       obj.maskDataPath = fullfile(obj.dataPath,'mask');
       obj.contourDataPath = fullfile(obj.dataPath,'contour');
 
       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
     end
 
 
     function output = getValues(obj)
       output = cellfun(@(property) obj.(property), properties(obj),...
                        'UniformOutput',false);
     end
 
     function measure(obj,segmentationSet)
       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;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
index 48f59d0..7efa616 100644
--- a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
+++ b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
@@ -1,113 +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'])));
     end
 
     function output = getSegmentedImage(obj)
       load(fullfile(obj.muscle.container.slicesDataPath,[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']));
       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 8a3d16e..f2ee3cf 100644
--- a/ShoulderCase/@MusclesContainer/MusclesContainer.m
+++ b/ShoulderCase/@MusclesContainer/MusclesContainer.m
@@ -1,107 +1,121 @@
 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');
 
       if not(isfolder(obj.dataPath))
         mkdir(obj.dataPath);
       end
       if not(isfolder(obj.slicesDataPath))
         mkdir(obj.slicesDataPath);
       end
 
       if obj.isEmpty
         return 
       end
       
       structfun(@(muscle) muscle.propagateDataPath, obj.list);
     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,segmentationSet)
       if obj.isEmpty
         return 
       end
 
       muscleNames = fields(obj.list);
       for i = 1:length(fields(obj.list))
         try
           obj.list.(muscleNames{i}).measure(segmentationSet);
         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/@Plane/Plane.m b/ShoulderCase/@Plane/Plane.m
index b74c8cd..3633ae3 100644
--- a/ShoulderCase/@Plane/Plane.m
+++ b/ShoulderCase/@Plane/Plane.m
@@ -1,81 +1,86 @@
 classdef Plane < handle
+% A plane is decribed by a point and a normal.
+%
+% Plane can define the desired plane's normal and point
+% by fitting an array of points.
+% Plane has a method to project points on the current plane.
 
 
 
   properties
     normal
     point
     fitPerformance
   end
 
 
 
   methods
 
     function obj = Plane()
       obj.normal        = [];
       obj.point         = [];
 
       obj.fitPerformance.points     = [];
       obj.fitPerformance.residuals  = [];
       obj.fitPerformance.RMSE       = [];
       obj.fitPerformance.R2         = [];
     end
 
     function plot(obj,varargin)
       % Plot a circle in the plane. The circle's radius can be given in argument
       radius = 100;
       if nargin == 2
         assert(isnum)
         radius = varargin{1};
       end
       circleInPlane = Circle(obj.point,radius,obj.normal,36);
       XYZ = circleInPlane.points;
       patch(XYZ(:,1),XYZ(:,2),XYZ(:,3),'black','FaceAlpha',0.3); % Transparent plane
     end
 
     function fit(obj,points)
       % Fit current plane to given points
       [coeff,score,roots] = pca(points);
       normal = cross(coeff(:,1),coeff(:,2));
       normal = normal/norm(normal);
 
       meanPoint = mean(points,1);
       estimatedPoints = repmat(meanPoint,size(points,1),1) +...
                                score(:,1:2)*coeff(:,1:2)';
       residuals = points - estimatedPoints;
 
       error =  vecnorm(residuals,2,2);
       RMSE = norm(error)/sqrt(size(points,1));
 
       sumSquareError = sum(error.^2);
       total = vecnorm(points - meanPoint,2,2)';
       sumSquareTotal = sum(total.^2);
       R2 = 1-(sumSquareError/sumSquareTotal); %http://en.wikipedia.org/wiki/Coefficient_of_determination
 
       obj.normal        = normal';
       obj.point         = meanPoint;
       obj.setFitPerformance(points,residuals,RMSE,R2);
     end
 
     function setFitPerformance(obj,points,residuals,RMSE,R2)
       obj.fitPerformance.points = points;
       obj.fitPerformance.residuals = residuals;
       obj.fitPerformance.RMSE = RMSE;
       obj.fitPerformance.R2 = R2;
     end
 
     function setPointOnPlane(obj,point)
       obj.point = point;
     end
 
     function projectedPoints = projectOnPlane(obj,points)
       N2 = obj.normal.'*obj.normal;
       projectedPoints = points*(eye(3)-N2)+repmat(obj.point*N2,size(points,1),1);
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@SCaseIDParser/SCaseIDParser.m b/ShoulderCase/@SCaseIDParser/SCaseIDParser.m
index 7b54d7e..b7e569a 100644
--- a/ShoulderCase/@SCaseIDParser/SCaseIDParser.m
+++ b/ShoulderCase/@SCaseIDParser/SCaseIDParser.m
@@ -1,74 +1,80 @@
 classdef SCaseIDParser < handle
+% Used to validate shoulder case ID.
+% Check the shoulder case type (first character of the ID).
+% Check the shoulder case number (last characters)
+%
+% Can give the shoulder case ID with filled digits.
+% This has been implemented to give ShoulderCase.id4c.
 
   properties (Access = private)
     maxDigits
     SCaseValidTypes
     rawID
   end
 
 
 
   methods
 
     function obj = SCaseIDParser(rawID)
       obj.rawID = rawID;
       assert(ischar(rawID),'Input argument must be a char array.');
 
       obj.maxDigits = ConfigFileExtractor.getVariable('maxSCaseIDDigits');
       assert(not(isempty(obj.maxDigits)),'maxDigits property is empty');
 
       obj.SCaseValidTypes = ConfigFileExtractor.getVariable('SCaseIDValidTypes');
       assert(not(isempty(obj.SCaseValidTypes)),'SCaseValidTypes property is empty');
     end
 
     function output = isValidID(obj)
       types = obj.SCaseValidTypes;
       maxDigits = obj.maxDigits;
       output = obj.SCaseIDHasTypesAndMaxDigits(types,maxDigits);
     end
 
     function output = isNormalCase(obj)
       types = 'N';
       maxDigits = obj.maxDigits;
       output = obj.SCaseIDHasTypesAndMaxDigits(types,maxDigits);
     end
 
     function output = isPathologicalCase(obj)
       types = 'P';
       maxDigits = obj.maxDigits;
       output = obj.SCaseIDHasTypesAndMaxDigits(types,maxDigits);
     end
 
     function output = getID(obj)
       output = obj.rawID;
     end
 
     function output = getIDWithNumberOfDigits(obj,size)
       type = obj.getCaseType;
       number = obj.getCaseNumber;
       fillingZeros = repmat('0',1,size-length(type)-length(number));
       output = [type fillingZeros number];
     end
 
   end
 
 
   methods (Access = private)
 
     function output = SCaseIDHasTypesAndMaxDigits(obj,types,maxDigits)
       expression = ['[' types ']\d{1,' num2str(maxDigits) '}'];
       output = textMatchesExpression(obj.rawID,expression);
     end
 
     function output = getCaseType(obj)
       output = obj.rawID(1);
     end
 
     function output = getCaseNumber(obj)
       output = obj.rawID(2:end);
     end
 
   end
 
 
 end
diff --git a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
index 4779918..58d37b7 100644
--- a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
+++ b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
@@ -1,93 +1,96 @@
 classdef ScapulaAuto < Scapula
+% To be used with ShoulderAuto data.
+% The load() method requires a specific implementation.
+% Instanciate a GlenoidAuto object.
 
   methods
 
     function createGlenoid(obj)
       obj.glenoid = GlenoidAuto(obj);
     end
 
     function outputArg = load(obj)
       % Load 2 files obtained by SSM auto segmentation
       % scapulaSurfaceAuto{L/R}.ply
       % scapulaLandmarksAuto{L/R}.mat
 
       outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise
 
       SCase = obj.shoulder.SCase;
       SCaseId4C = SCase.id4C;
       matlabDir = SCase.dataMatlabPath;
 
       side = obj.shoulder.side; % We might no have it yet !!
 
       if isempty(side)
         side = 'R';
       end
 
       % Try loading auto segmentation from matlab dir
       fileName = ['scapulaSurfaceAuto' side '.ply'];
       fileName = [matlabDir '/' fileName];
 
       if ~exist(fileName,'file')
         side = 'L';
         fileName = ['scapulaSurfaceAuto' side '.ply'];
         fileName = [matlabDir '/' fileName];
       end
 
 
       %side = SCase.shoulder.side;
 
       if exist(fileName,'file') == 2
           try
               face_index_start = 1;
               [points,faces] = loadPly(fileName,face_index_start); % % Load ply file
 
               % ptCloud = pcread(fileName); % load the pointCloud object in the ply file
               % points = ptCloud.Location; % get the array of (x,y,z) points
 
               obj.surface.points = points;
               obj.surface.faces = faces;
               obj.segmentation = 'A'; % Auto segmentation from matlab
           catch
               obj.segmentation = 'E'; % Error on loading
               outputArg = 0;
           end
       else
           obj.segmentation = 'N'; % No segmentation file
           outputArg = 0;
       end
 
       % Try loading auto scapula landmarks from matlab dir
       filename = ['scapulaLandmarksAuto' side '.mat'];
       filename = [matlabDir '/' filename];
       if exist(filename,'file') == 2
           try
               % Load mat file with scapula landmarks
               load(filename, 'ScapulaLandmarks');
               % Set loaded landmarks to scapulaAuto (obj)
               obj.angulusInferior       = ScapulaLandmarks.angulusInferior;
               obj.trigonumSpinae        = ScapulaLandmarks.trigonumSpinae;
               obj.processusCoracoideus  = ScapulaLandmarks.processusCoracoideus;
               obj.acromioClavicular     = ScapulaLandmarks.acromioClavicular;
               obj.angulusAcromialis     = ScapulaLandmarks.angulusAcromialis;
               obj.spinoGlenoidNotch     = ScapulaLandmarks.spinoGlenoidNotch;
               obj.pillar                = ScapulaLandmarks.pillar;
               obj.groove                = ScapulaLandmarks.groove;
           catch
               disp(SCaseId4C); % For debug (2 cases)
               outputArg = 0;
           end
       else
           outputArg = 0;
       end
       try
         obj.measurePlane;
       catch
         outputArg = 0;
       end
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@ScapulaManual/ScapulaManual.m b/ShoulderCase/@ScapulaManual/ScapulaManual.m
index ee4388a..efa2da9 100644
--- a/ShoulderCase/@ScapulaManual/ScapulaManual.m
+++ b/ShoulderCase/@ScapulaManual/ScapulaManual.m
@@ -1,122 +1,125 @@
 classdef ScapulaManual < Scapula
+% To be used with ShoulderManual data.
+% The load() method requires a specific implementation.
+% Instanciate a GlenoidManual object.
 
   methods
 
     function createGlenoid(obj)
       obj.glenoid = GlenoidManual(obj);
     end
 
     function outputArg = load(obj)
       % LOAD Load segmented surface and landmnarks
       %   Load the segmented scapula surface (if exist) and the
       %   scapula landmarks (6 scapula, 5 groove, 5 pillar) from
       %   amira directory.
 
       SCase = obj.shoulder.SCase;
       SCaseId4C = SCase.id4C;
       amiraDir = SCase.dataAmiraPath;
       matlabDir = SCase.dataMatlabPath;
       outputArg = 1; % Should report on loading result
 
       % Scapula (6) landmarks
       fileName = ['ScapulaLandmarks' SCaseId4C '.landmarkAscii'];
       fileName = [amiraDir '/' fileName];
 
       if exist(fileName,'file') == 2
           importedData = importdata(fileName, ' ', 14);
 
       else
           disp(['MATLAB:rmpath:DirNotFound',...
               'Could not find file: ',fileName]);
           outputArg = 0;
           return;
       end
       % Check that imported data are well formed
       try
           landmarks = importedData.data;
       catch
           outputArg = 0;
           return;
       end
 
       obj.angulusInferior      = landmarks(1,:);
       obj.trigonumSpinae       = landmarks(2,:);
       obj.processusCoracoideus = landmarks(3,:);
       obj.acromioClavicular    = landmarks(4,:);
       obj.angulusAcromialis    = landmarks(5,:);
       obj.spinoGlenoidNotch    = landmarks(6,:);
 
       % Scapula suprapinatus groove (5) landmarks
       fileName = ['ScapulaGrooveLandmarks' SCaseId4C '.landmarkAscii'];
       fileName = [amiraDir '/' fileName];
       if exist(fileName,'file') == 2
           importedData = importdata(fileName, ' ', 14);
           % Check that imported data are well formed
           try
               landmarks = importedData.data;
           catch
               outputArg = 0;
               return;
           end
       else
           disp(['MATLAB:rmpath:DirNotFound',...
               'Could not find file: ',fileName]);
           outputArg = 0;
       end
       % Check groove landmarks validity
       if length(landmarks) > 5
           landmarks = landmarks(1:5,:);
           warning(['More than 5 groove landmarks(' SCase.id ')']);
       elseif length(landmarks) < 5
           error('Less than 5 groove landmarks');
       end
 
       obj.groove = landmarks;
 
       % Scapula pillar (5) landmarks
       fileName = ['ScapulaPillarLandmarks' SCaseId4C '.landmarkAscii'];
       fileName = [amiraDir '/' fileName];
       if exist(fileName,'file') == 2
           importedData = importdata(fileName, ' ', 14);
           % Check that imported data are well formed
           try
               landmarks = importedData.data;
           catch
               outputArg = 0;
               return;
           end
       else
           disp(['MATLAB:rmpath:DirNotFound',...
               'Could not find file: ',fileName]);
           outputArg = 0;
       end
       obj.pillar = landmarks;
 
       % Import scapula surface
       % Try manual segmentation in amira folder
       fileName = ['scapula_' SCaseId4C '.stl'];
       fileName = [amiraDir '/' fileName];
       if exist(fileName,'file') == 2
           try
               [points,faces,~]=loadStl(fileName,1);
               obj.surface.points = points;
               obj.surface.faces = faces;
               obj.segmentation = 'M'; % Manual segmentation from amira
           catch
               obj.segmentation = 'E'; % Error in loading
           end
       else
           obj.segmentation = 'N'; % No segmentation
       end
       try
         obj.measurePlane;
       catch
         outputArg = 0;
       end
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m
index 23e9dea..f53b0aa 100644
--- a/ShoulderCase/@Shoulder/Shoulder.m
+++ b/ShoulderCase/@Shoulder/Shoulder.m
@@ -1,226 +1,229 @@
 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
     CTscan
     comment
   end
 
 
 
   properties (Hidden = true)
     dataPath
     SCase
   end
 
 
 
   methods (Abstract, Access = protected)
     createScapula(obj);
   end
 
   methods (Abstract)
     propagateDataPath(obj);
   end
 
   methods
     function obj = Shoulder(SCase)
       obj.side    = '';
       obj.humerus = Humerus(obj);
 
       obj.muscles = MusclesContainer(obj);
 
       obj.CTscan  = '';
       obj.comment = '';
 
       obj.SCase   = SCase;
 
       obj.createScapula;
     end
 
     function output = plot(obj)
       % Create a figure and plot inside:
       % the scapula surface, coordSys, groove points and landmarks;
       % the glenoid surface, sphere and centerLine;
       % the humerus sphere and centerLine.
 
       figure;
       hold on;
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %% Scapula
       % Plot the segmented surface
       if not(isempty(obj.scapula.surface))
           points = obj.scapula.surface.points;
           faces  = obj.scapula.surface.faces;
           x = points(:,1);
           y = points(:,2);
           z = points(:,3);
           trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud');
       end
 
       % Plot scapula groove
       if not(isempty(obj.scapula.groove))
         points = obj.scapula.groove;
         x = points(:,1);
         y = points(:,2);
         z = points(:,3);
         scatter3(x,y,z,100,'blue');
       end
 
       % Plot scapula markers
       points = [...
           obj.scapula.angulusInferior;
           obj.scapula.trigonumSpinae;
           obj.scapula.processusCoracoideus;
           obj.scapula.acromioClavicular;
           obj.scapula.angulusAcromialis;
           obj.scapula.spinoGlenoidNotch...
           ];
       x = points(:,1);
       y = points(:,2);
       z = points(:,3);
       scatter3(x,y,z,100,'blue','LineWidth',2);
 
 
       % Plot scapular plane
       plot(obj.scapula.plane);
 
       % Plot scapular axis
       pt1 = obj.scapula.coordSys.origin;
       axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]);
       pt2 = pt1 - obj.scapula.coordSys.ML * axisLength;
       x = [pt1(1), pt2(1)];
       y = [pt1(2), pt2(2)];
       z = [pt1(3), pt2(3)];
       plot3(x,y,z,'Color','blue','LineWidth', 5);
 
       % Plot scapular coordinate system
       axisLength = 50; % Length of the coord. sys. axes
       pt0 = obj.scapula.coordSys.origin;
       pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X
       pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y
       pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z
       x = [pt0(1), pt1(1)];
       y = [pt0(2), pt1(2)];
       z = [pt0(3), pt1(3)];
       plot3(x,y,z,'Color','red','LineWidth', 5);
       x = [pt0(1), pt2(1)];
       y = [pt0(2), pt2(2)];
       z = [pt0(3), pt2(3)];
       plot3(x,y,z,'Color','green','LineWidth', 5);
       x = [pt0(1), pt3(1)];
       y = [pt0(2), pt3(2)];
       z = [pt0(3), pt3(3)];
       plot3(x,y,z,'Color','blue','LineWidth', 5);
 
       % plot scapular axis on glenoid center
       axisLength = 50;
       pt1 = obj.scapula.glenoid.center;
       pt2 = pt1 + obj.scapula.coordSys.ML*axisLength;
       x = [pt1(1), pt2(1)];
       y = [pt1(2), pt2(2)];
       z = [pt1(3), pt2(3)];
       plot3(x,y,z,'Color','magenta','LineWidth', 5);
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Glenoid
 
       % Plot the glenoid surface
       % We might move them in the direction of th glenoid center
       % instead of the ML axis
 
       points = obj.scapula.glenoid.surface.points +...
           obj.scapula.coordSys.ML * 0.2;
       faces  = obj.scapula.glenoid.surface.faces;
       x = points(:,1);
       y = points(:,2);
       z = points(:,3);
       trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none');
 
 
       % plot glenoid centerline
       pt1 = obj.scapula.glenoid.center;
       pt2 = pt1 + obj.scapula.glenoid.centerLine;
       x = [pt1(1), pt2(1)];
       y = [pt1(2), pt2(2)];
       z = [pt1(3), pt2(3)];
       plot3(x,y,z,'Color','magenta','LineWidth', 5);
 
       % We might add spherical cap
 
       % plot glenoid sphere
       n = 20; % number of lines
       [X,Y,Z] = sphere(n);
       X = X*obj.scapula.glenoid.radius;
       Y = Y*obj.scapula.glenoid.radius;
       Z = Z*obj.scapula.glenoid.radius;
       Xc = obj.scapula.glenoid.center(1);
       Yc = obj.scapula.glenoid.center(2);
       Zc = obj.scapula.glenoid.center(3);
       mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Humerus
 
       % plot humeral head sphere
       n = 20; % number of lines
       [X,Y,Z] = sphere(n);
       X = X*obj.humerus.radius;
       Y = Y*obj.humerus.radius;
       Z = Z*obj.humerus.radius;
       Xc = obj.humerus.center(1);
       Yc = obj.humerus.center(2);
       Zc = obj.humerus.center(3);
       mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
 
 
       % plot humeral head centerline
       pt1 = obj.scapula.glenoid.center;
       pt2 = obj.humerus.center;
       x = [pt1(1), pt2(1)];
       y = [pt1(2), pt2(2)];
       z = [pt1(3), pt2(3)];
       plot3(x,y,z,'Color','magenta','LineWidth', 5);
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Graph properties
 
       axis off; % remove axis
       lightPosition = obj.scapula.glenoid.center + ...
           (...
           obj.scapula.coordSys.ML + ...
           obj.scapula.coordSys.IS + ...
           obj.scapula.coordSys.PA ...
           ) * 100;
       light('Position',lightPosition,'Style','local');
       grid on;
       xlabel('X (CT)');
       ylabel('Y (CT)');
       zlabel('Z (CT)');
       axis square;
       axis vis3d; % fixed limits and scaling
       viewPoint = obj.scapula.coordSys.ML + ...
                   obj.scapula.coordSys.IS + ...
                   obj.scapula.coordSys.PA;
       viewPoint = obj.scapula.coordSys.PA;
       view(viewPoint); % view from lateral side
       % Align IS axis with window vertical axis
       camup(obj.scapula.coordSys.IS);
 
       figureHandle = gca;
       figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis
       figureHandle.Box = 'On';
       set(gcf,'color','w'); % Set background color to white
 
       output = 1;
     end
   end
 end
diff --git a/ShoulderCase/@ShoulderAuto/ShoulderAuto.m b/ShoulderCase/@ShoulderAuto/ShoulderAuto.m
index 974a7bc..e084e21 100644
--- a/ShoulderCase/@ShoulderAuto/ShoulderAuto.m
+++ b/ShoulderCase/@ShoulderAuto/ShoulderAuto.m
@@ -1,31 +1,33 @@
 classdef ShoulderAuto < Shoulder
+% Shoulder using automatically measured data.
+% Instanciate a ScapulaAuto object.
 
   methods (Access = protected)
 
     function createScapula(obj)
       obj.scapula = ScapulaAuto(obj);
     end
 
   end
 
 
 
   methods
 
     function propagateDataPath(obj)
       % Update current dataPath
       % Propagate the path to objects in properties
 
       assert(not(isempty(obj.SCase.dataMatlabPath)),'SCase.dataMatlabPath is empty')
       obj.dataPath = fullfile(obj.SCase.dataMatlabPath,'shoulder_auto');
 
       if not(isfolder(obj.dataPath))
         mkdir(obj.dataPath);
       end
 
       obj.muscles.propagateDataPath;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m b/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m
index fcbdc80..2e9830c 100644
--- a/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m
+++ b/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m
@@ -1,6 +1,7 @@
 function output = getSmoothDicomInfo(obj)
   % Return the result of dicominfo() with the first dicom file found
+  
   dicomFiles = dir(fullfile(obj.getSmoothDicomPath,'*.dcm'));
   assert(not(isempty(dicomFiles)),'No dicom file found there %s',obj.getSmoothDicomPath);
   output = dicominfo(fullfile(dicomFiles(1).folder,dicomFiles(1).name));
 end
diff --git a/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m b/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m
index 862c90d..5ee1de6 100644
--- a/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m
+++ b/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m
@@ -1,6 +1,8 @@
 function output = getSmoothDicomPath(obj)
+  % Return the path to the dicom folder containing 'Smooth' files
+  
   output = [];
   smoothDicomPath = fullfile([obj.dataCTPath(1:end-1) '2'],'dicom');
   assert(isfolder(smoothDicomPath),'No smooth dicom for this SCase. %s is not a valid folder',smoothDicomPath);
   output = smoothDicomPath;
 end
diff --git a/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m b/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m
index 245e1c2..7057cf7 100644
--- a/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m
+++ b/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m
@@ -1,144 +1,167 @@
 classdef ShoulderCaseLoader < handle
+% Class to interact with the database.
+%
+% Use CasesFinder and SCaseIDParser classes to 
+% create a list of valid cases with their location 
+% thanks to a given data root path.
+%
+% Because of cross-systems utilisation of the database
+% the paths stored in the ShoulderCase instances may 
+% vary from one execution to the other. 
+% Thus, this class evolved to be the access point to
+% the ShoulderCase constructor, to avoid misusing
+% the paths. This was also the incentive to create the
+% ShoulderCase.propagateDataPath() method.
+%
+% The main purpose of the ShoulderCaseLoader is to
+% easily load the cases.
+%
+% example:
+%   database = ShoulderCaseLoader() % if not given to the contructor
+%                                   % as an argument, the data root
+%                                   % path will be searched in the
+%                                   % config file.
+%   SCase = database.loadCase('P500');
 
   properties (Access = private)
     normal
     pathological
   end
 
 
 
   methods
 
     function obj = ShoulderCaseLoader(varargin)
       dataRootPath = obj.getDataRootPath(varargin);
       obj.normal = CasesFinderNormal(dataRootPath);
       obj.pathological = CasesFinderPathological(dataRootPath);
     end
 
     function output = containsCase(obj,SCase)
       output = or( obj.normal.containsCase(SCase), obj.pathological.containsCase(SCase));
     end
 
     function output = getNumberOfCases(obj)
       output = obj.normal.getNumberOfCases + obj.pathological.getNumberOfCases;
     end
 
     function output = isEmpty(obj)
       output = not((obj.getNumberOfCases > 0));
     end
 
     function output = getAllNormalCasesIDs(obj)
       output = obj.normal.getAllCasesIDs;
     end
 
     function output = getAllPathologicalCasesIDs(obj)
       output = obj.pathological.getAllCasesIDs;
     end
 
     function output = getAllCasesIDs(obj)
       normalIDs = obj.normal.getAllCasesIDs;
       pathologicalIDs = obj.pathological.getAllCasesIDs;
       output = {};
       for i = 1:obj.normal.getNumberOfCases
         output{end+1} = normalIDs{i};
       end
       for i = 1:obj.pathological.getNumberOfCases
         output{end+1} = pathologicalIDs{i};
       end
     end
 
     function output = loadCase(obj,SCaseID)
       % Load existing ShoulderCase
 
       assert(obj.containsCase(SCaseID),'SCase not found')
 
       % Find ShoulderCase
       container = obj.getContainerOf(SCaseID);
       pathToCase = container.getPathTo(SCaseID);
 
       % Try to find saved ShoulderCase data
       caseMatlabFile = fullfile(pathToCase,'matlab','SCase.mat');
       assert(logical(exist(caseMatlabFile,'file')),'No archive found for this SCase');
 
       loaded = load(caseMatlabFile);
       assert(isfield(loaded,'SCase'),'No SCase field found in the archive');
 
       SCase = loaded.SCase;
       assert(isequal(class(SCase),'ShoulderCase'),'The found SCase is not a ShoulderCase instance');
 
       % Updated output with found archive
       output = SCase;
       % Update data paths
       output.dataCTPath = pathToCase;
       try
         output.propagateDataPath;
       catch ME 
         warning('Impossible to propagate dataPath. The object version might be deprecated.');
         warning(ME.message);
       end
     end
 
     function output = createEmptyCase(obj,SCaseID)
       % Create empty ShoulderCase
 
       assert(obj.containsCase(SCaseID),'Not enough data found to create the SCase. Must contains dicom files.')
 
       % Find ShoulderCase
       container = obj.getContainerOf(SCaseID);
       pathToCase = container.getPathTo(SCaseID);
       SCase = ShoulderCase(SCaseID,pathToCase);
 
       % Updated output with found archive
       output = SCase;
     end
 
     function output = loadAllCases(obj)
       output = {};
       caseList = obj.getAllCasesIDs();
       for i = 1:length(caseList)
         try
           output{end+1} = obj.loadCase(caseList{i});
         catch
           output{end+1} = obj.createEmptyCase(caseList{i});
         end
       end
     end
 
   end
 
 
 
   methods (Access = private, Hidden = true)
 
     function output = getDataRootPath(obj,arguments)
       output = [];
 
       if isempty(arguments)
         TextFileOrPath = 'config.txt';
       elseif not(ischar(arguments{1}))
         return
       else
         TextFileOrPath = arguments{1};
       end
 
       if endsWith(TextFileOrPath,'.txt')
         output = FileExtractor.getVariableFromTextFile('dataDir',TextFileOrPath);
       else
         output = TextFileOrPath;
       end
     end
 
     function output = getContainerOf(obj,SCase)
       output = [];
       if obj.normal.containsCase(SCase)
         output = obj.normal;
       elseif obj.pathological.containsCase(SCase)
         output = obj.pathological;
       end
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/@ShoulderManual/ShoulderManual.m b/ShoulderCase/@ShoulderManual/ShoulderManual.m
index 55bb67e..da707f9 100644
--- a/ShoulderCase/@ShoulderManual/ShoulderManual.m
+++ b/ShoulderCase/@ShoulderManual/ShoulderManual.m
@@ -1,31 +1,33 @@
 classdef ShoulderManual < Shoulder
+% Shoulder using manually measured data.
+% Instanciate a ScapulaManual object.
 
   methods (Access = protected)
 
     function createScapula(obj)
       obj.scapula = ScapulaManual(obj);
     end
 
   end
 
 
 
   methods
 
     function propagateDataPath(obj)
       % Update current dataPath
       % Propagate the path to objects in properties
         
       assert(not(isempty(obj.SCase.dataMatlabPath)),'SCase.dataMatlabPath is empty')
       obj.dataPath = fullfile(obj.SCase.dataMatlabPath,'shoulder_manual');
 
       if not(isfolder(obj.dataPath))
         mkdir(obj.dataPath);
       end
 
       obj.muscles.propagateDataPath;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@Sphere/Sphere.m b/ShoulderCase/@Sphere/Sphere.m
index 3c9873f..23a5975 100644
--- a/ShoulderCase/@Sphere/Sphere.m
+++ b/ShoulderCase/@Sphere/Sphere.m
@@ -1,30 +1,35 @@
 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 (Access = ?Glenoid)
     function obj = Sphere()
       obj.center            = [];
       obj.radius            = [];
       obj.residuals         = [];
       obj.R2                = [];
       obj.RMSE              = [];
     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
   end
 
 end
diff --git a/ShoulderCase/@Timer/Timer.m b/ShoulderCase/@Timer/Timer.m
index 3be5866..485e15c 100644
--- a/ShoulderCase/@Timer/Timer.m
+++ b/ShoulderCase/@Timer/Timer.m
@@ -1,60 +1,92 @@
 classdef Timer
+% Handy timer class to start and stop mutliple timer in parallel.
+%
+% It's basically a linked list with the timers' starting timestamp 
+% stored in the data.
+% 
+% Timer does not need to be instanciated. Simply use its start() and
+% stop() static methods.
+%
+% example:
+%   
+%   Timer.start();
+%   % run some code here
+%   Timer.stop()
+%
+% will return the running time between the two Timer's methods call.
+%
+% example:
+%
+%   timer1 = Timer.start();
+%   % run some code here
+%     timer2 = Timer.start();
+%     % run some other code here
+%     Timer.stop(timer2);
+%   % run some other code here
+%   Timer.stop(timer1)
+%
+% Stop the timers when wanted. If no argument is given to Timer.stop()
+% it will return the running time between Timer.stop() and the
+% last time Timer.start() has been called.
+%
+% The Timer.format property is not used yet but it would be a nice
+% feature.
 
   properties (Constant, Hidden = true)
     data = LinkedList;
     format = struct('string','%sh %sm %ss %sms');
   end
 
 
 
   methods (Static)
 
     function output = start
       output = Timer.data.append(datetime('now'));
     end
 
     function output = stop(varargin)
       if (Timer.data.length == 0)
         error('There is no timer currently running.');
       end
       if (nargin == 1)
         assert(ischar(varargin{1}),'The ID of a timer must be a char array.');
         timerStartTime = Timer.data.pop(varargin{1});
       else
         timerStartTime = Timer.data.pop;
       end
       output = Timer.getElapsedTime(datetime('now')-timerStartTime);
     end
 
     function clearAll
       Timer.data.clear;
     end
 
   end
 
 
   methods (Static, Hidden = true)
 
     function output = getElapsedTime(duration)
       elapsedTime = split(char(duration),':');
       elapsedMilliseconds = Timer.getElapsedMilliseconds(duration);
       output = [elapsedTime{1} ' hours ',...
                 elapsedTime{2} ' minutes ',...
                 elapsedTime{3} ' seconds ',...
                 elapsedMilliseconds ' milliseconds'];
       output = sprintf(Timer.format.string,elapsedTime{1},elapsedTime{2},elapsedTime{3},elapsedMilliseconds);
     end
 
     function output = getElapsedMilliseconds(duration)
       durationInSeconds = num2str(seconds(duration));
       output = extractAfter(durationInSeconds,'.');
       if isempty(output);
         output = '000';
       end
     end
 
   end
 
 
 
 end
diff --git a/ShoulderCase/textMatchesExpression.m b/ShoulderCase/textMatchesExpression.m
index 2f7f46c..94c9931 100644
--- a/ShoulderCase/textMatchesExpression.m
+++ b/ShoulderCase/textMatchesExpression.m
@@ -1,14 +1,15 @@
 function output = textMatchesExpression(text,expression)
+  % Apparently only used in SCaseIDParser. Could be moved there.
   if or( not(ischar(text)), not(ischar(expression)) )
     output = false;
     return;
   end
 
   matchingResult = regexp(text,expression,'match','once');
 
   if isempty(matchingResult)
     output = false;
   else
     output = true;
   end
 end
diff --git a/measureSCase.m b/measureSCase.m
index a1f510f..6831a6f 100755
--- a/measureSCase.m
+++ b/measureSCase.m
@@ -1,670 +1,671 @@
 function [status,message] = measureSCase(varargin)
 % MEASURESCASE Update the entire SCase DB with anatomical measurements.
 
 % The function can be run from (lbovenus) server with the following shell command
 % cd /home/shoulder/methods/matlab/database; matlab -nodisplay -nosplash -r "measureSCase;quit"
 
 % Progress can be checked through log file with following shell command
 % cd /home/shoulder/methods/matlab/database;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.getAllCasesIDs();
 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(contains(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;
     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'})))) )
       orderMorphologyMeasurementManual = true;
       orderDegenerationMeasurementManual = true;
       SCase = loadedSCase;
     end
     if and(calcMusclesDegeneration, not(all(table2array(musclesAuto.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');
 
   try
     shoulder.muscles.createRotatorCuffSlices();
+    shoulder.muscles.createRotatorCuffSlicesOld();
     shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuff');
     shoulder.muscles.measureAllMuscles('auto');
   catch ME
     fprintf(logFileID,ME.message);
   end
 
   fprintf(logFileID, ': OK');
 end