diff --git a/muscleDegeneration.m b/muscleDegeneration.m index a429b73..ea114ca 100644 --- a/muscleDegeneration.m +++ b/muscleDegeneration.m @@ -1,190 +1,170 @@ 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 +% Load the reference 2D CT slice and associated muscle contour 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 +contouredMuscle = load(sprintf('%s%s',inputPath,inputFile)); -myMuscle(~A.binaryImage1) = 0; +% --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)); -% 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); +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 -% segmentation osteochondroma -% select pixels above threshold (exclude muscle, include osteochondroma) -myOsteo = imbinarize(blackMaskedImage,normalizedThresholdOsteochondroma); -myOsteo = imfill(myOsteo, 'holes'); -myOsteo(~myMuscle) = 0; +% --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); -% 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; +% --Osteochondroma segmentation-- -%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; +%Segment using osteochondroma threshold +myOsteo = imbinarize(maskedImage,normalizedThresholdOsteochondroma); +myOsteo = imfill(myOsteo, 'holes'); % Fill holes -% % Calculate infiltration -Ri = Si/Stot; -% % Calculate osteochondroma -Ro = So/Stot; -% % Calculate degeneration -Rd = (Stot-Sm)/Stot; +% --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') - imshow(A.inputImage,visualizationWindow); - [r,c,~] = size(A.inputImage); + + % --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 - % 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. + % 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(x3, y3,'color','red'); + plot(x, y,'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. + % 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(x3, y3,'color','magenta'); + 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); - 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 \ No newline at end of file