diff --git a/ShoulderCase/@Muscle/Muscle.m b/ShoulderCase/@Muscle/Muscle.m
new file mode 100644
index 0000000..fb7dbb0
--- /dev/null
+++ b/ShoulderCase/@Muscle/Muscle.m
@@ -0,0 +1,49 @@
+classdef Muscle < handle
+
+  properties
+    name = '';
+    segmentation = '';
+    sliceName = '';
+    degeneration = [];
+    atrophy = [];
+    osteochondroma = [];
+    PCSA = [];
+  end
+
+
+
+  properties (Hidden = true)
+    dataPath = '';
+    container = [];
+  end
+
+
+  methods
+
+    function obj = Muscle(musclesContainer,muscleName,sliceWithMuscleFilename)
+      obj.container = musclesContainer;
+      obj.name = muscleName;
+      obj.sliceName = sliceWithMuscleFilename;
+    end
+
+    function propagateDataPath(obj)
+      obj.dataPath = fullfile(obj.container.dataPath,obj.name);
+
+      if not(isfolder(obj.dataPath))
+        mkdir(obj.dataPath);
+      end
+    end
+
+
+    function output = getValues(obj)
+      output = cellfun(@(property) obj.(property), properties(obj),...
+                       'UniformOutput',false);
+    end
+
+    function measure(obj,segmentationSet)
+      
+    end
+
+  end
+
+end
diff --git a/ShoulderCase/@Muscle/measureDegenerationFromSegmentationSet.m b/ShoulderCase/@Muscle/measureDegenerationFromSegmentationSet.m
new file mode 100644
index 0000000..fdd07d0
--- /dev/null
+++ b/ShoulderCase/@Muscle/measureDegenerationFromSegmentationSet.m
@@ -0,0 +1,128 @@
+function measureDegenerationFromSegmentationSet(obj,segmentationSetName)
+  measurementsResults = struct('area',[],...
+                               'degeneration',[],...
+                               'PCSA',[]);
+  SCase = obj.shoulder.SCase;
+  segmentationSetPNGFiles = dir(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,'*.png'));
+  segmentedAreasName = extractBefore({segmentationSetPNGFiles.name},'.png');
+  load(fullfile(SCase.dataCTPath,'muscles','imagesPixelSpacings.mat'));
+  obj.slicesPixelSpacings = imagesPixelSpacings;
+
+  for i = 1:length(segmentedAreasName)
+      segmentedImage = getSegmentedImage(obj,segmentationSetName,segmentedAreasName{i});
+      try
+        areas = getThresholdedAreas(segmentedImage);
+      catch ME
+        disp(ME.message);
+        continue;
+      end
+      degeneration = getDegenerationFromAreas(areas);
+      PCSA = getPCSAFromAreas(obj,areas);
+      measurementsResults.area = segmentedAreasName{i};
+      measurementsResults.degeneration = [measurementsResults.degeneration degeneration];
+      measurementsResults.PCSA = [measurementsResults.PCSA PCSA];
+      obj.(segmentedAreasName{i}) = degeneration;
+  end
+  measurementsResults.area = struct2table(measurementsResults);
+  save(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,'measurementsResults.mat'),'measurementsResults');
+end
+
+
+function output = getSegmentedImage(obj,segmentationSetName,muscleName)
+  SCase = obj.shoulder.SCase;
+  load(fullfile(SCase.dataCTPath,'muscles','imageForDegeneration.mat'));
+  normalizedImage = getNormalizedImage(imageForDegeneration);
+  mask = logical(imread(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,[muscleName '.png'])));
+
+  segmentedImage = normalizedImage.*mask;
+
+  output = segmentedImage;
+end
+
+
+
+function output = getNormalizedImage(image)
+  % Normalize image and thresholds
+  image(image < -1000) = -1000;
+  image(image > 1000) = 1000;
+  rangeHU = double([-1000 1000]);
+  output = (double(image)-rangeHU(1))/(rangeHU(2)-rangeHU(1));
+end
+
+
+
+function output = getThresholdedAreas(image)
+  thresholds = getThresholds();
+
+  % muscle
+  muscleArea = imbinarize(image,thresholds.muscle);
+  muscleArea = imfill(muscleArea,'holes');
+  % Keep biggest island only
+  islands = bwconncomp(muscleArea);
+  islandsSizes = cellfun(@numel,islands.PixelIdxList);
+  [~,biggestIslandIndex] = max(islandsSizes);
+  muscleArea = false(size(muscleArea));
+  muscleArea(islands.PixelIdxList{biggestIslandIndex}) = true;
+
+
+  muscleImage = image.*muscleArea;
+
+  % fat
+  fatArea = muscleArea & not(imbinarize(muscleImage,thresholds.fat));
+
+  % osteochondroma
+  osteochondromaArea = imbinarize(muscleImage,thresholds.osteochondroma);
+  osteochondromaArea = imfill(osteochondromaArea,'holes');
+
+  % healthy muscle
+  healthyMuscleArea = muscleArea & not(fatArea) & not(osteochondromaArea);
+
+  output.segmented = logical(image);
+  output.muscle = muscleArea;
+  output.healthyMuscle = healthyMuscleArea;
+  output.fat = fatArea;
+  output.osteochondroma = osteochondromaArea;
+end
+
+
+
+function output = getThresholds()
+  % MATLAB does not accept to threshold an image with a negative value, as is
+  % the case when using HU scale. CT scanners usually have a CT range from
+  % -1000 HU to 1000 HU or -2000 HU to 2000 HU for modern scanners. In this
+  % case, we're not using HU values lower than -1000 and higher than +1000 so
+  % all pixel intensity values of the input image< -1000 HU are set to -1000
+  % HU, and values > 1000 HU are set to 1000 HU. Then, the input image is
+  % mapped to the [0 1] range -1000->0 1000->1, as are the thresholds for
+  % muscle and osteochondroma
+
+  rangeHU = double([-1000 1000]);
+  HUThresholdMuscle = 0;
+  HUThresholdFat = 30;
+  HUThresholdOsteochondroma = 166;
+
+  normalizedThresholdMuscle = ((HUThresholdMuscle - 1)-rangeHU(1))/(rangeHU(2)-rangeHU(1));
+  normalizedThresholdFat = ((HUThresholdFat - 1)-rangeHU(1))/(rangeHU(2)-rangeHU(1));
+  normalizedThresholdOsteochondroma = ((HUThresholdOsteochondroma - 1)-rangeHU(1))/(rangeHU(2)-rangeHU(1));
+
+  output.muscle = normalizedThresholdMuscle;
+  output.fat = normalizedThresholdFat;
+  output.osteochondroma = normalizedThresholdOsteochondroma;
+end
+
+
+
+function output = getDegenerationFromAreas(areas)
+  degeneration = bwarea(areas.segmented & not(areas.healthyMuscle)) / bwarea(areas.segmented);
+
+  output = degeneration;
+end
+
+
+
+function output = getPCSAFromAreas(obj,areas)
+  pixelSurface = (obj.slicesPixelSpacings(1) * obj.slicesPixelSpacings(2));
+  PCSA = bwarea(areas.healthyMuscle) * pixelSurface;
+
+  output = PCSA;
+end
diff --git a/ShoulderCase/@MusclesContainer/MusclesContainer.m b/ShoulderCase/@MusclesContainer/MusclesContainer.m
new file mode 100644
index 0000000..9523724
--- /dev/null
+++ b/ShoulderCase/@MusclesContainer/MusclesContainer.m
@@ -0,0 +1,64 @@
+classdef MusclesContainer < handle
+
+  properties
+    summary = [];
+    list = {};
+  end
+
+
+
+  properties (Hidden = true)
+    dataPath = '';
+    shoulder
+  end
+
+
+
+  methods
+
+    function obj = MusclesContainer(shoulder)
+      obj.shoulder = shoulder;
+    end
+
+    function propagateDataPath(obj)
+      % Update current dataPath
+      % Propagate the path to objects in properties
+
+      obj.dataPath = fullfile(obj.shoulder.dataPath,'muscles');
+
+      if not(isfolder(obj.dataPath))
+        mkdir(obj.dataPath);
+      end
+
+      cellfun(@(muscle) muscle.propagateDataPath, obj.list);
+    end
+
+    function addMuscle(obj,muscleName,sliceWithMuscleFilename)
+      newMuscle = Muscle(obj,muscleName,sliceWithMuscleFilename);
+      newMuscle.propagateDataPath;
+      obj.list{end+1} = newMuscle;
+
+      obj.updateSummary;
+    end
+
+    function updateSummary(obj)
+      musclesValues = obj.getAllMusclesValues();
+      obj.summary = cell2table(musclesValues,...
+                               'VariableNames',properties('Muscle'));
+    end
+
+    function output = getAllMusclesValues(obj)
+      output = {};
+      for i = 1:length(obj.list)
+        output(end+1,:) = obj.list{i}.getValues';
+      end
+    end
+
+    function measureAllMuscles(obj,segmentationSet)
+      cellfun(@(muscle) muscle.measure(segmentationSet), obj.list);
+      obj.updateSummary;
+    end
+
+  end
+
+end