基于字典學習的低劑量X-ray CT圖像去噪

朱永成 陳 陽 ，2，3 羅立民，2，3 Toumoulin Christine，2，3

(1東南大學影像科學與技術實驗室，南京210096)

(2中法生物醫(yī)學信息研究中心，法國雷恩35000)

(3法國雷恩大學信號與圖像處理實驗室，法國雷恩35042)

計算機斷層成像(computed tomography，CT)技術自1970年初面世以來在醫(yī)學診斷、治療中的應用越來越廣泛.與其他常規(guī)影像檢查手段(如X射線)相比，CT技術能夠提供具有較高分辨率的斷層解剖學圖像，但是其較高的輻射劑量會對人體帶來一定的傷害.低劑量 CT(low-dose CT，LDCT)圖像一般可以通過降低X射線管電壓或管電流的方法得到.由于塊狀噪聲的影響，圖像質(zhì)量會顯著下降.因此，如何有效抑制LDCT圖像中的噪聲成為CT技術在臨床應用中的關鍵［1］.改進CT圖像質(zhì)量的方法可概括為以下3類:①對重建算法進行改進，減小噪聲對重建過程的影響.該類算法的研究有賴于投影數(shù)據(jù)及其相關信息的獲得.②對投影數(shù)據(jù)進行去噪，再利用傳統(tǒng)的濾波反投影算法進行重建.此技術的關鍵在于理解不同劑量CT投影數(shù)據(jù)的噪聲統(tǒng)計特性，并據(jù)此設計有效算法以抑制投影數(shù)據(jù)中的噪聲［2］.③ 在圖像空間中對重建圖像進行噪聲抑制.由于CT成像過程是經(jīng)過反投影重建得到的，其噪聲統(tǒng)計特性比較復雜，常規(guī)的圖像空間處理方法難以取得較好的去噪效果［3］.

近年來，在自然圖像分析領域，基于自適應稀疏表示理論的方法得到越來越多的關注.該理論起源于Olshausen等［4］提出的一種通過學習得到過完備字典的算法，該算法從一組自然圖像中提取小圖像塊作為訓練樣本，通過計算圖像塊稀疏表示的方式得到圖像字典.字典中的元素具有類似于簡單細胞反應的特性，因此可以利用這種字典來表示自然圖像.隨著機器學習理論的發(fā)展，基于圖像樣本學習構造自適應圖像字典并將其用于圖像去噪［5］、壓縮［6］和恢復［7］等領域的自適應稀疏方法得到了迅速且廣泛的應用.1999年，Engan等［8］采用最優(yōu)方向算法(MOD)來解決稀疏表示中的字典設計問題.在這一優(yōu)化過程中，系數(shù)更新和字典更新交替進行，其中系數(shù)的更新通過現(xiàn)有的各種稀疏優(yōu)化算法逐列進行，而字典的更新則采用簡單的整體求導方式.算法每一次迭代可以得到一個局部最小解，幾十次迭代后即可達到收斂.這一算法的缺點在于涉及矩陣求逆運算，計算復雜度較大.2005年，Aharon等［9］在MOD算法基礎上提出了K-SVD算法，該算法采用逐列更新字典的方式，解決MOD算法中矩陣求逆的問題，通過對部分圖像集進行迭代學習，并在修正樣本中稀疏編碼的同時不斷更新字典，使其更好地符合樣本集的固有特性，有效降低了計算的復雜度.除此之外，近年來還涌現(xiàn)出了很多參數(shù)化和結構化的字典，如平移不變字典［10］、多分辨率字典［11］等.

本文首先將K-SVD算法應用于體模圖像去噪中，通過比較去噪結果圖像的PSNR值，確定該算法中的參數(shù)設定和字典訓練時的樣本選擇;然后將該算法應用于病人圖像的去噪研究中，以證明算法的有效性.實驗結果表明，該算法在低劑量CT圖像去噪中取得了較好的效果，得到了對臨床診斷有意義的去噪結果.

1 基于字典學習方法的算法

K-SVD算法可分為以下3步:

①字典的訓練.旨在尋找最優(yōu)字典來對訓練樣本進行稀疏表示，使基于此字典的稀疏表示對于樣本的逼近均方差達到最小.

②求解稀疏表示.基于步驟①得到的字典，利用某一貪婪算法(如正交匹配跟蹤算法)求解噪聲圖像的稀疏表示.

③圖像的重建.將步驟①得到的字典與步驟②得到的稀疏表示相乘，得到結果圖像，即可完成圖像的重建.

K-SVD算法中最關鍵的步驟是圖像字典的訓練.要達到此目的，需構造如下的目標函數(shù):

式(2)中L的取值決定于LDCT圖像的質(zhì)量，此處根據(jù)LDCT圖像的標準差來決定其取值.對于處理后的圖像，~yi=D~xi.在適當?shù)淖值湎?，無噪聲圖像符合一定稀疏表示的特點，而噪聲破壞了這種稀疏表示性;式(2)中定義的模型通過估計無噪聲圖像的稀疏表示系數(shù)來恢復圖像，從而達到抑制LDCT圖像中噪聲的目的.

在字典更新階段，首先計算式(2)中定義的模型.因為該模型是NP難問題，故采用正交匹配跟蹤、FOCUSS或基跟蹤等算法求近似解，計算殘差后對殘差應用奇異值分解，更新字典和系數(shù)矩陣，迭代更新以獲得更優(yōu)的字典.在去噪階段，字典是已知的，因此再次計算噪聲圖像在該字典下的稀疏表示，即得到相應解.需要注意的是，算法中的計算都是針對圖像中指定大小的圖像塊，且這些圖像塊來源于不斷平移得到的重疊的圖像塊集合，因此最終的實驗結果是將這一圖像塊集合中所有元素的計算結果進行加權平均得到的.

2 實驗與結果分析

2.1 實驗數(shù)據(jù)

本實驗中所使用的LDCT圖像是管電流為60 mA時得到的 CT圖像，而在臨床中使用的NDCT圖像是管電流為240 mA時得到的CT圖像.該LDCT圖像來自于醫(yī)院采集的病人信息，包括一組肺窗數(shù)據(jù)和一組縱膈窗數(shù)據(jù).此外，為了定量分析該算法的去噪效果以及參數(shù)對于實驗結果的影響，對一組體模數(shù)據(jù)進行去噪處理，并計算去噪前后的峰值信噪比PSNR.PSNR作為一種評價圖像的客觀標準，計算公式如下:

式中，Mp為圖像的最大可能像素值，由圖像的位數(shù)決定，例如對于一幅8 bit的灰度圖像，Mp=28-1=255;MS為標準圖像與處理后圖像之間的均方誤差.對于2幅大小為p×q的圖像P和Q，MS的計算公式為

本實驗中，用于計算均方誤差的對比圖像是與LDCT圖像相對應的、管電流為240 mA時得到的NDCT圖像.由于病人在2次掃描過程中不可避免會因位置移動或呼吸等動作，造成NDCT圖像難以與LDCT圖像完全一致，因此為了確保PSNR定量分析的準確性，需針對體模數(shù)據(jù)進行處理.CT掃描器的參數(shù)設置如下:管電壓為120 kVp;掃描片厚度為2 mm;旋轉時間為0.5 s;肺窗和縱膈窗的重建算法核分別為B70f和B31f;CT劑量和管電流的大小成線性正比關系，240 mA條件下記錄的劑量為9.36 mGy，60 mA 條件下為 2.34 mGy.所處理的CT圖像數(shù)據(jù)以DICOM文件存儲，用于處理的主機配置如下:Intel CoreTM2 Duo CPU;主頻為3.00 GHz，內(nèi)存為2 GB;顯卡為 NVIDIA GeForce 8 400 GS.

2.2 實驗參數(shù)設置

本算法要求在字典訓練和稀疏表示過程中僅針對較小圖像塊進行處理，且塊的大小受到處理機內(nèi)存的制約.此處選取圖像塊大小為8×8.根據(jù)冗余字典理論，在實驗中訓練得到的字典大小設定為64×441，其中64為圖像塊的大小，即每個圖像塊排成列向量的行數(shù);K=441表示字典的原子個數(shù).本算法對每個像素點進行O(nKLJ)次操作，其中J為字典訓練時的迭代次數(shù)，n為字典中圖像塊大小.圖像噪聲未知，需要對其進行估計.最常用的方式是局部標準差法，即通過合理的數(shù)學模型進行估計.當標準差σ=10時，L的平均值為2.96;當σ=20 時，L 的平均值為 1.12.Elad 等［12］通過實驗將拉格朗日參數(shù) λ取為30/σ，噪聲增益 C取為1.15;本文采用相同取值.

字典訓練時，有2種源數(shù)據(jù)的選擇方式［6］:①選取噪聲圖像作為訓練的源數(shù)據(jù)，根據(jù)待處理圖像的不同特點來調(diào)整字典，得到每幅圖像的自適應字典;②選取一組與LDCT噪聲圖像掃描部位相近的NDCT圖像集合作為訓練數(shù)據(jù).

2.3 實驗結果

分別以待處理的LDCT圖像和一組先前保存的NDCT圖像作為訓練源數(shù)據(jù)，根據(jù)2.2節(jié)中的參數(shù)設置，分別得到肺窗和縱膈窗圖像顯示下的去噪結果.算法中參數(shù)取不同值時，字典訓練時間不同，兩者之間的關系見表1.產(chǎn)生這一現(xiàn)象的原因在于，L的取值與噪聲水平相關，圖像的噪聲水平越高，L的取值則應盡量小，理論上需要的迭代次數(shù)應盡量多.訓練字典的去噪過程是對每一個圖像塊進行稀疏表示，其運行時間取決于所選取圖像塊以及噪聲圖像整體的大小.本實驗中，圖像塊大小為8，圖像大小為512 ×512，所需時間約為0.25 h.

表1 字典訓練的時間 h

2.3.1 字典的訓練

假定L=3，迭代次數(shù)為20，初始字典設置為如圖1(a)所示的DCT字典，由此便可得到如圖1(b)和(c)所示的圖像字典.

實驗所用參數(shù)中，有2個參數(shù)是與字典相關的，即字典大小K和字典中圖像塊大小n.實驗中首先保證其他參數(shù)值不變，分別改變這2個參數(shù)值，計算體模數(shù)據(jù)(縱膈窗和肺窗)的PSNR，結果見表2.由表可知，以NDCT圖像作為字典訓練的數(shù)據(jù)源，字典越大且字典內(nèi)圖像塊越大時，去噪效果越好.這是因為圖像塊和字典尺寸越大，可包含的圖像信息越多，從而能夠更好地對噪聲圖像進行稀疏表示.以LDCT圖像自身作為訓練數(shù)據(jù)源時，圖像塊越大，去噪效果越好;但是如果字典太大時，其去噪效果反而降低了.其原因在于，LDCT圖像中有偽影或噪聲信息，字典越大，則會導致越多的信息被引入到字典中，從而使得去噪效果變差.

圖1 初始以及訓練中得到的圖像字典

表2 不同參數(shù)下結果圖像的PSNR dB

2.3.2 體模數(shù)據(jù)實驗結果

圖2和圖3分別是體模肺窗和縱膈窗圖像的去噪結果.圖中，方式1和方式2分別是指以LDCT圖像和NDCT圖像作為訓練樣本訓練字典.實驗結果驗證了不同的訓練數(shù)據(jù)源及不同參數(shù)對去噪效果的影響.根據(jù)表2數(shù)據(jù)分析可知:本文算法可以有效地提高LDCT圖像的質(zhì)量;同時，以NDCT圖像作為訓練數(shù)據(jù)源時，可以取得更好的去噪結果，這是因為NDCT圖像中含有較少的噪聲或者偽影，構建的字典能夠更好地去除LDCT圖像中的噪聲.

2.3.3 病人真實數(shù)據(jù)實驗結果

圖4和圖5分別為病人肺窗和縱膈窗圖像的去噪結果.由圖可知，K-SVD算法可以有效去除LDCT圖像中的噪聲，并保留了被掃描部位人體組織的特性，有利于臨床上做出正確的診斷.

圖2 體模肺窗圖像的去噪結果

圖3 體?？v膈窗圖像的去噪結果

圖4 病人肺窗圖像的去噪結果

圖5 病人縱膈窗圖像的去噪結果

3 結語

本文采用K-SVD算法，對體模圖像和病人圖像分別進行去噪處理和定量分析.在體模圖像實驗中，研究了K-SVD算法中的相關參數(shù)對結果的影響.同時，以不同的數(shù)據(jù)源作為字典訓練樣本進行試驗，結果表明，以NDCT圖像作為訓練樣本的去噪結果優(yōu)于以LDCT圖像作為訓練樣本的結果.病人圖像的去噪結果表明，本文算法可以在有效去噪的同時盡量保留圖像細節(jié)，有利于臨床上依據(jù)結果圖像作出診斷.下一步工作主要是對該算法進行進一步優(yōu)化，提高算法的執(zhí)行速度;利用近年來提出的參數(shù)化字典理論，訓練得到結構更優(yōu)的圖像字典，從而進一步提高結果圖像的圖像質(zhì)量.

References)

［1］Yazdi M，Beaulieu L.Artifacts in spiral X-ray CT scanners:problems and solutions［J］.Proc World Acad Sci Eng Technol，2007，26(3):376-380.

［2］Wang J，Lu H，Wen J，et al.Multiscale penalized weighted least-squares sinogram restoration for low-dose X-Ray computed tomography［J］.IEEE Transactions on Biomedical Engineering，2008，55(3):1022-1031.

［3］Rust G F，Aurich V，Reiser M.Noise/dose reduction and image improvements in screening virtual colonoscopy with tube currents of 20mAs with nonlinear Gaussian filter chains［J］.Proc SPIE Int Soc Opt Eng，2002，4683:186-197.

［4］Olshausen B A，F(xiàn)ield D J.Emergence of simple-cell receptive field properties by learning a sparse code for natural images［J］.Nature，1996，381(6583):607-609.

［5］Starck J L，Candes E J，Donoho D L.The curvelet transform for image denoising［J］.IEEE Transactions on Image Processing，2002，11(6):670-684.

［6］Wakin M B，Romberg J K，Choi H，et al.Wavelet-domain approximation and compression ofpiecewise smooth images［J］.IEEE Transactions on Image Processing，2006，15(5):1071-1087.

［7］Needell D，Tropp J A.CoSaMP:iterative signal recovery from incomplete and inaccurate samples［J］.Applied and Computational Harmon Analysis.2009，26(3):301-321.

［8］Engan K，Aase S O，Hankon H J.Method of optimal directions for frame design［C］//Proceedings of 1999 IEEE International Conference on Acoustics，Speech and Signal Processing.Phoenix，AZ，USA，1999:2443-2446.

［9］Aharon M，Elad M，Bruckstein A.The K-SVD:an algorithm for designing overcompletedictionariesfor sparse representation［J］.IEEE Transactions on Signal Processing，2006，54(11):4311-4322.

［10］Blumensath T，Davies M.Sparse and shift-invariant representations of music［J］.IEEE Transactions on Speech Audio Processing，2006，14(1):50-55.

［11］Rubinstein R，Zibulevsky M，Elad M.Double sparsity:learning sparse dictionaries for sparse signal approximation［J］.IEEE Transactions on Signal Processing，2010，58(3):1553-1564.

［12］Elad M，Aharon M.Image denoising via sparse and redundantrepresentations overlearned dictionaries［J］.IEEE Transactions on Image Processing，2006，15(12):3736-3745.