基于小波和非局部的全變差中值先驗重建算法

張 芳，張 權，崔學英，董嬋嬋，劉 祎，孫未雅，白云蛟，桂志國，2＋

（1.中北大學 電子測試技術國家重點實驗室，山西 太原030051；2.中北大學 儀器科學與動態測試教育部重點實驗室，山西 太原030051）

0 引 言

低劑量CT［1］投影數據噪聲模型的數據特點已經被廣泛研究，其中Wang等通過對多個體模進行反復實驗和分析，最終得出低劑量CT 投影數據經過對數變換后的均值和方差之間呈非線性遞增關系，近似服從非平穩高斯分布的結論。針對投影數據統計特性的分析，由于最大似然期望最大算法 （MLEM）在重建過程中考慮了觀測數據的統計特性，且具有非負性、全局收斂性和計數保持的特點，重建圖像的效果得到了較好的改善，故而被廣泛的使用。然而在實際中，當迭代次數達到一定后，隨著迭代次數的增多重建圖像的質量會出現棋盤效應，從而造成圖像失真。基于Bayesian理論的最大后驗 （maximum a posterior，MAP）法可以有效解決上述問題［2］，該方法既考慮了低劑量CT投影數據的統計特性，又通過對先驗分布加入先驗信息，使進行多次迭代后仍可以很好地抑制噪聲且克服MLEM 重建算法收斂慢的缺點。MAP 重建的思想是在傳統的MLEM 算法的基礎上加入先驗約束，從而達到抑制噪聲、平滑圖像、增強圖像邊緣的目的［3］。然而，基于傳統的貝葉斯方法提供的先驗信息是有限的，往往會使低劑量重建圖像出現階梯狀偽影和過平滑現象［4］。本文將從此問題出發對低劑量CT 重建進行進一步的研究。

由于利用圖像稀疏的先驗知識能夠在低劑量的條件下很好的重建圖像，因此在不完備CT 投影數據的重建過程中對得到的梯度圖像求最小化的l1范數，即運用全變差可以改善低劑量CT 重建圖像的質量，該方法雖然提高了重建圖像的抗噪聲性能，但重建圖像依然會出現圖像邊緣不清晰且存在塊狀偽影的問題。小波收縮由于其很好的時頻局部特性，可以在去除噪聲的同時保持圖像的邊緣和細節信息［5］，非局部也有很好的保持邊緣細節的能力，本文受其啟發，提出一種基于小波收縮和非局部的TV 中值先驗重建算法。該算法在每次迭代過程中，用基于TV 的MP重建算法對低劑量投影數據直接進行重建，然后再通過小波變換用小波收縮和非局部對重建后的結果圖進行處理。實驗結果表明，本文算法明顯改善了圖像的質量，同時圖像的信噪比也得到了很好的提高。

1 中值先驗重建算法

Alenius等在1997年提出MRP （median root prior）算法。該算法具有邊緣保持效果，其基本思想為把中值濾波器應用于重建迭代過程中，從而使得圖像的像素值逐漸接近其鄰域像素值的中值［5，6］。但是由于該算法不是真正意義上的MAP算法，只是一種經驗公式。2003年，Hsiao等在MRP的基礎上，提出一種具有優異邊緣保持能力的真正MAP算法。

Hsicao定義的中值先驗的目標函數如下式

式中：Φ（y｜f）——對數型的似然函數，R（f，m）——一種新的先驗分布的目標函數。y 和f 分別表示觀測數據向量和圖像向量，m 是與f 具有相同維數的輔助向量。

先驗分布的目標函數如下

2 基于小波收縮和非局部的TV 中值先驗重建算法

2.1 基于TV的MP重建算法

由于MAP方法引入了圖像的先驗分布信息，將重建結果約束在正則空間，故改善了重建圖像的質量，提高了圖像的信噪比，在一定程度上保持了圖像的細節信息［7］。但是基于傳統的貝葉斯方法只可以提供有限局部先驗信息，因此重建圖像會過平滑，而且會出現階梯狀邊緣的偽影。

由于利用圖像稀疏的先驗知識能夠在欠采樣的條件下很好的重建圖像。而將圖像進行離散梯度變換可得到稀疏圖像，因此在低劑量CT 投影數據的重建過程中對得到的梯度圖像求最小化的l1范數可以改善低劑量CT 重建圖像的質量。而梯度的l1范數即為所謂的全變差。受此啟發，本文將全變差［8－10］引入到MAP的目標函數中，從而提高圖像的抗噪聲性能，基于TV 的MP重建算法的目標函數為

式中：全變差 （total variation，TV）的連續形式為

對于N×M 大小的圖像，其離散形式為

對目標函數進行求解可得

則圖像第j個元素fj的TV 偏導數為

2.2 基于小波收縮和非局部的TV中值先驗重建算法

盡管基于TV 的MP 重建算法取得了很好的效果，但是應用TV 最小化這一先驗依然存在一些問題，對于低劑量的CT 投影數據，重建圖像會出現邊緣不清晰現象且存在著塊狀偽影，為了解決此不足，本文在基于TV 的MP重建算法的基礎上，提出了一種基于小波收縮和非局部的TV 中值先驗重建算法，從而進一步提高了重建圖像的質量。

小波閾值去噪又稱小波收縮［11］，是在最小均方誤差下達到近似最優時對小波系數進行統一處理。圖像進行小波變換后，圖像的大部分能量集中在低頻部分，高頻部分包含了圖像的噪聲和邊緣信息，且只集中了少部分能量，在文獻 ［12］中提出噪聲對應的小波系數幅值在噪聲處較小，在邊緣處較大，因此采用小波收縮的方法處理圖像可以得到令人滿意的效果。和正交小波變換相比，平穩小波變換具有平移不變性從而更適合用于圖像的處理［13，14］。因此本文選用平穩小波變換對圖像進行處理。其中2j為相應的尺度。由此分析，本文對基于TV 的MP 重建算法的每次迭代中，再進行3級離散平穩小波分解，然后對集中較多噪聲的高頻部分進行小波收縮處理，對低頻成分的處理選取可較好保持細節和邊緣的方法，從而保證圖像的細節信息。

Buades等在2005 年利用圖像包含眾多相似結構的特性，提出了一種非局部均值 （nonlocal means，NLM）的圖像去噪算法。該算法的基本思想是在全局范圍內搜索當前像素點所在圖像塊相應的相似塊，然后進行加權平均達到去噪的目的。其具體的計算方法如下

上式是以像素k、像素j 為中心像素點的鄰域內所有像素灰度值的加權Euclidean距離，它為像素i和像素j 的鄰域的一種相似度測量，其中，r 表示離散化噪聲圖像，為某向量的模值；r（Vi）和r（Vj）為灰度值向量，即表示像素i和j 周圍的局部子塊的像素集合。兩者之間的相似性決定像素i和j 之間的相似性。

受其啟發，本文選取非局部作為小波變換后低頻成分的處理方法，實驗表明該算法可以在降噪的同時有效的保持圖像的邊緣和細節，改善圖像的信噪比，提高圖像的抗噪聲性能。

通過上述的分析，本文具體重建算法如下：

（1）基于TV 的MP重建

（2）在基于TV 的MP 重建算法的每次迭代中進行小波收縮處理：對低頻成分采用保持圖像邊緣和細節很好的非局部進行降噪

對高頻成分采用軟閾值收縮，本文選取的閾值為長度對數閾值。最后使用平穩小波逆變換進行重構，得到優化的圖像。

（3）重復以上過程一定次數后，得到最終的重建圖像。

3 實驗結果與分析

3.1 重建圖像比較

本文算法及所有比較算法的實驗仿真環境為：Windows 7旗艦版32位SP1 （DirectX 11）操作系統，英特爾Celeron （賽揚）E3300＠2.50GHz雙核處理器，2GB內存。編程工具使用MATLAB7.6.0 （R2008a）。本文選取大小為128mm×128mm 的Sheep－Logan頭部剖面圖模型作為實驗對象，如圖1 （a）所示。本文選取大小為 （128×128）×（128×128）的系統矩陣，本實驗采用平行投影方式，且在180個角度中均勻選取128個投影方向，每個投影方向分布128對探測器。本文按照下式向理想投影數據加入期望和方差成非線性關系的高斯噪聲來仿真低劑量CT 投影數據

其中，i＝1，2…N 表示探測器信道，N 表示信道總數，λi表示第i個探測器獲得投影數據的均值，σi2表示第i個探測器獲得投影數據的方差，ki表示第i 個探測器的參數，T 表示用于描述掃描系統校準過程的系統參數，對于給定的CT采集系統，ki與T 是給定的。本文參數選取：ki＝200，T ＝12000。為了驗證本文算法的有效性，將本文算法與傳統的MLEM，OS－PML－OSL （ordered subsets－penalized maximum likelihood－one step late），MRP （median root prior），和在每次重建迭代中使用基于方差的各項異性擴散降噪的算法進行了比較。

圖1 各種算法的對比結果

其中圖1 （a）為原始圖像。圖1 （b）采用傳統的MLEM 重建算法。圖1 （c）采用基于有序子集的OSPML－OSL重建算法。圖1 （d）使用本文介紹的MRP 重建算法。為了把本文算法和降噪效果很好的各項異性擴散算法進行比較，圖1 （e）在每次迭代中使用基于方差的各項異性擴散降噪算法對重建圖像進行降噪。圖1 （f）為本文提出的算法。各種算法中涉及到的各種參數以及迭代次數，均為反復實驗后得到的最優值。從結果圖中可得，傳統的MLEM 重建圖像的質量最差；OS－PML－OSL 得到的圖像雖然可以達到降噪的效果，卻也模糊了圖像的邊緣和細節；MRP算法利用中值先驗對重建圖像的噪聲進行了一定程度的消除且可獲得比較清晰的圖像，但是圖像中存在一些明顯的塊狀陰影；基于方差的各項異性擴散降噪處理雖然可達到不錯的效果，但是依然存在少量的塊狀陰影；本文算法不僅可以有效的解決低劑量CT 圖像的噪聲問題，而且較好的保持了圖像的邊緣和細節信息，從視覺上分析，重建效果優于其它幾種算法，初步說明該算法的有效性。

3.2 重建精度比較

從上述分析可知，本文提出的算法對低劑量CT 的重建圖像有很好的去噪以及保持邊緣和細節能力，為了定量描述本文算法的可行性，接下來采用歸一化均方距離、均方絕對誤差、信噪比對重建圖像進行定量描述。其定義如下：

（1）歸一化均方距離 （normalized mean square distance，NMSD）

（2）均方絕對誤差 （mean absolute error，MAE）

（3）信噪比 （signal to noise ratio，SNR）

其中，J 表示圖像像素點的總數，Fi和fi分別表示重建圖像與原始圖像的第i個像素的灰度值，Mi和mi分別表示重建圖像與原始圖像的平均值。這些指標分別從不同的方面評價重建圖像與原始圖像的接近程度以及重建圖像的質量，表1為本文算法與其它幾種比較算法的客觀評價結果。

表1 各種算法的客觀評價

從表1可得出本文算法的信噪比最大，其它指標值均比其它的幾種比較算法小。該結論說明本文算法的重建圖像最接近原始圖像。因此在定量評價方面，同樣表明本算法在CT 重建中是可行的。

圖2給出了本文所用的Sheep－Logan模型的原始理想圖像與在以上各種算法的重建圖像的側面輪廓線的比較圖。從圖中可以看出本文算法的重建圖像與原始圖像的吻合度最高，最接近理想圖像，具有最小的噪聲波動，可以在保持圖像細節信息的同時較好的消除低劑量CT 重建圖像的噪聲。

4 結束語

圖2 各種算法第65行側面輪廓線的對比結果

本文提出了一種基于小波收縮和非局部的全變差中值先驗重建算法。該算法先在中值先驗MP 算法的基礎上，引入降噪性能優異的TV 方法，對目標函數進行修訂，形成基于TV 的MP 重建算法；該算法可以很好的改善圖像的質量，但是依然會有一些塊狀偽影；小波收縮可以在去除噪聲的同時較好的保持圖像的細節信息，非局部也可以很好的保持圖像的細節，結合這幾種方法的優勢，提出了本文的算法，即在基于TV 的MP重建算法的每次迭代中，在小波變換后的小波域再進行小波收縮和非局部降噪。實驗結果表明，該算法改善了圖像的質量，提高了圖像的抗噪聲性能，在主觀效果和客觀效果來看，均說明了該算法的可行性。

［1］Bian J G，Yang K，Boone J M，et al.Investigation of iterative image reconstruction in low－dose breast CT ［J］.Physics in Medicine and Biology，2014，59 （11）：2659－2685.

［2］Ma J H，Zhang H，Gao Y，et al.Iterative image reconstruction for cerebral perfusion CT using apre－contrast scan induced edge－preserving prior ［J］.Physics in Medicine and Biology，2012，57 （22）：7519－7542.

［3］HE Lingjun，PAN Jinxiao，KONG Huihua.Adaptive regularized MAP of CT image reconstruction method ［J］.Computer Engineering and Applications，2011，47 （28）：198－200 （in Chinese）. ［何玲君，潘晉孝，孔慧華.自適應正則MAP 的CT 圖像重建方法研究 ［J］.計算機工程與應用，2011，47（28）：198－200.］

［4］LI Xiaohong，ZHANG Quan，LIU Yi，et al.High quality median prior image reconstruction algorithm based on wavelet shrinkage and forward－and－backward diffusion ［J］.Journal of Computer Applications，2012，32 （12）：3357－3360 （in Chinese）.［李曉紅，張權，劉祎，等.基于小波收縮和正逆擴散結合的優質中值先驗圖像重建算法 ［J］.計算機應用，2012，32 （12）：3357－3360.］

［5］Hsiao I T，Huang H M，Lin K J.A fast OS－type Bayesian reconstruction with an edge－preserving median prior［J］.Journal of Instrumentation，2009，4 （7）：P07012.

［6］Karali E，Koutsouris D.Towards novel regularization approaches to PET image reconstruction ［J］.Journal of Biosciences and Medicines，2013，1 （2）：6－9.

［7］Tang J，Rahmim A.Bayesian PET image reconstruction incorporating anato－functional joint entropy ［J］.Physics in Medicine and Biology，2009，54 （23）：7063－7075.

［8］Liu X W，Huang L H.A new nonlocal total variation regularization algorithm for image denoising ［J］.Mathematics and Computers in Simulation，2014，97：224－233.

［9］Bhadauria H S，Dewal M L.Medical image denoising using adaptive fusion of curvelet transform and total variation ［J］.Computers ＆ Electrical Engineering，2013，39 （5）：1451－1460.

［10］Jovanovski O.Convergence bound in total variation for an image restoration model［J］.Statistics ＆ Probability Letters，2014，90：11－16.

［11］Lee K，Vidakovic B.Semi－supervised wavelet shrinkage［J］.Computational Statistics ＆ Data Analysis，2012，56 （6）：1681－1691.

［12］Fierro M，Ha H G，Ha Y H.Noise reduction based on partialreference，dual－tree complex wavelet transform shrinkage ［J］.IEEE Trans Image Process，2013，22 （5）：1859－1872.

［13］Deng J M，Yue H Z，Zhuo Z Z，et al.A stationary wavelet transform based approach to registration of planning CT and setup cone beam－CT images in radiotherapy ［J］.Journal of Medical Systems，2014，38 （5）：1－9.

［14］Abdullah A J.Denoising of an image using discrete stationary wavelet transform and various thresholding techniques ［J］.Journal of Signal and Information Processing，2013，4 （1）：33－41.