关闭

基于代数距离的椭圆拟合

标签: 椭圆拟合
1412人阅读 评论(2) 收藏 举报
分类:

问题

给定离散点集Xi=(xi,yi),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》。

椭圆拟合的难点

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

Mini=1Nf(xi,E)

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

a1x2+a2xy+a3y2+a4x+a5y+a6=0

与直线类似(直线表达为:a1x+a2y+a3=0),上式的一组参数也是齐次量,即只能定义到一个比例因子。而且表达式也能表达双曲线和抛物线。而椭圆
通常要求:a224a1a4<0。最好通过令a224a1a4=1,即可同时解决这两个问题。

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

一般的求解过程

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

f(a,(x,y))=Da=0

这里D=(x2,xy,y2,x,y,1),a=(a1,a2,a3,a4,a5,a6).
于是给定N个离散点Xi=(xi,yi),i=1...N,通过最小化代数距离:

MinΔ(a,X)=i=1N(f(a,Xi))2

重写上式,即有:

MinΔ(a,X)=i=1NaTDTiDia=aTSa

这里S=DTiDi为6x6矩阵。
另外,我们可以把各个Di合并起来,则有

D^=D1D2DN

因此:
S=D^TD^

此外,我们还有约束条件:
a224a1a4<0,可以改写为:

aTCa=1

其中CR6×6,且大部分为0,只有C1,3=C3,1=2,C2,2=1
因此综合上来即为,求解如下的最优化问题:

Min   aTSa

s.t.   aTCa=1

公式和拉格朗日乘子法

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

L(a)=aTSaλ(aTCa1)

然后求解aL(a)=0可以求解得到:

Sa=λCa

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

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

当我们求解得到了圆锥曲线的系数,即a,我们一般需要获得椭圆的标准方程,也就是获得椭圆的如下参数:
中心(cx,cy)、长短半轴r1,r2,旋转角度θ.

这里我们可以给出结果,也可以自己结合椭圆的标准形式与圆锥曲线的方程去推导.
image
参考文献: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一样的效果。其中最关键的是广义特征值的求解。我们使用了Eigenclapack库。其中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

实验结果

image

参考文献

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

椭圆拟合

论文作者,著名学者主页:

0
0
查看评论
发表评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场

基于几何距离的椭圆拟合

问题给定离散点集Xi=(xi,yi)X_i=(x_i,y_i),我们希望找到最好的椭圆去拟合这些离散点。方法通常我们使用最小二乘法求解如下的最优化问题:Min∑i=1Nf(xi,E)2 Min \su...
  • xiamentingtao
  • xiamentingtao
  • 2017-02-08 23:24
  • 1597

图像处理中的椭圆拟合(一)

图像处理中的椭圆检测用处还是挺多的,找到这里来的同学大多是想用椭圆检测来解决某些实际问题吧,所以我就不做介绍,直奔主题。我研究这块也有一段时间了,也查找了挺多资料,貌似通用的椭圆算法还没有,不排除我孤...
  • zhiyuanzhe007
  • zhiyuanzhe007
  • 2014-03-30 21:18
  • 2167

最小二乘法椭圆拟合

感谢原作者的分享 直接上代码: .h文件为: //在这里,我实现了两种算法,一种是 //http://wenku.baidu.com/link?url=7kIrC8LoOMCtlmAH8yqkp...
  • xiao_lxl
  • xiao_lxl
  • 2015-07-02 14:13
  • 3129

基于几何距离的椭圆拟合

问题给定离散点集Xi=(xi,yi)X_i=(x_i,y_i),我们希望找到最好的椭圆去拟合这些离散点。方法通常我们使用最小二乘法求解如下的最优化问题:Min∑i=1Nf(xi,E)2 Min \su...
  • xiamentingtao
  • xiamentingtao
  • 2017-02-08 23:24
  • 1597

最小二乘拟合圆及椭圆

opencv-最小二乘法拟合圆及椭圆
  • qq_32171677
  • qq_32171677
  • 2017-11-30 15:14
  • 167

opencv中的椭圆拟合

参考:http://blog.sina.com.cn/s/blog_662c785901011i7z.html
  • suky520
  • suky520
  • 2014-01-21 12:09
  • 8149

OpencV使用fitEllipse拟合椭圆后,获取椭圆参数

OpencV使用fitEllipse拟合椭圆后,获取椭圆参数
  • sinat_31425585
  • sinat_31425585
  • 2017-07-14 12:04
  • 1455

OpenCV:椭圆拟合

#include #include int slider_pos = 70; IplImage *image02 = 0, *image03 = 0, *image04 = 0; void p...
  • u010002704
  • u010002704
  • 2015-02-11 15:33
  • 2703

椭圆拟合算法总结

椭圆拟合算法 椭圆拟合在医疗图像、工程设计、交通识别、生物识别、人工智能、林火蔓延等领域具有广泛应用,本文总结前人研究综合出基于MATLAB的三种椭圆拟合的算法:直接计算方法,二次多项式拟合方法,标...
  • shumu_xzw
  • shumu_xzw
  • 2015-07-31 15:26
  • 5577

最小二乘法用于直线,多项式,圆,椭圆的拟合及程序实现

参考: http://blog.csdn.net/jairuschan/article/details/7517773/
  • piaoxuezhong
  • piaoxuezhong
  • 2017-03-31 13:05
  • 1610
    个人资料
    • 访问:284118次
    • 积分:4121
    • 等级:
    • 排名:第8689名
    • 原创:117篇
    • 转载:82篇
    • 译文:8篇
    • 评论:106条
    个人网站
    最新评论