问题

给定离散点集$X_i=(x_i,y_i)$，我们希望找到最好的椭圆去拟合这些离散点。

方法

通常我们使用最小二乘法求解如下的最优化问题：

$$Min\sum _{i=1}^{N}f(x_i,E{)}^2$$
这里$f(x_i,E)$表示点$x_i$到E（指待拟合的椭圆）的最小距离。

通常我们有两种方法来表达$f(x_i,E)$，分别是：几何距离和代数距离。前面我们已经描述了基于代数距离的椭圆拟合，下面我们将重点介绍基于几何距离的椭圆拟合方法，这种方法也是椭圆拟合方法中最鲁棒、精度最高的拟合方法。

我们主要参考论文《Least-squares orthogonal distances fitting of circle,
sphere, ellipse, hyperbola, and parabola》。

论文的核心思路

显然，由于离散点到椭圆的几何距离是非线性的，因此椭圆拟合是一种非线性优化的问题，需要通过多次迭代求解。

因此下面我们阐述论文的思路时，将分成两个片段来讲解。第一段，主要介绍基于高斯-牛顿法的非线性优化方法。第二段，具体介绍这种方法在椭圆拟合中的应用。

非线性最小二乘拟合

问题

考虑一般的非线性最小二乘问题如下：
$$\underset{a}{\min }||X-F(a)|{|}^2（1）$$
这里$a\in R^q$为优化参数，F为非线性函数，$X\in R^p$是已知向量。如下我们给出基于梯度的迭代公式：
$$a_{k+1}=a_k+\lambda A{J}_F^T(X-F(a_k))（2）$$
这里$\lambda$为步长，A为缩放因子，J_F为F在当前参数a下的Jacobian矩阵。

各种优化方法不同，主要是A的选择不同。具体如下：

• 当$A=H^{-1}$时，为牛顿迭代法。
• 当$A=(J_F^T{J}_F{)}^{-1}$时，为高斯-牛顿迭代法。
• 当$A=I$时，为梯度下降法。
  这里我们选择使用高斯-牛顿迭代法。

我们将A带入（2）,即可以改写为：
$$\{\begin{array}{c}{a}_{k+1}=a_k+\lambda \Delta a(3)\\ J_F\Delta a=X-F(a_k)(4)\end{array}$$
这里$J_F=(\frac{\partial F_i}{\partial a_j}),i=1,2,...,p,j=1,2,...,q$
并且，我们根据（1）可以写出迭代算法的性能表现指数，其实就是参差的表达式：
$${\sigma }_0^2=||X-F(a)|{|}^2=[X-F(a){]}^T[X-F(a)](5)$$

求解

接下来，我们需要求解出$\Delta a$,才能进行迭代。

事实上求解$\Delta a$只需要求解方程组（4)即可.我们可以使用奇异值分解求解方程组。

基于几何距离的椭圆拟合

平面上的椭圆可以使用5个参数唯一地表达，即圆心$(X_c,Y_c)$、主轴$a,b$,角度$\alpha (-\pi /2<\alpha \le \pi /2)$.

对于椭圆的最小二乘正交距离拟合，我们将使用第一段所讲的方法。此时，我们定义
$$\hat{a}=(X_c,Y_c,a,b,\alpha {)}^T$$
$X$可以看成给定的离散点$X_i$，而自然$F(\hat{a})$就是定点$X_i$在椭圆上的最近点$X_i^{\prime }$(也称为正交关联点)。迭代公式就变成了如下：
$$\{\begin{array}{ccc}{\hat{a}}_{k+1}& =& {\hat{a}}_k+\lambda \Delta \hat{a}(6)\\ J_{X_i^{\prime },\hat{a}}\Delta \hat{a}& =& X_i-X_i^{\prime }，i=1,2....m(7)\end{array}$$
这里m指给定的离散点的个数。各个离散点均满足（7）,因此关于（7）可以联立。实际上关于（7）的方程是2m个，而未知数a的维数是5，显然2m>>5,因此解是不唯一的，我们需要求出最小二乘解即可。

显然我们必须计算$X_i^{\prime }$以及在$X_i^{\prime }$点上的Jacobian矩阵$J_{X_i^{\prime },\hat{a}}$。下面我们将阐述如何求解这两个量。

坐标系转换

为了便于后面的求解，我们需要考虑,将原来的坐标系XY，旋转$\alpha$角，建立一个位于$(0,0)$的临时坐标系xy。见Fig.3

则有
$$x=R(X-X_c)(8)$$
or
$$X=R^{-1}x+X_c(9)$$
这里
$$R=\left(\begin{array}{cc}C& S\\ -S& C\end{array}\right)和R^{-1}=R^T$$
$C=\cos (\alpha ),S=\sin (\alpha )$

椭圆上的正交关联点

经过坐标轴转换，5个椭圆参数变成了2个，仅仅包含了主轴$a,b$。椭圆上的点可以用标准方程描述如下：
$$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1(10)$$
从Fig.3上，我们可以看到位于xy坐标系上的正交关联点$（x_i^{\prime },y_i^{\prime }）$的切向量与两点（即$（x_i,y_i）、（x_i^{\prime },y_i^{\prime }）$）的连线是正交的。因此，有：
$$\frac{dy}{dx}\cdot \frac{y_i-y}{x_i-x}=\frac{-b^{2}x}{a^{2}y}\cdot \frac{y_i-y}{x_i-x}=-1(11)$$
重写（10,11）有：

$$\begin{gathered} f_1(x,y)=\frac{1}{2}(a^2{y}^2+b^2{x}^2-a^2{b}^2)=0(12) \\ {f}_2(x,y)=b^{2}x(y_i-y)-a^{2}y(x_i-x)=0(13) \end{gathered}$$
上面两式称为正交关联条件,我们将使用广义牛顿法来求解上面方程组的解。
即使用如下的迭代公式求解：
$$\{\begin{array}{c}{x}_{k+1}=x_k+\Delta x\\ Q_k\Delta x=-f(x_k)\end{array}(14)$$
$$Q=\left(\begin{array}{cc}\frac{\partial f_1}{\partial x}& \frac{\partial f_1}{\partial y}\\ \frac{\partial f_2}{\partial x}& \frac{\partial f_2}{\partial y}\end{array}\right)=\left(\begin{array}{cc}{b}^{2}x& a^{2}y\\ (a^2-b^2)y+b^2{y}_i& (a^2-b^2)x-a^2{x}_i\end{array}\right)(15)$$
对于迭代公式（14）,我们提供初始解$x_0$,如下：(见Fig.3)
$$x_0=\frac{1}{2}({\hat{x}}_{k1}+{\hat{x}}_{k2})$$
这里
$${\hat{x}}_{k1}=\left(\begin{array}{c}{x}_{k1}\\ y_{k1}\end{array}\right)=\left(\begin{array}{c}{x}_i\\ y_i\end{array}\right)ab/\sqrt{b^2{x}_i^2+a^2{y}_i}$$

一般经过3-4轮迭代就可以提供了一个很高精度的解。

我们先把$X_i$通过转换公式(8)转换成xy坐标系的$x_i$,然后通过求解广义的牛顿法求解得到正交关联点$x_i^{\prime }$,最后再通过转换公式（9）转换成XY坐标系的$X_i^{\prime }$,最后给出正交误差距离向量为：
$$X_i^{\prime \prime }=X_i-X_i^{\prime }(16)$$

椭圆上的正交关联点的Jacobian矩阵

下面我们直接给出椭圆上的正交关联点的Jacobian矩阵，（注：本部分推导比较复杂，贴出本人的详细推导过程，需要的可以参考推导1,2,3,4,5）
$$J_{X_i^{\prime },\hat{a}}=(R^{-1}{Q}^{-1}B){|}_{x=x_i^{\prime }}(17)$$
这里$Q$见（15），$B=(B_1,B_2,B_3,B_4,B_5)$
此时

椭圆的正交距离拟合

给定m个二维点，我们可以利用(16)、（17）给定的误差距离向量$X_i^{\prime \prime }$和Jacobian矩阵$J_{X_i^{\prime },\hat{a}}$,构造p(=2m)个线性方程组。如下：

,当我们使用奇异值分解求解出上述方程组的解时，就可以进行高斯-牛顿迭代。
此时我们还需要一个初值，一般可以选择基于代数拟合的椭圆或者使用圆作为初始值来进行迭代。通常我们推荐使用圆来作为初始值参数。

显然，从圆拟合中获得的初始参数为：
$(X_c,Y_c{)}_{ellipse}=(X_c,Y_c{)}_{circle},a=b=R,\alpha =0$
并且在迭代过程中，如果$a<b$，则需要令$\alpha \leftarrow \alpha -sign(\alpha )\frac{\pi }{2}$(保持a为主轴长)

标准坐标轴下的椭圆拟合

有时候，我们也需要为椭圆拟合增加一些约束，经典的就是要求在标准坐标轴下进行椭圆拟合（$\alpha =0或者\pi /2$），或者要求增加椭圆面积约束.此时，我们只需要在原来的第（7）式增加如下约束即可。

此时$w_3,w_4$为权重参数，一般可以取$1.0\times 10^6$.

实验结果

例1.

我们对Table 7的离散点进行椭圆拟合，其中初始参数$a_0$源自基于几何距离的圆拟合。
取$\lambda =1.2$，在经过19次迭代后，最终的校正向量的范数为$||\Delta a||=4.2\times 10^{-6}$,并且得到误差指数${\sigma }_0=1.1719$.见如下：

为了更好地比较，我们也考虑了初始参数$a_0$源自基于代数距离的椭圆拟合，而达到与上述相同的估计，需要进行21次迭代（$||\Delta a||=1.1\times 10^{-6}$）.
我们也比较了我们的结果与Gander的结果，而后者达到相同的估计，需要进行71次迭代（$||\Delta a||=1.0\times 10^{-6}$）.具体可参考Table 8 的III,IV.

Gander的算法有一个缺点就是不能使用圆的结果作为初始参数，为了克服这个困难，一般取$ a=R, b=R/2$.为了更加直接验证我们算法的鲁棒性，我们也考虑了使用两种设置作为初始值，分别见Table 8的II，III。即：
II：
$$a=R,b=R/2,\alpha =0$$
III:
$$a=R,b=R/2,\alpha =-1.211$$.
我们分别使用17次和23次迭代达到了相同的估计，其中校正向量的范数分别是$1.2\times 10^{-6},5.2\times 10^{-6}$.

从以上结果来看，我们的算法是鲁棒的，即使初始值不够好，最终也能迭代到较好的效果。