diff --git a/ShoulderCase/@Muscle/Muscle.m b/ShoulderCase/@Muscle/Muscle.m
index 4132f85..2585515 100644
--- a/ShoulderCase/@Muscle/Muscle.m
+++ b/ShoulderCase/@Muscle/Muscle.m
@@ -1,87 +1,103 @@
 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
+    container = nan;
     name = "";
-    segmentationSet = "";
+    segmentationName = "";
     sliceName = "";
     
     PCSA = nan;
     atrophy = nan;
     fat = nan;
     osteochondroma = nan;
     degeneration = nan;
 
-    container = nan;
-    insertions = nan;
+    centroid = nan;
+    forceApplicationPoint = nan;
     forceVector = nan;
   end
 
   methods
     function obj = Muscle(musclesContainer,muscleName)
       obj.container = musclesContainer;
       obj.name = muscleName;
 
       [~, ~] = mkdir(obj.dataPath);
       [~, ~] = mkdir(obj.dataMaskPath);
-      [~, ~] = mkdir(obj.dataContourPath);
     end
 
     function output = dataPath(obj)
       output = fullfile(obj.container.dataPath, obj.name);
     end
 
     function output = dataMaskPath(obj)
       output = fullfile(obj.dataPath, "mask");
     end
 
-    function output = dataContourPath(obj)
-      output = fullfile(obj.dataPath, "contour");
+    function setSliceName(obj, value)
+      obj.sliceName = value;
     end
 
-    function output = getValues(obj)
-      output = cellfun(@(property) obj.(property), properties(obj),...
-        'UniformOutput',false);
+    function setSegmentationName(obj, value)
+      obj.segmentationName = value;
     end
 
-    function output = summary(obj)
-      summary.name = string(obj.name);
-      summary.segmentationSet = string(obj.segmentationSet);
-      summary.sliceName = string(obj.sliceName);
-      summary.PCSA = obj.PCSA;
-      summary.atrophy = obj.atrophy;
-      summary.fat = obj.fat;
-      summary.osteochondroma = obj.osteochondroma;
-      summary.degeneration = obj.degeneration;
-      output = struct2table(summary);
-    end
-
-    function measureDegeneration(obj,sliceName, segmentationSet)
-      obj.sliceName = sliceName;
-      obj.segmentationSet = segmentationSet;
-
-      try
-        measurer = MuscleMeasurer(obj);
-        obj.PCSA = measurer.getPCSA;
-        obj.atrophy = measurer.getRatioAtrophy;
-        obj.fat = measurer.getRatioFat;
-        obj.osteochondroma = measurer.getRatioOsteochondroma;
-        obj.degeneration = measurer.getRatioDegeneration;
-      end
+    function output = loadMask(obj, maskSuffix)
+      output = logical(imread(fullfile(obj.dataMaskPath,...
+        obj.segmentationName+"_"+maskSuffix+".png")));
     end
 
-    function output = loadMask(obj)
-      output = imread(fullfile(obj.dataMaskPath,obj.segmentationSet + "_Mask.png")) > 0;
+    function saveMask(obj, mask, maskSuffix)
+      imwrite(mask, fullfile(obj.dataMaskPath,...
+        obj.segmentationName+"_"+maskSuffix+".png"));
     end
 
     function output = loadPixelCoordinates(obj)
       output = load(fullfile(obj.container.dataSlicesPath,...
         obj.sliceName + "_PixelCoordinates.mat")).imagesPixelCoordinates;
     end
+
+
+
+    function output = measureFirst(obj)
+      % Call methods that can be run after morphology() methods has been run 
+      % by all ShoulderCase objects.
+      obj.setSliceName("rotatorCuffMatthieu");
+      obj.setSegmentationName("autoMatthieu");
+      success = Logger.timeLogExecution(...
+        obj.name+" measure and save contour: ",...
+        @(obj) obj.measureAndSaveContour(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" measure and save centroid: ",...
+        @(obj) obj.measureAndSaveCentroid(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" PCSA: ",...
+        @(obj) obj.measurePCSA(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" atrophy: ",...
+        @(obj) obj.measureAtrophy(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" fat infiltration: ",...
+        @(obj) obj.measureFatInfiltration(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" osteochondroma: ",...
+        @(obj) obj.measureOsteochondroma(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" degeneration: ",...
+        @(obj) obj.measureDegeneration(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" humeral head contact point: ",...
+        @(obj) obj.measureHumerusContactPoint(), obj);
+      success = success | Logger.timeLogExecution(...
+        obj.name+" force: ",...
+        @(obj) obj.measureForceVector(), obj);
+      output = success;
+    end
   end
 
 end
diff --git a/ShoulderCase/@Muscle/getCentroid.m b/ShoulderCase/@Muscle/getCentroid.m
deleted file mode 100644
index c1291ab..0000000
--- a/ShoulderCase/@Muscle/getCentroid.m
+++ /dev/null
@@ -1,25 +0,0 @@
-function output = getCentroid(obj)
-  muscleMask = obj.loadMask();
-  slicePixelCoordinates = obj.loadPixelCoordinates();
-
-  output.mask = getCentroidMask(muscleMask);
-  output.indices = getCentroidIndices(muscleMask);
-  output.coordinates = getCentroidCoordinates(muscleMask,slicePixelCoordinates);
-end
-
-function output = getCentroidMask(muscleMask)
-  centroidIndices = getCentroidIndices(muscleMask);
-  centroidMask = false(size(muscleMask));
-  centroidMask(centroidIndices(1),centroidIndices(2)) = true;
-  output = centroidMask;
-end
-
-function output = getCentroidCoordinates(muscleMask,pixelCoordinates)
-  centroidIndices = getCentroidIndices(muscleMask);
-  centroidCoordinates = pixelCoordinates{centroidIndices(1),centroidIndices(2)} ;
-  output = centroidCoordinates;
-end
-
-function output = getCentroidIndices(muscleMask)
-  output = round(regionprops(muscleMask,"Centroid").Centroid);
-end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/getContour.m b/ShoulderCase/@Muscle/getContour.m
deleted file mode 100644
index e4f03f6..0000000
--- a/ShoulderCase/@Muscle/getContour.m
+++ /dev/null
@@ -1,29 +0,0 @@
-function output = getContour(obj)
-  muscleMask = obj.loadMask();
-  slicePixelCoordinates = obj.loadPixelCoordinates();
-
-  output.mask = getContourMask(muscleMask);
-  output.indices = getContourIndices(muscleMask);
-  output.coordinates = getContourCoordinates(muscleMask,slicePixelCoordinates);
-end
-
-function output = getContourMask(muscleMask)
-  contourMask = bwmorph(muscleMask,"remove");
-  output = contourMask;
-end
-
-function output = getContourCoordinates(muscleMask,pixelCoordinates)
-  contourIndices = getContourIndices(muscleMask);
-  contourCoordinates = [];
-  for row = 1:size(contourIndices,1)
-    contourCoordinates = [contourCoordinates; pixelCoordinates{contourIndices(row,1),contourIndices(row,2)}];
-  end
-  output = contourCoordinates;
-end
-
-function output = getContourIndices(muscleMask)
-  contourMask = getContourMask(muscleMask);
-  [row, column] = ind2sub(size(contourMask),find(contourMask));
-  contourIndices = [column row];
-  output = contourIndices;
-end
diff --git a/ShoulderCase/@Muscle/getMaskCoordinates.m b/ShoulderCase/@Muscle/getMaskCoordinates.m
new file mode 100644
index 0000000..ddac808
--- /dev/null
+++ b/ShoulderCase/@Muscle/getMaskCoordinates.m
@@ -0,0 +1,10 @@
+function output = getMaskCoordinates(obj, maskSuffix)
+  mask = obj.loadMask(maskSuffix);
+  maskIndices = find(mask);
+  slicePixelCoordinates = obj.loadPixelCoordinates();
+  maskCoordinates = [...
+    slicePixelCoordinates.x(maskIndices),...
+    slicePixelCoordinates.y(maskIndices),...
+    slicePixelCoordinates.z(maskIndices)];
+  output = maskCoordinates;
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureAndSaveCentroid.m b/ShoulderCase/@Muscle/measureAndSaveCentroid.m
new file mode 100644
index 0000000..d04222c
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureAndSaveCentroid.m
@@ -0,0 +1,11 @@
+function measureAndSaveCentroid(obj)
+  muscleContour = obj.loadMask("Contour");
+  muscleCentroid = false(size(muscleContour));
+  centroidIndices = round(regionprops(muscleContour).Centroid);
+  % centroidIndices(2) is the vertical top->down index
+  % centroidIndices(1) is the horizontal left->right index
+  % kind of opposite to the matlab array indexing...
+  muscleCentroid(centroidIndices(2), centroidIndices(1)) = true;
+  obj.saveMask(muscleCentroid, "Centroid");
+  obj.centroid = obj.getMaskCoordinates("Centroid");
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureAndSaveContour.m b/ShoulderCase/@Muscle/measureAndSaveContour.m
new file mode 100644
index 0000000..e645430
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureAndSaveContour.m
@@ -0,0 +1,5 @@
+function measureAndSaveContour(obj)
+  muscleSegmentation = obj.loadMask("Segmentation");
+  muscleContour = bwmorph(muscleSegmentation, "remove");
+  obj.saveMask(muscleContour, "Contour");
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureAtrophy.m b/ShoulderCase/@Muscle/measureAtrophy.m
new file mode 100644
index 0000000..1a0afa7
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureAtrophy.m
@@ -0,0 +1,4 @@
+function measureAtrophy(obj)
+  measurer = MuscleMeasurer(obj);
+  obj.atrophy = measurer.getRatioAtrophy();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureDegeneration.m b/ShoulderCase/@Muscle/measureDegeneration.m
new file mode 100644
index 0000000..7db9589
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureDegeneration.m
@@ -0,0 +1,4 @@
+function measureDegeneration(obj)
+  measurer = MuscleMeasurer(obj);
+  obj.degeneration = measurer.getRatioDegeneration();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureFatInfiltration.m b/ShoulderCase/@Muscle/measureFatInfiltration.m
new file mode 100644
index 0000000..a5467e5
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureFatInfiltration.m
@@ -0,0 +1,4 @@
+function measureFatInfiltration(obj)
+  measurer = MuscleMeasurer(obj);
+  obj.fat = measurer.getRatioFat();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureForceVector.m b/ShoulderCase/@Muscle/measureForceVector.m
index bf00aeb..790273e 100644
--- a/ShoulderCase/@Muscle/measureForceVector.m
+++ b/ShoulderCase/@Muscle/measureForceVector.m
@@ -1,5 +1,5 @@
 function measureForceVector(obj)
-  force = Vector(obj.insertions, obj.getCentroid().coordinates);
-  force = force*obj.PCSA / force.norm;
+  force = Vector(obj.forceApplicationPoint, obj.centroid);
+  force = (force / force.norm) * obj.PCSA;
   obj.forceVector = force;
 end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureHumerusContactPoint.m b/ShoulderCase/@Muscle/measureHumerusContactPoint.m
new file mode 100644
index 0000000..703b0c9
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureHumerusContactPoint.m
@@ -0,0 +1,15 @@
+function measureHumerusContactPoint(obj)
+  humerus = obj.container.shoulder.humerus;
+  humeralHeadContactPointFinder = SphereContactPoint(...
+    Sphere(humerus.center, humerus.radius));
+
+  % muscles humeral head insertion is estimated to be the most lateral point of
+  % the humeral head fitted sphere 
+  MLAxis = obj.container.shoulder.scapula.coordSys.ML;
+  muscleInsertion = humerus.center + humerus.radius * MLAxis;
+  humeralHeadContactPointFinder.setAnchorPoint(muscleInsertion);
+
+  humeralHeadContactPointFinder.setPolarPoint(obj.centroid);
+
+  obj.forceApplicationPoint = humeralHeadContactPointFinder.getContactPoint();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measureOsteochondroma.m b/ShoulderCase/@Muscle/measureOsteochondroma.m
new file mode 100644
index 0000000..e05e4a3
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureOsteochondroma.m
@@ -0,0 +1,4 @@
+function measureOsteochondroma(obj)
+  measurer = MuscleMeasurer(obj);
+  obj.osteochondroma = measurer.getRatioOsteochondroma();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Muscle/measurePCSA.m b/ShoulderCase/@Muscle/measurePCSA.m
new file mode 100644
index 0000000..8df5381
--- /dev/null
+++ b/ShoulderCase/@Muscle/measurePCSA.m
@@ -0,0 +1,4 @@
+function measurePCSA(obj)
+  measurer = MuscleMeasurer(obj);
+  obj.PCSA = measurer.getPCSA();
+end
\ No newline at end of file
diff --git a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
index 762cab7..b4d1f54 100644
--- a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
+++ b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m
@@ -1,119 +1,119 @@
 classdef MuscleMeasurer < handle
 % Used to perfom PCSA and degeneration measurements of a Muscle object.
 % Slices and mask images are expected to be found at some specific places
 % and have specific names.
 %
 % This class is the results of extracting measurements methods from the
 % project of Nathan Donini.  
 
   properties
     muscle
     segmentedArea
     segmentedImage
 
     rangeHU = double([-1000 1000]);
     muscleThreshold = 0;
     fatThreshold = 30;
     osteochondromaThreshold = 166;
   end
 
 
 
   methods
 
     function obj = MuscleMeasurer(muscle)
       obj.muscle = muscle;
       obj.segmentedArea = obj.getSegmentedArea;
       obj.segmentedImage = obj.getSegmentedImage;
 
       obj.normalizeThresholds;
     end
 
     function output = getSegmentedArea(obj)
-      output = logical(imread(fullfile(obj.muscle.dataMaskPath,obj.muscle.segmentationSet + "_Mask.png")));
+      output = logical(imread(fullfile(obj.muscle.dataMaskPath,obj.muscle.segmentationName + "_Segmentation.png")));
     end
 
     function output = getSegmentedImage(obj)
       load(fullfile(obj.muscle.container.dataSlicesPath,obj.muscle.sliceName + "_ForMeasurements.mat"));
       normalizedImage = obj.getNormalizedImage(imageForMeasurements);
 
       segmentedImage = normalizedImage.*obj.segmentedArea;
 
       output = segmentedImage;
     end
 
     function output = getNormalizedImage(obj,image)
       image(image < obj.rangeHU(1)) = obj.rangeHU(1);
       image(image > obj.rangeHU(2)) = obj.rangeHU(2);
       output = (double(image)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
     end
 
     function normalizeThresholds(obj)
       obj.muscleThreshold = ((obj.muscleThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
       obj.fatThreshold = ((obj.fatThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
       obj.osteochondromaThreshold = ((obj.osteochondromaThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1));
     end
 
     function output = getPCSA(obj)
       load(fullfile(obj.muscle.container.dataSlicesPath,obj.muscle.sliceName + "_PixelSpacings.mat"));
       pixelSurface = (imagesPixelSpacings(1) * imagesPixelSpacings(2)) / 100; % cm^2
       output = pixelSurface * bwarea(obj.segmentedArea);
     end
 
     function output = getRatioAtrophy(obj)
       output = obj.getRatioAreas(obj.getAreaAtrophy,obj.segmentedArea);
     end
 
     function output = getRatioFat(obj)
       output = obj.getRatioAreas(obj.getAreaFat,obj.segmentedArea);
     end
 
     function output = getRatioOsteochondroma(obj)
       output = obj.getRatioAreas(obj.getAreaOsteochondroma,obj.segmentedArea);
     end
 
     function output = getRatioDegeneration(obj)
       output = obj.getRatioAtrophy + obj.getRatioFat + obj.getRatioOsteochondroma;
     end
 
     function output = getRatioAreas(obj,partialArea,totalArea)
       output = bwarea(partialArea)/bwarea(totalArea);
     end
 
     function output = getAreaMuscle(obj)
       % The area "Muscle" is the area of the segmented image that is not atrophied.
       % Looking for a better name.
       areaMuscle = imbinarize(obj.segmentedImage,obj.muscleThreshold);
       areaMuscle = imfill(areaMuscle,'holes');
       % Keep biggest island only
       islands = bwconncomp(areaMuscle);
       islandsSizes = cellfun(@numel,islands.PixelIdxList);
       [~,biggestIslandIndex] = max(islandsSizes);
       areaMuscle = false(size(areaMuscle));
       areaMuscle(islands.PixelIdxList{biggestIslandIndex}) = true;
 
       output = areaMuscle;
     end
 
     function output = getAreaAtrophy(obj)
       output = obj.segmentedArea & not(obj.getAreaMuscle);
     end
 
     function output = getAreaFat(obj)
       muscleImage = obj.segmentedImage.*obj.getAreaMuscle;
       areaFat = obj.getAreaMuscle & not(imbinarize(muscleImage,obj.fatThreshold));
 
       output = areaFat;
     end
 
     function output = getAreaOsteochondroma(obj)
       muscleImage = obj.segmentedImage.*obj.getAreaMuscle;
       areaOsteochondroma = imbinarize(muscleImage,obj.osteochondromaThreshold);
       areaOsteochondroma = imfill(areaOsteochondroma,'holes');
 
       output = areaOsteochondroma;
     end
 
   end
 
 end
diff --git a/ShoulderCase/@RotatorCuff/RotatorCuff.m b/ShoulderCase/@RotatorCuff/RotatorCuff.m
index 7a65b4a..bfa4014 100644
--- a/ShoulderCase/@RotatorCuff/RotatorCuff.m
+++ b/ShoulderCase/@RotatorCuff/RotatorCuff.m
@@ -1,110 +1,58 @@
 classdef RotatorCuff < handle
 
   properties
     SC
     SS
     IS
     TM
     anteroPosteriorImbalance = nan;
 
     shoulder;
   end
 
   methods
     function obj = RotatorCuff(shoulder)
       obj.shoulder = shoulder;
 
       obj.SC = Muscle(obj, "SC");
       obj.SS = Muscle(obj, "SS");
       obj.IS = Muscle(obj, "IS");
       obj.TM = Muscle(obj, "TM");
 
       [~, ~] = mkdir(obj.dataPath);
       [~, ~] = mkdir(obj.dataSlicesPath);
     end
 
     function output = dataPath(obj)
       output = fullfile(obj.shoulder.dataPath, "muscles");
     end
 
     function output = dataSlicesPath(obj)
       output = fullfile(obj.dataPath, "oblique_slices");
     end
 
     function output = summary(obj)
-      summary = [];
-      summary = [summary; obj.SC.summary];
-      summary = [summary; obj.SS.summary];
-      summary = [summary; obj.IS.summary];
-      summary = [summary; obj.TM.summary];
+      summary = [...
+        getTabulatedProperties(obj.SC);...
+        getTabulatedProperties(obj.SS);...
+        getTabulatedProperties(obj.IS);...
+        getTabulatedProperties(obj.TM)...
+        ];
       output = summary;
     end
 
     function sliceAndSegment(obj, varargin)
       obj.createSliceMatthieu();
       obj.segmentMuscles("rotatorCuffMatthieu", "autoMatthieu");
     end
-
-    function output = measure(obj, sliceName, segmentationSet)
-      try
-        obj.measureDegenerations(sliceName, segmentationSet);
-        obj.measureInsertions();
-        obj.measureAnteroPosteriorImbalance();
-      end
-      output = obj.summary;
-    end
-
-    function measureDegenerations(obj, sliceName, segmentationSet)
-      for muscleName = ["SC", "SS", "IS", "TM"]
-        obj.(muscleName).measureDegeneration(sliceName, segmentationSet);
-      end
-    end
-
-    function measureInsertions(obj)
-      humerus = obj.shoulder.humerus;
-      humeralHead = Sphere(humerus.center, humerus.radius);
-
-      % IS humeral insertion estimation
-      obj.IS.measureInsertionTangentToSphere(...
-        humeralHead,...
-        obj.shoulder.scapula.coordSys.IS,...
-        -obj.shoulder.scapula.coordSys.PA);
-      obj.IS.measureForceVector();
-
-      % TM humeral insertion estimation
-      obj.TM.measureInsertionTangentToSphere(...
-        humeralHead,...
-        obj.shoulder.scapula.coordSys.IS,...
-        -obj.shoulder.scapula.coordSys.PA);
-      obj.TM.measureForceVector();
-
-      % SC humeral insertion estimation
-      obj.SC.measureInsertionTangentToSphere(...
-        humeralHead,...
-        obj.shoulder.scapula.coordSys.IS,...
-        obj.shoulder.scapula.coordSys.PA);
-      obj.SC.measureForceVector();
-    end
-
-    function measureAnteroPosteriorImbalance(obj)
-      forceSC = obj.SC.forceVector;
-      forceIS = obj.IS.forceVector;
-      forceTM = obj.TM.forceVector;
-      scapulaCS = obj.shoulder.scapula.coordSys;
-
-      forceResultant = Vector(scapulaCS.origin, scapulaCS.origin) + ...
-        forceSC + forceIS + forceTM; 
-
-      % Take infero-superior orthogonal component <=> project on scapular axial plane
-      inferoSuperior = Vector(scapulaCS.origin, scapulaCS.origin + scapulaCS.IS);
-      imbalanceVector = inferoSuperior.orthogonalComplementToPoint(forceResultant.target);
-
-      lateroMedial = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.ML);
-      anteroPosterior = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.PA);
-      obj.anteroPosteriorImbalance = ...
-        sign(dot(imbalanceVector, anteroPosterior)) *...
-        rad2deg(lateroMedial.angle(imbalanceVector));
-    end
+    
+    function output = measureSecond(obj)
+      % Call methods that can be run after measureFirst() methods has been run 
+      % by all ShoulderCase objects.
+      success = Logger.timeLogExecution(...
+        "Rotator cuff antero-posterior imbalance: ",...
+        @(obj) obj.measureAnteroPosteriorImbalance(), obj);
+    end    
   end
 
 end
\ No newline at end of file
diff --git a/ShoulderCase/@RotatorCuff/measureAnteroPosteriorImbalance.m b/ShoulderCase/@RotatorCuff/measureAnteroPosteriorImbalance.m
new file mode 100644
index 0000000..e2875ab
--- /dev/null
+++ b/ShoulderCase/@RotatorCuff/measureAnteroPosteriorImbalance.m
@@ -0,0 +1,21 @@
+function measureAnteroPosteriorImbalance(obj)
+  scapulaCS = obj.shoulder.scapula.coordSys;
+
+  forceResultant = Vector(scapulaCS.origin, scapulaCS.origin)...
+    + obj.SC.forceVector...
+    + obj.SS.forceVector...
+    + obj.IS.forceVector...
+    + obj.TM.forceVector;
+
+  % Take infero-superior orthogonal component <=> project on scapular axial plane
+  inferoSuperior = Vector(scapulaCS.origin, scapulaCS.origin + scapulaCS.IS);
+  imbalanceVector = inferoSuperior.orthogonalComplementTo(forceResultant);
+
+  lateroMedial = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.ML);
+  anteroPosterior = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.PA);
+
+  % posterior imbalance is positive 
+  obj.anteroPosteriorImbalance = ...
+    sign(dot(imbalanceVector, anteroPosterior)) *...
+    rad2deg(lateroMedial.angle(imbalanceVector));
+end
\ No newline at end of file
diff --git a/ShoulderCase/@RotatorCuff/segmentMuscles.m b/ShoulderCase/@RotatorCuff/segmentMuscles.m
index f6343dc..d0d3a98 100644
--- a/ShoulderCase/@RotatorCuff/segmentMuscles.m
+++ b/ShoulderCase/@RotatorCuff/segmentMuscles.m
@@ -1,66 +1,66 @@
 function segmentMuscles(obj, sliceName, maskName)
   % sliceName is the part before the '_ForSegmentation.png' part of the name of
   % the file that is sent to rcseg.
-  % maskName is the part before the '_Mask.png' part of the name of
+  % maskName is the part before the '_Segmentation.png' part of the name of
   % the file that is saved in the SCase folder.
   cleanRotatorCuffSegmentationWorkspace(obj);
   sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName);
   callRotatorCuffSegmentation(obj);
   saveSegmentationResults(obj,maskName);
   cleanRotatorCuffSegmentationWorkspace(obj);
 end
 
 
 
 function cleanRotatorCuffSegmentationWorkspace(obj)
   rotatorCuffSegmentationPath = getConfig().pythonDir;
   delete(fullfile(rotatorCuffSegmentationPath,'input','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'IS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SC','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'SS','*'));
   delete(fullfile(rotatorCuffSegmentationPath,'TM','*'));
 end
 
 
 
 function sendImageToRotatorCuffSegmentationWorkspace(obj,sliceName)
   SCase = obj.shoulder.SCase;
   rotatorCuffSegmentationPath = getConfig().pythonDir;
   imageForSegmentationPath = fullfile(...
     obj.dataSlicesPath,sliceName + "_ForSegmentation.png");
   copyfile(imageForSegmentationPath,fullfile(...
     rotatorCuffSegmentationPath,'input',SCase.id+".png"));
 end
 
 
 
 function callRotatorCuffSegmentation(obj)
   %% Apply UNIBE method to automatically segment normal muscle
 
   pythonDir = getConfig("convertTextToString", false).pythonDir;
 
   pythonCommandActivateEnvironment= ['source ' fullfile(pythonDir,'venv','bin','activate') ';'];
   pythonCommandMoveToSegmentationWorkspace = ['cd ' pythonDir ';'];
   pythonCommandExecuteSegmentation = [fullfile(pythonDir,'rcseg.py') ' segment input' ';'];
 
   [~, pythonCommandOutput] = system([pythonCommandActivateEnvironment,...
     pythonCommandMoveToSegmentationWorkspace,...
     pythonCommandExecuteSegmentation]);
   disp(pythonCommandOutput);
 end
 
 
 
 function saveSegmentationResults(obj,maskName)
   SCaseID = obj.shoulder.SCase.id;
   rotatorCuffSegmentationPath = getConfig().pythonDir;
   
   for muscleName = ["SC", "SS", "IS", "TM"]
     try
       copyfile(fullfile(rotatorCuffSegmentationPath,muscleName,SCaseID + ".png"),...
-        fullfile(obj.(muscleName).dataMaskPath,maskName + "_Mask.png") );
+        fullfile(obj.(muscleName).dataMaskPath,maskName + "_Segmentation.png") );
     catch ME
       warning(ME.message)
     end
   end
 end
diff --git a/ShoulderCase/@RotatorCuff/sliceAndSegment.m b/ShoulderCase/@RotatorCuff/sliceAndSegment.m
new file mode 100644
index 0000000..687966c
--- /dev/null
+++ b/ShoulderCase/@RotatorCuff/sliceAndSegment.m
@@ -0,0 +1,4 @@
+function sliceAndSegment(obj)
+  obj.createSliceMatthieu();
+  obj.segmentMuscles("rotatorCuffMatthieu", "autoMatthieu");
+end
\ No newline at end of file
diff --git a/setup.m b/setup.m
index f9d7091..99d1de4 100644
--- a/setup.m
+++ b/setup.m
@@ -1,42 +1,39 @@
 function setup()
     % to execute once after cloning the repo.
     createDefaultConfigFile();
 
     addpath(genpath(pwd));
 end
 
 
 
 function createDefaultConfigFile()
     defaultConfig.maxSCaseIDDigits = 3;
     defaultConfig.SCaseIDValidTypes = ["N", "P"];
     defaultConfig.pythonDir = "/shoulder/methods/python/rcseg";
     defaultConfig.dataDir = "/shoulder/dataDev";
     
     defaultConfig.shouldersToMeasure.rightAuto = true;
     defaultConfig.shouldersToMeasure.rightManual = true;
     defaultConfig.shouldersToMeasure.leftAuto = true;
     defaultConfig.shouldersToMeasure.leftManual = true;
     
     defaultConfig.standardMeasurementsToRun.loadData = true;
     defaultConfig.standardMeasurementsToRun.morphology = true;
     defaultConfig.standardMeasurementsToRun.measureFirst = true;
     defaultConfig.standardMeasurementsToRun.measureSecond = true;
     defaultConfig.standardMeasurementsToRun.measureThird = true;
 
     defaultConfig.specialMeasurementsToRun.sliceAndSegment = false;
-    defaultConfig.specialMeasurementsToRun.measureDegenerations = false;
     defaultConfig.specialMeasurementsToRun.calcDensity = false;
 
     defaultConfig.specialMeasurementsArguments.sliceAndSegment = [];
-    defaultConfig.specialMeasurementsArguments.measureDegenerations =...
-        ["rotatorCuffMatthieu", "autoMatthieu"];
     defaultConfig.specialMeasurementsArguments.calcDensity = [];
 
     defaultConfig.overwriteMeasurements = false;
     defaultConfig.saveMeasurements = false;
 
     fid = fopen('config.json','w');
     fprintf(fid, replace(jsonencode(defaultConfig), ",", ",\n"));
     fclose(fid);
 end
\ No newline at end of file
diff --git a/utils/@Sphere/Sphere.m b/utils/@Sphere/Sphere.m
index 2e7bd92..b25c686 100644
--- a/utils/@Sphere/Sphere.m
+++ b/utils/@Sphere/Sphere.m
@@ -1,59 +1,44 @@
 classdef Sphere < handle
 % Simple sphere object.
 % Used to be fitted to glenoid points.
 % 
 % This class could be used elsewhere and the constructor
 % access might be widen.
 
   properties
     center = [];
     radius = [];
     residuals = [];
     R2 = [];
     RMSE = [];
   end
 
   methods
     function obj = Sphere(varargin)
       if nargin == 2
         obj.center = varargin{1};
         obj.radius = varargin{2};
       end
     end
 
     function fitTo(obj,points)
       [center,radius,residuals,R2] = fitSphere(points);
       obj.center            = center';
       obj.radius            = radius;
       obj.residuals         = residuals;
       obj.R2                = R2;
       obj.RMSE              = norm(residuals)/sqrt(length(residuals));
     end
 
     function output = isempty(obj)
       output = any(cellfun(@(field)isempty(obj.(field)),fields(obj)));
     end
 
-    function output = polarTangentPoint(obj, polePoint, slicePlaneNormal,...
-      polarPointOrientation)
-      poleToCenter = Vector(polePoint, obj.center);
-      assert(poleToCenter.norm > obj.radius, "Pole point is inside sphere and doesn't have any tangent.");
-
-      poleToCenterPerpendicular = poleToCenter.cross(slicePlaneNormal).orientToward(polarPointOrientation).direction;
-      tangentPointAngle = acos(obj.radius/poleToCenter.norm);
-
-      polarTangentPoint = obj.center + ...
-        obj.radius*sin(tangentPointAngle)*poleToCenterPerpendicular + ...
-        obj.radius*cos(tangentPointAngle)*(-poleToCenter.direction);
-      
-      output = polarTangentPoint;
-    end
-
     function exportPly(obj, filename)
       sphereToExport = BinarySphere(round(obj.radius), round(obj.center));
       sphereToExport.exportPly(filename);
     end
  
   end
 
 end
diff --git a/utils/@SphereContactPoint/SphereContactPoint.m b/utils/@SphereContactPoint/SphereContactPoint.m
new file mode 100644
index 0000000..d12e893
--- /dev/null
+++ b/utils/@SphereContactPoint/SphereContactPoint.m
@@ -0,0 +1,95 @@
+classdef SphereContactPoint < handle 
+
+  properties (Access = private)
+    sphere 
+    anchorPoint = []
+    polarPoint = []
+  end
+
+  methods 
+    function obj = SphereContactPoint(sphere)
+      obj.sphere = sphere;
+    end
+
+    function set.sphere(obj, value)
+      assert(isequal(class(value), "Sphere"), "Must be a Sphere instance");
+      obj.sphere = value;
+    end
+
+    function set.anchorPoint(obj, value)
+      assert(isnumeric(value) & isequal(size(value), [1 3]),...
+        "Must be a 1x3 numeric array")
+      assert(not(isequal(value, obj.sphere.center)),...
+        "Anchor point can't be superposed with sphere center");
+      obj.anchorPoint = value;
+    end
+
+    function set.polarPoint(obj, value)
+      assert(isnumeric(value) & isequal(size(value), [1 3]),...
+        "Must be a 1x3 numeric array")
+      assert(not(Vector(obj.sphere.center, value).norm <= obj.sphere.radius),...
+        "Polar point must be outside the sphere");
+      obj.polarPoint = value;
+    end
+
+    function setSphere(obj, sphere)
+      obj.sphere = sphere;
+    end
+    
+    function setAnchorPoint(obj, point)
+      obj.anchorPoint = point;
+    end
+    
+    function setPolarPoint(obj, point)
+      obj.polarPoint = point;
+    end
+
+    function output = getContactPoint(obj)
+      assert(not(isempty(obj.anchorPoint)), "Anchor point missing");
+      assert(not(isempty(obj.polarPoint)), "Polar point missing");
+      assert(obj.sphereIsBetweenAnchorAndPolarPoints,...
+        "There is no contact point");
+
+      poleToCenter = Vector(obj.polarPoint, obj.sphere.center);
+      anchorToCenter = Vector(obj.anchorPoint, obj.sphere.center);
+      sliceNormal = poleToCenter.cross(anchorToCenter);
+
+      poleToCenterPerpendicular = poleToCenter.cross(sliceNormal);
+      % angle at vertex sphere.center between contact point and polar point
+      contactPointsAngle = acos(obj.sphere.radius/poleToCenter.norm);
+
+      firstContactPoint = obj.sphere.center...
+        + obj.sphere.radius*sin(contactPointsAngle)*poleToCenterPerpendicular.direction...
+        + obj.sphere.radius*cos(contactPointsAngle)*(-poleToCenter.direction);
+
+      secondContactPoint = obj.sphere.center...
+        + obj.sphere.radius*sin(contactPointsAngle)*(-poleToCenterPerpendicular.direction)...
+        + obj.sphere.radius*cos(contactPointsAngle)*(-poleToCenter.direction);
+
+      output = obj.getContactPointClosestToAnchor(...
+        firstContactPoint, secondContactPoint);
+    end
+
+    function output = getContactPointClosestToAnchor(obj,...
+        firstContactPoint, secondContactPoint)
+      anchorToFirstContactPoint = Vector(obj.anchorPoint, firstContactPoint);
+      anchorToSecondContactPoint = Vector(obj.anchorPoint, secondContactPoint);
+
+      if anchorToFirstContactPoint.norm < anchorToSecondContactPoint.norm
+        output = firstContactPoint;
+      else
+        output = secondContactPoint;
+      end
+    end
+
+    function output = sphereIsBetweenAnchorAndPolarPoints(obj)
+      poleToCenter = Vector(obj.polarPoint, obj.sphere.center);
+      polarApertureAngle = atan(obj.sphere.radius / poleToCenter.norm);
+      poleToAnchor = Vector(obj.polarPoint, obj.anchorPoint);
+      sphereCenterToAnchorAngle = angle(poleToCenter, poleToAnchor);
+
+      output = sphereCenterToAnchorAngle < polarApertureAngle;
+    end
+  end
+
+end
\ No newline at end of file