高 華，傅德勝，李瀟陽(yáng)

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

0 引 言

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

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

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

這些方法都取得了較高的信噪比，但是目前第一類方法進(jìn)展緩慢，第二類算法仍處起步階段，而第三類雖研究較早但實(shí)際應(yīng)用不多。這主要是因?yàn)榈谌愃惴ㄟ\(yùn)算開(kāi)銷較大，且有部分需要經(jīng)驗(yàn)設(shè)置。

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

1 噪聲分析與模型

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

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

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

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

然而實(shí)際數(shù)據(jù)中也發(fā)現(xiàn)，在一些低角度的數(shù)據(jù)中也會(huì)出現(xiàn)異值，分析原因是:X 射線的路徑并不完全是直線，即使是在穿透距離較短的情況下仍然會(huì)有大量光子跑到其他傳感器上，因此，就造成了即使不在異值重點(diǎn)區(qū)域也會(huì)出現(xiàn)異值的特點(diǎn);同時(shí)也會(huì)因?yàn)閭鞲衅骼匣斐晒庾咏邮軘?shù)值有較大偏差，如圖2 所示。

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

2 算法描述

2.1 算法分析

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

為方便檢測(cè)到異值點(diǎn)，設(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)計(jì)學(xué)原理，選擇2σgm為95%置信區(qū)間。由于可能存在期望值相近，為此，再引入標(biāo)準(zhǔn)離差

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

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

本文采用了一種改進(jìn)的運(yùn)作方式來(lái)判斷和更新像素中的異值點(diǎn)，首先相同的就是選擇

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

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

在確定異值點(diǎn)后，并更新圖像后單角度上的數(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 實(shí)驗(yàn)材料與仿真結(jié)果

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

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

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

可以看出，本文使用的方法較好地保存了邊緣信息，清晰度比PWLS 要好，但是也發(fā)現(xiàn)在下方的3 個(gè)連續(xù)的橢圓上邊緣有模糊，分析結(jié)果證明:由于本文算法在緩和投影數(shù)據(jù)間的差距時(shí)，也將一些特別小的物體的投影數(shù)據(jù)當(dāng)作是真實(shí)數(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 個(gè)差值圖像可以看出直接FBP 效果(圖(d))最差，在其邊緣、內(nèi)部平滑區(qū)域都存在嚴(yán)重的噪點(diǎn)，PWLS(圖(e))與原始圖像差距減小，但也發(fā)現(xiàn)其在邊緣和在窄邊緣發(fā)生的復(fù)制現(xiàn)象，相比較而言，PWLS 抑制了一部分的內(nèi)部平滑區(qū)域噪聲，但在窄地區(qū)抑制效果一般。而本文方法可以明顯看出，其與原圖在內(nèi)部平滑區(qū)域相差不大，但由于其采用直接對(duì)投影域數(shù)據(jù)處理的方法，并沒(méi)有在后續(xù)重建中對(duì)邊緣進(jìn)行處理，因此，可以看到邊緣的差異。

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

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

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

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

4 結(jié)束語(yǔ)

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

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

［2］ 劉圓圓，張 麗，陳志強(qiáng).改進(jìn)的單程分裂合并分割方法在雙能CT 欠采樣方法中的應(yīng)用［C］∥全國(guó)射線數(shù)字成像與CT新技術(shù)研討會(huì)論文集，上海: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 的對(duì)照研究［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］ 王 旭，陳志強(qiáng)，熊 華，等.聯(lián)合代數(shù)重建算法中基于像素的投影計(jì)算方法［J］.核電子學(xué)與探測(cè)技術(shù)，2006，25(6):785－788.

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

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

［10］張?jiān)疲瑥堒娪?，盧虹冰.低劑量CT 投影圖像噪聲分析及去噪算法研究［J］.光電子.激光，2010，21(7):1073－1078.

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