24小时热门版块排行榜    

北京石油化工学院2026年研究生招生接收调剂公告
查看: 4741  |  回复: 9

zhhhero

铜虫 (初入文坛)

[求助] matllab中空间圆如何拟合?

求大侠,帮忙,如何在matlab中最小二乘法拟合空间圆,求出圆心,半径。最好有matlab程序,谢谢
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

学到老,活到老
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

change0618

铁杆木虫 (著名写手)

方丈大师

【答案】应助回帖

★ ★
感谢参与,应助指数 +1
xiegangmai: 金币+2, 谢谢参与! 2013-05-02 23:10:37
可以去matlab App中心下载
CODE:
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));

CODE:
%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
  

CODE:
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);

2楼2013-05-02 09:22:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhhhero

铜虫 (初入文坛)

对了,我是有10组空间点三维坐标,现在想拟合空间圆。
学到老,活到老
3楼2013-05-02 09:58:24
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

【答案】应助回帖


感谢参与,应助指数 +1
xiegangmai: 金币+1, 谢谢参与! 2013-05-02 23:10:45
空间圆的点必然共面,将空间坐标转到平面坐标下,拟合后再转回去。
The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
4楼2013-05-02 10:23:35
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhhhero

铜虫 (初入文坛)

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

有matlab的代码吗
学到老,活到老
5楼2013-05-02 10:25:58
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

【答案】应助回帖


xiegangmai: 金币+1, 谢谢参与! 2013-05-02 23:10:53
我做过空间凸多边形的做大内切圆计算,有代码。不过空间坐标转到平面坐标这个比较简单,你可以自己写个。
The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
6楼2013-05-02 10:28:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

dbb627

荣誉版主 (著名写手)

【答案】应助回帖

★ ★
xiegangmai: 金币+2, 谢谢参与! 2013-05-02 23:11:00
CODE:
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];

The more you learn, the more you know, the more you know, and the more you forget. The more you forget, the less you know. So why bother to learn.
7楼2013-05-02 11:45:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhhhero

铜虫 (初入文坛)

引用回帖:
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 ...

对了,我是有10组空间点三维坐标,现在想拟合空间圆
学到老,活到老
8楼2013-05-02 22:22:14
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zhhhero

铜虫 (初入文坛)

引用回帖:
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 ...

大侠,给的程序是三点确定空间圆,能否对8个点,拟合空间圆给出思路呢。谢谢
学到老,活到老
9楼2013-05-02 22:49:09
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

等稻草的人

金虫 (正式写手)

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

楼主还在吗  ,现在我也在做空间圆拟合的问题,不知道你做好了没有 能教我下吗
10楼2014-06-04 16:19:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 zhhhero 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考研] 生物与医药273求调剂 +4 荔题南墙 2026-04-05 4/200 2026-04-05 20:17 by barlinike
[考研] 315求调剂 +5 &123456789 2026-04-05 5/250 2026-04-05 19:55 by nepu_uu
[考研] 材料工程302分求调剂 +8 zyx上岸! 2026-04-04 8/400 2026-04-05 18:48 by 蓝云思雨
[考研] 085410人工智能 初试316分 求调剂 +6 残星拂曙 2026-03-31 6/300 2026-04-05 12:15 by rainbow11
[考研] 270求调剂 +9 小杰pp 2026-03-31 11/550 2026-04-05 11:02 by 风雨无晴
[考研] 0703化学 +12 妮妮ninicgb 2026-04-04 13/650 2026-04-05 10:46 by 啊俊!
[考研] 272求调剂 +4 电气李 2026-04-05 4/200 2026-04-05 10:41 by lbsjt
[考研] 323求调剂(计算机视觉和大模型项目经历) +3 chaoxiicy 2026-03-31 3/150 2026-04-05 10:33 by zhq0425
[考博] 申博 +7 IQwQl 2026-04-04 7/350 2026-04-04 23:32 by mumin1990
[考研] 309分085801求调剂 +11 MY_angel 2026-03-31 11/550 2026-04-04 19:11 by 蓝云思雨
[考研] 一志愿北京科技大学材料工程085601,求调剂 +17 cdyw 2026-04-02 18/900 2026-04-04 11:14 by w_xuqing
[考研] 357求调剂 +13 1050389037 2026-04-03 13/650 2026-04-03 22:27 by 无际的草原
[考研] 22408 266求调剂 +3 masss11222 2026-04-02 3/150 2026-04-02 18:11 by 笔落锦州
[考研] 298求B区调剂 +4 zzz,,r 2026-04-02 5/250 2026-04-02 12:17 by 土木硕士招生
[考研] 一志愿厦门大学化学工程(专硕)-数二英二406分-求调剂 +5 厦大化工 2026-04-01 5/250 2026-04-02 10:03 by jp9609
[考研] 279求调剂 +6 学而思兮知 2026-04-01 6/300 2026-04-02 09:16 by vgtyfty
[考研] 298求调剂 +4 什么是胖头鱼 2026-03-30 6/300 2026-04-01 22:06 by 客尔美德
[考研] 求0861交通运输专硕or材料专硕调剂 +4 勒布朗@ 2026-03-31 4/200 2026-04-01 09:54 by 一只好果子?
[考研] 301求调剂 +8 axibli 2026-04-01 8/400 2026-04-01 09:51 by 我的船我的海
[考研] 293求调剂 +3 末未mm 2026-03-30 5/250 2026-03-30 17:23 by 王保杰33
信息提示
请填处理意见