diff --git a/ShoulderCase/Scapula.m b/ShoulderCase/Scapula.m index b7dc42b..82f4dd9 100644 --- a/ShoulderCase/Scapula.m +++ b/ShoulderCase/Scapula.m @@ -1,338 +1,338 @@ classdef 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 surface % surface points and triangles of the scapula segmentation % either none 'N', manual 'M' or automatic 'A' glenoid % Glenoid class acromion % Acromion class comment % any comment end properties (Hidden = true) SCase end methods (Access = ?Shoulder) % Only Shoulder is allowed to construct a Scapula function obj = Scapula(SCase) %SCAPULA Construct an instance of this class % Constructor of the scapula object, sets all properties to % zero. obj.SCase = SCase; % All definition below should be replaced by [] obj.angulusInferior = []; obj.angulusInferiorLocal = []; obj.trigonumSpinae = []; obj.processusCoracoideus = []; obj.acromioClavicular = []; obj.angulusAcromialis = []; obj.spinoGlenoidNotch = []; obj.pillar = []; obj.groove = []; obj.surface = []; obj.segmentation = 'N'; obj.coordSys.origin = []; obj.coordSys.PA = []; % X postero-anterior axis obj.coordSys.IS = []; % Y infero-superior axis obj.coordSys.ML = []; % Z medio-lateral axis obj.glenoid = Glenoid(SCase); obj.acromion = Acromion(SCase); obj.comment = ''; end end methods function outputArg = load(obj) % LOAD Load segmented surface and landmnarks % Load the segmented scapula surface (if exist) and the % scapula landmarks (6 scapula, 5 groove, 5 pillar) from % amira directory. SCase = obj.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; matlabDir = SCase.dataMatlabPath; outputArg = 1; % Should report on loading result % Scapula (6) landmarks fileName = ['ScapulaLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); else - disp('MATLAB:rmpath:DirNotFound',... - 'Could not find file: %s',fileName); + disp(['MATLAB:rmpath:DirNotFound',... + 'Could not find file: ',fileName]); outputArg = 0; return; end % Check that imported daat are well formed try landmarks = importedData.data; catch outputArg = 0; return; end obj.angulusInferior = landmarks(1,:); obj.trigonumSpinae = landmarks(2,:); obj.processusCoracoideus = landmarks(3,:); obj.acromioClavicular = landmarks(4,:); obj.angulusAcromialis = landmarks(5,:); obj.spinoGlenoidNotch = landmarks(6,:); % Scapula suprapinatus groove (5) landmarks fileName = ['ScapulaGrooveLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); landmarks = importedData.data; else - disp('MATLAB:rmpath:DirNotFound',... - 'Could not find file: %s',fileName); + disp(['MATLAB:rmpath:DirNotFound',... + 'Could not find file: ',fileName]); outputArg = 0; end % Check groove landmarks validity if length(landmarks) > 5 landmarks = landmarks(1:5,:); warning(['More than 5 groove landmarks(' SCase.id ')']); elseif length(landmarks) < 5 error('Less than 5 groove landmarks'); end obj.groove = landmarks; % Scapula pillar (5) landmarks fileName = ['ScapulaPillarLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); landmarks = importedData.data; else - disp('MATLAB:rmpath:DirNotFound',... - 'Could not find file: %s',fileName); + disp(['MATLAB:rmpath:DirNotFound',... + 'Could not find file: ',fileName]); outputArg = 0; end obj.pillar = landmarks; % Import scapula surface % Try manual segmentation in amira folder fileName = ['scapula_' SCaseId4C '.stl']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 try [points,faces,~]=loadStl(fileName,1); obj.surface.points = points; obj.surface.faces = faces; obj.segmentation = 'M'; % Manual segmentation from amira catch obj.segmentation = 'E'; % Error in loading end else obj.segmentation = 'N'; % No segmentation end end function outputArg = loadAuto(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.SCase; SCaseId4C = SCase.id4C; matlabDir = SCase.dataMatlabPath; side = SCase.shoulder.side; % We might no have it yet !! % Try loading auto segmentation from matlab dir fileName = ['scapulaSurfaceAuto' side '.ply']; fileName = [matlabDir '/' fileName]; 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]; 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 end function outputArg = coordSysSet(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. % Scapular plane is fitted on 3 points (angulusInferior, % trigonumSpinae, most laretal scapular goove landmark). % When the scapular coordinate is set, we then express the % local scapula landmarks in this system. SCase = obj.SCase; % Define local variables for landmarks aiLand = obj.angulusInferior; tsLand = obj.trigonumSpinae; pcLand = obj.processusCoracoideus; % acLand = obj.acromioClavicular; % not used here aaLand = obj.angulusAcromialis; snLand = obj.spinoGlenoidNotch; % spLand = obj.pillar; % not used here sgLand = obj.groove; % Test validity of groove if length(sgLand) ~= 5 outputArg = 0; return end % Define approximate anatomical axes (not orthogonal) postAnt = pcLand - aaLand; % ~ X infSup = snLand - aiLand; % ~ Y medLat = snLand - tsLand; % ~ Z % Find scapula side (right/left) and set to shoulder.side infSupRight = cross(medLat, postAnt); % inf-sup for right if dot(infSup, infSupRight) > 0 SCase.shoulder.side = 'R'; % Right scapula else SCase.shoulder.side = 'L'; % Left scapula end % Fit scapular plane to aiLand, tsLand and most distal scapula % groove landmarks % Determine position of most lateral landmark of the scapula groove % by calculating the distance between ts landmark and each of the 5 % landmarks of the scapula groove. dist = zeros(length(sgLand),1); for i=1:length(sgLand); dist(i) = sqrt(... (tsLand(1)-sgLand(i,1))^2 + ... (tsLand(2)-sgLand(i,2))^2 + ... (tsLand(3)-sgLand(i,3))^2 ... ); end maxDist = max(dist); % Maximum distance pos = dist == maxDist; % Location of point with max distance in the array sgLatLan = sgLand(pos,:); % Coordnates of most lateral point of the scapula groove sgLatLan = sgLatLan(1,:); % If there are more than 1, take the 1st % Fit scapular plane on aiLand, tsLand and sgLat [scapPlane, PlaneMean, ~, PlaneRMSE, ~] = ... fitPlane([aiLand; tsLand; sgLatLan]); scapPlane = scapPlane'; % Unit normal row vector % Set scapulaPlaneNormal oriented posterior-anterior (as x axis) if dot(scapPlane, postAnt) < 0 scapPlane = -scapPlane; end PAaxis = scapPlane; % PA (X) axis correspond to scapular plane normal % Scapular transverse (medio-lateral) axis % 1) Project groove points to scapular plane % 2) Fit axis to projected points % 3) Set medio-lateral direction sgProj = project2Plane(sgLand,scapPlane,sgLatLan,size(sgLand,1)); [grooveAxis, grooveMean, ~, grooveRMSE, R2Line] = fitLine(sgProj); transAxis = (grooveAxis/norm(grooveAxis))'; % Set transvere axis direction to med-lat if dot(transAxis, medLat) < 0 transAxis = -transAxis; end MLaxis = transAxis; % ML (Z) axis along scapular transverse axis % Infero-superio axis ISaxis = cross(PAaxis, MLaxis); if dot(ISaxis, infSup) < 0 % Check orientation of axis ISaxis = -ISaxis; end % !! Code below is not fully correct: grooveMean is not in % scapPlane nor in scapAxis % Set the origin of the coordinate system to the spino-gemoid notch % projected on the scapular axis transAxisMean = project2Plane(grooveMean,scapPlane,PlaneMean,1); % Project spinoGlenoidNotch to transverseAxis to set Origin origin = transAxisMean + ... dot((snLand - transAxisMean),transAxis)*transAxis; obj.coordSys.origin = origin; obj.coordSys.PA = PAaxis; obj.coordSys.IS = ISaxis; obj.coordSys.ML = MLaxis; % Express the scapula landmarks in the scapular coordinate % system R = [PAaxis' ISaxis', MLaxis']; % Rotation matrix to align to scapula coordinate system obj.angulusInferiorLocal = (obj.angulusInferior - obj.coordSys.origin) * R; % Not used % outputStuct = struct('PlaneRMSE',PlaneRMSE, 'GrooveRMSE', grooveRMSE, 'R2Line', R2Line); outputArg = 1; end end end