Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/ShoulderCase/Scapula.m
502abc4d511dBashardoust_SystCoord

Last Change
Raw File

Scapula.m
View Options

classdef Scapula < handle
%SCAPULA Summary of this class goes here
% This class defines the scapula. Landmarks are used to define its
% coordinate system. It includes the glenoid object.
properties
angulusInferior % Landmark Angulus Inferior read from amira
angulusInferiorLocal % Same as above, expressed the scapular coordinate system rather than the CT coordinate system
trigonumSpinae % Landmark Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine read from amira
processusCoracoideus % Landmark Processus Coracoideus (most lateral point) read from amira
acromioClavicular % Landmark Acromio-clavicular joint (most lateral part on acromion)read from amira
angulusAcromialis % Landmark Angulus Acromialis (most laterodorsal point) read from amira
spinoGlenoidNotch % Landmark Spino-glenoid notch read from amira
pillar % 5 landmarks on the pillar
groove % 5 landmarks on the scapula (supraspinatus) groove
coordSys % Scapular coordinate system
surface % surface points and triangles of the scapula
segmentation % either none 'N', manual 'M' or automatic 'A'
glenoid % Glenoid class
acromion % Acromion class
comment % any comment
end
properties (Hidden = true)
SCase
end
methods (Access = ?Shoulder) % Only Shoulder is allowed to construct a Scapula
function obj = Scapula(SCase)
%SCAPULA Construct an instance of this class
% Constructor of the scapula object, sets all properties to
% zero.
obj.SCase = SCase;
% All definition below should be replaced by []
obj.angulusInferior = [];
obj.angulusInferiorLocal = [];
obj.trigonumSpinae = [];
obj.processusCoracoideus = [];
obj.acromioClavicular = [];
obj.angulusAcromialis = [];
obj.spinoGlenoidNotch = [];
obj.pillar = [];
obj.groove = [];
obj.surface = [];
obj.segmentation = 'N';
obj.coordSys.origin = [];
obj.coordSys.PA = []; % X postero-anterior axis
obj.coordSys.IS = []; % Y infero-superior axis
obj.coordSys.ML = []; % Z medio-lateral axis
obj.glenoid = Glenoid(SCase);
obj.acromion = Acromion(SCase);
obj.comment = '';
end
end
methods
function outputArg = load(obj)
% LOAD Load segmented surface and landmnarks
% Load the segmented scapula surface (if exist) and the
% scapula landmarks (6 scapula, 5 groove, 5 pillar) from
% amira directory.
SCase = obj.SCase;
SCaseId4C = SCase.id4C;
amiraDir = SCase.dataAmiraPath;
matlabDir = SCase.dataMatlabPath;
outputArg = 1; % Should report on loading result
% Scapula (6) landmarks
fileName = ['ScapulaLandmarks' SCaseId4C '.landmarkAscii'];
fileName = [amiraDir '/' fileName];
if exist(fileName,'file') == 2
importedData = importdata(fileName, ' ', 14);
else
disp(['MATLAB:rmpath:DirNotFound',...
'Could not find file: ',fileName]);
outputArg = 0;
return;
end
% Check that imported data are well formed
try
landmarks = importedData.data;
catch
outputArg = 0;
return;
end
obj.angulusInferior = landmarks(1,:);
obj.trigonumSpinae = landmarks(2,:);
obj.processusCoracoideus = landmarks(3,:);
obj.acromioClavicular = landmarks(4,:);
obj.angulusAcromialis = landmarks(5,:);
obj.spinoGlenoidNotch = landmarks(6,:);
% Scapula suprapinatus groove (5) landmarks
fileName = ['ScapulaGrooveLandmarks' SCaseId4C '.landmarkAscii'];
fileName = [amiraDir '/' fileName];
if exist(fileName,'file') == 2
importedData = importdata(fileName, ' ', 14);
% Check that imported data are well formed
try
landmarks = importedData.data;
catch
outputArg = 0;
return;
end
else
disp(['MATLAB:rmpath:DirNotFound',...
'Could not find file: ',fileName]);
outputArg = 0;
end
% Check groove landmarks validity
if length(landmarks) > 5
landmarks = landmarks(1:5,:);
warning(['More than 5 groove landmarks(' SCase.id ')']);
elseif length(landmarks) < 5
error('Less than 5 groove landmarks');
end
obj.groove = landmarks;
% Scapula pillar (5) landmarks
fileName = ['ScapulaPillarLandmarks' SCaseId4C '.landmarkAscii'];
fileName = [amiraDir '/' fileName];
if exist(fileName,'file') == 2
importedData = importdata(fileName, ' ', 14);
% Check that imported data are well formed
try
landmarks = importedData.data;
catch
outputArg = 0;
return;
end
else
disp(['MATLAB:rmpath:DirNotFound',...
'Could not find file: ',fileName]);
outputArg = 0;
end
obj.pillar = landmarks;
% Import scapula surface
% Try manual segmentation in amira folder
fileName = ['scapula_' SCaseId4C '.stl'];
fileName = [amiraDir '/' fileName];
if exist(fileName,'file') == 2
try
[points,faces,~]=loadStl(fileName,1);
obj.surface.points = points;
obj.surface.faces = faces;
obj.segmentation = 'M'; % Manual segmentation from amira
catch
obj.segmentation = 'E'; % Error in loading
end
else
obj.segmentation = 'N'; % No segmentation
end
end
function outputArg = loadAuto(obj)
% Load 2 files obtained by SSM auto segmentation
% scapulaSurfaceAuto{L/R}.ply
% scapulaLandmarksAuto{L/R}.mat
outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise
SCase = obj.SCase;
SCaseId4C = SCase.id4C;
matlabDir = SCase.dataMatlabPath;
side = SCase.shoulder.side; % We might no have it yet !!
% Try loading auto segmentation from matlab dir
fileName = ['scapulaSurfaceAuto' side '.ply'];
fileName = [matlabDir '/' fileName];
if exist(fileName,'file') == 2
try
face_index_start = 1;
[points,faces] = loadPly(fileName,face_index_start); % % Load ply file
% ptCloud = pcread(fileName); % load the pointCloud object in the ply file
% points = ptCloud.Location; % get the array of (x,y,z) points
obj.surface.points = points;
obj.surface.faces = faces;
obj.segmentation = 'A'; % Auto segmentation from matlab
catch
obj.segmentation = 'E'; % Error on loading
outputArg = 0;
end
else
obj.segmentation = 'N'; % No segmentation file
outputArg = 0;
end
% Try loading auto scapula landmarks from matlab dir
filename = ['scapulaLandmarksAuto' side '.mat'];
filename = [matlabDir '/' filename];
if exist(filename,'file') == 2
try
% Load mat file with scapula landmarks
load(filename, 'ScapulaLandmarks');
% Set loaded landmarks to scapulaAuto (obj)
obj.angulusInferior = ScapulaLandmarks.angulusInferior;
obj.trigonumSpinae = ScapulaLandmarks.trigonumSpinae;
obj.processusCoracoideus = ScapulaLandmarks.processusCoracoideus;
obj.acromioClavicular = ScapulaLandmarks.acromioClavicular;
obj.angulusAcromialis = ScapulaLandmarks.angulusAcromialis;
obj.spinoGlenoidNotch = ScapulaLandmarks.spinoGlenoidNotch;
obj.pillar = ScapulaLandmarks.pillar;
obj.groove = ScapulaLandmarks.groove;
catch
disp(SCaseId4C); % For debug (2 cases)
outputArg = 0;
end
else
outputArg = 0;
end
end
function outputArg = coordSysSet(obj)
% 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.
% Scapular plane is fitted on 3 points (angulusInferior,
% trigonumSpinae, most laretal scapular goove landmark).
% When the scapular coordinate is set, we then express the
% local scapula landmarks in this system.
SCase = obj.SCase;
% Define local variables for landmarks
aiLand = obj.angulusInferior;
tsLand = obj.trigonumSpinae;
pcLand = obj.processusCoracoideus;
% acLand = obj.acromioClavicular; % not used here
aaLand = obj.angulusAcromialis;
snLand = obj.spinoGlenoidNotch;
% spLand = obj.pillar; % not used here
sgLand = obj.groove;
% Test validity of groove
if length(sgLand) ~= 5
outputArg = 0;
return
end
% Define approximate anatomical axes (not orthogonal)
postAnt = pcLand - aaLand; % ~ X
infSup = snLand - aiLand; % ~ Y
medLat = snLand - tsLand; % ~ Z
% Find scapula side (right/left) and set to shoulder.side
infSupRight = cross(medLat, postAnt); % inf-sup for right
if dot(infSup, infSupRight) > 0
SCase.shoulder.side = 'R'; % Right scapula
else
SCase.shoulder.side = 'L'; % Left scapula
end
% Fit scapular plane to aiLand, tsLand and most distal scapula
% groove landmarks
% Determine position of most lateral landmark of the scapula groove
% by calculating the distance between ts landmark and each of the 5
% landmarks of the scapula groove.
dist = zeros(length(sgLand),1);
for i=1:length(sgLand);
dist(i) = sqrt(...
(tsLand(1)-sgLand(i,1))^2 + ...
(tsLand(2)-sgLand(i,2))^2 + ...
(tsLand(3)-sgLand(i,3))^2 ...
);
end
maxDist = max(dist); % Maximum distance
pos = dist == maxDist; % Location of point with max distance in the array
sgLatLan = sgLand(pos,:); % Coordnates of most lateral point of the scapula groove
sgLatLan = sgLatLan(1,:); % If there are more than 1, take the 1st
% Fit scapular plane on aiLand, tsLand and sgLat
[scapPlane, PlaneMean, ~, PlaneRMSE, ~] = ...
fitPlane([aiLand; tsLand; sgLatLan]);
scapPlane = scapPlane'; % Unit normal row vector
% Set scapulaPlaneNormal oriented posterior-anterior (as x axis)
if dot(scapPlane, postAnt) < 0
scapPlane = -scapPlane;
end
PAaxis = scapPlane; % PA (X) axis correspond to scapular plane normal
% Scapular transverse (medio-lateral) axis
% 1) Project groove points to scapular plane
% 2) Fit axis to projected points
% 3) Set medio-lateral direction
sgProj = project2Plane(sgLand,scapPlane,sgLatLan,size(sgLand,1));
[grooveAxis, grooveMean, ~, grooveRMSE, R2Line] = fitLine(sgProj);
transAxis = (grooveAxis/norm(grooveAxis))';
% Set transvere axis direction to med-lat
if dot(transAxis, medLat) < 0
transAxis = -transAxis;
end
MLaxis = transAxis; % ML (Z) axis along scapular transverse axis
% Infero-superio axis
ISaxis = cross(PAaxis, MLaxis);
if dot(ISaxis, infSup) < 0 % Check orientation of axis
ISaxis = -ISaxis;
end
% !! Code below is not fully correct: grooveMean is not in
% scapPlane nor in scapAxis
% Set the origin of the coordinate system to the spino-gemoid notch
% projected on the scapular axis
transAxisMean = project2Plane(grooveMean,scapPlane,PlaneMean,1);
% Project spinoGlenoidNotch to transverseAxis to set Origin
origin = transAxisMean + ...
dot((snLand - transAxisMean),transAxis)*transAxis;
obj.coordSys.origin = origin;
obj.coordSys.PA = PAaxis;
obj.coordSys.IS = ISaxis;
obj.coordSys.ML = MLaxis;
% Express the scapula landmarks in the scapular coordinate
% system
R = [PAaxis' ISaxis', MLaxis']; % Rotation matrix to align to scapula coordinate system
obj.angulusInferiorLocal = (obj.angulusInferior - obj.coordSys.origin) * R;
% Not used
% outputStuct = struct('PlaneRMSE',PlaneRMSE, 'GrooveRMSE', grooveRMSE, 'R2Line', R2Line);
outputArg = 1;
end
end
end
c4science · Help