椭圆拟合过程中的椭圆倾角计算问题

%% 测试验证椭圆形的随机离散方法


% 按照文献中所给算法，但是感觉不准

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