function   [xc,yc,R,a] = circfit(x,y)
%
%   [xc yx R] = circfit(x,y)
%
%   fits a circle  in x,y plane in a more accurate
%   (less prone to ill condition )
%  procedure than circfit2 but using more memory
%  x,y are column vector where (x(i),y(i)) is a measured point
%
%  result is center point (yc,xc) and radius R
%  an optional output is the vector of coeficient a
% describing the circle's equation
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
%
%  By:  Izhak bucher 25/oct /1991,
    x=x(:); y=y(:);
   a=[x y ones(size(x))]\[-(x.^2+y.^2)];
   xc = -.5*a(1);
   yc = -.5*a(2);
   R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));

%try_circ_fit
%
% IB
%
% revival of a 13 years old code


  % Create data for a circle + noise
  
  th = linspace(0,2*pi,20)';
  R=1.1111111;
  sigma = R/10;
  x = R*cos(th)+randn(size(th))*sigma;
  y = R*sin(th)+randn(size(th))*sigma;
  
   plot(x,y,'o'), title(' measured points')
   pause(1)
   
   % reconstruct circle from data
   [xc,yc,Re,a] = circfit(x,y);
      xe = Re*cos(th)+xc; ye = Re*sin(th)+yc;
   
     plot(x,y,'o',[xe;xe(1)],[ye;ye(1)],'-.',R*cos(th),R*sin(th)),
     title(' measured fitted and true circles')
      legend('measured','fitted','true')
      text(xc-R*0.9,yc,sprintf('center (%g , %g );  R=%g',xc,yc,Re))
     xlabel x, ylabel y
     axis equal
  

function [center,rad,v1n,v2nb] = circlefit3d(p1,p2,p3)
% circlefit3d: Compute center and radii of circles in 3d which are defined by three points on the circumference
% usage: [center,rad,v1,v2] = circlefit3d(p1,p2,p3)
%
% arguments: (input)
%  p1, p2, p3 - vectors of points (rowwise, size(p1) = [n 3])
%               describing the three corresponding points on the same circle.
%               p1, p2 and p3 must have the same length n.
%
% arguments: (output)
%  center - (nx3) matrix of center points for each triple of points in p1,  p2, p3
%
%  rad    - (nx1) vector of circle radii.
%           if there have been errors, radii is a negative scalar ( = error code)
%
%  v1, v2 - (nx3) perpendicular vectors inside circle plane
%
% Example usage:
%
%  (1)
%      p1 = rand(10,3);
%      p2 = rand(10,3);
%      p3 = rand(10,3);
%      [center, rad] = circlefit3d(p1,p2,p3);
%      % verification, result should be all (nearly) zero
%      result(:,1)=sqrt(sum((p1-center).^2,2))-rad;
%      result(:,2)=sqrt(sum((p2-center).^2,2))-rad;
%      result(:,3)=sqrt(sum((p3-center).^2,2))-rad;
%      if sum(sum(abs(result))) < 1e-12,
%       disp('All circles have been found correctly.');
%      else,
%       disp('There had been errors.');
%      end
%
%
% (2)
%       p1=rand(4,3);p2=rand(4,3);p3=rand(4,3);
%       [center,rad,v1,v2] = circlefit3d(p1,p2,p3);
%       plot3(p1(:,1),p1(:,2),p1(:,3),'bo');hold on;plot3(p2(:,1),p2(:,2),p2(:,3),'bo');plot3(p3(:,1),p3(:,2),p3(:,3),'bo');
%       for i=1:361,
%           a = i/180*pi;
%           x = center(:,1)+sin(a)*rad.*v1(:,1)+cos(a)*rad.*v2(:,1);
%           y = center(:,2)+sin(a)*rad.*v1(:,2)+cos(a)*rad.*v2(:,2);
%           z = center(:,3)+sin(a)*rad.*v1(:,3)+cos(a)*rad.*v2(:,3);
%           plot3(x,y,z,'r.');
%       end
%       axis equal;grid on;rotate3d on;
%
%
% Author: Johannes Korsawe
% E-mail: johannes.korsawe@volkswagen.de
% Release: 1.0
% Release date: 26/01/2012

% Default values
center = [];rad = 0;v1n=[];v2nb=[];

% check inputs
% check number of inputs
if nargin~=3,
    fprintf('??? Error using ==> cirlefit3d\nThree input matrices are needed.\n');rad = -1;return;
end
% check size of inputs
if size(p1,2)~=3 || size(p2,2)~=3 || size(p3,2)~=3,
    fprintf('??? Error using ==> cirlefit3d\nAll input matrices must have three columns.\n');rad = -2;return;
end
n = size(p1,1);
if size(p2,1)~=n || size(p3,1)~=n,
    fprintf('??? Error using ==> cirlefit3d\nAll input matrices must have the same number or rows.\n');rad = -3;return;
end
% more checks are to follow inside calculation

% Start calculation
% v1, v2 describe the vectors from p1 to p2 and p3, resp.
v1 = p2 - p1;v2 = p3 - p1;
% l1, l2 describe the lengths of those vectors
l1 = sqrt((v1(:,1).*v1(:,1)+v1(:,2).*v1(:,2)+v1(:,3).*v1(:,3)));
l2 = sqrt((v2(:,1).*v2(:,1)+v2(:,2).*v2(:,2)+v2(:,3).*v2(:,3)));
if find(l1==0) | find(l2==0), %#ok<OR2>
    fprintf('??? Error using ==> cirlefit3d\nCorresponding input points must not be identical.\n');rad = -4;return;
end
% v1n, v2n describe the normalized vectors v1 and v2
v1n = v1;for i=1:3, v1n(:,i) = v1n(:,i)./l1;end
v2n = v2;for i=1:3, v2n(:,i) = v2n(:,i)./l2;end
% nv describes the normal vector on the plane of the circle
nv = [v1n(:,2).*v2n(:,3) - v1n(:,3).*v2n(:,2) , v1n(:,3).*v2n(:,1) - v1n(:,1).*v2n(:,3) , v1n(:,1).*v2n(:,2) - v1n(:,2).*v2n(:,1)];
if find(sum(abs(nv),2)<1e-5),
    fprintf('??? Warning using ==> cirlefit3d\nSome corresponding input points are nearly collinear.\n');
end
% v2nb: orthogonalization of v2n against v1n
dotp = v2n(:,1).*v1n(:,1) + v2n(:,2).*v1n(:,2) + v2n(:,3).*v1n(:,3);
v2nb = v2n;for i=1:3,v2nb(:,i) = v2nb(:,i) - dotp.*v1n(:,i);end
% normalize v2nb
l2nb = sqrt((v2nb(:,1).*v2nb(:,1)+v2nb(:,2).*v2nb(:,2)+v2nb(:,3).*v2nb(:,3)));
for i=1:3, v2nb(:,i) = v2nb(:,i)./l2nb;end

% remark: the circle plane will now be discretized as follows
%
% origin: p1                    normal vector on plane: nv
% first coordinate vector: v1n  second coordinate vector: v2nb

% calculate 2d coordinates of points in each plane
% p1_2d = zeros(n,2); % set per construction
% p2_2d = zeros(n,2);p2_2d(:,1) = l1; % set per construction
p3_2d = zeros(n,2); % has to be calculated
for i = 1:3,
    p3_2d(:,1) = p3_2d(:,1) + v2(:,i).*v1n(:,i);
    p3_2d(:,2) = p3_2d(:,2) + v2(:,i).*v2nb(:,i);
end

% calculate the fitting circle
% due to the special construction of the 2d system this boils down to solving
% q1 = [0,0], q2 = [a,0], q3 = [b,c] (points on 2d circle)
% crossing perpendicular bisectors, s and t running indices:
% solve [a/2,s] = [b/2 + c*t, c/2 - b*t]
% solution t = (a-b)/(2*c)

a = l1;b = p3_2d(:,1);c = p3_2d(:,2);
t = 0.5*(a-b)./c;
scale1 = b/2 + c.*t;scale2 = c/2 - b.*t;

% centers
center = zeros(n,3);
for i=1:3,
    center(:,i) = p1(:,i) + scale1.*v1n(:,i) + scale2.*v2nb(:,i);
end

% radii
rad = sqrt((center(:,1)-p1(:,1)).^2+(center(:,2)-p1(:,2)).^2+(center(:,3)-p1(:,3)).^2);

4楼: Originally posted by dbb627 at 2013-05-02 10:23:35
空间圆的点必然共面，将空间坐标转到平面坐标下，拟合后再转回去。

function [T,dp]=Coordinate3DTranfer(A,B,C)
%% 函数使用说明
%以不共线的三点ABC（xyz坐标）确定平面为ABC，建立新平面坐标XYZ
%新平面坐标XYZ 的原点是xyz坐标原点O到平面ABC垂足P
% new OZ 方向余弦(OP) 即ABC 平面法向量（x y z）
% new OX 方向余弦(PA)
% new OY 方向余弦     即OPA 平面法向量(ly,my,ny)
%
% 将空间平面xoy（xyz坐标）变换为new XOY平面（XYZ坐标）（Z=inv(T)*dp'）
% inv(T)*point
% [X1 X2 X3]             [x1 x2 x3]
% [Y1 Y2 Y3]=inv(T)*[y1 y2 y3]
% [Z   Z   Z ]              [z1 z2  z3]
% 逆变换,new XOY平面上2D坐标为X Y时Z=Z(3),变换到原xoy (xyz)
% T*point
% [x1 x2 x3]      [X1 X2 X3]
% [y1 y2 y3]=T*[Y1 Y2 Y3]
% [z1 z2 z3]       [Z  Z  Z ]
%% Example
% A=[1 3 2];B=[4 2 3];C=[2 1 4];
% [T,dp]=Coordinate3DTranfer(A,B,C)
% Z=inv(T)*dp'
% point=(inv(T)*[A' B' C'])'% 变换
% T*point'                  

x1=A(:,1);y1=A(:,2);z1=A(:,3);
x2=B(:,1);y2=B(:,2);z2=B(:,3);
x3=C(:,1);y3=C(:,2);z3=C(:,3);
%% ABC 平面法向量（x y z）
x=(y1.*z2 - y2.*z1 - y1.*z3 + y3.*z1 + y2.*z3 - y3.*z2);
y =-(x1.*z2 - x2.*z1 - x1.*z3 + x3.*z1 + x2.*z3 - x3.*z2);
z =(x1.*y2 - x2.*y1 - x1.*y3 + x3.*y1 + x2.*y3 - x3.*y2);
%% 平面方程 x(X-x1)+y(Y-Y1)+z(Z-z1)=0 O到平面ABC垂足P（xp yp zp）平移参数
%[xp,yp,zp]=solve('xp*(x2-x1)+yp*(y2-y1)+zp*(z2-z1)=0','xp*(x3-x1)+yp*(y3-y1)+zp*(z3-z1)=0','x*(xp-x1)+y*(yp-y1)+z*(zp-z1)','xp,yp,zp')
PP=(x.*y1.*z2 - x.*y2.*z1 - x1.*y.*z2 + x1.*y2.*z + x2.*y.*z1 - x2.*y1.*z - x.*y1.*z3 +...
    x.*y3.*z1 + x1.*y.*z3 - x1.*y3.*z - x3.*y.*z1 + x3.*y1.*z + x.*y2.*z3 - x.*y3.*z2 - x2.*y.*z3 +...
    x2.*y3.*z + x3.*y.*z2 - x3.*y2.*z);
xp =(y.*y1.^2.*z2 - y.*y1.^2.*z3 - y2.*z.*z1.^2 + y3.*z.*z1.^2 + x.*x1.*y1.*z2 - x.*x1.*y2.*z1 -...
    x.*x1.*y1.*z3 + x.*x1.*y3.*z1 + x.*x1.*y2.*z3 - x.*x1.*y3.*z2 - y.*y1.*y2.*z1 + y.*y1.*y3.*z1 +...
    y.*y1.*y2.*z3 - y.*y1.*y3.*z2 + y1.*z.*z1.*z2 - y1.*z.*z1.*z3 + y2.*z.*z1.*z3 - y3.*z.*z1.*z2)./PP;
yp =-((x.*x1 + y.*y1 + z.*z1).*(x1.*z2 - x2.*z1 - x1.*z3 + x3.*z1 + x2.*z3 - x3.*z2))./PP;
zp =((x.*x1 + y.*y1 + z.*z1).*(x1.*y2 - x2.*y1 - x1.*y3 + x3.*y1 + x2.*y3 - x3.*y2))./PP;
%% new OZ 方向余弦(OP) ABC 平面法向量（x y z）
lz=x./sqrt(x.^2+y.^2+z.^2);mz=y./sqrt(x.^2+y.^2+z.^2);nz=z./sqrt(x.^2+y.^2+z.^2);
%% new Ox 方向余弦(PA)
d=sqrt((x1-xp).^2+(y1-yp).^2+(z1-zp).^2);
lx=(xp-x1)./d;mx=(yp-y1)./d;nx=(zp-z1)./d;
%% new Oy 方向余弦  OPA 平面法向量(ly,my,ny)
xx=(y1.*zp - yp.*z1);
yy =-(x1.*zp - xp.*z1);
zz =(x1.*yp - xp.*y1);
ly=xx./sqrt(xx.^2+yy.^2+zz.^2);my=yy./sqrt(xx.^2+yy.^2+zz.^2);ny=zz./sqrt(xx.^2+yy.^2+zz.^2);
T=[lx ly lz;mx my mz;nx ny nz];
dp=[xp yp zp];

2楼: Originally posted by change0618 at 2013-05-02 09:22:48
可以去matlab App中心下载


function    = circfit(x,y)
%
%    = circfit(x,y)
%
%   fits a circle  in x,y plane in a more accurate
%   (less prone to ill condition )
%  procedure than circfit2 ...

2楼: Originally posted by change0618 at 2013-05-02 09:22:48
可以去matlab App中心下载


function    = circfit(x,y)
%
%    = circfit(x,y)
%
%   fits a circle  in x,y plane in a more accurate
%   (less prone to ill condition )
%  procedure than circfit2 ...

3楼: Originally posted by zhhhero at 2013-05-02 09:58:24
对了，我是有10组空间点三维坐标，现在想拟合空间圆。