Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/ShoulderCase/@Glenoid/calcDensity.m
5e81a03b9817Bashardoust_SystCoord

Last Change
Raw File

calcDensity.m
View Options

classdef Glenoid < handle
%GLENOID Summary of this class goes here
% Detailed explanation goes here
% Need to load surface points from Amira (or auto-segment)
% Need to set anatomical calculation, and use the scapular coord.Syst
% TODO
% loadMatlab()
properties
surface % points of the glenoid surface load glenoid surface
center % Glenoid center in CT coordinate system
centerLocal % Glenoid center in scapula coordinate system
radius
sphereRMSE
centerLine
depth
width
height
versionAmpl
versionOrient
version
inclination
density
walch
comment
end
properties (Hidden = true)
SCase
end
methods (Access = ?Scapula) % Only Scapula is allowed to construct a Glenoid
function obj = Glenoid(SCase)
%GLENOID Construct an instance of this class
% Detailed explanation goes here
obj.surface = [];
obj.center = [];
obj.centerLocal.x = [];
obj.centerLocal.y = [];
obj.centerLocal.z = [];
obj.radius = [];
obj.sphereRMSE = [];
obj.depth = [];
obj.width = [];
obj.height = [];
obj.versionAmpl = [];
obj.versionOrient = [];
obj.version = [];
obj.inclination = [];
obj.density = [];
obj.walch = [];
obj.comment = [];
obj.SCase = SCase;
end
end
methods
function outputArg = load(obj)
% LOAD Summary of this method goes here
% Load the points of the glenoid surface from the amira
% directoty.
SCase = obj.SCase;
SCaseId4C = SCase.id4C;
amiraDir = SCase.dataAmiraPath;
% Import glenoid surface points
fileName = ['ExtractedSurface' SCaseId4C '.stl'];
fileName = [amiraDir '/' fileName];
if exist(fileName,'file') == 2
[points,faces,~]=loadStl(fileName,1);
obj.surface.points = points;
obj.surface.faces = faces;
else
error(['Could not find file: ' fileName])
end
outputArg = 1; % Should report on loading result
end
function outputArg = loadAuto(obj)
% LOADAUTO Load genoid surface from auto segmentation
% Load the points of the glenoid surface from matlab dir
SCase = obj.SCase;
SCaseId4C = SCase.id4C;
matlabDir = SCase.dataMatlabPath;
side = SCase.shoulder.side;
% Import glenoid surface points
fileName = ['glenoidSurfaceAuto' side '.ply'];
fileName = [matlabDir '/' fileName];
if exist(fileName,'file') == 2
[points,faces]=loadPly(fileName,1);
obj.surface.points = points;
obj.surface.faces = faces;
else
outputArg = 0;
return
% error(['Could not find file: ' fileName])
end
outputArg = 1; % Should report on loading result
end
function outputArg = morphology(obj)
% Calculate morphological properties of the glenoid from its surface:
% center, radius, depth, width, height,
% versionAmpl, versionOrient, version, inclination
SCase = obj.SCase;
glenSurf = obj.surface.points;
% Get local coordinate system from parent scapula object
origin = SCase.shoulder.scapula.coordSys.origin;
xAxis = SCase.shoulder.scapula.coordSys.PA;
yAxis = SCase.shoulder.scapula.coordSys.IS;
zAxis = SCase.shoulder.scapula.coordSys.ML;
% Glenoid center is the closet point between the mean of the extracted
% surface and all the points of the extracted surface.
% This might be improved by find the "real" center, instead of
% the closest one from the surface nodes.
glenCenter = mean(glenSurf);
vect = zeros(size(glenSurf));
dist = zeros(size(glenSurf,1),1);
for i=1:size(glenSurf,1)
vect(i,:) = glenSurf(i,:) - glenCenter;
dist(i) = norm(vect(i,:));
end
glenCenterNode = find(dist == min(dist));
glenCenter = glenSurf(glenCenterNode,:);
% Fit sphere on glenoid surface
[sphereCenter,sphereRadius,sphereResiduals, R2Glenoid] = fitSphere(glenSurf);
sphereCenter = sphereCenter'; % transform a vertical to horizotal vector
glenRadius = sphereRadius;
% Root Mean Square Error
sphereRMSE = norm(sphereResiduals) / sqrt(length(sphereResiduals));
% Glenoid centerline
centerLine = sphereCenter - glenCenter;
% Depth of the glenoid surface
% Project the points of the surface on the centerline
glenSurfProj = zeros(size(glenSurf));
for i = 1:size(glenSurf,1)
glenSurfProj(i,:) = glenCenter + (dot(centerLine,(glenSurf(i,:) - glenCenter)) / ...
dot(centerLine,centerLine)) * centerLine;
dist(i) = norm(glenCenter - glenSurfProj(i,:));
end
depth = max(dist);
% Glenoid width and height
% This part should be rewritten
width = '';
height = '';
% Glenoid version 3D, defined by the amplitude and orientation of the angle
% between the glenoid centerline and the scapular transverse axis (z axis).
% The orientation is the angle between the x axis and the glenoid
% centerline projected in the xy plane. Zero orientation correspond to -x
% axis (posterior side), 90 to superior, 180 to anterior orientaion, and
% -90 to inferior. Glenoid version is also reported as version2D (>0 for
% retroversion) and inclination (>0 for superior tilt).
versionAmplRot = vrrotvec(zAxis,centerLine); % angle between centerLine and scapular axis
versionAmpl = rad2deg(versionAmplRot(4));
xProj = dot(centerLine, xAxis);
yProj = dot(centerLine, yAxis);
zProj = dot(centerLine, zAxis);
if xProj > 0 % anterior side
if yProj > 0 % ant-sup quadran
versionOrient = pi - atan(yProj/xProj);
else % ant-inf quadran
versionOrient = -pi - atan(yProj/xProj);
end
elseif xProj < 0 % posterior side
versionOrient = - atan(yProj/xProj);
else % xProj = 0
versionOrient = pi/2 - pi * (yProj < 0);
end
versionOrient = rad2deg(versionOrient);
if zProj > 0
version = atan(xProj/zProj);
inclination = atan(yProj/zProj);
else % not realistic situation
version = NaN;
inclination = NaN;
end
version = rad2deg(version);
inclination = rad2deg(inclination);
% Angle between plane of maximum wear and CT axis [0,0,1]
maxWearPlane = cross(centerLine, zAxis);
wearPlaneAngle = vrrotvec(maxWearPlane, [0 0 1]);
wearPlaneAngle = rad2deg(wearPlaneAngle(4));
if wearPlaneAngle > 90
wearPlaneAngle = 180 - wearPlaneAngle;
end
% Glenoid center in scapula coordinate system
CT2scap = [xAxis' yAxis' zAxis'];
glenCenterLocal = (glenCenter - origin) * CT2scap;
if length(glenCenterLocal) ~= 3
warning(['Variable glenCenterLocal not defined (' SCase.id ')']);
outputArg = 0;
return
end
% Set object properties
obj.center = glenCenter;
obj.centerLocal.x = glenCenterLocal(1);
obj.centerLocal.y = glenCenterLocal(2);
obj.centerLocal.z = glenCenterLocal(3);
obj.radius = glenRadius;
obj.sphereRMSE = sphereRMSE;
obj.depth = depth;
obj.centerLine = centerLine;
obj.width = width;
obj.height = height;
obj.versionAmpl = versionAmpl;
obj.versionOrient = versionOrient;
obj.version = version;
obj.inclination = inclination;
outputArg = 1;
end
function outputArg = calcDensity(obj)
outputArg = 1;
end
end
end
c4science · Help