各向異性加權(quán)先驗(yàn)?zāi)Ｐ蚆AP投影域降噪

李安迪，劉 祎，張 權(quán)，桂志國

1.中北大學(xué) 山西省生物醫(yī)學(xué)成像與影像大數(shù)據(jù)重點(diǎn)實(shí)驗(yàn)室，太原 030051

2.中北大學(xué) 信息與通信工程學(xué)院，太原 030051

1 引言

計(jì)算機(jī)斷層掃描技術(shù)（Computed Tomography，CT）的出現(xiàn)極大地促進(jìn)了臨床醫(yī)學(xué)診斷的準(zhǔn)確性，已成為一種有效而重要的成像手段。但常規(guī)劑量的CT輻射對人體有一定的傷害，受檢者需要承擔(dān)未知的健康風(fēng)險(xiǎn)。因此如何在保證CT圖像低噪聲高分辨率的同時(shí)減少射線劑量是CT領(lǐng)域的一個(gè)研究方向。降低X射線劑量的主要措施有降低管電流或管電壓，增大準(zhǔn)直寬度，增大螺距等，最簡單有效的措施為降低管電流。然而由于管電流降低，探測到的光子數(shù)會(huì)隨之減少，投影數(shù)據(jù)會(huì)被大量的量子噪聲污染，使得重建圖像質(zhì)量嚴(yán)重下降，干擾臨床醫(yī)學(xué)診斷的準(zhǔn)確性[1]。目前對這一問題，主要有三種解決途徑：（1）圖像后處理，對重建出的圖像直接去噪；（2）對噪聲特性進(jìn)行建模，用迭代法進(jìn)行重建；（3）投影域處理，對經(jīng)過降噪處理后的投影域數(shù)據(jù)進(jìn)行直接FBP（Filter Back Projection）重建[2]。

圖像后處理的優(yōu)點(diǎn)在于不需要獲得投影數(shù)據(jù)而可直接對重建圖像降噪。在這一方向，Li等人[3]提出了一種改進(jìn)的RSF（Region-Scalable Fitting）模型，將其與穩(wěn)健統(tǒng)計(jì)方法相結(jié)合，獲得了良好的去噪效果，且收斂速度更快；Chen等人[4]提出針對低劑量腹部CT圖像的大尺度加權(quán)強(qiáng)度平均（Weighted Intensity Averaging over Large-Scale Neighborhoods，WIA-LN）算法，利用在大尺度鄰域中具有與目標(biāo)像素相似周圍結(jié)構(gòu)的像素加權(quán)強(qiáng)度平均值來更新目標(biāo)像素，取得了較好的視覺效果；后于2014年將字典學(xué)習(xí)應(yīng)用于低劑量CT圖像域處理中，提出 ASDL（Artifact Suppressed Dictionary Learning）算法，利用偽影的方向和尺度信息訓(xùn)練出偽影原子，再與組織特征原子結(jié)合建立三個(gè)判別字典，可有效去除條形偽影[5]。第二種方法的優(yōu)點(diǎn)在于考慮了投影數(shù)據(jù)的統(tǒng)計(jì)特性。針對此，Wu等人[6]提出一種新的迭代CT重建算法，該算法利用K-稀疏自動(dòng)編碼器（K-Sparse Auto-Encoder，K SAE）從FBP重建出的常規(guī)CT劑量圖像中學(xué)習(xí)一種非線性稀疏先驗(yàn)，基于此先驗(yàn)進(jìn)行迭代重建；受壓縮感知理論和TV最小化的影響，Xu等人[7]在統(tǒng)計(jì)迭代重建框架中，將冗余字典中的系數(shù)約束納入目標(biāo)函數(shù)中，在降低噪聲和保護(hù)細(xì)節(jié)結(jié)構(gòu)特征方面都取得了較好的效果；Makeev等人[8]將基于懲罰似然目標(biāo)函數(shù)和Huber先驗(yàn)的統(tǒng)計(jì)迭代重建算法應(yīng)用于胸部CT中，取得了較好的結(jié)果。但迭代算法計(jì)算量大，成本昂貴。針對第三種方法，Lauzier等人[9]研究了先驗(yàn)圖像約束壓縮感知在CT劑量減少時(shí)的應(yīng)用，得到了較好的去噪效果；Chen等人[1]提出一種應(yīng)用于低劑量CT圖像重建的自適應(yīng)加權(quán)非局部先驗(yàn)?zāi)Ｐ停勺赃m應(yīng)選擇利用圖像全局信息，在分辨率保持和噪聲去除之間取得較好的平衡。Zhang等人[10]提出一種低劑量CT聯(lián)合去噪算法，原始投影數(shù)據(jù)經(jīng)過廣義Anscombe變換（Generalized Anscombe Transform，GAT）后，所含噪聲可視為具有同一方差的加性獨(dú)立高斯白噪聲，用BM3D對其濾波，對所得到的投影數(shù)據(jù)作GAT反變換，F(xiàn)DK重建后再進(jìn)行RNLM濾波處理，所得結(jié)果與原始圖像較為接近；2012年，Zhang[11]等人提出一種基于各向異性加權(quán)先驗(yàn)?zāi)Ｐ偷淖畲蠛篁?yàn)（Maximum A Posteriori，MAP）正弦圖平滑算法，將各向異性擴(kuò)散系數(shù)與傳統(tǒng)的二階PM先驗(yàn)?zāi)Ｐ拖嘟Y(jié)合構(gòu)造出新的正則項(xiàng)，最后用高斯-賽德爾法來求解基于此先驗(yàn)?zāi)Ｐ偷腗AP正弦圖優(yōu)化模型。該算法雖然能依據(jù)不同區(qū)域自適應(yīng)地調(diào)節(jié)平滑程度，有效減少條形偽影，但視覺上偽影依然明顯且邊緣已經(jīng)模糊。

近年來，直覺模糊集在圖像處理中取得了廣泛的應(yīng)用，Atanassov[12]首次提出了直覺模糊集的概念。直覺模糊熵作為直覺模糊集的測度，能夠較好地描述信息的確定性與不確定性。圖像過渡區(qū)域本身存在模糊性，即不確定某一像素是否屬于平坦區(qū)域，給圖像處理帶來不便。對此Chaira[13]提出利用直覺模糊熵檢測圖像邊緣，區(qū)分平滑區(qū)域和細(xì)節(jié)區(qū)域；上官宏[14]定義了正弦圖的局部直覺模糊熵，并將其應(yīng)用于正弦圖平滑算法中。受上述算法啟發(fā)，本文提出一種融合了直覺模糊熵的各向異性加權(quán)先驗(yàn)?zāi)Ｐ停⑵渑cMAP優(yōu)化估計(jì)算法框架結(jié)合，實(shí)現(xiàn)對投影數(shù)據(jù)不同區(qū)域進(jìn)行不同強(qiáng)度的降噪處理。

2 算法原理

2.1 MAP降噪算法模型

MAP是通過最大化后驗(yàn)概率分布函數(shù)來求解目標(biāo)最優(yōu)值。將經(jīng)系統(tǒng)校正和對數(shù)變換后的低劑量投影數(shù)據(jù)看作q，去噪后投影數(shù)據(jù)看作 f，其可建模成一個(gè)馬爾科夫場。根據(jù)Bayesian理論可得后驗(yàn)概率分布為：

由此實(shí)現(xiàn) f的估計(jì)，可建立如下優(yōu)化模型：

近年來確定待估計(jì)參數(shù) f的先驗(yàn)分布P(f)一直是眾多學(xué)者研究的熱點(diǎn)，也是MAP估計(jì)框架的關(guān)鍵點(diǎn)。將 f建模成一個(gè)馬爾科夫隨機(jī)場（Markov Random Field，MRF），根據(jù) Hammersley-Clifford定理給出的MRF與吉布斯隨機(jī)場等效的條件，可得到起正則化作用的先驗(yàn)分布項(xiàng)：

其中，Z為配分函數(shù)，通常設(shè)置為常數(shù)1；U(f)為能量函數(shù)；Uj(f)代表在去噪投影數(shù)據(jù)中像素 j處的能量函數(shù)值；β為一全局超參數(shù)，用來控制先驗(yàn)正則項(xiàng)在優(yōu)化過程中所起作用的大小。

對式（3）兩邊取負(fù)對數(shù)變換得到后驗(yàn)?zāi)芰亢瘮?shù)如下：

文獻(xiàn)[14]研究表明，經(jīng)系統(tǒng)校正和對數(shù)變換后的投影數(shù)據(jù)q近似服從非平穩(wěn)高斯分布，該投影數(shù)據(jù)方差為：

其中，rj、fj為與掃描系統(tǒng)相關(guān)的參數(shù)。采用統(tǒng)計(jì)迭代法來估計(jì)去噪后的投影數(shù)據(jù)，有如下似然函數(shù)形式：

可得到式（4）中負(fù)對數(shù)似然函數(shù)為：

其中，Σ是以σj為組成元素的對角矩陣。最大后驗(yàn)概率模型可轉(zhuǎn)化為求解如下函數(shù)，即通過最小化后驗(yàn)?zāi)芰亢瘮?shù)來實(shí)現(xiàn)投影數(shù)據(jù)的降噪：

2.2 先驗(yàn)?zāi)Ｐ?/h3>

一般情況下，能量函數(shù)通過鄰域中像素間差異的加權(quán)和求得：

其中，Nj表示中心像素為 j的鄰域；權(quán)重ωij表示鄰域像素i與中心像素 j間的相關(guān)關(guān)系。傳統(tǒng)的二次方先驗(yàn)?zāi)Ｐ?u(t)=t2/2)對各個(gè)像素?zé)o差別對待，由此會(huì)造成重建圖像邊緣模糊和平坦區(qū)域過平滑等問題。Huber先驗(yàn)?zāi)Ｐ涂筛鶕?jù)圖像中心像素與鄰域像素的灰度差值進(jìn)行自適應(yīng)性平滑，因此本文選用Huber先驗(yàn)?zāi)Ｐ?，其形式如下?/p>

Huber先驗(yàn)?zāi)Ｐ屠瞄撝礹th來控制圖像的平滑程度，當(dāng)像素差值大于hth時(shí)，該區(qū)域?qū)儆谶吘壖?xì)節(jié)區(qū)域的可能性更大，采用平滑程度較輕的線性形式；當(dāng)像素差值小于hth時(shí)，該區(qū)域?qū)儆谄教箙^(qū)域的可能性更大，采用二次項(xiàng)形式加大平滑程度。

權(quán)重ωij表示中心像素與鄰域內(nèi)其他像素間的相關(guān)關(guān)系，通常情況下定義為兩個(gè)像素間歐式距離的反比。考慮到二維正弦圖中探測器探元方向和角度方向上對鄰域內(nèi)像素的影響，且探元方向影響大于角度方向，水平方向權(quán)重取1.00，垂直方向權(quán)重取0.25，如圖1所示。

圖1 鄰域內(nèi)權(quán)重值

2.3 改進(jìn)的各向異性加權(quán)先驗(yàn)?zāi)Ｐ?/h3>

傳統(tǒng)二次先驗(yàn)?zāi)Ｐ筒荒苡行У貐^(qū)分邊緣細(xì)節(jié)區(qū)域和平坦區(qū)域，不能實(shí)現(xiàn)對不同區(qū)域的自適應(yīng)平滑，重建結(jié)果中易導(dǎo)致邊緣細(xì)節(jié)區(qū)域過平滑，平坦區(qū)域平滑力度不足等情況，同時(shí)考慮到權(quán)重ωij是固定不變的，提出一種融合了各向異性擴(kuò)散系數(shù)的Huber先驗(yàn)?zāi)Ｐ停?/p>

常用的各向異性擴(kuò)散系數(shù)有如下形式：

正弦圖中平坦區(qū)域和邊緣細(xì)節(jié)區(qū)域區(qū)分不明顯，即本身具有模糊性。根據(jù)直覺模糊集理論，將正弦圖看作直覺模糊集合，某一像素j對該直覺模糊集合的隸屬度φ(j,γ)、非隸屬度 η(j,γ)和猶豫度 ρ(j,γ)定義如下：

構(gòu)造出正弦圖的局部直覺模糊熵：

K衡量像素點(diǎn) j在正弦圖中處于平坦區(qū)域還是細(xì)節(jié)區(qū)域，K值越大，像素點(diǎn) j處于細(xì)節(jié)或邊緣區(qū)域的確定性越大；反之K值越小，像素點(diǎn) j處于平坦區(qū)域的確定性越大。

綜合考慮正弦圖的梯度模、局部方差與局部直覺模糊熵，可構(gòu)造一種新的各向異性擴(kuò)散系數(shù)：

根據(jù)以上各向異性擴(kuò)散系數(shù)的擴(kuò)散特點(diǎn)，此先驗(yàn)?zāi)Ｐ涂蓪φ覉D進(jìn)行自適應(yīng)平滑。若像素點(diǎn) j處于平坦區(qū)域，其對應(yīng)的K值、梯度模均較小值變大，平滑程度加大；若像素點(diǎn) j處于邊緣區(qū)域，其對應(yīng)的K值、梯度模均較大，值變小，平滑程度減弱。因此可有效地在去除噪聲的同時(shí)，較好地保持圖像邊緣細(xì)節(jié)信息。

綜上所述，本文提出的新的基于各向異性擴(kuò)散系數(shù)的先驗(yàn)?zāi)Ｐ蜑椋?/p>

根據(jù)本文提出的各向異性加權(quán)Huber先驗(yàn)?zāi)Ｐ?，式?）所描述的算法優(yōu)化模型為：

3 實(shí)驗(yàn)結(jié)果與分析

本文采用三種體模來驗(yàn)證基于改進(jìn)的各向異性加權(quán)先驗(yàn)?zāi)Ｐ偷腗AP投影域降噪算法的有效性：體模1為數(shù)字骨盆模型，體模2為Shepp-Logan頭模型，體模3為數(shù)字胸腔模型（thorax phantom）。采用三種對比算法：經(jīng)典FBP直接重建算法，Wang等人[15]提出的懲罰重加權(quán)最小二乘優(yōu)化算法PRWLS，Zhang等人[11]提出的各向異性加權(quán)先驗(yàn)MAP平滑算法。各模型原始圖像如圖2所示，體模大小均為512×512。實(shí)驗(yàn)所用編程工具為Matlab，版本為R2016a。計(jì)算機(jī)操作系統(tǒng)為Windows 10 64位操作系統(tǒng)，處理器為Intel?CoreTMi7-7700 CPU@3.60 GHz。

圖2 各模型的原始圖像

實(shí)驗(yàn)中所用模擬投影數(shù)據(jù)通過等角扇束掃描方式得到，探測器探元弧形排列為888個(gè)，頂點(diǎn)距離X射線源949 mm，距離旋轉(zhuǎn)中心408 mm，984個(gè)投影角度均勻分布在360°圓周上，因此正弦圖大小為984×888。實(shí)驗(yàn)過程中所求正弦圖方差相關(guān)參數(shù)為rj=150,η=42 000,迭代次數(shù)設(shè)置為15次。超全局參數(shù)β控制正則項(xiàng)在平滑過程中的作用，本文所提算法對Huber先驗(yàn)不同表達(dá)式設(shè)置不同的β值，當(dāng)像素差值大于閾值時(shí)，增大β值以加大先驗(yàn)正則項(xiàng)的作用，控制了平滑程度，反之像素差值小于閾值時(shí)，β值較小，先驗(yàn)項(xiàng)在對投影數(shù)據(jù)處理過程中作用較小，平滑程度加大。通過手動(dòng)調(diào)節(jié)其值以得到視覺效果最優(yōu)的重建結(jié)果圖，各算法β值設(shè)置如表1所示。

表1 參數(shù)設(shè)置

將FBP直接重建算法、PRWLS算法、文獻(xiàn)[12]算法和本文算法分別應(yīng)用于體模1，調(diào)整參數(shù)得到各個(gè)算法對應(yīng)的最優(yōu)效果，如圖3所示。從結(jié)果圖中可看出，直接使用FBP算法得到的重建圖像被條形偽影覆蓋，視覺效果最差；經(jīng)PRWLS算法重建出的圖像視覺效果得到顯著提升，但在中心部位條形偽影依然嚴(yán)重；文獻(xiàn)[12]算法一定程度上減少了中心條形偽影；本文算法在有效減少中心條形偽影的同時(shí)，更好地保持了邊緣與細(xì)節(jié)信息，在左右兩側(cè)嵌套的黑白兩色圓形處邊緣保持更明顯。

圖3 體模1實(shí)驗(yàn)結(jié)果圖

同樣的方法應(yīng)用于體模2和體模3，圖4（a1）～（e1）和圖5（a1）～（e1）分別為體模2和體模3用四種方法重建得到的結(jié)果，圖4（a2）～（e2）和圖5（a2）～（e2）為其對應(yīng)的中心放大圖。從體模2的中心放大區(qū)域可看出三個(gè)白色橢圓區(qū)域噪聲明顯減少，且兩個(gè)黑色區(qū)域部分更平滑；體模3中整體背景部分更平滑，放大部分中可看出在白色橢圓部分尤其是中間橢圓，前三種方法處理后還有可見的黑色噪聲點(diǎn)，而采用本文算法則有效去除了噪聲點(diǎn)。由此可證明，文獻(xiàn)[12]算法在去除噪聲方面有明顯效果，但本文算法不僅能夠有效地去除噪聲，且在保持細(xì)節(jié)與邊緣信息方面有更好的效果。以上實(shí)驗(yàn)充分說明了了本文算法的有效性。

圖4 體模2實(shí)驗(yàn)結(jié)果圖

圖5 體模3實(shí)驗(yàn)結(jié)果圖

經(jīng)過仿真實(shí)驗(yàn)驗(yàn)證，本文算法在去除噪聲和保護(hù)邊緣信息方面都有提升。這是由于相對于其他算法，本文算法在對投影數(shù)據(jù)的處理中結(jié)合了局部直覺模糊熵，采用Huber先驗(yàn)?zāi)Ｐ?，更有效地區(qū)分了平坦區(qū)域和邊緣細(xì)節(jié)區(qū)域，使在平坦區(qū)域平滑力度加大以有效去除噪聲，在邊緣細(xì)節(jié)區(qū)域減小平滑力度，保留相關(guān)信息。為了進(jìn)一步說明算法有效性，選取體模中像素值變化較多的行或列進(jìn)行比較，圖6～圖8分別給出了體模1在第245行的像素值，體模2在第220列的像素值，體模3在第300行的像素值相對于四種算法的對比。

圖6 體模1第245行像素值

由像素值對比圖可看出本文提出的基于改進(jìn)的各向異性擴(kuò)散Huber先驗(yàn)平滑算法所重建出的圖像在邊緣和平滑區(qū)域與原始干凈無噪圖像像素值更加接近。為了進(jìn)一步說明這一問題，圖9～圖11選取像素值變化平緩區(qū)域進(jìn)行像素值比較?？煽闯霰疚乃惴ㄖ亟ǔ龅膱D像整體尤其是在中心區(qū)域更接近原始圖像。

圖7 體模2第220列像素值

圖8 體模3第300行像素值

圖9 體模1第310行像素值

圖10 體模2第100列像素值

圖11 體模3第80列像素值

為了客觀說明本文算法的有效性，采用信噪比定量描述各種算法，同時(shí)記錄了各算法運(yùn)行時(shí)間。計(jì)算信噪比所用的公式如下：

表2 各算法客觀評價(jià)指標(biāo)

其中，y代表經(jīng)過算法處理后的圖像像素值。表2給出了各種算法的客觀評價(jià)指標(biāo)。

由表2可以看出本文算法在信噪比上有明顯提升，重建圖像質(zhì)量優(yōu)于其他算法，且計(jì)算速度與文獻(xiàn)[12]算法相當(dāng)，從客觀上說明了本文算法的有效性。

4 結(jié)束語

本文提出了一種改進(jìn)的基于各向異性加權(quán)先驗(yàn)?zāi)Ｐ偷腗AP投影域降噪算法。針對原始算法中使用的二次方先驗(yàn)?zāi)Ｐ蛯Ω鱾€(gè)像素?zé)o差別對待問題，本文采用Huber先驗(yàn)?zāi)Ｐ?，針對不同的區(qū)域進(jìn)行不同程度的平滑處理。同時(shí)考慮到局部直覺模糊熵可有效地區(qū)分正弦圖中的平滑區(qū)域和細(xì)節(jié)邊緣區(qū)域，將其與傳統(tǒng)的各向異性擴(kuò)散系數(shù)相結(jié)合，并采用正弦圖局部方差來實(shí)現(xiàn)該系數(shù)的自適應(yīng)性調(diào)節(jié)。最后將其融合于最大后驗(yàn)概率估計(jì)模型中，實(shí)現(xiàn)對整個(gè)正弦圖的平滑處理。實(shí)驗(yàn)證明，本文算法相對于其他三種算法能夠有效去除圖像噪聲，同時(shí)可以有效地保持圖像邊緣與細(xì)節(jié)信息，處理得到的圖像與原始圖像更接近。但同時(shí)算法中一些參數(shù)也需要根據(jù)經(jīng)驗(yàn)手動(dòng)調(diào)整，可以作為下一步研究的方向。