diff --git a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m index aedc83c..3ed8664 100644 --- a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m +++ b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m @@ -1,38 +1,38 @@ classdef GlenoidAuto < Glenoid % To be used with ShoulderAuto data. % The load() method requires a specific implementation. methods function outputArg = load(obj) % LOADAUTO Load genoid surface from auto segmentation % Load the points of the glenoid surface from matlab dir SCase = obj.scapula.shoulder.SCase; SCaseId4C = SCase.id4C; matlabDir = SCase.dataMatlabPath; side = SCase.shoulderAuto.side; % Import glenoid surface points - fileName = ['glenoidSurfaceAuto' side '.ply']; - fileName = [matlabDir '/' fileName]; + fileName = "glenoidSurfaceAuto" + side + ".ply"; + fileName = fullfile(matlabDir, fileName); if exist(fileName,'file') == 2 [points,faces]=loadPly(fileName,1); obj.surface.points = points; obj.surface.meanPoint = mean(points); obj.surface.faces = faces; else outputArg = 0; return % error(['Could not find file: ' fileName]) end outputArg = 1; % Should report on loading result end end end diff --git a/ShoulderCase/@GlenoidManual/GlenoidManual.m b/ShoulderCase/@GlenoidManual/GlenoidManual.m index cc14619..45b3e7e 100644 --- a/ShoulderCase/@GlenoidManual/GlenoidManual.m +++ b/ShoulderCase/@GlenoidManual/GlenoidManual.m @@ -1,39 +1,39 @@ classdef GlenoidManual < Glenoid % To be used with ShoulderManual data. % The load() method requires a specific implementation. methods function outputArg = load(obj) % LOAD Summary of this method goes here % Load the points of the glenoid surface from the amira % directoty. SCase = obj.scapula.shoulder.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; % Import glenoid surface points - fileName = ['ExtractedSurface' SCaseId4C '.stl']; - fileName = [amiraDir '/' fileName]; + fileName = "ExtractedSurface" + SCaseId4C + ".stl"; + fileName = fullfile(amiraDir, fileName); if exist(fileName,'file') == 2 [points,faces,~]=loadStl(fileName,1); obj.surface.points = points; obj.surface.meanPoint = mean(points); obj.surface.faces = faces; else - warning(['Could not find file: ' fileName]); + warning("Could not find file: " + fileName); outputArg = 0; return; % error(['Could not find file: ' fileName]) end outputArg = 1; % Should report on loading result end end end diff --git a/ShoulderCase/@Humerus/Humerus.m b/ShoulderCase/@Humerus/Humerus.m index 906b9d8..337bd21 100644 --- a/ShoulderCase/@Humerus/Humerus.m +++ b/ShoulderCase/@Humerus/Humerus.m @@ -1,257 +1,225 @@ classdef Humerus < handle %HUMERUS Properties and methods associated to the humerus % Detailed explanation goes here % Author: Alexandre Terrier, EPFL-LBO % Creation date: 2018-07-01 % Revision date: 2019-06-29 % % TO DO: % Local coordinate system properties landmarks % 5 3D points center % Center of the humeral head (sphere fit on 5 points radius % Radius of the humeral head (sphere fit on 5 points jointRadius % Radius of cuvature of the articular surface (todo) SHSAngle % Scapulo-humeral subluxation angle SHSPA % Scapulo-humeral subluxation angle in the postero-anterior direction (posterior is negative, as for glenoid version) SHSIS % Scapulo-humeral subluxation angle in the infero-superior direction SHSAmpl % Scapulo-humeral subluxation (center ofset / radius) SHSOrient % Scapulo-humral subluxation orientation (0 degree is posterior) GHSAmpl % Gleno-humeral subluxation (center ofset / radius) GHSOrient % Gleno-humral subluxation orientation (0 degree is posterior) subluxationIndex3D shoulder end methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Humerus function obj = Humerus(shoulder) %HUMERUS Construct an instance of this class % Detailed explanation goes here obj.landmarks = []; obj.center = []; obj.radius = []; obj.jointRadius = []; obj.SHSAngle = []; obj.SHSPA = []; obj.SHSIS = []; obj.SHSAmpl = []; obj.SHSOrient = []; obj.GHSAmpl = []; obj.GHSOrient = []; obj.shoulder = shoulder; end end methods function outputArg = load(obj) - %LOAD Load 5 humeral head landmarks - % TODO: Should be check for 'old' vs 'new' file definition - SCase = obj.shoulder.SCase; - SCaseId4C = SCase.id4C; - amiraDir = SCase.dataAmiraPath; - - % Try 1st with new defiunition of HHLandmarks - fileName = "newHHLandmarks" + SCaseId4C + ".landmarkAscii"; - fileName = fullfile(amiraDir, fileName); - - if exist(fileName,'file') == 2 - importedData = importdata(fileName, ' ', 14); - else % Check for old version of HHLandmarks - fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii']; - fileName = [amiraDir '/' fileName]; - if exist(fileName,'file') == 2 - importedData = importdata(fileName, ' ', 14); - warning([SCaseId4C ' is using old humeral landmarks definition']) - else - outputArg = 0; - return - end + if exist(fullfile(SCase.dataSlicerPath, "HH_landmarks.mrk.json"), "file") + success = obj.loadSlicerLandmarks(); + elseif not(isempty(dir(fullfile(SCase.dataAmiraPath, "*HH*.landmarksAscii")))) + success = obj.loadAmiraLandmarks(); + else + success = 0; end - % with other methods - - % Check validity of data in file - % Variable points belos should be renamed landmarks, for consistency - try - points = importedData.data; - catch - % Cant get landarks from file - warning([SCaseId4C ' cant load humeral landmarks']) - outputArg = 0; - return; - end - % We should add a test for validity of points - - - obj.landmarks = points; - outputArg = 1; % Should report on loading result + outputArg = success; % Should report on loading result end function outputArg = subluxation(obj) % Calcul of humeral head subluxation % The function also get center and radius of the humeral head. % TODO: % Calcul would be simpler in the scapular cordinate system, as % done in Acromion object. % Radius of curvature of humeral head articular surface % Might be removed anatomy for concistency SCase = obj.shoulder.SCase; scapula = obj.shoulder.scapula; %% Get data from parent shoulder object % Shoulder side side = obj.shoulder.side; %% Get data from parent scapula object % Scapula coordinate system xAxis = scapula.coordSys.PA; yAxis = scapula.coordSys.IS; zAxis = scapula.coordSys.ML; %% Get data from parent glenoid object % Glenoid center glenCenter = scapula.glenoid.center; % Glenoid centerLine from parent glenoid object glenCenterline = scapula.glenoid.centerLine; %% Get data from humerus % humeral head landmarks (5 or ?) humHead = obj.landmarks; % Humeral head landmarks %% Anatomical measurements % Fit sphere on humeral head landmarks [humHeadCenter,humHeadRadius,HHResiduals, R2HH] = fitSphere(humHead); humHeadCenter = humHeadCenter'; % transform vertical to horizantal vector %% Scapulo-humeral subluxation (SHS) % Vector from glenoid center to humeral head center SHSvect = humHeadCenter - glenCenter; % Vector from humeral head center to scapular (z) axis, % perpendicular to z-axis --> similar to projection in xy plane % This might be rewritten more clearly SHSvectxy = SHSvect - (dot(SHSvect,zAxis)*zAxis); % SHS amplitude (ratio rather than amplitude) SHSAmpl = norm(SHSvectxy) / humHeadRadius; % SHS orientation xProj = dot(SHSvectxy, xAxis); yProj = dot(SHSvectxy, yAxis); if xProj > 0 % anterior side if yProj > 0 % ant-sup quadran SHSOrient = pi - atan(yProj/xProj); else % ant-inf quadran SHSOrient = -pi - atan(yProj/xProj); end elseif xProj < 0 % posterior side SHSOrient = - atan(yProj/xProj); else % xProj = 0 SHSOrient = pi/2 - pi * (yProj < 0); end SHSOrient = rad2deg(SHSOrient); % SHSAngle: angle between glenHumHead and zAxis SHSAngle = vrrotvec(SHSvect,zAxis); SHSAngle = SHSAngle(4); SHSAngle = rad2deg(SHSAngle); % SHSAP is angle between zxProjection and zAxis % Project glenHumHead on yAxis yProjVect = dot(SHSvectxy, yAxis)*yAxis; % Remove projection from glenHumHead --> zxProjection glenHumHeadzx = SHSvect - yProjVect; a = glenHumHeadzx; b = zAxis; rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b rotVect3 = rotVect4(1:3); rotAngle = rotVect4(4); rotSign = sign(dot(rotVect3,yAxis)); % Sign of rotation SHSPA = rotAngle * rotSign; SHSPA = -SHSPA; % Convention that posterior is negative if strcmp(side,'L') % Change signe for left-handed scapula coordinate system SHSPA = -SHSPA; end SHSPA = rad2deg(SHSPA); % Angle in degrees % SHSIS is angle between yzProjection and zAxis % Project glenHumHead on xAxis xProjVect = dot(SHSvectxy, xAxis)*xAxis; % Remove projection from glenHumHead --> yzProjection glenHumHeadyz = SHSvect - xProjVect; a = glenHumHeadyz; b = zAxis; rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b rotVect3 = rotVect4(1:3); rotAngle = rotVect4(4); rotSign = sign(dot(rotVect3,xAxis)); % Sign of rotation SHSIS = rotAngle * rotSign; SHSIS = SHSIS; % Convention that superior is positive if strcmp(side,'L') % Change signe for left-handed scapula coordinate system SHSIS = -SHSIS; end SHSIS = rad2deg(SHSIS); % Angle in degrees %% Angle between plane of maximum SH subluxation and CT axis [0,0,1] % Not used anymore % SHSPlane = cross(glenHumHead, zAxis); % G = vrrotvec(SHSPlane, [0 0 1]); % SHSPlaneAngle = rad2deg(G(4)); % if SHSPlaneAngle > 90 % SHSPlaneAngle = 180 - SHSPlaneAngle; % end %% Gleno-humeral subluxation (GHS) % Vector from glenoid sphere centre to glenoid centerline (perpendicular) glenCenterlineNorm = glenCenterline/norm(glenCenterline); GHS = SHSvect - dot(SHSvect, glenCenterlineNorm)*glenCenterlineNorm; % GHS amplitude GHSAmpl = norm(GHS) / humHeadRadius; % GHS orientation xProj = dot(GHS, -xAxis); yProj = dot(GHS, yAxis); if xProj ~= 0 if yProj ~= 0 GHSOrient = atan(yProj/xProj); else GHSOrient = pi * (xProj > 0); end else GHSOrient = pi/2 - pi * (yProj > 0); end GHSOrient = rad2deg(GHSOrient); %% Angle between plane of maximum GH subluxation and CT axis [0,0,1] % Not used anymore % GHSPlane = cross(glenCenterline, glenHumHead); % GHSPlaneAngle = vrrotvec(GHSPlane, [0 0 1]); % GHSPlaneAngle = rad2deg(GHSPlaneAngle(4)); % if GHSPlaneAngle > 90 % GHSPlaneAngle = 180 - GHSPlaneAngle; % end %% Set humerus properties obj.center = humHeadCenter; obj.radius = humHeadRadius; obj.SHSAngle = SHSAngle; obj.SHSPA = SHSPA; obj.SHSIS = SHSIS; obj.SHSAmpl = SHSAmpl; obj.SHSOrient = SHSOrient; obj.GHSAmpl = GHSAmpl; obj.GHSOrient = GHSOrient; % Commented ince not used % outputArgStruct = struct('HHResiduals',HHResiduals,'R2HH',R2HH); % Useful? outputArg = 1; end end end diff --git a/ShoulderCase/@Humerus/loadAmiraLandmarks.m b/ShoulderCase/@Humerus/loadAmiraLandmarks.m new file mode 100644 index 0000000..302afb3 --- /dev/null +++ b/ShoulderCase/@Humerus/loadAmiraLandmarks.m @@ -0,0 +1,44 @@ +function outputArg = loadAmiraLandmarks(obj) + %LOAD Load 5 humeral head landmarks + % TODO: Should be check for 'old' vs 'new' file definition + + SCase = obj.shoulder.SCase; + SCaseId4C = SCase.id4C; + amiraDir = SCase.dataAmiraPath; + + % Try 1st with new defiunition of HHLandmarks + fileName = "newHHLandmarks" + SCaseId4C + ".landmarkAscii"; + fileName = fullfile(amiraDir, fileName); + + if exist(fileName,'file') == 2 + importedData = importdata(fileName, ' ', 14); + else % Check for old version of HHLandmarks + fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii']; + fileName = [amiraDir '/' fileName]; + if exist(fileName,'file') == 2 + importedData = importdata(fileName, ' ', 14); + warning([SCaseId4C ' is using old humeral landmarks definition']) + else + outputArg = 0; + return + end + end + % with other methods + + % Check validity of data in file + % Variable points belos should be renamed landmarks, for consistency + try + points = importedData.data; + catch + % Cant get landarks from file + warning([SCaseId4C ' cant load humeral landmarks']) + outputArg = 0; + return; + end + % We should add a test for validity of points + + + obj.landmarks = points; + + outputArg = 1; % Should report on loading result +end \ No newline at end of file diff --git a/ShoulderCase/@Humerus/loadSlicerLandmarks.m b/ShoulderCase/@Humerus/loadSlicerLandmarks.m index 33251ed..9a8e24d 100644 --- a/ShoulderCase/@Humerus/loadSlicerLandmarks.m +++ b/ShoulderCase/@Humerus/loadSlicerLandmarks.m @@ -1,7 +1,12 @@ -function loadSlicerLandmarks(obj) - landmarksFilename = fullfile(... - obj.shoulder.SCase.dataSlicerPath,... - "HH_landmarks.mrk.json"); - fileData = jsondecode(fileread(landmarksFilename)); - obj.landmarks = [fileData.markups.controlPoints.position]'; +function output = loadSlicerLandmarks(obj) + try + landmarksFilename = fullfile(... + obj.shoulder.SCase.dataSlicerPath,... + "HH_landmarks.mrk.json"); + fileData = jsondecode(fileread(landmarksFilename)); + obj.landmarks = [fileData.markups.controlPoints.position]'; + output = 1; + catch + output = 0; + end end \ No newline at end of file diff --git a/ShoulderCase/@RotatorCuff/RotatorCuff.m b/ShoulderCase/@RotatorCuff/RotatorCuff.m index 6163329..2e95829 100644 --- a/ShoulderCase/@RotatorCuff/RotatorCuff.m +++ b/ShoulderCase/@RotatorCuff/RotatorCuff.m @@ -1,105 +1,110 @@ classdef RotatorCuff < handle properties SC SS IS TM anteroPosteriorImbalance = nan; shoulder; end methods function obj = RotatorCuff(shoulder) obj.shoulder = shoulder; obj.SC = Muscle(obj, "SC"); obj.SS = Muscle(obj, "SS"); obj.IS = Muscle(obj, "IS"); obj.TM = Muscle(obj, "TM"); [~, ~] = mkdir(obj.dataPath); [~, ~] = mkdir(obj.dataSlicesPath); end function output = dataPath(obj) output = fullfile(obj.shoulder.dataPath, "muscles"); end function output = dataSlicesPath(obj) output = fullfile(obj.dataPath, "oblique_slices"); end function output = summary(obj) summary = []; summary = [summary; obj.SC.summary]; summary = [summary; obj.SS.summary]; summary = [summary; obj.IS.summary]; summary = [summary; obj.TM.summary]; output = summary; end + function sliceAndSegment(obj) + obj.createRotatorCuffSliceMatthieu(); + obj.segmentMuscles("rotatorCuffMatthieu", "autoMatthieu"); + end + function output = measure(obj, sliceName, segmentationSet) try obj.measureDegenerations(sliceName, segmentationSet); obj.measureInsertions(); obj.measureAnteroPosteriorImbalance(); end output = obj.summary; end function measureDegenerations(obj, sliceName, segmentationSet) for muscleName = ["SC", "SS", "IS", "TM"] obj.(muscleName).measureDegeneration(sliceName, segmentationSet); end end function measureInsertions(obj) humerus = obj.shoulder.humerus; humeralHead = Sphere(humerus.center, humerus.radius); % IS humeral insertion estimation obj.IS.measureInsertionTangentToSphere(... humeralHead,... obj.shoulder.scapula.coordSys.IS,... -obj.shoulder.scapula.coordSys.PA); obj.IS.measureForceVector(); % TM humeral insertion estimation obj.TM.measureInsertionTangentToSphere(... humeralHead,... obj.shoulder.scapula.coordSys.IS,... -obj.shoulder.scapula.coordSys.PA); obj.TM.measureForceVector(); % SC humeral insertion estimation obj.SC.measureInsertionTangentToSphere(... humeralHead,... obj.shoulder.scapula.coordSys.IS,... obj.shoulder.scapula.coordSys.PA); obj.SC.measureForceVector(); end function measureAnteroPosteriorImbalance(obj) forceSC = obj.SC.forceVector; forceIS = obj.IS.forceVector; forceTM = obj.TM.forceVector; scapulaCS = obj.shoulder.scapula.coordSys; forceResultant = Vector(scapulaCS.origin, scapulaCS.origin) + ... forceSC + forceIS + forceTM; % Take infero-superior orthogonal component <=> project on scapular axial plane inferoSuperior = Vector(scapulaCS.origin, scapulaCS.origin + scapulaCS.IS); imbalanceVector = inferoSuperior.orthogonalComplementToPoint(forceResultant.target); lateroMedial = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.ML); anteroPosterior = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.PA); obj.anteroPosteriorImbalance = ... sign(dot(imbalanceVector, anteroPosterior)) *... rad2deg(lateroMedial.angle(imbalanceVector)); end end end \ No newline at end of file diff --git a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m index a9eba51..c02cfdc 100644 --- a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m +++ b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m @@ -1,103 +1,102 @@ classdef ScapulaAuto < Scapula % To be used with ShoulderAuto data. % The load() method requires a specific implementation. % Instanciate a GlenoidAuto object. methods function createGlenoid(obj) obj.glenoid = GlenoidAuto(obj); end function outputArg = load(obj) % Load 2 files obtained by SSM auto segmentation % scapulaSurfaceAuto{L/R}.ply % scapulaLandmarksAuto{L/R}.mat outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise SCase = obj.shoulder.SCase; SCaseId4C = SCase.id4C; matlabDir = SCase.dataMatlabPath; % Set side equal to manual side or 'R' if manual side is not available - manualSide = obj.shoulder.SCase.shoulderManual.side; - if isequal(manualSide,'L') - side = manualSide; - otherSide = 'R'; - elseif isequal(manualSide,'R') - side = manualSide; - otherSide = 'L'; + if not(isempty(obj.shoulder.side)) + side = obj.shoulder.side; + elseif not(isempty(SCase.shoulderManual.side)) + side = SCase.shoulderManual.side; else - side = 'R'; - otherSide = 'L'; + side = "R"; end + assert(ismember(side, ["R" "L"]), "Shoulder side must be R or L"); + otherSide = char((side == "R")*'L' + (side == "L")*'R'); + % Try loading auto segmentation from matlab dir - fileName = ['scapulaSurfaceAuto' side '.ply']; + fileName = "scapulaSurfaceAuto" + side + ".ply"; fileName = fullfile(matlabDir, fileName); if ~exist(fileName,'file') side = otherSide; - fileName = ['scapulaSurfaceAuto' side '.ply']; + fileName = "scapulaSurfaceAuto" + side + ".ply"; fileName = fullfile(matlabDir, fileName); end %side = SCase.shoulder.side; if exist(fileName,'file') == 2 try face_index_start = 1; [points,faces] = loadPly(fileName,face_index_start); % % Load ply file % ptCloud = pcread(fileName); % load the pointCloud object in the ply file % points = ptCloud.Location; % get the array of (x,y,z) points obj.surface.points = points; obj.surface.faces = faces; obj.segmentation = 'A'; % Auto segmentation from matlab catch obj.segmentation = 'E'; % Error on loading outputArg = 0; end else obj.segmentation = 'N'; % No segmentation file outputArg = 0; end % Try loading auto scapula landmarks from matlab dir - filename = ['scapulaLandmarksAuto' side '.mat']; - filename = [matlabDir '/' filename]; + filename = "scapulaLandmarksAuto" + side + ".mat"; + filename = matlabDir + "/" + filename; if exist(filename,'file') == 2 try % Load mat file with scapula landmarks load(filename, 'ScapulaLandmarks'); % Set loaded landmarks to scapulaAuto (obj) obj.angulusInferior = ScapulaLandmarks.angulusInferior; obj.trigonumSpinae = ScapulaLandmarks.trigonumSpinae; obj.processusCoracoideus = ScapulaLandmarks.processusCoracoideus; obj.acromioClavicular = ScapulaLandmarks.acromioClavicular; obj.angulusAcromialis = ScapulaLandmarks.angulusAcromialis; obj.spinoGlenoidNotch = ScapulaLandmarks.spinoGlenoidNotch; obj.pillar = ScapulaLandmarks.pillar; obj.groove = ScapulaLandmarks.groove; catch disp(SCaseId4C); % For debug (2 cases) outputArg = 0; end else outputArg = 0; end try obj.measurePlane; catch outputArg = 0; end end end end diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m index c1a7025..31ae096 100644 --- a/ShoulderCase/@Shoulder/Shoulder.m +++ b/ShoulderCase/@Shoulder/Shoulder.m @@ -1,24 +1,46 @@ classdef (Abstract) Shoulder < handle % Contains all the shoulder parts (bones) and their measurements. % % Can be used to plot an overview of the case. properties side = ""; % R or L scapula humerus rotatorCuff CTscan = ""; comment = ""; SCase end methods function obj = Shoulder(SCase) obj.SCase = SCase; obj.humerus = Humerus(obj); obj.rotatorCuff = RotatorCuff(obj); end + + function loadData(obj) + obj.scapula.load(); + obj.scapula.glenoid.load(); + obj.humerus.load(); + end + + function measureMorphology(obj) + obj.scapula.coordSysSet(); + obj.scapula.glenoid.morphology(); + obj.humerus.subluxation(); + obj.scapula.acromion.morphology(); + end + + function measureDensity(obj) + obj.scapula.genoid.calcDensity(); + end + + function measureMuscles(obj) + obj.scapula.rotatorCuff.sliceAndSegment(); + obj.scapula.rotatorCuff.measure("rotatorCuffMatthieu", "autoMatthieu"); + end end end