Matlab在離散點擬合橢圓及極值距離計算中的應(yīng)用

吳美容，王建國

(南京太亞科技有限責(zé)任公司，江蘇 南京 210061)

Matlab在離散點擬合橢圓及極值距離計算中的應(yīng)用

吳美容，王建國

(南京太亞科技有限責(zé)任公司，江蘇 南京 210061)

目前關(guān)于根據(jù)離散點建立橢圓的一般二次曲線方程，并由此推導(dǎo)橢圓的標(biāo)準(zhǔn)方程以及計算離散點距離橢圓極值方面的文獻(xiàn)較少，文章重點介紹了利用Matlab軟件進(jìn)行公式推導(dǎo)的詳細(xì)過程以及計算極值的方法。極值計算方法介紹了幾何作圖法、迭代計算法、窮舉法、人工智能算法。實例驗證表明了思路與算法的正確性和實用性，可以為解決其他二次曲線和二次曲面的問題提供借鑒。

橢圓；擬合；距離極值；算法；人工智能；最優(yōu)化

0 引言

橢圓、雙曲線、拋物線等都是平面二次曲線，本文針對橢圓進(jìn)行討論。張元元將橢圓擬合應(yīng)用在工程隧道形體檢測[1]，閆蓓將橢圓擬合應(yīng)用在醫(yī)學(xué)瞳孔中心定位系統(tǒng)[2]，應(yīng)該說現(xiàn)實世界中橢圓形狀的物體較多，經(jīng)常會碰到需要確定這些橢圓實體的幾何參數(shù)(中心、長短半軸、旋轉(zhuǎn)角度)，通過采集橢圓上一系列離散點進(jìn)行橢圓擬合是較為常用的處理方法。目前較多的文獻(xiàn)都是直接給出公式，側(cè)重于橢圓擬合的應(yīng)用，卻沒有給出數(shù)據(jù)處理的整體思路以及公式推導(dǎo)的方法，讀者不清楚公式是如何得出的，有的公式甚至還是錯誤的，直接套用會發(fā)現(xiàn)結(jié)果不正確，如文獻(xiàn)[2]公式(4)、(5)。本文著重介紹Matlab在橢圓擬合數(shù)據(jù)處理中的具體應(yīng)用。此外，目前討論離散點距擬合橢圓極值距離的文獻(xiàn)較少，本文也在這方面進(jìn)行了探討。

1 橢圓方程的建立

橢圓的一般二次曲線方程可表示為：

F(X，Y)=AX2+BXY+CY2+DX+EY+F=0

(1)

式(1)有6個未知參數(shù)，其中，F(xiàn)為常數(shù)項，因此，至少需要5個以上的點來確定橢圓的二次曲線方程。為了數(shù)據(jù)處理的方便，一般需要對原始數(shù)據(jù)X、Y進(jìn)行正規(guī)化處理，代碼描述：

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；

正規(guī)化數(shù)據(jù)對應(yīng)的二次曲線方程記為：

f(x，y)=ax2+bxy+cy2+dx+ey+f=0

(2)

對于計算式(2)中的6個參數(shù)，可以調(diào)用Matlab自帶的regress()，該函數(shù)為基于最小二乘原理多元線性回歸，構(gòu)造矩陣：

輸入命令：regress(Y，X)即可算出系數(shù)。該步驟也可以輸入如下命令：inv(X′*X)*X′*Y。實際計算時可取f=1。通過以上步驟，6個參數(shù)abcdef都已經(jīng)計算出來了，為已知值?，F(xiàn)在要反歸一化，通過abcdef反算ABCDEF。推導(dǎo)過程借助Matlab符號運算，代碼描述：

syms x y X Y mx my sx sy a b c d e f；

x =(X-mx)/sx；

y =(Y-my)/sy；

f=a*x^2+b*x*y+c*y^2+d*x+e*y+f；

pretty(f)

結(jié)果：

分別將X^2，XY，Y^2，X，Y系數(shù)提取出來，可以輔助用collect()命令簡化提取工作。

通過以上步驟，ABCDEF已經(jīng)計算出來了，為已知值。截至目前，橢圓的二次曲線方程已知。下面推導(dǎo)如何轉(zhuǎn)換為幾何標(biāo)準(zhǔn)形式方程，進(jìn)而計算出5個參數(shù)。

首先假設(shè)式(1)中的各系數(shù)已知?，F(xiàn)在把坐標(biāo)軸逆時針轉(zhuǎn)動θ角度，其中(x′，y′)表示點(x，y)在坐標(biāo)軸變動后的新坐標(biāo)：

(3)

把式(3)代入式(1)，設(shè)新方程為：

a′x′2+b′x′y′+c′y′2+d′x′+e′y′+f′=0

(4)

容易求出：

a′=Acos2θ+Bcosθsinθ+Csin2θ

b′=(C-A)sin2θ+Bcos2θ

c′=Asin2θ+Bcosθsinθ+Ccos2θ

d′=Dcosθ+Esinθ

e′=-Dsinθ+Ecosθ

f′=F

為了使坐標(biāo)軸變換后，方程不再出現(xiàn)x′y′一項，即b′=0，則有：

(5)

為了讓x′y′項的系數(shù)為0，需令

(6)

(7)

則有：

(8)

此時的U和V即橢圓的兩半軸長，再利用式(3)將目前的橢圓中心點坐標(biāo)(x0′，y0′)轉(zhuǎn)變?yōu)樽鴺?biāo)軸變換之前的(x0，y0)。以上符號公式的推導(dǎo)仍然借助于Matlab符號運算。

于是，從橢圓一般方程的系數(shù)得到了橢圓的兩半軸長、中心坐標(biāo)以及U半軸關(guān)于X坐標(biāo)軸的角度(逆時針為正向)。

2 離散點距擬合橢圓極值距離計算

方法1(幾何作圖法)：

離散點到曲線最小距離，采用文獻(xiàn)[3]里面的第1種思路也是可以的，但是手工計算存在效率較低的問題，可作為一種輔助計算手段。

方法2(解非線性方程組法)：

文獻(xiàn)[3]里面的第2種思路也可以，k1×k2=-1，根據(jù)該思路可建立如下2個方程：

(9)

Ax02+Bx0y0+Cy02+Dx0+Ey0+F=0

(10)

具體如下：

將式(9)、式(10)聯(lián)立計算離散點(px，py)對應(yīng)于擬合橢圓上的最近點(x0，y0)，反算得到距離極值。求解上面的方程組得出關(guān)于x0、y0的解析式較麻煩，考慮迭代算法。直接利用fsolve()函數(shù)，以點(112.6614，-165.52)為例。首先建立M文件：

function f=myfun(x)

A=4.89939E-05；B=5.62288E-08；

C=0.000115288；D=-0.008437616；

E=0.016865307；F=-0.031413295；

px=112.6614；py=-165.52；

f1=(-(2*A*x(1)+B*x(2)+D)/(B*x(1)+2*C*x(2)+E)*(py-x(2))/(px-x(1)))+1；

f2=(A*x(1)^2+B*x(1)*x(2)+C*x(2)^2+D*x(1)+E*x(2)+F)；

f=[f1 f2]；

end

然后輸入如下代碼：

fsolve(@myfun，[112 -165]，optimset(‘MaxFunEvals’，100，‘MaxIter’，100))%最小距離對應(yīng)的橢圓點

結(jié)果如下：

ans =112.6287 -165.2517

輸入如下代碼：

fsolve(@myfun，[14 -26]，optimset(‘MaxFunEvals’，100，‘MaxIter’，100))%最大距離對應(yīng)的橢圓點

結(jié)果如下：

ans =-35.8312 -23.6046

采用該方法能很好計算出正確的結(jié)果，坐標(biāo)反算最小最大距離分別為：0.270 3，205.402 1。針對該函數(shù)計算時較依賴初值，容易陷入局部極小，采取如下對策：計算最短距離時，初值可直接取離散點數(shù)值的整數(shù)部分或所在象限中以橢圓長短半軸為兩邊所構(gòu)成的長方形區(qū)域(近似)中心點的數(shù)值；計算最大距離時，初值可取對角象限中以橢圓長短半軸為兩邊所構(gòu)成的長方形區(qū)域(近似)中心點的數(shù)值。

方法3(傳統(tǒng)最優(yōu)化算法)：

>>syms t；px′=26.5498；py′=-92.3433；

x=143.7082651857737*cos(t)；y=93.6831595721793*sin(t)；

f=sqrt((x-px′).^2+(y-py′).^2)；

>> [x，minf]=minNewton(f，0)

計算結(jié)果：

x =2.5848

minf =205.4007

>> [x，minf]=minPWX(f，0，4)

計算結(jié)果：

x =4.8980

minf =0.2707

以上分別用基本牛頓法和拋物線法計算了極大值和極小值。該方法較實用，如調(diào)用Matlab自帶的函數(shù)fminsearch()、fminunc()也可以。

方法4(人工智能算法)：

群智能算法SIA，本文以粒子群算法和魚群算法為例進(jìn)行介紹，經(jīng)筆者測試其他智能算法，如遺傳算法ga()、模擬退火算法saa()等也都能在很短的時間內(nèi)獲得準(zhǔn)確的結(jié)果。

粒子群算法：

首先建立目標(biāo)函數(shù)文件myfunpso()：

function f=myfunpso(t)

px1=26.5498；py1=-92.3433；

x=143.7082651857737*cos(t)；

y=93.6831595721793*sin(t)；

f=sqrt((x-px1).^2+(y-py1).^2)；

end

命令窗口輸入：

>> [xm，fv]=PSO(@myfunpso，40，2，2，0.5，100，1)

結(jié)果如下：

xm =-1.385213897691377 fv =0.270784374867295

計算最大值，將目標(biāo)函數(shù)文件改為：f=-sqrt((x-px1).^2+(y-py1).^2)；

輸入命令后，結(jié)果如下：

xm =2.584791470317521 fv =-2.054007419382480e+002

魚群算法迭代過程如圖1所示。

圖1 魚群算法迭代過程Fig.1 Iterative process of fish-swarm algorithm

對比上面方法3用傳統(tǒng)的最優(yōu)化算法，結(jié)果是一致的。群智能算法屬于隨機(jī)算法，不要求函數(shù)連續(xù)、可導(dǎo)，應(yīng)用時需要設(shè)置恰當(dāng)?shù)膮?shù)，多次運行程序以提高可靠性。

方法5(遍歷法)：

本文再介紹一種經(jīng)典的遍歷算法來計算極值。思路：第一，將離散點坐標(biāo)轉(zhuǎn)換到標(biāo)準(zhǔn)橢圓坐標(biāo)系；第二，判斷離散點位于哪個象限；第三，離散點與屬于該象限橢圓上的點按步長計算距離并比較取最小值，此外，一、二象限，離散點y值大于橢圓點y值，距離取正，否則取負(fù)；三、四象限，離散點y值小于橢圓點y值，距離取正，否則取負(fù)；第四，離散點與所屬象限的對角象限上的橢圓點按步長計算距離并比較取最大值。算法流程圖如圖2所示。

圖2 算法流程圖Fig.2 Algorithm flow chart

程序采用的思路，根據(jù)點位落在的象限區(qū)間，縮小搜索范圍，不需要每個點都要全區(qū)域搜索，提高了運行的效率。

3 具體應(yīng)用

測試用例采用文獻(xiàn)[4]，具體數(shù)據(jù)如表1所示，程序截圖如圖3、圖4、圖5所示。

表1 離散點數(shù)據(jù)

圖3 讀取數(shù)據(jù)文件Fig.3 To read the data files

圖4 計算結(jié)果Fig.4 Calculation results

圖5 離散點與擬合橢圓Fig.5 Discrete points and fitting ellipse

4 結(jié)論

通過應(yīng)用分析，筆者可得到如下結(jié)論：

1)程序只需讀取原始測量數(shù)據(jù)TXT文件即可自動輸出擬合橢圓及離散點分布圖，直觀明了；能夠自動輸出橢圓5個參數(shù)以及每個離散點對應(yīng)于擬合橢圓的最大距離與最短距離；文獻(xiàn)[4]橢圓中心坐標(biāo)值86.105 7錯誤，應(yīng)該改為86.150 7，最小值輸出結(jié)果與文獻(xiàn)[4]結(jié)果相差2 mm內(nèi)。程序結(jié)果正確，操作簡單，自動化程度高。

2)程序橢圓擬合算法基于最小二乘原理，最小二乘原理本身不具備較強(qiáng)的抗差能力，若加上穩(wěn)健估計或粗差剔除功能模塊，將更加提高程序的自動化程度。

3)無論平面坐標(biāo)轉(zhuǎn)換還是空間三維坐標(biāo)轉(zhuǎn)換都是工程測量和工業(yè)測量數(shù)據(jù)處理過程中最基礎(chǔ)也是必要的步驟。

4)文章介紹了5種計算離散點到擬合橢圓極值距離的數(shù)據(jù)處理方法，實踐中綜合運用以上5種方法可以取得很好的計算效果。計算極值問題，若采用人工智能算法，如：模擬退火算法等，具有很好全局最優(yōu)解的搜索能力，屬于隨機(jī)算法，對初值也不敏感，若設(shè)置好冷卻進(jìn)度表相關(guān)參數(shù)，計算效果會很好；若采用傳統(tǒng)算法，為避免陷入局部極小，采用作圖法作大致判斷也是很有必要。各種算法的優(yōu)勢互補(bǔ)，算法集成是很有意義的研究方向。例如，將模擬退火算法、遺傳算法、粒子群算法與文獻(xiàn)[5]的BP神經(jīng)網(wǎng)絡(luò)算法進(jìn)行融合，優(yōu)化神經(jīng)網(wǎng)絡(luò)連接權(quán)值和閾值。

5)程序的思路可以為點到一般二次曲線(如：雙曲線，拋物線等)及點到一般二次曲面的極值計算問題提供借鑒。

6)文獻(xiàn)[6]較詳細(xì)的介紹了溢流堰表面形體檢測工作中涉及的檢測點距離直線、多圓心圓弧、拋物線、緩和曲線等線元最小值的計算方法，實際上溢流堰圖紙也經(jīng)常設(shè)計有四分之一橢圓，本文針對點距離該線元極值的計算提出了比較實用的解決思路。

[1] 張元元，楊國東，王鳳艷.丹麥法穩(wěn)健估計在隧道橢圓擬合中的應(yīng)用[J].世界地質(zhì)，2012，31(1)：199-203.

[2] 閆蓓，王斌，李媛.基于最小二乘法的橢圓擬合改進(jìn)算法[J].北京航空航天大學(xué)學(xué)報，2008，34(3)：295-298.

[3] 王建國.迭代計算方法在測繪領(lǐng)域中的具體應(yīng)用[J].城市勘測，2011(2)：153-154

[4] 王解先，季凱敏.工業(yè)測量擬合[M].北京：測繪出版社，2008：57-61.

[5] 王建國，吳美容.運用BP人工神經(jīng)網(wǎng)絡(luò)設(shè)計變形預(yù)報模型[J].礦山測量，2010(4)：73-75.

[6] 吳美容，王建國.溢流堰表面形體檢測研究[J].現(xiàn)代測繪，2016，39(2)：54-56.

Application of Matlab in Fitting Ellipse and Calculating Extremum Distance of Discrete Points to Ellipse

WU Mei-rong，WANG Jian-guo

(NanjingPaciaTechnology&ScienceCo.,Ltd.，NanjingJiangsu210061，China)

According to the fact that there is no enough documents about establishing general quadratic curve equation based on discrete points，deducing standard geometric symbolic equation，and calculating the extremum distance with the use of mathematical computer language(Matlab).This article focuses on the detailed process of formula derivation and the method of calculating the extreme value with Matlab.The article presents the following methods：geometric drawing method，iterative method，exhaustive method，and AI algorithm.The practical example verifies the correctness and practicability of the method and the algorithm，which can be used for reference to solve the problems of other conic curves and quadric surfaces.

ellipse；fitting；distance extremum；algorithm；artificial intelligence；optimization

2016-07-05

TP 301.6

：A

：1007-9394(2016)04-0020-04

吳美容(1982～)，女，江蘇海安人，學(xué)士，工程師，現(xiàn)主要從事地理信息系統(tǒng)數(shù)據(jù)處理及項目管理方面的工作。