同心圓亞像素中心定位方法

張方程, 李曉天

(中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所， 吉林 長(zhǎng)春 130033)

同心圓亞像素中心定位方法

張方程, 李曉天

(中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所， 吉林 長(zhǎng)春 130033)

基于橢圓曲線擬合和曲率濾波,提出同心圓投影中心定位方法。采用單應(yīng)矩陣變換和坐標(biāo)平移法生成了大量的理想測(cè)試圖像，并與角點(diǎn)檢測(cè)法的測(cè)量精度進(jìn)行了比較，結(jié)果表明，文中方法的測(cè)量精度在有無(wú)噪聲情況下均高于角點(diǎn)檢測(cè)方法。當(dāng)信噪比介于35 dB和25 dB之間時(shí)，角點(diǎn)檢測(cè)法的測(cè)量誤差突然急劇增大，文中方法的測(cè)量誤差變化仍比較平緩，測(cè)量誤差小于0.15 pixel。

機(jī)器視覺(jué)； 圖像處理； 同心圓定標(biāo)； 橢圓擬合

0 引 言

平面同心圓的投影中心可用作對(duì)攝像機(jī)進(jìn)行內(nèi)外參數(shù)標(biāo)定，彌補(bǔ)角點(diǎn)檢測(cè)法和圓心法的定位不準(zhǔn)確問(wèn)題。一個(gè)圓的投影圓心與其在圖像中的橢圓中心不具備仿射不變對(duì)應(yīng)關(guān)系[1-3]，因此采用圖像橢圓中心來(lái)代替投影圓心的成像坐標(biāo)方法將帶來(lái)一定的定位誤差。根據(jù)小孔成像和仿射變換原理，平面內(nèi)同心圓在經(jīng)過(guò)攝像機(jī)拍攝后的數(shù)字圖像中表現(xiàn)形式為兩個(gè)中心分離的橢圓，當(dāng)且僅當(dāng)攝像機(jī)光軸與同心圓所在平面嚴(yán)格垂直時(shí)，同心圓所成的像為兩個(gè)中心重合的圓。2001年Kim和Kweon[1]指出，按照仿射變換原理，平面同心圓的投影中心與兩個(gè)成像橢圓中心共線，并提出采用交比不變的方法對(duì)同心圓的投影中心進(jìn)行定位。隨后，Xianghua Ying[4],Jun-Sik Kim[2],Ashutosh Morde[3],Jun Peng Xue[5]等學(xué)者對(duì)同心圓投影中心定位方法進(jìn)行了研究，但主要圍繞采用不同幾何點(diǎn)進(jìn)行仿射不變方法實(shí)現(xiàn)同心圓投影中心定位及其在攝像機(jī)定標(biāo)中的應(yīng)用進(jìn)行探討，沒(méi)有對(duì)同心圓的投影中心的提取方法及其提取精度進(jìn)行分析。

文中提出采用經(jīng)過(guò)曲率濾波的曲線擬合方法來(lái)實(shí)現(xiàn)同心圓投影中心的亞像素定位。首先提取兩橢圓邊界的亞像素定位坐標(biāo)，并采用曲率濾波方法剔除掉誤差較大的橢圓邊界點(diǎn)對(duì)兩個(gè)橢圓方程進(jìn)行最小二乘法擬合。然后根據(jù)兩個(gè)橢圓中心亞像素級(jí)坐標(biāo)以及兩個(gè)橢圓的方程，采用仿射變換方法對(duì)同心圓的投影中心進(jìn)行定位和逼近。最后在有無(wú)噪聲的情況下對(duì)文中方法和角點(diǎn)檢測(cè)法[6]的圖像點(diǎn)提取精度進(jìn)行對(duì)比仿真分析，并給出了實(shí)驗(yàn)測(cè)量結(jié)果。

1 數(shù)學(xué)建模

1.1圖像預(yù)處理

攝像機(jī)采集到的圖像存在著大量的噪聲，因此在進(jìn)行邊緣檢測(cè)前要先對(duì)圖像進(jìn)行平滑去噪處理。這里采用的是具有一定的邊緣保持功能的二維高斯濾波器

(1)

然后通過(guò)閾值化方法對(duì)投影同心圓質(zhì)心進(jìn)行粗定位,并提取出包含同心圓的感興趣區(qū)域(ROI)，后續(xù)的計(jì)算只對(duì)該區(qū)域所包含的256色灰度圖像進(jìn)行操作，如圖1所示。

1.2雙線性插值與橢圓邊緣的粗定位

根據(jù)方形孔徑原理，由于光學(xué)元器件的卷積以及衍射作用，圖像邊緣灰度值變換表征為高斯分布，高斯曲線的頂點(diǎn)對(duì)應(yīng)的位置即為真實(shí)邊緣點(diǎn)的位置。為提高計(jì)算速度，文中采用一維高斯邊緣檢測(cè)的方法對(duì)橢圓邊緣進(jìn)行粗定位。首先通過(guò)質(zhì)心(xg，yg)的θ角在0到 π之間的各條直線所對(duì)應(yīng)的像素坐標(biāo)進(jìn)行鄰近點(diǎn)雙線性插值，然后提取一維亞像素級(jí)高斯邊緣，從而獲得兩個(gè)橢圓的邊緣點(diǎn)圖像坐標(biāo)。

圖1 ROI圖像的像素級(jí)邊緣提取

1.2.1 雙線性插值

在進(jìn)行雙線性插值前，要先將直線L1離散化，此時(shí)要根據(jù)L1斜率的不同，而選擇不同的離散方法。當(dāng)|tanθ|≤1時(shí)，令x在直線L1被圖1所包含的范圍內(nèi)以1個(gè)像素等間隔進(jìn)行遞增變化，根據(jù)直線L1方程對(duì)相應(yīng)的y值進(jìn)行求解；當(dāng)|tanθ|>1時(shí)，令y在直線L1被圖1所包含的范圍內(nèi)以1個(gè)像素等間隔進(jìn)行遞增變化，根據(jù)直線L1方程對(duì)相應(yīng)的x值進(jìn)行求解。

雙線性插值方法認(rèn)為圖像的灰度在橫向和縱向方向上都是線性變化的，這與光學(xué)成像的漸變?cè)硐嘁恢耓6-8]。設(shè)直線L1上任一點(diǎn)P的坐標(biāo)為(xp,yp)，若xp和yp皆為整數(shù)，則P點(diǎn)的灰度值IP可直接從同心圓圖像中讀取，記為I(xp,yp)。若xp和yp不全為整數(shù)，則P點(diǎn)的雙線性插值灰度值IP為:

(2)

式中:a,b----分別與xp,yp最近整數(shù);

I(a,b),I(a+1,b),I(a,b+1)，I(a+1,b+1)----分別為對(duì)應(yīng)的整數(shù)點(diǎn)坐標(biāo)的灰度值。

1.2.2 一階高斯離散變換

設(shè)離散高斯變換的長(zhǎng)度f(wàn)為奇數(shù)，則離散的高斯函數(shù)一階導(dǎo)數(shù)的表達(dá)式為

(3)

沿著直線L1的方向?qū)D像進(jìn)行雙線性插值并獲得一維圖像數(shù)據(jù)后,采用離散的高斯函數(shù)一階導(dǎo)數(shù)對(duì)一維圖像灰度數(shù)據(jù)進(jìn)行卷積運(yùn)算，可獲得曲線如圖2所示。

圖2 ROI圖像某θ角方向的一維一階高斯變換曲線

由于圖2曲線峰值附近各點(diǎn)呈一維高斯分布，而對(duì)一維高斯函數(shù)的等式兩端取自然對(duì)數(shù)，獲得的新函數(shù)變換規(guī)律符合拋物線方程，因此,可根據(jù)曲線頂點(diǎn)附近3點(diǎn)的灰度差值(M軸的坐標(biāo)值)的自然對(duì)數(shù)求得拋物線的頂點(diǎn)，該點(diǎn)所對(duì)應(yīng)的圖像位置即為邊緣點(diǎn)亞像素級(jí)位置。

由圖2中的虛線圓所包含的頂峰上的頂點(diǎn)坐標(biāo)(d0,M0)以及其左右兩點(diǎn)坐標(biāo)(d0-1,ML)和(d0+1,MR)，對(duì)拋物線頂點(diǎn)坐標(biāo)進(jìn)行計(jì)算，求得邊緣點(diǎn)的坐標(biāo)為

(4)

將d0sub根據(jù)tan(θ)的取值不同，代入到直線L1的方程中，即可獲得亞像素級(jí)邊緣點(diǎn)的坐標(biāo)。

1.3曲率濾波

由于圖像中噪聲的存在，使得少部分邊緣點(diǎn)的計(jì)算位置與實(shí)際位置存在較大偏差。因此,在進(jìn)行橢圓方程曲線擬合前，需將這些計(jì)算誤差較大的邊緣點(diǎn)過(guò)濾掉，否則將影響橢圓方程的擬合精度。

設(shè)(x1,y1),(x2,y2),(x3,y3)為橢圓邊緣上的任意連續(xù)3點(diǎn)，則(x2,y2)點(diǎn)處的曲率kr計(jì)算公式為

(5)

其中

根據(jù)橢圓方程可知，橢圓的曲率在長(zhǎng)軸的端點(diǎn)處取最大值，在短軸端點(diǎn)處取最小值，曲率變化呈漸變趨勢(shì)，因此可根據(jù)式(5)將曲率值發(fā)生突變的邊緣點(diǎn)過(guò)濾掉。文中對(duì)經(jīng)過(guò)邊緣檢測(cè)的兩橢圓的邊緣點(diǎn)加入了噪聲，進(jìn)行曲率濾波前后的邊緣點(diǎn)曲率變化曲線分別如圖3和圖4所示。

(a) 濾波前

(b) 濾波后

(a) 濾波前

(b) 濾波后

1.4橢圓邊緣點(diǎn)的精確提取以及同心圓投影中心定位

在對(duì)橢圓邊緣點(diǎn)進(jìn)行濾波后，對(duì)橢圓進(jìn)行最小二乘擬合，設(shè)獲得橢圓的初始方程為:

(6)

根據(jù)橢圓方程可求得橢圓上任意點(diǎn)P(xp,yp)的梯度斜率kg為:

(7)

在上述的橢圓邊緣檢測(cè)中，沒(méi)有嚴(yán)格在邊緣點(diǎn)的梯度方向上對(duì)邊緣點(diǎn)進(jìn)行定位，所以將帶來(lái)一定的邊緣定位誤差。為了獲得更高精度的橢圓方程，需在橢圓邊緣點(diǎn)附近鄰域內(nèi)沿著式(7)所確定的梯度方向?qū)E圓邊緣點(diǎn)重新進(jìn)行亞像素級(jí)邊緣檢測(cè)和曲率濾波，然后對(duì)兩個(gè)橢圓的方程重新進(jìn)行最小二乘曲線擬合。

設(shè)擬合后的兩橢圓方程分別為F1(x,y)和F2(x,y)，兩橢圓的幾何中心為O1和O2。過(guò)兩橢圓的幾何中心O1和O2作一條直線L，該直線與兩橢圓的交點(diǎn)分別為A1,A2,B2,B1，如圖5所示。

圖5 同心圓投影中心定位原理圖

則同心圓的投影中心O可由仿射變換交比不變?cè)砬蟮?

(8)

2 同心圓投影中心定位實(shí)驗(yàn)

為了對(duì)文中算法的精度進(jìn)行分析，生成了同心圓測(cè)試圖像。

首先繪制帶有同心圓環(huán)和角點(diǎn)模板，然后用3*3的濾波器對(duì)該圖像進(jìn)行平滑，用以模擬攝像機(jī)拍攝中所產(chǎn)生的漸變效應(yīng)，從而生成標(biāo)準(zhǔn)測(cè)試圖像，如圖6所示。

圖6 同心圓環(huán)標(biāo)準(zhǔn)測(cè)試圖像

然后對(duì)該圖像進(jìn)行不同的單應(yīng)矩陣變換和坐標(biāo)平移，用以模擬不同角度的攝像機(jī)拍攝圖像，如圖7所示。

圖7 經(jīng)單應(yīng)變換的同心圓環(huán)測(cè)試圖像

圖7圖像的同心圓環(huán)投影中心實(shí)際值可由圖6圖像的圓環(huán)中心實(shí)際值經(jīng)單應(yīng)矩陣變換和坐標(biāo)平移得到。

文中首先分析了高斯平滑濾波器、一階高斯邊緣檢測(cè)算子對(duì)測(cè)量精度的影響，分別如圖8～圖10所示。

圖8 不同GSF尺寸的誤差分析 (SNR: 41 dB)

圖9 GSF不同σ的誤差分析 (SNR: 41 dB)

圖10 不同GEF尺寸的誤差分析

由圖8可知，當(dāng)高斯平滑濾波器的尺寸在4～10 pixels之間時(shí)，測(cè)量誤差較小；由圖9可知，當(dāng)σ≤4.2時(shí)，測(cè)量誤差較??；由圖10可知，邊緣檢測(cè)算子尺寸在11～13 pixels之間時(shí)，測(cè)量誤差較小，而且隨著噪聲的增加，測(cè)量誤差有明顯的上升趨勢(shì)。

最后，分別用文中的同心圓法、角點(diǎn)定位法對(duì)測(cè)試圖像進(jìn)行了定位誤差分析，如圖11所示。

圖11 文中方法與角點(diǎn)檢測(cè)法的像素誤差對(duì)比

由圖11可知，當(dāng)圖像中存在噪聲時(shí)，在信噪比大于35 dB時(shí)，文中方法與角點(diǎn)法的測(cè)量誤差均小于0.1 pixel;當(dāng)信噪比介于35 dB和25 dB之間時(shí)，角點(diǎn)檢測(cè)法的測(cè)量誤差突然急劇增大，而文中方法的測(cè)量誤差變化仍比較平緩，測(cè)量誤差小于0.15 pixel。

3 結(jié) 語(yǔ)

提出了采用經(jīng)過(guò)曲率濾波的橢圓曲線擬合法與射影交比不變法相結(jié)合來(lái)實(shí)現(xiàn)對(duì)同心圓的投影中心進(jìn)行的亞像素定位。并采用單應(yīng)矩陣變換和坐標(biāo)平移法生成了大量的理想測(cè)試圖像，將文中方法與角點(diǎn)檢測(cè)法的測(cè)量精度進(jìn)行了比較分析，結(jié)果表明，文中方法的測(cè)量精度在有無(wú)噪聲情況下均高于角點(diǎn)檢測(cè)方法。

[1] Jun-Sik Kim, In-So Kweon. A New Camera Calibration Method for Robotic Applications[C]// International Conference on Intelligent Robots and Systems.2001:778-783.

[2] Ashutosh Morde, Mourad Bouzit. Lawrence Rabiner[C]//A Novel Approach to planar camera calibration, " International Conference on Computer Vision Theory and Applications.2006:87-92.

[3] X H Ying, H B Zha. Efficient detection of projected concentric circles using four intersection points on a secant Line[C]//International Conference on Pattern Recognition.2008：1-4.

[4] Junpeng Xue, Xianyu Su, Liqun Xiang, et al. Using concentric circles and wedge grating for camera calibration[J]. Applied Optics,2012,51(17):3811-3816.

[5] Chris Harris, Mike Stephens. A combined corner and edge detector[C]//Proceeding of the 4th Alvey Vision Conference.1988:147-151.

[6] Li L F, Feng Z R, He K L. A randomized algorithm for detecting multiple ellipses based on least square approach [J]. Opto-electronics Review,2005,13(1):61-67.

[7] Halir R, Flusser J. Numerically stable direct least squares fitting of ellipses[C]// International Conference in Central Europe on Computer Graphics and Visualization.1998:125-132.

[8] 王毅，徐成杰.基于多目標(biāo)規(guī)劃的飛機(jī)使用計(jì)劃模型研究[J].長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2009,30(1):73-77.

Concentric circle sub-pixel center positioning method

ZHANG Fang-cheng, LI Xiao-tian

(Changchun Institute of Optics and Fine Mechanics and Physics, Chinese Academy of Sciences, Changchun 130033， China)

A location method for concentric circles' projection center based on ellipses’curve fitting and curvature filter is put forward. The test images are generated with homography transform and coordinate shift, and the method is compared with the corner detection. The results show that the accuracy of the method is better than that of corner detection no matter how the images have noise or not. When images’ S/N ratio is between 25～35 dB, the error of the corner detection increases sharply while that of concentric circle method is only up to 0.15 pixel.

machine vision; image processings; concentric circle calibration; ellipse fitness.

2014-08-22

國(guó)家重大科研裝備研制基金資助項(xiàng)目(ZBYZ2008-1)

張方程(1958-)，男，漢族，河北寧津人，中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所工程師，主要從事衍射光柵鍍膜工藝及檢測(cè)技術(shù)方向研究,E-mail:zhangfc_ciomp@163.com.

TP 311.5

A

1674-1374(2014)05-0500-06