基于多元統計的線陣CT圖像環形偽影去除方法

，，，，

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

基于多元統計的線陣CT圖像環形偽影去除方法

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

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

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

工業CT；環形偽影；噪聲去除；多元統計

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

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

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

1 基于多元統計的環形偽影校正

1.1 環形偽影成因

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

圖1 環形偽影形成原理示意與鋁合金CT圖像環形偽影示例

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

1.2 多元直方圖分類統計

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

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

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

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

1.3 主元消減法

根據坐標系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所在位置的灰度值。對于灰度域中兩端點數據的處理方法為，當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為正整數，范圍0≤r≤N/2，計算過程中四舍五入至最近正整數；θ數值最小間隔為360/(π×N)，計算過程中四舍五入至最近間隔位。

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

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

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

圖2 插值優先示意

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

1.5 偽影均值統計

1.6 偽影校正

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

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

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

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

2 試驗設備及對象

2.1 試驗設備

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

圖3 6 MeV高能工業CT檢測系統

2.2 試驗對象

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

圖4 鋁合金試塊CT圖

3 試驗結果與分析

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

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

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

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

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

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

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

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

圖9 零件實物及其CT圖

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

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

4 結語

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

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

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

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

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

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

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

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

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

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