Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/muscleDegeneration.m
260f40b8f635aterrier-plotUpdate

Last Change
Raw File

muscleDegeneration.m
View Options

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,output)
%
% 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 contours
structBoundaries = bwboundaries(myFat);
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','red','DisplayName','Fat infiltration');
end
% --Osteochondroma segmentation--
%Segment using osteochondroma threshold
myOsteo = imbinarize(maskedImage,normalizedThresholdOsteochondroma);
myOsteo = imfill(myOsteo, 'holes'); % Fill holes
% Plot osteochondroma contours
structBoundaries = bwboundaries(myOsteo);
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','magenta','DisplayName','Osteochondroma');
end
% --Measure atrophy and degeneration so that Stot = Sa + Sm + Si + So--
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];
legend(p);
outputFile = sprintf('%sDegeneration%s_%d',char(muscleName),char(caseID),HUThresholdMuscle);
saveas(gcf,sprintf('%s%s.png','C:\Users\vmalfroy\Google Drive\DegenerationCuffStudy_VMC\improve_segmentation\21-Jun-2017_ActiveContours1\',outputFile))
close all
end
c4science · Help