Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/muscleDegeneration.m
4e0ab13607fcaterrier-plotUpdate

Last Change
Raw File

muscleDegeneration.m
View Options

function [Stot,Sm,Sa,Si,So,Ra,Ri,Ro,Rd] = muscleDegeneration(subject,muscleName,HUThresholdMuscle,HUThresholdOsteochondroma,output)
%%MUSCLEDEGENERATION computes degeneration parameters from muscle contours
%
% This function is called by "degenerationAll.m" and returns muscle atrophy,
% fat infiltration, osteochondroma and degeneration values.
%
% USE: muscleDegeneration(subject,muscleName,HUThresholdMuscle,HUThresholdOsteochondroma,output)
%
% INPUT subject: patient number
% muscleName: name of the rotator cuff muscle being measured (SS,IS, SC or TM)
% HUThresholdMuscle: HU threshold used for muscle segmentation (muscle-fat limit)
% HUThresholdMuscle:HU threshold used for osteochondroma segmentation (muscle-osteochondroma limit)
% output: type of output returned (image)
%
% OUTPUT Stot: total area of the muscle fossa
% Sm: area of the degenerated muscle, without atrophy, fat infiltration and osteochondroma
% Sa: area of the atrophy (Stot - Sm)
% Si: area of fat infiltration
% So: area of osteochondroma
% Ra: ratio of muscle atrophy (area atrophy over total area of muscle fossa)
% Ri: ratio of fat infiltration in the muscle (area fat infiltration over total area of muscle fossa)
% Ro: ratio of osteochondroma in the muscle (area osteochondroma over total area of muscle fossa)
% Rd: the ratio of degeneration of the muscle (sum of atrophy, fat infiltration and osteochondroma over total area of muscle fossa)
%
% REMARKS The degeneration values are saved in .dat files.
%
% created with MATLAB ver.: 9.2.0.556344 (R2017a) on Windows 7
% Author: EPFL-LBO-VMC
% Date: 07-Jun-2017
%
close all
fontSize = 14;
% visualization window
visualizationWindow = [-100 200];
% Load data
inputPath = sprintf('results/%s/%s/', muscleName, subject);
inputFile = sprintf('ssContours%s.mat',subject);
A = load(sprintf('%s%s',inputPath,inputFile));
% muscle segmentation
imageWindow = double([min(A.inputImage(:)) max(A.inputImage(:))]);
normalizedThresholdMuscle = (HUThresholdMuscle-imageWindow(1))/(imageWindow(2)-imageWindow(1));
normalizedThresholdOsteochondroma = (HUThresholdOsteochondroma-imageWindow(1))/(imageWindow(2)-imageWindow(1));
% normalize image
normalizedInputImage = double(A.inputImage)-imageWindow(1);
normalizedInputImage = normalizedInputImage/(imageWindow(2)-imageWindow(1));
blackMaskedImage = normalizedInputImage;
blackMaskedImage(~A.binaryImage1) = 0; % If outside contours -> black
segmented = imbinarize(blackMaskedImage,normalizedThresholdMuscle);
myMuscle = imfill(segmented, 'holes'); % Fill holes
myMuscle(~A.binaryImage1) = 0;
% segmentation mask for fat infiltration - select pixels above threshold
% (exclude fat, include muscle)
blackMaskedImage = normalizedInputImage;
blackMaskedImage(~myMuscle) = 0; % If outside muscle selection -> 0
binaryImage3 = imbinarize(blackMaskedImage,normalizedThresholdMuscle);
% segmentation osteochondroma
% select pixels above threshold (exclude muscle, include osteochondroma)
myOsteo = imbinarize(blackMaskedImage,normalizedThresholdOsteochondroma);
myOsteo = imfill(myOsteo, 'holes');
myOsteo(~myMuscle) = 0;
% Measure atrophy and degeneration
% Define different zones
% Stot = manually delimited contours
% Sa = atrophy zone (around muscle)
% Sm = remaining muscle (muscle only)
% Si = infiltration
% So = osteochondroma
% so that Stot = Sa + Sm + Si + So
Stot = bwarea(A.binaryImage1);
Sa = bwarea(A.binaryImage1)-bwarea(myMuscle);
So = bwarea(myOsteo);
Sm = bwarea(binaryImage3)-bwarea(myOsteo);
Si = bwarea(myMuscle)-Sm-So;
%Old definitions
% % Calculate SS atrophy
% ssAtrophy = 100*(1-bwarea(myMuscle)/bwarea(A.binaryImage1));
%
% % Calculate SS infiltration
% ssInfiltration = 100*(1-bwarea(binaryImage3)/bwarea(myMuscle));
%
% % Calculate Osteochondroma
% ssOsteochondroma = 100*(bwarea(myOsteo)/bwarea(myMuscle));
%New definitions
% % Calculate atrophy
Ra = Sa/Stot;
% % Calculate infiltration
Ri = Si/Stot;
% % Calculate osteochondroma
Ro = So/Stot;
% % Calculate degeneration
Rd = (Stot-Sm)/Stot;
% Report results
if exist('output','var') == 1
if strcmp(output,'image')
imshow(A.inputImage,visualizationWindow);
[r,c,~] = size(A.inputImage);
title(sprintf('Muscle threshold: %.0fHU, Muscle atrophy: %.2f%%, Muscle infiltration: %.2f%%, Osteochondroma: %.2f%%', ...
HUThresholdMuscle,Ra*100,Ri*100,Ro*100), 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize'), 'Visible', 'on'); % Maximize figure.
hold on
% Get coordinates and plot of the boundaries drawn regions.
structBoundaries1 = bwboundaries(A.binaryImage1);
xy1=structBoundaries1{1}; % Get n by 2 array of x,y coordinates.
x1 = xy1(:, 2); % Columns.
y1 = xy1(:, 1); % Rows.
plot(x1, y1, 'color','green');
% Plot the selected region (fat+muscle) within the drawn regions
structBoundaries2 = bwboundaries(myMuscle);
xy2=structBoundaries2{1}; % Get n by 2 array of x,y coordinates.
x2 = xy2(:, 2); % Columns.
y2 = xy2(:, 1); % Rows.
plot(x2, y2, 'color','blue');
% Plot the muscle only within the fat+muscle region
structBoundaries3 = bwboundaries(binaryImage3);
for i = 1:size(structBoundaries3,1)
xy3 = structBoundaries3{i}; % Get n by 2 array of x,y coordinates.
x3 = xy3(:, 2); % Columns.
y3 = xy3(:, 1); % Rows.
%subplot(2, 3, 3); % Plot over original image.
hold on; % Keep the image.
plot(x3, y3,'color','red');
end
% Plot the osteochondroma within the muscle only region
structBoundaries4 = bwboundaries(myOsteo);
for i = 1:size(structBoundaries4,1)
xy3 = structBoundaries4{i}; % Get n by 2 array of x,y coordinates.
x3 = xy3(:, 2); % Columns.
y3 = xy3(:, 1); % Rows.
%subplot(2, 3, 3); % Plot over original image.
hold on; % Keep the image.
plot(x3, y3,'color','magenta');
end
outputFile = sprintf('ssDegeneration%s.%d',subject,HUThresholdMuscle);
saveas(gcf,sprintf('%s%s.png',inputPath,outputFile))
set(gca,'Units','normalized','Position',[0 0 1 1])
set(gcf,'Units','pixels','Position',[200 200 c r]);
f = getframe(gcf);
imwrite(f.cdata,sprintf('%s%s.jpg',inputPath,outputFile));
dlmwrite(sprintf('%s%s.dat',inputPath,outputFile),...
[HUThresholdMuscle; Ra*100; Ri*100; Ro*100; Rd*100],'precision',6);
disp(fprintf('%s%s.png\n%s%s.dat'...
,inputPath,outputFile,inputPath,outputFile))
close all
end
end
Ra = num2cell(Ra);
Ri = num2cell(Ri);
Ro = num2cell(Ro);
Rd = num2cell(Rd);
end
c4science · Help