平面二维任意椭圆数据拟合算法推导及程序实现详解

         在刚刚过去的2017全国大学生数学建模比赛中，笔者有幸指导了一组本科学生参赛。对于赛题A《 CT系统参数标定及成像》中的CT系统参数标定，经过将问题进一步的提炼，问题最终变成了在平面二维空间中对任意椭圆进行拟合的问题，笔者花了大概四个小时的时间建立了该问题的数学公式表达、并推导出了求解该问题的算法、同时编写了实现该算法的Matlab程序。整个过程一气呵成，没有一丝的错误，当时就把我震惊了，本来以为会出很多的问题，但是当我编写完代码按下F5的那一刻，出来的结果和我之前预想的一模一样，为了更进一步的验证算法的正确性后续又做了很多的测试，结果表明算法完全正确。那一刻我真的感受到了自己平时公式推导的严谨态度和规范代码编写带来的巨大效率的提升。

公式推导

      好了，闲话不多说了，进入正题，我们通常见到的椭圆表达式是这样的：

      

    这种表达式称为椭圆的标准表达式，即其长短轴分别和坐标轴平行，且中心在坐标原点，但是更一般的椭圆其长短轴不是和坐标轴平行的，且中心也不在原点。我们将上式中的椭圆以坐标原点为旋转中心逆时针旋转θ，得到如下表达式：

展开得：

令：

则有：

此表达式表明此时的椭圆的长短轴已不再和坐标轴平行，但其中心仍在原点，接下来我们将其中心移至点(x_0,y_0)得到：

展开得：

令：

则得到椭圆的一般式：

    从这里我们就可以看到，我们只需要确定参数g, c, d, e, f就可以得到椭圆的中心(x_0,y_0)、旋转角度θ、及长短轴a和b。
    我们可以看到当知道了椭圆上的五个点的坐标后，我们就可以列出一个五元一次方程组，然后求得g, c, d, e, f。然而在实际应用中我们往往得到的数据点是大于5的因此列出的方程是一个超定方程，没有唯一解。但是我们可以求出它的最小二乘解。具体的过程如下：
    设已知椭圆上N个点的坐标的检测值(x_i,y_i) (含有检测误差),i=1,2,•••，N，将每组检测值(x_i,y_i)代入上式，则有方程误差ε_i,即：

    方程误差ε_i是由于对椭圆上点的坐标的检测误差引起的，通常N远大于5。最小二乘法原理就是用极小化方程误差的平方和来确定未知椭圆模型参数(g, c, d, e, f),即他们的极小化性能指标

由极值原则，置

即：

令：

带入上式化简得：

通过求解上述方程从而求得最小二乘意义上的参数(g, c, d, e, f)，进而通过关系式：

求得椭圆的中心(x_0,y_0)、旋转角度θ、及长短轴a和b。

Matlab仿真

%椭圆拟合测试程序

clc;
close all;
clear all

alpha = (1:360)/180*pi;  %标准椭圆参数方程角度取值
a = 15;                  %短轴长
b = 40;                  %长轴长
x0 = 5;                  %中心点x坐标
y0 = -10;                %中心点y坐标
Thita = 30/180*pi;       %旋转角度
SNR = 20;                %信噪比

%*********************************生成椭圆数据******************************
%得到标准形式的椭圆数据
x = a*cos(alpha);
y = b*sin(alpha);
%旋转角度Thita
x_r = x*cos(Thita)-y*sin(Thita);
y_r = x*sin(Thita)+y*cos(Thita);
%加入白噪声
x_r = x_r + 10^(-SNR/20)*a*randn(size(x_r));
y_r = y_r + 10^(-SNR/20)*a*randn(size(x_r));
%中心平移到（x0，y0）
x_r_s = x_r + x0;
y_r_s = y_r + y0;

[a_f,b_f,x0_f,y0_f,Thita_f]=Ellipse_fitting(x_r_s,y_r_s);    %椭圆拟合
%%
%**********利用拟合的结果生成拟合后的椭圆数据***********
x = a_f*cos(alpha);
y = b_f*sin(alpha);
x_r = x*cos(Thita_f/180*pi)-y*sin(Thita_f/180*pi);
y_r = x*sin(Thita_f/180*pi)+y*cos(Thita_f/180*pi);
x_r_s_f = x_r + x0_f;
y_r_s_f = y_r + y0_f;
%% 绘图
figure;
plot(x_r_s,y_r_s,'b',x0,y0,'b*');
xlabel('x');
ylabel('y');
axis auto equal
grid on;
hold on;
plot(x_r_s_f,y_r_s_f,'r',x0_f,y0_f,'r*');
legend('原始数据','仿真的椭圆中心','拟合后的','拟合后的椭圆中心');

fprintf('椭圆拟合结果如下(信噪比 SNR=%ddB)：\n',SNR);
fprintf('参数 a=%f,b=%f\n',a_f,b_f);
fprintf('中心坐标（x0=%f,y0=%f）\n',x0_f,y0_f);
fprintf('旋转角度 Thita=%f°\n',Thita_f);

椭圆拟合子函数

function [a,b,x0,y0,xita]=Ellipse_fitting(x,y)
%椭圆拟合子函数
%参数说明
%  a：   x轴上的轴半径
%  b：   x轴上的轴半径
%  x0:   椭圆的中心x轴坐标
%  y0:   椭圆的中心y轴坐标
%  xita：椭圆的旋转角度

N = length(x);   %求总点数

avr_x = sum(x)/N;
avr_y = sum(y)/N;
avr_xx = sum(x.*x)/N;
avr_xy = sum(x.*y)/N;
avr_yy = sum(y.*y)/N;

avr_xxx = sum(x.*x.*x)/N;
avr_xxy = sum(x.*x.*y)/N;
avr_xyy = sum(x.*y.*y)/N;
avr_yyy = sum(y.*y.*y)/N;

avr_xxxy = sum(x.*x.*x.*y)/N;
avr_xxyy = sum(x.*x.*y.*y)/N;
avr_xyyy = sum(x.*y.*y.*y)/N;
avr_yyyy = sum(y.*y.*y.*y)/N;

%求系数矩阵
AA = [avr_xxyy,avr_xyyy,avr_xxy,avr_xyy,avr_xy;...
     avr_xyyy,avr_yyyy,avr_xyy,avr_yyy,avr_yy;...
     avr_xxy, avr_xyy, avr_xx, avr_xy, avr_x;...
     avr_xyy, avr_yyy, avr_xy, avr_yy, avr_y;...
     avr_xy,  avr_yy,  avr_x,  avr_y,  1];
 
bb = [-avr_xxxy;...
      -avr_xxyy;...
      -avr_xxx;...
      -avr_xxy;...
      -avr_xx];
   
out = inv(AA)*bb;
%得到一般式的系数   
g = out(1);
c = out(2);
d = out(3);
e = out(4);
f = out(5);
%得到椭圆中心坐标
AAA = [-2, -g;...
       -g, -2*c];
bbb = [d;e];
out1 = inv(AAA)*bbb;
x0 = out1(1);
y0 = out1(2);
%得到过渡系数
A = 1/(x0^2 + c*y0^2 + g*x0*y0 - f);
B = c*A;
C = g*A;
%求得旋转角度
tan_2xita = C/(A-B);
xita = atan(tan_2xita)/2;
sin_xita = sin(xita);
cos_xita = cos(xita);
%求得轴半径
AAAA = [cos_xita^2,sin_xita^2;...
        sin_xita^2,cos_xita^2];
bbbb =[A;B];
out2 = inv(AAAA)*bbbb;
a = sqrt(1/out2(1));
b = sqrt(1/out2(2));

xita = xita/pi*180;       %转化成°

运行一下该程序得到如下结果：

matlab命令窗口输出：

可见拟合的结果与仿真设定的值很接近。