基于自動抑制的傳感器欠采樣CT 影像降噪方法

高 華，傅德勝，李瀟陽

(1.南京信息工程大學(xué) 計算機與軟件學(xué)院，江蘇 南京210044;2.南京信息工程大學(xué) 環(huán)境科學(xué)與工程學(xué)院，江蘇 南京210044)

0 引 言

現(xiàn)今，CT 掃描已被作為一種常見診療手段而得到廣泛應(yīng)用，但相比其他X 射線檢查它既有優(yōu)點也有明顯缺點。其中，CT 的大量輻射問題是最為突出的，因此，國內(nèi)外研究者都在研究如何有效降低CT 的掃描輻射劑量。實際使用中，醫(yī)療工作者主要通過降低X 射線管的電壓或者管電流來達(dá)到降低輻射劑量的目的，但這會使重建后的影像具有較低的信噪比，嚴(yán)重影響了診斷。

目前研究者已經(jīng)提出了很多方法有效降低患者的受輻射劑量，主要有三方面研究重點:加快硬件掃描速度、采用感興趣區(qū)域掃描和稀疏重建、低劑量CT 掃描降噪重建。

第一類加快硬件掃描速度是各大廠商研究的重點，例如:目前推出的雙源螺旋多層掃描等［1］，其主要方式是提高單位時間內(nèi)掃描的速度，間接縮短時間從而降低輻射劑量。第二類是目前剛興起的研究領(lǐng)域，它主要是通過稀疏數(shù)據(jù)進行數(shù)學(xué)重建。因為采用了局部感興趣和稀疏采樣降低了輻射劑量［2］，典型算法有基于CS 壓縮感知［3］的稀疏重建和欠采樣的ROI 重建［4］。第三類是使用降低掃描電流或管電壓的方式獲得的低輻射投影數(shù)據(jù)［5］，并在此基礎(chǔ)上采用一些算法降噪重建，典型算法如:PWLS［6］、投影域像素降噪［7］、TV 正則降噪等算法［8］。

這些方法都取得了較高的信噪比，但是目前第一類方法進展緩慢，第二類算法仍處起步階段，而第三類雖研究較早但實際應(yīng)用不多。這主要是因為第三類算法運算開銷較大，且有部分需要經(jīng)驗設(shè)置。

本文提出了一種基于投影域數(shù)據(jù)分析抑制的改進的自適應(yīng)算法，目的是糾正由傳感器老化和欠采樣等造成的嚴(yán)重噪聲，在其抑制過程中自動修正抑制參數(shù)并提高其抑制效果，實驗結(jié)果表明:該算法抑制噪聲取得了較好的效果。

1 噪聲分析與模型

CT 噪聲主要由2 個方面造成:1)機械產(chǎn)生的噪聲，包括了震動、傳感器老化等，這一類是無法避免的;2)人工修改機器參數(shù)，導(dǎo)致傳感器發(fā)生光子饑餓。早先國內(nèi)學(xué)者盧宏冰等人［9］采用固定角度多次掃描數(shù)據(jù)，并分析其統(tǒng)計特性，發(fā)現(xiàn)其校正和對數(shù)變換后的數(shù)據(jù)近似服從加性噪聲式(1)污染且滿足式(2)的總體非線性非平穩(wěn)高斯分布

其中，f 為理想數(shù)據(jù)，n 為非平穩(wěn)高斯噪聲。

但實際投影數(shù)據(jù)中經(jīng)常會有一些數(shù)值較高的點，在文獻(xiàn)［10］中張元科博士對其進行論述，并且給出了異值點集中區(qū)域。但該方法是基于一種肯定假設(shè):即大部分?jǐn)?shù)據(jù)是高于并限定在一定的數(shù)值范圍內(nèi)，并且假設(shè)異值點較為集中在這些區(qū)域。圖1 給出了示例。

圖1 異值點集中區(qū)域圖形Fig 1 Different values in concentrated area

然而實際數(shù)據(jù)中也發(fā)現(xiàn)，在一些低角度的數(shù)據(jù)中也會出現(xiàn)異值，分析原因是:X 射線的路徑并不完全是直線，即使是在穿透距離較短的情況下仍然會有大量光子跑到其他傳感器上，因此，就造成了即使不在異值重點區(qū)域也會出現(xiàn)異值的特點;同時也會因為傳感器老化而造成光子接受數(shù)值有較大偏差，如圖2 所示。

因此，在去除明顯異噪集中區(qū)域的噪聲外，也需要對非集中區(qū)的異值點進行抑制。而異值點通常是明顯高于周邊傳感器的數(shù)值的，而且其在單一方向上的震動幅度是滿足一定范圍的，因此，這就提供了一個較好的判斷依據(jù)。

2 算法描述

2.1 算法分析

圖2 短穿透投影距離下存在異值Fig 2 Different value exists in short penetrate distance(zoomed area，arrow point to)

為方便檢測到異值點，設(shè)投影圖像為p(x，y)，則g(x，y)為以(x，y)為中心的n×n 窗口，gm(x，y)為窗口均值，則d=g－gm，w(x，y)為窗口內(nèi)8 領(lǐng)域，t(x，y)為該角度上的投影數(shù)據(jù)的左右各n+1 方向。示意圖如圖3 所示。

圖3 濾波窗口示例Fig 3 Example of filtering window

選擇T1=2σgm，T2=gm，T3=tm，T4=σtm，其中，tm為某角度上的以(x，y)為中心的投影數(shù)據(jù)條內(nèi)均值，σgm為窗口內(nèi)的方差，根據(jù)統(tǒng)計學(xué)原理，選擇2σgm為95%置信區(qū)間。由于可能存在期望值相近，為此，再引入標(biāo)準(zhǔn)離差

比較σ1，σ2，σ2，若σ1較小，選擇方形窗口內(nèi)的中值濾波;相反，則選條狀區(qū)域的均值。

在確定自適應(yīng)濾波方式后，如何確定異值點，在文獻(xiàn)［10］中給出了方法，該方法是一種比較可靠的方法，它為感興趣的區(qū)域設(shè)置閾值，按照其規(guī)則確認(rèn)異值點，但該方法需要人工設(shè)置參數(shù)，且在一些角度上數(shù)據(jù)全部偏大時進行了過度濾值，并且它是以8 領(lǐng)域方形窗口為主要濾波窗口，這忽略了在單一掃描方向上的數(shù)值非平穩(wěn)連續(xù)特性。

本文采用了一種改進的運作方式來判斷和更新像素中的異值點，首先相同的就是選擇

其中，σ'為兩標(biāo)準(zhǔn)離差較小的那個所對應(yīng)的窗口方差。

然后選擇更新圖像中的點，不過這里選用了2 個窗口中的數(shù)值，其主要目的是進行保守的異值點更新，其更新公式為

在確定異值點后，并更新圖像后單角度上的數(shù)據(jù)更加趨于緩和。處理后結(jié)果如圖4、圖5 所示。

圖4 整體抑制Fig 4 Global inhibition

圖5 短穿透抑制Fig 5 Short penetration suppression

2.2 算法總結(jié)

1)設(shè)置窗口大小為3，標(biāo)準(zhǔn)離差和窗口方差、投影方向均值均為0;

4)重復(fù)直至全部掃描完成;

5)更新g(x，y)=g'(x，y);

6)其他經(jīng)典重建算法重建。

3 實驗材料與仿真結(jié)果

實驗材料選用256×256 的Sheep－Logan 頭部，其圖形如圖6 所示。為了模擬CT 掃描使用扇形束投影fanbeam。旋轉(zhuǎn)角度是180°，984 個采樣方向，探測器一共有888 個，探測器距離是1.029 3，X 射線到旋轉(zhuǎn)中心距離是541，中心到探測器的距離是408。在這些參數(shù)下仿真出CT 投影數(shù)據(jù)。

然后利用公式(2)和圖6 的數(shù)據(jù)進行加性噪聲的添加，之后，再隨機添加異值點，獲得低劑量情況下的模擬數(shù)據(jù)，其中，β 取22 000。

下面給出傳統(tǒng)FBP 直接重建［11］、PWLS、以及本文給定的方法濾波后再進行PWLS 方法得到的重建。最后比對重建后與原始圖差值，其結(jié)果如圖7(f)所示。

可以看出，本文使用的方法較好地保存了邊緣信息，清晰度比PWLS 要好，但是也發(fā)現(xiàn)在下方的3 個連續(xù)的橢圓上邊緣有模糊，分析結(jié)果證明:由于本文算法在緩和投影數(shù)據(jù)間的差距時，也將一些特別小的物體的投影數(shù)據(jù)當(dāng)作是真實數(shù)據(jù)造成。

圖6 Sheep-Logan 模型Fig 6 Sheep-Logan model

圖7 三種方法重建圖(a)～(c)，比較差值圖(d)～(f)Fig 7 Reconstruction images(show in a～c)，difference value comparison show in d～f use 3 methods

從圖7(d)～(f)3 個差值圖像可以看出直接FBP 效果(圖(d))最差，在其邊緣、內(nèi)部平滑區(qū)域都存在嚴(yán)重的噪點，PWLS(圖(e))與原始圖像差距減小，但也發(fā)現(xiàn)其在邊緣和在窄邊緣發(fā)生的復(fù)制現(xiàn)象，相比較而言，PWLS 抑制了一部分的內(nèi)部平滑區(qū)域噪聲，但在窄地區(qū)抑制效果一般。而本文方法可以明顯看出，其與原圖在內(nèi)部平滑區(qū)域相差不大，但由于其采用直接對投影域數(shù)據(jù)處理的方法，并沒有在后續(xù)重建中對邊緣進行處理，因此，可以看到邊緣的差異。

作為主觀判斷可以得到一定的結(jié)果，但為了更好地表示數(shù)據(jù)分析的特征，這里引入客觀的信噪比分析工具信—噪功率比(PSNR)，其計算公式

式中 xij為重建后圖像，fij為理想圖像，M，N 為圖像大小，i=1，2，3，…，M;j=1，2，3，…，N;M，N∈Z+，它反映理想圖像與處理后圖像的信噪比情況。計算結(jié)果如表1 所示。

表1 仿真重建后的數(shù)據(jù)計算PSNRTab 1 PSNR computed by datas after simulation reconstruction

從計算結(jié)果與主觀判斷看，它提高了圖像信號所占比例，且圖像的質(zhì)量是可以接受的。

4 結(jié)束語

本文分析CT 投影數(shù)據(jù)的特點，介紹了一些研究成果，提出其不足。分析了投影數(shù)據(jù)在單角度上的數(shù)學(xué)特性，發(fā)現(xiàn)投影數(shù)據(jù)去異值點后服從一定的非平穩(wěn)高斯噪聲分布。相鄰角度數(shù)據(jù)相關(guān)，且低穿透區(qū)域也存在異值，解釋了產(chǎn)生的可能原因。在波動較大時，同一投影角度上震動幅度大體相當(dāng)，并針對該特性提出算法改進，最后進行了仿真實驗。實驗結(jié)果表明:濾除噪聲后的數(shù)據(jù)趨于非平穩(wěn)高斯分布，重建效果較好，克服了傳感器老化等因素，且自動化程度提高，免除了人工設(shè)定異值范圍的麻煩。

［1］ 張宗軍，盧光明.雙源CT 及其臨床應(yīng)用［J］.醫(yī)學(xué)研究生學(xué)報，2007，20(4):416－418.

［2］ 劉圓圓，張 麗，陳志強.改進的單程分裂合并分割方法在雙能CT 欠采樣方法中的應(yīng)用［C］∥全國射線數(shù)字成像與CT新技術(shù)研討會論文集，上海:2009.

［3］ 練秋生，郝鵬鵬.基于壓縮傳感和代數(shù)重建法的CT 圖像重建［J］.光學(xué)技術(shù)，2009(3):422－425.

［4］ Xu Q，Mou X Q，Wang G，et al.Statistical interior tomography［J］.IEEE Transactions on Medical Imaging，2011，30(5):1116－1128.

［5］ 朱天照，唐光健，蔣學(xué)祥.低劑量螺旋CT 肺結(jié)節(jié)掃描與常規(guī)劑量CT 的對照研究［J］.中華放射學(xué)雜志，2004，38(4):428－431.

［6］ Wang Q X，Li H，Lam K Y.Development of a new meshless－point weighted least－squares(PWLS)method for computational mechanics［J］.Computational Mechanics，2005，35(3):170－181.

［7］ 王 旭，陳志強，熊 華，等.聯(lián)合代數(shù)重建算法中基于像素的投影計算方法［J］.核電子學(xué)與探測技術(shù)，2006，25(6):785－788.

［8］ 張瀚銘，閆 鑌，李 磊.一種基于TV 正則化的有限角度CT圖像重建算法［C］∥綿陽:全國射線數(shù)字成像與CT 新技術(shù)研討會論文集，2012.

［9］ 盧虹冰，李 響，蕭穎聰，等.利用K－L 域懲罰加權(quán)最小均方平滑對低劑量CT 投影數(shù)據(jù)進行噪聲濾除［J］.中國醫(yī)學(xué)物理學(xué)雜志，2002，19(4):200－204.

［10］張元科，張軍英，盧虹冰.低劑量CT 投影圖像噪聲分析及去噪算法研究［J］.光電子.激光，2010，21(7):1073－1078.

［11］錢姍姍，黃 靜，馬建華，等.基于投影數(shù)據(jù)非單調(diào)性全變分恢復(fù)的低劑量CT 重建［J］.電子學(xué)報，2011，39(7):1702－1707.