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

張方程, 李曉天

(中國科學院長春光學精密機械與物理研究所， 吉林 長春 130033)

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

張方程, 李曉天

(中國科學院長春光學精密機械與物理研究所， 吉林 長春 130033)

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

機器視覺； 圖像處理； 同心圓定標； 橢圓擬合

0 引 言

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

文中提出采用經過曲率濾波的曲線擬合方法來實現同心圓投影中心的亞像素定位。首先提取兩橢圓邊界的亞像素定位坐標，并采用曲率濾波方法剔除掉誤差較大的橢圓邊界點對兩個橢圓方程進行最小二乘法擬合。然后根據兩個橢圓中心亞像素級坐標以及兩個橢圓的方程，采用仿射變換方法對同心圓的投影中心進行定位和逼近。最后在有無噪聲的情況下對文中方法和角點檢測法[6]的圖像點提取精度進行對比仿真分析，并給出了實驗測量結果。

1 數學建模

1.1圖像預處理

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

(1)

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

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

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

圖1 ROI圖像的像素級邊緣提取

1.2.1 雙線性插值

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

雙線性插值方法認為圖像的灰度在橫向和縱向方向上都是線性變化的，這與光學成像的漸變原理相一致[6-8]。設直線L1上任一點P的坐標為(xp,yp)，若xp和yp皆為整數，則P點的灰度值IP可直接從同心圓圖像中讀取，記為I(xp,yp)。若xp和yp不全為整數，則P點的雙線性插值灰度值IP為:

(2)

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

I(a,b),I(a+1,b),I(a,b+1)，I(a+1,b+1)----分別為對應的整數點坐標的灰度值。

1.2.2 一階高斯離散變換

設離散高斯變換的長度f為奇數，則離散的高斯函數一階導數的表達式為

(3)

沿著直線L1的方向對圖像進行雙線性插值并獲得一維圖像數據后,采用離散的高斯函數一階導數對一維圖像灰度數據進行卷積運算，可獲得曲線如圖2所示。

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

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

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

(4)

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

1.3曲率濾波

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

設(x1,y1),(x2,y2),(x3,y3)為橢圓邊緣上的任意連續3點，則(x2,y2)點處的曲率kr計算公式為

(5)

其中

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

(a) 濾波前

(b) 濾波后

(a) 濾波前

(b) 濾波后

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

在對橢圓邊緣點進行濾波后，對橢圓進行最小二乘擬合，設獲得橢圓的初始方程為:

(6)

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

(7)

在上述的橢圓邊緣檢測中，沒有嚴格在邊緣點的梯度方向上對邊緣點進行定位，所以將帶來一定的邊緣定位誤差。為了獲得更高精度的橢圓方程，需在橢圓邊緣點附近鄰域內沿著式(7)所確定的梯度方向對橢圓邊緣點重新進行亞像素級邊緣檢測和曲率濾波，然后對兩個橢圓的方程重新進行最小二乘曲線擬合。

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

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

則同心圓的投影中心O可由仿射變換交比不變原理求得:

(8)

2 同心圓投影中心定位實驗

為了對文中算法的精度進行分析，生成了同心圓測試圖像。

首先繪制帶有同心圓環和角點模板，然后用3*3的濾波器對該圖像進行平滑，用以模擬攝像機拍攝中所產生的漸變效應，從而生成標準測試圖像，如圖6所示。

圖6 同心圓環標準測試圖像

然后對該圖像進行不同的單應矩陣變換和坐標平移，用以模擬不同角度的攝像機拍攝圖像，如圖7所示。

圖7 經單應變換的同心圓環測試圖像

圖7圖像的同心圓環投影中心實際值可由圖6圖像的圓環中心實際值經單應矩陣變換和坐標平移得到。

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

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

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

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

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

最后，分別用文中的同心圓法、角點定位法對測試圖像進行了定位誤差分析，如圖11所示。

圖11 文中方法與角點檢測法的像素誤差對比

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

3 結 語

提出了采用經過曲率濾波的橢圓曲線擬合法與射影交比不變法相結合來實現對同心圓的投影中心進行的亞像素定位。并采用單應矩陣變換和坐標平移法生成了大量的理想測試圖像，將文中方法與角點檢測法的測量精度進行了比較分析，結果表明，文中方法的測量精度在有無噪聲情況下均高于角點檢測方法。

[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] 王毅，徐成杰.基于多目標規劃的飛機使用計劃模型研究[J].長春工業大學學報：自然科學版，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

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

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

TP 311.5

A

1674-1374(2014)05-0500-06