Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (master) Browse

R8218/ShoulderCase/@Muscle/measureDegenerationFromSegmentationSet.m
ed6300bff59fmaster

Last Change
Raw File

measureDegenerationFromSegmentationSet.m
View Options

function measureDegenerationFromSegmentationSet(obj,segmentationSetName)
measurementsResults = struct('area',[],...
'degeneration',[],...
'PCSA',[]);
SCase = obj.shoulder.SCase;
segmentationSetPNGFiles = dir(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,'*.png'));
segmentedAreasName = extractBefore({segmentationSetPNGFiles.name},'.png');
load(fullfile(SCase.dataCTPath,'muscles','imagesPixelSpacings.mat'));
obj.slicesPixelSpacings = imagesPixelSpacings;
for i = 1:length(segmentedAreasName)
segmentedImage = getSegmentedImage(obj,segmentationSetName,segmentedAreasName{i});
try
areas = getThresholdedAreas(segmentedImage);
catch ME
disp(ME.message);
continue;
end
degeneration = getDegenerationFromAreas(areas);
PCSA = getPCSAFromAreas(obj,areas);
measurementsResults.area = segmentedAreasName{i};
measurementsResults.degeneration = [measurementsResults.degeneration degeneration];
measurementsResults.PCSA = [measurementsResults.PCSA PCSA];
obj.(segmentedAreasName{i}) = degeneration;
end
measurementsResults.area = struct2table(measurementsResults);
save(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,'measurementsResults.mat'),'measurementsResults');
end
function output = getSegmentedImage(obj,segmentationSetName,muscleName)
SCase = obj.shoulder.SCase;
load(fullfile(SCase.dataCTPath,'muscles','imageForDegeneration.mat'));
normalizedImage = getNormalizedImage(imageForDegeneration);
mask = logical(imread(fullfile(SCase.dataCTPath,'muscles','segmentations',segmentationSetName,[muscleName '.png'])));
segmentedImage = normalizedImage.*mask;
output = segmentedImage;
end
function output = getNormalizedImage(image)
% Normalize image and thresholds
image(image < -1000) = -1000;
image(image > 1000) = 1000;
rangeHU = double([-1000 1000]);
output = (double(image)-rangeHU(1))/(rangeHU(2)-rangeHU(1));
end
function output = getThresholdedAreas(image)
thresholds = getThresholds();
% muscle
muscleArea = imbinarize(image,thresholds.muscle);
muscleArea = imfill(muscleArea,'holes');
% Keep biggest island only
islands = bwconncomp(muscleArea);
islandsSizes = cellfun(@numel,islands.PixelIdxList);
[~,biggestIslandIndex] = max(islandsSizes);
muscleArea = false(size(muscleArea));
muscleArea(islands.PixelIdxList{biggestIslandIndex}) = true;
muscleImage = image.*muscleArea;
% fat
fatArea = muscleArea & not(imbinarize(muscleImage,thresholds.fat));
% osteochondroma
osteochondromaArea = imbinarize(muscleImage,thresholds.osteochondroma);
osteochondromaArea = imfill(osteochondromaArea,'holes');
% healthy muscle
healthyMuscleArea = muscleArea & not(fatArea) & not(osteochondromaArea);
output.segmented = logical(image);
output.muscle = muscleArea;
output.healthyMuscle = healthyMuscleArea;
output.fat = fatArea;
output.osteochondroma = osteochondromaArea;
end
function output = getThresholds()
% 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
rangeHU = double([-1000 1000]);
HUThresholdMuscle = 0;
HUThresholdFat = 30;
HUThresholdOsteochondroma = 166;
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));
output.muscle = normalizedThresholdMuscle;
output.fat = normalizedThresholdFat;
output.osteochondroma = normalizedThresholdOsteochondroma;
end
function output = getDegenerationFromAreas(areas)
degeneration = bwarea(areas.segmented & not(areas.healthyMuscle)) / bwarea(areas.segmented);
output = degeneration;
end
function output = getPCSAFromAreas(obj,areas)
pixelSurface = (obj.slicesPixelSpacings(1) * obj.slicesPixelSpacings(2));
PCSA = bwarea(areas.healthyMuscle) * pixelSurface;
output = PCSA;
end
c4science · Help