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