基于代数距离的椭圆拟合

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

问题

给定离散点集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

椭圆拟合

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

查看评论

数学中几种常用的距离

数学中有很多不同种类的距离,常用于几何、高等代数等数学研究。多种多样的距离在数学建模、计算机学习中有着不小的应用。比如,A*搜索时的评估函数。比如,在机器学习中,做分类时常常需要估算不同样本之间的相似...
  • u013007900
  • u013007900
  • 2016年03月24日 11:01
  • 3810

机器视觉halcon软件:获取圆形的半径和圆度

1、要测量一个圆的直径,可以先通过阈值筛选出所需部分:threshold(); 2、然后将选中部分组合connection (),选择其中类圆度高的部分(根据选中部分与圆的相似度筛选出圆型部分)...
  • LeasonQ
  • LeasonQ
  • 2017年05月31日 18:17
  • 3091

Halcon学习备忘六(拟合圆的一种方法)

对于圆形区域,我们想拟合出它的圆形
  • haohaoxuexi320
  • haohaoxuexi320
  • 2014年08月11日 17:15
  • 8799

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

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

Halcon12.0 椭圆拟合

  • 2017年04月25日 16:40
  • 8KB
  • 下载

【Halcon】轮廓处理

 Coutours 属于XLD结构 a) 图  边缘覆盖在原始图像上 b)缩放 a)中的 矩形区域,凸显控制点 c)轮廓可以分割成线。圆。椭圆等,可以得出它们的 角度 中心 半径等 ...
  • huixingshao
  • huixingshao
  • 2015年05月08日 09:17
  • 23937

点到椭圆的最短距离

Minimum Distance to an EllipseDate: 06/10/99 at 11:18:23 From: Christian Gauthier Subject: Geometr...
  • stereohomology
  • stereohomology
  • 2015年09月02日 10:19
  • 1861

OpenCV 基于轮廓提取的二值图像分析与连通区域标记算法

简单介绍了OpenCV 提供的可用于二值图像分析的算法接口,并在指出现有算法接口不足的基础上,给出了一个带有标记结果输出的 连通区域标记算法接口。...
  • kewei9
  • kewei9
  • 2017年07月07日 16:28
  • 3459

椭圆曲线算法:入门(1)

https://www.jianshu.com/p/2e6031ac3d50?from=groupmessage 很多人都听说过加密算法,包括ECC、ECDH或者ECDSA。ECC是Elliptic...
  • PZ0605
  • PZ0605
  • 2018年01月11日 09:46
  • 516

椭圆曲线多倍点运算的实现:

 椭圆曲线多倍点运算的实现:有多种方法,下面只说二进制展开法 输入:点P,l比特的整数k[kl-1,…k0] 输出:Q=[k]P 1:置...
  • qq_24375053
  • qq_24375053
  • 2017年07月14日 15:29
  • 672
    个人资料
    持之以恒
    等级:
    访问量: 34万+
    积分: 4718
    排名: 7692
    个人网站
    最新评论