基于中值濾波的口腔CT圖像環(huán)形偽影去除算法

王闖 李詩 吳少海 彭曉龍 廣東省醫(yī)療器械質(zhì)量監(jiān)督檢驗所 （廣東 廣州 510663）

錐束CT系統(tǒng)，采用錐形束X射線和平板探測器，圍繞物體進行圓周旋轉(zhuǎn)，采集到的一系列投影數(shù)據(jù)進行計算機圖像重建，得到連續(xù)多層CT圖像，在口腔臨床醫(yī)學成像方面具有廣泛的應用前景和較高的臨床意義[1]。然而，平板探測器上分布有按矩陣形式排列的大量探測元，由于加工工藝精度的限制，無法保證它們的響應是一致的，這就導致重建出的CT圖像含有以圖像中心為圓心的一系列環(huán)形偽影，這些偽影降低了圖像質(zhì)量，限制了后期的圖像處理精度。

CT圖像環(huán)形偽影去除算法主要分為兩類。第一類是前處理方法，即在投影數(shù)據(jù)中利用圖像識別技術(shù)識別出探測器上響應不一致的像素點，然后利用臨近的像素點插值得到該像素點的估計值，最后利用CT圖像重建算法重建出無環(huán)形偽影的CT圖像[2-5]。該類方法的難度在于識別探測器像素點和像素點值插值，若處理不當，重建后CT圖像中的環(huán)形偽影無法全部去除，且部分細節(jié)信息被破壞。第二類是后處理方法，即對含有環(huán)形偽影的CT圖像進行處理[6-8]。一般地，將CT圖像從笛卡爾坐標系轉(zhuǎn)換到極坐標系，環(huán)形偽影從環(huán)形變?yōu)橹本€型，再利用圖像直線檢測技術(shù)和插值/濾波技術(shù)對偽影進行定位和信息恢復，最后變換回笛卡爾坐標系，得到最終圖像。該類方法由于需要坐標系的變換，而且圖像上像素點是離散的，容易造成信息失真。本文旨在提出一種偽影去除效果好且計算復雜度較低的口腔CT圖像環(huán)形偽影去除算法，提出基于中值濾波的口腔CT圖像環(huán)形偽影去除算法，實驗結(jié)果表明，新方法能夠在偽影去除和結(jié)構(gòu)信息保持方面取得滿意的效果。

1.方法

1.1 一維和二維中值濾波算法

在一維信號處理中，一維中值濾波算法常用來對信號中的突變點進行處理[9]。設(shè)一維離散信號為x(i)，則第i個信號點的中值濾波數(shù)學公式表示為：

其中，濾波前的信號記為xold，濾波后的信號記為xnew。i表示信號的第i個信號點。r為整數(shù)，決定以(i)為中心的矩形窗鄰域的范圍，其取值范圍為-R≤r≤R。median(.)函數(shù)表示對xold(i+r)像素點構(gòu)成的集合進行升或降排序，然后取中值作為xnew(i)。

在二維信號處理（即圖像）中，二維中值濾波算法常用來對圖像中椒鹽噪聲進行處理[10]。設(shè)一幅圖像矩陣為x，則像素點(i,j)的中值濾波數(shù)學公式表示為：

其中，濾波前的圖像記為xold，濾波后的圖像記為xnew。(i,j)表示圖像的第i行第j列的像素點。r1和r2為整數(shù)，決定以(i,j)為中心的矩形窗鄰域的范圍，其取值范圍為[-R,R]。median(.)函數(shù)表示對xold(i+r1,j+r2)像素點構(gòu)成的集合進行升或降排序，然后取中值作為xnew(i,j)。圖1所示為二維中值濾波示意圖，窗寬為3×3，對x1～x9的像素值排序后的中間值，即為最終值。

1.2 中值濾波的使用和不足

將CT圖像從笛卡爾坐標系轉(zhuǎn)換到極坐標系，環(huán)形偽影從環(huán)形變?yōu)橹本€型，若濾除這些直線偽影，一維濾波的游走方向垂直于直線偽影，這等價于未做坐標變換前沿著圖像中心向四周進行徑向的一維濾波，避免了因坐標系變換帶來的信息失真問題。實驗發(fā)現(xiàn)這種處理方法有效去除了環(huán)形偽影，但其他的正常結(jié)構(gòu)邊緣容易產(chǎn)生“毛刺”效應。

對CT圖像直接進行二維中值濾波，正常結(jié)構(gòu)邊緣得到比較好的保持，但是環(huán)形偽影不能得到完全消除。

1.3 本文方法

圖1.窗寬為3×3的中值濾波示意圖

圖2.帶有環(huán)形偽影的CT圖像和二維中值濾波結(jié)果

鑒于一維和二維中值濾波去除環(huán)形偽影和正常結(jié)構(gòu)邊緣保持的利弊，融合各自優(yōu)點形成無偽影且正常信息結(jié)構(gòu)得以保持的CT圖像，提出如下方法：

①輸入帶有環(huán)形偽影的口腔CT圖像；

②對①中圖像進行一維中值濾波處理得到圖像A；

③對①中圖像進行二維中值濾波處理得到圖像B；

④用圖像A減去圖像B，得到差值圖像Idiff；

⑤對圖像Idiff進行歸一化處理，使圖像中的每個像素值范圍在0～1：

⑥計算權(quán)重因子：

⑦對圖像A和B加權(quán)求和，得到最終圖像I new：

2.實驗

實驗環(huán)境為CPU：Intel?Core?i3-4160，內(nèi)存：5.0G，操作系統(tǒng)為Windows 10 64bits，編程軟件：Matlab 2009。實驗數(shù)據(jù)為廣州華端科技有限公司的口腔CT系統(tǒng)拍攝的頭部模型CT圖像，仿真產(chǎn)生兩圈明顯的環(huán)形偽影，見圖2a。

圖2b為窗寬為3×3的二維中值濾波處理結(jié)果，可見，圖中多處箭頭處仍含有殘留的環(huán)形偽影；圖2c為窗寬為5×5的二維中值濾波處理結(jié)果，雖然較圖2b有更明顯的改善，但圖中箭頭處的殘留偽影仍然無法去除，采用更大的窗寬同樣無法去除。

圖3a是窗寬為3的一維徑向中值濾波處理結(jié)果，環(huán)形偽影得到完全的去除，但是對于骨組織的邊緣結(jié)構(gòu)處出現(xiàn)了明顯的“毛刺”結(jié)構(gòu)，如圖3a箭頭處，放大圖更顯而易見；圖3b是窗寬為5的一維徑向中值濾波處理結(jié)果，環(huán)形偽影完全消除，但是如圖中箭頭所示的骨組織的邊緣發(fā)生了像素值的改變，表現(xiàn)為骨線不連續(xù)，這是因為在一維徑向中值濾波處理中，誤將骨線當作偽影線進行處理，且窗寬更大其效應越嚴重。圖3c是利用本文方法處理得到的結(jié)果，參數(shù)a=0.4，可見，綜合了圖2c與圖3a的各自優(yōu)勢融合成無偽影且邊緣結(jié)構(gòu)組織得到有效保持的CT圖像。

圖3.一維中值濾波結(jié)果和本文方法處理結(jié)果

3.結(jié)語

本研究在一維和二維中值濾波處理含有環(huán)形偽影的口腔CT圖像基礎(chǔ)上，通過引入權(quán)重函數(shù)模型，結(jié)合不同濾波手段后的圖像的優(yōu)勢，融合成為無環(huán)形偽影且邊緣結(jié)構(gòu)組織得到有效保持的CT圖像。方法通過口腔CT圖像數(shù)據(jù)進行了證實，能夠得到滿意的效果，證明該算法的可行性。

[1]Jaffray DA, Siewerdsen JH.Cone-beam computed tomography with a flat-panel imager: initial performance characterization[J].Medical physics,2000,27(6):1311-1323.

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

[3]余曉鍔,羅君芳,陳武凡.基于弦圖的CT圖像環(huán)形偽影校正[J].醫(yī)學爭鳴,2009,30(3):207-209.

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

[5]王玨,黃蘇紅,蔡玉芳.工業(yè)CT圖像環(huán)形偽影校正[J].光學精密工程,2010,18(5):1226-1233.

[6]齊宏亮,洪虹,徐圓,等.一種快速去除CT環(huán)形偽影的方法[J].南方醫(yī)科大學學報,2012,32(12):1748-1751.

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

[8]Kyriakou Y, Prell D, Kalender WA.Ring artifact correction for high-resolution micro CT[J].Physics in Medicine & Biology,2009,54(17):N385.

[9]Veiga AL, Grunfeld CM.Hardware implementation of a single-cycle one-dimensional median filter[C].Circuits & Systems.IEEE,2015:1-4.

[10]Huang T, Yang G, Tang G.A fast two-dimensional median filtering algorithm[J].IEEE Transactions on Acoustics Speech & Signal Processing,1979,27(1):13-18.