Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/ShoulderCase/CoordinateSystem.m
ceac7cb096efaterrier-plotUpdate

Last Change
Raw File

CoordinateSystem.m
View Options

classdef CoordinateSystem < handle
properties
origin
xAxis
yAxis
zAxis
end
properties (Hidden = true)
rotationMatrix
% enable retro-compatibility with codes that uses PA, IS, and ML axes
PA
IS
ML
end
methods
function obj = CoordinateSystem()
obj.origin = [0 0 0];
obj.xAxis = [1 0 0];
obj.yAxis = [0 1 0];
obj.zAxis = [0 0 1];
end
function set.xAxis(obj,value)
% enable retro-compatibility with codes that uses a PA axis
obj.xAxis = value;
obj.PA = value;
obj.rotationMatrix(:,1) = value';
end
function set.yAxis(obj,value)
% enable retro-compatibility with codes that uses an IS axis
obj.yAxis = value;
obj.IS = value;
obj.rotationMatrix(:,2) = value';
end
function set.zAxis(obj,value)
% enable retro-compatibility with codes that uses a ML axis
obj.zAxis = value;
obj.ML = value;
obj.rotationMatrix(:,3) = value';
end
function points = express(obj,points)
% Give the coordinates of points in present coordinate system
if (size(points,2) ~= 3)
error('Points have to be a Nx3 double array.')
return;
end
points = (points-obj.origin) * obj.rotationMatrix;
end
function output = isOrthogonal(obj)
% The axes comparison can't be done looking for strict equality due to
% rounded values. Therefore the values must be evaluated with a given
% tolerance
minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of
% a millionth of a radian
if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle)
output = false;
return
end
if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle)
output = false;
return
end
if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle)
output = false;
return
end
output = true;
end
function output = isRightHanded(obj)
% The comparison between a theoretical right-handed axis and the actual
% z Axis can't be done looking for strict vector equality due to rounded
% values. Therefore the norm of the sum of the two vectors is evaluated.
if not(obj.isOrthogonal)
output = false;
return
end
rightHandedAxis = cross(obj.xAxis,obj.yAxis);
if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis)
output = false;
return
end
output = true;
end
function output = isLeftHanded(obj)
if not(obj.isOrthogonal)
output = false;
return
end
if obj.isRightHanded
output = false;
return
end
output = true;
end
function setSystemWithScapulaLandmarks(obj,scapula)
% 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.
if (length(scapula.groove) ~= 5)
error(['Can''t set scapular coordinate system with actual',...
'scapula.groove']);
return;
end
% The four following methods have to be called in this specific order
obj.setXAxisWithScapulaLandmarks(scapula);
obj.setZAxisWithScapulaLandmarks(scapula);
obj.setYAxisWithScapulaLandmarks(scapula);
obj.setOriginWithScapulaLandmarks(scapula);
end
function setSystemWith(obj)
%
%
%
%
%
%
%
%
%
%
%
%
%
% New method to set a coordinate system can be written here
%
%
%
%
%
%
%
%
%
%
%
%
end
end
methods (Access = private, Hidden = true)
function setOriginWithScapulaLandmarks(obj,scapula)
% Set the origin of the coordinate system to the spino-gemoid notch
% projected on the scapular axis
spinoGlenoid = scapula.spinoGlenoidNotch;
grooveMeanProjection = mean(scapula.plane.projectOnPlane(scapula.groove));
obj.origin = grooveMeanProjection + ...
dot((spinoGlenoid - grooveMeanProjection),obj.zAxis)*obj.zAxis;
end
function setXAxisWithScapulaLandmarks(obj,scapula)
% Set the postero-anterior axis to be the scapula plane normal
obj.xAxis = scapula.plane.normal;
end
function setYAxisWithScapulaLandmarks(obj,scapula)
% Set the infero-superior axis to be orthogonal with the two other axes
superior = scapula.spinoGlenoidNotch;
inferior = scapula.angulusInferior;
obj.yAxis = cross(obj.xAxis,obj.zAxis);
obj.yAxis = orientVectorToward(obj.yAxis,(superior-inferior));
end
function setZAxisWithScapulaLandmarks(obj,scapula)
% Set the media-lateral axis to be the line fitted to the projection of
% the groove points on the scapula plane
lateral = scapula.spinoGlenoidNotch;
medial = scapula.trigonumSpinae;
groovePointsProjection = scapula.plane.projectOnPlane(scapula.groove);
grooveAxis = fitLine(groovePointsProjection);
grooveAxis = (grooveAxis/norm(grooveAxis))';
obj.zAxis = orientVectorToward(grooveAxis,(lateral-medial));
end
end
end
c4science · Help