Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/ShoulderCase/Humerus.m
492db744a284Bashardoust_SystCoord

Last Change
Raw File

Humerus.m
View Options

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)
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
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);
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.SCase;
%% Get data from parent shoulder object
% Shoulder side
side = SCase.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
c4science · Help