diff --git a/ShoulderCase/@Humerus/Humerus.m b/ShoulderCase/@Humerus/Humerus.m index 1eccad8..a4e4c43 100644 --- a/ShoulderCase/@Humerus/Humerus.m +++ b/ShoulderCase/@Humerus/Humerus.m @@ -1,258 +1,259 @@ 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 end properties (Hidden = true) 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 = [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 function outputArg = subluxation(obj,scapula) % 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; %% 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/measureSubluxationIndex3D.m b/ShoulderCase/@Humerus/measureSubluxationIndex3D.m new file mode 100644 index 0000000..822979a --- /dev/null +++ b/ShoulderCase/@Humerus/measureSubluxationIndex3D.m @@ -0,0 +1,15 @@ +function measureSubluxationindex3D(obj) +% - GL: (Glenoid Line) Line that joins the most anterior glenoid rim point and +% the most posterior glenoid rim point. +% - HHv: (Humeral Head vector) Vector from the glenoid center to the humeral +% head’s fitted sphere center. +% - HHp: (Humeral Head projection) Projection of HHv on GL. +% - HHr: (Humeral Head radius) Humeral head’s fitted sphere radius. +% +% The 3D humeral subluxation index is the ratio of ||HHp|| + ||HHr|| over 2*||HHr||. + GL = obj.shoulder.scapula.glenoid.posteroAnteriorLine; + HHv = Vector(obj.shoulder.scapula.glenoid.center, obj.center); + HHp = GL.project(HHv); + HHr = obj.radius; + obj.subluxationIndex3D = (HHp.norm + HHr) / (2*HHr); +end \ No newline at end of file