matlab練習程序（螺線擬合）

這里待擬合的螺線我們選擇阿基米德螺線，對數螺線類似。

螺線的笛卡爾坐標系方程為：

 

螺線從笛卡爾坐標轉為極坐標方程為：

 

阿基米德螺線在極坐標系下極徑r和極角theta為線性關系，方程為：

 

計算步驟如下：

1.通常我們首先得到螺線在笛卡爾坐標下的一些點x,y。

2.然后根據x,y計算出r和theta。

3.最后擬合的目標就是計算出a和b，這一步可以用最小二乘。

 

擬合結果：

下圖藍色線為原始線（這里可能看不到），紅色線為擬合線，紅色點為測量點。

放大看一下：

不過有時候擬合也會失敗（這時候就可以祭出ransac大法了）：

matlab代碼如下：

clear all;
close all;
clc;

%%生成阿基米德螺線
a=6.34;
b=4.23;
theta=0:0.01:5*pi;
r = a+b*theta;

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

%%生成待擬合數據
ind = randperm(length(x),50);
dat=[x(ind)' y(ind)'] + rand(50,2)/5;
hold on;
plot(dat(:,1),dat(:,2),'r.');

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

%%因為T是周期為pi循環數列，因此需要根據不同圈數加pi
D=[R T];
D=sortrows(D);
E=D;
n = 0;
for i=2:length(D)
    if D(i,2)-D(i-1,2)<0 && D(i,2)<0
       n=n+1;
    end
    E(i,2) = E(i,2) + n*pi;   
end

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

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


%%生成對數螺線
a=1.34;
b=2.23;
theta=0:0.01:5*pi;
r = a*exp(b*theta/10);

x = r.*cos(theta);
y = r.*sin(theta);
figure;
plot(x,y)

ind = randperm(length(x),50);
dat=[x(ind)' y(ind)'] + rand(50,2)/5;
hold on;
plot(dat(:,1),dat(:,2),'r.');

最后還生成了對數螺線，大家可以自行嘗試擬合一下哈。