Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/muscles/muscleDegeneration3D.m
7ac612db8877aterrier-plotUpdate

Last Change
Raw File

muscleDegeneration3D.m
View Options

function muscle = muscleDegeneration3D(caseID,muscleName,HUThresholdMuscle,HUThresholdFat,HUThresholdOsteochondroma,muscleDirectory)
%%MUSCLEDEGENERATION3D computes degeneration parameters from muscle contours
%
% This function returns muscle atrophy, fat infiltration, osteochondroma
% and degeneration values computed from a CT slice for a given muscle
%
% USE: muscleDegeneration3D(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: 29-Jun-2017
%
% Get slice number
segments = textscan(muscleDirectory,'%s','Delimiter','\');
segments = segments{1};
sliceNum = segments(end);
% Find contours
inputFile = ['ssContours' caseID '_' sliceNum{1} '.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.SliceNumber = str2double(sliceNum);
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 ', Slice #: ' sliceNum{1}]; [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
sliceNum = sprintf(num2str(str2double(sliceNum{1}),'%03i'));
outputFile = sprintf('%sDegeneration%s_%s',muscleName,caseID,sliceNum);
saveas(gcf,sprintf('%s%s.png',[muscleDirectory '\'],outputFile));
close all
end
c4science · Help