diff --git a/ShoulderCase/@DicomVolume/DicomVolume.m b/ShoulderCase/@DicomVolume/DicomVolume.m new file mode 100644 index 0000000..b4d1bf8 --- /dev/null +++ b/ShoulderCase/@DicomVolume/DicomVolume.m @@ -0,0 +1,55 @@ +classdef (Abstract) DicomVolume < handle + + properties + volume = [] + spatial = [] + dicomInfo = [] + dicomFolderPath = [] + end + + + + methods + function set.volume(obj,volume) + assert(length(size(volume))>=3,'Volume must be a 3D array') + obj.volume = squeeze(volume); + end + + function set.spatial(obj,spatial) + requiredStructureField = {'PatientPositions','PixelSpacings','PatientOrientations'}; + assert(all(contains(fields(spatial),requiredStructureField)),'Invalid spatial structure'); + obj.spatial = spatial; + end + + function loadDataFromFolder(obj,dicomFolderPath) + assert(obj.pathContainsDicomFiles(dicomFolderPath)); + obj.dicomFolderPath = dicomFolderPath; + [obj.volume,obj.spatial,~] = dicomreadVolume(dicomFolderPath); + obj.loadDicomInfoFromFolder(dicomFolderPath); + end + + function loadDicomInfoFromFolder(obj,dicomFolderPath) + assert(obj.pathContainsDicomFiles(dicomFolderPath)); + dicomFiles = dir(fullfile(dicomFolderPath,'*.dcm')); + obj.dicomInfo = dicominfo(fullfile(dicomFiles(1).folder,dicomFiles(1).name)); + end + + function output = pathContainsDicomFiles(obj,path) + output = false; + assert(isfolder(path),'Provided argument is not a valid folder path'); + dicomFiles = dir(fullfile(path,'*.dcm')); + assert(not(isempty(dicomFiles)),'No dicom file found there %s',path); + output = true; + end + + function output = getPointIndexInVolume(obj,point) + % the output are the three indices of the given point in the given volume. + % i,j,k: output indices + k = find(abs(obj.spatial.PatientPositions(:,3)-point(3)) == min(abs(obj.spatial.PatientPositions(:,3)-point(3)))); + i = round( (point(1) - obj.spatial.PatientPositions(k,1)) / obj.spatial.PixelSpacings(k,1) ); + j = round( (point(2) - obj.spatial.PatientPositions(k,2)) / obj.spatial.PixelSpacings(k,2) ); + assert(all([i j k] < size(obj.volume)),'Indices found are outside the volume boundaries'); + output = [i,j,k]; + end + end +end diff --git a/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m new file mode 100644 index 0000000..906df6c --- /dev/null +++ b/ShoulderCase/@DicomVolumeNormaliser/DicomVolumeNormaliser.m @@ -0,0 +1,55 @@ +classdef DicomVolumeNormaliser < handle & DicomVolume + + properties + resolution = []; + normalisedVolume = []; + end + + + + methods + + function obj = DicomVolumeNormaliser(varargin) + if (nargin == 1) + obj.loadDataFromFolder(varargin{1}); + end + end + + function setMinimalResolution(obj) + zSpacing = obj.spatial.PatientPositions(2:end,3)-obj.spatial.PatientPositions(1:end-1,3); + minZ = min(zSpacing,[],'all'); + minXY = min(obj.spatial.PixelSpacings,[],'all'); + obj.resolution = min([minXY minZ],[],'all')/2; + end + + function normalise(obj) + normalisedSizeZ = round(abs(obj.spatial.PatientPositions(1,3) - obj.spatial.PatientPositions(end,3))/obj.resolution); + normalisedSizeX = round(max(obj.spatial.PixelSpacings(:,1),[],'all')*size(obj.volume,1)/obj.resolution); + normalisedSizeY = round(max(obj.spatial.PixelSpacings(:,2),[],'all')*size(obj.volume,2)/obj.resolution); + obj.normalisedVolume = imresize3(obj.volume,[normalisedSizeX normalisedSizeY normalisedSizeZ]); + end + + function OLDnormalise(obj) + normalisedSizeZ = round(abs(obj.spatial.PatientPositions(1,3) - obj.spatial.PatientPositions(end,3))/obj.resolution); + normalisedSizeX = round(max(obj.spatial.PixelSpacings(:,1),[],'all')*size(obj.volume,1)/obj.resolution); + normalisedSizeY = round(max(obj.spatial.PixelSpacings(:,2),[],'all')*size(obj.volume,2)/obj.resolution); + obj.normalisedVolume = zeros(normalisedSizeX,normalisedSizeY,normalisedSizeZ); + for x = 1:size(obj.normalisedVolume,1) + disp(sprintf('%d', round(100*x/size(obj.normalisedVolume,1)))) + for y = 1:size(obj.normalisedVolume,2) + for z = 1:size(obj.normalisedVolume,3) + pointPosition = obj.spatial.PatientPositions(1,:)+obj.resolution*[x y z]; + try + pointIndices = obj.getPointIndexInDicomVolume(pointPosition); + catch + [x y z] + pointPosition + end + obj.normalisedVolume(x,y,z) = obj.volume(pointIndices(1),pointIndices(2),pointIndices(3)); + end + end + end + end + end + +end diff --git a/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m new file mode 100644 index 0000000..bee14b4 --- /dev/null +++ b/ShoulderCase/@DicomVolumeSlicer/DicomVolumeSlicer.m @@ -0,0 +1,30 @@ +classdef DicomVolumeSlicer < handle & DicomVolume & DicomVolumeNormaliser + + properties + sliced + slicedPixelSpacing + end + + + + methods + + function obj = DicomVolumeSlicer(varargin) + if (nargin == 1) + obj.loadDataFromFolder(varargin{1}); + end + end + + function slice(obj,point,normalVector) + pointIndices = obj.getPointIndexInVolume(point); + vectorIndices = pointIndices - obj.getPointIndexInVolume(point+normalVector); + + [obj.sliced,x,y,z] = obliqueslice(obj.volume,pointIndices,vectorIndices,... + 'FillValues',min(obj.volume,[],'all')); + obj.sliced = obj.sliced'; + obj.sliced = obj.sliced*obj.dicomInfo.RescaleSlope + obj.dicomInfo.RescaleIntercept; + end + + end + +end diff --git a/ShoulderCase/@MethodsMonitor/MethodsMonitor.m b/ShoulderCase/@MethodsMonitor/MethodsMonitor.m new file mode 100644 index 0000000..bcec1ae --- /dev/null +++ b/ShoulderCase/@MethodsMonitor/MethodsMonitor.m @@ -0,0 +1,116 @@ +classdef (Abstract) MethodsMonitor < handle + + properties (SetObservable) + methodsMonitor + end + + + methods + function startMethodsMonitoring(obj) + objectMethods = obj.getObjectMethods; + + obj.methodsMonitor = table(); + for i = 1:length(objectMethods) + obj.methodsMonitor = [obj.methodsMonitor, {'';'';''}]; + end + obj.methodsMonitor.Properties.VariableNames = objectMethods; + obj.methodsMonitor.Properties.RowNames = {'state';'runtime';'comment'}; + + addlistener(obj,'methodsMonitor','PostSet',@obj.launchOrderedMethods); + obj.setAllMethodsUnavailable; + end + + function orderMethod(obj,method) + obj.setMethodState(method,'ordered') + end + + function setAllMethodsAvailable(obj) + obj.setAllMethodsState('available'); + end + + function setMethodAvailable(obj,method) + obj.setMethodState(method,'available') + end + + function setAllMethodsUnavailable(obj) + obj.setAllMethodsState('unavailable'); + end + + function setMethodUnavailable(obj,method) + obj.setMethodState(method,'unavailable') + end + + function setMethodComment(obj,method,comment) + obj.setMethodVariableValue(method,'comment',comment); + end + + function output = isMethodAvailable(obj,method) + output = isequal(obj.(method)('state'),{'available'}); + end + end + + + + methods (Access = protected) + function output = getObjectMethods(obj) + objectMethods = methods(obj); + handleMethods = methods('handle'); + MethodsMonitorMethods = methods('MethodsMonitor'); + objectMethods = cellArraysDifference(objectMethods,handleMethods,MethodsMonitorMethods); + + output = objectMethods; + end + + function launchOrderedMethods(obj,~,~) + method = obj.getNextOrderedMethod; + while not(isempty(method)) + obj.runMethod(method); + method = obj.getNextOrderedMethod; + end + end + + function setAllMethodsState(obj,state) + objectMethods = obj.getObjectMethods; + for i = 1:length(objectMethods) + obj.setMethodState(objectMethods{i},state); + end + end + + function setMethodState(obj,method,state) + obj.setMethodVariableValue(method,'state',state); + end + + function setMethodRuntime(obj,method,runtime) + obj.setMethodVariableValue(method,'runtime',runtime); + end + + function setMethodVariableValue(obj,method,variable,value) + obj.methodsMonitor.(method)(variable) = {value}; + end + + function runMethod(obj,method) + Timer.start; + try + eval(['obj.' method ';']); + obj.setMethodState(method,'done'); + obj.setMethodComment(method,''); + catch ME + obj.setMethodState(method,'failed'); + obj.setMethodComment(method,ME.message) + end + obj.setMethodRuntime(method,Timer.stop); + end + + function output = getNextOrderedMethod(obj) + output = []; + objectMethods = obj.methodsMonitor.Properties.VariableNames; + for i = 1:length(objectMethods) + if isequal(obj.methodsMonitor.(objectMethods{i})('state'),{'ordered'}) + output = objectMethods{i}; + return; + end + end + end + end + +end 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 diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m index 4852a0f..9fad905 100644 --- a/ShoulderCase/@Shoulder/Shoulder.m +++ b/ShoulderCase/@Shoulder/Shoulder.m @@ -1,39 +1,221 @@ classdef (Abstract) Shoulder < handle properties side % R or L scapula humerus muscles CTscan comment end properties (Hidden = true) SCase end methods (Abstract, Access = protected) createScapula(obj); end methods function obj = Shoulder(SCase) obj.side = ''; obj.humerus = Humerus(obj); obj.muscles = Muscles(obj); obj.CTscan = ''; obj.comment = ''; obj.SCase = SCase; obj.createScapula; end + + function output = plot(obj) + % Create a figure and plot inside: + % the scapula surface, coordSys, groove points and landmarks; + % the glenoid surface, sphere and centerLine; + % the humerus sphere and centerLine. + + figure; + hold on; + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %% Scapula + % Plot the segmented surface + if not(isempty(obj.scapula.surface)) + points = obj.scapula.surface.points; + faces = obj.scapula.surface.faces; + x = points(:,1); + y = points(:,2); + z = points(:,3); + trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud'); + end + + % Plot scapula groove + if not(isempty(obj.scapula.groove)) + points = obj.scapula.groove; + x = points(:,1); + y = points(:,2); + z = points(:,3); + scatter3(x,y,z,100,'blue'); + end + + % Plot scapula markers + points = [... + obj.scapula.angulusInferior; + obj.scapula.trigonumSpinae; + obj.scapula.processusCoracoideus; + obj.scapula.acromioClavicular; + obj.scapula.angulusAcromialis; + obj.scapula.spinoGlenoidNotch... + ]; + x = points(:,1); + y = points(:,2); + z = points(:,3); + scatter3(x,y,z,100,'blue','LineWidth',2); + + + % Plot scapular plane + plot(obj.scapula.plane,100); + + % Plot scapular axis + pt1 = obj.scapula.coordSys.origin; + axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]); + pt2 = pt1 - obj.scapula.coordSys.ML * axisLength; + x = [pt1(1), pt2(1)]; + y = [pt1(2), pt2(2)]; + z = [pt1(3), pt2(3)]; + plot3(x,y,z,'Color','blue','LineWidth', 5); + + % Plot scapular coordinate system + axisLength = 50; % Length of the coord. sys. axes + pt0 = obj.scapula.coordSys.origin; + pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X + pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y + pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z + x = [pt0(1), pt1(1)]; + y = [pt0(2), pt1(2)]; + z = [pt0(3), pt1(3)]; + plot3(x,y,z,'Color','red','LineWidth', 5); + x = [pt0(1), pt2(1)]; + y = [pt0(2), pt2(2)]; + z = [pt0(3), pt2(3)]; + plot3(x,y,z,'Color','green','LineWidth', 5); + x = [pt0(1), pt3(1)]; + y = [pt0(2), pt3(2)]; + z = [pt0(3), pt3(3)]; + plot3(x,y,z,'Color','blue','LineWidth', 5); + + % plot scapular axis on glenoid center + axisLength = 50; + pt1 = obj.scapula.glenoid.center; + pt2 = pt1 + obj.scapula.coordSys.ML*axisLength; + x = [pt1(1), pt2(1)]; + y = [pt1(2), pt2(2)]; + z = [pt1(3), pt2(3)]; + plot3(x,y,z,'Color','magenta','LineWidth', 5); + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Glenoid + + % Plot the glenoid surface + % We might move them in the direction of th glenoid center + % instead of the ML axis + + points = obj.scapula.glenoid.surface.points +... + obj.scapula.coordSys.ML * 0.2; + faces = obj.scapula.glenoid.surface.faces; + x = points(:,1); + y = points(:,2); + z = points(:,3); + trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none'); + + + % plot glenoid centerline + pt1 = obj.scapula.glenoid.center; + pt2 = pt1 + obj.scapula.glenoid.centerLine; + x = [pt1(1), pt2(1)]; + y = [pt1(2), pt2(2)]; + z = [pt1(3), pt2(3)]; + plot3(x,y,z,'Color','magenta','LineWidth', 5); + + % We might add spherical cap + + % plot glenoid sphere + n = 20; % number of lines + [X,Y,Z] = sphere(n); + X = X*obj.scapula.glenoid.radius; + Y = Y*obj.scapula.glenoid.radius; + Z = Z*obj.scapula.glenoid.radius; + Xc = obj.scapula.glenoid.center(1); + Yc = obj.scapula.glenoid.center(2); + Zc = obj.scapula.glenoid.center(3); + mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud'); + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Humerus + + % plot humeral head sphere + n = 20; % number of lines + [X,Y,Z] = sphere(n); + X = X*obj.humerus.radius; + Y = Y*obj.humerus.radius; + Z = Z*obj.humerus.radius; + Xc = obj.humerus.center(1); + Yc = obj.humerus.center(2); + Zc = obj.humerus.center(3); + mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud'); + + + % plot humeral head centerline + pt1 = obj.scapula.glenoid.center; + pt2 = obj.humerus.center; + x = [pt1(1), pt2(1)]; + y = [pt1(2), pt2(2)]; + z = [pt1(3), pt2(3)]; + plot3(x,y,z,'Color','magenta','LineWidth', 5); + + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Graph properties + + axis off; % remove axis + lightPosition = obj.scapula.glenoid.center + ... + (... + obj.scapula.coordSys.ML + ... + obj.scapula.coordSys.IS + ... + obj.scapula.coordSys.PA ... + ) * 100; + light('Position',lightPosition,'Style','local'); + grid on; + xlabel('X (CT)'); + ylabel('Y (CT)'); + zlabel('Z (CT)'); + axis square; + axis vis3d; % fixed limits and scaling + viewPoint = obj.scapula.coordSys.ML + ... + obj.scapula.coordSys.IS + ... + obj.scapula.coordSys.PA; + viewPoint = obj.scapula.coordSys.PA; + view(viewPoint); % view from lateral side + % Align IS axis with window vertical axis + camup(obj.scapula.coordSys.IS); + + figureHandle = gca; + figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis + figureHandle.Box = 'On'; + set(gcf,'color','w'); % Set background color to white + + output = 1; + end end end diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m index 968d14e..07f325c 100755 --- a/ShoulderCase/@ShoulderCase/ShoulderCase.m +++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m @@ -1,486 +1,298 @@ classdef ShoulderCase < handle % Properties and methods associated to the SoulderCase object. % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2018-12-31 % TODO: % Need method to load properties (from amira, matlab, excel, SQL) % Need method to save properties (in matlab, excel, SQL) % Remove datapath from constants % properties id % id of the shoulder case, as it appears in the database (Pnnn) diagnosis treatment outcome patient shoulderManual shoulderAuto study comment end % properties (Constant, Hidden = true) % dataPath = '/Volumes/shoulder/dataDev'; % path to data folder from package folder % % This should be done differently. Check where we are, and where we % % should go. % end properties (Hidden = true) dataPath dataCTPath dataAmiraPath dataMatlabPath id4C % SCaseId with P/N followed by 3 digits --> 4 char end methods function obj = ShoulderCase(SCaseID) %SHOULDERCASE Construct an instance of this class % If called with argument SCaseId, set id and path of CT rawSCaseID = SCaseIDParser(SCaseID); assert(rawSCaseID.isValidID,'The input argument is not a valid SCaseID.') obj.id = SCaseID; % Initialsation obj.diagnosis = ''; obj.treatment = ''; obj.outcome = ''; obj.patient = Patient(obj); obj.shoulderManual = ShoulderManual(obj); obj.shoulderAuto = ShoulderAuto(obj); obj.study = ''; obj.comment = ''; obj.dataCTPath = ''; obj.dataAmiraPath = ''; obj.dataMatlabPath = ''; obj.id4C = ''; % Set path of the SCase % obj.path; To be set only if varargin ~isempty end function outputArg = path(obj) % Should be private % PATH Summary of this method goes here % Detailed explanation goes here SCase = obj; SCaseId = SCase.id; % The validity of the format should be checked Pnnn or Nnnn. if (numel(regexp(SCaseId,'^[PN]\d{1,3}$')) == 0) error(['Invalid format of SCaseId: ' SCaseId{1i} '. CaseID must start with "P" or "N" and be followed by 1 to 3 digits.']); end % Set the folder of the SCaseId SCaseDirLevelPN = SCaseId(1); % Either 'P' or 'N' strLengthSCaseId = strlength(SCaseId(2:end)); if (strLengthSCaseId < 2) SCaseDirLevel2 = '0'; % Hunderets SCaseDirLevel1 = '0'; % Dozent elseif (strLengthSCaseId < 3) SCaseDirLevel2 = '0'; % Hunderets SCaseDirLevel1 = SCaseId(2); % Dozent else SCaseDirLevel2 = SCaseId(2); % Hunderets SCaseDirLevel1 = SCaseId(3); % Dozent end % Set SCaseId with fixed 3 digits after P/N (needed when id < % 10 inloading of data in amira directory. SCaseId4C = [SCaseDirLevelPN SCaseDirLevel2 SCaseDirLevel1 SCaseId(end)]; obj.id4C = SCaseId4C; % Check if a (!unique! to be done) directory exists for this SCaseId dataDir = obj.dataPath; FindSCaseIdFolder = dir([dataDir '/' SCaseDirLevelPN '/' SCaseDirLevel2 '/' SCaseDirLevel1 '/' SCaseId '*']); if (isempty(FindSCaseIdFolder)) % No directory for this SCaseId error(['Missing directory for SCaseId: ' SCaseId]); end SCaseDirLevel0 = FindSCaseIdFolder.name; SCaseDir = [dataDir '/' SCaseDirLevelPN '/' SCaseDirLevel2 '/' SCaseDirLevel1 '/' SCaseDirLevel0]; %{ %8 lines below to remove when -1 consrain avoided % Check if this SCaseId has a CT directory with '-1' postfix (preoperative % sharp CT FindCTFolder = dir([SCaseDir '/' 'CT*-1']); if (isempty(FindCTFolder)) % No CT directory for this SCaseId error(['Missing directory for SCaseId: ' SCaseId]); % ! should exist script here ! end CTDir = [SCaseDir '/' FindCTFolder.name]; %} CTdirList=dir([SCaseDir '/CT-*']); % List directories with CT iCTdirAmira = 0; iCTdir2use = 0; for iCTdirList = length(CTdirList):-1:1 % Loop from last to first (4 to 1) CTdir = CTdirList(iCTdirList); % Check that dir name ends with -1,-2,-3,-4 dirName = CTdir.name; if strcmp(dirName(end-1),'-') % Exclude postoperative 'p' CT CTnum = str2num(dirName(end)); if CTnum <= 4 % Exlude non shoulder (elbow) CT % Check that the dir contains a dicom dir CTdir = [CTdir.folder '/' CTdir.name]; dicomDir = [CTdir '/dicom']; if exist(dicomDir, 'dir') % Chech if amira dir exist amiraDir = [CTdir '/amira']; if exist(amiraDir, 'dir') % This is the CT folder to use iCTdirAmira = iCTdirList; end iCTdir2use = iCTdirList; end end end end if iCTdir2use == 0 error(['\n', SCaseID, ' has no valid dicom directory']); else if iCTdirAmira % If amira dir exist, use it iCTdir2use = iCTdirAmira; end CTdir = CTdirList(iCTdir2use); CTdirPath = [CTdir.folder '/' CTdir.name]; end CTDir = CTdirPath; amiraDir = [CTDir '/amira']; % we should check for existance matlabDir = [CTDir '/matlab']; % we should check for existance obj.dataCTPath = CTDir; obj.dataAmiraPath = amiraDir; obj.dataMatlabPath = matlabDir; outputArg = 1; % Should report on directories existence if exist(amiraDir, 'dir') ~= 7 outputArg = 0; end end function outputArg = output(obj, varargin) % Below is copy/past from previous version % This function is used to output the variable Results, used % by scapula_measurePKG the modified funcion scapula_measure. % inputArg could be as in scapula_calculation % 'density', 'References', 'obliqueSlice', 'display' Result.SCase_id = obj.id; % 1 Result.glenoid_Radius = obj.shoulder.scapula.glenoid.radius; % 2 Result.glenoid_radiusRMSE = obj.shoulder.scapula.glenoid.sphereRMSE; % 3 % 3 %Result.glenoidSphericity = ' '; %Result.glenoid_biconcave = ' '; Result.glenoid_depth = obj.shoulder.scapula.glenoid.depth; % 4 Result.glenoid_width = obj.shoulder.scapula.glenoid.width; % 5 Result.glenoid_height = obj.shoulder.scapula.glenoid.height; % 6 Result.glenoid_center_PA = obj.shoulder.scapula.glenoid.centerLocal(1); % 7 Result.glenoid_center_IS = obj.shoulder.scapula.glenoid.centerLocal(2); % 8 Result.glenoid_center_ML = obj.shoulder.scapula.glenoid.centerLocal(3); % 9 Result.glenoid_version_ampl = obj.shoulder.scapula.glenoid.versionAmpl; % 10 Result.glenoid_version_orient = obj.shoulder.scapula.glenoid.versionOrient; % 11 Result.glenoid_Version = obj.shoulder.scapula.glenoid.version; % 12 Result.glenoid_Inclination = obj.shoulder.scapula.glenoid.inclination; % 13 Result.humerus_joint_radius = ' '; % 14 Result.humeral_head_radius = obj.shoulder.humerus.radius; % 15 % Result.humerus_GHsublux_2D = ' '; % Result.humerus_SHsublux_2D = ' '; Result.humerus_GHsubluxation_ampl = obj.shoulder.humerus.GHSAmpl; % 16 Result.humerus_GHsubluxation_orient = obj.shoulder.humerus.GHSOrient; % 17 Result.humerus_SHsubluxation_ampl = obj.shoulder.humerus.SHSAmpl; % 18 Result.humerus_SHsubluxation_orient = obj.shoulder.humerus.SHSOrient; % 19 Result.scapula_CSA = obj.shoulder.scapula.acromion.criticalShoulderAngle; % radCSA; % 5 Lines below should be updated Result.scapula_CTangle = 0; %CTorientation; %20 Result.scapula_CTangleVersion = 0; % WearPlaneAngle; %21 Result.scapula_CTangleSHS = 0; % SHSPlaneAngle; % 22 Result.scapula_CTangleGHS = 0; % GHSPlaneAngle; % 23 Result.scapula_PlaneRMSE = 0; %PlaneRMSE;%24 Result.scapula_AI = obj.shoulder.scapula.acromion.acromionIndex; outputArg = Result; end function outputArg = saveMatlab(obj) % Save SCase to matlab file dir = obj.dataMatlabPath; % Create dir if not exist try if ~exist(dir, 'dir') % create directory if it does not exist mkdir(dir); end catch fprintf('Error creating the matlab directory \n'); % Should be in log end % Save SCase in matlab directoty, in a file named SCaseCNNN.m filename = 'SCase'; filename = [dir '/' filename '.mat']; try SCase = obj; save(filename, 'SCase'); outputArg = 1; catch fprintf('Error creating SCase matlab file \n'); % Should be in log outputArg = -1; end end function outputArg = appendToCSV(obj,filename) % Save SCase to csv file logFid = fopen('log/measureSCase.log', 'a'); dataDir = openConfigFile('config.txt', logFid); xlsDir = [dataDir '/Excel/xlsFromMatlab']; fid = fopen([xlsDir '/' filename],'a'); fprintf(fid,[... obj.id ','... % SCase_id obj.shoulder.side ','... % shoulder_side num2str(obj.shoulder.scapula.glenoid.radius) ','... % glenoid_radius num2str(obj.shoulder.scapula.glenoid.sphereRMSE) ','... % glenoid_sphereRMSE num2str(obj.shoulder.scapula.glenoid.depth) ','... % glenoid_depth num2str(obj.shoulder.scapula.glenoid.width) ','... % glenoid_width num2str(obj.shoulder.scapula.glenoid.height) ','... % glenoid_height num2str(obj.shoulder.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA num2str(obj.shoulder.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS num2str(obj.shoulder.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML num2str(obj.shoulder.scapula.glenoid.versionAmpl) ','... % glenoid_versionAmpl num2str(obj.shoulder.scapula.glenoid.versionOrient) ','... % glenoid_versionOrient num2str(obj.shoulder.scapula.glenoid.version) ','... % glenoid_version num2str(obj.shoulder.scapula.glenoid.inclination) ','... % glenoid_inclination num2str(obj.shoulder.humerus.jointRadius) ','... % humerus_jointRadius num2str(obj.shoulder.humerus.radius) ','... % humerus_headRadius num2str(obj.shoulder.humerus.GHSAmpl) ','... % humerus_GHSAmpl num2str(obj.shoulder.humerus.GHSOrient) ','... % humerus_GHSOrient num2str(obj.shoulder.humerus.SHSAmpl) ','... % humerus_SHSAmpl num2str(obj.shoulder.humerus.SHSOrient) ','... % humerus_SHSOrient num2str(obj.shoulder.humerus.SHSAngle) ','... % humerus_SHSAgle num2str(obj.shoulder.humerus.SHSPA) ','... % humerus_SHSPA num2str(obj.shoulder.humerus.SHSIS) ','... % humerus_SHSIS num2str(obj.shoulder.scapula.acromion.AI) ','... % acromion_AI num2str(obj.shoulder.scapula.acromion.CSA) ','... % acromion_CSA num2str(obj.shoulder.scapula.acromion.PSA) ','... % acromion_PSA num2str(obj.shoulder.scapula.acromion.AAA) '\n'... % acromion_AAA ]); fclose(fid); fclose(logFid); outputArg = 1; end function outputArg = saveExcel(~) % Save SCase to Excel file outputArg = 1; end function outputArg = saveSQL(~) % Save SCase to MySQL database outputArg = 1; end - - %{ - function outputArg = plot(obj) - % code moved to file plot.m - - % Plot SCase - figure; % Create figure - hold on; % Superpose objects on the same figure - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Scapula - - % Plot the scapula surface as points - % If exists, plot the segemented surface (manual) - if strcmp(obj.shoulder.scapula.segmentation, 'M') - points = obj.shoulder.scapula.surface.points; - faces = obj.shoulder.scapula.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',1,'EdgeColor','none','FaceLighting','gouraud'); - %scatter3(x,y,z,1,'black'); % Point points - end - % If exists, plot the segemented surface (auto) - if strcmp(obj.shoulder.scapulaAuto.segmentation, 'A') - points = obj.shoulder.scapulaAuto.surface.points; - faces = obj.shoulder.scapulaAuto.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - trisurf(faces,x,y,z,'Facecolor','red','FaceAlpha',0.25,'EdgeColor','none','FaceLighting','gouraud'); - end - - % Plot scapula groove (manual) - points = obj.shoulder.scapula.groove; - x = points(:,1); - y = points(:,2); - z = points(:,3); - scatter3(x,y,z,100,'blue'); - % If exists, plot segments between manual and auto - - % Plot scapula markers - points = [... - obj.shoulder.scapula.angulusInferior; - obj.shoulder.scapula.trigonumSpinae; - obj.shoulder.scapula.processusCoracoideus; - obj.shoulder.scapula.acromioClavicular; - obj.shoulder.scapula.angulusAcromialis; - obj.shoulder.scapula.spinoGlenoidNotch... - ]; - x = points(:,1); - y = points(:,2); - z = points(:,3); - scatter3(x,y,z,100,'blue','LineWidth',2); - - % Plot scapular plane - axisLength = 50; - pt0 = obj.shoulder.scapula.coordSys.origin; - pt1 = pt0 + (obj.shoulder.scapula.coordSys.ML + obj.shoulder.scapula.coordSys.IS) * axisLength ; - pt2 = pt1 - obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]); - pt3 = pt2 - obj.shoulder.scapula.coordSys.ML * pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae] ); - pt4 = pt3 + obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]); - X = [pt1(1) pt2(1) pt3(1) pt4(1)]; - Y = [pt1(2) pt2(2) pt3(2) pt4(2)]; - Z = [pt1(3) pt2(3) pt3(3) pt4(3)]; - patch(X,Y,Z,'black','FaceAlpha',0.3); % Transparent plane - - - % Plot scapular axis - pt1 = obj.shoulder.scapula.coordSys.origin; - axisLength = 10 + pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae]); - pt2 = pt1 - obj.shoulder.scapula.coordSys.ML * axisLength; - x = [pt1(1), pt2(1)]; - y = [pt1(2), pt2(2)]; - z = [pt1(3), pt2(3)]; - plot3(x,y,z,'Color','blue','LineWidth', 5); - - % Plot scapular coordinate system - axisLength = 50; % Length of the cood. sys. axes - pt0 = obj.shoulder.scapula.coordSys.origin; - pt1 = pt0 + obj.shoulder.scapula.coordSys.PA*axisLength; % X - pt2 = pt0 + obj.shoulder.scapula.coordSys.IS*axisLength; % Y - pt3 = pt0 + obj.shoulder.scapula.coordSys.ML*axisLength; % Z - x = [pt0(1), pt1(1)]; - y = [pt0(2), pt1(2)]; - z = [pt0(3), pt1(3)]; - plot3(x,y,z,'Color','red','LineWidth', 5); - x = [pt0(1), pt2(1)]; - y = [pt0(2), pt2(2)]; - z = [pt0(3), pt2(3)]; - plot3(x,y,z,'Color','green','LineWidth', 5); - % Not plotting ML axis to avoid its superposition with scapular axis on glenoid center -% x = [pt0(1), pt3(1)]; -% y = [pt0(2), pt3(2)]; -% z = [pt0(3), pt3(3)]; -% plot3(x,y,z,'Color','blue','LineWidth', 5); - - % plot scapular axis on glenoid center - axisLength = 50; - pt1 = obj.shoulder.scapula.glenoid.center; - pt2 = pt1 + obj.shoulder.scapula.coordSys.ML*axisLength; - x = [pt1(1), pt2(1)]; - y = [pt1(2), pt2(2)]; - z = [pt1(3), pt2(3)]; - plot3(x,y,z,'Color','blue','LineWidth', 5); - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Glenoid - - % Plot the glenoid surface - % We might move them in the direction of th glenoid center - % instead of the ML axis - points = obj.shoulder.scapula.glenoid.surface.points +... - obj.shoulder.scapula.coordSys.ML * 0.2; - faces = obj.shoulder.scapula.glenoid.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - %scatter3(x,y,z,1,'blach'); - trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','yellow','FaceLighting','none'); - - % plot glenoid centerline - pt1 = obj.shoulder.scapula.glenoid.center; - pt2 = pt1 + obj.shoulder.scapula.glenoid.centerLine; - x = [pt1(1), pt2(1)]; - y = [pt1(2), pt2(2)]; - z = [pt1(3), pt2(3)]; - plot3(x,y,z,'Color','yellow','LineWidth', 5); - - % We might add spherical cap - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Humerus - - % plot humeral head sphere - n = 20; % number of lines - [X,Y,Z] = sphere(n); - X = X*obj.shoulder.humerus.radius; - Y = Y*obj.shoulder.humerus.radius; - Z = Z*obj.shoulder.humerus.radius; - Xc = obj.shoulder.humerus.center(1); - Yc = obj.shoulder.humerus.center(2); - Zc = obj.shoulder.humerus.center(3); - mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud'); - % hidden off; % Transparent sphere - % check 'FaceLighting','gouraud', - - % plot humeral head centerline - pt1 = obj.shoulder.scapula.glenoid.center; - pt2 = obj.shoulder.humerus.center; - x = [pt1(1), pt2(1)]; - y = [pt1(2), pt2(2)]; - z = [pt1(3), pt2(3)]; - plot3(x,y,z,'Color','magenta','LineWidth', 5); - - - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - % Graph properties - - % axis off; % remove axis - lightPosition = obj.shoulder.scapula.glenoid.center + ... - (... - obj.shoulder.scapula.coordSys.ML + ... - obj.shoulder.scapula.coordSys.IS + ... - obj.shoulder.scapula.coordSys.PA ... - ) * 100; - light('Position',lightPosition,'Style','local'); - grid on; - xlabel('X (CT)'); - ylabel('Y (CT)'); - zlabel('Z (CT)'); - axis square; - axis vis3d; % fixed limits and scaling - view(obj.shoulder.scapula.coordSys.IS); % view from top - % Align ML axis with window horizontal axis !to do - % obj.shoulder.scapula.coordSys.ML - %angleDegrees = rad2deg(obj.shoulder.scapula.coordSys.ML - % camroll(angleDegrees); - - figureHandle = gca; - figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis - figureHandle.Box = 'On'; - - outputArg = 1; - - end - %} end end diff --git a/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m b/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m new file mode 100644 index 0000000..fcbdc80 --- /dev/null +++ b/ShoulderCase/@ShoulderCase/getSmoothDicomInfo.m @@ -0,0 +1,6 @@ +function output = getSmoothDicomInfo(obj) + % Return the result of dicominfo() with the first dicom file found + dicomFiles = dir(fullfile(obj.getSmoothDicomPath,'*.dcm')); + assert(not(isempty(dicomFiles)),'No dicom file found there %s',obj.getSmoothDicomPath); + output = dicominfo(fullfile(dicomFiles(1).folder,dicomFiles(1).name)); +end diff --git a/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m b/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m new file mode 100644 index 0000000..862c90d --- /dev/null +++ b/ShoulderCase/@ShoulderCase/getSmoothDicomPath.m @@ -0,0 +1,6 @@ +function output = getSmoothDicomPath(obj) + output = []; + smoothDicomPath = fullfile([obj.dataCTPath(1:end-1) '2'],'dicom'); + assert(isfolder(smoothDicomPath),'No smooth dicom for this SCase. %s is not a valid folder',smoothDicomPath); + output = smoothDicomPath; +end diff --git a/ShoulderCase/@ShoulderCase/plot.m b/ShoulderCase/@ShoulderCase/plot.m index 9af258c..78afa71 100644 --- a/ShoulderCase/@ShoulderCase/plot.m +++ b/ShoulderCase/@ShoulderCase/plot.m @@ -1,241 +1,23 @@ -function outputArg = plot(obj) -%%PLOT Plot a SCase -figure; % Create figure -hold on; % Superpose objects on the same figure - -%% Scapula -% If exists, plot the segemented surface (manual) -if strcmp(obj.shoulder.scapula.segmentation, 'M') - points = obj.shoulder.scapula.surface.points; - faces = obj.shoulder.scapula.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud'); - %scatter3(x,y,z,1,'black'); % Point points -end -% If exists, plot the segemented surface (auto) -%{ -if strcmp(obj.shoulder.scapulaAuto.segmentation, 'A') - points = obj.shoulder.scapulaAuto.surface.points; - faces = obj.shoulder.scapulaAuto.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - trisurf(faces,x,y,z,'Facecolor','red','FaceAlpha',0.25,'EdgeColor','none','FaceLighting','gouraud'); -end -%} - -% Plot scapula groove (manual) - -points = obj.shoulder.scapula.groove; -x = points(:,1); -y = points(:,2); -z = points(:,3); -scatter3(x,y,z,100,'blue'); - - -% If exists, plot segments between manual and auto -%{ -if ~isempty(obj.shoulder.scapulaAuto.groove) - pointsAuto = obj.shoulder.scapulaAuto.groove; - for iPts=1:length(pointsAuto) - x = [points(iPts,1), pointsAuto(iPts,1)]; - y = [points(iPts,2), pointsAuto(iPts,2)]; - z = [points(iPts,3), pointsAuto(iPts,3)]; - plot3(x,y,z,'Color','blue','LineWidth', 5); - end -end -%} - -% Plot scapula markers - -points = [... - obj.shoulder.scapula.angulusInferior; - obj.shoulder.scapula.trigonumSpinae; - obj.shoulder.scapula.processusCoracoideus; - obj.shoulder.scapula.acromioClavicular; - obj.shoulder.scapula.angulusAcromialis; - obj.shoulder.scapula.spinoGlenoidNotch... - ]; -x = points(:,1); -y = points(:,2); -z = points(:,3); -scatter3(x,y,z,100,'blue','LineWidth',2); - - -% If exists, plot segments between manual and auto -%{ -if ~isempty(obj.shoulder.scapulaAuto.angulusInferior) - pointsAuto = [... - obj.shoulder.scapulaAuto.angulusInferior; - obj.shoulder.scapulaAuto.trigonumSpinae; - obj.shoulder.scapulaAuto.processusCoracoideus; - obj.shoulder.scapulaAuto.acromioClavicular; - obj.shoulder.scapulaAuto.angulusAcromialis; - obj.shoulder.scapulaAuto.spinoGlenoidNotch... - ]; - for iPts=1:length(pointsAuto) - x = [points(iPts,1), pointsAuto(iPts,1)]; - y = [points(iPts,2), pointsAuto(iPts,2)]; - z = [points(iPts,3), pointsAuto(iPts,3)]; - plot3(x,y,z,'Color','blue','LineWidth', 5); - end +function output = plot(obj,varargin) + % Call a ShoulderCase.shoulder.plot function. + % + % input: (''), ('manual'), ('auto') + % input arguments are not case sensitive. + % + % output: call obj.shoulderManual.plot (default) or obj.shoulderAuto.plot. + + assert(nargin <= 2,'Too many arguments given') + if nargin == 1 + output = obj.shoulderManual.plot; + return + end + assert(ischar(varargin{1}),'Given argument must be a char array') + shoulderToPlot = lower(varargin{1}); % given argument is not case sensitive + if isequal(shoulderToPlot,'manual') + output = obj.shoulderManual.plot; + elseif isequal(shoulderToPlot,'auto') + output = obj.shoulderAuto.plot; + else + error('Wrong type of shoulder given. Input must be ''manual'' or ''auto''.'); + end end -%} - -% Plot scapular plane -axisLength = 50; -pt0 = obj.shoulder.scapula.coordSys.origin; -pt1 = pt0 + (obj.shoulder.scapula.coordSys.ML + obj.shoulder.scapula.coordSys.IS) * axisLength ; -pt2 = pt1 - obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]); -pt3 = pt2 - obj.shoulder.scapula.coordSys.ML * pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae] ); -pt4 = pt3 + obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]); -X = [pt1(1) pt2(1) pt3(1) pt4(1)]; -Y = [pt1(2) pt2(2) pt3(2) pt4(2)]; -Z = [pt1(3) pt2(3) pt3(3) pt4(3)]; -patch(X,Y,Z,'black','FaceAlpha',0.3); % Transparent plane - -% Plot scapular axis -pt1 = obj.shoulder.scapula.coordSys.origin; -axisLength = 10 + pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae]); -pt2 = pt1 - obj.shoulder.scapula.coordSys.ML * axisLength; -x = [pt1(1), pt2(1)]; -y = [pt1(2), pt2(2)]; -z = [pt1(3), pt2(3)]; -plot3(x,y,z,'Color','blue','LineWidth', 5); - -% Plot scapular coordinate system -axisLength = 50; % Length of the cood. sys. axes -pt0 = obj.shoulder.scapula.coordSys.origin; -pt1 = pt0 + obj.shoulder.scapula.coordSys.PA*axisLength; % X -pt2 = pt0 + obj.shoulder.scapula.coordSys.IS*axisLength; % Y -%pt3 = pt0 + obj.shoulder.scapula.coordSys.ML*axisLength; % Z -x = [pt0(1), pt1(1)]; -y = [pt0(2), pt1(2)]; -z = [pt0(3), pt1(3)]; -plot3(x,y,z,'Color','red','LineWidth', 5); -x = [pt0(1), pt2(1)]; -y = [pt0(2), pt2(2)]; -z = [pt0(3), pt2(3)]; -plot3(x,y,z,'Color','green','LineWidth', 5); -% Not plotting ML axis to avoid its superposition with scapular axis on glenoid center -% x = [pt0(1), pt3(1)]; -% y = [pt0(2), pt3(2)]; -% z = [pt0(3), pt3(3)]; -% plot3(x,y,z,'Color','blue','LineWidth', 5); - -% plot scapular axis on glenoid center -axisLength = 50; -pt1 = obj.shoulder.scapula.glenoid.center; -pt2 = pt1 + obj.shoulder.scapula.coordSys.ML*axisLength; -x = [pt1(1), pt2(1)]; -y = [pt1(2), pt2(2)]; -z = [pt1(3), pt2(3)]; -plot3(x,y,z,'Color','blue','LineWidth', 5); - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Glenoid - -% Plot the glenoid surface -% We might move them in the direction of th glenoid center -% instead of the ML axis - -points = obj.shoulder.scapula.glenoid.surface.points +... - obj.shoulder.scapula.coordSys.ML * 0.2; -faces = obj.shoulder.scapula.glenoid.surface.faces; -x = points(:,1); -y = points(:,2); -z = points(:,3); -%scatter3(x,y,z,1,'blach'); -trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none'); - - -% plot glenoid centerline -pt1 = obj.shoulder.scapula.glenoid.center; -pt2 = pt1 + obj.shoulder.scapula.glenoid.centerLine; -x = [pt1(1), pt2(1)]; -y = [pt1(2), pt2(2)]; -z = [pt1(3), pt2(3)]; -plot3(x,y,z,'Color','magenta','LineWidth', 5); - -% We might add spherical cap - -% plot glenoid sphere -%{ -n = 20; % number of lines -[X,Y,Z] = sphere(n); -X = X*obj.shoulder.scapula.glenoid.radius; -Y = Y*obj.shoulder.scapula.glenoid.radius; -Z = Z*obj.shoulder.scapula.glenoid.radius; -Xc = pt2(1); -Yc = pt2(2); -Zc = pt2(3); -mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud'); -%} - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Humerus - -% plot humeral head sphere -%{ -n = 20; % number of lines -[X,Y,Z] = sphere(n); -X = X*obj.shoulder.humerus.radius; -Y = Y*obj.shoulder.humerus.radius; -Z = Z*obj.shoulder.humerus.radius; -Xc = obj.shoulder.humerus.center(1); -Yc = obj.shoulder.humerus.center(2); -Zc = obj.shoulder.humerus.center(3); -mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud'); -% hidden off; % Transparent sphere -% check 'FaceLighting','gouraud', -%} - - -% plot humeral head centerline -pt1 = obj.shoulder.scapula.glenoid.center; -pt2 = obj.shoulder.humerus.center; -x = [pt1(1), pt2(1)]; -y = [pt1(2), pt2(2)]; -z = [pt1(3), pt2(3)]; -plot3(x,y,z,'Color','magenta','LineWidth', 5); - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% Graph properties - -axis off; % remove axis -lightPosition = obj.shoulder.scapula.glenoid.center + ... - (... - obj.shoulder.scapula.coordSys.ML + ... - obj.shoulder.scapula.coordSys.IS + ... - obj.shoulder.scapula.coordSys.PA ... - ) * 100; -light('Position',lightPosition,'Style','local'); -%{ -grid on; -xlabel('X (CT)'); -ylabel('Y (CT)'); -zlabel('Z (CT)'); -axis square; -%} -axis vis3d; % fixed limits and scaling -viewPoint = obj.shoulder.scapula.coordSys.ML + ... - obj.shoulder.scapula.coordSys.IS + ... - obj.shoulder.scapula.coordSys.PA; -viewPoint = obj.shoulder.scapula.coordSys.PA; -view(viewPoint); % view from lateral side -% Align IS axis with window vertical axis -camup(obj.shoulder.scapula.coordSys.IS); - - - -figureHandle = gca; -figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis -figureHandle.Box = 'On'; -set(gcf,'color','w'); % Set background color to white - -outputArg = 1; - -end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m b/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m index f41f516..c7ea718 100644 --- a/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m +++ b/ShoulderCase/@ShoulderCaseLoader/ShoulderCaseLoader.m @@ -1,111 +1,113 @@ classdef ShoulderCaseLoader < handle properties (Access = private) normal pathological end methods function obj = ShoulderCaseLoader(varargin) dataRootPath = obj.getDataRootPath(varargin); obj.normal = CasesFinderNormal(dataRootPath); obj.pathological = CasesFinderPathological(dataRootPath); end function output = containsCase(obj,SCase) output = or( obj.normal.containsCase(SCase), obj.pathological.containsCase(SCase)); end function output = getNumberOfCases(obj) output = obj.normal.getNumberOfCases + obj.pathological.getNumberOfCases; end function output = isEmpty(obj) output = not((obj.getNumberOfCases > 0)); end function output = getAllNormalCasesIDs(obj) output = obj.normal.getAllCasesIDs; end function output = getAllPathologicalCasesIDs(obj) output = obj.pathological.getAllCasesIDs; end function output = getAllCasesIDs(obj) normalIDs = obj.normal.getAllCasesIDs; pathologicalIDs = obj.pathological.getAllCasesIDs; output = {}; for i = 1:obj.normal.getNumberOfCases output{end+1} = normalIDs{i}; end for i = 1:obj.pathological.getNumberOfCases output{end+1} = pathologicalIDs{i}; end end function output = loadCase(obj,SCase) if not(obj.containsCase(SCase)) output = []; return end output = ShoulderCase(SCase); container = obj.getContainerOf(SCase); pathToCase = container.getPathTo(SCase); caseMatlabFile = fullfile(pathToCase,'matlab','SCase.mat'); if exist(caseMatlabFile,'file') loaded = load(caseMatlabFile); if not(isfield(loaded,'SCase')) return end if not(isequal(class(loaded.SCase),'ShoulderCase')) return end + loaded.SCase.dataCTPath = pathToCase; + loaded.SCase.dataPath = obj.getDataRootPath(''); output = loaded.SCase; return end end end methods (Access = private, Hidden = true) function output = getDataRootPath(obj,arguments) output = []; if isempty(arguments) TextFileOrPath = 'config.txt'; elseif not(ischar(arguments{1})) return else TextFileOrPath = arguments{1}; end if endsWith(TextFileOrPath,'.txt') output = FileExtractor.getVariableFromTextFile('dataDir',TextFileOrPath); else output = TextFileOrPath; end end function output = getContainerOf(obj,SCase) output = []; if obj.normal.containsCase(SCase) output = obj.normal; elseif obj.pathological.containsCase(SCase) output = obj.pathological; end end end end diff --git a/ShoulderCase/@Timer/Timer.m b/ShoulderCase/@Timer/Timer.m index 0d6e238..3be5866 100644 --- a/ShoulderCase/@Timer/Timer.m +++ b/ShoulderCase/@Timer/Timer.m @@ -1,58 +1,60 @@ classdef Timer properties (Constant, Hidden = true) data = LinkedList; + format = struct('string','%sh %sm %ss %sms'); end methods (Static) function output = start output = Timer.data.append(datetime('now')); end function output = stop(varargin) if (Timer.data.length == 0) error('There is no timer currently running.'); end if (nargin == 1) assert(ischar(varargin{1}),'The ID of a timer must be a char array.'); timerStartTime = Timer.data.pop(varargin{1}); else timerStartTime = Timer.data.pop; end output = Timer.getElapsedTime(datetime('now')-timerStartTime); end function clearAll Timer.data.clear; end end methods (Static, Hidden = true) function output = getElapsedTime(duration) elapsedTime = split(char(duration),':'); elapsedMilliseconds = Timer.getElapsedMilliseconds(duration); output = [elapsedTime{1} ' hours ',... elapsedTime{2} ' minutes ',... elapsedTime{3} ' seconds ',... elapsedMilliseconds ' milliseconds']; + output = sprintf(Timer.format.string,elapsedTime{1},elapsedTime{2},elapsedTime{3},elapsedMilliseconds); end function output = getElapsedMilliseconds(duration) durationInSeconds = num2str(seconds(duration)); output = extractAfter(durationInSeconds,'.'); if isempty(output); output = '000'; end end end end diff --git a/ShoulderCase/Circle.m b/ShoulderCase/Circle.m new file mode 100644 index 0000000..b89b8c6 --- /dev/null +++ b/ShoulderCase/Circle.m @@ -0,0 +1,40 @@ +classdef Circle < handle + + properties + points + numberOfPoints + radius + center + normal + end + + + + methods + function obj = Circle(center,radius,normal,numberOfPoints) + obj.center = center; + obj.radius = radius; + obj.numberOfPoints = numberOfPoints; + obj.normal = normal; + + obj.createPoints; + end + + function createPoints(obj) + obj.points = obj.radius*[cos(2*pi*(1:obj.numberOfPoints)/obj.numberOfPoints)',... + sin(2*pi*(1:obj.numberOfPoints)/obj.numberOfPoints)',... + zeros(obj.numberOfPoints,1)]; + obj.rotate(vrrotvec([0 0 1],obj.normal)); + obj.translate(obj.center); + end + + function rotate(obj,rotationVector) + obj.points = (vrrotvec2mat(rotationVector)*obj.points')'; + end + + function translate(obj,vector) + obj.points = obj.points + vector; + end + end + +end diff --git a/ShoulderCase/CoordinateSystem.m b/ShoulderCase/CoordinateSystem.m index 8b3bf95..74d483c 100644 --- a/ShoulderCase/CoordinateSystem.m +++ b/ShoulderCase/CoordinateSystem.m @@ -1,232 +1,235 @@ classdef CoordinateSystem < handle properties origin xAxis yAxis zAxis end properties (Hidden = true) rotationMatrix % enable retro-compatibility with codes that uses PA, IS, and ML axes PA IS ML end methods function obj = CoordinateSystem() obj.origin = [0 0 0]; obj.xAxis = [1 0 0]; obj.yAxis = [0 1 0]; obj.zAxis = [0 0 1]; end function set.xAxis(obj,value) % enable retro-compatibility with codes that uses a PA axis obj.xAxis = value; obj.PA = value; obj.rotationMatrix(:,1) = value'; end function set.yAxis(obj,value) % enable retro-compatibility with codes that uses an IS axis obj.yAxis = value; obj.IS = value; obj.rotationMatrix(:,2) = value'; end function set.zAxis(obj,value) % enable retro-compatibility with codes that uses a ML axis obj.zAxis = value; obj.ML = value; obj.rotationMatrix(:,3) = value'; end function output = getRotationMatrix(obj) output = obj.rotationMatrix; end - function points = express(obj,points) - % Give the coordinates of points in present coordinate system - if (size(points,2) ~= 3) - error('Points have to be a Nx3 double array.') - return; - end - points = (points-obj.origin) * obj.rotationMatrix; + function localPoints = express(obj,globalPoints) + % Give the coordinates of global cartesian points in present coordinate system + assert(size(globalPoints,2) == 3,'Points have to be a Nx3 double array.') + localPoints = (globalPoints-obj.origin) * obj.rotationMatrix; + end + + function globalPoints = getGlobalCoordinatesOfLocalPoints(obj,localPoints) + % Give the coordinates of local points in global cartesian coordinate system + assert(size(localPoints,2) == 3,'Points have to be a Nx3 double array.') + globalPoints = (localPoints * obj.rotationMatrix') + obj.origin; end function projected = projectOnXYPlane(obj,points) projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.xAxis);(obj.origin + obj.yAxis)]); projected = projectionPlane.projectOnPlane(points); end function projected = projectOnYZPlane(obj,points) projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.yAxis);(obj.origin + obj.zAxis)]); projected = projectionPlane.projectOnPlane(points); end function projected = projectOnZXPlane(obj,points) projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.zAxis);(obj.origin + obj.xAxis)]); projected = projectionPlane.projectOnPlane(points); end function output = isOrthogonal(obj) % The axes comparison can't be done looking for strict equality due to % rounded values. Therefore the values must be evaluated with a given % tolerance minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of % a millionth of a radian if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle) output = false; return end if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle) output = false; return end if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle) output = false; return end output = true; end function output = isRightHanded(obj) % The comparison between a theoretical right-handed axis and the actual % z Axis can't be done looking for strict vector equality due to rounded % values. Therefore the norm of the sum of the two vectors is evaluated. if not(obj.isOrthogonal) output = false; return end rightHandedAxis = cross(obj.xAxis,obj.yAxis); if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis) output = false; return end output = true; end function output = isLeftHanded(obj) if not(obj.isOrthogonal) output = false; return end if obj.isRightHanded output = false; return end output = true; end function setSystemWithScapulaLandmarks(obj,scapula) % Calculate the (EPFL) scapular coordinate system % The EPFL scapular coordinate system is defined for a right scapula see % paper (10.1302/0301-620X.96B4.32641). The X axis is % antero-posterior, the Y axis is infero-superior, the Z axis % is meddio-lateral. This system is right-handed. % This orientation corresponds approximatively to the ISB % recommendadtion (doi:10.1016/j.jbiomech.2004.05.042). % For left scapula we keep the x axis as postero-anterior and % get a left-handed coordinate system. if (length(scapula.groove) < 2) error(['Can''t set scapular coordinate system with actual',... ' scapula.groove']); return; end % The four following methods have to be called in this specific order obj.setXAxisWithScapulaLandmarks(scapula); obj.setZAxisWithScapulaLandmarks(scapula); obj.setYAxisWithScapulaLandmarks(scapula); obj.setOriginWithScapulaLandmarks(scapula); end function setSystemWith(obj) % % % % % % % % % % % % % % New method to set a coordinate system can be written here % % % % % % % % % % % % end end methods (Access = private, Hidden = true) function setOriginWithScapulaLandmarks(obj,scapula) % Set the origin of the coordinate system to the spino-gemoid notch % projected on the scapular axis spinoGlenoid = scapula.spinoGlenoidNotch; grooveMeanProjection = mean(scapula.plane.projectOnPlane(scapula.groove)); obj.origin = grooveMeanProjection + ... dot((spinoGlenoid - grooveMeanProjection),obj.zAxis)*obj.zAxis; end function setXAxisWithScapulaLandmarks(obj,scapula) % Set the postero-anterior axis to be the scapula plane normal obj.xAxis = scapula.plane.normal; end function setYAxisWithScapulaLandmarks(obj,scapula) % Set the infero-superior axis to be orthogonal with the two other axes superior = scapula.spinoGlenoidNotch; inferior = scapula.angulusInferior; obj.yAxis = cross(obj.xAxis,obj.zAxis); obj.yAxis = orientVectorToward(obj.yAxis,(superior-inferior)); end function setZAxisWithScapulaLandmarks(obj,scapula) % Set the media-lateral axis to be the line fitted to the projection of % the groove points on the scapula plane lateral = scapula.spinoGlenoidNotch; medial = scapula.trigonumSpinae; groovePointsProjection = scapula.plane.projectOnPlane(scapula.groove); grooveAxis = fitLine(groovePointsProjection); grooveAxis = (grooveAxis/norm(grooveAxis))'; obj.zAxis = orientVectorToward(grooveAxis,(lateral-medial)); end end end diff --git a/ShoulderCase/Plane.m b/ShoulderCase/Plane.m index 4c37d13..43815ed 100644 --- a/ShoulderCase/Plane.m +++ b/ShoulderCase/Plane.m @@ -1,69 +1,76 @@ classdef Plane < handle properties normal point fitPerformance end methods function obj = Plane() obj.normal = []; obj.point = []; obj.fitPerformance.points = []; obj.fitPerformance.residuals = []; obj.fitPerformance.RMSE = []; obj.fitPerformance.R2 = []; end + function plot(obj,radius) + % Plot a circle in the plane + circleInPlane = Circle(obj.point,radius,obj.normal,36); + XYZ = circleInPlane.points; + patch(XYZ(:,1),XYZ(:,2),XYZ(:,3),'black','FaceAlpha',0.3); % Transparent plane + end + function fit(obj,points) % Fit current plane to given points [coeff,score,roots] = pca(points); normal = cross(coeff(:,1),coeff(:,2)); normal = normal/norm(normal); meanPoint = mean(points,1); estimatedPoints = repmat(meanPoint,size(points,1),1) +... score(:,1:2)*coeff(:,1:2)'; residuals = points - estimatedPoints; error = vecnorm(residuals,2,2); RMSE = norm(error)/sqrt(size(points,1)); sumSquareError = sum(error.^2); total = vecnorm(points - meanPoint,2,2)'; sumSquareTotal = sum(total.^2); R2 = 1-(sumSquareError/sumSquareTotal); %http://en.wikipedia.org/wiki/Coefficient_of_determination obj.normal = normal'; obj.point = meanPoint; obj.setFitPerformance(points,residuals,RMSE,R2); end function setFitPerformance(obj,points,residuals,RMSE,R2) obj.fitPerformance.points = points; obj.fitPerformance.residuals = residuals; obj.fitPerformance.RMSE = RMSE; obj.fitPerformance.R2 = R2; end function setPointOnPlane(obj,point) obj.point = point; end function projectedPoints = projectOnPlane(obj,points) N2 = obj.normal.'*obj.normal; projectedPoints = points*(eye(3)-N2)+repmat(obj.point*N2,size(points,1),1); end end end diff --git a/ShoulderCase/cellArraysDifference.m b/ShoulderCase/cellArraysDifference.m new file mode 100644 index 0000000..d1802c2 --- /dev/null +++ b/ShoulderCase/cellArraysDifference.m @@ -0,0 +1,26 @@ +function output = cellArraysDifference(evaluatedArray,substractedArray,varargin) + % Remove all occurence of substracted arrays' elements from the evaluated array. + % Any number of substracted arrays can be given in arguments with at least done + % substracted array. + + assert(iscell(evaluatedArray),'Arguments must be cell arrays.') + assert(iscell(substractedArray),'Arguments must be cell arrays.') + output = {}; + + for i = 1:numel(evaluatedArray) + evaluatedElement = evaluatedArray{i}; + for j = 1:numel(substractedArray) + if isequal(evaluatedElement,substractedArray{j}) + evaluatedElement = []; + break; + end + end + if not(isempty(evaluatedElement)) + output{end+1} = evaluatedElement; + end + end + + if (nargin > 2) + output = cellArraysDifference(output,varargin{1},varargin{2:end}); + end +end diff --git a/plotSCase.m b/plotSCase.m index 423e18e..2487992 100755 --- a/plotSCase.m +++ b/plotSCase.m @@ -1,85 +1,86 @@ -function SCase = plotSCase(SCaseId) -%PLOTSCASE Plot a SCase. - -% Input: id char of shoulder case (e.g. 'P315') +function SCase = plotSCase(SCaseId,varargin) +% Call a ShoulderCase.plot function +% +% Inputs: id char of shoulder case (e.g. 'P315') +% (optional) 'manual' or 'auto' (default is 'manual'). See ShouldeCase.plot +% docstring for further explanation % % Output: Corresponding ShoulderCase object % % Example: SCase = plotSCase('P315'); % % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2018-12-30 % % TODO % Chech that the SCase exist and can be ploted %% Open log file logFileID = openLogFile('plotSCase.log'); %% Set the data directory from the configuration file config.txt -dataDir = openConfigFile('config.txt', logFileID); +dataDir = ConfigFileExtractor.getVariable('dataDir'); %% Add path of ShoulderCase classes addpath('ShoulderCase'); %% Instance of a ShoulderCase object if (exist('SCase','var') > 0) clear SCase; % Check for delete method end -SCase = ShoulderCase; % Instanciate a ShoulderCase object +SCase = ShoulderCase(SCaseId); % Instanciate a ShoulderCase object % Load the SCase % This should be a method of ShoulderCase class % The validity of the format should be checked Pnnn or Nnnn. if (numel(regexp(SCaseId,'^[PN]\d{1,3}$')) == 0) error('Invalid format of SCaseId argument. SCaseID must start with "P" or "N" and be followed by 1 to 3 digits.'); end % Directory of SCaseId SCaseDirLevel0 = SCaseId(1); % Either 'P' or 'N' strLengthSCaseId = strlength(SCaseId(2:end)); if (strLengthSCaseId < 2) SCaseDirLevel1 = '0'; % Hunderets SCaseDirLevel2 = '0'; % Dozent elseif (strLengthSCaseId < 3) SCaseDirLevel1 = '0'; % Hunderets SCaseDirLevel2 = SCaseId(2); % Dozent else SCaseDirLevel1 = SCaseId(2); % Hunderets SCaseDirLevel2 = SCaseId(3); % Dozent end % Check if a (!unique! to be done) directory exists for this SCaseId FindSCaseIdFolder = dir([dataDir '/' SCaseDirLevel0 '/' SCaseDirLevel1 '/' SCaseDirLevel2 '/' SCaseId '*']); if (isempty(FindSCaseIdFolder)) % No directory for this SCaseId error(['Missing directory for SCaseId: ' SCaseId]); end SCaseDirLevel3 = FindSCaseIdFolder.name; SCaseDir = [dataDir '/' SCaseDirLevel0 '/' SCaseDirLevel1 '/' SCaseDirLevel2 '/' SCaseDirLevel3]; % Check if this SCaseId has a CT directory with '-1' postfix (preoperative FindCTFolder = dir([SCaseDir '/' 'CT*-1']); if (isempty(FindCTFolder)) % No CT directory for this SCaseId error('Missing CT directory for SCaseId'); % ! should exist script here ! end CTDir = [SCaseDir '/' FindCTFolder.name]; matlabDir = [CTDir '/matlab']; % we should check for existance matlabFile = [matlabDir '/SCase.mat']; if exist(matlabFile, 'file') % SCaseL = load(matlabFile); % Load SCase dat file % SCase.id = SCaseL.sCase.id; % SCase.shoulder = SCaseL.sCase.shoulder; load(matlabFile); %SCase=sCase; - SCase.plot; % Plot SCase + SCase.plot(varargin{:}); % Plot SCase fprintf(logFileID, ['\n' SCaseId]); else error('Not matlab file for this SCase'); end fclose(logFileID); % Close log file end -