diff --git a/ShoulderCase/Humerus.m b/ShoulderCase/Humerus.m old mode 100644 new mode 100755 index 7525770..89623b0 --- a/ShoulderCase/Humerus.m +++ b/ShoulderCase/Humerus.m @@ -1,245 +1,258 @@ 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: 2018-08-10 + % Revision date: 2019-06-29 % - % TO DO: + % 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) end - + properties (Hidden = true) SCase end - + methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Humerus function obj = Humerus(SCase) %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.SCase = SCase; end end - + methods function outputArg = load(obj) %LOAD Load 5 humeral head landmarks - % TODO: Should be check for 'old' vs 'new' file definition - + % TODO: Should be check for 'old' vs 'new' file definition + SCase = obj.SCase; SCaseId4C = SCase.id4C; amiraDir = SCase.dataAmiraPath; - + + % Try 1st with new defiunition of HHLandmarks fileName = ['newHHLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; - + if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); - points = importedData.data; else % Check for old version of HHLandmarks fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii']; fileName = [amiraDir '/' fileName]; if exist(fileName,'file') == 2 importedData = importdata(fileName, ' ', 14); - points = importedData.data; - warning([SCaseId4C ' is using old humeral landmarks definition']) + warning([SCaseId4C ' is using old humeral landmarks definition']) else error('MATLAB:rmpath:DirNotFound',... 'Could not find file: %s',fileName) 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 - + 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 + % Radius of curvature of humeral head articular surface % Might be removed anatomy for concistency - + SCase = obj.SCase; - + %% Get data from parent shoulder object % Shoulder side side = SCase.shoulder.side; - + %% Get data from parent scapula object % Scapula coordinate system xAxis = SCase.shoulder.scapula.coordSys.PA; yAxis = SCase.shoulder.scapula.coordSys.IS; zAxis = SCase.shoulder.scapula.coordSys.ML; - + %% Get data from parent glenoid object % Glenoid center glenCenter = SCase.shoulder.scapula.glenoid.center; % Glenoid centerLine from parent glenoid object glenCenterline = SCase.shoulder.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 + % 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 + 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? + % outputArgStruct = struct('HHResiduals',HHResiduals,'R2HH',R2HH); % Useful? outputArg = 1; end end end -