基于代数距离的椭圆拟合

问题

给定离散点集$X_i=(x_i,y_i),i=1,2,...N$，我们希望找到误差最小的椭圆去拟合这些离散点。

方法

由于椭圆的形式可以给定， 自然我们将使用最小二乘法来求解椭圆。主要依据论文《Direct least squares fitting of ellipsees, Fitzgibbon, Pilu and Fischer in Fitzgibbon, A.W., Pilu, M., and Fischer R.B.,Proc. of the 13th Internation Conference on Pattern Recognition, pp 253–257, Vienna, 1996》。

椭圆拟合的难点

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

$$Min\sum _{i=1}^{N}f(x_i,E)$$
这里$f(x_i,E)$表示点$x_i$到E（指待拟合的椭圆）的最小距离。一般称为几何距离。但是我们很难直接给定几何距离的解析表达，因此很难求出。因此我们退而求其次，我们采用：为椭圆写下隐式方程，然后将点的坐标带入此隐式方程就得到了点到椭圆的距离。这种方法对于直线拟合、圆拟合，返回的就是到其的真实距离。而对于椭圆拟合，它返回的值是与距离有类似的属性，但不是一个真正的距离值。因此这个距离被称为代数距离。
椭圆可以用下面的隐式方程表达：

$$a_1{x}^2+a_{2}xy+a_3{y}^2+a_{4}x+a_{5}y+a_6=0$$

与直线类似（直线表达为：$a_{1}x+a_{2}y+a_3=0$）,上式的一组参数也是齐次量，即只能定义到一个比例因子。而且表达式也能表达双曲线和抛物线。而椭圆
通常要求:$a_2^2-4a_1{a}_4<0$。最好通过令$a_2^2-4a_1{a}_4=-1$,即可同时解决这两个问题。

请参考维基百科椭圆定义和wolfram MathWorld 椭圆定义

一般的求解过程

对上面的椭圆方程改写为：

$$f(a,(x,y))=D\cdot a=0$$

这里$D=(x^2,xy,y^2,x,y,1)$,$a=(a_1,a_2,a_3,a_4,a_5,a_6)$.
于是给定N个离散点$X_i=(x_i,y_i),i=1...N$,通过最小化代数距离：

$$Min\Delta (a,X)=\sum _{i=1}^N(f(a,X_i){)}^2$$

重写上式，即有：

$$Min\Delta (a,X)=\sum _{i=1}^N{a}^T{D}_i^T{D}_{i}a=a^{T}Sa$$

这里$S=\sum D_i^T{D}_i$为6x6矩阵。
另外，我们可以把各个$D_i$合并起来，则有
$$\hat{D}=\left(\begin{array}{c}{D}_1\\ D_2\\ ⋮\\ D_N\end{array}\right)$$
因此：
$$S={\hat{D}}^T\hat{D}$$

此外，我们还有约束条件：
$a_2^2-4a_1{a}_4<0$，可以改写为：

$$a^{T}Ca=-1$$

其中$C\in R^{6\times 6}$,且大部分为0，只有$C_{1,3}=C_{3,1}=-2,C_{2,2}=1$。
因此综合上来即为，求解如下的最优化问题：

$$Min{a}^{T}Sa$$
$$s.t.a^{T}Ca=-1$$

公式和拉格朗日乘子法

通过使用拉格朗日乘子法 ，我们可以得到

$$L(a)=a^{T}Sa-\lambda (a^{T}Ca+1)$$

然后求解${\partial }_{a}L(a)=0$可以求解得到：
$$Sa=\lambda Ca$$

上式是经典的广义特征值问题。我们可以直接求解，其中$\lambda$为广义特征值，而a为广义特征向量
论文 中指出有且仅有一个负的特征值，且其对应的特征向量即为我们需要的椭圆方程的系数，即这里的a.

一般的圆锥方程到标准椭圆方程的转化。

当我们求解得到了圆锥曲线的系数，即a,我们一般需要获得椭圆的标准方程，也就是获得椭圆的如下参数：
中心$(c_x,c_y)$、长短半轴$r_1,r_2$,旋转角度$\theta$.

这里我们可以给出结果，也可以自己结合椭圆的标准形式与圆锥曲线的方程去推导.

参考文献：http://mathworld.wolfram.com/Ellipse.html.
##编程实现
下面描述matlab与c++实现。

matlab

% fitellip gives the 6 parameter vector of the algebraic circle fit
% to a(1)x^2 + a(2)xy + a(3)y^2 + a(4)x + a(5)y + a(6) = 0
% X & Y are lists of point coordinates and must be column vectors.(或者行向量)
function a = fitellip(X,Y)

   % normalize data
   mx = mean(X);
   my = mean(Y);
   sx = (max(X)-min(X))/2;
   sy = (max(Y)-min(Y))/2;
   x = (X-mx)/sx;
   y = (Y-my)/sy;
   % Force to column vectors
   x = x(:);
   y = y(:);

   % Build design matrix
   D = [ x.*x  x.*y  y.*y  x  y  ones(size(x)) ];

   % Build scatter matrix
   S = D'*D;

   % Build 6x6 constraint matrix
   C(6,6) = 0; C(1,3) = -2; C(2,2) = 1; C(3,1) = -2;

   % Solve eigensystem
   [gevec, geval] = eig(S,C);

   % Find the negative eigenvalue
   [NegR, NegC] = find(geval < 0 & ~isinf(geval));

   % Extract eigenvector corresponding to positive eigenvalue
   A = gevec(:,NegC);

   % unnormalize
   a = [
        A(1)*sy*sy,   ...
        A(2)*sx*sy,   ...
        A(3)*sx*sx,   ...
        -2*A(1)*sy*sy*mx - A(2)*sx*sy*my + A(4)*sx*sy*sy,   ...
        -A(2)*sx*sy*mx - 2*A(3)*sx*sx*my + A(5)*sx*sx*sy,   ...
        A(1)*sy*sy*mx*mx + A(2)*sx*sy*mx*my + A(3)*sx*sx*my*my   ...
                - A(4)*sx*sy*sy*mx - A(5)*sx*sx*sy*my   ...
                + A(6)*sx*sx*sy*sy   ...
       ]';

c++

我们参考了网上的版本，并修复了他的问题，最终也做出了和matlab一样的效果。其中最关键的是广义特征值的求解。我们使用了Eigen与clapack库。其中Eigen易于表达矩阵，和matlab用法类似，是个强大的C++线性代数库。而CLAPACK是线性代数包Lapack面向C/c++的接口。里面包含了很丰富的线性代数算法，包括广义特征值求解接口，而且速度很快。我们希望将二者结合起来使用。

Eigen的安装

Eigen直接以源代码的方式提供给用户，因此我们从官网上下载下后，直接在工程中包含其头文件路径即可。具体可参考：http://blog.csdn.net/abcjennifer/article/details/7781936

clapack的安装

请查看官网，里面包含了详细的使用与安装步骤。

也可以使用我们已经编译了的vc2010和vc2013的库，可以点击下载。

尽管clapack面向c语言，因此需要我们在包含头文件的时候，记得加上extern “C”.但是最新的版本（比如CLAPACK 3.2.1）已经为我们在头文件中加上了这些限制符，因此最新的版本可以兼容c和c++，所以直接在项目包含头文件即可。

比如像下面一样：

//Eigen
#include <Eigen/Dense>
#include <Eigen/Core>
#include <iostream>

//clapack,必须放在Eigen后面

#include <f2c.h>
#include <clapack.h>

而且应该注意Eigen与CLAPACK混合使用的时候，CLAPACK的头文件要加在Eigen的后面。否则会出错。

代码

请查看Github主页：https://github.com/xiamenwcy/EllipseFitting

实验结果

参考文献

Eigen、LAPACK

1. Interfacing Eigen with LAPACK
2. Using BLAS/LAPACK from Eigen
3. CLAPACK for Windows
4. 如何在VC中调用CLAPACK
5. 使用Lapack库的一些重要规则
6. 在VS中用CLAPACK解决广义特征值问题
7. LAPACK 3.7.0
8. C＋＋编程：C++中同时使用Eigen和CLAPACK
9. CLAPACK在vc++6.0中成功调用
10. CLAPACK动态调用
11. Leading dimension
12. leading dimension 2

椭圆拟合

• Fitting an Ellipse to a Set of Data Points(python实现)
• 二次曲线Quadratic Curve
• fitellipse.cpp(opencv)
• Creating Bounding rotated boxes and ellipses for contours
• OpenCV’s fitEllipse() sometimes returns completely wrong ellipses
• Fitting an Ellipse to a Set of Data Points
• Direct Least Square Fitting of Ellipses (论文)
• best-fitting-line-circle-and-ellipse
• opencv中的椭圆拟合
• Python and C implementations of an ellipse fitting algorithm in double and fixed point precision.
• C++ code for circle fitting algorithms
• OriginDelight/EllipseFit (c++源码)

论文作者，著名学者主页：

• Bob Fisher： http://homepages.inf.ed.ac.uk/rbf/
• Fitzgibbon, Pilu, Fisher: Direct Least Squares Fitting of Ellipses：http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/FITZGIBBON/ELLIPSE/