网站加入搜索引擎怎么做,网站备案怎么找人备,手机开发游戏,wordpress弹出小提示问题
给定离散点集Xi(xi,yi),i1,2,...NX_i(x_i,y_i) ,i1,2,...NXi(xi,yi),i1,2,...N#xff0c;我们希望找到误差最小的椭圆去拟合这些离散点。
方法
由于椭圆的形式可以给定#xff0c; 自然我们将使用最小二乘法来求解椭圆。主要依据论文《Direct least squares f…问题
给定离散点集Xi(xi,yi),i1,2,...NX_i(x_i,y_i) ,i1,2,...NXi(xi,yi),i1,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∑i1Nf(xi,E)Min \sum_{i1}^N f(x_i,E)Mini1∑Nf(xi,E) 这里f(xi,E)f(x_i,E)f(xi,E)表示点xix_ixi到E指待拟合的椭圆的最小距离。一般称为几何距离。但是我们很难直接给定几何距离的解析表达因此很难求出。因此我们退而求其次我们采用为椭圆写下隐式方程然后将点的坐标带入此隐式方程就得到了点到椭圆的距离。这种方法对于直线拟合、圆拟合返回的就是到其的真实距离。而对于椭圆拟合它返回的值是与距离有类似的属性但不是一个真正的距离值。因此这个距离被称为代数距离。 椭圆可以用下面的隐式方程表达
a1x2a2xya3y2a4xa5ya60a_1x^2a_2xya_3y^2a_4xa_5ya_60a1x2a2xya3y2a4xa5ya60
与直线类似直线表达为a1xa2ya30a_1xa_2ya_30a1xa2ya30,上式的一组参数也是齐次量即只能定义到一个比例因子。而且表达式也能表达双曲线和抛物线。而椭圆 通常要求:a22−4a1a4lt;0a_2^2-4a_1a_4lt;0a22−4a1a40。最好通过令a22−4a1a4−1a_2^2-4a_1a_4-1a22−4a1a4−1,即可同时解决这两个问题。
请参考维基百科椭圆定义和wolfram MathWorld 椭圆定义
一般的求解过程
对上面的椭圆方程改写为
f(a,(x,y))D⋅a0f(a,(x,y))D\cdot a0f(a,(x,y))D⋅a0
这里D(x2,xy,y2,x,y,1)D(x^2, xy, y^2, x, y, 1)D(x2,xy,y2,x,y,1),a(a1,a2,a3,a4,a5,a6)a(a_1, a_2, a_3, a_4, a_5, a_6)a(a1,a2,a3,a4,a5,a6). 于是给定N个离散点Xi(xi,yi),i1...NX_i(x_i,y_i),i1...NXi(xi,yi),i1...N,通过最小化代数距离
MinΔ(a,X)∑i1N(f(a,Xi))2Min \Delta(a,X)\sum_{i1}^{N}(f(a,X_i))^2MinΔ(a,X)i1∑N(f(a,Xi))2
重写上式即有
MinΔ(a,X)∑i1NaTDiTDiaaTSaMin \Delta(a,X)\sum_{i1}^{N}a^TD_i^TD_iaa^TSaMinΔ(a,X)i1∑NaTDiTDiaaTSa
这里S∑DiTDiS\sum D_i^T D_iS∑DiTDi为6x6矩阵。 另外我们可以把各个DiD_iDi合并起来则有 D^(D1D2⋮DN)\hat{D}\begin{pmatrix} D_1\\ D_2\\ \vdots \\ D_N \end{pmatrix}D^⎝⎜⎜⎜⎛D1D2⋮DN⎠⎟⎟⎟⎞ 因此 SD^TD^S\hat{D}^T\hat{D}SD^TD^
此外我们还有约束条件 a22−4a1a4lt;0a_2^2-4a_1a_4lt;0a22−4a1a40可以改写为
aTCa−1a^TCa-1aTCa−1
其中C∈R6×6C\in R^{6\times 6}C∈R6×6,且大部分为0只有C1,3C3,1−2,C2,21C_{1,3}C_{3,1}-2,C_{2,2}1C1,3C3,1−2,C2,21。 因此综合上来即为求解如下的最优化问题
MinaTSaMin~~~ a^TSa Min aTSa s.t.aTCa−1s.t. ~~~a^TCa-1s.t. aTCa−1
公式和拉格朗日乘子法
通过使用拉格朗日乘子法 我们可以得到
L(a)aTSa−λ(aTCa1)L(a)a^TSa-\lambda (a^TCa1)L(a)aTSa−λ(aTCa1)
然后求解∂aL(a)0\partial_aL(a)0∂aL(a)0可以求解得到 SaλCaSa\lambda CaSaλCa
上式是经典的广义特征值问题。我们可以直接求解其中λ\lambdaλ为广义特征值而a为广义特征向量 论文 中指出有且仅有一个负的特征值且其对应的特征向量即为我们需要的椭圆方程的系数即这里的a.
一般的圆锥方程到标准椭圆方程的转化。
当我们求解得到了圆锥曲线的系数即a,我们一般需要获得椭圆的标准方程也就是获得椭圆的如下参数 中心(cx,cy)(c_x,c_y)(cx,cy)、长短半轴r1,r2r_1,r_2r1,r2,旋转角度θ\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 datamx 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 vectorsx x(:);y y(:);% Build design matrixD [ x.*x x.*y y.*y x y ones(size(x)) ];% Build scatter matrixS D*D;% Build 6x6 constraint matrixC(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 eigenvalueA gevec(:,NegC);% unnormalizea [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 LAPACKUsing BLAS/LAPACK from EigenCLAPACK for Windows 如何在VC中调用CLAPACK使用Lapack库的一些重要规则在VS中用CLAPACK解决广义特征值问题LAPACK 3.7.0C编程C中同时使用Eigen和CLAPACKCLAPACK在vc6.0中成功调用CLAPACK动态调用Leading dimensionleading dimension 2
椭圆拟合
Fitting an Ellipse to a Set of Data Points(python实现)二次曲线Quadratic Curvefitellipse.cpp(opencv)Creating Bounding rotated boxes and ellipses for contoursOpenCV’s fitEllipse() sometimes returns completely wrong ellipsesFitting an Ellipse to a Set of Data PointsDirect Least Square Fitting of Ellipses (论文)best-fitting-line-circle-and-ellipseopencv中的椭圆拟合Python and C implementations of an ellipse fitting algorithm in double and fixed point precision.C code for circle fitting algorithmsOriginDelight/EllipseFit (c源码)
论文作者著名学者主页
Bob Fisher http://homepages.inf.ed.ac.uk/rbf/Fitzgibbon, Pilu, Fisher: Direct Least Squares Fitting of Ellipseshttp://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/FITZGIBBON/ELLIPSE/
