Three dimensional coordinate point fitting sphere (MATLAB and C)

  3D coordinate point reading .csv file ; Fit a ball , Output ball center and radius

% Least square method for fitting 
clear all;
close all
clc;
R = 2;         % Spherical radius 
x0 = 100;      % Spheric Center x coordinate 
y0 = 1;        % Spheric Center y coordinate 
z0 = 76;       % Spheric Center z coordinate 

points = load( '2.csv');
[mm,nn]=size(points);
x=points(:,1);
y=points(:,2);
z=points(:,3);
% %  The mean value added is 0 Gaussian distribution noise  
% amp = 0.5;
% x = x + amp*rand(mm,1);
% y = y + amp*rand(mm,1);
% z = z + amp*rand(mm,1);


% Spherical fitting algorithm 
num_points = mm;
x_avr = sum(x)/num_points;
y_avr = sum(y)/num_points;
z_avr = sum(z)/num_points;

xx_avr = sum(x.*x)/num_points;
yy_avr = sum(y.*y)/num_points;
zz_avr = sum(z.*z)/num_points;
xy_avr = sum(x.*y)/num_points;
xz_avr = sum(x.*z)/num_points;
yz_avr = sum(y.*z)/num_points;

xxx_avr = sum(x.*x.*x)/num_points;
xxy_avr = sum(x.*x.*y)/num_points;
xxz_avr = sum(x.*x.*z)/num_points;
xyy_avr = sum(x.*y.*y)/num_points;
xzz_avr = sum(x.*z.*z)/num_points;
yyy_avr = sum(y.*y.*y)/num_points;
yyz_avr = sum(y.*y.*z)/num_points;
yzz_avr = sum(y.*z.*z)/num_points;
zzz_avr = sum(z.*z.*z)/num_points;
% Calculate the coefficient matrix for solving linear equations 
A = [xx_avr - x_avr*x_avr,xy_avr - x_avr*y_avr,xz_avr - x_avr*z_avr;
     xy_avr - x_avr*y_avr,yy_avr - y_avr*y_avr,yz_avr - y_avr*z_avr;
     xz_avr - x_avr*z_avr,yz_avr - y_avr*z_avr,zz_avr - z_avr*z_avr];
b = [xxx_avr - x_avr*xx_avr + xyy_avr - x_avr*yy_avr + xzz_avr - x_avr*zz_avr;
     xxy_avr - y_avr*xx_avr + yyy_avr - y_avr*yy_avr + yzz_avr - y_avr*zz_avr;
     xxz_avr - z_avr*xx_avr + yyz_avr - z_avr*yy_avr + zzz_avr - z_avr*zz_avr];
b = b/2;

resoult = A\b

x00 = resoult(1) ;  % Fitted x coordinate 
y00 = resoult(2);     % Fitted y coordinate 
z00 = resoult(3) ;    % Fitted z coordinate 
r = sqrt(xx_avr-2*x00*x_avr+x00*x00 + yy_avr-2*y00*y_avr+y00*y00 + zz_avr-2*z00*z_avr+z00*z00)   % The fitted sphere radius r

C Version code ： use ransac Algorithm , Enter a series of points , Output radius and Center ; Here is 2 Vader , It can be modified to three-dimensional , Three dimensional center and radius calculation can refer to the algorithm of four points to determine a circle ;

https://github.com/xiapengchng/Ransac-Fit-a-circle