李佳慶,雷 蕾
(1.山西經(jīng)濟管理干部學院,山西 太原 030024;2.山西城市職業(yè)技術(shù)學院,山西 太原 030027)
X-CT 低劑量薄層掃描,即“薄掃”,能更精準地檢出肺癌早期關(guān)鍵病灶。通過X 射線斷層掃描和低劑量計算機斷層掃描(Computed Tomography,CT)技術(shù),對病灶進行篩查,發(fā)現(xiàn)病灶后再進行局部高分辨率掃描的檢查。分析時,把電子探測器接收信號轉(zhuǎn)化為數(shù)字信號輸入計算機,再由計算機轉(zhuǎn)化為圖像。這種掃描已成為肺癌早期診斷最有效、最直接的影像學方式。
本文研究低輻射劑量、薄層掃描技術(shù)的CT 肺部序列結(jié)構(gòu)影像,對連續(xù)高密掃描影像進行檢測分析。由于成像設(shè)備和獲取條件等因素的影響,需要對部分出現(xiàn)質(zhì)量問題的序列圖像進行快速去噪,利用圖像去噪技術(shù)改善CT 圖像的質(zhì)量,在最大可能地保留和增強醫(yī)學細微特征和信息的基礎(chǔ)上,解決不完備投射帶來的影像數(shù)據(jù)質(zhì)量問題;去除低劑量方式帶來的圖像噪聲,在去噪過程中避免將微小結(jié)節(jié)誤清除。
由于醫(yī)學圖像去噪的特殊性,本文擬在研究圖像去噪算法的基礎(chǔ)上,研究保留醫(yī)學細微特征信息,盡可能地去除掉圖像的干擾成分,同時保證原始圖像有用信息的完整度。本文研究算法實現(xiàn)快速有效的圖像去噪處理方法來支持肺部病灶的檢出,實現(xiàn)CT 薄層掃描下序列醫(yī)學圖像快速去噪。
目前,有許多學者在圖像去噪方面做出了研究。MUKHOPADHYAY 等人[1]提出了一種利用遺傳算法的隨機化技術(shù)對閾值進行設(shè)定與優(yōu)化的圖像去噪方法,利用該方法將含噪聲的圖像分成固定大小的小塊,通過小波變換對每個小塊圖像進行去噪。YANG H 等人[2]針對圖像去噪時對邊緣和紋理保護的問題,提出了一種基于雙支持向量機的圖像去噪方法,在對圖像進行去噪的過程中充分利用Shearlet 變換的多分辨率分析和最優(yōu)逼近的特點,對圖像進行去噪的同時很好地保留了圖像的邊緣特征。PROTTER 等人[3]提出了一個通過稀疏和冗余表示對序列圖像去噪的方法,比循環(huán)地處理單張圖像有更好的峰值信噪比(Peak Signal to Noise Ratio,PSNR)。LI H 等人[4]根據(jù)小波變換、稀疏表示和冗余表示,提出了一個單尺度小波K-SVD 算法,對于強噪聲圖像處理有很好的峰值信噪比。
本文考慮薄層掃描切片及切片上肺結(jié)節(jié)的分析,在分析肺部CT 序列上下層切片相關(guān)性的基礎(chǔ)上,結(jié)合專家提供的經(jīng)驗知識,設(shè)計非局部信息稀疏和冗余表示的圖像去噪方法,將圖像序列轉(zhuǎn)化成一個既包含時間相關(guān)度又包含空間相關(guān)度的體素,從而將圖像的去噪問題轉(zhuǎn)換成一個基于冗余字典的最優(yōu)稀疏表示的二次尋優(yōu)問題。通過序列圖像信息與稀疏表示方法對圖像進行去噪,從而將圖像的去噪問題轉(zhuǎn)換成基于冗余字典的最優(yōu)稀疏表示的尋優(yōu)問題[5]。
時域和頻域是信號的基本性質(zhì),為信號分析提供了不同的角度。對于需要處理的信號來說,在時域上本身就具有稀疏性的信號很少。為了尋找更好的稀疏,總能找到某種變換,使得在某種變換域中稀疏,即能夠用較少的系數(shù)來表達。
常見的變換有正交變換、余弦基變換(Discrete Cosine Transformation,DCT)、小 波 變 換(Wavelet Transformation,WT)、多分辨率分析(小波變換)、多 尺 度 幾 何 分 析[6](包 括Ridgelet、Curvelet、Bandelet、Contourlet 等一系列多尺度分析工具)。這些工具符合人類視覺皮層對圖像有效表示的要求,通過變換的數(shù)學分析方法把圖像分解在不同的尺度(不同的精細程度)上來處理。
含噪的圖像是有一定的結(jié)構(gòu)的,可以從信號中提取出能表達圖像性質(zhì)的原子,即將信號表示為字典和稀疏的乘積,將信號提取出來,即為去噪后的圖像。而噪聲是隨機、不相關(guān)的,沒有結(jié)構(gòu)特性,因此提取不出來。通過這樣的原理可以將圖像提取出來,以達到圖像去噪的目的。
稀疏表示[7]是獲得對信號的良好近似,優(yōu)化模型是從信號重建的角度來設(shè)計的,將信號表示為字典的線性表示。即,y=Ax,其中y是待處理的信號,A是字典,x是稀疏系數(shù)。稀疏(性)是指x的大部分元素為0 或接近0,只有少部分元素不為0。其中,x滿足稀疏性。稀疏表示的目的是選擇最少的系數(shù),從而重構(gòu)信號y。已知量是輸入信號y,未知量是字典A與稀疏系數(shù)x。對信號稀疏表示的目的是尋找一個字典使得信號在該字典下的系數(shù)(表達)最稀疏。
稀疏表示中稀疏的原因是:字典是在一個足夠大的訓練樣本空間內(nèi),對于一個類別的物體,可以大致由訓練樣本中的樣本子空間線性表示,因此當該物體有整個樣本空間表示時,其表示的系數(shù)是稀疏的[8]。
稀疏表示的兩個主要的問題是字典的生成和稀疏分解(編碼)。
2.3.1 字典的生成
字典是一個矩陣,一個原子(基)是字典A矩陣中的一列。由這樣一列一列的原子排列成的矩陣就是一個字典。字典表達信號的能力取決于信號的特征是否與字典中基的特征相吻合。字典是基的組合,有表達更多、更復雜的信號的能力。
字典有固定字典和自適應字典兩種。常用的固定字典有超完備字典DCT、小波字典Contourlet 和曲波字典Wavelet 等。固定字典不能隨著信號的變化做出相應的變化,一經(jīng)選定,對所有的信號一視同仁,信號的表達形式單一;而學習字典是通過大量圖像數(shù)據(jù)學習得到的,通過訓練學習得到的字典可以根據(jù)輸入信號的不同做出自適應的變化,能夠更好地適應不同的圖像數(shù)據(jù)。
2.3.2 稀疏分解(編碼)
求信號在字典上的稀疏表示,即求出信號y在字典A上的稀疏表示x,常用的方法是正交匹配追蹤 算 法(Orthogonal Matching Pursuit,OMP)。主要目標是找出x中最主要的k個分量和值。假設(shè)x=(x0,0,0,x1,x2,…)。為尋求與y最為接近的x,首先需計算信號與原子(字典中的每一列)之間的內(nèi)積。內(nèi)積可理解為計算矢量在某一指定方向上的投影長度。其次,確定最大值對應的原子,并利用最小二乘法來確立其稀疏系數(shù)。在接下來的步驟中,將余下的原子與字典中相應的列進行內(nèi)積計算,以識別出與最大值匹配的字典列,不斷迭代,直至達到算法的收斂狀態(tài)。
稀疏分解算法的概念最早由MALLAT 提出,名為匹配追蹤算法(Matching Pursuit,MP),以其簡潔和易于實現(xiàn)的特性得到了廣泛的應用[9]。此后,PATI 等學者在MP 算法的基礎(chǔ)上提煉出OMP,其特點是具備更大的收斂速度。隨后的研究進一步優(yōu)化了OMP 算法,孕育出如壓縮采樣匹配追蹤算法(Compressive Sampling Matching Pursuit,CoSaMP)、正則化正交匹配追蹤算法(Regularized Orthogonal Matching Pursuit,ROMP)、分段式正交匹配追蹤算法(Stagewise OMP,StOMP)以及子空間追蹤算法(Subspace Pursuit,SP)等多種變種。
2.3.3 優(yōu)化問題
2.3.3.1 凸優(yōu)化問題
字典A是一個N×M的矩陣。由于M遠小于N,用M個方程求解N個未知數(shù),即用少的方程求解多的未知數(shù),因此y=A*x是個欠定方程(有無窮多解),說明使用字典的線性組合表達信號將是不唯一的,所以求解此問題需要找到最優(yōu)解。
解優(yōu)化問題,如果加上適當?shù)南薅l件即有唯一解。那么需要解決的問題就是尋找是否存在一種最好的表達方式,使得在字典的表示下系數(shù)最稀疏。||x||0為x的稀疏度,表示x中非零系數(shù)的個數(shù)。要求解x盡可能地稀疏,即||x||0盡可能地小。0 范數(shù)的優(yōu)化問題很難求解,而1 范數(shù)優(yōu)化問題容易求解(存在且唯一)。在滿足一定條件下將問題轉(zhuǎn)換為1范數(shù)的優(yōu)化問題(最小化問題),此時優(yōu)化問題變成一個凸優(yōu)化問題[10],應用已有的理論就可以解決。
2.3.3.2 最大期望算法
優(yōu)化的模型是兩個參數(shù)相乘,需要同時估計字典A和稀疏表示系數(shù)x,基本思想是:第一步,將字典A固定,求出x的值,即求信號y在A上稀疏表示;第二步,使用上一步得到的x來更新字典A,即字典更新,使得y和Ax最接近。如此反復迭代幾次,即可得到優(yōu)化的A和x。
在稀疏表示里,已知y和A求x的方法是OMP。已知y,求一個比較好的A的方法是Module(MOD)和K-SVD 等。字典的更新方法主要有最 大 似 然(Maximum Likelihood,ML)、最 大 后驗 概 率(Maximum A posteriori Probability,MAP)、MOD、FOCUSS 字典學習算法、主成分分析(Principal Component Analysis,PCA)、K-SVD 以 及Online 等。本文以應用最廣泛的K-SVD[11]算法為例做對比實驗。
輸入含噪序列圖像后,K-SVD[12]訓練字典用于圖像去噪的算法步驟如下。
(1)選擇字典作為過完備字典,訓練字典;選取字典集,并更新字典。
(2)將含噪圖像矩陣以m×m為一塊,將一個大的矩陣劃分為多個子矩陣。對于每個子矩陣,將元素按列連接為一個列向量,一個m×m的塊轉(zhuǎn)換為m2×1 的列向量,生成一個新的矩陣B。最后生成的矩陣B就是由A個塊矩陣生成的A列矩陣。矩陣B就作為字典學習的輸入。
(3)先取B中的M列,更新其中的每一個元素,B[i][j]B[i][j]-means{B[:,j]}B[i][j]=B[i][j]-means{B[:,j]},式中B[i][j]表示B[i][j]矩陣中第i行第j列元素,means{B[:,j]}表示元素(i,j)所在列(第j列)的平均值。
(4)固定字典,利用OMP 算法,求解該B在字典下面的稀疏系數(shù)cofes。
(5)更新B的每一列:B[:][j]=B[:][j]+means{B[:,j]},在此上下文中,B[:][j]第j列可被視為原子的代表。K-SVD 策略涉及逐一地更新原子(即字典的某一特定列)及其相應的稀疏系數(shù)。只有當所有的原子都經(jīng)歷了更新過程后,該算法才進入下一迭代。經(jīng)過多次的精細迭代,最終能夠獲得經(jīng)過優(yōu)化的字典及其稀疏系數(shù)。
(6)取下一個M列,重復步驟(3)~(5),直到B矩陣的所有列都處理完成。
(7)將B的每一列(m2×1)轉(zhuǎn)換為m×m的小塊,逆運算得到去噪后的圖像。
最終,輸出去噪后的序列圖像。
本文采用模擬退火算法[13]得到最優(yōu)值,算法是一種隨機性組合優(yōu)化方法,來源于對固體退火降溫的過程。即將固體加溫至充分高,再讓其溫度慢慢退卻。固體加熱時,固體內(nèi)粒子隨溫度的升高而變?yōu)闊o序狀,而冷卻時,粒子逐漸趨于有序。這一過程中,每個溫度下都達到平衡狀態(tài),不但能夠解優(yōu)化問題,而且能夠優(yōu)化字典中原子的組合,使字典具備更好的表達信息的能力,從而更好地去噪。
算法是從整體上考慮的,概率與系統(tǒng)能量E、溫度T密切相關(guān)。在高溫下,接受與當前狀態(tài)能量差較大的新狀態(tài);在低溫下,接受與當前狀態(tài)能量差較小的新狀態(tài)。這樣避免陷入局部最優(yōu),有的放矢地高效找出最優(yōu)解。所提方法描述如圖1 所示。模擬退火算法[14]稀疏表示去噪算法具體步驟如下。

圖1 方法描述
首先輸入y=A*x中A與x的無窮多個解,限定條件||x||0和||x||1盡可能小(稀疏)。
(1)初始化狀態(tài)信息,隨機產(chǎn)生初始解x0(即在溫度T條件下的解為x0)。
(2)設(shè)目標(能量)函數(shù)E=y-y*,其中y*=Ax,計算差值ΔE=Ej-Ei,i為當前狀態(tài),j為下一狀態(tài)。若ΔE<0,則接受j為當前狀態(tài),即更接近最優(yōu)解(內(nèi)能變小更趨于有序狀態(tài)),否則以一定的概率來接 受 狀 態(tài)。min{1,exp(-ΔE/KT)}>random[0,1],K是Metropolis 接受準則,是以概率接受新狀態(tài)。溫度T為控制參數(shù)。對于exp(A)函數(shù),當A<0,其值小于1;當A>0,其值大于1。ΔE<0,-ΔE>0,-Δf/TK>0,exp(Δf/TK)>1,min{1,exp(-Δf/TK)}取1,接受全部新值。當解為差解時ΔE>0,-ΔE<0,-ΔE/KT<0,exp(-ΔE/KT)<1,min{1,exp(-ΔE/KT)}中取exp(-ΔE/KT),最后通過概率p=exp(-ΔE/KT)>random[0,1]來接受新值,即概率p=exp[-(Ej-Ei)/KT]若大于[0,1]區(qū)間的隨機數(shù),接受狀態(tài)j為當前狀態(tài);若不成立,則保留狀態(tài)i為當前狀態(tài)。
(3)不斷產(chǎn)生新解,計算目標函數(shù)差,判斷是否接受新解。經(jīng)過大量的解變化后,可以求得給定控制參數(shù)的優(yōu)化問題的相對最優(yōu)解。
(4)設(shè)定閾值n是最優(yōu)的n個解。保留這些解,尋找與圖像最接近的若干解,并對這些解進行標記。利用模擬退火算法來優(yōu)化自適應字典,以找到每個解的最佳自適應字典。最終,在所有解中選擇最優(yōu)解。
輸出最優(yōu)解,得到最優(yōu)的字典A和稀疏系數(shù)x,逆運算Ax=y*得到去噪后的圖像,從而完成圖像去噪。
與廣泛使用的K-SVD 算法相比,K-SVD 算法需要先更新每一個元素再更新每一列,更新速度慢,迭代次數(shù)多。而本文算法是從宏觀的角度動態(tài)自適應地更新協(xié)調(diào)全部的值,更加整體和智能,在對比實驗中得到的去噪效果較好。
算法的實驗環(huán)境是Matlab R2015b 軟件,計算機處理器為Intel Core i7-3770,主頻3.40 GHz,內(nèi)存8 GB。用數(shù)學軟件進行算法研究和數(shù)據(jù)分析,從評價圖像質(zhì)量的指標峰值信噪比PSNR 和算法性能的時間復雜度兩方面進行實驗。實驗中,將提出的基于模擬退火算法的稀疏表示圖像去噪算法與K-SVD 算法進行去噪效果對比。實驗數(shù)據(jù)為采集的山西省人民醫(yī)院患者數(shù)據(jù),采集設(shè)備是醫(yī)院現(xiàn)有X-Ray 薄掃設(shè)備。依次標記20 組序列圖像,將每組序列圖像數(shù)據(jù)集作為一組測試集?,F(xiàn)列出一組序列圖像的實驗數(shù)據(jù)來說明算法[15]。
訓練的字典如圖2 所示。圖2(a)是K-SVD算法的自適應字典,圖2(b)是本文算法所得出的字典。

圖2 訓練字典
序列圖像中其中一幅圖像去噪后的實驗結(jié)果對比如圖3 所示。從圖3 可以看出,本文算法的實驗結(jié)果有較高的峰值信噪比,從直觀感受方面看得到了很好的去噪效果。

圖3 噪聲圖像與去噪后的圖片對比
從整體看,一組序列圖像的實驗數(shù)據(jù)如圖4、圖5 所示。

圖4 本文方法與K-SVD 方法的信噪比對比

圖5 本文方法和K-SVD 方法時間復雜度的對比
從實驗數(shù)據(jù)可以看出,對于同一序列圖像,傳統(tǒng)的K-SVD 方法在序列圖像的信噪比上有間歇性的不穩(wěn)定,對于某些圖像的去噪效果不好,但本文方法有很好的信噪比與穩(wěn)定性。傳統(tǒng)的K-SVD 方法時間復雜度較高,而本文方法在時間復雜度上有所降低[16]。
接下來,擴展至更多的數(shù)據(jù)集,使實驗更具廣泛性和穩(wěn)定性[17-18]。以下是一組序列圖像的實驗結(jié)果,運行程序共運行了20 組實驗數(shù)據(jù),現(xiàn)列出5組數(shù)據(jù)平均值,以說明實驗結(jié)果及去噪效果。將本文算法與現(xiàn)在廣泛使用的K-SVD 的相關(guān)數(shù)據(jù)對比,結(jié)果如表1、表2 所示。表1 給出了其中5 組序列圖像的去噪峰值信噪比平均值??梢钥闯?,在對薄掃圖像進行去噪時,本文提出的方法具有較高的信噪比。

表1 峰值信噪比平均值對比 單位:dB
時間復雜度是評估算法性能的重要參數(shù)。實驗中使用的算法對樣本的時間復雜度如表2 所示。在實驗比較的算法中,本文算法的速度比K-SVD 算法快。在測試階段,本文算法識別一個手勢圖像序列消耗的平均時間為99.39 s。
在將快速去噪算法應用于整個肺部結(jié)節(jié)的檢出實驗中,根據(jù)先驗知識與診斷結(jié)果,對是否檢出病灶的實驗數(shù)據(jù)進行比對,進而算出結(jié)節(jié)識別率進行對比,對于病灶不明顯的圖像具有一定的作用,實驗數(shù)據(jù)存在偶然性,就目前所具有的數(shù)據(jù),識別的結(jié)果有存在不一致的情況,這樣降低漏診率,在預處理方面對此課題有一定的輔助作用[19]。相對于傳統(tǒng)的K-SVD 方法的性能優(yōu)勢的描述,包括信噪比的提高、穩(wěn)定性的增強以及時間復雜度的降低。這些描述強調(diào)了本文方法在圖像去噪任務中的優(yōu)越性,并提供了實驗數(shù)據(jù)支持,進一步強調(diào)了本文方法的創(chuàng)新性和實用性,使讀者更容易理解本文方法是一個有價值的改進,為圖像去噪帶來了重要的進展。因此,本文提出的算法更適用于稀疏表示去噪,是一個行之有效的方法,在效果和性能方面都有一定的提高。
將稀疏表示用于去噪,分為多尺度幾何分析、稀疏表示和優(yōu)化問題。其中,稀疏表示的兩大任務是字典的生成和信號的稀疏分解。運用多尺度幾何分析將圖像轉(zhuǎn)換到某個域中,使稀疏盡可能地稀疏,能夠達到好的效果。本文運用優(yōu)化方法進行解優(yōu)化問題。在字典生成中,下一步的工作是用online 子空間聚類去批量去噪。對稀疏表示去噪還會繼續(xù)進行研究,以為醫(yī)學診斷做好堅實的前期工作。肺癌死亡率高的主要原因是其早期的細小病變特征很難被發(fā)現(xiàn),臨床漏診率高。肺癌早期的檢出對提高治愈率起著至關(guān)重要的作用,能夠使達到醫(yī)學圖像診斷誤差率盡可能小。