Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/muscleDegeneration.m
5233bc7a4d91aterrier-plotUpdate

Last Change
Raw File

muscleDegeneration.m
View Options

%% muscleDegeneration
%
% This script is used by "degenerationAll.m" and returns the atrophy, the infiltration, osteo and degeneration values.
%
% Author: EPFL-LBO
% Date: 2016-09-05
%
%%
% Source is 'all', 'new' or 'specific' (for only one case, the subject is required as second argument)
% Subject is the patient number (only in case of 'specific')
function [Stot,Sm,Sa,Si,So,Ra,Ri,Ro,Rd] = muscleDegeneration(subject,muscleName,HUThresholdMuscle,HUThresholdOsteochondroma,output)
close all
filtreOsteo = 1;
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));
filteredInputImage = medfilt2(normalizedInputImage);
% isolate SS muscle
bigSE = strel('disk', 5);
smallSE = strel('disk', 3);
blackMaskedImage = filteredInputImage;
blackMaskedImage(~A.binaryImage1) = 0; % If outside contours -> black
segmented = im2bw(blackMaskedImage,normalizedThresholdMuscle);
myMuscle = medfilt2(segmented); % Remove isolated points
myMuscle = imfill(myMuscle, 'holes'); % Fill holes
% Control
tempMuscle = myMuscle; small = 0;
for i=1 : 2
tempMuscle = imerode(tempMuscle, bigSE); % Erosion - Refine selection
end
% tempMuscle = imfill(tempMuscle, 'holes');
if max(max(tempMuscle)) == 0
tempMuscle = myMuscle;
tempMuscle = imerode(tempMuscle, bigSE);
% tempMuscle = imfill(tempMuscle, 'holes');
small = 1;
if max(max(tempMuscle)) == 0
tempMuscle = myMuscle;
tempMuscle = imerode(tempMuscle, smallSE);
% tempMuscle = imfill(tempMuscle, 'holes');
small = 2;
end
end
myMuscle = tempMuscle;
island = sum(sum(myMuscle)) / 10;
myMuscle = bwareaopen(myMuscle, double(uint16(island))); % Remove islands smaller than 1/10
for i=1 : 2
if small == 2
myMuscle = imdilate(myMuscle, smallSE);
break
end
myMuscle = imdilate(myMuscle, bigSE); % Dilatation/Grow image
if small == 1 break
end
end
myMuscle = bwmorph(myMuscle, 'majority'); % Smooth contours, remove small peaks ? Smoothing
myMuscle = imerode(myMuscle, smallSE);
myMuscle = imdilate(myMuscle, smallSE); % mod on 25.04
myMuscle = imfill(myMuscle, 'holes');
myMuscle(~A.binaryImage1) = 0;
%% segmentation mask for fat infiltration - select pixels above threshold
% (exclude fat, include muscle)
blackMaskedImage = filteredInputImage;
blackMaskedImage(~myMuscle) = 0; % If outside muscle selection -> 0
binaryImage3 = im2bw(blackMaskedImage,normalizedThresholdMuscle);
%% segmentation osteochondroma
% select pixels above threshold (exclude muscle, include osteochondroma)
blackMaskedImage_O = filteredInputImage;
blackMaskedImage_O(~myMuscle) = 0;
if filtreOsteo == 1
segmented_O = im2bw(blackMaskedImage_O,normalizedThresholdOsteochondroma);
myOsteo = medfilt2(segmented_O); % Remove isolated points
myOsteo = imfill(myOsteo, 'holes'); % Fill holes
% Control
tempOsteo = myOsteo; small = 0;
% for i=1 : 2
tempOsteo = imerode(tempOsteo, smallSE); % Erosion - Refine selection
%tempOsteo = imerode(tempOsteo, bigSE);
% end
if max(max(tempOsteo)) == 0
tempOsteo = myOsteo;
tempOsteo = imerode(tempOsteo, smallSE);
%tempOsteo = imerode(tempOsteo, bigSE);
small = 1;
if max(max(tempOsteo)) == 0
tempOsteo = myOsteo;
tempOsteo = imerode(tempOsteo, smallSE);
small = 2;
end
end
myOsteo = tempOsteo;
island = sum(sum(myOsteo)) / 20;
myOsteo = bwareaopen(myOsteo, double(uint16(island))); % Remove islands smaller than 1/10
for i=1 : 1
if small == 2
myOsteo = imdilate(myOsteo, smallSE);
break
end
myOsteo = imdilate(myOsteo, bigSE); % Dilatation/Grow image
if small == 1 break
end
end
myOsteo = bwmorph(myOsteo, 'majority'); % Smooth contours, remove small peaks ? Smoothing
myOsteo = imerode(myOsteo, smallSE);
myOsteo = imdilate(myOsteo, smallSE);
myOsteo = imfill(myOsteo, 'holes');
myOsteo(~myMuscle) = 0;
else
myOsteo = im2bw(blackMaskedImage,normalizedThresholdOsteochondroma);
end
%% 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,d] = 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
if strcmp(output,'OriginalImage');
imshow(A.inputImage,visualizationWindow);
[r,c,d] = size(A.inputImage);
title(sprintf('Subject %s',subject), 'FontSize', fontSize);
set(gcf, 'Position', get(0,'Screensize'), 'Visible', 'off'); % Maximize figure.
oPath = ':\Users\Raph\Desktop\TestContours';
oName = sprintf('Picture_%s',subject);
saveas(gcf,sprintf('%s%s.png',oPath,oName))
end
end
Ra = num2cell(Ra);
Ri = num2cell(Ri);
Ro = num2cell(Ro);
Rd = num2cell(Rd);
end
c4science · Help