diff --git a/ShoulderCase/Muscles.m b/ShoulderCase/Muscles.m new file mode 100644 index 0000000..f2e67be --- /dev/null +++ b/ShoulderCase/Muscles.m @@ -0,0 +1,597 @@ +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) + SCase + end + + methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Muscles + function obj = Muscles(SCase) + %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.SCase = SCase; + end + end + + methods + function outputArg = 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 + SCase = obj.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 = SCase.shoulder.scapula.coordSys.origin; % do not mix up with origin_image, "origin" is the origin of the Scapula Coordinate System + xAxis = SCase.shoulder.scapula.coordSys.PA; + yAxis = SCase.shoulder.scapula.coordSys.IS; + zAxis = SCase.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 SCase.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 SCase.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 + + Folder_input_py='/shoulder/methods/python/rcseg/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 + + pythonDir = '/shoulder/methods/python/rcseg'; + 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); + 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.m b/ShoulderCase/Shoulder.m index 039dfe6..0783c5b 100644 --- a/ShoulderCase/Shoulder.m +++ b/ShoulderCase/Shoulder.m @@ -1,37 +1,36 @@ classdef Shoulder < handle %SHOULDER is a class to define the shoulder object, child of the %SouderCase object. % It defines the side (R/L), the scapula, humerus, muscles and CT. - + properties side % R or L scapula scapulaAuto % Same as scapula but obtained by SSM segmentation humerus muscles CTscan comment - + end properties (Hidden = true) SCase end - + methods (Access = ?ShoulderCase) % Only ShoulderCase is allowed to construct a Shoulder function obj = Shoulder(SCase) %SHOULDER Construct an instance of this class % Detailed explanation goes here obj.SCase = SCase; - + obj.side = ''; - obj.muscles = ''; + obj.muscles = Muscles(SCase); obj.CTscan = ''; obj.comment = ''; - + obj.scapula = Scapula(SCase); obj.scapulaAuto = Scapula(SCase); obj.humerus = Humerus(SCase); end end end - diff --git a/measureSCase.m b/measureSCase.m index 334c8b7..46462c3 100755 --- a/measureSCase.m +++ b/measureSCase.m @@ -1,483 +1,625 @@ function [status,message] = measureSCase(varargin) % MEASURESCASE Update the entire SCase DB with anatomical measurements. % The function can be run from (lbovenus) server with the following shell command % cd /home/shoulder/methods/matlab/database; matlab -nodisplay -nosplash -r "measureSCase;quit" % Progress can be checked through log file with following shell command % cd /home/shoulder/methods/matlab/database;tail -f log/measureSCase.log % It take ~30 min for 700 cases when run from server. % After run from server, the output file % /shoulder/data/Excel/xlsFromMatlab/anatomy.csv % must be open from Excel and saved % as xls at the same location. Excel file % /shoulder/data/Excel/ShoulderDataBase.xlsx should also be open, updated, and saved. % The script measureSCase was replacing 3 previous functions % 1) rebuildDatabaseScapulaMeasurements.m --> list of all SCases --> execute scapula_measure.m % 2) scapula_measure.m --> execute scapula_calculation and save results in Excel and SQL % 3) scapula_calculation.m --> perform anatomical analysis --> results varaiable % Input arguments can be one or several among: % {Empty,'GlenoidDensity', 'update', '[NP][1-999]', '[NP]'} % % If there is no argument the function will measure every new case. % The last 'N' or 'P' given as an argument will force the function to be % executed on all the corresponding 'N' or 'P' cases. This will overwrite % every other ShoulderCase argument '[NP][1-999]' given before. % % Examples: % measureSCase('GlenoidDensity') % is a valid instruction where every new case is measured % including its glenoid density. % measureSCase('P400','N520') % is a valid instruction where the cases 'P400' and 'N520' are % measured only if they have not already been measured yet (new % cases). % measureSCase('P400','N520','update') % is a valid instruction where the cases 'P400' and 'N520' are % measured whether they have already been measured or not. % measureSCase('GlenoidDensity','P400','N520','update','P') % is a valid instruction where all the 'P' cases are measured % including their glenoid density whether they have been % measured or not. % measureSCase('Z1590','trollingArgumentLol') % is a valid instruction which raises a warning about the % arguments and will result in the same as executing % measureSCase() i.e. every new case is calculated. % Output % File /shoulder/data/Excel/xlsFromMatlab/anatomy.csv % File /shoulder/data/Excel/xlsFromMatlab/anatomy.xls (only windows) % File /home/shoulder/data/matlab/SCaseDB.mat % File {SCaseId}.mat in each {SCaseId}/{CT}/matlab directory % logs in /log/measureSCase.log % Example: measureSCase % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2019-06-29 % TO DO: % Read first last anatomy.csv and/or SCaseDB.mat % Fill with patient and clinical data, from excel or MySQL % use argument to redo or update % Save csv within the main loop (not sure it's a good idea) % Add more message in log % Add more test to assess validity of data % Update MySQL % Rename all objects with a capital % Make it faster by not loading scapula surface if already in SCaseDB.mat % Add argument 'load' to force reload files even if they are alerady in DB % Add argument 'measure' to force re-measuring event if measurement is already in DB % Add "scapula_addmeasure.m" --> should be coded in ShoulderCase class - -tic; % Start stopwatch timer +startTime = datetime('now'); %% Open log file logFileID = openLogFile('measureSCase.log'); %% Set the data directory from the configuration file config.txt dataDir = openConfigFile('config.txt', logFileID); %% Location of the XLS ShoulderDatabase xlsDir = [dataDir '/Excel/xlsFromMatlab']; matlabDir = [dataDir '/matlab']; %% Get list of all SCase % Get list of SCases from varargin -% Delault value for inputArg +% Delault value for varargin calcGlenoidDensity = false; updateCalculations = false; specificCases = true; fprintf(logFileID, '\n\nList of SCase'); cases = {}; for argument = 1:nargin - switch varargin{argument} % Test every argument in varargin. - case 'GlenoidDensity' + switch lower(varargin{argument}) % Test every argument in varargin. + case 'glenoiddensity' calcGlenoidDensity = true; case 'update' updateCalculations = true; - case regexp(varargin{argument},'[NP]','match') + case regexp(lower(varargin{argument}),'[np]','match') specificCases = false; % Detect if argument is 'N' or 'P'. cases = {varargin{argument}}; % ends here. - case regexp(varargin{argument},'[NP][1-9][0-9]{0,2}','match') % Test for a correct argument [NP][1-999] + case regexp(lower(varargin{argument}),'[np][1-9][0-9]{0,2}','match') % Test for a correct argument [NP][1-999] if specificCases - cases{end+1} = varargin{argument}; % - for element = 1:length(cases)-1 % Check if the case has already been added to the case list before - if varargin{argument} == cases{element} % - cases = cases(1:end-1); % Delete the last element added to cases if it is already in the list + cases{end+1} = upper(varargin{argument}); % + for element = 1:length(cases)-1 % Check if the case has already been added to the case list before + if and( (upper(varargin{argument}(1))==cases{element}(1)),... % + (upper(varargin{argument}(end-2:end))==cases{element}(end-2:end)) ) + cases = cases(1:end-1); % Delete the last element added to cases if it is already in the list end end end otherwise warning(['\n "%s" is not a valid argument and has not been added to the'... - ' list of cases.\n measureSCase arguments format must be "[NP]",'... - ' "[NP][1-999]", "GlenoidDensity", or "update"'],varargin{argument}) + ' list of cases.\n measureSCase arguments format must be "[npNP]",'... + ' "[npNP][1-999]", "GlenoidDensity", or "update"'],varargin{argument}) end end SCaseList = listSCase(cases); %% Load current xls database (for clinical data) -fprintf(logFileID, '\nLoad xls database'); +fprintf(logFileID, '\nLoad xls database\n\n'); addpath('XLS_MySQL_DB'); filename = [dataDir '/Excel/ShoulderDataBase.xlsx']; excelRaw = rawFromExcel(filename); % cell array excelSCaseID = excelRaw(2:end, 1); excelDiagnosis = excelRaw(2:end, 4); excelPatientHeight = excelRaw(2:end, 23); excelGlenoidWalch = excelRaw(2:end, 55); % Lines below adapted from XLS_MySQL_DB/MainExcel2XLS % Excel.diagnosisList = diagnosisListFromExcel(Excel.Raw); % Excel.treatmentList = treatmentListFromExcel(Excel.Raw); % ExcelData.Patient = patientFromExcel(excelRaw); % Structure with patient data % Excel.shoulder = shoulderFromExcel(Excel.patient, Excel.Raw); % [Excel.SCase, Excel.diagnosis, Excel.treatment, Excel.outcome, Excel.study] = SCaseFromExcel(... % Excel.shoulder, Excel.patient, Excel.diagnosisList, Excel.treatmentList, Excel.Raw); % fprintf(logFileID, ': OK'); %% Add path to ShoulderCase class 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 SCaseDB = ShoulderCase; SCase.dataPath = dataDir; % Set dataDir for SCase + +%%%% definition of analysis variables +% +% % manual +% +ok = {}; +scapulaLoad = {}; +scapulaCoord = {}; +degeneration = {}; +glenoidLoad = {}; +glenoidMorphology = {}; +glenoidDensity = {}; +humerusLoad = {}; +humerusSubluxation = {}; +acromionMorphology = {}; +% +% % auto +% +autoOk = {}; +autoScapulaLoad = {}; +autoScapulaCoord = {}; +%autoDegeneration = {}; +autoGlenoidLoad = {}; +autoGlenoidMorphology = {}; +%autoGlenoidDensity = {}; +%autoHumerusLoad = {}; +%autoHumerusSubluxation = {}; +%autoAcromionMorphology = {}; +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %% Start loop over SCases in SCaseList + nSCase = length(SCaseList); % Number of SCases for iSCaseID = 1:nSCase SCaseID = SCaseList(iSCaseID).id; SCase(iSCaseID).id = SCaseID; SCase(iSCaseID).dataPath = dataDir; % Set dataDir for SCase output = SCase(iSCaseID).path; % Set data path of this SCase percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f'); - fprintf(logFileID, ['\n\nSCaseID: ' SCaseID ' (' percentProgress '%%)']); + fprintf(logFileID, ['\n___________________________________\n|\n| SCaseID: ' SCaseID ' (' percentProgress '%%)']); % There are 3 parts within this SCase loop: % 1) Load & analyses manual data from amira % 2) Load & analyses auto data from matlab (Statistical Shape Model) % 3) Load clinical data from Excel database, and set them to SCase % 4) Save SCase in file SCaseID/matlab/SCase.mat % 1) Load & analyses manual data from amira - if (~exist([SCase(iSCaseID).dataMatlabPath '/SCase.mat'],'file') || updateCalculations) % New calculation of SCase.mat or update the current calculation - fprintf(logFileID, '\n Segmentation manual '); - if output % Continue if amira dir exists in SCase - fprintf(logFileID, '\n Load scapula surface and landmarks'); - output = SCase(iSCaseID).shoulder.scapula.load; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Set coord. syst.'); - output = SCase(iSCaseID).shoulder.scapula.coordSysSet; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Load glenoid surface'); - output = SCase(iSCaseID).shoulder.scapula.glenoid.load; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Measure glenoid anatomy'); - output = SCase(iSCaseID).shoulder.scapula.glenoid.morphology; - if output - fprintf(logFileID, ': OK'); - if calcGlenoidDensity - fprintf(logFileID, '\n Measure glenoid density'); - output = SCase(iSCaseID).shoulder.scapula.glenoid.calcDensity; - if output - fprintf(logFileID, ': OK'); - else - fprintf(logFileID, ': Could no read dicom'); - end - end - fprintf(logFileID, '\n Load humerus data'); - output = SCase(iSCaseID).shoulder.humerus.load; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Measure humerus subluxation'); - output = SCase(iSCaseID).shoulder.humerus.subluxation; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Measure acromion anatomy'); - % Should be run after SCase.shoulder.humerus (for AI) - output = SCase(iSCaseID).shoulder.scapula.acromion.morphology; - if output - fprintf(logFileID, ': OK'); - else - fprintf(logFileID, '\n: Measure acromion anatomy error'); - end % Acromion anatomy - else - fprintf(logFileID, '\n: Measure humerus subluxation error'); - end % Subluxation humerus - else - fprintf(logFileID, '\n: Load humerus data error'); - end % Load humerus - else - fprintf(logFileID, '\n: Measure glenoid anatomy error'); - end % Anatomy glenoid - else - fprintf(logFileID, '\n: Density calculus error'); - end % Load glenoid - else - fprintf(logFileID, '\n Scapula coord syst error '); - end % Set coord. syst. - - else - fprintf(logFileID, '\n Scapula loading error'); - end % Load scapula - - else - fprintf(logFileID, '\n No amira folder'); - end % Amira path exist + + if (exist([SCase(iSCaseID).dataMatlabPath '/SCase.mat'],'file') && ~updateCalculations) % New calculation of SCase.mat or update the current calculation + continue + end + + while true + + tic; + + fprintf(logFileID, '\n| Segmentation manual '); + + if ~output % Continue if amira dir exists in SCase + fprintf(logFileID, '\n| No amira folder'); + break + end + + fprintf(logFileID, '\n| Load scapula surface and landmarks'); + output = SCase(iSCaseID).shoulder.scapula.load; + + if ~output + fprintf(logFileID, ': ERROR'); + scapulaLoad{end+1} = SCaseID; + break + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Set coord. syst.'); + output = SCase(iSCaseID).shoulder.scapula.coordSysSet; + + if ~output + fprintf(logFileID, ': ERROR'); + scapulaCoord{end+1} = SCaseID; + break; + end + + fprintf(logFileID, '\n| Measure muscle''s degeneration'); + output = SCase(iSCaseID).shoulder.muscles.degeneration(logFileID); + + switch output % output is reversed here to differenciate different types of degeneration errors + case 1 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'dicomRead() error'; + break; + case 2 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'slice() error'; + break; + case 3 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'slice is to small'; + break; + case 4 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'crop error'; + break; + case 5 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'imwrite() gave black image'; + break; + case 6 + degeneration{1,end+1} = SCaseID; + degeneration{2,end} = output; + degeneration{3,end} = 'automatic segmentation error'; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Load glenoid surface'); + output = SCase(iSCaseID).shoulder.scapula.glenoid.load; + + if ~output + fprintf(logFileID, '\n| glenoid surface loading error'); + glenoidLoad{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure glenoid anatomy'); + output = SCase(iSCaseID).shoulder.scapula.glenoid.morphology; + + if ~output + fprintf(logFileID, ': ERROR'); + glenoidMorphology{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + if calcGlenoidDensity + fprintf(logFileID, '\n| Measure glenoid density'); + output = SCase(iSCaseID).shoulder.scapula.glenoid.calcDensity; + + if ~output + fprintf(logFileID, ': Density calculus error. Could not read dicom'); + glenoidDensity{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + end + + fprintf(logFileID, '\n| Load humerus data'); + output = SCase(iSCaseID).shoulder.humerus.load; + + if ~output + fprintf(logFileID, ': ERROR'); + humerusLoad{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure humerus subluxation'); + output = SCase(iSCaseID).shoulder.humerus.subluxation; + + if ~output + fprintf(logFileID, ': ERROR'); + humerusSubluxation{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure acromion anatomy'); + % Should be run after SCase.shoulder.humerus (for AI) + output = SCase(iSCaseID).shoulder.scapula.acromion.morphology; + + if ~output + fprintf(logFileID, ': ERROR'); + acromionMorphology{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + ok{end+1} = SCaseID; + % 2) Set clinical data from Excel database to SCase - fprintf(logFileID, '\n Set clinical data from Excel'); + fprintf(logFileID, '\n| Set clinical data from Excel'); % Get idx of excelRaw for SCaseID % Use cell2struct when dots are removed from excel headers idx = find(strcmp(excelRaw, SCaseID)); % Index of SCaseID in excelRaw SCase(iSCaseID).diagnosis = excelRaw{idx,4}; SCase(iSCaseID).treatment = excelRaw{idx,9}; SCase(iSCaseID).patient.gender = excelRaw{idx,19}; SCase(iSCaseID).patient.age = excelRaw{idx,21}; SCase(iSCaseID).patient.height = excelRaw{idx,23}; SCase(iSCaseID).patient.weight = excelRaw{idx,24}; SCase(iSCaseID).patient.BMI = excelRaw{idx,25}; SCase(iSCaseID).shoulder.scapula.glenoid.walch = excelRaw{idx,55}; % Patient data --> TO DO % find patient id (in SQL) corresponding to iSCaseID % % % idxSCase = find(contains(ExcelData, SCaseID)); % SCase(iSCaseID).patient.id = ExcelData(idxSCase).patient.patient_id; % SQL patient id % SCase(iSCaseID).patient.idMed = ExcelData(idxSCase).patient.IPP; % IPP to be replace by coded (anonymized) IPP from SecuTrial % SCase(iSCaseID).patient.gender = ExcelData(idxSCase).patient.gender; % SCase(iSCaseID).patient.age = []; % Age (years) at treatment, or preop CT (we might also add birthdate wi day et at 15 of the month) % SCase(iSCaseID).patient.ethnicity = ExcelData(idxSCase).patient.IPP; % SCase(iSCaseID).patient.weight = ExcelData(idxSCase).patient.weight; % SCase(iSCaseID).patient.height = ExcelData(idxSCase).patient.height; % SCase(iSCaseID).patient.BMI = ExcelData(idxSCase).pateint.BMI; % SCase(iSCaseID).patient.comment = ExcelData(idxSCase).patient.comment; % 3) Save SCase in file SCaseID/matlab/SCase.mat - fprintf(logFileID, '\n Save SCase.mat'); + fprintf(logFileID, '\n| Save SCase.mat'); output = SCase(iSCaseID).saveMatlab; if output fprintf(logFileID, ': OK'); else - fprintf(logFileID, ': Error'); + fprintf(logFileID, ': ERROR'); end + break; end + fprintf(logFileID, '\n| Elapsed time (s): %s\n|', num2str(toc)); - % 4) Load & analyses auto data from matlab (Statistical Shape Model) - % Load scapula surface and landmarks - if (~exist([SCase(iSCaseID).dataMatlabPath '/scapulaLandmarksAutoR.mat'],'file') || updateCalculations) % New calculation of SCase.mat or update the current calculation - fprintf(logFileID, '\n Segmentation auto '); - fprintf(logFileID, '\n Load scapula surface and landmarks'); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if (exist([SCase(iSCaseID).dataMatlabPath '/scapulaLandmarksAutoR.mat'],'file') && ~updateCalculations) + continue + end + + while true + + tic; + + fprintf(logFileID, '\n| Segmentation auto '); + fprintf(logFileID, '\n| Load scapula surface and landmarks'); output = SCase(iSCaseID).shoulder.scapulaAuto.loadAuto; % Load auto data from matlab directory - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Set coord. syst.'); - output = SCase(iSCaseID).shoulder.scapulaAuto.coordSysSet; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Load glenoid surface'); - output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.loadAuto; - if output - fprintf(logFileID, ': OK'); - fprintf(logFileID, '\n Measure glenoid anatomy'); - output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.morphology; - if output - fprintf(logFileID, ': OK'); - else - fprintf(logFileID, '\n: Error in measuring glenoid anatomy'); - end - else - fprintf(logFileID, '\n: Error in loading glenoid surface'); - end - else - fprintf(logFileID, '\n: Error in setting the coordonnate System'); - end - else - fprintf(logFileID, '\n: Error in loading scapula surface and landmarks'); - end + + if ~output + fprintf(logFileID, ': ERROR'); + autoScapulaLoad{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Set coord. syst.'); + output = SCase(iSCaseID).shoulder.scapulaAuto.coordSysSet; + + if ~output + fprintf(logFileID, ': ERROR'); + autoScapulaCoord{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Load glenoid surface'); + output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.loadAuto; + + if ~output + fprintf(logFileID, ': ERROR'); + autoGlenoidLoad{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + fprintf(logFileID, '\n| Measure glenoid anatomy'); + output = SCase(iSCaseID).shoulder.scapulaAuto.glenoid.morphology; + + if ~output + fprintf(logFileID, ': ERROR'); + autoGlenoidMorphology{end+1} = SCaseID; + break; + end + + fprintf(logFileID, ': OK'); + autoOk{end+1} = SCaseID; + break; + end + fprintf(logFileID, '\n| Elapsed time (s): %s\n|___________________________________\n\n', num2str(toc)); + + % 4) Load & analyses auto data from matlab (Statistical Shape Model) + % Load scapula surface and landmarks + save('measureSCase_analysis.mat',... + 'ok','scapulaLoad','scapulaCoord',... + 'degeneration','glenoidLoad','glenoidMorphology','glenoidDensity',... + 'humerusLoad','humerusSubluxation','acromionMorphology',... + 'autoOk','autoScapulaLoad','autoScapulaCoord','autoGlenoidLoad',... + 'autoGlenoidMorphology') end % End of loop on SCaseList %% Write csv file % This might be a function (SCase.csvSave()). The input would be a filename and a structure % data % Replace header and data by structure. Currently not working %{ txtFilename = [xlsDir, '/anatomy.txt']; % Name of the txt file DataStruc = struc(... 'SCase_id', SCase.id, ... 'shoulder_side', SCase.shoulder.side,... 'glenoid_radius', SCase.shoulder.scapula.glenoid.radius,... 'glenoid_sphereRMSE', SCase.shoulder.scapula.glenoid.sphereRMSE,... 'glenoid_depth', SCase.shoulder.scapula.glenoid.depth,... 'glenoid_width', SCase.shoulder.scapula.glenoid.width,... 'glenoid_height', SCase.shoulder.scapula.glenoid.height,... 'glenoid_centerPA', SCase.shoulder.scapula.glenoid.centerLocal(1),... 'glenoid_centerIS', SCase.shoulder.scapula.glenoid.centerLocal(2),... 'glenoid_centerML', SCase.shoulder.scapula.glenoid.centerLocal(3),... 'glenoid_versionAmpl', SCase.shoulder.scapula.glenoid.versionAmpl,... 'glenoid_versionOrient', SCase.shoulder.scapula.glenoid.versionOrient,... 'glenoid_version', SCase.shoulder.scapula.glenoid.version,... 'glenoid_inclination', SCase.shoulder.scapula.glenoid.inclination,... 'humerus_jointRadius', SCase.shoulder.humerus.jointRadius,... 'humerus_headRadius', SCase.shoulder.humerus.radius,... 'humerus_GHSAmpl', SCase.shoulder.humerus.GHSAmpl,... 'humerus_GHSOrient', SCase.shoulder.humerus.GHSOrient,... 'humerus_SHSAmpl', SCase.shoulder.humerus.SHSAmpl,... 'humerus_SHSOrient', SCase.shoulder.humerus.SHSOrient,... 'humerus_SHSAngle', SCase.shoulder.humerus.SHSAngle,... 'humerus_SHSPA', SCase.shoulder.humerus.SHSPA,... 'humerus_SHSIS', SCase.shoulder.humerus.SHSIS,... 'acromion_AI', SCase.shoulder.scapula.acromion.AI,... 'acromion_CSA', SCase.shoulder.scapula.acromion.CSA,... 'acromion_PS', SCase.shoulder.scapula.acromion.PS... ); DataTable = strct2table(Data); writetable(DataTable,filename); %} % Header of the csv file dataHeader = [... 'SCase_id,' ... 'shoulder_side,' ... 'glenoid_radius,' ... 'glenoid_sphereRMSE,' ... 'glenoid_depth,' ... 'glenoid_width,' ... 'glenoid_height,' ... 'glenoid_centerPA,' ... 'glenoid_centerIS,' ... 'glenoid_centerML,' ... 'glenoid_versionAmpl,' ... 'glenoid_versionOrient,' ... 'glenoid_version,' ... 'glenoid_inclination,' ... 'humerus_jointRadius,' ... 'humerus_headRadius,' ... 'humerus_GHSAmpl,' ... 'humerus_GHSOrient,' ... 'humerus_SHSAmpl,' ... 'humerus_SHSOrient,' ... 'humerus_SHSAngle,' ... 'humerus_SHSPA,' ... 'humerus_SHSIS,' ... 'acromion_AI,' ... 'acromion_CSA,' ... 'acromion_PSA,'... 'acromion_AAA\n'... ]; updateList = listSCase; % The list of SCase to use to construct anatomy.csv and SCaseDB.mat fprintf(logFileID, '\n\nSave anatomy.csv file'); fid = fopen([xlsDir, '/anatomy.csv'],'w'); %Last csv is deleted and a new one, updated is being created fprintf(fid,dataHeader); fclose(fid); % The following lines re-create the SCaseDB.mat and contruct the anatomy.csv file S = ShoulderCase; for Case=1:length(updateList) S(Case).id = updateList(Case).id; % S(Case).dataPath = dataDir; % S(Case).path; % Compute the SCase.dataMatlabPath loaded = load([S(Case).dataMatlabPath '/SCase.mat']); % Load SCase.mat SCaseDB(Case) = loaded.SCase; SCaseDB(Case).appendToCSV('anatomy.csv'); end fprintf(logFileID, ': OK'); %% Save the entire SCaseDB array as a matlab file fprintf(logFileID, '\n\nSave SCase database'); filename = 'SCaseDB'; filename = [matlabDir '/' filename '.mat']; try save(filename, 'SCaseDB'); fprintf(logFileID, ': OK'); catch fprintf(logFileID, ': Error'); end %fprintf(logFileID, [': ' csvFilename]); %% Write xls file (Windows only) % [~,message] = csv2xlsSCase(csvFilename); % fprintf(logFileID, message); % If run from non-windows system, only csv will be produced. The xls can be % produced by opening the csv from Excel and save as xls. % Could be a function with 1 input (cvsFilenames %{ xlsFilename = strrep(csvFilename, '.csv', '.xls'); % Replace .csv by .xls if ispc % Check if run from windows! % Read cvs file as table and write the table as xls cvsTable = readtable(csvFilename); % Get table from csv (only way to read non-numeric data) cvsCell = table2cell(cvsTable); % Tranform table cell array dataHeader = cvsTable.Properties.VariableNames; % Get header cvsCell = [dataHeader; cvsCell]; % Add header sheet = 'anatomy'; % Sheet name of the xls file [status,message] = xlswrite(xlsFilename, cvsCell, sheet); % Save xls if status fprintf(logFileId, ['\nSaved ' xlsFilename]); else fprintf(logFileId, ['\nCould not save ' xlsFilename]); fprintf(logFileId, ['\n' message.identifier]); fprintf(logFileId, ['\n' message.message]); end else fprintf(logFileId, ['\n' xlsFilename ' not saved. Open and save as xls from Excel']); end %} %% Close log file fprintf(logFileID, '\n\nSCase measured: %s', num2str(length(SCaseList))); % Number of SCases measured -fprintf(logFileID, '\nElapsed time (s): %s', num2str(toc)); % stopwatch timer + +elapsedTime = char(datetime('now')-startTime); +elapsedTime = [elapsedTime(1:2) ' hours ' elapsedTime(4:5) ' minutes ' elapsedTime(7:8) ' seconds']; + +fprintf(logFileID, ['\nTotal Elapsed time: ' elapsedTime ]); fprintf(logFileID, '\n'); fclose(logFileID); % Close log file %% Output of the SCase if required by argument % SCase.output % inputArg = abaqus/density/references/ % update mySQL % SCase.sqlSave status = 1; message = 'OK'; end diff --git a/openConfigFile.m b/openConfigFile.m index 7032fca..5f18436 100644 --- a/openConfigFile.m +++ b/openConfigFile.m @@ -1,48 +1,47 @@ function dataDir = openConfigFile(configFile, logFileID) %OPENCONFIGFILE open file config.txt and read dataDir % Config file must follow specific/restrictive rules -% The file config.txt should be writen as below. +% The file config.txt should be writen as below. % One line should contain dataDir = %{ % This configuration file contains the definition of variable dataDir for database access. % It is used by functions: listSCase.m, measureSCase.m, plotScase.m % /home/shoulder/data <-- acces data dir from lbovenus or lbomars (default) % /home/shoulder/dataDev <-- acces dataDev dir from lbovenus or lbomars % /Volumes/shoulder/data <-- acces data dir from lbovenus/lbomars mounted on macos % /Volumes/shoulder/dataDev <-- acces dataDev dir from lbovenus mounted on macos % Z:\data <-- acces data dir from lbovenus/lbomars mounted on windows % Z:\dataDev <-- acces dataDev dir from lbovenus mounted on windows dataDir = /Volumes/shoulder/dataDev %} % Author: AT % Date: 2019-01-15 % Modified by: JSM, 2019-02-13 % TODO: Less restrict rules for writting config.txt dataDir = []; if nargin==2 if exist(configFile, 'file') fileID = fopen(configFile,'r'); dataDir = fscanf(fileID, '%s'); % Read config file without spaces fclose(fileID); k = strfind(dataDir, 'dataDir='); % Find definition of dataDir k = k + 8; % Remove 'dataDir=' % Check if dataDir exists dataDir = dataDir(k:end); % Assume that config file ends with dataDir content if ~exist(dataDir, 'dir') fprintf(logFileID, ['Data directory not found, check ' configFile]); % error(['Data directory not found, check ' configFile]); end end end if ~exist(dataDir, 'dir') - errordlg(sprintf("%s%s%s%s",... - "Please select manually the data directory",... - "Error: ConfigFile not provided or not found",... - "dataDir not found = ",dataDir)) - dataDir = uigetdir(pwd,"Please select the directory 'data or dataDev'"); + errordlg(sprintf('%s%s%s%s',... + 'Please select manually the data directory',... + 'Error: ConfigFile not provided or not found',... + 'dataDir not found = ',dataDir)) + %dataDir = uigetdir(pwd,'Please select the directory ''data or dataDev'''); end end -