改進的Zernike 矩工業(yè)CT 圖像邊緣檢測

陶 李，王 玨 ，鄒永寧，伍立芬，王慧倩

(重慶大學(xué) 自動化學(xué)院 光電技術(shù)及系統(tǒng)教育部重點實驗室ICT 研究中心，重慶400044)

1 引 言

近年來，隨著工業(yè)CT 技術(shù)的迅速發(fā)展，工業(yè)CT 系統(tǒng)已不僅用于對物體的無損探傷和對物體內(nèi)部結(jié)構(gòu)的定性分析和質(zhì)量評價，還廣泛應(yīng)用于物體內(nèi)部結(jié)構(gòu)尺寸和材料密度的定量測量，因此提高測量精度和測量速度以及自動化程度已成為工業(yè)CT 領(lǐng)域當前研究的重點之一。通過工業(yè)CT圖像測量物體的幾何尺寸主要包括: 點與點之間的距離、目標面積、表面積和體積等。

通過工業(yè)CT 圖像實現(xiàn)對被測物體幾何尺寸實時精確的測量，首先要提高工業(yè)CT 圖像的邊緣定位精度和速度，這直接影響后續(xù)測量結(jié)果的精度和速度。提高圖像測量精度和速度主要通過:(1) 提高工業(yè)CT 系統(tǒng)的硬件性能，如提高掃描系統(tǒng)速度和探測系統(tǒng)精度等; ( 2) 利用軟件方法，如提高圖像邊緣定位算法的精度和速度等。前者設(shè)備成本高、周期長而且在精度和速度的提高上受限; 后者設(shè)備成本低、速度快、易于更新。如果能夠?qū)T 圖像邊緣定位在亞像素級別，則可直接提高測量系統(tǒng)精度。

亞像素邊緣定位算法具有邊緣定位精度高等優(yōu)點，常用于亞像素邊緣檢測當中，方法主要包括: 基于一階微分期望值、小波變換［2］、空間矩［3］、插值、擬合［4］等，其中，Haralick［5］等人提出的基于Facet 模型的亞像素級邊緣檢測算法，通過3 次多項式擬合Facet 模型，再根據(jù)方向?qū)?shù)信息求得二階導(dǎo)數(shù)過零點，最終確定亞像素邊緣位置。該方法穩(wěn)定可靠，能定位圖像微小細節(jié)，從而提高檢測精度，但是因為運算量巨大而耗時［6］，不能同時滿足精度和速度的要求［7］。利用基于Zernike 正交復(fù)數(shù)矩［8］的亞像素邊緣檢測算法在單位圓內(nèi)具有正交性，可減少亞像素邊緣檢測所需要的模板數(shù)目、降低函數(shù)的階數(shù)、增強抗干擾能力，從而提高邊緣定位精度［9］，但是運算耗時不能滿足實時在線測量要求。因此，為了提高工業(yè)CT 圖像邊緣檢測算法的精度和速度，本文研究了改進的Zernike 矩邊緣檢測算法，該方法不僅能夠高精度定位工業(yè)CT 圖像邊緣，而且能夠大幅度提高定位速度。通過高精度定位圖像的亞像素邊緣，保證了后續(xù)圖像測量的精度;在保證精度的前提下，算法大幅減少運算時間，提高檢測速度，實現(xiàn)了實時在線工業(yè)CT 圖像的高精度測量。

2 Zernike 矩邊緣檢測原理

2．1 Zernike 矩原理

根據(jù)復(fù)數(shù)域Zernike 多項式的定義，可得到具有正交、旋轉(zhuǎn)不變特性的Zernike 矩。Zernike 的n階多項式定義為［10］:

式中:n≥0，且n- |m|為偶數(shù)，實值多項式Rnm由下式給出:

Zernike 多項式在單位圓x2+y2=1 內(nèi)是正交的，即:

當且僅當n=p，m=q，公式( 3) 成立。其中* 表示復(fù)共軛。

一幅連續(xù)圖像f(x，y) 的n階m次Zernike 矩定義為:

式中:Vmn( ρ，θ) =Rmn( ρ) eimθ為積分核函數(shù)，* 為復(fù)共軛，ρ 為圓點到點(x，y) 的矢量長度，θ 角為矢量長度ρ 與x軸之間逆時針方向的夾角。Rmn( ρ) 是徑向多項式，定義為:

對離散的數(shù)字圖像而言，其表示形式為:

Zernike 矩的最重要性質(zhì)是旋轉(zhuǎn)不變性［10］，一幅圖像旋轉(zhuǎn)θ 角度后，A'mn=Amne-imθ，即圖像旋轉(zhuǎn)前后的Zernike 矩只改變相位角，而幅度保持不變。

2．2 Zernike 矩邊緣檢測算法

Ghosal［11］等人提出基于Zernike 矩的邊緣檢測算法，其基本原理是通過建立Zernike 矩與圖像邊緣檢測理想模型的4 個邊緣參數(shù)的關(guān)系，分別求解矩而得到圖像模型4 個邊緣參數(shù)，在將參數(shù)與預(yù)設(shè)的閾值進行比較判定為邊緣點之后，通過計算來精確定位圖像邊緣位置。利用Zernike 矩的旋轉(zhuǎn)不變性旋轉(zhuǎn)圖像后，可以容易地計算出4 個參數(shù)，進而實現(xiàn)對邊緣的精確定位。Zernike 矩是積分型算子，對噪聲不敏感，抗噪能力較強。

二維平面亞像素邊緣檢測理想邊緣模型如圖1所示。理想的采樣區(qū)域為單位圓，直線L被單位圓包含的部分代表理想邊緣，圓內(nèi)L兩側(cè)的灰度值分別為H和H+k，k為兩側(cè)的灰度差，l表示原點到邊緣的理論距離值，θ 表示線段l與x軸之間的夾角。

圖1 二維平面亞像素邊緣檢測理想模型Fig.1 Ideal model of two-dimensional sub-pixel edge detection

若把圖像邊緣順時針旋轉(zhuǎn)θ 角，線段l將與x軸平行，則:

式中，f'(x，y) 是圖像經(jīng)過順時針旋轉(zhuǎn)后的得到的邊緣函數(shù)。利用Zernike 矩進行邊緣定位時，需要用到3 個不同階次的Zernike 矩，分別是A00，A11，A20。這些不同階次的Zernike 矩相對應(yīng)的積分核函數(shù)分別為:V00=1，V11=x+iy，V20=2x2+2y2-1。

旋轉(zhuǎn)后圖像的Zernike 矩與原圖像的Zernike矩之間的關(guān)系為:

由于式(7) 是A＇11的虛部，因此有:

式中Re［A11］，Im［A11］分別為A11的實部和虛部。因此可以求出θ 角，即:

式中，(xs，ys) 是邊緣的亞像素坐標。

綜上所述，利用3 個Zernike 矩A00，A11，A20可以求解理想邊緣模型的4 個參數(shù)θ，l，k，h，并可以根據(jù)4 個參數(shù)計算出邊緣點的亞像素坐標。

2.2.1 Zernike 矩算子的7 ×7 模板

采用Zernike 矩算子進行圖像邊緣檢測，需要對每個邊緣點計算得到4 個參數(shù)，用來判斷該像素點是否為邊緣點。計算不同階次的Zernike 矩模板是對圖像進行亞像素邊緣檢測的前提條件。原Ghosal 算法僅僅推導(dǎo)了Zernike 矩的5 ×5 模板中A00，A10，A20的系數(shù)，在Ghosal 等人的5 ×5 模板基礎(chǔ)上，高世一等人計算了Zernike 矩的7 ×7 模板系數(shù)［13］，這是因為實際CT 圖像邊緣不是理想階躍邊緣，邊緣存在灰度過渡，而采用7 ×7 模板能夠細化邊緣，從而更好地抑制CT 系統(tǒng)帶來的噪聲。

2.2.2 消除模板效應(yīng)

仔細分析Ghosal 算法能夠發(fā)現(xiàn)，Ghosal的Zernike矩亞像素定位精度雖然較高，卻忽略了模板效應(yīng)導(dǎo)致的邊緣坐標系數(shù)計算不精確問題。模板效應(yīng)是指Zernike 矩算子根據(jù)實際需求要選擇不同大小的模板，而使邊緣坐標計算產(chǎn)生偏差［14］。假設(shè)所選的模板大小為N×N，由于邊緣點坐標公式的推導(dǎo)和計算都在單位圓內(nèi)進行，而在實際邊緣檢測過程中，模板是從左到右、從上到下移動著和圖像進行卷積計算，這時模板每次移動覆蓋了被檢圖像上N2個像素點，單位圓的半徑則變?yōu)镹/2，因此需要把計算出來的到邊緣的垂直距離l放大為Nl/2，相應(yīng)的式(12) 變換為:

3 算法改進

由于Zernike 矩邊緣算法需要將圖像和模板進行卷積計算，運算量隨著圖像的大小而增大。雖然它檢測精度高，但算法的運行時間較慢，不能滿足對工業(yè)CT 圖像邊緣實時在線檢測的要求，因此，本文提出了綜合Sobel 的改進Zernike 矩工業(yè)CT 圖像邊緣檢測算法。主要流程如圖2 所示。

3．1 Sobel 算子邊緣定位

首先利用Sobel 算子粗定位整個CT 圖像可能的邊緣點;

圖2 算法流程圖Fig.2 Algorithm flowchart

Sobel 邊緣檢測的主要目的是找出全部邊緣可能存在的位置，在Sobel 算子粗定位之前，需要對CT 圖像進行預(yù)處理，用來減少粗檢時產(chǎn)生多余的偽邊緣。

利用Sobel 邊緣算子檢測圖像邊緣的步驟如下:

(1) 對CT 圖像上灰度值為f(x，y) 的點(x，y)在x，y方向上分別計算灰度值的偏導(dǎo)數(shù)f'x與f'y，公式如下:

式中:f'x(x，y) ，f'y(x，y) 分別表示x，y方向的一階微分，G［f(x，y) ］為梯度。

(3) 設(shè)閾值為t，當G［f(x，y) ］＞t時，即可判定該點即為圖像可能的邊緣點。然后將這些可能的邊緣點保存在數(shù)組中。

3．2 Zernike 矩邊緣精確定位

在對整個圖像使用Sobel 算子進行粗檢后，利用Zernike 矩算子對保存在數(shù)組中可能的邊緣點圖像重新精確定位圖像邊緣，步驟如下:

(1) 將像素點和計算的7 ×7 模板進行卷積;

(2) 計算每一個像素點的邊緣角度θ，k，l;

(3) 如果像素點的參數(shù)滿足判定條件k≥kt且l≤lt;則利用式(13) 計算亞像素邊緣坐標。

3.3 邊緣點判據(jù)

從3.2 步驟(3) 可以看出，在Zernike 矩算子中利用式(9) 和式(11) 求出4 個參數(shù)之后，需要對參數(shù)判斷是否滿足k≥kt且l≤lt，其中kt和lt分別是階躍強度閾值和圓心到邊緣的距離閾值，只有當上面的條件成立時，才能將像素點判定為邊緣點并進行計算。通過分析可知，若強度閾值太小，則產(chǎn)生的偽邊緣較多，強度閾值太大雖然能抑制偽邊緣，但是會丟失部分真實邊緣，效果也不理想。由式(11) 可知，k是關(guān)于A＇11和l的函數(shù)，因此可以得到k的關(guān)系式:

而其中:

令A(yù)'11=amp，則k=1．5(1 -l2)-3/2amp=λ·amp

式中，λ =1．5( 1 -l2)-3/2，當l取值較小時，則有:

為了細化邊緣圖像，通常取l為＜0.3 的正數(shù)，l與λ 的關(guān)系見表1。因此，使用amp代替k沒有損失圖像邊緣檢測精度，其主要目的是減少算法的運算量。

表1 l 與λ 之間的關(guān)系Tab．1 Relationship between l and λ

4 試驗與結(jié)果分析

4.1 改進算法的邊緣檢測精度

利用Matlab2008Ra 工具驗證加速算法不損失圖像邊緣定位精度。實驗所用的為一幅800 pixel×800 pixel 工業(yè)CT 圖像［16］，如圖3 所示。

圖3 工業(yè)CT 圖像邊緣檢測結(jié)果比較Fig.3 Comparison of industrial CT image edge detection results

從圖3( b) 可以看出，Canny 邊緣算子檢測出的邊緣圖像為像素級，經(jīng)過非極大值抑制后的圖像邊緣為單像素，但是精度達不到亞像素級別，滿足不了測量需求; 圖3( c) 利用改進的Zernike 矩檢測出亞像素邊緣，但將亞像素邊緣圖像局部放大后( 見圖3( d) ) ，由于邊緣坐標為浮點數(shù)，可能會有小的斷裂［17］。

通過對邊緣圖像進行邊緣跟蹤連接［18］，能夠計算出圓半徑參數(shù)。表2、表3 列出了圖4 中圓半徑參數(shù)測量結(jié)果以及測量值與真實值的絕對誤差和相對誤差。

圖4 運算速度檢測圖像Fig.4 Images of processing speed detection

表2 Sobel+Zernike 矩工業(yè)CT 圖像中圓半徑參數(shù)的測量Tab．2 Circle radius parameter measurement of Sobel+Zernike moment industrial CT images

表3 Zernike 矩工業(yè)CT 圖像中圓半徑參數(shù)測量Tab．3 Circle radius parameter measurement of Zernike moment industrial CT images

由表2、表3 可以看出，測量值較真實值偏大，使用Zernike 矩的半徑測量值絕對誤差＜0.24 pixel，相對誤差＜0.64%; 使用改進的Zernike 矩半徑測量值絕對誤差＜0.24，相對誤差＜0.55%;可以看出改進的Zernike 矩算法在測量精度上和原Zernike 矩算法一致，沒有降低精度。

4.2 改進算法的運算速度

為了驗證加速算法的有效性，對加速算法的運算時間進行了測量。運行環(huán)境為個人計算機( 配 置CPU 為DualCore AMD Athlon II X 2 2 5 0 3GHz，內(nèi)存為2G) 和Matlab2008。仿真了4 個圖像中不同的邊緣點數(shù)量，在不同圖像尺寸( 256 pixel×256 pixel、512 pixel × 512 pixel、1 024 pixel×1 024 pixel、2 048 pixel×2 048 pixel)中測量了運算時間，分析了邊緣點數(shù)量對算法運算速度的影響。

從表4 中可以得出，Zernike 矩算子的邊緣檢測時間隨著圖像尺寸大小而劇增，改進的Zernike 矩算子能夠大幅度減少運算時間。

表4 Sobel、zernike 矩和改進的Zernike 矩運算時間Tab．4 Computation time of Sobel，Zernike moments and improved Zernike moments (s)

圖5 圖像尺寸、邊緣數(shù)量與加速比Fig.5 Image size，edge quantity and acceleration ratio

由圖5 可以看出，圖像尺寸越大，改進算法的加速比越大;圖像邊緣點數(shù)量多，加速比越小;邊緣點數(shù)量越少，改進算法的加速比越大，加速效果也會越明顯。在實際CT 圖像中，邊緣數(shù)量比例小于5%，因此改進算法可提高實際CT 圖像的邊緣檢測速度。

改進的Zernike 矩算子運算時間較Zernike 矩算子快的主要原因在于該算子先利用快速的Sobel算子對整個圖像進行了粗檢測，縮小Zernike矩算子檢測范圍，再利用7 ×7 大模板的Zernike矩算子精確定位邊緣。

綜上所述，圖像中邊緣點個數(shù)決定了整個算法的運算時間。在一幅圖像中當邊緣點比例為n≈87%時，兩個算子的運行時間才相等。然而在任何情況下邊緣點比例n≤50%?？紤]到Soble算子粗檢出可能的邊緣點會多于實際邊緣點數(shù)量，如果Sobel 算子粗檢出的可能邊緣點數(shù)量為真實邊緣點數(shù)量的2 倍，則n≈48%。在實際工業(yè)CT 圖像邊緣檢測中，由于存在噪聲和灰度過渡，Sobel 算子會粗檢出比真實邊緣多的邊緣點，驗證了改進后的Zernike 矩算法能夠大幅提高運算速度。

5 結(jié) 論

研究了一種改進的Zernike 矩的工業(yè)CT 亞像素邊緣檢測算法，并進行了工業(yè)CT 圖像實驗。實驗結(jié)果表明，該方法的邊緣定位精度無損失，運算速度快。運用該算法能夠精確定位工業(yè)CT 圖像邊緣亞像素位置，在定位精度上和Zernike 矩保持一致，半徑測量值絕對誤差＜0.1 mm，相對誤差＜0.5%;在速度上比單獨使用Zernike 矩算子快，在實際CT 圖像中，檢測速度提高了70%。因此該算法適用于對檢測精度和實時在線性要求較高的工業(yè)CT 圖像測量。

［1］ 莊天戈.CT 原理與算法［M］.上海:上海交通大學(xué)出版社，1992.ZHUANG T G.CT Theory and Arithmetic［M］. Shanghai:Shanghai Jiaotong University Press，1992.( in Chinese)

［2］ 林曉梅，李琳娜，牛剛，等.基于小波邊緣檢測的圖像去噪方法［J］.光學(xué) 精密工程，2004，12(1) :88-93.LIN X M，LI L N，NIU G，et al.. Image de-noising based on wavelet edge detection［J］.Opt． Precision Eng，2004，12(1) :88-93.( in Chinese)

［3］ 郭玉波，姚郁，遆曉光.基于空間矩的圓標記中心亞像素定位算法［J］.吉林大學(xué)學(xué)報，2009，39(1) :160-163.GUO Y B，YAO Y，DI X G. Subpixel location algorithm for circle target center based on spatial moment［J］.J． Jilin University( Engineering and Technology Edition) ，2009，39(1) :160-163.( in Chinese)

［4］ 尚雅層，陳靜，田軍委.高斯擬合亞像素邊緣檢測算法［J］.計算機應(yīng)用，2011，31(1) :179-181.SHANG Y C，CHEN J，TIAN J W. Sub-pixel edge detection algorithm based on Gauss fitting［J］.J． Computer Appl.，2011，31(1) :179-181.( in Chinese)

［5］ HARALICKR M. Digital step edges from zerocrossing of second directional derivatives［J］.IEEE T． Pattern Anal．，1984，6(1) :58-68.

［6］ HU Z F，DANG H SH，LI X R. A novel fast subpixel edge location method based on sobel-OFMM［C］//Proc. of the IEEE International Conference on Automation and Logistics，Sept 1-3，2008，Qingdao，China，2008

［7］ 王云慧，馬軍山，孫軍.基于邊緣檢測的七階矩方法［J］.計算機應(yīng)用，2010，30( S1) :159-161.WANG Y H，MA J SH，SUN J. Seven-order moments based on edge detection［J］.J． Computer Appl．，2010，30( S1) :159-161.( in Chinese)

［8］ 魏本征，趙志敏，華晉.基于改進形態(tài)學(xué)梯度和Zernike 矩的亞像素邊緣檢測方法［J］.儀器儀表學(xué)報，2010，31(4) :838-844.WEI B ZH，ZHAO ZH M，HUA J. Sub-pixel edge detection method based on improved morphological gradient and Zernike moment［J］.Chinese J． Sci． Instrument，2010，31(4) :838-844.( in Chinese)

［9］ XIANG L M，WAN T，HE C ZH. A stable least squares ellipses fitting algorithm based on Zernike moments［J］.Advanced Mater． Res．，2011，142:199-203.

［10］ GHOSAL S，MEHROTRA R. Orthogonal moment operators for subpixel edge detection［J］.Pattern Recognition，1993，26(2) :295-306.

［11］ Al-RAWI M. Fast Zernike moments［J］.Real-Time Image Proc．，2008，3(1-2) :89-96.

［12］ LIANG B，DONG M L，WANG J，et al. Sub-pixel location of center of target based on Zernike moment［J］.SPIE，2010，7544:75443A.

［13］ 高世一，趙明揚，張雷，等.基于Zernike 正交矩的圖像亞像素邊緣檢測算法改進［J］. 自動化學(xué)報，2008，34( 9) :1163-1168.GAO SH Y，ZHAO M Y，ZHANG L，et al.. Improved algorithm about subpixel edge detection of image based on Zernike orthogonal moments［J］.Acta Automatica Sinica，2008，34(9) :1163-1168.( in Chinese)

［14］ FANG Y H，XIAO H M，WANG B L. A sub-pixel feature extraction method for CT image based on Zernike moments［C］//2nd IEEE International Conference on Advanced Computer Control( ICACC 2010) ，Mar 27-29，2010，Shenyang，China，2010:539-542.

［15］ 曲迎東，李榮德，白彥華，等.高速的9 ×9 尺寸模板Zernike 矩邊緣算子［J］.光電子·激光，2010，21( 11) :1683-1687.QU Y D，LI R D，BAI Y H，et al.. A high-speed Zernike moments edge operator based on 9 ×9 masks［J］.J．Optoelectronics·Laser，2010，21(11) :1683-1687.( in Chinese)

［16］ 曲迎東，崔成松，陳善本，等.Soble-Zernike 矩亞像素邊緣算子的影響因素分析［J］.光電工程，2005，32(7) :71-73.QU Y D，CUI CH S，CHEN SH B，et al.. Analysis of the influencing factors on subpixel Sobel-Zernike moments edge operator［J］.Opto-electronic Eng.，2005，32(7) :71-73.( in Chinese)

［17］ 劉豐林，喬桂鋒，鄒斌.工業(yè)CT 圖像圓精確測量［J］.光學(xué) 精密工程，2009，11(17) :2842-2848.LIU F L，QIAO G F，ZOU B. Precise measurement of circles in industrial computed tomographic images［J］.Opt． Precision Eng.，2009，11(17) :2842-2848.( in Chinese)

［18］ 向才兵，曾理.綜合CV 和Facet 模型的裂紋邊緣檢測［J］.計算機工程與應(yīng)用，2011，47(2) :153-155.XIANG C B，ZENG L. Crack edge detection using comprehensive CV and facet model［J］.Computer Eng． Appl.，2011，47(2) :153-155.( in Chinese)

［19］ 馬睿，曾理，盧艷平.工業(yè)CT 圖像的亞像素級面積測量［J］.計算機工程與應(yīng)用，2009，45(15) :233-236.MA R，ZENG L，LU Y P. Sub-pixel area measurement method of ICT image［J］.Computer Eng． Appl.，2009，45(15) :233-236.( in Chinese)

［20］ AGNIESZKA L. Multiscale moments-based edge detection［C］//Eurocon 2009. IEEE，May 18-23，2009，St.-petersburg，USA，2009:1414-1419.

［21］ 陳娟，陳乾輝，師路歡，等.圖像跟蹤中的邊緣檢測技術(shù)［J］.中國光學(xué)與應(yīng)用光學(xué)，2009，2(1) :46-52.CHEN J，CHEN Q H，SHI L H，et al.. Edge detection technology in imaging tracking［J］.Chinese J． Optics and Appl．Optics，2009，2(1) :46-52.( in Chinese)

［22］ 袁理，葉露，賈建祿.基于Hough 變換的橢圓檢測算法［J］.中國光學(xué)與應(yīng)用光學(xué)，2010，3(4) :379-384.YUAN L，YE L，JIA J L. Ellipse detection algorithm based on Hough transforms［J］.Chinese J． Opt． and Appl． Optics，2010，3(4) :379-384.( in Chinese)