Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (pezhman_excelExport) Browse

R8218/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
02f485383659pezhman_excelExport

Last Change
Raw File

createRotatorCuffSlicesNathan.m
View Options

function createRotatorCuffSlicesNathan(obj,sliceName)
% This function and its results might be reffered as
% Nathan algorithm and Nathan slices in some studies.
%
%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
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;
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;
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)))
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
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
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;
imwrite(HU_slice_8_bit,fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']));
imageForMeasurements = HU_slice;
save(fullfile(obj.slicesDataPath,[sliceName '_ForMeasurements.mat']),'imageForMeasurements');
imagesPixelSpacings = [NewPixSize_x NewPixSize_y];
save(fullfile(obj.slicesDataPath,[sliceName '_PixelSpacings.mat']),'imagesPixelSpacings');
end
c4science · Help