Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (master) Browse

R8218/muscleDegeneration.m
8d2978dd4d8cmaster

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
%
% Load the reference 2D CT slice and associated muscle contour
CTDatabaseLocation = 'Z://data'; % Location of the CT database
listDirInCurrDir = dir([CTDatabaseLocation '/' subject(1) '/' subject(2) '/' subject(3)]);
listDirInCurrDir = listDirInCurrDir(3:end);
for i=1:length(listDirInCurrDir)
if ~isempty(regexp(listDirInCurrDir(i).name,['^[PN]\d{0,2}[' subject(4) '][-]\w*'], 'once'))
listDirInCaseFolder = dir([listDirInCurrDir(i).folder '\' listDirInCurrDir(i).name]);
listDirInCaseFolder = listDirInCaseFolder(3:end);
% In all CT folders in the current Case Folder, find CT #1
for j=1:length(listDirInCaseFolder)
if ~isempty(regexp(listDirInCaseFolder(j).name,'^[C][T][-][PN]\d{1,3}[-]\d*[-][1]', 'once'))
% If 'muscles' exists, save iput path
if exist([listDirInCaseFolder(j).folder '\' listDirInCaseFolder(j).name '\muscles'],'dir') == 7
inputPath = [listDirInCaseFolder(j).folder '\' listDirInCaseFolder(j).name '\muscles\' muscleName '\'];
inputFile = ['ssContours' subject '.mat'];
else
error(['Missing "muscles" directory for CaseID: ' subject]);
end
end
end
end
end
contouredMuscle = load(sprintf('%s%s',inputPath,inputFile));
% --Atrophied muscle segmentation--
% Normalize image and thresholds
rangeHU = double([min(contouredMuscle.inputImage(:)) max(contouredMuscle.inputImage(:))]);
normalizedThresholdMuscle = (HUThresholdMuscle-rangeHU(1))/(rangeHU(2)-rangeHU(1));
normalizedThresholdOsteochondroma = (HUThresholdOsteochondroma-rangeHU(1))/(rangeHU(2)-rangeHU(1));
normalizedInputImage = double(contouredMuscle.inputImage)-rangeHU(1);
normalizedInputImage = normalizedInputImage/(rangeHU(2)-rangeHU(1));
% Mask image using muscle fossa contours
maskedImage = normalizedInputImage;
maskedImage(~contouredMuscle.binaryImage1) = 0;
% Segment using muscle threshold
myMuscle = imbinarize(maskedImage,normalizedThresholdMuscle);
myMuscle = imfill(myMuscle, 'holes'); % Fill holes
% --Atrophied muscle without fat infiltration segmentation--
% Mask image using atrophied muscle segmentation
maskedImage = normalizedInputImage;
maskedImage(~myMuscle) = 0;
% Segment using muscle threshold
myMuscleNoFat = imbinarize(maskedImage,normalizedThresholdMuscle);
% --Osteochondroma segmentation--
%Segment using osteochondroma threshold
myOsteo = imbinarize(maskedImage,normalizedThresholdOsteochondroma);
myOsteo = imfill(myOsteo, 'holes'); % Fill holes
% --Measure atrophy and degeneration so that Stot = Sa + Sm + Si + So--
Stot = bwarea(contouredMuscle.binaryImage1); % Stot = manually delimited contours
Sa = bwarea(contouredMuscle.binaryImage1) - bwarea(myMuscle); % Sa = atrophy zone (around muscle)
So = bwarea(myOsteo); % So = osteochondroma
Sm = bwarea(myMuscleNoFat) - bwarea(myOsteo); % Sm = remaining muscle (muscle only)
Si = bwarea(myMuscle) - Sm - So; % Si = fat infiltration
Ra = Sa/Stot; % Atrophy ratio
Ro = So/Stot; % Osteochondroma ratio
Ri = Si/Stot; % Fat infiltration ratio
Rd = (Stot-Sm)/Stot; % Degeneration
% Report results
if exist('output','var') == 1
if strcmp(output,'image')
% --Plot CT slice with all contours--
fontSize = 14;
visualizationWindow = [-100 200];
imshow(contouredMuscle.inputImage,visualizationWindow);
[r,c,~] = size(contouredMuscle.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
% 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.
plot(x, y, 'color','green');
% Plot atrophied muscle contours
structBoundaries = bwboundaries(myMuscle);
xy=structBoundaries{1}; % Get n by 2 array of x,y coordinates.
x = xy(:, 2); % Columns.
y = xy(:, 1); % Rows.
plot(x, y, 'color','blue');
% Plot fat infiltration contours
structBoundaries = bwboundaries(myMuscleNoFat);
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.
plot(x, y,'color','red');
end
% Plot osteochondroma contours
structBoundaries = bwboundaries(myOsteo);
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.
plot(x, y,'color','magenta');
end
% --Save images--
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));
close all
% --Write to file--
dlmwrite(sprintf('%s%s.dat',inputPath,outputFile),...
[HUThresholdMuscle; Ra*100; Ri*100; Ro*100; Rd*100],'precision',6);
end
end
Ra = num2cell(Ra);
Ri = num2cell(Ri);
Ro = num2cell(Ro);
Rd = num2cell(Rd);
end
c4science · Help