diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m index 497deb1..f236b6b 100644 --- a/ShoulderCase/@ShoulderCase/ShoulderCase.m +++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m @@ -1,463 +1,464 @@ classdef ShoulderCase < handle % Properties and methods associated to the SoulderCase object. % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 - % Revision date: 2018-08-11 + % 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 shoulder study comment end - properties (Constant, Hidden = true) - dataPath = '../../../data'; % path to data folder from package folder - % This should be done differently. Check where we are, and where we - % should go. - 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 if nargin > 0 obj.id = SCaseId; %{ % SCaseId must be P or N followed by 1-3 digit. if ischar(SCaseId) if strlength(SCaseId) == 4 obj.id = SCaseId; else error('id must have 4 char'); end else % error('id must be a string'); end %} else % error('id is required'); end % Initialsation obj.diagnosis = ''; obj.treatment = ''; obj.outcome = ''; obj.patient = Patient(obj); obj.shoulder = Shoulder(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 = saveCSV(~) % Save SCase to csv file 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/XLS_MySQL_DB/rawFromExcel.m b/XLS_MySQL_DB/rawFromExcel.m index d979e21..74db504 100644 --- a/XLS_MySQL_DB/rawFromExcel.m +++ b/XLS_MySQL_DB/rawFromExcel.m @@ -1,36 +1,36 @@ function [ rawExcel ] = rawFromExcel(filename) %RAWFROMEXCERL read the Excel file filename and output content in raw % Set the sheet and the range % Returns the raw data fprintf('\nGet rawExcel from Excel... '); % Set Excel file name & directory -%directory = '../../../../data/Excel/'; -%filename = 'ShoulderDataBase.xlsx';filename = strcat(directory,filename); +% directory = '../../../../data/Excel/'; +% filename = 'ShoulderDataBase.xlsx';filename = strcat(directory,filename); sheet = 'SCase'; % 'Normal' & 'TSA' xlRange = 'A1:CX999'; % range of column and row to be imported [~,~,rawExcel] = xlsread(filename,sheet,xlRange,'basic'); % Get column of variables header = rawExcel(1,:); SCase_id_col = find(strcmp([header],'SCase_ID')); % Delete rows without SCase row_idx = 2; [rowN, ~] = size(rawExcel); while row_idx <= rowN row = rawExcel(row_idx, :); % get the entire row SCase_id = row{SCase_id_col}; if isnan(SCase_id) rawExcel(row_idx, :) = []; rowN = rowN - 1; else row_idx = row_idx + 1; end end fprintf('Done\n'); -end \ No newline at end of file +end diff --git a/dicominfoSCase.m b/dicominfoSCase.m index c51be93..d102c94 100755 --- a/dicominfoSCase.m +++ b/dicominfoSCase.m @@ -1,514 +1,526 @@ function outputArg = dicominfoSCase(varargin) %DICOMINFOSCASE Get info from dicom fields and save in csv, xls, mat, and SQL % This function should be started from /shoulder/methods/matlab/database % A logfile is written during execution in % /shoulder/methods/matlab/database/log/dicominfoSCase.log % If the function is called without arguments the csv file is overwritten, % otherwise it is updated. The same applies for xls, mat and SQL. % The fonction can be run from server (lbovenus) with % cd /home/shoulder/methods/matlab/database; matlab -nojvm -nodisplay -nosplash -r "dicominfoSCase;quit" % For 700 cases, it takes less than 60 seconds (when executed from server) % Progress can be be checked throuhth log file with following shell command % cd /home/shoulder/methods/matlab/database;tail -f log/dicominfoSCase.log % This funtion replaces % /home/shoulder/methods/matlab/database/CTinfoCT.m % which only works matlab R2015a on a Windows system % Syntax: outputArg = dicominfoSCase(varargin) % % Input: % varargin: can be nothing (all cases), 'N', 'P', or SCase(s) % % Output: % outputArg: 1 if ok and 0 otherwise % % Example: CTdicomInfo, CTdicominfo('P'), CTdicominfo('P315') -% Other M-files required: -% Subfunctions: listSCase +% Other M-files required: listSCase +% Subfunctions: % See also: measureSCase % % Author: Alexandre Terrier % EPFL-LBO % Creation date: 2018-07-01 -% Revision date: 2018-08-29 +% Revision date: 2018-12-30 % % TO DO: % Re-check for double slash in paths % We should implement the ShoulderCase class % We should save data in SCase/matlab directory % We might add family name, given name (to check anonymity) % We should first read csv and then update % We should update MySQL % Check 'Index exceeds matrix dimensions' when 'using dicomdisp' % Structure SCaseDicom should match csv (and xlx aand sql) firld names % rather than dicom - -%% Initialization -% This should be in a function call -dataDir = '../../../data/'; % location of data -xlsDataDir = strcat(dataDir, 'Excel/xlsFromMatlab/'); % Location of the XLS ShoulderDatabase -logDir = 'log'; % log file in current directory - %% Open log file -filename = 'dicominfoSCase.log'; -filename = [logDir '/' filename]; -logFileId = fopen(filename, 'w'); +% Check if log dir exists in working directory +logDir = 'log'; +if ~exist('log', 'dir') % If not create dir log + mkdir('log') +end +logFile = [logDir '/dicominfoSCase.log']; +logFileID = fopen(logFile, 'w'); +fprintf(logFileID, 'dicominfoSCase log file'); +fprintf(logFileID, '\nDate: '); +fprintf(logFileID, datestr(datetime)); -fprintf(logFileId, 'dicominfoSCase log file\n'); -fprintf(logFileId, '\nDate: '); -fprintf(logFileId, datestr(datetime)); -fprintf(logFileId, '\n'); -fprintf(logFileId, '\nGet dicom info of SCase'); -fprintf(logFileId, '\n'); +%% Set the data directory from the configuration file config.txt +configFile = 'config.txt'; +if ~exist(configFile, 'file') + error([configFile, ' required']); +end +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=' +dataDir = dataDir(k:end); % Assume that config file ends with dataDir content +% Check if dataDir exists +if ~exist(dataDir, 'dir') + fprintf(logFileID, ['Data directory not found, check ' configFile]); + error(['Data directory not found, check ' configFile]); +end +% Location of the XLS ShoulderDatabase +xlsDataDir = strcat(dataDir, '/Excel/xlsFromMatlab'); %% Get list of SCases if nargin == 0 % dicominfoSCase called without argument - SCaseList = listSCase; %listSCase is an external function which provides a list with all the SCases in Excel DB (commented by JSM) + SCaseList = listSCase; % All SCases in data dir with valid dicom CT dir update = 0; elseif nargin == 1 update = 1; inputArg = varargin{1}; - SCaseList = listSCase(inputArg); + SCaseList = listSCase(inputArg); % SCase related to inputArg ('N', 'P', or list) end -%% Read dicom info +%% Read dicom info tic; % Start stopwatch timer -fprintf(logFileId, '\n\nRead dicom meta data\n'); % log file message +fprintf(logFileID, '\n\nRead dicom meta data\n'); % log file message % Struct below correponds to regular dicom field names, except SCase_ID SCaseDicom = struct(... 'SCase_id' ,[],... 'PatientID' ,[],... 'PatientSex' ,[],... 'PatientAge' ,[],... 'PatientBirthDate' ,[],... 'PatientSize' ,[],... 'PatientWeight' ,[],... 'PatientFamilyName' ,[],... 'PatientGivenName' ,[],... 'AcquisitionDate' ,[],... 'ConvolutionKernel' ,[],... 'PixelSpacing' ,[],... 'SpacingBetweenSlices' ,[],... 'SliceThickness' ,[],... 'KVP' ,[],... 'XrayTubeCurrent' ,[],... 'InstitutionName' ,[],... 'StudyDescription' ,[],... 'SeriesDescription' ,[],... % Might describe shoulder side 'Manufacturer' ,[],... 'ManufacturerModelName' ,[],... 'ProtocolName' ,[]... ); % Loop over all SCases collected above for SCase_idx=1:1:numel(SCaseList) SCase_id = SCaseList(SCase_idx).id; dicomDir = [SCaseList(SCase_idx).dir '/dicom']; dicomFiles = dir([dicomDir, '/*.dcm']); dicomFile1 = [dicomDir, '/', dicomFiles(10).name]; % seems better to check 10th than 1st ! dicomInfoError = false; - fprintf(logFileId, '\n'); - fprintf(logFileId, SCaseList(SCase_idx).id); - fprintf(logFileId, ': '); + fprintf(logFileID, '\n'); + fprintf(logFileID, SCaseList(SCase_idx).id); + fprintf(logFileID, ': '); % Try first with matlab function dicominfo try dicomInfo1 = dicominfo(dicomFile1, 'UseVRHeuristic', true); - fprintf(logFileId, 'dicominfo ok'); + fprintf(logFileID, 'dicominfo ok'); % If dicominfo didn't worked, try with matlab function dicomdisp catch ME dicomInfoError = true; - fprintf(logFileId, '\n'); - fprintf(logFileId, 'dicominfo error --> using dicomdisp'); - fprintf(logFileId, ['\n' ME.message]); + fprintf(logFileID, '\n'); + fprintf(logFileID, 'dicominfo error --> using dicomdisp'); + fprintf(logFileID, ['\n' ME.message]); dicomdispStr = evalc('dicomdisp(dicomFile1)'); end % If dicominfo didn't worked, continue with output of dicomdisp if dicomInfoError % get dicom info from dicomdisp() PatientID = extractAfter(dicomdispStr,'PatientID'); PatientID = extractAfter(PatientID,'['); PatientID = extractBefore(PatientID,']'); PatientSex = extractAfter(dicomdispStr,'PatientSex'); PatientSex = extractAfter(PatientSex,'['); PatientSex = extractBefore(PatientSex,']'); PatientAge = extractAfter(dicomdispStr,'PatientAge'); PatientAge = extractAfter(PatientAge,'['); PatientAge = extractBefore(PatientAge,']'); PatientBirthDate = extractAfter(dicomdispStr,'PatientBirthDate'); PatientBirthDate = extractAfter(PatientBirthDate,'['); PatientBirthDate = extractBefore(PatientBirthDate,']'); PatientSize = extractAfter(dicomdispStr,'PatientSize'); PatientSize = extractAfter(PatientSize,'['); PatientSize = extractBefore(PatientSize,']'); PatientSize = str2double(PatientSize); PatientWeight = extractAfter(dicomdispStr,'PatientWeight'); PatientWeight = extractAfter(PatientWeight,'['); PatientWeight = extractBefore(PatientWeight,']'); PatientWeight = str2double(PatientWeight); AcquisitionDate = extractAfter(dicomdispStr,'AcquisitionDate'); AcquisitionDate = extractAfter(AcquisitionDate,'['); AcquisitionDate = extractBefore(AcquisitionDate,']'); ConvolutionKernel = extractAfter(dicomdispStr,'ConvolutionKernel'); ConvolutionKernel = extractAfter(ConvolutionKernel,'['); ConvolutionKernel = extractBefore(ConvolutionKernel,']'); PixelSpacing = extractAfter(dicomdispStr,'PixelSpacing'); PixelSpacing = extractAfter(PixelSpacing,'['); PixelSpacing = extractBefore(PixelSpacing,'\'); SpacingBetweenSlices = extractAfter(dicomdispStr,'SpacingBetweenSlices'); SpacingBetweenSlices = extractAfter(SpacingBetweenSlices,'['); SpacingBetweenSlices = extractBefore(SpacingBetweenSlices,']'); SpacingBetweenSlices = str2double(SpacingBetweenSlices); SliceThickness = extractAfter(dicomdispStr,'SliceThickness'); SliceThickness = extractAfter(SliceThickness,'['); SliceThickness = extractBefore(SliceThickness,']'); SliceThickness = str2double(SliceThickness); KVP = extractAfter(dicomdispStr,'KVP'); KVP = extractAfter(KVP,'['); KVP = extractBefore(KVP,']'); KVP = str2double(KVP); XrayTubeCurrent = extractAfter(dicomdispStr,'XrayTubeCurrent'); XrayTubeCurrent = extractAfter(XrayTubeCurrent,'['); XrayTubeCurrent = extractBefore(XrayTubeCurrent,']'); XrayTubeCurrent = str2double(XrayTubeCurrent); InstitutionName = extractAfter(dicomdispStr,'InstitutionName'); InstitutionName = extractAfter(InstitutionName,'['); InstitutionName = extractBefore(InstitutionName,']'); Manufacturer = extractAfter(dicomdispStr,'Manufacturer'); Manufacturer = extractAfter(Manufacturer,'['); Manufacturer = extractBefore(Manufacturer,']'); ManufacturerModelName = extractAfter(dicomdispStr,'ManufacturerModelName'); ManufacturerModelName = extractAfter(ManufacturerModelName,'['); ManufacturerModelName = extractBefore(ManufacturerModelName,']'); ProtocolName = extractAfter(dicomdispStr,'ProtocolName'); ProtocolName = extractAfter(ProtocolName,'['); ProtocolName = extractBefore(ProtocolName,']'); SeriesDescription = extractAfter(dicomdispStr,'SeriesDescription'); SeriesDescription = extractAfter(SeriesDescription,'['); SeriesDescription = extractBefore(SeriesDescription,']'); else % dicomInfo exists if isfield(dicomInfo1,'PatientID') PatientID = dicomInfo1.PatientID; else PatientID = ''; end if isfield(dicomInfo1,'PatientSex') PatientSex = dicomInfo1.PatientSex; else PatientSex = ''; end if isfield(dicomInfo1,'PatientBirthDate') PatientBirthDate = dicomInfo1.PatientBirthDate; else PatientBirthDate = ''; end if isfield(dicomInfo1,'PatientAge') PatientAge = dicomInfo1.PatientAge; else PatientAge = ''; end if isfield(dicomInfo1,'PatientSize') PatientSize = dicomInfo1.PatientSize; else PatientSize = ''; end if isfield(dicomInfo1,'PatientWeight') PatientWeight = dicomInfo1.PatientWeight; else PatientWeight=''; end if isfield(dicomInfo1,'AcquisitionDate') AcquisitionDate = dicomInfo1.AcquisitionDate; else AcquisitionDate = ''; end if isfield(dicomInfo1,'ConvolutionKernel') ConvolutionKernel = dicomInfo1.ConvolutionKernel; else ConvolutionKernel = ''; end if isfield(dicomInfo1,'PixelSpacing') PixelSpacing = dicomInfo1.PixelSpacing(1); PixelSpacing = num2str(PixelSpacing); % PixelSpacing = [,PixelSpacing,]; % Commented because I (AT) don't % understand the effect else PixelSpacing=''; end if isfield(dicomInfo1,'SpacingBetweenSlices') SpacingBetweenSlices = dicomInfo1.SpacingBetweenSlices; else % SpacingBetweenSlices estimated from 2 slices % fprintf(logFileID, '\n'); - fprintf(logFileId, ', SpacingBetweenSlices estimated from 2 slices'); + fprintf(logFileID, ', SpacingBetweenSlices estimated from 2 slices'); dicomFile2 = [dicomDir, '/', dicomFiles(11).name]; dicomInfo2 = dicominfo(dicomFile2, 'UseVRHeuristic', true); ImagePositionPatient1 = dicomInfo1.ImagePositionPatient(3); ImagePositionPatient2 = dicomInfo2.ImagePositionPatient(3); posdDiff = abs(ImagePositionPatient2 - ImagePositionPatient1); InstanceNumber1 = dicomInfo1.InstanceNumber; InstanceNumber2 = dicomInfo2.InstanceNumber; InstanceNumberDiff = abs(InstanceNumber2 - InstanceNumber1); if InstanceNumberDiff == 0 InstanceNumberDiff=1; end SpacingBetweenSlices = posdDiff/InstanceNumberDiff; end if isfield(dicomInfo1,'SliceThickness') SliceThickness = dicomInfo1.SliceThickness; else SliceThickness=''; end if isfield(dicomInfo1,'KVP') KVP = dicomInfo1.KVP; else KVP= ''; end if isfield(dicomInfo1,'XrayTubeCurrent') XrayTubeCurrent = dicomInfo1.XrayTubeCurrent; else XrayTubeCurrent = ''; end if isfield(dicomInfo1,'InstitutionName') InstitutionName = dicomInfo1.InstitutionName; else InstitutionName = ''; end if isfield(dicomInfo1,'Manufacturer') Manufacturer = dicomInfo1.Manufacturer; else Manufacturer =''; end if isfield(dicomInfo1,'ManufacturerModelName') ManufacturerModelName = dicomInfo1.ManufacturerModelName; else ManufacturerModelName = ''; end if isfield(dicomInfo1,'ProtocolName') ProtocolName = dicomInfo1.ProtocolName; else ProtocolName = ''; end if isfield(dicomInfo1,'SeriesDescription') SeriesDescription = dicomInfo1.SeriesDescription; else SeriesDescription = ''; end % End of reading dimcom info end %% Format data PatientBirthDate = datetime(PatientBirthDate,'InputFormat','yyyyMMdd'); PatientBirthDate = datestr(PatientBirthDate, 'yyyy-mm-dd'); PatientAge = str2double(strtok(PatientAge,'Y')); PatientSize = round(100 * PatientSize); if PatientSize == 0 PatientSize = ''; end PatientWeight = round(PatientWeight); if PatientWeight == 0 PatientWeight = ''; end AcquisitionDate = datetime(AcquisitionDate,'InputFormat','yyyyMMdd'); AcquisitionDate = datestr(AcquisitionDate, 'yyyy-mm-dd'); PixelSpacing = str2double(PixelSpacing); PixelSpacing = round(PixelSpacing, 3); SpacingBetweenSlices = round(SpacingBetweenSlices, 3); SliceThickness = round(SliceThickness, 3); KVP = round(KVP); XrayTubeCurrent = round(XrayTubeCurrent); % Replace any ',' by ' ' to avoid problem in csv InstitutionName = strrep(InstitutionName,',',' '); Manufacturer = strrep(Manufacturer,',',' '); ManufacturerModelName = strrep(ManufacturerModelName,',',' '); ProtocolName = strrep(ProtocolName,',',' '); SeriesDescription = strrep(SeriesDescription,',',' '); %% Add data in SCaseDicom structure % This could be directly in SCase.CT object SCaseDicom(SCase_idx).id = SCase_id; SCaseDicom(SCase_idx).PatientID = PatientID; SCaseDicom(SCase_idx).PatientSex = PatientSex; SCaseDicom(SCase_idx).PatientBirthDate = PatientBirthDate; SCaseDicom(SCase_idx).PatientAge = PatientAge; SCaseDicom(SCase_idx).PatientSize = PatientSize; SCaseDicom(SCase_idx).PatientWeight = PatientWeight; SCaseDicom(SCase_idx).AcquisitionDate = AcquisitionDate; SCaseDicom(SCase_idx).ConvolutionKernel = ConvolutionKernel; SCaseDicom(SCase_idx).PixelSpacing = PixelSpacing; SCaseDicom(SCase_idx).SpacingBetweenSlices = SpacingBetweenSlices; SCaseDicom(SCase_idx).SliceThickness = SliceThickness; SCaseDicom(SCase_idx).KVP = KVP; SCaseDicom(SCase_idx).XrayTubeCurrent = XrayTubeCurrent; SCaseDicom(SCase_idx).InstitutionName = InstitutionName; SCaseDicom(SCase_idx).Manufacturer = Manufacturer; SCaseDicom(SCase_idx).ManufacturerModelName = ManufacturerModelName; SCaseDicom(SCase_idx).ProtocolName = ProtocolName; SCaseDicom(SCase_idx).SeriesDescription = SeriesDescription; end -fprintf(logFileId, '\n\nElapsed time (s): %s', num2str(toc)); % stopwatch timer +fprintf(logFileID, '\n\nElapsed time (s): %s', num2str(toc)); % stopwatch timer %% Write csv file % This might be a function. The input would be a filename and a structure % data dataHeader = [... 'SCase_ID,' ... 'patient_IPP,' ... 'patient_gender,' ... 'patient_birthDate,' ... 'patient_age,' ... 'patient_height,' ... 'patient_weight,' ... 'CT_date,' ... 'CT_kernel,' ... 'CT_pixelSize,' ... 'CT_sliceSpacing,' ... 'CT_sliceThickness,' ... 'CT_tension,' ... 'CT_current,' ... 'CT_institution,' ... 'CT_manufacturer,' ... 'CT_model,' ... 'CT_protocol,' ... 'CT_SeriesDescription' ... ]; -csvFilename = strcat(xlsDataDir, 'CTdicomInfo.csv'); +csvFilename = strcat(xlsDataDir, '/CTdicomInfo.csv'); % If dicominfoSCase is called with argument, we have to update the csv. if update % Read the csv and put data in a structure like SCaseDicom % We should chech that file exists csvTable = readtable(csvFilename,'DatetimeType', 'text'); % If 'DatetimeType' is specified as 'text', then the type for imported % date and time data depends on the value specified in the 'TextType' % parameter: If 'TextType' is 'char', then readtable returns dates as a % cell array of character vectors. If 'TextType' is 'string', then % readtable returns dates as an array of strings. % Set csv data in temporary SCaseDicomTmp structure SCaseDicom_csv = table2struct(csvTable); % Update the SCaseDicomTmp structure with new data for iSCaseDicom = 1:length(SCaseDicom) % loop on new SCase(s) - iSCaseDicomTmp = strcmp({SCaseDicom_csv.SCase_id}, SCaseDicom(iSCaseDicom).id); % 0/1 array + iSCaseDicomTmp = strcmp({SCaseDicom_csv.SCase_ID}, SCaseDicom(iSCaseDicom).id); % 0/1 array iSCaseDicomTmp = find(iSCaseDicomTmp); % index with non-zero, correponding to SCaseID % We should test iSCaseDicomTmp SCaseDicom_csv(iSCaseDicomTmp).patient_IPP = SCaseDicom(iSCaseDicom).PatientID; SCaseDicom_csv(iSCaseDicomTmp).patient_gender = SCaseDicom(iSCaseDicom).PatientSex; SCaseDicom_csv(iSCaseDicomTmp).patient_birthDate = SCaseDicom(iSCaseDicom).PatientBirthDate; SCaseDicom_csv(iSCaseDicomTmp).patient_age = SCaseDicom(iSCaseDicom).PatientAge; SCaseDicom_csv(iSCaseDicomTmp).patient_height = SCaseDicom(iSCaseDicom).PatientSize; SCaseDicom_csv(iSCaseDicomTmp).patient_weight = SCaseDicom(iSCaseDicom).PatientWeight; SCaseDicom_csv(iSCaseDicomTmp).CT_date = SCaseDicom(iSCaseDicom).AcquisitionDate; SCaseDicom_csv(iSCaseDicomTmp).CT_kernel = SCaseDicom(iSCaseDicom).ConvolutionKernel; SCaseDicom_csv(iSCaseDicomTmp).CT_pixelSize = SCaseDicom(iSCaseDicom).PixelSpacing; SCaseDicom_csv(iSCaseDicomTmp).CT_sliceSpacing = SCaseDicom(iSCaseDicom).SpacingBetweenSlices; SCaseDicom_csv(iSCaseDicomTmp).CT_sliceThickness = SCaseDicom(iSCaseDicom).SliceThickness; SCaseDicom_csv(iSCaseDicomTmp).CT_tension = SCaseDicom(iSCaseDicom).KVP; SCaseDicom_csv(iSCaseDicomTmp).CT_current = SCaseDicom(iSCaseDicom).XrayTubeCurrent; SCaseDicom_csv(iSCaseDicomTmp).CT_institution = SCaseDicom(iSCaseDicom).InstitutionName; SCaseDicom_csv(iSCaseDicomTmp).CT_manufacturer = SCaseDicom(iSCaseDicom).Manufacturer; SCaseDicom_csv(iSCaseDicomTmp).CT_model = SCaseDicom(iSCaseDicom).ManufacturerModelName; SCaseDicom_csv(iSCaseDicomTmp).CT_protocol = SCaseDicom(iSCaseDicom).ProtocolName; SCaseDicom_csv(iSCaseDicomTmp).CT_SeriesDescription = SCaseDicom(iSCaseDicom).SeriesDescription; end % SCaseDicom must now be updated for iSCaseDicom = 1:length(SCaseDicom) - SCaseDicom(iSCaseDicom).id = SCaseDicom_csv(iSCaseDicom).SCase_id; + SCaseDicom(iSCaseDicom).id = SCaseDicom_csv(iSCaseDicom).SCase_ID; SCaseDicom(iSCaseDicom).PatientID = SCaseDicom_csv(iSCaseDicom).patient_IPP; SCaseDicom(iSCaseDicom).PatientSex = SCaseDicom_csv(iSCaseDicom).patient_gender; SCaseDicom(iSCaseDicom).PatientBirthDate = SCaseDicom_csv(iSCaseDicom).patient_birthDate; SCaseDicom(iSCaseDicom).PatientAge = SCaseDicom_csv(iSCaseDicom).patient_age; SCaseDicom(iSCaseDicom).PatientSize = SCaseDicom_csv(iSCaseDicom).patient_height; SCaseDicom(iSCaseDicom).PatientWeight = SCaseDicom_csv(iSCaseDicom).patient_weight; SCaseDicom(iSCaseDicom).AcquisitionDate = SCaseDicom_csv(iSCaseDicom).CT_date; SCaseDicom(iSCaseDicom).ConvolutionKernel = SCaseDicom_csv(iSCaseDicom).CT_kernel; SCaseDicom(iSCaseDicom).PixelSpacing = SCaseDicom_csv(iSCaseDicom).CT_pixelSize; SCaseDicom(iSCaseDicom).SpacingBetweenSlices = SCaseDicom_csv(iSCaseDicom).CT_sliceSpacing; SCaseDicom(iSCaseDicom).SliceThickness = SCaseDicom_csv(iSCaseDicom).CT_sliceThickness; SCaseDicom(iSCaseDicom).KVP = SCaseDicom_csv(iSCaseDicom).CT_tension; SCaseDicom(iSCaseDicom).XrayTubeCurrent = SCaseDicom_csv(iSCaseDicom).CT_current; SCaseDicom(iSCaseDicom).InstitutionName = SCaseDicom_csv(iSCaseDicom).CT_institution; SCaseDicom(iSCaseDicom).Manufacturer = SCaseDicom_csv(iSCaseDicom).CT_manufacturer; SCaseDicom(iSCaseDicom).ManufacturerModelName = SCaseDicom_csv(iSCaseDicom).CT_model; SCaseDicom(iSCaseDicom).ProtocolName = SCaseDicom_csv(iSCaseDicom).CT_protocol; SCaseDicom(iSCaseDicom).SeriesDescription = SCaseDicom_csv(iSCaseDicom).CT_SeriesDescription; end end % Write csv file csvFileId = fopen(csvFilename, 'w'); fprintf(csvFileId,dataHeader); for SCase_idx=1:1:numel(SCaseDicom) % Loop on all SCase provided by listSCase() fprintf(csvFileId,'\n%s,%s,%s,%s,%f,%f,%f,%s,%s,%f,%f,%f,%f,%f,%s,%s,%s,%s,%s',... SCaseDicom(SCase_idx).id,... SCaseDicom(SCase_idx).PatientID,... SCaseDicom(SCase_idx).PatientSex,... SCaseDicom(SCase_idx).PatientBirthDate,... SCaseDicom(SCase_idx).PatientAge,... SCaseDicom(SCase_idx).PatientSize,... SCaseDicom(SCase_idx).PatientWeight,... SCaseDicom(SCase_idx).AcquisitionDate,... SCaseDicom(SCase_idx).ConvolutionKernel,... SCaseDicom(SCase_idx).PixelSpacing,... SCaseDicom(SCase_idx).SpacingBetweenSlices,... SCaseDicom(SCase_idx).SliceThickness,... SCaseDicom(SCase_idx).KVP,... SCaseDicom(SCase_idx).XrayTubeCurrent,... SCaseDicom(SCase_idx).InstitutionName,... SCaseDicom(SCase_idx).Manufacturer,... SCaseDicom(SCase_idx).ManufacturerModelName,... SCaseDicom(SCase_idx).ProtocolName,... SCaseDicom(SCase_idx).SeriesDescription... ); end fclose(csvFileId); -fprintf(logFileId, ['\n\nSaved csv file: ' csvFilename]); +fprintf(logFileID, ['\n\nSaved csv file: ' csvFilename]); %% Write xls file (Windows only) [~,message] = csv2xlsSCase(csvFilename); -fprintf(logFileId, message); +fprintf(logFileID, message); %% Update MySQL % To be done %{ addpath('XLS_MySQL_DB'); MainExcel2SQL; %} %% Final log -fprintf(logFileId, '\n\nSCases treated: %s', num2str(numel(SCaseDicom))); % Number of SCases treated -fprintf(logFileId, '\nElapsed time (s): %s', num2str(toc)); % stopwatch timer -fclose(logFileId); % Close log file +fprintf(logFileID, '\n\nSCases treated: %s', num2str(numel(SCaseDicom))); % Number of SCases treated +fprintf(logFileID, '\nElapsed time (s): %s', num2str(toc)); % stopwatch timer +fclose(logFileID); % Close log file outputArg = 1; % Might be extended end \ No newline at end of file diff --git a/listSCase.m b/listSCase.m index feda8b6..3234604 100755 --- a/listSCase.m +++ b/listSCase.m @@ -1,264 +1,230 @@ function SCaseList = listSCase(varargin) -%LISTSCASE output list of sCases (id and path), accordind to inputList and +%LISTSCASE output list of sCases (id and path), accordind to inputList or %filter arguments. % Each sCase directory should at least contain at % directory with a name staring with 'CT' and ending with either % '-1' : preop/shoulder/sharp/noPhantom % '-2' : preop/shoulder/smooth/noPhantom % '-3' : preop/shoulder/sharp/phantom % '-4' : preop/shoulder/smooth/phantom % and containing a 'dicom' subdirectory. The priority is '1', '2', '3', '4', % unless there is a amira folder already present in one of these diretories. % -% The optional input argument inputList if an array of SCaseID. If empty, -% listSCase will search for all SCaseID. The optional input argument -% filter can restrict le list. +% The optional input argument inputList if an array of SCaseID, or 'N' or +% 'P'. If empty, listSCase will search for all SCaseID. The optional +% input argument filter can restrict le list. % The fonction can be run from server (lbovenus) with % cd /home/shoulder/methods/matlab/database; matlab -nojvm -nodisplay -nosplash -r "listSCase;quit" % For 700 cases, it takes less than 60 seconds (when executed from server) % Progress can be be checked throuhth log file with following shell command % cd /home/shoulder/methods/matlab/database;tail -f log/listSCase.log - - % Input: % inputArg: can be nothing (=all), 'N', 'P', or SCase(s) % % Output: % outputArg: 1 if ok and 0 otherwise % % Example: listSCase, listSCase('P'), listSCase('P001','P002') -% Other M-files required: -% Subfunctions: -% See also: +% Other M-files required: None +% Subfunctions: None +% See also: Used by dicominfoSCase, measureSCase, % % Author: Alexandre Terrier % EPFL-LBO % Creation date: 2018-07-01 -% Revision date: 2018-08-09 -% -% TO DO: -% If no dir CT*-1, check for other, wit the priority with one having amira -% directory +% Revision date: 2018-12-30 +% +% TODO +% Two main looops might be merged in one, to avoid duplicagte post +% treatment -%% Initialisation of paths -% Should be replaced by a function -currentDir = pwd; -rootDir = strfind(currentDir, 'methods/matlab'); % Check it for windows -rootDir = currentDir(1:rootDir-2); -% We should check here that currentDir contains 'Dev' to associate dataDev -dataDir = [rootDir '/data']; % location of data -%dataDir = '../../../data'; % location of data --> remove if above ok -logDir = 'log'; % log file in xls folder - %% Open log file -logFile = [logDir '/listSCase.log']; -logFileId = fopen(logFile, 'w'); -fprintf(logFileId, 'listSCase log file'); -fprintf(logFileId, '\nDate: '); -fprintf(logFileId, datestr(datetime)); +logFileID = openLogFile('listSCase.log'); + +%% Set the data directory from the configuration file config.txt +dataDir = openConfigFile('config.txt', logFileID); %% Struct contains id and associated directory SCaseList = struct(... 'id' ,[],... 'dir' ,[],... 'comment' ,[]... % Might be used later ); %% Check varargin dirL0Array = ''; SCaseListIn = {}; if nargin == 0 dirL0Array = ['N';'P']; % All SCase, default when no argument - fprintf(logFileId, '\nGet list of all SCase with a valid dicom folder'); + fprintf(logFileID, '\nGet list of all SCase with a valid dicom folder'); elseif nargin == 1 if strcmp(varargin,'N') % All normal SCase dirL0Array = 'N'; - fprintf(logFileId, '\nGet list of all normal SCase with a valid dicom folder'); + fprintf(logFileID, '\nGet list of all normal SCase with a valid dicom folder'); elseif strcmp(varargin,'P') % All pathological SCase dirL0Array = 'P'; - fprintf(logFileId, '\nGet list of all pathological SCase with a valid dicom folder'); + fprintf(logFileID, '\nGet list of all pathological SCase with a valid dicom folder'); else % A specific SCase - fprintf(logFileId, '\nGet list of a SCase with a valid dicom folder'); + fprintf(logFileID, '\nGet list of a SCase with a valid dicom folder'); SCaseListIn = varargin; end else % list of multiple SCase given in varargin - fprintf(logFileId, '\nGet list of a SCase with a valid dicom folder'); + fprintf(logFileID, '\nGet list of a SCase with a valid dicom folder'); SCaseListIn = varargin; end %% Only one of the two loops below are performed SCaseListIdx = 0; % Used in one of two loops below %% Loop over SCase (N, P, or both), if length(dirL0Array) > 1 for dirL0idx = 1:length(dirL0Array) % N or P - for dirL1idx = 0:9 % Hunderts - for dirL2idx = 0:9 % Tens + for dirL1idx = 0:9 % Hunderts dir + for dirL2idx = 0:9 % Tens dir dirL0Str = dirL0Array(dirL0idx); dirL1Str = num2str(dirL1idx); dirL2Str = num2str(dirL2idx); dirTens = [dataDir '/' dirL0Str '/' dirL1Str '/' dirL2Str]; - dirTensList = dir(dirTens); % List all directories in tens directory + dirTensList = dir(dirTens); % List all directories in tens dir dirTensList = dirTensList(~ismember({dirTensList.name},{'.','..'})); % Remove . and .. % We might add here a filter (only diretories, {N/P} followed % by digits for dirL3idx = 1:numel(dirTensList) sCaseDirName = dirTensList(dirL3idx).name; % Root directory of this case sCaseDirPath = dirTensList(dirL3idx).folder; - % Check that this sCase directory contains a directory with the - % name [CT?*-1] and with a dicom directory in it. + % Check if this SCaseID has a CT dir [CT?*-?] with a valid preoprative CT SCaseID = strtok(sCaseDirName, '-'); % Get chars before '-' - - % Check if this SCaseID has a CT directory with a valid preoprative CT CTdirList = dir([sCaseDirPath '/' sCaseDirName '/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 - fprintf(logFileId, ['\n', SCaseID, ' has no valid dicom directory']); + fprintf(logFileID, ['\n', SCaseID, ' has no valid dicom directory']); else if iCTdirAmira % If amira dir exist, use it iCTdir2use = iCTdirAmira; end CTdir = CTdirList(iCTdir2use); SCaseListIdx = SCaseListIdx + 1; SCaseList(SCaseListIdx).id = SCaseID; CTdirPath = [CTdir.folder '/' CTdir.name]; SCaseList(SCaseListIdx).dir = CTdirPath; if CTnum ~= 1 - fprintf(logFileId, ['\n', SCaseID, ' dicom directory is CT-' num2str(CTnum)]); + fprintf(logFileID, ['\n', SCaseID, ' dicom directory is CT-' num2str(CTnum)]); end - end - - - - -% % Listing of (CT) folders in sCaseDirPath staring with CT -% % and ending with '-1' -% CTdir = dir([sCaseDirPath '/' sCaseDirName '/CT*-1']); -% if ~isempty(CTdir) -% CTdirPath = [CTdir.folder '/' CTdir.name]; -% dicomDirName = [CTdirPath '/dicom']; -% if exist(dicomDirName, 'dir') -% SCaseListIdx = SCaseListIdx + 1; -% SCaseList(SCaseListIdx).id = SCaseID; -% SCaseList(SCaseListIdx).dir = CTdirPath; -% % We might get here the readme.txt if present -% else -% fprintf(logFileId, ['\n', SCaseID, ' has no valid dicom directory']); -% end -% else -% fprintf(logFileId, ['\n', SCaseID, ' has no CT*-1 directory']); -% end end end end end %% Loop over input SCase list (given as input argument), if length(SCaseListIn) > 0 for iSCaseListIn = 1:length(SCaseListIn) % Check if there is a valid dicom directory for this SCaseID SCaseID = char(SCaseListIn{iSCaseListIn}); % Build directory from SCaseID (as in function ShoulderCase.path) % The validity of the format should be either Pnnn or Nnnn. if (numel(regexp(SCaseID,'^[PN]\d{1,3}$')) == 0) error(['Invalid format of SCaseID: ' SCaseID{1} '. 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)]; % Check if a (!unique! to be done) directory exists for this SCaseID 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]; % Check if this SCaseID has a CT directory with a valid preoprative CT 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 - fprintf(logFileId, ['\n', SCaseID, ' has no valid dicom directory']); + fprintf(logFileID, ['\n', SCaseID, ' has no valid dicom directory']); else if iCTdirAmira % If amira dir exist, use it iCTdir2use = iCTdirAmira; end CTdir = CTdirList(iCTdir2use); SCaseListIdx = SCaseListIdx + 1; SCaseList(SCaseListIdx).id = SCaseID; CTdirPath = [CTdir.folder '/' CTdir.name]; SCaseList(SCaseListIdx).dir = CTdirPath; + if CTnum ~= 1 + fprintf(logFileID, ['\n', SCaseID, ' dicom directory is CT-' num2str(CTnum)]); + end end end -fprintf(logFileId, ['\nNumber of SCase: ' num2str(length(SCaseList))]); +fprintf(logFileID, ['\nNumber of SCase: ' num2str(length(SCaseList))]); +fclose(logFileID); % Close log file end diff --git a/measureSCase.m b/measureSCase.m index dfb1398..14bc9b0 100755 --- a/measureSCase.m +++ b/measureSCase.m @@ -1,433 +1,450 @@ function [status,message] = measureSCase(varargin) % MEASURESCASE Update the entire SCase DB with anatimical 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. % The script measureSCase will replace 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 % The current version produces correctly the entire csv when run % without arguments only (for the entire database) % Input % Either empty, 'N', 'P', or a SCase ('P315') % 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: 2018-09-18 +% Revision date: 2018-12-31 % 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 -%% Initialisation of paths -% Should be replaced by a function -dataDir = '../../../data'; % location of data -logDir = 'log'; % log file in xls folder -xlsDir = [dataDir '/Excel/xlsFromMatlab']; -matlabDir = [dataDir '/matlab']; +%% Open log file +logFileID = openLogFile('measureSCase.log'); -%% Add path to ShoulderCase class -addpath('ShoulderCase'); +%% Set the data directory from the configuration file config.txt +dataDir = openConfigFile('config.txt', logFileID); -%% Open log file -filename = 'measureSCase.log'; -filename = [logDir '/' filename]; -logFileID = fopen(filename, 'w'); -fprintf(logFileID, 'Log file of measureSCase'); -fprintf(logFileID, ['\nDate: ' datestr(datetime)]); -%% Instance of a ShoulderCase object -if (exist('SCase','var') > 0) - clear SCase; % Check for delete method -end -SCase = ShoulderCase; % Instanciate a ShoulderCase object + +% configFile = 'config.txt'; +% if ~exist(configFile, 'file') +% error([configFile, ' required']); +% end +% 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=' +% dataDir = dataDir(k:end); % Assume that config file ends with dataDir content +% % Check if dataDir exists +% if ~exist(dataDir, 'dir') +% fprintf(logFileID, ['Data directory not found, check ' configFile]); +% error(['Data directory not found, check ' configFile]); +% end + + +%% Location of the XLS ShoulderDatabase +xlsDir = [dataDir '/Excel/xlsFromMatlab']; +matlabDir = [dataDir '/matlab']; %% Get list of all SCase % Get list of SCases from varargin fprintf(logFileID, '\n\nList of SCase'); if nargin == 0 % findScapulaLandmarksSCase called without argument SCaseList = listSCase; elseif nargin == 1 inputArg = varargin{1}; SCaseList = listSCase(inputArg); else error('inputArg can be empty, N, P, or a SCaseID') end %% Load current xls database (for clinical data) fprintf(logFileID, '\nLoad xls database'); 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 +SCase.dataPath = dataDir; % Set dataDir for SCase + %% 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 percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f'); fprintf(logFileID, ['\n\nSCaseID: ' 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 fprintf(logFileID, '\n Segmentation manual '); output = SCase(iSCaseID).path; % Set data path of this SCase 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'); fprintf(logFileID, '\n Measure glenoid density'); output = SCase(iSCaseID).shoulder.scapula.glenoid.calcDensity; if output fprintf(logFileID, ': OK'); 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, ': Error'); end % Acromion anatmy else fprintf(logFileID, ': Error'); end % Subluxation humerus else fprintf(logFileID, ': Error'); end % Load humerus else fprintf(logFileID, ': Error'); end % Anatomy glenoid else fprintf(logFileID, ': Error'); end % Load glenoid else fprintf(logFileID, '\n No amira folder'); end % Set coord. syst. end % Load scapula end % Amira path exist % 2) Load & analyses auto data from matlab (Statistical Shape Model) fprintf(logFileID, '\n Segmentation auto '); % Load scapula surface and landmarks 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, ': Error'); end else fprintf(logFileID, ': Error'); end else fprintf(logFileID, ': Error'); end else fprintf(logFileID, ': Error'); end % 3) Set clinical data from Excel database to SCase 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; % 4) Save SCase in file SCaseID/matlab/SCase.mat fprintf(logFileID, '\n Save SCase'); output = SCase(iSCaseID).saveMatlab; if output fprintf(logFileID, ': OK'); else fprintf(logFileID, ': Error'); end end % End of loop on SCaseList %% Save the entire SCase array as a matlab file fprintf(logFileID, '\n\n Save SCase database'); filename = 'SCaseDB'; filename = [matlabDir '/' filename '.mat']; try SCaseDB = SCase; save(filename, 'SCaseDB'); fprintf(logFileID, ': OK'); catch fprintf(logFileID, ': Error'); end %% Write csv file % This might be a function (SCase.csvSave()). The input would be a filename and a structure % data fprintf(logFileID, '\n Save csv database'); % 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'... ]; csvFilename = [xlsDir, '/anatomy.csv']; % Name of the csv file csvFileId = fopen(csvFilename, 'w'); fprintf(csvFileId,dataHeader); % Write header for iSCaseID = 1:length(SCaseList) % Loop on SCaseList to write data in csv disp(SCase(iSCaseID).id); % For debug fprintf(csvFileId,['\n'... SCase(iSCaseID).id ','... % SCase_id SCase(iSCaseID).shoulder.side ','... % shoulder_side num2str(SCase(iSCaseID).shoulder.scapula.glenoid.radius) ','... % glenoid_radius num2str(SCase(iSCaseID).shoulder.scapula.glenoid.sphereRMSE) ','... % glenoid_sphereRMSE num2str(SCase(iSCaseID).shoulder.scapula.glenoid.depth) ','... % glenoid_depth num2str(SCase(iSCaseID).shoulder.scapula.glenoid.width) ','... % glenoid_width num2str(SCase(iSCaseID).shoulder.scapula.glenoid.height) ','... % glenoid_height num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS num2str(SCase(iSCaseID).shoulder.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML num2str(SCase(iSCaseID).shoulder.scapula.glenoid.versionAmpl) ','... % glenoid_versionAmpl num2str(SCase(iSCaseID).shoulder.scapula.glenoid.versionOrient) ','... % glenoid_versionOrient num2str(SCase(iSCaseID).shoulder.scapula.glenoid.version) ','... % glenoid_version num2str(SCase(iSCaseID).shoulder.scapula.glenoid.inclination) ','... % glenoid_inclination num2str(SCase(iSCaseID).shoulder.humerus.jointRadius) ','... % humerus_jointRadius num2str(SCase(iSCaseID).shoulder.humerus.radius) ','... % humerus_headRadius num2str(SCase(iSCaseID).shoulder.humerus.GHSAmpl) ','... % humerus_GHSAmpl num2str(SCase(iSCaseID).shoulder.humerus.GHSOrient) ','... % humerus_GHSOrient num2str(SCase(iSCaseID).shoulder.humerus.SHSAmpl) ','... % humerus_SHSAmpl num2str(SCase(iSCaseID).shoulder.humerus.SHSOrient) ','... % humerus_SHSOrient num2str(SCase(iSCaseID).shoulder.humerus.SHSAngle) ','... % humerus_SHSAgle num2str(SCase(iSCaseID).shoulder.humerus.SHSPA) ','... % humerus_SHSPA num2str(SCase(iSCaseID).shoulder.humerus.SHSIS) ','... % humerus_SHSIS num2str(SCase(iSCaseID).shoulder.scapula.acromion.AI) ','... % acromion_AI num2str(SCase(iSCaseID).shoulder.scapula.acromion.CSA) ','... % acromion_CSA num2str(SCase(iSCaseID).shoulder.scapula.acromion.PSA) ','... % acromion_PSA num2str(SCase(iSCaseID).shoulder.scapula.acromion.AAA)... % acromion_AAA ]); end fclose(csvFileId); 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 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 new file mode 100644 index 0000000..a02eef6 --- /dev/null +++ b/openConfigFile.m @@ -0,0 +1,40 @@ +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. +% 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 +% TODO: Less restrict rules for writting config.txt + +if ~exist(configFile, 'file') + error([configFile, ' required']); +end +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=' +dataDir = dataDir(k:end); % Assume that config file ends with dataDir content +% Check if dataDir exists +if ~exist(dataDir, 'dir') + fprintf(logFileID, ['Data directory not found, check ' configFile]); + error(['Data directory not found, check ' configFile]); +end + +end + diff --git a/openLogFile.m b/openLogFile.m new file mode 100644 index 0000000..662408e --- /dev/null +++ b/openLogFile.m @@ -0,0 +1,20 @@ +function logFileID = openLogFile(filename) +%SETLOGFILE Creates a log file using filename in a log directory, within +%current working directory. +% The directory log is created if not present. + +% Check if log dir exists in working directory +logDir = 'log'; +if ~exist('log', 'dir') % If not create dir log + mkdir('log') +end +logFile = [logDir '/' filename]; +logFileID = fopen(logFile, 'w'); +fprintf(logFileID, filename); % Write name of file +fprintf(logFileID, '\nDate: '); % Write date +fprintf(logFileID, datestr(datetime)); % Write time + +output = 1; +message = 'OK'; +end + diff --git a/plotSCase.m b/plotSCase.m index 65e8185..423e18e 100755 --- a/plotSCase.m +++ b/plotSCase.m @@ -1,81 +1,85 @@ function SCase = plotSCase(SCaseId) %PLOTSCASE Plot a SCase. % Input: id char of shoulder case (e.g. 'P315') % % Output: Corresponding ShoulderCase object % % Example: SCase = plotSCase('P315'); % % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 -% Revision date: 2018-08-13 +% Revision date: 2018-12-30 % % TODO % Chech that the SCase exist and can be ploted -% Initialisation of paths -% Should be replaced by a function -dataDir = '../../../data'; % location of data +%% Open log file +logFileID = openLogFile('plotSCase.log'); -% Add path of ShoulderCase classes +%% Set the data directory from the configuration file config.txt +dataDir = openConfigFile('config.txt', logFileID); + +%% 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 % 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; + % 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 + fprintf(logFileID, ['\n' SCaseId]); else error('Not matlab file for this SCase'); end +fclose(logFileID); % Close log file end