diff --git a/ShoulderCase/Muscles.m b/ShoulderCase/@Muscles/Muscles.m similarity index 99% rename from ShoulderCase/Muscles.m rename to ShoulderCase/@Muscles/Muscles.m index a047dea..8ef6838 100644 --- a/ShoulderCase/Muscles.m +++ b/ShoulderCase/@Muscles/Muscles.m @@ -1,597 +1,599 @@ 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_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) 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_man = []; obj.SC_man = []; obj.IS_man = []; obj.TM_man = []; obj.SS_auto = []; obj.SC_auto = []; obj.IS_auto = []; obj.TM_auto = []; obj.shoulder = shoulder; end end methods - function outputArg = degeneration(obj,logFileID) % obj is the output of the previous function + 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 HUThresholdFat = 30; %#ok HUThresholdOsteochondroma = 166; %#ok % 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