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
    segmentationMask
    segmentationDensity
    slicePixelSpacings

    rangeHU = double([-1000 1000]);
    muscleThreshold = 0;
    fatThreshold = 30;
    osteochondromaThreshold = 166;
  end



  methods

    function obj = MuscleMeasurer(segmentationMask, sliceForMeasurement, slicePixelSpacings)
      obj.segmentationMask = segmentationMask;
      obj.segmentationDensity =...
        obj.getNormalizedImage(sliceForMeasurement).*segmentationMask;
      obj.slicePixelSpacings = slicePixelSpacings;

      obj.normalizeThresholds;
    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)
      pixelSurface = (obj.slicePixelSpacings(1) * obj.slicePixelSpacings(2)) / 100; % cm^2
      output = pixelSurface * sum(obj.segmentationMask, "all");
    end

    function output = getRatioAtrophy(obj)
      output = obj.getRatioAreas(obj.getAreaAtrophy,obj.segmentationMask);
    end

    function output = getRatioFat(obj)
      output = obj.getRatioAreas(obj.getAreaFat,obj.segmentationMask);
    end

    function output = getRatioOsteochondroma(obj)
      output = obj.getRatioAreas(obj.getAreaOsteochondroma,obj.segmentationMask);
    end

    function output = getRatioDegeneration(obj)
      output = obj.getRatioAtrophy + obj.getRatioFat + obj.getRatioOsteochondroma;
    end

    function output = getRatioAreas(obj,partialArea,totalArea)
      output = sum(partialArea, "all")/sum(totalArea, "all");
    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.segmentationDensity,obj.muscleThreshold);
      areaMuscle = imfill(areaMuscle,'holes');
      % Keep biggest island only
      islands = bwconncomp(areaMuscle);
      areaMuscle = false(size(areaMuscle));
      if not(isempty(islands.PixelIdxList))
        islandsSizes = cellfun(@numel,islands.PixelIdxList);
        [~,biggestIslandIndex] = max(islandsSizes);
        areaMuscle(islands.PixelIdxList{biggestIslandIndex}) = true;
      end

      output = areaMuscle;
    end

    function output = getAreaAtrophy(obj)
      output = obj.segmentationMask & not(obj.getAreaMuscle);
    end

    function output = getAreaFat(obj)
      muscleImage = obj.segmentationDensity.*obj.getAreaMuscle;
      areaFat = obj.getAreaMuscle & not(imbinarize(muscleImage,obj.fatThreshold));

      output = areaFat;
    end

    function output = getAreaOsteochondroma(obj)
      muscleImage = obj.segmentationDensity.*obj.getAreaMuscle;
      areaOsteochondroma = imbinarize(muscleImage,obj.osteochondromaThreshold);
      areaOsteochondroma = imfill(areaOsteochondroma,'holes');

      output = areaOsteochondroma;
    end

  end

end