Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (master) Browse

R8218/ShoulderCase/@Glenoid/measureCenterWithEllipsoid.m
8f8aa39e7bd4master

Last Change
Raw File

measureCenterWithEllipsoid.m
View Options

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
c4science · Help