姚佳麗　 李中源　 吳華珍　 李　光　 羅守華

(東南大學生物科學與醫(yī)學工程學院, 南京 210096)

?

基于字典學習的超分辨率顯微CT圖像重建

姚佳麗 李中源 吳華珍 李光 羅守華

(東南大學生物科學與醫(yī)學工程學院, 南京 210096)

為提高顯微CT重建圖像的空間分辨率,提出了一種基于字典學習的超分辨率圖像重建算法．首先,將重建圖像進行網(wǎng)格細化,并使用面積權(quán)值模型實現(xiàn)對投影過程的精確建模．然后,選擇高質(zhì)量的圖像作為訓練樣本,采用K-SVD算法構(gòu)建圖像字典．基于該圖像字典,利用正交匹配追蹤算法實現(xiàn)對重建圖像的稀疏表達,并以此作為稀疏項約束引入到重建算法的目標函數(shù)中．最后,使用梯度下降法求解目標函數(shù)．實驗結(jié)果表明:與傳統(tǒng)的基于插值的超分辨率重建算法相比,所提算法的超分辨率結(jié)果在圖像對比度、邊緣保持方面具有優(yōu)勢,并且保留了更多的圖像高頻信息,從而有效提高了重建圖像的空間分辨率．

超分辨率重建;字典學習;面積權(quán)值;微計算機斷層掃描技術(shù)

微計算機斷層掃描技術(shù)(又稱顯微CT)是一種非破壞性、非介入式、高分辨率的成像技術(shù),能夠在不破壞樣本的情況下了解樣本的內(nèi)部顯微結(jié)構(gòu),從而被廣泛應(yīng)用于醫(yī)學、藥學、材料學、工業(yè)無損探測等領(lǐng)域．在保證X射線安全劑量的前提下,重建圖像的空間分辨率是各種應(yīng)用中所關(guān)注的最重要的指標之一,它的提高對于顯微CT在各領(lǐng)域中的進一步應(yīng)用具有顯著意義[1]．然而,在實際成像過程中,重建圖像空間分辨率的提高往往受限于設(shè)備的物理條件和實際應(yīng)用要求,如X射線輻射劑量、掃描時間、探測器配置以及掃描的幾何結(jié)構(gòu)等．因此,提高顯微CT成像的空間分辨率成為研究的熱點和難點．

目前的解決方案可分為以下3類[2]:① 通過縮短探測器孔隙寬度及探測器單元間距,提高探測器的分辨率;但該方案會受到物理硬件條件的限制[3],若探測器單元大小超過衍射極限,還需要考慮衍射造成的影響．② 通過局部成像,將感興趣區(qū)域放大,以提高投影的空間分辨率[3];但該方案會導(dǎo)致投影數(shù)據(jù)橫向截斷,理論上不能實現(xiàn)精確重建,重建圖像往往存在杯狀偽影,從而導(dǎo)致圖像質(zhì)量下降．③ 改變掃描方式以提高采樣頻率,具體實現(xiàn)方式包括探測器中心偏置1/4像素掃描技術(shù)[4]、射線源焦點擺動法[5]及像素錯位技術(shù)[6]等;但該方案對機械精度要求較高,并且多次曝光會使X射線的輻射劑量加倍．

近年來,壓縮感知理論在圖像放大、去噪及恢復(fù)等領(lǐng)域得到越來越多的關(guān)注,在該理論基礎(chǔ)上建立的基于學習的超分辨率重建算法成為相關(guān)研究的熱點．Freeman等[7]提出了基于樣本圖像學習的超分辨率重建算法,充分利用自然圖像集,實現(xiàn)了基于學習的超分辨率重建算法的擴展和開發(fā); Wang等[8]提出了一種基于實例的超分辨率重建算法,通過對圖像塊之間稀疏性關(guān)聯(lián)的學習，重建了高分辨率圖像;Jeong等[9]利用K均值算法學習構(gòu)建了分類字典,并利用分類字典實現(xiàn)了單幅圖像的重建．上述算法通過對高、低分辨率訓練樣本的學習,建立了一種高低分辨率的對應(yīng)關(guān)系,相比于傳統(tǒng)的基于插值的超分辨率重建算法[10]可獲得更多的細節(jié)信息．然而,目前關(guān)于這類算法的研究主要集中于自然圖像．本文將字典學習的方法應(yīng)用到CT超分辨率圖像重建中．通過將重建圖像網(wǎng)格細化和使用精度更高的系數(shù)矩陣,采用字典學習的方法對重建圖像進行稀疏表達,并作為重建算法的稀疏項約束,利用直接獲得的投影數(shù)據(jù)進行超分辨率圖像重建．

1　超分辨重建模型的建立

圖像離散化后,CT圖像重建可以轉(zhuǎn)化為如下的代數(shù)方程:

(1)

式中,vj為重建圖像V的像素值,其中,j=1,2,…,N,N=Nwid×Nhei為重建圖像的像素點總數(shù),Nwid,Nhei分別為重建圖像的寬度和高度,即V∈RNwid×Nhei;投影數(shù)據(jù)矩陣P的元素pi表示第i條投影射線的投影值,其中,i=1,2,…,M,M=Mproj×Mbin為投影射線總數(shù),Mproj,Mbin分別為投影的角度數(shù)和探測器的寬度,即P∈RMproj×Mbin;系數(shù)矩陣W的元素wij為投影系數(shù),表示第j個像素對第i個投影值的貢獻．式(1)可寫成矩陣相乘的形式,即

WV=P

(2)

傳統(tǒng)采樣過程中,將重建圖像V離散化為Nwid×Nhei的離散圖像(見圖1(a)),每個像素的大小與探測器像元的大小相同．王麗艷等[11]從采樣定律和壓縮感知理論出發(fā),利用稀疏表示的先驗知識,證明可通過原有的低數(shù)據(jù)量測量進行高精度圖像重建．由此提出了一種超分辨重建模型,將重建圖像V進一步細分為ηNwid×ηNhei,其中η∈Z為超分辨因子,且η>1．超分辨因子為2時重建圖像的超分辨率網(wǎng)格的示意圖見圖1(b),重新離散化后像素的大小為探測器像元大小的1/η．

(a) 重建圖像網(wǎng)格

(b) 重建圖像的超分辨率網(wǎng)格

根據(jù)提出的重建模型,CT圖像重建轉(zhuǎn)化為新的代數(shù)方程:

AX=P

(3)

式中,X∈RηNwid×ηNhei為圖像V經(jīng)過超分辨率網(wǎng)格細化后的重建圖像;A為圖像V經(jīng)過超分辨率網(wǎng)格細化后的系數(shù)矩陣,其元素aij為投影系數(shù)．在壓縮感知理論下,假設(shè)重建圖像X經(jīng)過某種變換φ作用后稀疏化,利用原有的投影數(shù)據(jù)P,可以重建出空間分辨率提高了η2的圖像．重建過程為求解稀疏優(yōu)化的問題,即

s.t.AX=P

(4)

超分辨重建的目的在于盡可能保持圖像的高頻信息,即需要對投影過程進行精確建模并對系數(shù)矩陣精確求解,因此aij的獲取是重建算法中的關(guān)鍵．常用方法為長度權(quán)值法(見圖2(a)),即將射線看作沒有寬度的直線,將其與圖像像素相交的長度作為該像素對射線貢獻的權(quán)值(即圖2(a)中箭頭標注的線段長度);該方法簡單易行,但精度較低,難以滿足本文算法的要求．為此,本文采用如圖2(b)所示的面積權(quán)值法來確定權(quán)重,即將射線視作具有寬度的射束,將其與像素方塊重疊的面積作為貢獻的權(quán)值(即圖中陰影區(qū)域的面積)．與長度權(quán)值法相比,面積權(quán)值法雖然計算復(fù)雜,但更符合工程實際,重建圖像質(zhì)量更優(yōu)．

(a) 長度權(quán)值法

(b) 面積權(quán)值法

以等角扇束掃描方式為例,在二維平面上,當射線源旋轉(zhuǎn)至某一投影角度時,計算相鄰2條投影射線在像素邊界的交點坐標,以及被穿過的像素在整個重建圖像中的索引號,并將其保存到數(shù)組中;然后,從第1個交點坐標開始,根據(jù)相鄰射線在像素邊界的交點情況,確定射線束與像素格相交的多邊形形狀．圖3(a)為射線束穿過重建區(qū)域的示意圖．圖中,tim表示第i條投影射線與像素網(wǎng)格的交點,m=1,2,…,10為交點的索引號．

(a) 射線束穿過重建區(qū)域的示意圖

(b) 射線束與像素格相交產(chǎn)生多邊形組合的一般情況

相鄰射線可能穿過不同的像素格,產(chǎn)生不同的多邊形組合(見圖3(b))．由圖可知,情況1下,射線束穿過2個像素,與像素格相交的區(qū)域分別形成2個不規(guī)則多邊形;情況2下,射線束穿過1個像素,與像素格相交的區(qū)域形成平行四邊形;情況3和情況4下,射線束穿過2個像素,與像素格相交的區(qū)域分別形成三角形和不規(guī)則多邊形．根據(jù)交點坐標值計算面積大小,并將其作為該像素對該射線束的貢獻值,即權(quán)值aij．繼續(xù)在各自的交點數(shù)組中遍歷下一個交點坐標,若一個交點數(shù)組優(yōu)先遍歷結(jié)束,另一個交點數(shù)組需要繼續(xù)遍歷,直到最后一個交點．

除上述討論的一般情況外,對于一些特殊的個別情況,例如射線斜率為±1、相鄰射線從重建區(qū)域的不同邊界射入等,討論方法類似．

2　基于字典學習的方法

2.1字典表達

根據(jù)稀疏表示理論,信號可通過有限的字典原子的線性組合來表示．因此,圖像塊可表示為

(5)

給定字典后,需要找到表示圖像塊x的稀疏表達,即確保稀疏表示系數(shù)α中非零個數(shù)最少,等價于求解如下的優(yōu)化問題:

(6)

式中,ε≥0為誤差限．

引入拉格朗日常數(shù)δ,可將式(6)改寫成如下的非約束形式:

(7)

一般采用貪婪算法(如正交匹配追蹤算法[12])來求解式(6)或式(7)．即在迭代過程中,選擇與信號殘差最匹配的原子,將其正交投影到被選擇原子所構(gòu)成的空間中,重新計算殘差信息,直至殘差滿足預(yù)先設(shè)定的誤差限．

2.2字典學習

(8)

式(8)等價于求解以下任一目標函數(shù):

s.t.?s,‖αs‖0≤T0

(9)

(10)

式中,T0為稀疏度．

本文使用高劑量圖像作為訓練樣本,利用K-SVD算法[13]獲取合適的字典．該算法對字典D中的原子逐個進行更新,因此可將式(9)細化到每一個信號,即

s.t.‖αs‖≤T0

(11)

K-SVD算法的目的在于使得式(11)值最小,聯(lián)合式(9)和式(11)將目標函數(shù)修改為下式進行求解:

(12)

3　基于字典學習的圖像重建算法

CT圖像重建算法中的迭代重建算法抗噪性能較好,偽影抑制能力較強,并且在無完備投影數(shù)據(jù)的情況下仍能重建出質(zhì)量較好的圖像,因此得到了廣泛的應(yīng)用．另外,字典學習的方法不僅引入了先驗知識,減少了方程數(shù)量,還包含了圖像的結(jié)構(gòu)信息．將字典學習的方法與圖像的迭代重建算法相結(jié)合,便可有效實現(xiàn)超分辨率圖像重建．圖像重建的過程即為求解以下目標函數(shù):

(13)

式中,Es為從離散圖像X中提取圖像塊的矩陣函數(shù),且Es∈RK×N;λ為正則化系數(shù)．式(13)中第1項表示最小化計算投影與測量投影之間的差異,為保真項;第2項為使用字典進行圖像稀疏分解的先驗懲罰項．式(13)的目標函數(shù)可分解為如下2個函數(shù)交替求解:

1) 固定D和X求解αs,即

(14)

2) 根據(jù)字典D和αs更新X,即

(15)

式(15)為二次凸函數(shù)形式,可利用梯度下降法求解．

4　實驗

利用本文算法進行超分辨率圖像重建,以驗證超分辨模型和本文算法的有效性．投影數(shù)據(jù)由蘇州海斯菲德信息科技有限公司HiscanM-1000型小動物CT拍攝小鼠獲得．CT掃描管電壓為40kV,管電流為200μA．用于重建的投影數(shù)據(jù)大小為992×360,其中992表示探測器的寬度,360表示投影的角度數(shù)．原始重建圖像的大小為886×886像素,按照本文模型重建后的圖像大小為1 772×1 772像素．為獲得高質(zhì)量的圖像以作為訓練圖像字典的樣本,在360°范圍內(nèi)按圓形軌跡掃描方式均勻采集了1 800個角度的投影圖,并使用FBP算法重建得到高質(zhì)量的圖像．以8×8像素大小的圖像塊對該高質(zhì)量圖像進行采樣,得到訓練集,字典大小設(shè)定為64×256,其中64為圖像塊按行重排為列向量后的行數(shù),256為字典的原子個數(shù)．訓練時選取的稀疏度為3,迭代次數(shù)為20．訓練得到的圖像字典如圖4所示．

圖4　基于高質(zhì)量圖像訓練獲得的圖像字典

實驗中,設(shè)置圖像初始值X0=0．計算機配置為Inteli5-2400,3.10GHzCPU．為驗證本文算法的性能,將其與以下2種傳統(tǒng)的基于插值的超分辨率重建算法進行比較:① 雙線性插值后處理算法(簡稱為算法1),即利用基于面積權(quán)值的迭代重建算法得到原始圖像,再使用雙線性插值算法將圖像插值為1 772×1 772像素的圖像;② 投影數(shù)據(jù)一維徑向線性插值前處理算法(簡稱為算法2),即將投影數(shù)據(jù)在探測器寬度上進行線性插值算法處理,得到1 984×360大小的投影數(shù)據(jù),再使用FBP算法重建為1 772×1 772像素大小的圖像．

選取不同λ進行對比實驗后得出,當λ=1.0時本文算法的結(jié)果相對更好．CT圖像的數(shù)據(jù)位數(shù)為12,顯示窗寬為2 483,窗位為1 434．圖5所示為3種算法的重建結(jié)果．

圖5　3種算法下小動物CT的重建圖像

為更好地觀察各算法的重建結(jié)果,分別選取圖5中代表性區(qū)域Ⅰ,Ⅱ,Ⅲ進行放大,結(jié)果見圖6．圖6(a)為原始圖像直接放大后的結(jié)果,模糊程度較大,很難觀察到圖像的細節(jié)信息,且邊緣有鋸齒狀偽影．圖6(b)和(c)分別為基于算法1和算法2的重建圖像,具有一定的超分辨率效果,能夠觀察到被放大后的圖像細節(jié),且邊緣基本不存在鋸齒狀偽影;但圖像對比度不高,細節(jié)信息的層次感不強．圖6(d)為基于本文算法的重建圖像,與圖6(b)和(c)相比,在圖像對比度、細節(jié)層次的豐富性和邊緣保持方面更具有優(yōu)越性,超分辨率效果更好．

為客觀評價算法的性能,分別選取圖5中代表性區(qū)域Ⅰ,Ⅱ,Ⅲ來計算相鄰像素灰度方差(SMD),它是用于衡量圖像清晰程度的最常用的評價函數(shù)之一[14]．計算結(jié)果見表1．由表可見,基于本文算法的重建圖像的SMD參數(shù)明顯高于其他2種算法,表明所得的重建圖像清晰度最高．

(a) 原始圖像　(b) 算法1　(c) 算法2　(d) 本文算法

表1　重建圖像中3處代表性區(qū)域的相鄰像素灰度方差

此外,選取各重建圖像中多處具有代表性的邊緣,計算并繪制了調(diào)制傳遞函數(shù)(MTF)[15]曲線,結(jié)果見圖7．由圖可見,基于算法1和算法2的重建圖像的MTF曲線下降速度明顯快于采用本文算法時的情況,表明利用本文算法得到的重建圖像中高頻信息保留得最完整,算法的優(yōu)越性得到進一步驗證．

圖7　重建圖像的MTF曲線

下面利用由蘇州海斯菲德信息科技有限公司提供的Hiscan H-100型高分辨率CT拍攝木質(zhì)牙簽獲得的投影數(shù)據(jù),對3種重建算法的性能進行比較．圖8給出了3種算法的重建圖像．圖9為重建圖像中代表性區(qū)域Ⅰ,Ⅱ放大后的結(jié)果．由圖9中箭頭所指區(qū)域可見,采用本文算法時孔洞結(jié)構(gòu)保留得最清晰、完整．

綜上可知,本文算法較其他算法在圖像清晰度、細節(jié)的完整性和邊緣保持方面有比較明顯的優(yōu)勢．

圖8　高分辨率CT下3種算法的重建圖像

(a) 算法1　　　　(b) 算法2　　　　(c) 本文算法

5　結(jié)語

本文將壓縮感知理論和圖像迭代重建算法相結(jié)合,提出了一種基于字典學習的超分辨率重建算法．在無需改變CT的硬件系統(tǒng)、機械結(jié)構(gòu)以及不需要對投影數(shù)據(jù)插值的情況下,通過將重建圖像網(wǎng)格細化,使用字典學習加迭代的方法進行稀疏求解,并在迭代過程中將原始的長度權(quán)值轉(zhuǎn)換為精度更高的面積權(quán)值,從而獲得了分辨率更高的重建圖像．實驗結(jié)果表明,較之傳統(tǒng)的基于插值的超分辨率重建算法,本文算法在圖像高頻信息的保留、圖像細節(jié)的完整性和邊緣保持方面具有優(yōu)勢,從而有效提高了重建圖像的空間分辨率．本文算法是在空域上實現(xiàn)的,故而可以擴展到三維錐束模型的求解上,但由此帶來的大規(guī)模數(shù)據(jù)量問題和算法效率下降問題有待進一步研究．

References)

[1]Ritman E L. Vision 20/20: Increased image resolution versus reduced radiation exposure[J].MedicalPhysics, 2008, 35(6): 2502-2512.DOI:10.1118/1.2919112.

[2]Park S C, Park M K, Kang M G. Super-resolution image reconstruction: A technical overview[J].IEEESignalProcessingMagazine. 2003, 20(3): 21-36. DOI:10.1109/msp.2003.1203207.

[3]Ffohr T G, Stiemtoffer K, Suss C, et a1.Novel ultrahigh resolution data acquisition and image reconstruction for multi-detector row CT[J].MedicalPhysics, 2007, 34(5): 1712-1723. DOI:10.1118/1.2722872.

[4]魏東波, 傅健, 李斌, 等. 提高三維ICT成像空間分辨率的掃描方式及其重建算法[J]. 航空動力學報, 2007, 22(5): 768-772. DOI:10.3969/j.issn.1000-8055.2007.05.014.

Wei Dongbo, Fu Jian, Li Bin, et al. Scan mode and its reconstruction algorithm for improving spatial resolution of 3D-ICT[J].JournalofAerospacePower, 2007, 22(5): 768-772. DOI:10.3969/j.issn.1000-8055.2007.05.014. (in Chinese)

[5]Lonn A H R. Computed tomography system with translatable focal spot:USA, 5173852[P]. 1992-12-22.

[6]傅健, 譚忍泊, 趙峰. 一種用于CT空間分辨率增強的圖像重建方法[J]. 航空動力學報, 2013, 28(5): 971-976.

Fu Jian, Tan Renbo, Zhao Feng. An image reconstruction method for improving CT spatial resolution[J].JournalofAerospacePower, 2013, 28(5): 971-976. (in Chinese)

[7]Freeman W T, Jones T R, Pasztor E C. Example-based super-resolution[J].IEEEComputerGraphicsandApplications, 2002, 22(2):56-65. DOI:10.1109/38.988747.

[8]Wang J, Zhu S, Gong Y. Resolution enhancement based on learning the sparse association of image patches[J].PatternRecognitionLetters, 2010, 31(1):1-10. DOI:10.1016/j.patrec.2009.09.004.

[9]Jeong S C, Song B C. Fast super-resolution algorithm based on dictionary size reduction usingk-means clustering[J].ETRIJournal, 2010, 32(4): 596-602.

[10]Gilman A, Bailey D G, Marsland S R. Interpolation models for image super-resolution[C]//The4thIEEEInternationalSymposiumonElectronicDesign,Test&Applications. Hong Kong,China, 2008: 55-60. DOI:10.1109/delta.2008.104.

[11]王麗艷, 韋志輝, 羅守華, 等. 壓縮感知在Micro-CT圖像超分辨重建中的應(yīng)用[J]. 中國圖象圖形學報, 2012, 17(4): 487-493.

Wang Liyan, Wei Zhihui, Luo Shouhua, et al. Image super reconstruction for micro-CT based on compressed sensing[J].JournalofImageandGraphics, 2012, 17(4): 487-493. (in Chinese)

[12]Donoho D L, Tsaig Y, Drori I, et al. Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[J].IEEETransactionsonInformationTheory, 2012, 58(2): 1094-1121.

[13]Aharon M, Elad M, Bruckstein A. K-SVD: An algorithm for designing over-complete dictionaries for sparse representation[J].IEEETransactionsonSignalProcessing, 2006, 54(11): 4311-4322. DOI:10.1109/tsp.2006.881199.

[14]Jarvis R A. Focus optimization criteria for computer image processing[J].Microscope, 1976, 24(2): 163-180.

[15]American Society for Testing and Materials(ASTM). E1441-00(2005) Standard guide for computed tomography(CT) imaging[S]. West Conshohocker, Pennsylvania,USA: ASTM, 2005.

Super-resolution image reconstruction for micro-CT based on dictionary learning

Yao Jiali Li Zhongyuan Wu Huazhen Li Guang Luo Shouhua

(School of Biological Sciences and Medical Engineering, Southeast University, Nanjing 210096, China)

To improve the spatial resolution of reconstructed images for micro computed tomography(micro-CT), a super-resolution image reconstruction algorithm based on dictionary learning is proposed. First, the reconstructed image grid is refined and the area weight model is used to achieve accurate modeling of the projection process. Then, high quality images are selected as the training samples. The K-means singular value decomposition (K-SVD) algorithm is adopted to train the image dictionary. On the basis of the image dictionary, the orthogonal matching pursuit algorithm is used to implement sparse representation of the reconstructed image, which is introduced into the objective function of the reconstruction algorithm as a sparse constraint. Finally, the gradient descent method is adopted to solve the objective function. The experimental results show that compared to the conventional interpolation-based super-resolution reconstruction algorithms, the proposed algorithm has advantages on the image contrast and the edge preservation, and retains more high frequency information of images, effectively improving the spatial resolution of the reconstructed images.

super-resolution reconstruction; dictionary learning; area weight; micro computed tomography(micro-CT)

10.3969/j.issn.1001-0505.2016.05.010

2016-01-06.作者簡介: 姚佳麗(1991—)，女，碩士生；羅守華(聯(lián)系人)，男，博士，副教授，luoshouhua@seu.edu.cn.

國家自然科學基金資助項目(61127002，61179035).

TP391.1

A

1001-0505(2016)05-0957-07

引用本文: 姚佳麗,李中源,吳華珍,等.基于字典學習的超分辨率顯微CT圖像重建[J].東南大學學報(自然科學版),2016,46(5):957-963. DOI:10.3969/j.issn.1001-0505.2016.05.010.