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 % E-mail: N/A
0038 % Created: 2010
0039 % Copyright 2010-2024 Levente Hunyadi
0040 
0041 narginchk(1,3);
0042 
0043 switch nargin  % n x 3 matrix
0044     case 1
0045         n = size(x,1);
0046         validateattributes(x, {'numeric'}, {'2d','real','size',[n,3]});
0047         z = x(:,3);
0048         y = x(:,2);
0049         x = x(:,1);
0050     otherwise  % three x,y,z vectors
0051         n = length(x(:));
0052         x = x(:);  % force into columns
0053         y = y(:);
0054         z = z(:);
0055         validateattributes(x, {'numeric'}, {'real','size',[n,1]});
0056         validateattributes(y, {'numeric'}, {'real','size',[n,1]});
0057         validateattributes(z, {'numeric'}, {'real','size',[n,1]});
0058 end
0059 
0060 % need four or more data points
0061 if n < 4
0062    error('spherefit:InsufficientData', ...
0063        'At least four points are required to fit a unique sphere.');
0064 end
0065 
0066 % solve linear system of normal equations
0067 A = [x, y, z, ones(size(x))];
0068 b = -(x.^2 + y.^2 + z.^2);
0069 a = A \ b;
0070 
0071 % return center coordinates and sphere radius
0072 center = -a(1:3)./2;
0073 radius = realsqrt(sum(center.^2)-a(4));
0074 sphere = [center' radius];
0075 
0076 % calculate residuals
0077 residuals = radius - sqrt(sum(bsxfun(@minus,[x y z],center.').^2,2));
0078 end