设为首页收藏本站

EPS数据狗论坛

 找回密码
 立即注册

QQ登录

只需一步,快速开始

查看: 2418|回复: 0

matlab练习程序(螺线拟合)

[复制链接]

10

主题

101

金钱

149

积分

入门用户

发表于 2019-12-16 15:56:59 | 显示全部楼层 |阅读模式

这里待拟合的螺线我们选择阿基米德螺线,对数螺线类似。

螺线的笛卡尔坐标系方程为:
螺线从笛卡尔坐标转为极坐标方程为:
1.png
2.png

阿基米德螺线在极坐标系下极径r和极角theta为线性关系,方程为:

计算步骤如下:

1.通常我们首先得到螺线在笛卡尔坐标下的一些点x,y。

2.然后根据x,y计算出r和theta。

3.最后拟合的目标就是计算出a和b,这一步可以用最小二乘。

拟合结果:

下图蓝色线为原始线(这里可能看不到),红色线为拟合线,红色点为测量点。
3.png

放大看一下:
4.png
不过有时候拟合也会失败(这时候就可以祭出ransac大法了):
5.png
matlab代码如下:
  1. clear all;
  2. close all;
  3. clc;

  4. %%生成阿基米德螺线
  5. a=6.34;
  6. b=4.23;
  7. theta=0:0.01:5*pi;
  8. r = a+b*theta;

  9. x = r.*cos(theta);
  10. y = r.*sin(theta);
  11. plot(x,y,'b')

  12. %%生成待拟合数据
  13. ind = randperm(length(x),50);
  14. dat=[x(ind)' y(ind)'] + rand(50,2)/5;
  15. hold on;
  16. plot(dat(:,1),dat(:,2),'r.');

  17. T = atan(dat(:,2)./dat(:,1));
  18. R = sqrt(dat(:,1).^2+dat(:,2).^2);

  19. %%因为T是周期为pi循环数列,因此需要根据不同圈数加pi
  20. D=[R T];
  21. D=sortrows(D);
  22. E=D;
  23. n = 0;
  24. for i=2:length(D)
  25.     if D(i,2)-D(i-1,2)<0 && D(i,2)<0
  26.        n=n+1;
  27.     end
  28.     E(i,2) = E(i,2) + n*pi;   
  29. end

  30. X = [E(:,2) ones(length(dat),1)];
  31. Y = E(:,1);
  32. C = inv(X'*X)*X'*Y;     

  33. theta=0:0.01:5*pi;
  34. r = C(2)+C(1)*theta;
  35. x = r.*cos(theta);
  36. y = r.*sin(theta);
  37. plot(x,y,'r')


  38. %%生成对数螺线
  39. a=1.34;
  40. b=2.23;
  41. theta=0:0.01:5*pi;
  42. r = a*exp(b*theta/10);

  43. x = r.*cos(theta);
  44. y = r.*sin(theta);
  45. figure;
  46. plot(x,y)

  47. ind = randperm(length(x),50);
  48. dat=[x(ind)' y(ind)'] + rand(50,2)/5;
  49. hold on;
  50. plot(dat(:,1),dat(:,2),'r.');
复制代码


最后还生成了对数螺线,大家可以自行尝试拟合一下哈。
6.png
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

关闭

站长推荐上一条 /1 下一条

客服中心
关闭
在线时间:
周一~周五
8:30-17:30
QQ群:
653541906
联系电话:
010-85786021-8017
在线咨询
客服中心

意见反馈|网站地图|手机版|小黑屋|EPS数据狗论坛 ( 京ICP备09019565号-3 )   

Powered by BFIT! X3.4

© 2008-2028 BFIT Inc.

快速回复 返回顶部 返回列表