，，，，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 主元消減法

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

1.4 CT圖像坐標變換

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

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

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

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

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

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

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

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

1.6 偽影校正

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

由式(7)計算G(x,y)各個點的灰度值。

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

對所有I(x,y)全部點進行校正處理，獲得G(x,y)。

2 試驗設備及對象

2.1 試驗設備

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

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

2.2 試驗對象

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

圖4 鋁合金試塊CT圖

3 試驗結果與分析

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

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

圖6 鋁合金試塊的極坐標展開圖像

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

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

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

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

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

圖9 零件實物及其CT圖

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

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

4 結語

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

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

[2] 張俊哲.無損檢測技術及其應用[M].北京:科學出版社,2010.

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

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

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

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

[7] 郭宏,曾棟,張華,等.基于投影域小波濾波處理的CT圖像環(huán)形偽影去除方法[J].南方醫(yī)科大學學報,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理論與應用研究,2012,21(2):247-254.

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

[12] 楊俊,甄鑫,盧文婷,等.基于圓掃描軌跡的錐形束CT重建與環(huán)形偽影消除[J].南方醫(yī)科大學學報,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

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

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

齊子誠, nathan1984@qq.com

10.11973/wsjc201712005