Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (master) Browse

R8218/ShoulderCase/CoordinateSystem.m
a04f3ed203e2master

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 output = getRotationMatrix(obj)
output = obj.rotationMatrix;
end
function localPoints = express(obj,globalPoints)
% Give the coordinates of global cartesian points in present coordinate system
assert(size(globalPoints,2) == 3,'Points have to be a Nx3 double array.')
localPoints = (globalPoints-obj.origin) * obj.rotationMatrix;
end
function globalPoints = getGlobalCoordinatesOfLocalPoints(obj,localPoints)
% Give the coordinates of local points in global cartesian coordinate system
assert(size(localPoints,2) == 3,'Points have to be a Nx3 double array.')
globalPoints = (localPoints * obj.rotationMatrix') + obj.origin;
end
function projected = projectOnXYPlane(obj,points)
projectionPlane = Plane;
projectionPlane.fit([obj.origin;(obj.origin + obj.xAxis);(obj.origin + obj.yAxis)]);
projected = projectionPlane.projectOnPlane(points);
end
function projected = projectOnYZPlane(obj,points)
projectionPlane = Plane;
projectionPlane.fit([obj.origin;(obj.origin + obj.yAxis);(obj.origin + obj.zAxis)]);
projected = projectionPlane.projectOnPlane(points);
end
function projected = projectOnZXPlane(obj,points)
projectionPlane = Plane;
projectionPlane.fit([obj.origin;(obj.origin + obj.zAxis);(obj.origin + obj.xAxis)]);
projected = projectionPlane.projectOnPlane(points);
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) < 2)
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