Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/ShoulderCase/fitSphere.m
4bbb99750a2eaterrier-plotUpdate

Last Change
Raw File

fitSphere.m
View Options

function [center,radius,residuals,R2] = fitSphere(x,y,z)
%SPHEREFIT find least squares sphere
%
% Fit a sphere to a set of xyz data points
% [center,radius,residuals] = shperefit(X)
% [center,radius,residuals] = spherefit(x,y,z);
% Input
% x,y,z Cartesian data, n x 3 matrix or three vectors (n x 1 or 1 x n)
% Output
% center: least squares sphere center coordinates, == [xc yc zc]
% radius: radius of curvature
% residuals: residuals in the radial direction
%
% Fit the equation of a sphere in Cartesian coordinates to a set of xyz
% data points by solving the overdetermined system of normal equations,
% ie, x^2 + y^2 + z^2 + a*x + b*y + c*z + d = 0
% The least squares sphere has radius R = sqrt((a^2+b^2+c^2)/4-d) and
% center coordinates (x,y,z) = (-a/2,-b/2,-c/2)
error(nargchk(1,3,nargin)); % check input arguments
if nargin == 1 % n x 3 matrix
if size(x,2) ~= 3
error ('input data must have three columns')
else
z = x(:,3); % save columns as x,y,z vectors
y = x(:,2);
x = x(:,1);
end
elseif nargin == 3 % three x,y,z vectors
x = x(:); % force into columns
y = y(:);
z = z(:);
if ~isequal(length(x),length(y),length(z)) % same length ?
error('input vectors must be same length');
end
else % must have one or three inputs
error('invalid input, n x 3 matrix or 3 n x 1 vectors expected');
end
% need four or more data points
if length(x) < 4
error('must have at least four points to fit a unique sphere');
end
% solve linear system of normal equations
A = [x y z ones(size(x))];
b = -(x.^2 + y.^2 + z.^2);
a = A \ b;
% return center coordinates and sphere radius
center = -a(1:3)./2;
radius = sqrt(sum(center.^2)-a(4));
% calculate residuals
if nargout > 2
residuals = radius - sqrt(sum(bsxfun(@minus,[x y z],center.').^2,2));
end
% whichstats = {'adjrsquare' 'rsquare'};
% stats = regstats(b,A(:,1:3),'linear',whichstats);
% R2 = stats.rsquare;
% R2adj = stats.adjrsquare;
sse = sum(residuals.^2);
meanX = [mean(x) mean(y) mean(z)];
for i=1:size(x)
X = [x(i) y(i) z(i)];
tot(i) = norm(meanX-X);
end
sst = sum(tot.^2);
R2 = 1-sse/sst;
end
c4science · Help