diff --git a/ShoulderCase/@Glenoid/Glenoid.m b/ShoulderCase/@Glenoid/Glenoid.m index 1bdaee4..42820dc 100644 --- a/ShoulderCase/@Glenoid/Glenoid.m +++ b/ShoulderCase/@Glenoid/Glenoid.m @@ -1,251 +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 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.fittedSphere = Sphere(); - obj.fittedEllipsoid = Ellipsoid(); obj.SCase = SCase; obj.scapula = []; 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 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; end end end diff --git a/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m b/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m deleted file mode 100644 index b4ba8b6..0000000 --- a/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m +++ /dev/null @@ -1,138 +0,0 @@ -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 diff --git a/ShoulderCase/Ellipsoid.m b/ShoulderCase/Ellipsoid.m deleted file mode 100644 index f589838..0000000 --- a/ShoulderCase/Ellipsoid.m +++ /dev/null @@ -1,43 +0,0 @@ -classdef Ellipsoid < handle - - properties - polynomial - center - radii - axes - RMSE - end - - methods ( Access = ?Glenoid) - - function obj = Ellipsoid() - obj.polynomial = []; - obj.center = []; - obj.axes = []; - obj.R = []; - end - - function fitTo(obj,points) - obj.polynomial = ellipsoidfit_direct(points(:,1),points(:,2),points(:,3)); - [center,radii,~,R] = ellipsoid_im2ex(obj.p); - obj.center = center'; - obj.radii = radii; - obj.axes = R'; - obj.measureRMSE(points); - end - end - - methods ( Access = private, Hidden = true) - - function measureRMSE(obj,points) - [x,y,p] = ellipsoidfit_residuals(points(:,1),points(:,2),points(:,3),... - obj.center,obj.radii,obj.axes'); - type = 2; - dimension = 2; - residualsEllipsoid = vecnorm((points-[x y z]),type,dimension); - obj.RMSE = norm(residualsEllipsoid)/sqrt(length(residualsEllipsoid)); - end - - end - -end