diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index da5b811..1bdaee4 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,249 +1,251 @@ 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 + fittedSphere + fittedEllipsoid end properties (Hidden = true) SCase scapula 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; - obj.scapula = []; + 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.fittedSphere = Sphere(); + obj.fittedEllipsoid = Ellipsoid(); + + obj.SCase = SCase; + obj.scapula = []; end end - methods (Access = private, Hidden = true) - output = measureGlenoidCenterWithEllipsoid(); - 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 - warning(['Could not find file: ' fileName]); - outputArg = 0; - return; - % error(['Could not find file: ' fileName]) + [points,faces,~]=loadStl(fileName,1); + obj.surface.points = points; + obj.surface.faces = faces; + else + 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.meanPoint = mean(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 + outputArg = 1; + SCase = obj.SCase; glenSurf = obj.surface.points; % Get local coordinate system from parent scapula object origin = obj.scapula.coordSys.origin; xAxis = obj.scapula.coordSys.PA; yAxis = obj.scapula.coordSys.IS; zAxis = obj.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 = -atan(xProj/zProj); inclination = -atan(yProj/zProj); 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