橢圓擬合過程中的橢圓傾角計算問題

　　matlab離散點擬合二維橢圓的時候，主要有兩種方法，一種是使用boundary函數拾取邊界點，直接擬合，一共五個參數，使用matlab的regress就可以，公式如下：

　　但是這種方法有一些缺點，就是在由於參數過多，擬合垂直情況和系數相差過大的異常橢圓情況會誤差會比較大，因此一般使用第二種方法，即將原始的離散點通過使用旋轉和平移變化移動到長軸和x軸重合，以原點為中心，最后通過平面變換回去，由於少了b,d,e三項，擬合的精度會好很多，可以解決一些過於扁平的離散點的橢圓擬合，但是這里面會遇到一個很重要的問題就是如何計算旋轉角，即通過離散點計算橢圓軸的傾角，由於不是一個一一對應的函數，所以計算起來會有些困難，這里參考http://mathworld.wolfram.com/Ellipse.html的計算方法，公式如下

 

代碼如下：

%% 測試驗證橢圓形的隨機離散方法


% 按照文獻中所給算法，但是感覺不准

for anglex=0:pi/90:2*pi
% 常規方法生成
a=20;b=5;
num=400;
x1=-a+2*a*rand(num,1);
y1=-b+2*b*rand(num,1);
xs=[];ys=[];

for i=1:length(x1)
    if (x1(i)^2/a^2+y1(i)^2/b^2)<=1
        xs=[xs;x1(i)];ys=[ys;y1(i)];
        hold on
    end
end
plot(xs,ys,'y*');
%% 生成經過平移和旋轉的離散點
deltax=3;deltay=6;tx=anglex;
pingyi=[1 0 deltax;0 1 deltay;0 0 1];
xuanzhuan=[cos(tx) -sin(tx) 0;sin(tx) cos(tx) 0;0 0 1]; %逆時針旋轉一定的角度
rota_xsys=pingyi*xuanzhuan*[xs';ys';ones(size(xs'))];
xs=rota_xsys(1,:);ys=rota_xsys(2,:);

% scatter(x1,y1,'filled');
 plot(xs,ys,'b*');
 axis equal
 hold on


%% 橢圓點邊界點擬合--顯然更加准確
kp=boundary(xs',ys',0.1);
xs1=[xs(kp)]';ys1=[ys(kp)]';
plot(xs1,ys1,'r');
[pxf,~,rp]=regress(-ones(size(xs1)),[xs1.^2 xs1.*ys1 ys1.^2 xs1 ys1]);
a=pxf(1);b=pxf(2);c=pxf(3);d=pxf(4);e=pxf(5);
xc=(b*e-2*c*d)/(4*a*c-b^2);yc=(b*d-2*a*e)/(4*a*c-b^2);  %中心
la=sqrt(2*(a*xc^2+c*yc^2+b*xc*yc-1)/(a+c+sqrt((a-c)^2+b^2))); %長半軸
lb=sqrt(2*(a*xc^2+c*yc^2+b*xc*yc-1)/(a+c-sqrt((a-c)^2+b^2))) ;%短半軸
if la<=lb
    lx=la;
    la=lb;
    lb=lx;
end

theta=linspace(0,2*pi,100);
xb=la*sin(theta);yb=lb*cos(theta);
plot(xb,yb,'g*');
%% 求橢圓的傾角
if b==0 && abs(a)<abs(c)
    rota=0;
    disp(['角度1=' num2str(rota)])
elseif b==0 && abs(a)>abs(c)
    rota=pi/2;
    disp(['角度2=' num2str(rota)])
elseif b~=0 && abs(a)<abs(c)
    ex=1/2*acot((a-c)/b);
    rota=ex;
    rotax=rota*180/pi;
    disp(['角度3=' num2str(rotax)])
elseif b~=0 && abs(a)>abs(c)
    ex=1/2*acot((a-c)/b);
    rota=pi/2+ex;
    rotax=rota*180/pi;
    disp(['角度4=' num2str(rotax)])
end




% disp(['角度' num2str(rotax)])
xuanzhuan=[cos(rota) -sin(rota) 0;sin(rota) cos(rota) 0;0 0 1];
pingyi=[1 0 xc;0 1 yc;0 0 1];

rota_xy=pingyi*xuanzhuan*[xb; yb;ones(size(xb))];
% rota_xy=[xc 1;1 yc]*[cos(rota) sin(rota);-sin(rota) cos(rota)]*[xb;yb]; %旋轉平移過去
xf1=rota_xy(1,:);yf1=rota_xy(2,:);
% zf1=(df-af.*xf1-bf.*yf1)./cf;
% fs=[xf1' yf1' zf1'];
hold on
plot(xf1,yf1,'k-');
title(['角度=' num2str(anglex*180/pi)]);

pause(1);
drawnow
clf

end

　　測試了360個角度的計算，擬合的橢圓和離散點符合地都比較好。