陳明惠,王 帆,張晨曦,李福剛,鄭 剛
(1. 上海理工大學(xué) 生物醫(yī)學(xué)光學(xué)與視光學(xué)研究所 教育部現(xiàn)代微創(chuàng)醫(yī)療器械及技術(shù)工程研究中心 上海介入醫(yī)療器械工程技術(shù)研究中心,上海 200093;2. 上海奧普生物醫(yī)藥有限公司,上海 201203)
光學(xué)相干斷層掃描(Optical Coherence Tomography,OCT)成像技術(shù)利用寬帶光源的低相干特性,以非侵入、無損傷的方式進行斷層成像,已廣泛應(yīng)用于醫(yī)學(xué)診斷和研究[1-2]。頻域OCT技術(shù)因其高靈敏度和高速性漸漸取代傳統(tǒng)時域OCT技術(shù),但是擁有這些優(yōu)點的同時意味著后續(xù)數(shù)據(jù)采集和處理系統(tǒng)的造價更高且負擔(dān)更重[3];同時,該成像技術(shù)需要對病變部位進行反復(fù)多角度多層次的深度掃描以獲得大量高分辨率圖像,高分辨率必然伴隨著高數(shù)據(jù)量,同時在對眼底部位進行成像時,眼球的抖動會造成運動模糊,所以迫切需要一種能降低硬件系統(tǒng)存儲和傳輸代價且花費時間較少并獲得較高圖像質(zhì)量的方式。
2006年,Donoho[4],E Candes[5-6]及華裔科學(xué)家Tao[7]等人提出了一種新的信息采集方式,即壓縮感知(Compressive Sensing,CS)理論。CS技術(shù)表明,如果信號(圖像)本身是稀疏的,或能在某個稀疏集上進行稀疏表示,則通過一定的測量矩陣進行測量,加以選取適當(dāng)?shù)膬?yōu)化算法,在獲得少量測量值的基礎(chǔ)上,能對原始信號進行有效地恢復(fù)重構(gòu)[8]。研究高性能的恢復(fù)重構(gòu)算法是CS領(lǐng)域的一個熱點,目前優(yōu)化算法主要分為貪婪算法、凸優(yōu)化算法和非凸算法三類[9]。研究者們在匹配追蹤算法(Matching Pursuit,MP)的基礎(chǔ)上獲得了更多的貪婪算法,Wei等人[10]提出的SAMP(Sparsity Adaptive Matching Pursuit)重構(gòu)算法通過更新支持集并逐漸增加稀疏度來逼近原始信號,實現(xiàn)了稀疏度未知的信號重構(gòu),但重構(gòu)精度不夠高,重構(gòu)后SNR最高僅為25 dB左右;Chen等人[11]利用小波變換獲得圖像的稀疏表示,使用改進的最佳正交匹配追蹤算法完成MRI圖像的重構(gòu),但算法計算效率不高,算法運行時間約為26 s;Ma等人[12]提出一種變步長SAMP算法,在一定程度上提高了重構(gòu)精度,但引入分塊思想后,圖像重構(gòu)出現(xiàn)明顯的人為噪聲,導(dǎo)致重構(gòu)質(zhì)量下降;石昊蘇等人[13]采用共軛梯度算法替換OMP算法中的最小二乘法求取估計值進行超聲圖像的重構(gòu),提高了重構(gòu)質(zhì)量,但PSNR值僅達到13.75 dB,算法運行時間約為16.6 s。第二類凸優(yōu)化算法通過將非凸問題轉(zhuǎn)換成凸問題進行求解,Qu等人[14]將基追蹤算法和非均勻快速傅立葉變換技術(shù)相結(jié)合,實現(xiàn)了觀測場景的成像重構(gòu),提高了成像質(zhì)量,但算法運行時間約為39 s,算法計算復(fù)雜度高,應(yīng)用范圍受限;鄭萬澤等人[15]基于Contourlet變換,采樣改進梯度投影算法(Gradient Projection for Sparse Reconstruction,GPSR)恢復(fù)稀疏處理后的系數(shù),實現(xiàn)了二維Lena和Barbara圖像的重構(gòu),采樣率為0.5時重構(gòu)圖像的PSNR值僅為30 dB左右且算法實時性較差。第三類非凸算法通過改變目標(biāo)函數(shù)求解優(yōu)化問題,該類算法克服了基追蹤算法中的計算量大問題,同時可以采用比正交匹配追蹤算法少的測量值達到同樣的重建效果,主要以FOCUSS算法[16]和稀疏貝葉斯學(xué)習(xí)[17]為代表。
上述文獻中改進的壓縮感知重構(gòu)算法在保證算法的實時性以及圖像最終重構(gòu)質(zhì)量之間總是存在矛盾,算法運行速度最快約為17 s,PSNR值最高為30 dB。為了在一定程度上實現(xiàn)頻域OCT圖像重構(gòu)效率和重構(gòu)質(zhì)量之間的平衡,本文以非凸算法中的FOCUSS算法為基礎(chǔ)進行改進,改進后算法運行時間僅耗時約2 s,PSNR值均在37 dB以上,具體從兩方面進行,一是根據(jù)頻域OCT圖像特點,結(jié)合分塊思想降低圖像的數(shù)據(jù)量,引入lp范數(shù)簡化運算過程,提高算法運行速度;二是從提高圖像重構(gòu)質(zhì)量進行考慮,對圖像進行分塊會造成嚴重的塊效應(yīng),因此在每次迭代求解過程中嵌入各向異性平滑濾波算子,抑制圖像塊效應(yīng),提高圖像質(zhì)量,實驗結(jié)果表明,壓縮感知技術(shù)可成功應(yīng)用于頻域OCT圖像,且改進的壓縮感知重構(gòu)算法在頻域OCT圖像重構(gòu)中能花費較短的時間獲得較高的重構(gòu)質(zhì)量。
傳統(tǒng)的信號獲取需滿足采樣頻率大于兩倍信號最高頻率的Nyquist采樣定理。而CS利用圖像(信號)的稀疏性,以遠低于傳統(tǒng)采樣方法獲取的數(shù)據(jù)量便可完成初始信號的完整恢復(fù),突破了傳統(tǒng)采樣定理的限制[4-8],因此本文將壓縮感知技術(shù)引入應(yīng)用到頻域OCT圖像的重構(gòu)上,下面將介紹CS理論的三個組成部分。
(1)稀疏表示
信號的稀疏表示是應(yīng)用壓縮感知的前提,假設(shè)φi是RN空間的一組向量基,構(gòu)成基矩陣Ψ=[φ1,φ2,...,φN],則對于任意的頻域OCT圖像x∈RN均可表示為:
(1)
式中:α是x在Ψ域的稀疏表示,若α中非零元數(shù)目比x少得多,則稱信號x是稀疏的。常見的稀疏變換基有DCT,DWT(Discrete Wavelet Transform),DDWT(Discrete Dyadic Wavelet Transform),Curvelet以及冗余字典基等,在下文仿真實驗中選擇DCT基獲得頻域OCT圖像的稀疏表示。
(2)觀測矩陣
壓縮感知理論指出,若采用一個與稀疏變換基Ψ不相關(guān)的測量矩陣Φ∈RM×N,對頻域OCT圖像x進行投影得到觀測值y,則可以通過求解優(yōu)化問題高概率精確地重構(gòu)出x,表達式為:
y=Φx=ΦΨα,
(2)
式中:A=ΦΨ為傳感矩陣,要求測量矩陣Φ和稀疏變換基Ψ不相關(guān),即在較高概率上滿足約束等距性條件,常見的觀測矩陣分為兩類,確定性矩陣和隨機矩陣,本文在仿真實驗中選擇第二類高斯隨機矩陣對頻域OCT圖像進行觀測。
(3)重構(gòu)算法
對式(2)求解可轉(zhuǎn)化為式(3):
(3)
求解上述l0范數(shù)是NP難解問題,需要列出α中所有非零項位置的線性組合才有可能得到精確解,計算過程復(fù)雜且不易實現(xiàn)。因此常將該問題轉(zhuǎn)換為其他替代模型進行求解。常用的重構(gòu)算法有三類,本文利用第三類非凸算法中的改進FOCUSS重構(gòu)算法進行求解,接下來將詳細闡述該改進算法的原理及實現(xiàn)。
利用CS可實現(xiàn)對傳統(tǒng)二維圖像的重構(gòu),但此方法是對整幅圖像進行重構(gòu),觀測矩陣所需的存儲空間大,重構(gòu)時間長。為了解決這個問題,本文提出分塊FOCUSS重構(gòu)算法,根據(jù)頻域OCT圖像特點,將圖像分塊理論和FOCUSS重構(gòu)算法相結(jié)合。具體做法是先將一幅大小為N×N的OCT圖像x分成n個B×B大小的圖像塊,設(shè)xi(i=1,2,...,n,n=(N/B)2)是第i個圖形塊的列向量形式,采用相同的高斯觀測矩陣ΦB,則對第i塊進行壓縮采樣可表示為:
yi=ΦB·xi.
(4)
本文中壓縮觀測矩陣ΦB大小為nB×B2。
傳統(tǒng)FOCUSS算法中由于偽逆的存在,往往導(dǎo)致能量分散,無法獲得精確的稀疏解[18]。將lp范數(shù)作為正則項引入到FOCUSS算法中,通過逐步迭代使解的能量局部化。對每一塊分別應(yīng)用FOCUSS算法,設(shè)置p=1,用以簡化運算過程,即求解以下優(yōu)化問題:
(5)
使用拉格朗日方法將上式中約束優(yōu)化問題轉(zhuǎn)化為非約束優(yōu)化問題,其中,拉格朗日函數(shù)被定義為:
L(xi,λ)=f(xi)+λT(ΦBxi-yi),
(6)
式中λ=(λ1,λ2,...,λmB)T構(gòu)成拉格朗日乘數(shù)。
(7)
(8)
為了找到目標(biāo)函數(shù)的梯度,將矩陣Ψ分解為其行向量φi得到:
(9)
分解之后,目標(biāo)函數(shù)表示為:
(10)
目標(biāo)函數(shù)相對于的xi的梯度為:
(11)
目標(biāo)函數(shù)相對于元素xj的偏導(dǎo)數(shù)為:
(12)
使用因子形式表示目標(biāo)函數(shù)的梯度向量即:
經(jīng)過半年多的運轉(zhuǎn),浮選車間運行狀況穩(wěn)定。浮選精煤灰分控制在8%以內(nèi),三班工作制,每班精煤壓濾機可卸3~4個循環(huán),按每個循環(huán)8 t計,則每天可生產(chǎn)浮選精煤80 t左右。每月可多洗精煤約2 500 t,按當(dāng)時精煤價格1 300元/t測算,每月可增加效益近350萬元。經(jīng)測算,實際運行4個月后,已收回投資成本。
(13)
式中Π(xψ)=diag(|xψ|-1)是大小為n×n的對角矩陣且xψ=Ψxi。
參考式(8)和式(13),得出:
(14)
參考式(5),得出:
(15)
綜合式(7),式(14),式(15),穩(wěn)定點滿足如下兩個表達式:
(16)
(17)
(18)
將式(18)代入式(17),得出:
(19)
將式(20)代入式(19),得出:
(20)
定義Aψφ=(ΦBΨ-1),經(jīng)過簡單的轉(zhuǎn)換,將(20)寫成:
(21)
(22)
(23)
通過修改拉格朗日函數(shù)來調(diào)整松弛方法以找到迭代解,得到如下方程:
(24)
(25)
對所有的小圖像塊進行相同的處理,再把得到的圖像塊加以組合就完成了整幅頻域OCT圖像的重構(gòu)。本文的算法流程圖如圖1所示。
圖1 改進的分塊壓縮感知FOCUSS重構(gòu)算法流程圖
本文選擇兩種頻域OCT圖像進行仿真實驗,其中指尖圖像來自實驗室Thorlabs公司VEGA系列的SS-OCT成像系統(tǒng),光源中心波長為1 310 nm,穿透深度為12 mm,軸向分辨率為16 μm,掃頻速率為100 kHz,系統(tǒng)靈敏度為101 dB。眼底圖像來自杜克大學(xué)眼科實驗數(shù)據(jù)采集中心[19],成像儀器是德國海德堡公司的Spectralis SD-OCT系統(tǒng),光源中心波長為870 nm,穿透深度為1.8 mm??v向分辨率和橫向分辨率分別為3.8 μm和14 μm。實驗室采集到的SS-OCT指尖圖像為455×512 pixel的RGB圖像,圖像細節(jié)較少,整體較平坦。杜克大學(xué)數(shù)據(jù)庫中的SD-OCT眼底圖像為496×512 pixel的灰度圖像,圖像細節(jié)較豐富。為了方便對其進行分塊處理,均將圖像大小調(diào)整為512×512 pixel。
為了驗證本文算法能否有效重構(gòu)頻域OCT圖像,除了直接從人眼視覺效果評價重構(gòu)圖像質(zhì)量,同時采用峰值信噪比(Peak Signal to Noise Ratio,PSNR),用RPSNR表示,結(jié)構(gòu)相似性(Structural Similarity Index,SSIM),用MSSIM表示,取值范圍0-1,以及重構(gòu)時間t定量評價算法的重構(gòu)性能。RPSNR值越大,表明與完整采樣圖像相比,重構(gòu)圖像與之越接近,重構(gòu)效果越好,重構(gòu)精度越高;MSSIM值越大,表明重構(gòu)圖像保留了越多的原始完整采樣圖像中的信息,重構(gòu)圖像質(zhì)量越好;重構(gòu)時間t越小,算法的重構(gòu)效率越高。
本文中所有實驗在64位Windows10系統(tǒng)下利用軟件MATLAB R2018a進行,采樣率(Sampling Rate,SR)是重構(gòu)過程中采樣值大小z與原始信號x大小的比值,首先對實驗室采集到的指尖OCT圖像進行壓縮感知重構(gòu)處理,利用DCT基對原始圖像稀疏表示,采用高斯隨機矩陣對圖像進行觀測,并運用5種優(yōu)化算法對指尖圖像進行恢復(fù)重構(gòu)。圖2給出SR為0.3,分塊大小為8×8時原始FOCUSS重構(gòu)算法、當(dāng)前的主流分塊重構(gòu)算法-平滑投影算法(SPL)以及改進FOCUSS重構(gòu)算法的恢復(fù)效果以及局部細節(jié)對比圖。
據(jù)圖2可知,當(dāng)圖像采樣率SR為0.3,分塊大小為8×8時,改進FOCUSS重構(gòu)算法對OCT指尖圖像進行重構(gòu)后能夠獲得更好的視覺重構(gòu)效果,圖像細節(jié)保持完整,圖像邊緣處容易分辨,與原圖更接近。圖2(b)表明直接運行原始FOCUSS重構(gòu)算法得到的重構(gòu)圖像有一定程度的信息缺失,圖像的邊緣處不光滑;從圖2(c)~圖2(e)可以看出3種主流的圖像壓縮感知分塊算法對圖像進行重構(gòu)后,重構(gòu)圖像的細節(jié)存在明顯的邊界效應(yīng);通過圖2(f)可以得出改進的FOCUSS算法顯著改善了圖像塊效應(yīng),獲得的重構(gòu)圖像細節(jié)清晰,輪廓明顯,滿足人眼視覺要求,成功將壓縮感知技術(shù)應(yīng)用到掃頻OCT系統(tǒng)產(chǎn)生的指尖圖像上,本文算法重構(gòu)細節(jié)較少,畫面平坦的彩色圖像效果較好。
為了便于定量比較實驗結(jié)果,將采樣率固定為0.3,重復(fù)試驗10次,取平均值作為最終仿真結(jié)果值,表1給出原始FOCUSS算法、使用3種不同稀疏變換基(DCT,DWT,DDWT)的分塊平滑投影重構(gòu)算法以及改進FOCUSS算法對指尖OCT圖像的壓縮感知重構(gòu)客觀評價指標(biāo)值。
圖2 采樣率為0.3時不同算法重構(gòu)效果及局部細節(jié)圖
表1SR為0.3時各算法對指尖圖像的重構(gòu)結(jié)果
Tab.1 Reconstruction results of fingertip images by each algorithm when the SR is 0.3
算法PSNR/dBSSIMt/s原始FOCUSS35.340.878 778.65分塊SPL-DCT35.010.825 68.64分塊SPL-DWT35.660.853 113.79分塊SPL-DDWT36.040.865 710.72改進FOCUSS37.710.938 31.89
據(jù)表1可知,與原始FOCUSS算法相比,改進FOCUSS算法對圖像進行重構(gòu)后提高了圖像質(zhì)量,PSNR值約高出2 dB,SSIM值提高到0.938 3,并且對圖像進行分塊后大大降低了算法運行時間,算法運行僅花費1.89 s,提高了算法運行效率,算法實時性強;與同樣結(jié)合了分塊思想的SPL算法相比,改進FOCUSS算法不僅改善了圖像的塊效應(yīng),能獲得更高的PSNR值和SSIM值,重構(gòu)圖像質(zhì)量高,與原始完整采樣OCT指尖圖像更接近且保留了圖像中的大部分信息,而且算法運行花費時間更短,說明該算法計算復(fù)雜度更低,運算效率高,更易推廣使用。同時,為了說明本文算法的普遍性,另外進行了4組仿真實驗,將圖像分塊大小分別固定為8×8,16×16,32×32,64×64 4種情況,仿真獲取SR分別為0.1,0.2,0.3,0.4和0.5時5種算法對指尖圖像的重構(gòu)結(jié)果圖和客觀評價指標(biāo)值,實驗圖和數(shù)據(jù)均表明,當(dāng)圖像分塊大小為8×8,16×16,32×32時,采樣率大于0.2時,改進的壓縮感知重構(gòu)算法較以往算法能獲得重構(gòu)速度和重構(gòu)精度的平衡,是具有一定優(yōu)勢的改進算法。
上述仿真采用的指尖圖像細節(jié)較少,為了進一步說明本文算法的有效性,同時具體分析采樣率與OCT圖像重構(gòu)質(zhì)量的關(guān)系,隨機選取一張眼底OCT切片,切片圖像的細節(jié)較豐富,將采樣率分別設(shè)置為0.1,0.2,0.3,0.4和0.5,圖3~圖7給出了當(dāng)采樣率為0.1,0.3和0.5時各算法對完整采樣眼底OCT圖像的重構(gòu)結(jié)果。據(jù)五組圖可以看出,隨著采樣率的增加,重構(gòu)圖像的細節(jié)、輪廓越來越清晰,視覺效果越來越好,與原始完整采樣圖像越來越接近。據(jù)圖3~圖7(b)可知當(dāng)采樣率為0.1時,各算法都不能很好地重構(gòu)原始完整采樣圖像,圖像信息缺失嚴重,圖像邊緣位置無法確定,圖像視覺效果差;據(jù)圖3~圖7(c)和圖7(d)可知,隨著采樣率增加至0.3以上時,重構(gòu)后的圖像更加清晰,邊緣細節(jié)保持完整,邊緣處較為光滑,圖像輪廓基本保持良好,能獲得較高的圖像質(zhì)量,容易分辨各視網(wǎng)膜層,與原始完整采樣圖像更接近,從視覺效果看并無明顯差別。所以采用壓縮感知技術(shù)可以用較少的采樣數(shù)據(jù)量精確重構(gòu)出譜域OCT系統(tǒng)產(chǎn)生的原始眼底OCT圖像,本文算法對細節(jié)豐富的灰度圖像的重構(gòu)效果也較好,適用范圍較大。
圖3 原始FOCUSS算法重構(gòu)結(jié)果
圖4 分塊SPL-DCT算法重構(gòu)結(jié)果
圖5 分塊SPL-DWT算法重構(gòu)結(jié)果
圖6 分塊SPL-DDWT算法重構(gòu)結(jié)果
圖7 改進FOCUSS算法重構(gòu)結(jié)果
為了直觀反映重構(gòu)算法的重構(gòu)性能和采樣率的關(guān)系,圖8給出了不同采樣率下各算法對眼底OCT切片圖像重構(gòu)后的客觀評價指標(biāo)值變化曲線圖??梢钥闯鲭S著采樣率的增加,各算法對圖像的重構(gòu)效果越來越好,重構(gòu)質(zhì)量越來越高,重構(gòu)時間并沒有大幅度變化,重構(gòu)性能越來越好。據(jù)圖8(a)可知,在采樣率為0.1和0.2時,各算法對圖像重構(gòu)后的PSNR值均在36 dB以下,重構(gòu)質(zhì)量不高,但是隨著采樣率增加至0.3以上時,本文改進的結(jié)合了分塊思想的FOCUSS重構(gòu)算法的PSNR值迅速高于其他4種對比算法,圖像重構(gòu)精度高;從圖8(b)可以看出改進的FOCUSS算法的SSIM值始終高于其他算法,對圖像進行重構(gòu)后保留了圖像中的大部分有用信息,圖像重構(gòu)質(zhì)量高。圖8(c)則可以看出原始FOCUSS重構(gòu)算法的重構(gòu)時間均在75 s以上,算法重構(gòu)效率不高,而改進的算法則將重構(gòu)時間控制在2 s以下,算法運行時間也遠遠低于3種SPL算法,優(yōu)勢顯著,性能良好,綜合來看,改進重構(gòu)算法能同時獲得較高的圖像重構(gòu)質(zhì)量和重構(gòu)效率。
本文改進的壓縮感知FOCUSS重構(gòu)算法結(jié)合了圖像分塊思想,為了進一步分析分塊大小對頻域OCT圖像重構(gòu)質(zhì)量和重構(gòu)效率的影響,本文將采樣率固定為0.5,將圖像塊大小分別設(shè)置為8×8,16×16,32×32,64×64對眼底圖像進行改進壓縮感知重構(gòu)算法的仿真實驗,重復(fù)10次,取平均值,最終結(jié)果如表2所示。
圖8 5種算法在不同采樣率下的PSNR值、SSIM值以及重構(gòu)時間
表2 改進FOCUSS算法在不同分塊大小下的仿真結(jié)果
Tab.2 Simulation results of improved FOCUSS algorithm under different block sizes
PSNR/dBSSIMt/s8×842.470.964 11.9616×1642.290.961 52.1632×3241.840.956 84.6164×6440.990.946 553.45
據(jù)表2可知,隨著OCT圖像分塊大小的增加,算法的運行時間加長,重構(gòu)效率降低,算法實時性變差,這是由于高斯隨機觀測矩陣對各個塊圖像進行測量是同時進行的,前一個圖像塊對后一個圖像塊沒有依賴性,圖像分塊越大,觀測矩陣越大,系統(tǒng)的存儲代價越大,算法復(fù)雜度越高。從PSNR值和SSIM值的變小可以看出圖像重構(gòu)后圖像質(zhì)量有小幅度降低,這是由于本文處理的OCT圖像并不是標(biāo)準測試圖像,自身不可避免含有散斑噪聲,并且在算法迭代過程中結(jié)合了各向異性平滑濾波算法,圖像塊尺寸越小,邊緣平滑得越好,散斑也得到抑制。綜合來看,當(dāng)圖像分塊大小為8×8時,重構(gòu)圖像的PSNR值高達42.47 dB,SSIM值高達0.964 1,接近于1,算法運行僅耗時1.96 s。所以綜合來看,頻域OCT圖像分塊大小為8×8時,在圖像進行重構(gòu)后可以同時獲得較高的重構(gòu)質(zhì)量和重構(gòu)效率。
本文首先介紹了壓縮感知的基本框架,通過分析頻域OCT圖像特點,成功將壓縮感知技術(shù)引入到頻域OCT圖像稀疏重構(gòu)上。通過對細節(jié)豐富程度不同的兩種頻域OCT圖像在不同采樣率下及不同分塊大小的仿真實驗得出,本文提出的改進FOCUSS重構(gòu)算法通過結(jié)合分塊思想、引入正則項lp范數(shù)以及在算法迭代求解過程中嵌入各向異性濾波算子,不僅降低了計算復(fù)雜度,同時有效抑制了重構(gòu)圖像的塊效應(yīng),提高了圖像質(zhì)量,對OCT圖像進行壓縮感知重構(gòu)后能獲得更好的視覺效果以及更短的計算時間,本文算法適應(yīng)范圍較廣。通過進一步仿真得到當(dāng)圖像分塊大小為8×8時能同時獲得較高的重構(gòu)質(zhì)量和重構(gòu)速度,PSNR值高達42.47 dB,遠高于文獻14中的30 dB,SSIM值可達到0.964 1,重構(gòu)時間僅花費1.96 s,比文獻12中的運行時間16.6 s要快得多,本文改進的壓縮感知重構(gòu)算法在一定程度上實現(xiàn)了重構(gòu)效率和重構(gòu)精確度的平衡,具有一定的優(yōu)勢。