function muscle = muscleDegeneration(caseID,muscleName,HUThresholdMuscle,HUThresholdFat,HUThresholdOsteochondroma,muscleDirectory)
%%MUSCLEDEGENERATION computes degeneration parameters from muscle contours
%
% This function returns muscle atrophy, fat infiltration, osteochondroma
% and degeneration values for a given muscle
%
% USE: muscleDegeneration(caseID,muscleName,HUThresholdMuscle,HUThresholdOsteochondroma,muscleDirectory)
%
% INPUT     caseID: patient number
%           muscleName: name of the rotator cuff muscle being measured (SS,IS, SC or TM)
%           HUThresholdMuscle: HU threshold used for muscle segmentation (atrophied muscle limit)
%           HUThresholdFat: HU threshold used for fat infiltration segmentation (fat-muscle limit)
%           HUThresholdOsteochondroma:HU threshold used for osteochondroma segmentation (muscle-osteochondroma limit)
%           muscleDirectory: Path to 'muscles' directory where contours and output images are stored
%
% OUTPUT    muscle: structure containing fields 'Stot', 'Satrophy', 
%           'Sinfiltration', 'Sosteochondroma' and 'Sdegeneration', 
%           'atrophy', 'infiltration', 'osteochondroma' and 'degeneration'
%
% REMARKS   An image of the segmentation of atrophied muscle, fat 
%           infiltration and osteochondroma is saved in muscles directory, 
%           in the 'muscleName' folder.
%
% created with MATLAB ver.: 9.2.0.556344 (R2017a) on Windows 7
% Author: EPFL-LBO-VMC
% Date: 27-Jun-2017
%

inputFile = ['ssContours' caseID '.mat'];
contouredMuscle = load([muscleDirectory '\' muscleName '\' inputFile]);

% --Plot CT slice with all contours--

fontSize = 14;
visualizationWindow = [-100 200];

imshow(contouredMuscle.inputImage,visualizationWindow);

hold on

% --Atrophied muscle segmentation--

% 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. 

% Normalize image and thresholds

contouredMuscle.inputImage(contouredMuscle.inputImage < -1000) = -1000;
contouredMuscle.inputImage(contouredMuscle.inputImage > 1000) = 1000;

rangeHU = double([-1000 1000]);

% Binarization of images works so that pixel value > threshold is kept. So
% if you condier that a msucle is normal for 30 HU and abnormal for 29 HU,
% the threshold should be set at 29 HU. 
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));

normalizedInputImage = (double(contouredMuscle.inputImage)-rangeHU(1))/(rangeHU(2)-rangeHU(1));

% Mask image using muscle fossa contours
maskedImage = normalizedInputImage;
maskedImage(~contouredMuscle.binaryImage1) = 0;

% Plot muscle fossa contours.
structBoundaries = bwboundaries(contouredMuscle.binaryImage1);
xy=structBoundaries{1}; % Get n by 2 array of x,y coordinates.
x = xy(:, 2); % Columns.
y = xy(:, 1); % Rows.
p = plot(x, y, 'color','green','DisplayName','Muscle fossa');
p_count = 1;

% Segment using muscle threshold
myMuscle = imbinarize(maskedImage,normalizedThresholdMuscle);
myMuscle = imfill(myMuscle, 'holes'); % Fill holes

% Keep only the biggest connected component (i.e. remove islands)
CC = bwconncomp(myMuscle);
numPixels = cellfun(@numel,CC.PixelIdxList);
[~,idx] = max(numPixels);
myMuscle(:,:) = 0;
myMuscle(CC.PixelIdxList{idx}) = 1;

% Plot atrophied muscle contours
structBoundaries = bwboundaries(myMuscle);

if size(structBoundaries,1) > 0
    p_count = p_count+1;
end

for i = 1:size(structBoundaries,1)
    xy = structBoundaries{i}; % Get n by 2 array of x,y coordinates.
    x = xy(:, 2); % Columns.
    y = xy(:, 1); % Rows.
    hold on; % Keep the image.
    p(p_count) = plot(x, y, 'color','blue','DisplayName','Atrophied Muscle');
end

% --Fat infiltration segmentation--

% Mask image using atrophied muscle segmentation
maskedImage = normalizedInputImage;
maskedImage(~myMuscle) = 0;

% Segment using fat threshold
myMuscleNoFat = imbinarize(maskedImage,normalizedThresholdFat);
myFat = ~myMuscleNoFat;
myFat(~myMuscle) = 0;

% Plot fat infiltration surfaces
structBoundaries = bwboundaries(myFat);

% The 'fill' function cannot handle hollow polygons (with holes = 2
% boundaries). To overcome this difficulty, fat infiltration is plotted as
% a coloured overlay image with transparency zero outside fat regions and 
% transparency 0.25 in fat regions.

myFatAlpha = double(myFat);
myFatAlpha(myFat) = 0.25;
red=zeros(size(myFat,1),size(myFat,2),3);
red(:,:,1)=1;
h=imshow(red);
set(h,'AlphaData',myFatAlpha);

if size(structBoundaries,1) > 0
    p_count = p_count+1;
    p(p_count) = fill(0, 0,'red','EdgeColor','none','FaceAlpha', 0.25,'DisplayName','Fat infiltration'); % plot a filled patch outside the plot region to add to the legend
end

% --Osteochondroma segmentation--

%Segment using osteochondroma threshold
myOsteo = imbinarize(maskedImage,normalizedThresholdOsteochondroma);
myOsteo = imfill(myOsteo, 'holes'); % Fill holes

% Plot osteochondroma surfaces
structBoundaries = bwboundaries(myOsteo);

% The 'fill' function cannot handle hollow polygons (with holes = 2
% boundaries). To overcome this difficulty, osteochondroma are plotted as
% a coloured overlay image with transparency zero outside osteochondroma  
% and transparency 0.25 in osteochondroma regions.

myOsteoAlpha = double(myOsteo);
myOsteoAlpha(myOsteo) = 0.5;
yellow=zeros(size(myOsteo,1),size(myOsteo,2),3);
yellow(:,:,1)=1;
yellow(:,:,2)=1;
h=imshow(yellow);
set(h,'AlphaData',myOsteoAlpha);

if size(structBoundaries,1) > 0
    p_count = p_count+1;
    p(p_count) = fill(0, 0,'yellow','EdgeColor','none','FaceAlpha', 0.5,'DisplayName','Osteochondroma'); % plot a filled patch outside the plot region to add to the legend
end

% --Measure atrophy and degeneration so that Stot = Sa + Sm + Si + So--

muscle.caseID = caseID;
muscle.Stot = bwarea(contouredMuscle.binaryImage1); % Stot: total area of the muscle fossa (manually delimited contours)
muscle.Satrophy = bwarea(contouredMuscle.binaryImage1) - bwarea(myMuscle); % Sa: area of the atrophy (Stot - Sm)
muscle.Sinfiltration = bwarea(myFat); % Si: area of fat infiltration
muscle.Sosteochondroma = bwarea(myOsteo); %  So: area of osteochondroma
muscle.Sdegeneration = muscle.Stot-muscle.Satrophy-muscle.Sinfiltration-muscle.Sosteochondroma; % Sm: area of the degenerated muscle, without atrophy, fat infiltration and osteochondroma
      
muscle.atrophy = muscle.Satrophy/muscle.Stot; % ratio of muscle atrophy (area atrophy over total area of muscle fossa)
muscle.infiltration = muscle.Sinfiltration/muscle.Stot; % ratio of fat infiltration in the muscle (area fat infiltration over total area of muscle fossa)
muscle.osteochondroma = muscle.Sosteochondroma/muscle.Stot; % ratio of osteochondroma in the muscle (area osteochondroma over total area of muscle fossa)
muscle.degeneration = (muscle.Stot-muscle.Sdegeneration)/muscle.Stot; % the ratio of degeneration of the muscle (sum of atrophy, fat infiltration and osteochondroma over total area of muscle fossa)


% --Save images--

title({['Case ID: ' caseID]; [muscleName ' Muscle threshold: ' sprintf('%d',HUThresholdMuscle) ' HU']; sprintf('Muscle atrophy: %.2f%%,  Muscle infiltration: %.2f%%, Osteochondroma: %.2f%%', ...
    muscle.atrophy*100,muscle.infiltration*100,muscle.osteochondroma*100)}, 'FontSize', fontSize);
fig = get(groot,'CurrentFigure');
fig.OuterPosition = [-1140 -1 871 1028];
hLegend = legend(p);
set(hLegend.BoxFace, 'ColorType','truecoloralpha', 'ColorData',uint8(255*[.8;.8;.8;.8]));  
hLegend.EdgeColor = 'none';

% % This piece of code can be used to make the rectangle patches in the
% legend have the correct alpha value (transparency) as used in the plot.
% Unfortunately, it doesnt work in MATLAB r2017a because calling the
% 'legend' function with multiple output arguments results in a graphics
% bug and prevent the rectangular patches to display. I keep this piece of
% code here, as this bug might be solved in a future version.
%
% VMC-29-Jun-2017
%
% % Get handles of all patch objects in current axis: 
% hla = findobj(gca,'type','patch'); 
% 
% % Get alpha properties of the current plotted patches: 
% hlaFA = get(hla,'FaceAlpha');
% hlaFC = get(hla,'FaceColor');
% 
% % Create a legend: 
% [hLegend,oLegend,~,~]= legend(p);
% set(hLegend.BoxFace, 'ColorType','truecoloralpha', 'ColorData',uint8(255*[.8;.8;.8;.8]));  
% hLegend.EdgeColor = 'none';
% 
% % Get handles of patch objects in the legend: 
% hlp = findobj(oLegend,'type','patch'); 
% 
% % Set alpha values: 
% for k = 1:length(hlp)
% set(hlp(k),'facec',hlaFC(k,:)); 
% set(hlp(k),'facea',hlaFA(k)); 
% end

outputFile = sprintf('%sDegeneration%s_%d',char(muscleName),char(caseID),HUThresholdMuscle);
% saveas(gcf,sprintf('%s%s.png',[muscleDirectory '/' muscleName '/'],outputFile))

close all

end