三维坐标点拟合球（matlab and C ）

 三维坐标点读取.csv文件；拟合一个球，输出球中心和半径

%最小二乘的方法进行拟合
clear all;
close all
clc;
R = 2;         %球面半径
x0 = 100;      %球心x坐标
y0 = 1;        %球心y坐标
z0 = 76;       %球心z坐标

points = load( '2.csv');
[mm,nn]=size(points);
x=points(:,1);
y=points(:,2);
z=points(:,3);
% % 加入均值为0的高斯分布噪声 
% amp = 0.5;
% x = x + amp*rand(mm,1);
% y = y + amp*rand(mm,1);
% z = z + amp*rand(mm,1);


%球面拟合算法
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;
%计算求解线性方程的系数矩阵
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) ;  %拟合出的x坐标
y00 = resoult(2);     %拟合出的y坐标
z00 = resoult(3) ;    %拟合出的z坐标
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)   %拟合出的球半径r

C版本代码：用ransac算法，输入一系列点，输出半径和圆心；下面是2维得，可以修改为三维的，三维的圆心与半径计算可以参考四点确定一个圆的算法；

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