Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/measureSCase.m
e504ac2a0db1aterrier-plotUpdate

Last Change
Raw File

measureSCase.m
View Options

function [status,message] = measureSCase(varargin)
% MEASURESCASE Update the entire SCase DB with anatomical measurements.
% The function should be run from the server (lbovenus), from the user account containing the git repository,
% with the following shell command
% matlab -nodisplay -nosplash -r "measureSCase;quit"
% For long calculations, it can be run in the background using tmux, followed
% by Ctrl+b d to detach the session and keep the process running.
% Progress can be checked through log file with following shell command
% 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 variable
% Input arguments can be one or several among:
% {Empty,'glenoidDensity', 'update', 'musclesDegeneration' '[NP][1-999]', '[NP]'}
%
% If there is no argument the function will measure every new case.
% The 'N' or 'P' given as an argument will load all the Normal or Pathological cases.
%
% The measurements are divided into three measurements: "morphology", "glenoidDensity",
% and "musclesDegeneration". The default measurement is the "morphology" one, the
% two others has to be given in arguments in order to be done.%
%
% Examples:
% measureSCase()
% is a valid instruction where on every new case is done the
% "morphology" measurement.
% measureSCase('glenoidDensity','musclesDegeneration')
% is a valid instruction where on every cases are done the
% "morphology","glenoidDensity", and "musclesDegeneration" measurements
% if one of these measurements is missing.
% measureSCase('P400','N520')
% is a valid instruction where on the cases 'P400' and 'N520' is
% done the "morphology" measurement only if they have not already
% been measured yet (new cases).
% measureSCase('P400','N520','update')
% is a valid instruction where on the cases 'P400' and 'N520' is
% done the "morphology" measurement whether they have already been
% measured yet or not.
% measureSCase('glenoidDensity','P400','N520','update','P')
% is a valid instruction where on all the 'P' cases and the case N520
% are done the "morphology" and "glenoidDensity" measurements whether
% they have been measured yet or not.
% 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
% Authors: Alexandre Terrier, EPFL-LBO
% Matthieu Boubat, EPFL-LBO
% Creation date: 2018-07-01
% Revision date: 2020-06-09
% 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
startTime = datetime('now');
%% Add environment to path
addpath(genpath(pwd));
%% Open log file
logFileID = openLogFile('measureSCase.log');
%% Set the data directory from the configuration file config.txt
dataDir = ConfigFileExtractor.getVariable('dataDir');
%% Location of the XLS ShoulderDatabase
xlsDir = [dataDir '/Excel/xlsFromMatlab'];
matlabDir = [dataDir '/matlab'];
%% Instanciate case loader
database = ShoulderCaseLoader();
%% Get list of all SCase
% Get list of SCases from varargin
% Delault value for varargin
calcGlenoidDensity = false;
calcMusclesDegeneration = false;
updateCalculations = false;
specificCases = true;
fprintf(logFileID, '\n\nList of SCase');
allCases = database.getAllCasesID();
assert(not(isempty(allCases)),'The database is empty, check that the provided data path is valid.');
requestedCases = {};
for argument = 1:nargin
switch lower(varargin{argument}) % Test every argument in varargin.
case 'glenoiddensity'
calcGlenoidDensity = true;
case 'musclesdegeneration'
calcMusclesDegeneration = true;
case 'update'
updateCalculations = true;
case regexp(lower(varargin{argument}),'[p]','match') % Detect if argument is 'P'.
requestedCases = [requestedCases database.getAllPathologicalCasesIDs()];
case regexp(lower(varargin{argument}),'[n]','match') % Detect if argument is 'N'.
requestedCases = [requestedCases database.getAllNormalCasesIDs()];
case regexp(lower(varargin{argument}),'[np][1-9][0-9]{0,2}','match') % Test for a correct argument [NP][1-999]
if database.containsCase(varargin{argument})
requestedCases = [requestedCases varargin{argument}];
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 "[npNP]",'...
' "[npNP][1-999]", "GlenoidDensity", or "update"'],varargin{argument})
end
end
SCaseList = allCases(find(ismember(allCases,requestedCases))); % Remove cases that appears multiple times
if isempty(SCaseList)
disp('No valid ShoulderCase have been given in argument. All the cases found in the database will be measured.');
SCaseList = allCases;
end
%% Load current xls database (for clinical data)
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
SCaseDB = {};
SCase = [];
% %%%% definition of analysis variables
% %
% % % manual
% %
% ok = {};
% amira = {};
% 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};
if not(updateCalculations)
try
SCase = database.loadCase(SCaseID);
catch ME
SCase = database.createEmptyCase(SCaseID);
warning(ME.message);
end
else
SCase = database.createEmptyCase(SCaseID);
end
% newSCase = ShoulderCase(SCaseID);
% SCase = newSCase;
% SCase.dataPath = dataDir; % Set dataDir for SCase
manualMeasurementAvailable = not(isempty(dir(SCase.dataAmiraPath))); % Set data path of this SCase
autoMeasurementsAvailable = not(isempty(dir([SCase.dataMatlabPath '/scapulaLandmarksAuto*.mat'])));
percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f');
fprintf(logFileID, ['\n___________________________________\n|\n| SCaseID: ' SCaseID ' (' percentProgress '%%)']);
%% The following section is a temporary solution to the arguments handling.
%% It splits the measurements process into six parts (shoulderMorphology,
%% glenoidDensity, musclesDegeneration for both manual and auto cases).
%% Each part are independently 'ordered' based on some conditions.
orderMorphologyMeasurementManual = false;
orderDensityMeasurementManual = false;
orderDegenerationMeasurementManual = false;
orderMorphologyMeasurementAuto = false;
orderDensityMeasurementAuto = false;
orderDegenerationMeasurementAuto = false;
if or(not(exist([SCase.dataMatlabPath '/SCase.mat'],'file')), updateCalculations)
orderMorphologyMeasurementManual = true;
orderMorphologyMeasurementAuto = true;
if (calcGlenoidDensity)
orderDensityMeasurementManual = true;
orderDensityMeasurementAuto = true;
end
if (calcMusclesDegeneration)
orderDegenerationMeasurementManual = true;
orderDegenerationMeasurementAuto = true;
end
else
loadedSCase = database.loadCase(SCase.id);
musclesManual = loadedSCase.shoulderManual.muscles;
musclesAuto = loadedSCase.shoulderAuto.muscles;
if and(calcGlenoidDensity, isempty(loadedSCase.shoulderManual.scapula.glenoid.density))
orderMorphologyMeasurementManual = true;
orderDensityMeasurementManual = true;
SCase = loadedSCase;
end
if and(calcGlenoidDensity, isempty(loadedSCase.shoulderAuto.scapula.glenoid.density))
orderMorphologyMeasurementAuto = true;
orderDensityMeasurementAuto = true;
SCase = loadedSCase;
end
if and(calcMusclesDegeneration, not(all(table2array(musclesManual.summary(:,{'PCSA'})))) )
orderMorphologyMeasurementManual = true;
orderDegenerationMeasurementManual = true;
SCase = loadedSCase;
end
if and(calcMusclesDegeneration, not(all(table2array(musclesAuto.summary(:,{'PCSA'})))) )
orderMorphologyMeasurementAuto = true;
orderDegenerationMeasurementAuto = true;
SCase = loadedSCase;
end
end
%%
%% The section above (and basically all the measureSCase function) will be removed
%% in the upcoming implementation of ShoulderCase measurements.
%% Matthieu Boubat, 06.03.2020
%%
% 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 manualMeasurementAvailable
tic;
fprintf(logFileID, '\n| Segmentation manual ');
try
if orderMorphologyMeasurementManual
measureShoulderMorphology(SCase.shoulderManual,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
try
if orderDensityMeasurementManual
measureGlenoidDensity(SCase.shoulderManual,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
try
if orderDegenerationMeasurementManual
measureMusclesDegeneration(SCase.shoulderManual,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
fprintf(logFileID, '\n| Elapsed time (s): %s\n|', num2str(toc));
end
% 2) Load & analyses auto data from matlab
if autoMeasurementsAvailable
tic;
fprintf(logFileID, '\n| Segmentation auto ');
try
if orderMorphologyMeasurementAuto
measureShoulderMorphology(SCase.shoulderAuto,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
try
if orderDensityMeasurementAuto
measureGlenoidDensity(SCase.shoulderAuto,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
try
if orderDegenerationMeasurementAuto
measureMusclesDegeneration(SCase.shoulderAuto,logFileID);
end
catch ME
fprintf(logFileID,ME.message);
end
fprintf(logFileID, '\n| Elapsed time (s): %s\n|', num2str(toc));
end
% 3) Set clinical data from Excel database to SCase
fprintf(logFileID, '\n|\n| Set clinical data from Excel');
% Get idx of excelRaw for SCaseID
% Use cell2struct when dots are removed from excel headers
try
idx = find(strcmp(excelRaw, SCaseID)); % Index of SCaseID in excelRaw
SCase.diagnosis = excelRaw{idx,4};
SCase.treatment = excelRaw{idx,9};
SCase.patient.gender = excelRaw{idx,19};
SCase.patient.age = excelRaw{idx,21};
SCase.patient.height = excelRaw{idx,23};
SCase.patient.weight = excelRaw{idx,24};
SCase.patient.BMI = excelRaw{idx,25};
SCase.shoulderAuto.scapula.glenoid.walch = excelRaw{idx,55};
SCase.shoulderManual.scapula.glenoid.walch = excelRaw{idx,55};
fprintf(logFileID, ': OK');
catch
fprintf(logFileID, ': ERROR');
end
% 4) Save SCase
try
fprintf(logFileID, '\n| Save SCase.mat');
SCase.saveMatlab;
fprintf(logFileID, ': OK');
catch
fprintf(logFileID, ': ERROR');
end
fprintf(logFileID, '\n|___________________________________\n\n');
% 4) Load & analyses auto data from matlab (Statistical Shape Model)
% Load scapula surface and landmarks
% save('measureSCase_analysis.mat',...
% 'ok','amira','scapulaLoad','scapulaCoord',...
% 'degeneration','glenoidLoad','glenoidMorphology','glenoidDensity',...
% 'humerusLoad','humerusSubluxation','acromionMorphology',...
% 'autoOk','autoScapulaLoad','autoScapulaCoord','autoGlenoidLoad',...
% 'autoGlenoidMorphology','autoGlenoidDensity','autoAcromionMorphology')
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'...
];
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
SCaseDB = database.loadAllCases();
cellfun(@(SCase)SCase.appendToCSV('anatomy.csv'),SCaseDB); % Call method appendToCSV() on all element of the SCaseDB cell array
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
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
function measureShoulderMorphology(shoulder,logFileID)
fprintf(logFileID, '\n| Load scapula surface and landmarks');
output = shoulder.scapula.load;
if ~output
fprintf(logFileID, ': ERROR');
% scapulaLoad{end+1} = SCaseID;
return
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Set coord. syst.');
output = shoulder.scapula.coordSysSet;
if ~output
fprintf(logFileID, ': ERROR');
% scapulaCoord{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Load glenoid surface');
output = shoulder.scapula.glenoid.load;
if ~output
fprintf(logFileID, '\n| glenoid surface loading error');
% glenoidLoad{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Measure glenoid anatomy');
output = shoulder.scapula.glenoid.morphology;
if ~output
fprintf(logFileID, ': ERROR');
% glenoidMorphology{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Load humerus data');
output = shoulder.humerus.load;
if ~output
fprintf(logFileID, ': ERROR');
% humerusLoad{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Measure humerus subluxation');
output = shoulder.humerus.subluxation(shoulder.scapula);
if ~output
fprintf(logFileID, ': ERROR');
% humerusSubluxation{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
fprintf(logFileID, '\n| Measure acromion anatomy');
% Should be run after SCase.shoulderManual.humerus (for AI)
output = shoulder.scapula.acromion.morphology;
if ~output
fprintf(logFileID, ': ERROR');
% acromionMorphology{end+1} = SCaseID;
return;
end
fprintf(logFileID, ': OK');
end
function measureGlenoidDensity(shoulder,logFileID)
fprintf(logFileID, '\n| Measure glenoid density');
output = shoulder.scapula.glenoid.calcDensity;
if ~output
fprintf(logFileID, ': Density calculus error. Could not read dicom');
return;
end
fprintf(logFileID, ': OK');
end
function measureMusclesDegeneration(shoulder,logFileID)
fprintf(logFileID, '\n| Measure muscle''s degeneration');
% auto segmentation with most recent slicing methods
try
shoulder.muscles.createRotatorCuffSlicesMatthieu('rotatorCuffMatthieu');
shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuffMatthieu','autoMatthieu');
catch ME
fprintf(logFileID,ME.message);
end
% auto segmentation with former slicing methods
try
shoulder.muscles.createRotatorCuffSlicesNathan('rotatorCuffNathan');
shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuffNathan','autoNathan');
catch ME
fprintf(logFileID,ME.message);
end
% muscles measurements
try
shoulder.muscles.measureAllMuscles('rotatorCuffMatthieu','autoMatthieu');
catch ME
fprintf(logFileID,ME.message);
end
fprintf(logFileID, ': OK');
end
c4science · Help