基于多元統(tǒng)計(jì)的線陣CT圖像環(huán)形偽影去除方法

，，，，

(1.中國兵器科學(xué)研究院寧波分院，寧波 315103;2.浙江工業(yè)大學(xué) 機(jī)械工程學(xué)院，杭州 310014)

基于多元統(tǒng)計(jì)的線陣CT圖像環(huán)形偽影去除方法

齊子誠1,2，倪培君1，李紅偉1，唐盛明1，郭智敏1

(1.中國兵器科學(xué)研究院寧波分院，寧波 315103;2.浙江工業(yè)大學(xué) 機(jī)械工程學(xué)院，杭州 310014)

針對(duì)線陣工業(yè)CT重建圖像中環(huán)形偽影干擾問題，設(shè)計(jì)一種基于多元統(tǒng)計(jì)的新型校正方法。首先采用灰度直方圖統(tǒng)計(jì)結(jié)合最大類間方差法對(duì)CT圖像中不同材料密度進(jìn)行分類，根據(jù)材料密度分布對(duì)原始CT圖像進(jìn)行材料去除，然后進(jìn)行極坐標(biāo)系轉(zhuǎn)換，對(duì)轉(zhuǎn)換后的極坐標(biāo)圖像進(jìn)行橫坐標(biāo)像素值統(tǒng)計(jì)，提取不同半徑r下偽影灰度幅值統(tǒng)計(jì)值，最后采用原始CT圖像消減偽影統(tǒng)計(jì)值的方式進(jìn)行偽影校正。采用6 MeV線陣探測(cè)器高能工業(yè)CT系統(tǒng)，對(duì)含有環(huán)形偽影的不同結(jié)構(gòu)CT圖像進(jìn)行校正試驗(yàn)。試驗(yàn)結(jié)果表明：該方法在有效去除CT圖像環(huán)形偽影的同時(shí)，還能較好地保持圖像細(xì)節(jié)和分辨率，顯著提高了處理速度。

工業(yè)CT；環(huán)形偽影；噪聲去除；多元統(tǒng)計(jì)

工業(yè)計(jì)算機(jī)層析技術(shù)能夠不受被檢測(cè)物體材料、形狀、表面狀況等的限制[1]，給出被檢測(cè)物體內(nèi)部結(jié)構(gòu)、組成、材料及缺損狀況的二維、三維圖像，是現(xiàn)代工業(yè)產(chǎn)品無損缺陷檢測(cè)的一項(xiàng)重要手段[2]。現(xiàn)階段應(yīng)用較為廣泛的是第三代掃描方法，其采用單射線源，具有大扇角、寬扇束、全包容被檢斷面的特點(diǎn)，且可控性強(qiáng)、效率高。環(huán)形偽影作為三代掃描模式中較為常見的一種偽影，不僅會(huì)影響圖像的質(zhì)量，還會(huì)對(duì)后續(xù)的處理、缺陷識(shí)別和尺寸測(cè)量造成干擾。因此，研究CT圖像環(huán)形偽影校正方法是十分必要的。

目前，少數(shù)CT系統(tǒng)選用高性能探測(cè)器抑制環(huán)形偽影，但圖像處理仍是主要的去環(huán)形偽影的技術(shù)手段[3]?；趫D像處理技術(shù)消除CT環(huán)形偽影的方法主要分為投影正弦圖校正法和極坐標(biāo)校正法兩類[4-10]。投影正弦圖校正法在投影域直接對(duì)投影圖像進(jìn)行濾波，然后對(duì)恢復(fù)后的投影數(shù)據(jù)采用濾波反投影(FBP)算法進(jìn)行重建，實(shí)現(xiàn)CT圖像環(huán)形偽影的校正。極坐標(biāo)校正法需要將笛卡爾坐標(biāo)下的環(huán)形偽影變換為極坐標(biāo)下的直線偽影，然后進(jìn)行水平方向平滑濾波，最后再變換到直角坐標(biāo)系中。上述校正方法對(duì)環(huán)形偽影都有不錯(cuò)的抑制效果，但也存在一定的不足。基于投影正弦圖的校正方法雖然在質(zhì)量和保真度上都要優(yōu)越一些，但是存在迭代算法處理等相關(guān)參數(shù)不易確定，閾值不易選擇的問題。極坐標(biāo)校正法由于兩次坐標(biāo)變換都要用到插值，并且需要針對(duì)各個(gè)不同的圖像設(shè)定不同的環(huán)形偽影濾波器的閾值，當(dāng)插值和閾值設(shè)定不妥當(dāng)時(shí)，容易造成圖像分辨率下降、細(xì)節(jié)丟失、運(yùn)算效率低下等問題。

針對(duì)上述方法的不足，筆者提出了基于多元統(tǒng)計(jì)的環(huán)形偽影去除方法，首先對(duì)原始圖像進(jìn)行直方圖統(tǒng)計(jì)，采用多項(xiàng)最大類間方差分割法確定多種密度分類，并確定每一類的主元，刪去圖像中各個(gè)主元灰度值，再進(jìn)行極坐標(biāo)變換，統(tǒng)計(jì)偽影線分布系數(shù)，對(duì)原始圖像進(jìn)行映射，然后用原始數(shù)據(jù)減去偽影系數(shù)實(shí)現(xiàn)對(duì)異常像素的校正。通過對(duì)實(shí)際CT圖像進(jìn)行試驗(yàn)驗(yàn)證，結(jié)果表明該方法無需進(jìn)行線條識(shí)別，在去除環(huán)形偽影的同時(shí)能保持圖像的細(xì)節(jié)，有效快速地校正CT圖像中的環(huán)形偽影。

1 基于多元統(tǒng)計(jì)的環(huán)形偽影校正

1.1 環(huán)形偽影成因

環(huán)形偽影是CT中常見的一種偽影，產(chǎn)生原因有X射線束硬化、重建過程中濾波器和相鄰探測(cè)器響應(yīng)不一致等。原始數(shù)據(jù)在經(jīng)過反投影算法后，將在每個(gè)視角方向上都產(chǎn)生一條直線，通過持續(xù)的累加會(huì)形成包絡(luò)，同時(shí)由于離散反投影需要插值運(yùn)算，這往往會(huì)加大此包絡(luò)的寬度，形成環(huán)形偽影，環(huán)形偽影形成原理示意與鋁合金CT圖像環(huán)形偽影示例如圖1所示。

圖1 環(huán)形偽影形成原理示意與鋁合金CT圖像環(huán)形偽影示例

在重建圖像上，CT環(huán)形偽影與被檢測(cè)對(duì)象往往是疊加的，表現(xiàn)為一系列的同心圓，并沿徑向半徑不斷擴(kuò)大。環(huán)形偽影的存在很大程度上降低了圖像的質(zhì)量，給基于圖像的測(cè)量、識(shí)別、噪聲處理以及圖像分割等進(jìn)一步的處理和分析帶來很大的困擾。

1.2 多元直方圖分類統(tǒng)計(jì)

多閾值最大類間方法的基本思想是在兩個(gè)波峰之間找到一個(gè)統(tǒng)計(jì)意義上的最佳波谷，用波谷劃分開兩個(gè)獨(dú)立的波峰，從而區(qū)分開目標(biāo)與背景。設(shè)包含環(huán)形偽影的CT圖像大小為N×N，該圖像上任意一點(diǎn)的像素值可以表示為I(x,y)，其中x，y分別為該點(diǎn)的橫坐標(biāo)和縱坐標(biāo)，0≤x，y≤N。圖像中包含M個(gè)灰度級(jí)(0,1…,M-1)，灰度值為i的象素點(diǎn)數(shù)為Ni，由式(1)計(jì)算出圖像總像素點(diǎn)數(shù)W。

由式(2)計(jì)算灰度值為i的圖像點(diǎn)數(shù)占比。

設(shè)CT圖像中由n種不同密度材料組成，通過不同灰度分布可把直方圖大致分割成L0,L1,…,Ln個(gè)類。當(dāng)鄰近的兩組數(shù)據(jù)的類間方差最大時(shí)，得到最佳分割閾值ti。分割后各類的均值、方差分別為ui,σi。由式(3)計(jì)算類間方差σi。

根據(jù)最大類間方差計(jì)算n個(gè)類之間的最佳分割閾值，重新分割各類，設(shè)新類為A(i)，0

1.3 主元消減法

根據(jù)坐標(biāo)系I(x,y)中，每個(gè)像素所在類進(jìn)行分類，減去該點(diǎn)所在類的灰度均值，記為I′(x,y)。設(shè)直角坐標(biāo)I(x,y)中任意一點(diǎn)為k，其可表示為I(xk,yk)，首先判斷k點(diǎn)的灰度屬于哪個(gè)類，設(shè)k∈(tm-1,tm)即屬于類A(m)，該類的均值為Tm，則I′(xk,yk)=I(xk,yk)-Tm替換原來k在極坐標(biāo)系θ-r所在位置的灰度值。對(duì)于灰度域中兩端點(diǎn)數(shù)據(jù)的處理方法為，當(dāng)m=1時(shí)，即類A(1)灰度范圍(0,t1)，則I′=I(xk,yk)-max[WM,M∈(0,t1)]。將直角坐標(biāo)系中的所有點(diǎn)進(jìn)行上述處理得到I′(x,y)。

1.4 CT圖像坐標(biāo)變換

設(shè)極坐標(biāo)圖像為I(r,θ)，θ∈(0,2π),r∈(0,N/2)。θ取值間隔為360/(π×N)。從CT圖像(直角坐標(biāo)系)映射到θ-r圖像(極坐標(biāo)系)，由式(5)計(jì)算。

式中：r為正整數(shù)，范圍0≤r≤N/2，計(jì)算過程中四舍五入至最近正整數(shù)；θ數(shù)值最小間隔為360/(π×N)，計(jì)算過程中四舍五入至最近間隔位。

當(dāng)CT圖像(直角坐標(biāo))中所有的I(x,y)運(yùn)算完成以后，在θ-r圖像(極坐標(biāo))中會(huì)出現(xiàn)以下三種情況：

(1)θ-r圖像(極坐標(biāo))中某一坐標(biāo)點(diǎn)(r1,θ1)被賦值1次，則當(dāng)前值即為該位置的值。

(2)θ-r圖像(極坐標(biāo))中某一坐標(biāo)點(diǎn)(r1,θ1)被賦值多次，采用多次累加得到的和除以賦值次數(shù)，求得平均值即為該坐標(biāo)點(diǎn)的值。上述方法采用的角度θ間隔為360/(π×N)，最大程度上減少了由于疊加平均而引起的圖像細(xì)節(jié)的丟失。

圖2 插值優(yōu)先示意

(3)θ-r圖像(極坐標(biāo))中某一坐標(biāo)點(diǎn)(r1,θ1)未被賦值，則當(dāng)前值需要通過在θ-r圖像(極坐標(biāo))上進(jìn)行，插值優(yōu)先級(jí)從大到小為(r2,θ2)>(r3,θ3)>…>(r8,θ8)>(r9,θ9)。如圖2所示，通過上述步驟實(shí)現(xiàn)從CT圖像(直角坐標(biāo))映射到θ-r圖像(極坐標(biāo))的轉(zhuǎn)換。

1.5 偽影均值統(tǒng)計(jì)

1.6 偽影校正

設(shè)去除環(huán)形偽影后的灰度值為G(x,y)，圖像中心位置為(N/2,N/2)，由式(6)計(jì)算該點(diǎn)(x,y)與圖像中心的距離Lx,y。

由式(7)計(jì)算G(x,y)各個(gè)點(diǎn)的灰度值。

式中：Sh為濾波參數(shù)Sr，其中h=N/2-Lx,y。

對(duì)所有I(x,y)全部點(diǎn)進(jìn)行校正處理，獲得G(x,y)。

2 試驗(yàn)設(shè)備及對(duì)象

2.1 試驗(yàn)設(shè)備

如圖3所示，試驗(yàn)設(shè)備為北京固鴻生產(chǎn)的6 MeV高能工業(yè)CT系統(tǒng)，其空間分辨率為2 lp·mm-1，重建圖像像素尺寸為4 096×4 096，重建視場(chǎng)半徑為100 mm。采用厚度為1 mm的切片進(jìn)行對(duì)比試驗(yàn)，微動(dòng)次數(shù)均10次。采用三代掃描方式。

圖3 6 MeV高能工業(yè)CT檢測(cè)系統(tǒng)

2.2 試驗(yàn)對(duì)象

試驗(yàn)對(duì)象為均質(zhì)鋁合金試塊，如圖4所示。由于工藝參數(shù)選擇不當(dāng)，工件CT圖像含有明顯環(huán)形偽影。選取局部圖像進(jìn)行放大處理，可見環(huán)形偽影和高密度夾雜。

圖4 鋁合金試塊CT圖

3 試驗(yàn)結(jié)果與分析

對(duì)含有環(huán)形偽影的CT圖像進(jìn)行偽影去除，如圖4(a)所示的CT圖像像素尺寸為4 096×4 096。圖5為直方圖統(tǒng)計(jì)后計(jì)算的各類的均值，進(jìn)行主元消減，從圖中可以發(fā)現(xiàn)含有較多的環(huán)形偽影，有多個(gè)相連的全環(huán)，環(huán)形偽影較原圖更加明顯。進(jìn)行極坐標(biāo)展開后的圖像如圖6所示，圖中橫線為偽影的位置，可以發(fā)現(xiàn)極坐標(biāo)圖中偽影信息基本都得到了加強(qiáng)。統(tǒng)計(jì)圖6中極坐標(biāo)系CT圖像每行偽影的灰度均值，建立不同中心距與偽影灰度幅值的對(duì)應(yīng)關(guān)系，如圖7所示。根據(jù)圖4(a)原圖中每個(gè)像素點(diǎn)的位置，對(duì)應(yīng)偽影灰度幅值，進(jìn)行減法校正形成最終偽影去除后的CT圖像，如圖8所示。

圖5 鋁合金試塊的主元消減后CT圖

圖6 鋁合金試塊的極坐標(biāo)展開圖像

圖7 鋁合金試塊的偽影灰度分布圖

圖8 處理后的鋁合金試塊CT圖

結(jié)合復(fù)雜結(jié)構(gòu)產(chǎn)品對(duì)文中的偽影去除效果進(jìn)行定量分析，采用高碳高鉻鋼零部件作為檢測(cè)對(duì)象，如圖9(a)所示。使用文中的多元統(tǒng)計(jì)方法、正弦投影方法[11]、極坐標(biāo)變換方法[12]對(duì)圖9(b)進(jìn)行環(huán)形偽影校正。并且，采用信噪比作為圖像質(zhì)量評(píng)價(jià)指標(biāo)，分別選取工件內(nèi)部、背景和中間部位三個(gè)局部區(qū)域(局部1～局部3)進(jìn)行效果比較。采用圖像的信噪比作為圖像質(zhì)量的定量評(píng)價(jià)指標(biāo)，如式(8)所示。

計(jì)算背景區(qū)域中的信噪比，得到不同校正方法灰度圖像的信噪比如表1所示。

從圖10和表1可看出，多元統(tǒng)計(jì)方法、正弦投影方法和極坐標(biāo)變換方法對(duì)CT圖像環(huán)形偽影的校正都有一定的抑制效果。但是，正弦投影方法的邊緣細(xì)節(jié)保護(hù)效果好，圖像信噪比低；極坐標(biāo)變換方法圖像視覺上模糊且分辨率低。綜合考慮，多元統(tǒng)計(jì)方法可減小均勻區(qū)域標(biāo)準(zhǔn)差，提高信噪比，具有較好的圖像細(xì)節(jié)保持能力，且只需進(jìn)行一次坐標(biāo)系變化，極大地縮短了處理時(shí)間。

圖9 零件實(shí)物及其CT圖

圖10 零件局部區(qū)域(局部1～局部3)的試驗(yàn)對(duì)比效果圖

灰度圖像局部1圖像信噪比局部2圖像信噪比局部3圖像信噪比原局部圖24.6628.0524.53多元統(tǒng)計(jì)方法24.9131.8824.84正弦投影方法24.7130.4524.75極坐標(biāo)變換法24.8230.1624.63

4 結(jié)語

提出了一種基于多元統(tǒng)計(jì)的線陣CT圖像環(huán)形偽影去除方法，根據(jù)CT圖像直方圖上多類密度區(qū)域分布的特點(diǎn)，采用多閾值類間方差法對(duì)CT圖像灰度區(qū)域進(jìn)行分割，以各個(gè)區(qū)域內(nèi)主元作為該區(qū)域像素的均值進(jìn)行減除，將減除后的CT圖像作為偽影統(tǒng)計(jì)的輸入圖像。理論分析和試驗(yàn)測(cè)試結(jié)果均驗(yàn)證了該方法對(duì)工業(yè)CT圖像環(huán)形偽影去除的有效性。相比較于正弦投影方法，該方法能減小均勻區(qū)域的標(biāo)準(zhǔn)差，提高信噪比，具有較好的圖像細(xì)節(jié)的保持能力；較傳統(tǒng)極坐標(biāo)變換法，其只需進(jìn)行一次坐標(biāo)系變化而極大地縮短了處理時(shí)間；在保證較高的尺寸測(cè)量精度的同時(shí)，較大地提高了處理效率。

[1] 張朝宗,郭志平,張朋，等. 工業(yè)CT技術(shù)和原理[M].北京:科學(xué)出版社,2009:32-80.

[2] 張俊哲.無損檢測(cè)技術(shù)及其應(yīng)用[M].北京:科學(xué)出版社,2010.

[3] GONZALEZ R C, WOODS R E. 數(shù)字圖像處理[M].北京:電子工業(yè)出版社, 2008.

[4] 馬繼明,宋巖,王群書,等.X射線CT環(huán)形偽影去除方法[J].強(qiáng)激光與粒子束,2014,26(12):177-182.

[5] 王玨,黃蘇紅,蔡玉芳,等.改進(jìn)Canny算法的CT圖像環(huán)形偽影校正[J].光學(xué) 精密工程,2011,19(11):2767-2773.

[6] 周意超,謝明元,楊玲,等.CT環(huán)狀偽影矯正方法的改進(jìn)研究[J].四川大學(xué)學(xué)報(bào)(醫(yī)學(xué)版),2016(3):420-424.

[7] 郭宏,曾棟,張華,等.基于投影域小波濾波處理的CT圖像環(huán)形偽影去除方法[J].南方醫(yī)科大學(xué)學(xué)報(bào),2015(9):1258-1262.

[8] POLUDNIOWSKI G, EVANS P M, HANSEN V N, et al. An efficient Monte Carlo-based algorithm for scatter correction in keV cone-beam CT[J]. Physics in Medicine and Biology,2009,54:3847-3864.

[10] 張華,黃魁東,史儀凱,等. 一種基于空氣掃描的錐束CT環(huán)形偽影校正方法[J].CT理論與應(yīng)用研究,2012,21(2):247-254.

[11] 李保磊,楊民,傅健,等.兩種CT成像環(huán)狀偽影校正方法[J].光學(xué)學(xué)報(bào),2009,29(7):1849-1853.

[12] 楊俊,甄鑫,盧文婷,等.基于圓掃描軌跡的錐形束CT重建與環(huán)形偽影消除[J].南方醫(yī)科大學(xué)學(xué)報(bào),2009,29(12):2379-2382.

IndustrialCTImagesRingArtifactsRemovalMethodBasedonMultivariateStatistics

QI Zicheng1,2, NI Peijun1, LI Hongwei1, TANG Shengming1, GUO Zhimin1

(1.The Ningbo Branch of Ordnance Science Institute of China, Ningbo 315103, China;2.College of Mechanical Engineering, Zhejiang University of Technology, Hangzhou 310014, China)

A multivariate statistical method was proposed to correct ring artifacts in industrial CT images. First, the grayscale histogram statistics and maximum interclass variance were used to classify the density in CT images. Then, the coordinate system was converted to the polar one and the abscissa pixel values were counted in the polar coordinate image and the artifacts of the gray scale amplitude were extracted at different radii. Finally, the correction was performed by subtracting artifacts from the original CT image. The CT images of different structures with ring artifacts were calibrated by using 6 MeV linear array detector high energy industrial CT system. The experimental results show that this method can effectively remove CT images from ring artifacts. The image details and resolution were well kept and the processing speed was significantly improved.

industrial CT; ring artifact; noise removal; multivariate statistics

TG115.28

A

1000-6656(2017)12-0020-05

2017-07-09

國家自然科學(xué)基金資助項(xiàng)目(61471411)；浙江省自然科學(xué)基金資助項(xiàng)目(LQ15E010003)；寧波國際科技合作資助項(xiàng)目(2015D10005)；寧波市自然科學(xué)基金資助項(xiàng)目(2016A610247)

齊子誠(1984-)，男，副研究員，主要從事無損檢測(cè)自動(dòng)化技術(shù)、圖像處理等方面的研究

齊子誠, nathan1984@qq.com

10.11973/wsjc201712005