各向異性加權先驗模型MAP投影域降噪

李安迪，劉 祎，張 權，桂志國

1.中北大學 山西省生物醫學成像與影像大數據重點實驗室，太原 030051

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

1 引言

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

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

近年來，直覺模糊集在圖像處理中取得了廣泛的應用，Atanassov[12]首次提出了直覺模糊集的概念。直覺模糊熵作為直覺模糊集的測度，能夠較好地描述信息的確定性與不確定性。圖像過渡區域本身存在模糊性，即不確定某一像素是否屬于平坦區域，給圖像處理帶來不便。對此Chaira[13]提出利用直覺模糊熵檢測圖像邊緣，區分平滑區域和細節區域；上官宏[14]定義了正弦圖的局部直覺模糊熵，并將其應用于正弦圖平滑算法中。受上述算法啟發，本文提出一種融合了直覺模糊熵的各向異性加權先驗模型，并將其與MAP優化估計算法框架結合，實現對投影數據不同區域進行不同強度的降噪處理。

2 算法原理

2.1 MAP降噪算法模型

MAP是通過最大化后驗概率分布函數來求解目標最優值。將經系統校正和對數變換后的低劑量投影數據看作q，去噪后投影數據看作 f，其可建模成一個馬爾科夫場。根據Bayesian理論可得后驗概率分布為：

由此實現 f的估計，可建立如下優化模型：

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

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

對式（3）兩邊取負對數變換得到后驗能量函數如下：

文獻[14]研究表明，經系統校正和對數變換后的投影數據q近似服從非平穩高斯分布，該投影數據方差為：

其中，rj、fj為與掃描系統相關的參數。采用統計迭代法來估計去噪后的投影數據，有如下似然函數形式：

可得到式（4）中負對數似然函數為：

其中，Σ是以σj為組成元素的對角矩陣。最大后驗概率模型可轉化為求解如下函數，即通過最小化后驗能量函數來實現投影數據的降噪：

2.2 先驗模型

一般情況下，能量函數通過鄰域中像素間差異的加權和求得：

其中，Nj表示中心像素為 j的鄰域；權重ωij表示鄰域像素i與中心像素 j間的相關關系。傳統的二次方先驗模型(u(t)=t2/2)對各個像素無差別對待，由此會造成重建圖像邊緣模糊和平坦區域過平滑等問題。Huber先驗模型可根據圖像中心像素與鄰域像素的灰度差值進行自適應性平滑，因此本文選用Huber先驗模型，其形式如下：

Huber先驗模型利用閾值hth來控制圖像的平滑程度，當像素差值大于hth時，該區域屬于邊緣細節區域的可能性更大，采用平滑程度較輕的線性形式；當像素差值小于hth時，該區域屬于平坦區域的可能性更大，采用二次項形式加大平滑程度。

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

圖1 鄰域內權重值

2.3 改進的各向異性加權先驗模型

傳統二次先驗模型不能有效地區分邊緣細節區域和平坦區域，不能實現對不同區域的自適應平滑，重建結果中易導致邊緣細節區域過平滑，平坦區域平滑力度不足等情況，同時考慮到權重ωij是固定不變的，提出一種融合了各向異性擴散系數的Huber先驗模型：

常用的各向異性擴散系數有如下形式：

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

構造出正弦圖的局部直覺模糊熵：

K衡量像素點 j在正弦圖中處于平坦區域還是細節區域，K值越大，像素點 j處于細節或邊緣區域的確定性越大；反之K值越小，像素點 j處于平坦區域的確定性越大。

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

根據以上各向異性擴散系數的擴散特點，此先驗模型可對正弦圖進行自適應平滑。若像素點 j處于平坦區域，其對應的K值、梯度模均較小值變大，平滑程度加大；若像素點 j處于邊緣區域，其對應的K值、梯度模均較大，值變小，平滑程度減弱。因此可有效地在去除噪聲的同時，較好地保持圖像邊緣細節信息。

綜上所述，本文提出的新的基于各向異性擴散系數的先驗模型為：

根據本文提出的各向異性加權Huber先驗模型，式（8）所描述的算法優化模型為：

3 實驗結果與分析

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

圖2 各模型的原始圖像

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

表1 參數設置

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

圖3 體模1實驗結果圖

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

圖4 體模2實驗結果圖

圖5 體模3實驗結果圖

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

圖6 體模1第245行像素值

由像素值對比圖可看出本文提出的基于改進的各向異性擴散Huber先驗平滑算法所重建出的圖像在邊緣和平滑區域與原始干凈無噪圖像像素值更加接近。為了進一步說明這一問題，圖9～圖11選取像素值變化平緩區域進行像素值比較。可看出本文算法重建出的圖像整體尤其是在中心區域更接近原始圖像。

圖7 體模2第220列像素值

圖8 體模3第300行像素值

圖9 體模1第310行像素值

圖10 體模2第100列像素值

圖11 體模3第80列像素值

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

表2 各算法客觀評價指標

其中，y代表經過算法處理后的圖像像素值。表2給出了各種算法的客觀評價指標。

由表2可以看出本文算法在信噪比上有明顯提升，重建圖像質量優于其他算法，且計算速度與文獻[12]算法相當，從客觀上說明了本文算法的有效性。

4 結束語

本文提出了一種改進的基于各向異性加權先驗模型的MAP投影域降噪算法。針對原始算法中使用的二次方先驗模型對各個像素無差別對待問題，本文采用Huber先驗模型，針對不同的區域進行不同程度的平滑處理。同時考慮到局部直覺模糊熵可有效地區分正弦圖中的平滑區域和細節邊緣區域，將其與傳統的各向異性擴散系數相結合，并采用正弦圖局部方差來實現該系數的自適應性調節。最后將其融合于最大后驗概率估計模型中，實現對整個正弦圖的平滑處理。實驗證明，本文算法相對于其他三種算法能夠有效去除圖像噪聲，同時可以有效地保持圖像邊緣與細節信息，處理得到的圖像與原始圖像更接近。但同時算法中一些參數也需要根據經驗手動調整，可以作為下一步研究的方向。