Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (pezhman_excelExport) Browse

R8218/ShoulderCase/@Muscles/Muscles.m
ed6300bff59fpezhman_excelExport

Last Change
Raw File

Muscles.m
View Options

classdef Muscles < handle
%MUSCLES This class defines properties and methods of the 4 rotator
%cuff muscles
% Detailed explanation goes here
% Author: Nathan Donini
% Creation date: 2019-02-08
% Revisor: Mathtieu Boubat
% Revision date: 2019-11-22
%
% TO DO:
% - multiple slices and segmentations to compute the average degeneration
properties
SS
SC
IS
TM
slicesPixelSpacings
SS_man % Degeneration ratio of the supraspinatus
SC_man % Degeneration ratio of the subscapularis
IS_man % Degeneration ratio of the infraspinatus
TM_man % Degeneration ratio of the teres minor
SS_auto % Degeneration ratio of the supraspinatus using automatic segmentation
SC_auto % Degeneration ratio of the subscapularis using automatic segmentation
IS_auto % Degeneration ratio of the infraspinatus using automatic segmentation
TM_auto % Degeneration ratio of the teres minors using automatic segmentation
end
properties (Hidden = true)
dataPath
shoulder
end
methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Muscles
function obj = Muscles(shoulder)
%MUSCLES Construct an instance of this class
% Detailed explanation goes here
obj.SS = [];
obj.SC = [];
obj.IS = [];
obj.TM = [];
obj.slicesPixelSpacings = [];
obj.shoulder = shoulder;
obj.SS_man = [];
obj.SC_man = [];
obj.IS_man = [];
obj.TM_man = [];
obj.SS_auto = [];
obj.SC_auto = [];
obj.IS_auto = [];
obj.TM_auto = [];
end
end
methods
function propagateDataPath(obj)
obj.dataPath = fullfile(obj.shoulder.dataPath,'muscles');
if not(isfolder(obj.dataPath))
mkdir(obj.dataPath);
end
end
function output = degeneration(obj,logFileID) % obj is the output of the previous function
%METHOD1
%One needs to create two empty files in the current directory :
%data_degenerescence.mat (a table) and Didnotwork.mat (a
%string)
% This script calculates the degeneration of the 4 muscles
% of the rotator cuff : SS, SC, IS, TM. It is composed of four
% different steps :
% - Load CT images (DICOM files) in order to create a 3D matrix
% consituted with HU (Hounsfield unit) values
% - Create a slice of this 3D matrix along a patient-specific
% plane. Do image processing and set saving format usable by
% automatic segmentation code (Python's code from ARTORG team
% (UNIBE))
% - Execute automatic segmentation code on the new created
% images
% - Apply degeneration code (from Laboratoire de biom�canique
% orthop�dique (EPFL)) on new images with automatic
% segmentation and on image with manual segmentation (from CHUV
% radiologist). Save degeneration values under the "muscles"
% properties of each case in the LBO database
% Some "try-catch" are present in order to catch the errors
% while keeping running the script on the rest of the database.
% They are summarized underneath :
% (DicomReadVolume) : do not manage to load the dicoms
% Nothing : Automatic segmentation does not work (generally for TM)
% * : slice function does not work
% ** : does not want to crop image
% *** : black image, don't want to save png
% Main code starts here !
%% Load dicom
output = {};
SCase = obj.shoulder.SCase;
dataCTPath = SCase.dataCTPath;
% Create "muscles" folder (and its 4 subfolders) if it does not exist
folder_muscle=[dataCTPath '/muscles'];
if isfolder(folder_muscle)
else
mkdir(dataCTPath,'muscles');
mkdir(folder_muscle,'IS');
mkdir(folder_muscle,'SS');
mkdir(folder_muscle,'SC');
mkdir(folder_muscle,'TM');
end
% Use folder "CT-2" (smoothed kernel) if it exists
dataCTPath2=dataCTPath;
dataCTPath2(1,end)='2';
if isfolder(dataCTPath2)
dataCTPath0=dataCTPath2;
else
dataCTPath0=dataCTPath;
end
dicomFolder = [dataCTPath0 '/dicom'];
%Extract the images as a 4D matrix where spatial is the
%positions of the Pixels in the CT coordinate system
% Return the case's number if it does not manage to load the DICOM files
try
[V,spatial,dim] = dicomreadVolume(dicomFolder);
catch
outputArg=1;
fprintf(logFileID, '\n| ERROR: dicomreadVolume()');
return;
end
%V needs to be 3D matrix containing the value in HU of the
%corresponding voxel
V = squeeze(V) ;
%Read the information from the middle dicom image (all image
%informations are the same for all dicom files) to avoid a problem
%if other files are located at the beginning or end of the
%folder
dicomList = dir(dicomFolder) ;
size_list = size(dicomList) ;
dicom_information = dicominfo([dicomFolder '/' dicomList(round(size_list(1)/2)).name]);
%Take into account that the data will need to be linearly transformed
%from stored space to memory to come back to true HU units afterwards:
Rescale_slope = dicom_information.RescaleSlope;
Rescale_intercept = dicom_information.RescaleIntercept;
%Convert to double type to be able to multiply it with a int64
%matrix afterwards
V = double(V);
% To filter hard kernel images
% if isfolder(dataCTPath2)
% else
% V = imgaussfilt3(V,1) ;
% end
size_V = size(V);
%Extract position of each first pixel of each image (upper left corner)
PatientPositions = spatial.PatientPositions ;
%Get the distance between each pixel in neighboring rows and
%columns
PixelSpacings = spatial.PixelSpacings;
% Get local coordinate system from parent scapula object(glenoid center)
origin = obj.shoulder.scapula.coordSys.origin; % do not mix up with origin_image, "origin" is the origin of the Scapula Coordinate System
xAxis = obj.shoulder.scapula.coordSys.PA;
yAxis = obj.shoulder.scapula.coordSys.IS;
zAxis = obj.shoulder.scapula.coordSys.ML;
%% Slice from the CT data set
%Get the position of the upper left pixel of the first image
%and calculate the position of each pixel in the x,y
%and z direction starting from the "origin" pixel
origin_image = PatientPositions(1,:); % coordinates of the first pixel
x_max = origin_image(1)+PixelSpacings(1,1)*(size_V(1)-1);
y_max = origin_image(2)+PixelSpacings(1,2)*(size_V(2)-1);
z_max = PatientPositions(end,3);
%Calculate the coefficient of the linear transformation (CT-scan
%coordinate system (x,y,z) to images coordinates system (i,j,k))
coefficients_i = polyfit([1 size_V(1)],[origin_image(1) , x_max],1);
a_i = coefficients_i(1);
b_i = coefficients_i(2);
coefficients_j = polyfit([1 size_V(2)],[origin_image(2) , y_max],1);
a_j = coefficients_j(1);
b_j = coefficients_j(2);
coefficients_k = polyfit([1 size_V(3)],[origin_image(3) , z_max],1);
a_k = coefficients_k(1);
b_k = coefficients_k(2);
Pixel_pos_x=(1:size_V(1))*a_i+b_i; % Create the position vector along x-axis
Pixel_pos_y=(1:size_V(2))*a_j+b_j;
Pixel_pos_z=(1:size_V(3))*a_k+b_k;
PixSize_z=(z_max-origin_image(3))/(size_V(3)-1); % Space between 2 slices
% The oblique slice to be created is perpendicular to the z
% axis of the scapula CS and passes through the origin point.
% We will apply the matlab function "slice", to create the
% oblique slice. One must first create a slice parallel to 2
% system's axis. Here we take a plane parallel to the (x,y)
% plane, passing through the point "origin". Then, we need to
% rotate this plane so that it becomes normal to the zAxis. V1
% is the vector normal to the horizontal slice and v2 the
% vector normal to the final oblique slice, parallel to the
% zAxis. We need to rotate around a vector v3 and from an angle
% given by "vrrotvec(v1,v2)".
v1=[0,0,1];
v2=zAxis/norm(zAxis);
rotation_info=vrrotvec(v1,v2);
v3=rotation_info(1:3);
theta=rotation_info(4)/pi*180; % Convert it in degrees
% Create the plane to be rotated
[x,y,z] = meshgrid(Pixel_pos_x,Pixel_pos_y,Pixel_pos_z);
s=1000; % Extension factor of the slice in order to cut the whole matrix.Set arbitrarily but sufficiently big so that it works for every patient
Pixel_pos_x_enlarged=[origin_image(1)+(-PixelSpacings(1,1)*s:PixelSpacings(1,1):-PixelSpacings(1,1)) Pixel_pos_x x_max+(PixelSpacings(1,1):PixelSpacings(1,1):PixelSpacings(1,1)*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
Pixel_pos_y_enlarged=[origin_image(2)+(-PixelSpacings(1,2)*s:PixelSpacings(1,2):-PixelSpacings(1,2)) Pixel_pos_y y_max+(PixelSpacings(1,2):PixelSpacings(1,2):PixelSpacings(1,2)*s)];
figure('visible','off');
hsp = surf(Pixel_pos_x_enlarged,Pixel_pos_y_enlarged,origin(3)*ones(length(Pixel_pos_y_enlarged),length((Pixel_pos_x_enlarged))));
rotate(hsp,v3,theta,origin);
xd_0 = hsp.XData; % x-coordinates of the plane's points
yd_0 = hsp.YData; % y-coordinates of the plane's points
zd_0 = hsp.ZData; % z-coordinates of the plane's points
% Slice the 3D matrix of CT-scan
try % Return case's number in case slice does not work
figure('visible','off');
hss=slice(x,y,z,V,xd_0,yd_0,zd_0,'Linear');
% Linear interpolation by default. 'Cubic' or 'Nearest'
catch
outputArg=2;
fprintf(logFileID,'\n| ERROR: slice()');
return
end
HU_0 = get(hss,'CData'); % HU values at each point of the plane
NaNplaces_0=find(isnan(HU_0));
HU_1=HU_0*Rescale_slope+Rescale_intercept; % getting true HU values
p_0=10e6; % value outside 3D matrix, set to a high distinguishable number
HU_1(NaNplaces_0)=p_0;
% If it's a right shoulder or left one, has to be turned in a
% particular way
if obj.shoulder.side=='R'
HU_1=flipud(HU_1');
else
HU_1=HU_1';
end
% Find the new pixel's sizes necessary to obtain a 1024*1024
% shoulder slice. The largest diagonal's length is set to 1024
% pixels
[left_1,left_2]=find(HU_1~=p_0,1);
left_r=left_1;
left_c=left_2;
[right_1,right_2]=find(rot90(HU_1,2)~=p_0,1);
right_r=right_1;
right_c=size(HU_1,2)-right_2+1;
[up_1,up_2]=find(rot90(HU_1)~=p_0,1);
up_c=size(HU_1)-up_1+1;
up_r=up_2;
[down_1,down_2]=find(rot90(HU_1,3)~=p_0,1);
down_r=size(HU_1,1)-down_2+1;
down_c=down_1;
% In case the slice does not cut the whole matrix (which
% normally never happens), ask to double "s" parameter
if or(or(right_c==size(HU_1,2),left_c==1),or(up_r==1,down_r==size(HU_1,1)))
fprintf(logFileID,'\n| ERROR: the slice doesn''t cut throughout the 3D CT-matrix');
outputArg = 3;
return
end
% Find which diagonal of the matrix' effective slice is the largest
w=right_c-left_c; % width
h=down_r-up_r; % height
if h>w % Set the largest diagonal at 1024 pixels
l=h;
else
l=w;
end
% Define new pixel's size
NewPixSize_x=PixelSpacings(1,1)*l/1023;
NewPixSize_y=PixelSpacings(1,2)*l/1023;
% New plane with the new pixel's sizes
Pixel_pos_x_enlarged_2=[origin_image(1)+(-NewPixSize_x*s:NewPixSize_x:-NewPixSize_x) origin_image(1):NewPixSize_x:x_max x_max+(NewPixSize_x:NewPixSize_x:NewPixSize_x*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
Pixel_pos_y_enlarged_2=[origin_image(2)+(-NewPixSize_y*s:NewPixSize_y:-NewPixSize_y) origin_image(2):NewPixSize_y:y_max y_max+(NewPixSize_y:NewPixSize_y:NewPixSize_y*s)];
figure('visible','off');
hsp2 = surf(Pixel_pos_x_enlarged_2,Pixel_pos_y_enlarged_2,origin(3)*ones(length(Pixel_pos_y_enlarged_2),length((Pixel_pos_x_enlarged_2))));
rotate(hsp2,v3,theta,origin);
xd = hsp2.XData; % x-coordinates of the rotated plane's points.
yd = hsp2.YData; % y-coordinates of the rotated plane's points.
zd = hsp2.ZData; % z-coordinates of the rotated plane's points.
% New slice with the good resolution (1024*1024)
figure('visible','off');
hss=slice(x,y,z,V,xd,yd,zd,'Linear');
% Linear interpolation by default. 'Cubic' or 'Nearest'
HU_2 = get(hss,'CData'); % HU values at each point of the plane
NaNplaces=find(isnan(HU_2));
HU_slice=HU_2*Rescale_slope+Rescale_intercept;
HU_slice_8_bit=HU_slice;
p_tif=0; % value outside 3D matrix, for the tif image
p_png=85; % for the png image
HU_slice(NaNplaces)=p_tif;
% Scaling to fit the training set (for the png only, grayscale)
HU_interval_png=[-100 160];
lin_coef=[HU_interval_png(1) 1 ; HU_interval_png(2) 1]^(-1)*[0 255]';
HU_slice_8_bit=lin_coef(1)*HU_slice_8_bit+lin_coef(2);
HU_slice_8_bit(NaNplaces)=p_png;
% If it a right shoulder or left one, has to be turned in a
% particular way
if obj.shoulder.side=='R'
HU_slice=flipud(HU_slice');
HU_slice_8_bit=flipud(HU_slice_8_bit');
else
HU_slice=HU_slice';
HU_slice_8_bit=HU_slice_8_bit';
end
% Determine the location of the 4 corners of the slice
[left_1,left_2]=find(HU_slice~=p_tif,1);
left_r=left_1;
left_c=left_2;
[right_1,right_2]=find(rot90(HU_slice,2)~=p_tif,1);
right_r=right_1;
right_c=size(HU_slice,2)-right_2+1;
[up_1,up_2]=find(rot90(HU_slice)~=p_tif,1);
up_c=size(HU_slice,2)-up_1+1;
up_r=up_2;
[down_1,down_2]=find(rot90(HU_slice,3)~=p_tif,1);
down_r=size(HU_slice,1)-down_2+1;
down_c=down_1;
Folder=[dataCTPath '/matlab'];
% Crop the image
try % If cropping does not work, return case's number
HU_slice=HU_slice(up_r:up_r+1023,left_c:left_c+1023);
HU_slice_8_bit=HU_slice_8_bit(up_r:up_r+1023,left_c:left_c+1023);
HU_slice=int16(HU_slice); % convert into 16-bit
catch
fprintf(logFileID,'\n| ERROR: crop');
outputArg=4;
return
end
% Create '.tif' file
file_name_tif=[SCase.id '.tif'];
t = Tiff(fullfile(Folder,file_name_tif),'w');
setTag(t,'Photometric',Tiff.Photometric.MinIsBlack);
setTag(t,'Compression',Tiff.Compression.None);
setTag(t,'BitsPerSample',16);
setTag(t,'SamplesPerPixel',1);
setTag(t,'ImageLength',1024);
setTag(t,'ImageWidth',1024);
setTag(t,'PlanarConfiguration',Tiff.PlanarConfiguration.Chunky);
setTag(t,'SampleFormat',Tiff.SampleFormat.Int);
write(t,HU_slice);
close(t);
% Create 8 bit '.png' file and save it at the right places
% (matlab folder and python's input folder)
file_name_png=[SCase.id '.png'];
HU_slice_8_bit=uint8(HU_slice_8_bit);
figure('visible','off');
image(HU_slice_8_bit);
colormap(gray(256));
try % return case's number if does not manage to write png in the "matlab" folder because it is a black image
imwrite(HU_slice_8_bit,[Folder '/' file_name_png]);
catch
fprintf(logFileID,'\n| ERROR: black image');
outputArg=5;
return
end
pythonDir = ConfigFileExtractor.getVariable('pythonDir');
Folder_input_py= [pythonDir '/input'];
imwrite(HU_slice_8_bit,[Folder_input_py '/' file_name_png]); % write png file in the python input folder
% imwrite(HU_slice_8_bit,['./rcseg/' file_name_png]); % write png file in local /rceg folder
% Create .txt file with pixel dimensions
file_name_info=[SCase.id '.info'];
text1='# Amira Stacked Slices' ;
text2=['pixelsize ', num2str(NewPixSize_x),' ',num2str(NewPixSize_y)];
fileID = fopen([Folder '/' file_name_info],'w');
fprintf(fileID,'%s\r\n\r\n',text1);
fprintf(fileID,'%s\n',text2);
fclose(fileID);
close all;
%% Apply UNIBE method to automatically segment normal muscle
pythonCmd1 = ['source ' pythonDir '/venv/bin/activate'];
pythonCmd2 = ['cd ' pythonDir];
pythonCmd3 = './rcseg.py segment';
pythonCmd3Arg = 'input';
pythonCmd = [pythonCmd1 ';' pythonCmd2 ';' pythonCmd3 ' ' pythonCmd3Arg];
[status, pythonCmdOut] = system(pythonCmd);
disp(pythonCmdOut);
% Transform output into logical and group with initial
% 16 bits file (HU_slice)
caseID=SCase.id;
muscleDirectory=[dataCTPath '/muscles'];
musclesList=char('IS','SC','SS','TM');
inputImage=HU_slice;
inputFile = ['ssContours' caseID '.mat'];
inputFile_auto=['ssContours' caseID '_auto.mat'];
% Save matlab file at the right place
for n=1:length(musclesList)
muscleName_0=musclesList(n,:);
Folder_pyth=[pythonDir '/' muscleName_0];
binaryImage1_8_bit=imread(fullfile(Folder_pyth,file_name_png));
binaryImage1=logical(binaryImage1_8_bit);
output{end+1} = binaryImage1;
Folder_muscles=[dataCTPath '/muscles/' muscleName_0];
complete_name=[Folder_muscles '/' inputFile_auto];
save(complete_name,'binaryImage1','inputImage');
end
%% Apply LBO methods to estimate degeneration from image and normal muscle segmentation
HUThresholdMuscle = 0; %#ok<NASGU>
HUThresholdFat = 30; %#ok<NASGU>
HUThresholdOsteochondroma = 166; %#ok<NASGU>
% Check if there is a manual segmentation
if isfile([muscleDirectory '/IS/' inputFile]) % Suppose if the file is in 'IS' folder, then it is present in the other muscles' folders
k=2;
else
k=1;
end
degeneration=zeros(k,4);
% Calculate the degeneration with automatic segmentation and
% with manual segmentation when it is present
for j=1:length(musclesList)
muscleName=musclesList(j,:);
for m=1:k
if m==1
inputFile_auto_man=inputFile_auto;
else
inputFile_auto_man=inputFile;
end
contouredMuscle = load([muscleDirectory '/' muscleName '/' inputFile_auto_man]);
% --Atrophied muscle segmentation--
% 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
% Normalize image and thresholds
contouredMuscle.inputImage(contouredMuscle.inputImage < -1000) = -1000;
contouredMuscle.inputImage(contouredMuscle.inputImage > 1000) = 1000;
rangeHU = double([-1000 1000]);
% Binarization of images works so that pixel value > threshold is kept. So
% if you condier that a muscle is normal for 30 HU and abnormal for 29 HU,
% the threshold should be set at 29 HU
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));
normalizedInputImage = (double(contouredMuscle.inputImage)-rangeHU(1))/(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
% Keep only the biggest connected component (i.e. remove islands)
CC = bwconncomp(myMuscle);
numPixels = cellfun(@numel,CC.PixelIdxList);
[~,idx] = max(numPixels);
myMuscle(:,:) = 0;
try % In case automatic segmentation does not work (has returned a black image for TM, most of the time)
myMuscle(CC.PixelIdxList{idx}) = 1;
catch
fprintf(logFileID,'\n| ERROR: automatic segmentation');
outputArg=6;
input_to_delete=[pythonDir '/input/' file_name_png];
delete(input_to_delete);
return
end
% --Fat infiltration segmentation--
% Mask image using atrophied muscle segmentation
maskedImage = normalizedInputImage;
maskedImage(~myMuscle) = 0;
% Segment using fat threshold
myMuscleNoFat = imbinarize(maskedImage,normalizedThresholdFat);
myFat = ~myMuscleNoFat;
myFat(~myMuscle) = 0;
% --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--
muscle.caseID = caseID;
muscle.Stot = bwarea(contouredMuscle.binaryImage1); % Stot: total area of the muscle fossa (manually delimited contours)
muscle.Satrophy = bwarea(contouredMuscle.binaryImage1) - bwarea(myMuscle); % Sa: area of the atrophy (Stot - Sm)
muscle.Sinfiltration = bwarea(myFat); % Si: area of fat infiltration
muscle.Sosteochondroma = bwarea(myOsteo); % So: area of osteochondroma
muscle.Sdegeneration = muscle.Stot-muscle.Satrophy-muscle.Sinfiltration-muscle.Sosteochondroma; % Sm: area of the degenerated muscle, without atrophy, fat infiltration and osteochondroma
muscle.atrophy = muscle.Satrophy/muscle.Stot; % ratio of muscle atrophy (area atrophy over total area of muscle fossa)
muscle.infiltration = muscle.Sinfiltration/muscle.Stot; % ratio of fat infiltration in the muscle (area fat infiltration over total area of muscle fossa)
muscle.osteochondroma = muscle.Sosteochondroma/muscle.Stot; % ratio of osteochondroma in the muscle (area osteochondroma over total area of muscle fossa)
muscle.degeneration = (muscle.Stot-muscle.Sdegeneration)/muscle.Stot; % the ratio of degeneration of the muscle (sum of atrophy, fat infiltration and osteochondroma over total area of muscle fossa)
degeneration(m,j)=muscle.degeneration;
atrophy(m,j)=muscle.atrophy;
infiltration(m,j)=muscle.infiltration;
osteochondroma(m,j)=muscle.osteochondroma;
end
end
% 4 or 8 outputs, depending if there exist a manual segmentation
obj.SS_auto = degeneration(1,1);
obj.SC_auto = degeneration(1,2);
obj.IS_auto = degeneration(1,3);
obj.TM_auto = degeneration(1,4);
if k==2 % if there exist a manual segmentation
obj.SS_man = degeneration(2,1);
obj.SC_man = degeneration(2,2);
obj.IS_man = degeneration(2,3);
obj.TM_man = degeneration(2,4);
end
close all
% Cleaning of the SCase matlab folder
delete([Folder '/' file_name_png]);
delete([Folder '/' file_name_tif]);
delete([Folder '/' file_name_info]);
% Avoid calculating several times the same case's
% degeneration values
input_to_delete=[pythonDir '/input/' file_name_png];
delete(input_to_delete);
outputArg = 0;
end
end
end
c4science · Help