diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index 578541f..da5b811 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,244 +1,249 @@ 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 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 = []; 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]) 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 = 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 diff --git a/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m b/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m new file mode 100644 index 0000000..b4ba8b6 --- /dev/null +++ b/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m @@ -0,0 +1,138 @@ +function output = measureCenterWithEllipsoid() + + % ** LH : Center of the glenoid by 2 other technics : + %1. Intersection sphere-MS axi, 2. Intersection ellipsoid-MS + %axis + + glenCenterOldScapulaCoord = scapulaCoord^-1*glenCenter'; % "old" center in scapula coordinates to compare with other technics + + % MS axis : Normalized axis of the sphere center and mean of glenoid points + meanGlen = mean(glenSurf); + MS = (sphereCenter-meanGlen)/norm(sphereCenter-meanGlen) ; + + %**LH : 1. Sphere intersection with MS axis + t=-(dot(MS,(meanGlen-sphereCenter))) ... + +sqrt((dot(MS,(meanGlen-sphereCenter)))^2 ... + -norm(meanGlen-sphereCenter)^2+sphereRadius^2); + % MBO : change 't' for 'sphereRadius' along with changing 'meanGlen' for 'sphereCenter' ? + glenCenterSph=meanGlen+t*MS; + + meanGlenScapulaCoord=scapulaCoord^-1*meanGlen'; + glenCenterSphScapulaCoord=scapulaCoord^-1*glenCenterSph'; % center in scapula coordinates to compare with other technics + + if glenCenterSphScapulaCoord(3)>meanGlenScapulaCoord(3) + t=-(dot(MS,(meanGlen-sphereCenter))) ... + -sqrt((dot(MS,(meanGlen-sphereCenter)))^2 ... + -norm(meanGlen-sphereCenter)^2+sphereRadius^2); + glenCenterSph=meanGlen+t*MS; + end + glenCenterSphScapulaCoord=scapulaCoord^-1*glenCenterSph'; + + % ** LH : 2. Sphere ellipso�d intersection with MS axis + p=ellipsoidfit_direct(glenSurf(:,1),glenSurf(:,2),glenSurf(:,3)); % p : analytic vector of the fit ellipsoid + try + [ellipsoidCenter,ellipsoidAxes,~,R] = ellipsoid_im2ex(p); % R : rotation matrix of the ellipsoid + catch ME + if(strcmp(ME.identifier,'ellipsoid_im2ex:InvalidArgumentValue')) + msg= ['No ellipsoid found']; + ellipsoidCenter=[NaN;NaN;NaN]; + R=[NaN,NaN,NaN;NaN,NaN,NaN;NaN,NaN,NaN]; + ellipsoidAxes=[NaN,NaN,NaN]; + end + end + ellipsoidAxesDirection=R^-1; + ellipsoidCenter=ellipsoidCenter'; % transform a vertical to horizotal vector + + % RMSE of ellipsoid fit + + [xp,yp,zp]=ellipsoidfit_residuals(glenSurf(:,1),glenSurf(:,2),glenSurf(:,3),ellipsoidCenter,ellipsoidAxes,R); + residualsEllipsoid=vecnorm((glenSurf-[xp,yp,zp])'); + ellipsoidRMSE=norm(residualsEllipsoid)/sqrt(length(residualsEllipsoid)); + + % Base matrix in the ellipsoid referential, in the order + % height-width-depth + px = []; % empty vector used to store the dot products between the ellipsoid axes and the scapula axes + py = []; + pz = []; + for i = 1:3 + px = [px abs(dot(ellipsoidAxesDirection(:,i),xAxis))]; + py = [py abs(dot(ellipsoidAxesDirection(:,i),yAxis))]; + pz = [pz abs(dot(ellipsoidAxesDirection(:,i),zAxis))]; + end + [~,a] = max(px); % Here, the width is mainly on the xAxis, + [~,b] = max(py); % height on the yAxis + [~,c] = max(pz); % and depth on the zAxis + + if a==b && a==c + a=1;b=2;c=3; + else + if a==b + if max(py)>max(px) + a=6-b-c; + else + b=6-a-c; + end + end + if a==c + if max(px)>max(pz) + c=6-a-b; + else + a=6-b-c; + end + end + if b==c + if max(py)>max(pz) + c=6-a-b; + else + b=6-a-c; + end + end + end + + NewEllipsoidAxes = [ellipsoidAxesDirection(:,b),ellipsoidAxesDirection(:,a), ... + ellipsoidAxesDirection(:,c)]; % New basis in ellipsoid referential, columns are then in the order "height-width-depth" + + % Intersection of ellipsoid and MS axis + + % ti are coefficients of the 2nd order polynom of the intersection between the ellipsoid (using analytic vector p) and MS + t2 = p(1).*MS(1).^2 + p(2).*MS(2).^2 + p(3).*MS(3).^2 + p(4).*MS(1).*MS(2) + p(5).*MS(1).*MS(3) + p(6).*MS(2).*MS(3); + + t1 = 2*p(1).*meanGlen(1).*MS(1) + 2*p(2).*meanGlen(2).*MS(2) + 2*p(3).*meanGlen(3).*MS(3) + ... + p(4).*(meanGlen(1).*MS(2) + meanGlen(2)*MS(1)) + p(5).*(meanGlen(1).*MS(3) + meanGlen(3).*MS(1)) + ... + p(6).*(meanGlen(2).*MS(3) + meanGlen(3).*MS(2)) + p(7).*MS(1) + p(8).*MS(2) + p(9).*MS(3); + + t0 = p(1).*meanGlen(1).^2 + p(2).*meanGlen(2).^2 + p(3).*meanGlen(3).^2 + ... + p(4).*meanGlen(1).*meanGlen(2) + p(5).*meanGlen(1).*meanGlen(3) + p(6).*meanGlen(2).*meanGlen(3) + ... + p(7)*meanGlen(1) + p(8).*meanGlen(2) + p(9).*meanGlen(3) + p(10); + + t=roots([t2 t1 t0]); + + if isreal(t) + if length(t)==2 + x=meanGlen(1) + t*MS(1); + y=meanGlen(2) + t*MS(2); + z=meanGlen(3) + t*MS(3); + + Center1=scapulaCoord^-1*[x(1);y(1);z(1)]; + Center2=scapulaCoord^-1*[x(2);y(2);z(2)]; + + if Center1(3)>Center2(3) %Take the intersection the farest from meanGlen in ML axis + glenCenterE=[x(2),y(2),z(2)]; + else + glenCenterE=[x(1),y(1),z(1)]; + end + + elseif length(t)==1 + x=meanGlen(1) + t*MS(1); + y=meanGlen(2) + t*MS(2); + z=meanGlen(3) + t*MS(3); + glenCenterE=[x,y,z]; + end + else + glenCenterE=[NaN,NaN,NaN]; + end + glenCenterEScapulaCoord=scapulaCoord^-1*glenCenterE'; + + %**LH : end of glenoid center + +end