diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index abc30b2..84b24c6 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,242 +1,245 @@ 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]) + warning(['Could not find file: ' fileName]); + outputArg = 0; + return; + % 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 end end