問(wèn)題

給定離散點(diǎn)集$X_i=(x_i,y_i),i=1,2,...N$，我們希望找到誤差最小的橢圓去擬合這些離散點(diǎn)。

方法

由于橢圓的形式可以給定， 自然我們將使用最小二乘法來(lái)求解橢圓。主要依據(jù)論文《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》。

橢圓擬合的難點(diǎn)

通常我們使用最小二乘法求解如下的最優(yōu)化問(wèn)題：

$$Min\sum _{i=1}^{N}f(x_i,E)$$
這里$f(x_i,E)$表示點(diǎn)$x_i$到E（指待擬合的橢圓）的最小距離。一般稱為幾何距離。但是我們很難直接給定幾何距離的解析表達(dá)，因此很難求出。因此我們退而求其次，我們采用：為橢圓寫下隱式方程，然后將點(diǎn)的坐標(biāo)帶入此隱式方程就得到了點(diǎn)到橢圓的距離。這種方法對(duì)于直線擬合、圓擬合，返回的就是到其的真實(shí)距離。而對(duì)于橢圓擬合，它返回的值是與距離有類似的屬性，但不是一個(gè)真正的距離值。因此這個(gè)距離被稱為代數(shù)距離。
橢圓可以用下面的隱式方程表達(dá)：

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

與直線類似（直線表達(dá)為：$a_{1}x+a_{2}y+a_3=0$）,上式的一組參數(shù)也是齊次量，即只能定義到一個(gè)比例因子。而且表達(dá)式也能表達(dá)雙曲線和拋物線。而橢圓
通常要求:$a_2^2?4a_1{a}_4<0$。最好通過(guò)令$a_2^2?4a_1{a}_4=?1$,即可同時(shí)解決這兩個(gè)問(wèn)題。

請(qǐng)參考維基百科橢圓定義和wolfram MathWorld 橢圓定義

一般的求解過(guò)程

對(duì)上面的橢圓方程改寫為：

$$f(a,(x,y))=D?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個(gè)離散點(diǎn)$X_i=(x_i,y_i),i=1...N$,通過(guò)最小化代數(shù)距離：

$$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矩陣。
另外，我們可以把各個(gè)$D_i$合并起來(lái)，則有
$$\hat{D}=?\begin{array}{c}{D}_1\\ D_2\\ ?\\ D_N\end{array}?$$
因此：
$$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$。
因此綜合上來(lái)即為，求解如下的最優(yōu)化問(wèn)題：

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

公式和拉格朗日乘子法

通過(guò)使用拉格朗日乘子法 ，我們可以得到

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

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

上式是經(jīng)典的廣義特征值問(wèn)題。我們可以直接求解，其中$\lambda$為廣義特征值，而a為廣義特征向量
論文 中指出有且僅有一個(gè)負(fù)的特征值，且其對(duì)應(yīng)的特征向量即為我們需要的橢圓方程的系數(shù)，即這里的a.

一般的圓錐方程到標(biāo)準(zhǔn)橢圓方程的轉(zhuǎn)化。

當(dāng)我們求解得到了圓錐曲線的系數(shù)，即a,我們一般需要獲得橢圓的標(biāo)準(zhǔn)方程，也就是獲得橢圓的如下參數(shù)：
中心$(c_x,c_y)$、長(zhǎng)短半軸$r_1,r_2$,旋轉(zhuǎn)角度$\theta$.

這里我們可以給出結(jié)果，也可以自己結(jié)合橢圓的標(biāo)準(zhǔn)形式與圓錐曲線的方程去推導(dǎo).

參考文獻(xiàn)：http://mathworld.wolfram.com/Ellipse.html.
##編程實(shí)現(xiàn)
下面描述matlab與c++實(shí)現(xiàn)。

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++

我們參考了網(wǎng)上的版本，并修復(fù)了他的問(wèn)題，最終也做出了和matlab一樣的效果。其中最關(guān)鍵的是廣義特征值的求解。我們使用了Eigen與clapack庫(kù)。其中Eigen易于表達(dá)矩陣，和matlab用法類似，是個(gè)強(qiáng)大的C++線性代數(shù)庫(kù)。而CLAPACK是線性代數(shù)包Lapack面向C/c++的接口。里面包含了很豐富的線性代數(shù)算法，包括廣義特征值求解接口，而且速度很快。我們希望將二者結(jié)合起來(lái)使用。

Eigen的安裝

Eigen直接以源代碼的方式提供給用戶，因此我們從官網(wǎng)上下載下后，直接在工程中包含其頭文件路徑即可。具體可參考：http://blog.csdn.net/abcjennifer/article/details/7781936

clapack的安裝

請(qǐng)查看官網(wǎng)，里面包含了詳細(xì)的使用與安裝步驟。

也可以使用我們已經(jīng)編譯了的vc2010和vc2013的庫(kù)，可以點(diǎn)擊下載。

盡管clapack面向c語(yǔ)言，因此需要我們?cè)诎^文件的時(shí)候，記得加上extern “C”.但是最新的版本（比如CLAPACK 3.2.1）已經(jīng)為我們?cè)陬^文件中加上了這些限制符，因此最新的版本可以兼容c和c++，所以直接在項(xiàng)目包含頭文件即可。

比如像下面一樣：

//Eigen #include <Eigen/Dense> #include <Eigen/Core> #include <iostream>//clapack,必須放在Eigen后面#include <f2c.h> #include <clapack.h>

而且應(yīng)該注意Eigen與CLAPACK混合使用的時(shí)候，CLAPACK的頭文件要加在Eigen的后面。否則會(huì)出錯(cuò)。

代碼

請(qǐng)查看Github主頁(yè)：https://github.com/xiamenwcy/EllipseFitting

實(shí)驗(yàn)結(jié)果

參考文獻(xiàn)

Eigen、LAPACK

Interfacing Eigen with LAPACK

Using BLAS/LAPACK from Eigen

CLAPACK for Windows

如何在VC中調(diào)用CLAPACK

使用Lapack庫(kù)的一些重要規(guī)則

在VS中用CLAPACK解決廣義特征值問(wèn)題

LAPACK 3.7.0

C＋＋編程：C++中同時(shí)使用Eigen和CLAPACK

CLAPACK在vc++6.0中成功調(diào)用

CLAPACK動(dòng)態(tài)調(diào)用

Leading dimension

leading dimension 2

橢圓擬合

• Fitting an Ellipse to a Set of Data Points(python實(shí)現(xiàn))
• 二次曲線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++源碼)

論文作者，著名學(xué)者主頁(yè)：

• 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/

總結(jié)

以上是生活随笔為你收集整理的基于代数距离的椭圆拟合的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。