fitSphere

FITSPHERE Fit a sphere to 3D points using the least squares approach.

   SPHERE = fitSphere(PTS)
   Fits the equation of a sphere in Cartesian coordinates to the N-by-3 
   array PTS using the least squares approach. The sphere is represented 
   by its center [xc yc zc] and its radius r: SPHERE = [xc yc zc r].

   SPHERE = fitSphere(X, Y, Z)
   Use three vectors X, Y and Z with the length N instead of a the 
   N-by-3 array PTS.
  
   [SPHERE, RESIDUALS] = fitSphere(...)
   Additionally outputs the residuals in the radial direction.

   Example:
     center=-100 + 200*rand(1,3);
     radius = randi([10 100]);
     [x,y,z]=drawSphere(center, radius);
     x=x+rand(size(x)); y=y+rand(size(y)); z=z+rand(size(z));
     sampleIdx = randi(numel(x),[1,randi([4, numel(x)])]);
     x=x(sampleIdx); y=y(sampleIdx); z=z(sampleIdx);
     sphere = fitSphere(x,y,z);
     figure('color','w'); hold on; axis equal tight; view(3)
     drawPoint3d(x,y,z)
     drawSphere(sphere,'FaceAlpha',0.5)
   
   See also:
     createSphere, drawSphere, intersectLineSphere, intersectPlaneSphere

   Source:
     Levente Hunyadi - Fitting quadratic curves and surfaces:
     https://de.mathworks.com/matlabcentral/fileexchange/45356

0001 function [sphere, residuals] = fitSphere(x,y,z)
0002 %FITSPHERE Fit a sphere to 3D points using the least squares approach.
0003 %
0004 %   SPHERE = fitSphere(PTS)
0005 %   Fits the equation of a sphere in Cartesian coordinates to the N-by-3
0006 %   array PTS using the least squares approach. The sphere is represented
0007 %   by its center [xc yc zc] and its radius r: SPHERE = [xc yc zc r].
0008 %
0009 %   SPHERE = fitSphere(X, Y, Z)
0010 %   Use three vectors X, Y and Z with the length N instead of a the
0011 %   N-by-3 array PTS.
0012 %
0013 %   [SPHERE, RESIDUALS] = fitSphere(...)
0014 %   Additionally outputs the residuals in the radial direction.
0015 %
0016 %   Example:
0017 %     center=-100 + 200*rand(1,3);
0018 %     radius = randi([10 100]);
0019 %     [x,y,z]=drawSphere(center, radius);
0020 %     x=x+rand(size(x)); y=y+rand(size(y)); z=z+rand(size(z));
0021 %     sampleIdx = randi(numel(x),[1,randi([4, numel(x)])]);
0022 %     x=x(sampleIdx); y=y(sampleIdx); z=z(sampleIdx);
0023 %     sphere = fitSphere(x,y,z);
0024 %     figure('color','w'); hold on; axis equal tight; view(3)
0025 %     drawPoint3d(x,y,z)
0026 %     drawSphere(sphere,'FaceAlpha',0.5)
0027 %
0028 %   See also:
0029 %     createSphere, drawSphere, intersectLineSphere, intersectPlaneSphere
0030 %
0031 %   Source:
0032 %     Levente Hunyadi - Fitting quadratic curves and surfaces:
0033 %     https://de.mathworks.com/matlabcentral/fileexchange/45356
0034 
0035 % ------
0036 % Author: Levente Hunyadi, oqilipo (minor adaptions for matGeom)
0037 % Created: 2010
0038 % Copyright 2010 Levente Hunyadi
0039 
0040 narginchk(1,3);
0041 
0042 switch nargin  % n x 3 matrix
0043     case 1
0044         n = size(x,1);
0045         validateattributes(x, {'numeric'}, {'2d','real','size',[n,3]});
0046         z = x(:,3);
0047         y = x(:,2);
0048         x = x(:,1);
0049     otherwise  % three x,y,z vectors
0050         n = length(x(:));
0051         x = x(:);  % force into columns
0052         y = y(:);
0053         z = z(:);
0054         validateattributes(x, {'numeric'}, {'real','size',[n,1]});
0055         validateattributes(y, {'numeric'}, {'real','size',[n,1]});
0056         validateattributes(z, {'numeric'}, {'real','size',[n,1]});
0057 end
0058 
0059 % need four or more data points
0060 if n < 4
0061    error('spherefit:InsufficientData', ...
0062        'At least four points are required to fit a unique sphere.');
0063 end
0064 
0065 % solve linear system of normal equations
0066 A = [x, y, z, ones(size(x))];
0067 b = -(x.^2 + y.^2 + z.^2);
0068 a = A \ b;
0069 
0070 % return center coordinates and sphere radius
0071 center = -a(1:3)./2;
0072 radius = realsqrt(sum(center.^2)-a(4));
0073 sphere = [center' radius];
0074 
0075 % calculate residuals
0076 residuals = radius - sqrt(sum(bsxfun(@minus,[x y z],center.').^2,2));
0077 end