diff --git a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m index 7efa616..5d3bc0b 100644 --- a/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m +++ b/ShoulderCase/@MuscleMeasurer/MuscleMeasurer.m @@ -1,119 +1,119 @@ classdef MuscleMeasurer < handle % Used to perfom PCSA and degeneration measurements of a Muscle object. % Slices and mask images are expected to be found at some specific places % and have specific names. % % This class is the results of extracting measurements methods from the % project of Nathan Donini. properties muscle segmentedArea segmentedImage rangeHU = double([-1000 1000]); muscleThreshold = 0; fatThreshold = 30; osteochondromaThreshold = 166; end methods function obj = MuscleMeasurer(muscle) obj.muscle = muscle; obj.segmentedArea = obj.getSegmentedArea; obj.segmentedImage = obj.getSegmentedImage; obj.normalizeThresholds; end function output = getSegmentedArea(obj) - output = logical(imread(fullfile(obj.muscle.maskDataPath,[obj.muscle.segmentationSet '_Mask.png']))); + output = logical(imread(fullfile(obj.muscle.maskDataPath,[obj.muscle.segmentationSet '.png']))); end function output = getSegmentedImage(obj) load(fullfile(obj.muscle.container.slicesDataPath,[obj.muscle.sliceName '_ForMeasurements.mat'])); normalizedImage = obj.getNormalizedImage(imageForMeasurements); segmentedImage = normalizedImage.*obj.segmentedArea; output = segmentedImage; end function output = getNormalizedImage(obj,image) image(image < obj.rangeHU(1)) = obj.rangeHU(1); image(image > obj.rangeHU(2)) = obj.rangeHU(2); output = (double(image)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1)); end function normalizeThresholds(obj) obj.muscleThreshold = ((obj.muscleThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1)); obj.fatThreshold = ((obj.fatThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1)); obj.osteochondromaThreshold = ((obj.osteochondromaThreshold - 1)-obj.rangeHU(1))/(obj.rangeHU(2)-obj.rangeHU(1)); end function output = getPCSA(obj) load(fullfile(obj.muscle.container.slicesDataPath,[obj.muscle.sliceName '_PixelSpacings.mat'])); pixelSurface = (imagesPixelSpacings(1) * imagesPixelSpacings(2)) / 100; % cm^2 output = pixelSurface * bwarea(obj.segmentedArea); end function output = getRatioAtrophy(obj) output = obj.getRatioAreas(obj.getAreaAtrophy,obj.segmentedArea); end function output = getRatioFat(obj) output = obj.getRatioAreas(obj.getAreaFat,obj.segmentedArea); end function output = getRatioOsteochondroma(obj) output = obj.getRatioAreas(obj.getAreaOsteochondroma,obj.segmentedArea); end function output = getRatioDegeneration(obj) output = obj.getRatioAtrophy + obj.getRatioFat + obj.getRatioOsteochondroma; end function output = getRatioAreas(obj,partialArea,totalArea) output = bwarea(partialArea)/bwarea(totalArea); end function output = getAreaMuscle(obj) % The area "Muscle" is the area of the segmented image that is not atrophied. % Looking for a better name. areaMuscle = imbinarize(obj.segmentedImage,obj.muscleThreshold); areaMuscle = imfill(areaMuscle,'holes'); % Keep biggest island only islands = bwconncomp(areaMuscle); islandsSizes = cellfun(@numel,islands.PixelIdxList); [~,biggestIslandIndex] = max(islandsSizes); areaMuscle = false(size(areaMuscle)); areaMuscle(islands.PixelIdxList{biggestIslandIndex}) = true; output = areaMuscle; end function output = getAreaAtrophy(obj) output = obj.segmentedArea & not(obj.getAreaMuscle); end function output = getAreaFat(obj) muscleImage = obj.segmentedImage.*obj.getAreaMuscle; areaFat = obj.getAreaMuscle & not(imbinarize(muscleImage,obj.fatThreshold)); output = areaFat; end function output = getAreaOsteochondroma(obj) muscleImage = obj.segmentedImage.*obj.getAreaMuscle; areaOsteochondroma = imbinarize(muscleImage,obj.osteochondromaThreshold); areaOsteochondroma = imfill(areaOsteochondroma,'holes'); output = areaOsteochondroma; end end end diff --git a/ShoulderCase/@MusclesContainer/plot.m b/ShoulderCase/@MusclesContainer/plot.m index 9872881..5d0407f 100644 --- a/ShoulderCase/@MusclesContainer/plot.m +++ b/ShoulderCase/@MusclesContainer/plot.m @@ -1,30 +1,36 @@ -function output = plot(obj) +function output = plot(obj,varargin) % For now this function is specific to plotting the rotator cuff segmentation % results + if nargin == 2 + maskName = [varargin{1} '.png']; + else + maskName = 'auto_Mask.png'; + end + SCase = obj.shoulder.SCase; muscleNames = fields(obj.list); rotatorCuffCrossSection = imread(fullfile(obj.slicesDataPath,... [obj.list.(muscleNames{1}).sliceName '_ForSegmentation.png'])); leg = {}; plotHandle(1) = imshow(rotatorCuffCrossSection); hold on colors = {'g','r','b','y'}; for i = 1:length(muscleNames) try - segmentedImage = imread(fullfile(obj.list.(muscleNames{i}).maskDataPath,'auto_Mask.png')); + segmentedImage = imread(fullfile(obj.list.(muscleNames{i}).maskDataPath,maskName)); catch continue end if (max(segmentedImage,[],'all') > 0) % if segmentation result is not empty leg{end+1} = muscleNames{i}; [~,plotHandle(i+1)] = contour(segmentedImage,'color',colors{i}); end end plotHandle(i+2) = legend(leg); output = plotHandle; end diff --git a/ShoulderCase/@Scapula/Scapula.m b/ShoulderCase/@Scapula/Scapula.m index 005bff9..51886c7 100644 --- a/ShoulderCase/@Scapula/Scapula.m +++ b/ShoulderCase/@Scapula/Scapula.m @@ -1,121 +1,126 @@ classdef (Abstract) Scapula < handle %SCAPULA Summary of this class goes here % This class defines the scapula. Landmarks are used to define its % coordinate system. It includes the glenoid object. properties angulusInferior % Landmark Angulus Inferior read from amira angulusInferiorLocal % Same as above, expressed the scapular coordinate system rather than the CT coordinate system trigonumSpinae % Landmark Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine read from amira processusCoracoideus % Landmark Processus Coracoideus (most lateral point) read from amira acromioClavicular % Landmark Acromio-clavicular joint (most lateral part on acromion)read from amira angulusAcromialis % Landmark Angulus Acromialis (most laterodorsal point) read from amira spinoGlenoidNotch % Landmark Spino-glenoid notch read from amira pillar % 5 landmarks on the pillar groove % 5 landmarks on the scapula (supraspinatus) groove coordSys % Scapular coordinate system plane % Plane class segmentation % either none 'N', manual 'M' or automatic 'A' surface % surface points and triangles of the scapula glenoid % Glenoid class acromion % Acromion class comment % any comment end properties (Hidden = true) shoulder end methods (Abstract) load(obj); createGlenoid(obj); end methods (Access = protected) function obj = Scapula(shoulder) %SCAPULA Construct an instance of this class % Constructor of the scapula object, sets all properties to % zero. obj.shoulder = shoulder; obj.angulusInferior = []; obj.angulusInferiorLocal = []; obj.trigonumSpinae = []; obj.processusCoracoideus = []; obj.acromioClavicular = []; obj.angulusAcromialis = []; obj.spinoGlenoidNotch = []; obj.pillar = []; obj.groove = []; obj.coordSys = CoordinateSystemScapula(); obj.plane = Plane(); obj.segmentation = 'N'; obj.surface = []; obj.createGlenoid(); obj.acromion = Acromion(obj); obj.acromion.scapula = obj; obj.comment = ''; end end methods function output = coordSysSet(obj) try obj.setCoordinateSystemWithLandmarks(); catch ME warning(ME.message); output = 0; end % To enable retro compatibility with scapula.coordSysSet former % implementation, the two following commands are called here obj.setShoulderSideWithCoordinateSystem; obj.angulusInferiorLocal = obj.coordSys.express(obj.angulusInferior); output = 1; end function measurePlane(obj) % Scapular plane is fitted on 3 points (angulusInferior, % trigonumSpinae, most laretal scapular goove landmark). inferior = obj.angulusInferior; medial = obj.trigonumSpinae; mostLateralGrooveIndex = findLongest3DVector(medial-obj.groove); mostLateralGroovePoint = obj.groove(mostLateralGrooveIndex(1),:); + + % trying alternative: the lateral point is spinoGlenoidNotch + mostLateralGroovePoint = obj.spinoGlenoidNotch; + % trying alternative: the lateral point is the mean of groove points + % mostLateralGroovePoint = mean(obj.groove); obj.plane.fit([inferior; medial; mostLateralGroovePoint]); anterior = obj.processusCoracoideus; posterior = obj.angulusAcromialis; obj.plane.normal = orientVectorToward(obj.plane.normal,(anterior-posterior)); end function setShoulderSideWithCoordinateSystem(obj) % The scapula coordinate system is right-handed for the right shoulder only if obj.coordSys.isRightHanded obj.shoulder.side = 'R'; return elseif obj.coordSys.isLeftHanded obj.shoulder.side = 'L'; return end warning(['Couldn''t find the shoulder side by evaluating the scapula',... 'coordinate system.']); end end end diff --git a/ShoulderCase/@Scapula/setCoordinateSystemWithLandmarks.m b/ShoulderCase/@Scapula/setCoordinateSystemWithLandmarks.m index 68ccdff..d1e4723 100644 --- a/ShoulderCase/@Scapula/setCoordinateSystemWithLandmarks.m +++ b/ShoulderCase/@Scapula/setCoordinateSystemWithLandmarks.m @@ -1,56 +1,60 @@ function setCoordinateSystemWithLandmarks(obj) % Calculate the (EPFL) scapular coordinate system % The EPFL scapular coordinate system is defined for a right scapula see % paper (10.1302/0301-620X.96B4.32641). The X axis is % antero-posterior, the Y axis is infero-superior, the Z axis % is meddio-lateral. This system is right-handed. % This orientation corresponds approximatively to the ISB % recommendadtion (doi:10.1016/j.jbiomech.2004.05.042). % For left scapula we keep the x axis as postero-anterior and % get a left-handed coordinate system. if (length(obj.groove) < 2) error(['Can''t set scapular coordinate system with actual',... ' obj.groove']); return; end % The four following methods have to be called in this specific order setPAAxis(obj); setMLAxis(obj); setISAxis(obj); setOrigin(obj); end function setPAAxis(obj) % Set the postero-anterior axis to be the scapula plane normal obj.coordSys.PA = obj.plane.normal; end function setMLAxis(obj) % Set the media-lateral axis to be the line fitted to the projection of % the groove points on the scapula plane lateral = obj.spinoGlenoidNotch; medial = obj.trigonumSpinae; groovePointsProjection = obj.plane.projectOnPlane(obj.groove); - grooveAxis = fitLine(groovePointsProjection); + + % trying alternative: add spinoGlenoidNotch and trigonumSpinae + grooveAxis = fitLine([mean(groovePointsProjection);lateral;medial]); grooveAxis = (grooveAxis/norm(grooveAxis))'; obj.coordSys.ML = orientVectorToward(grooveAxis,(lateral-medial)); end function setISAxis(obj) % Set the infero-superior axis to be orthogonal with the two other axes superior = obj.spinoGlenoidNotch; inferior = obj.angulusInferior; obj.coordSys.IS = cross(obj.coordSys.PA,obj.coordSys.ML); obj.coordSys.IS = orientVectorToward(obj.coordSys.IS,(superior-inferior)); end function setOrigin(obj) % Set the origin of the coordinate system to the spino-glenoid notch % projected on the scapular axis spinoGlenoid = obj.spinoGlenoidNotch; - grooveMeanProjection = mean(obj.plane.projectOnPlane(obj.groove)); + + % trying alternative: add spinoGelnoidNotch and trigonumSpinae + grooveMeanProjection = mean(obj.plane.projectOnPlane([mean(obj.groove);obj.trigonumSpinae;obj.spinoGlenoidNotch])); obj.coordSys.origin = grooveMeanProjection + ... dot((spinoGlenoid - grooveMeanProjection),obj.coordSys.ML)*obj.coordSys.ML; end diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m index 0b82305..0a6baea 100755 --- a/ShoulderCase/@ShoulderCase/ShoulderCase.m +++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m @@ -1,212 +1,212 @@ classdef ShoulderCase < handle % Properties and methods associated to the SoulderCase object. % Author: Alexandre Terrier, EPFL-LBO % Matthieu Boubat, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2020-07-23 % TODO: % Need method to load properties (from amira, matlab, excel, SQL) % Need method to save properties (in matlab, excel, SQL) % Remove datapath from constants % properties id = []; % id of the shoulder case, as it appears in the database (Pnnn) diagnosis = []; treatment = []; outcome = []; patient = []; shoulderManual = []; shoulderAuto = []; study = []; comment = []; end % properties (Constant, Hidden = true) % dataPath = '/Volumes/shoulder/dataDev'; % path to data folder from package folder % % This should be done differently. Check where we are, and where we % % should go. % end properties (Hidden = true) dataPath = []; dataCTPath = []; dataAmiraPath = []; dataMatlabPath = []; dataDicomPath = []; id4C = []; % SCaseId with P/N followed by 3 digits --> 4 char end - methods (Access = ?ShoulderCaseLoader) + methods (Access = public)%?ShoulderCaseLoader) function obj = ShoulderCase(SCaseID,dataCTPath) % SCaseID validation rawSCaseID = SCaseIDParser(SCaseID); assert(rawSCaseID.isValidID,'The input argument is not a valid SCaseID.') obj.id = SCaseID; obj.id4C = rawSCaseID.getIDWithNumberOfDigits(4); % path attribution obj.dataCTPath = dataCTPath; % Initialsation obj.patient = Patient(obj); obj.shoulderManual = ShoulderManual(obj); obj.shoulderAuto = ShoulderAuto(obj); obj.propagateDataPath; end end methods function propagateDataPath(obj) % Update current dataPath % Propagate the path to objects in properties obj.dataAmiraPath = fullfile(obj.dataCTPath,'amira'); obj.dataMatlabPath = fullfile(obj.dataCTPath,'matlab'); obj.dataDicomPath = fullfile(obj.dataCTPath,'dicom'); obj.shoulderManual.propagateDataPath; obj.shoulderAuto.propagateDataPath; 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.shoulderManual.scapula.glenoid.radius; % 2 Result.glenoid_radiusRMSE = obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE; % 3 % 3 %Result.glenoidSphericity = ' '; %Result.glenoid_biconcave = ' '; Result.glenoid_depth = obj.shoulderManual.scapula.glenoid.depth; % 4 Result.glenoid_width = obj.shoulderManual.scapula.glenoid.width; % 5 Result.glenoid_height = obj.shoulderManual.scapula.glenoid.height; % 6 Result.glenoid_center_PA = obj.shoulderManual.scapula.glenoid.centerLocal(1); % 7 Result.glenoid_center_IS = obj.shoulderManual.scapula.glenoid.centerLocal(2); % 8 Result.glenoid_center_ML = obj.shoulderManual.scapula.glenoid.centerLocal(3); % 9 Result.glenoid_version_ampl = obj.shoulderManual.scapula.glenoid.versionAmplitude; % 10 Result.glenoid_version_orient = obj.shoulderManual.scapula.glenoid.versionOrientation; % 11 Result.glenoid_Version = obj.shoulderManual.scapula.glenoid.version; % 12 Result.glenoid_Inclination = obj.shoulderManual.scapula.glenoid.inclination; % 13 Result.humerus_joint_radius = ' '; % 14 Result.humeral_head_radius = obj.shoulderManual.humerus.radius; % 15 % Result.humerus_GHsublux_2D = ' '; % Result.humerus_SHsublux_2D = ' '; Result.humerus_GHsubluxation_ampl = obj.shoulderManual.humerus.GHSAmpl; % 16 Result.humerus_GHsubluxation_orient = obj.shoulderManual.humerus.GHSOrient; % 17 Result.humerus_SHsubluxation_ampl = obj.shoulderManual.humerus.SHSAmpl; % 18 Result.humerus_SHsubluxation_orient = obj.shoulderManual.humerus.SHSOrient; % 19 Result.scapula_CSA = obj.shoulderManual.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.shoulderManual.scapula.acromion.acromionIndex; outputArg = Result; end function outputArg = saveMatlab(obj) % Save SCase to matlab file dir = obj.dataMatlabPath; % Create dir if not exist try if ~exist(dir, 'dir') % create directory if it does not exist mkdir(dir); end catch fprintf('Error creating the matlab directory \n'); % Should be in log end % Save SCase in matlab directoty, in a file named SCaseCNNN.m filename = 'SCase'; filename = [dir '/' filename '.mat']; try SCase = obj; save(filename, 'SCase'); outputArg = 1; catch fprintf('Error creating SCase matlab file \n'); % Should be in log outputArg = -1; end end function outputArg = appendToCSV(obj,filename) % Save SCase to csv file logFid = fopen('log/measureSCase.log', 'a'); dataDir = ConfigFileExtractor.getVariable('dataDir'); xlsDir = [dataDir '/Excel/xlsFromMatlab']; fid = fopen([xlsDir '/' filename],'a'); fprintf(fid,[... obj.id ','... % SCase_id obj.shoulderManual.side ','... % shoulder_side num2str(obj.shoulderManual.scapula.glenoid.radius) ','... % glenoid_radius num2str(obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE) ','... % glenoid_sphereRMSE num2str(obj.shoulderManual.scapula.glenoid.depth) ','... % glenoid_depth num2str(obj.shoulderManual.scapula.glenoid.width) ','... % glenoid_width num2str(obj.shoulderManual.scapula.glenoid.height) ','... % glenoid_height num2str(obj.shoulderManual.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA num2str(obj.shoulderManual.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS num2str(obj.shoulderManual.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML num2str(obj.shoulderManual.scapula.glenoid.versionAmplitude) ','... % glenoid_versionAmpl num2str(obj.shoulderManual.scapula.glenoid.versionOrientation) ','... % glenoid_versionOrientation num2str(obj.shoulderManual.scapula.glenoid.version) ','... % glenoid_version num2str(obj.shoulderManual.scapula.glenoid.inclination) ','... % glenoid_inclination num2str(obj.shoulderManual.humerus.jointRadius) ','... % humerus_jointRadius num2str(obj.shoulderManual.humerus.radius) ','... % humerus_headRadius num2str(obj.shoulderManual.humerus.GHSAmpl) ','... % humerus_GHSAmpl num2str(obj.shoulderManual.humerus.GHSOrient) ','... % humerus_GHSOrient num2str(obj.shoulderManual.humerus.SHSAmpl) ','... % humerus_SHSAmpl num2str(obj.shoulderManual.humerus.SHSOrient) ','... % humerus_SHSOrient num2str(obj.shoulderManual.humerus.SHSAngle) ','... % humerus_SHSAgle num2str(obj.shoulderManual.humerus.SHSPA) ','... % humerus_SHSPA num2str(obj.shoulderManual.humerus.SHSIS) ','... % humerus_SHSIS num2str(obj.shoulderManual.scapula.acromion.AI) ','... % acromion_AI num2str(obj.shoulderManual.scapula.acromion.CSA) ','... % acromion_CSA num2str(obj.shoulderManual.scapula.acromion.PSA) ','... % acromion_PSA num2str(obj.shoulderManual.scapula.acromion.AAA) '\n'... % acromion_AAA ]); fclose(fid); fclose(logFid); outputArg = 1; end function outputArg = saveExcel(~) % Save SCase to Excel file outputArg = 1; end function outputArg = saveSQL(~) % Save SCase to MySQL database outputArg = 1; end end end diff --git a/measureSCase.m b/measureSCase.m index 48feb62..a8b0d81 100755 --- a/measureSCase.m +++ b/measureSCase.m @@ -1,674 +1,674 @@ 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.getAllCasesIDs(); 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(contains(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'); try shoulder.muscles.createRotatorCuffSlices(); shoulder.muscles.createRotatorCuffSlicesOld(); shoulder.muscles.segmentRotatorCuffMuscles('rotatorCuff'); - shoulder.muscles.measureAllMuscles('auto'); + shoulder.muscles.measureAllMuscles('auto_Mask'); catch ME fprintf(logFileID,ME.message); end fprintf(logFileID, ': OK'); end