基于代数距离的椭圆拟合
来源:互联网 发布:centos制作u盘启动 编辑:程序博客网 时间:2024/05/29 07:32
问题
给定离散点集
方法
由于椭圆的形式可以给定, 自然我们将使用最小二乘法来求解椭圆。主要依据论文《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》。
椭圆拟合的难点
通常我们使用最小二乘法求解如下的最优化问题:
这里
椭圆可以用下面的隐式方程表达:
与直线类似(直线表达为:
通常要求:
请参考维基百科椭圆定义和wolfram MathWorld 椭圆定义
一般的求解过程
对上面的椭圆方程改写为:
这里
于是给定N个离散点
重写上式,即有:
这里
另外,我们可以把各个
因此:
此外,我们还有约束条件:
其中
因此综合上来即为,求解如下的最优化问题:
公式和拉格朗日乘子法
通过使用拉格朗日乘子法 ,我们可以得到
然后求解
上式是经典的广义特征值问题。我们可以直接求解,其中
论文 中指出有且仅有一个负的特征值,且其对应的特征向量即为我们需要的椭圆方程的系数,即这里的a.
一般的圆锥方程到标准椭圆方程的转化。
当我们求解得到了圆锥曲线的系数,即a,我们一般需要获得椭圆的标准方程,也就是获得椭圆的如下参数:
中心
这里我们可以给出结果,也可以自己结合椭圆的标准形式与圆锥曲线的方程去推导.
参考文献: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
- Interfacing Eigen with LAPACK
- Using BLAS/LAPACK from Eigen
- CLAPACK for Windows
- 如何在VC中调用CLAPACK
- 使用Lapack库的一些重要规则
- 在VS中用CLAPACK解决广义特征值问题
- LAPACK 3.7.0
- C++编程:C++中同时使用Eigen和CLAPACK
- CLAPACK在vc++6.0中成功调用
- CLAPACK动态调用
- Leading dimension
- 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/
- 基于代数距离的椭圆拟合
- 基于几何距离的椭圆拟合
- 最小二乘拟合&基于RANSAC的直线拟合&椭圆拟合
- 椭圆与圆的拟合
- 椭圆拟合
- 椭圆拟合
- 椭圆中心到椭圆切线的距离
- opencv 椭圆拟合
- 圆及椭圆拟合
- opencv中的椭圆拟合
- matlab 椭圆方程拟合
- 椭圆拟合fit_ellipse.m
- opencv 椭圆拟合
- 【OpenCV】椭圆拟合
- OpenCV:椭圆拟合
- 最小二乘法椭圆拟合
- 最小二乘法椭圆拟合
- 椭圆拟合算法总结
- Performing summary statistics and plots —— Python Data Science Cookbook
- 深入理解JVM(二)——揭开HotSpot对象创建的奥秘
- UOJ 73 [WC2015]未来程序
- Linux下" >/dev/null 2>&1 "相关知识说明
- Android获取联系人
- 基于代数距离的椭圆拟合
- Spring_22_基于配置文件的方式来配置 AOP
- U盘制作Apple系统安装盘
- android-BroadcastReceiver 发送有序广播
- 找老乡 (sdut oj)
- 数据结构实验之二叉树的建立与遍历
- CF699A
- 关于Jadepool3.0的使用
- post请求https接口