基于SIDWT和迭代自一致性的快速并行成像重建方法

段繼忠, 錢青青

(昆明理工大學(xué) 信息工程與自動(dòng)化學(xué)院,昆明 650504)

磁共振成像(Magnetic Resonance Imaging, MRI)利用人體組織中氫原子核的核磁共振現(xiàn)象,通過接收線圈采集k空間(即頻域)數(shù)據(jù),然后對(duì)采集到的k空間數(shù)據(jù)進(jìn)行傅里葉逆變換,從而重建出人體器官圖像[1].MRI具有無電離輻射、多角度成像、對(duì)人體組織無損傷等諸多優(yōu)點(diǎn),是臨床醫(yī)學(xué)和醫(yī)學(xué)科研中非常重要的檢測(cè)手段[2].

然而,MRI掃描時(shí)間較長(zhǎng),病人在掃描過程中易感到不適或發(fā)生移動(dòng)導(dǎo)致圖像偽影,影響病人的就醫(yī)體驗(yàn)和掃描圖像質(zhì)量.為加快MRI的掃描速度,通常對(duì)k空間數(shù)據(jù)進(jìn)行欠采樣后再使用算法重建出圖像,但欠采樣易使圖像質(zhì)量下降,影響醫(yī)生對(duì)患者病情的診斷.在算法模型中加入不同正則項(xiàng)可改善圖像重建質(zhì)量,但會(huì)導(dǎo)致重建速度較慢從而出診斷結(jié)果較慢,最終影響病人的就醫(yī)滿意度,限制了算法的推廣與應(yīng)用[3].因此,在保證MRI重建質(zhì)量的同時(shí),研究MRI的快速重建方法是非常實(shí)際且有意義的.

并行成像技術(shù)[4]和壓縮感知(Compressed Sensing, CS)[5]是減少M(fèi)RI掃描時(shí)間的兩種有效方法.其中,并行磁共振成像(Parallel MRI, PMRI)方法使用具有不同空間靈敏度信息的多通道接收線圈陣列同時(shí)采集磁共振信號(hào),利用通道間的關(guān)聯(lián)信息即可重建出未采集的信號(hào).而CS理論突破了奈奎斯特采樣定理,利用圖像的可壓縮性和稀疏性即可從很少的觀測(cè)值中重建出高質(zhì)量圖像.為進(jìn)一步提高M(jìn)RI的重建性能,通常將PMRI方法與CS理論結(jié)合應(yīng)用[6].

在PMRI技術(shù)中,Lustig等[7]提出一種針對(duì)任意k空間軌跡的并行成像重建方法,即迭代自一致性并行成像重建(Iterative Self-Consistent Parallel Imaging Reconstruction, SPIRiT)模型.該方法基于k空間數(shù)據(jù)自身的一致性,包括校正一致性和數(shù)據(jù)采集一致性,利用相鄰k空間點(diǎn)之間的相關(guān)性恢復(fù)丟失的信息,已被廣泛應(yīng)用于臨床實(shí)踐中[8].SPIRiT是一種基于自一致性的廣義PMRI重建框架,可以方便地與各種正則化方法結(jié)合,以實(shí)現(xiàn)更高性能的MRI重建.因此,許多研究者引入不同正則項(xiàng)對(duì)SPIRiT模型進(jìn)行一系列的研究和改進(jìn)[9-17].其中,段繼忠等[9]將聯(lián)合稀疏性正則項(xiàng)引入SPIRiT模型,提出一種高效的重建方法.之后,Duan等[12]又在SPIRiT模型中引入聯(lián)合全變分和聯(lián)合L1范數(shù)正則項(xiàng)約束,提出另一種基于SPIRiT的快速重建方法.這兩種方法都將復(fù)雜優(yōu)化問題進(jìn)行簡(jiǎn)化,然后利用算子分裂技術(shù)將其分解為易于計(jì)算或求解的子問題,最后使用快速迭代軟閾值算法(Fast Iterative Shrinkage-Thresholding Algorithm, FISTA)[18]進(jìn)行加速,縮短圖像的重建時(shí)間.最近,Zhang等[15]將L1范數(shù)正則項(xiàng)加入SPIRiT模型中,提出一種基于快速投影迭代軟閾值算法(projected FISTA, pFISTA)[19]的快速重建方法(pFISTA-SPIRiT).該方法允許在并行MR圖像重建中使用不同的緊框架,是一個(gè)廣義模型,已被用于人腦成像、肝臟動(dòng)態(tài)增強(qiáng)成像和精準(zhǔn)的腦腫瘤血管通透性成像等[20].

但pFISTA-SPIRiT方法對(duì)校正一致項(xiàng)與數(shù)據(jù)一致項(xiàng)都進(jìn)行梯度計(jì)算,導(dǎo)致收斂速度較慢.將重建模型中的復(fù)雜優(yōu)化問題進(jìn)行簡(jiǎn)化可以有效提升算法的收斂速度[9,12].因此,本文基于pFISTA-SPIRiT模型,提出一種基于平移不變離散小波變換(Shift-Invariant Discrete Wavelets Transform, SIDWT)和SPIRiT的快速并行磁共振成像重建方法(fast SIDWT-SPIRiT, fSIDWT-SPIRiT).與pFISTA-SPIRiT利用圖像域SPIRiT模型進(jìn)行重建不同,該方法利用頻域SPIRiT模型進(jìn)行重建.使用SIDWT和SPIRiT模型,將校正一致項(xiàng)和數(shù)據(jù)一致項(xiàng)合并為一項(xiàng)后再使用pFISTA技術(shù)進(jìn)行求解.實(shí)驗(yàn)結(jié)果表明,與其他兩種方法相比,該方法能夠有效提升算法的收斂速度,從而減少圖像的重建時(shí)間,并且重建質(zhì)量不變.

1 SPIRiT模型

SPIRiT是一種逐線圈、自動(dòng)校準(zhǔn)的并行MRI重建模型.設(shè)r為k空間數(shù)據(jù)點(diǎn)的位置,Rr為從k空間中選擇需要的數(shù)據(jù)點(diǎn)的算子;xj為第j個(gè)線圈上的全部k空間數(shù)據(jù),j=1,2,…,J,J為總線圈數(shù);Rrxj為以第j個(gè)線圈位置r為中心從所有線圈提取的k空間鄰域數(shù)據(jù)點(diǎn);xi(r)為第i個(gè)線圈上位置r處的k空間數(shù)據(jù),則xi(r)的重建如下:

(1)

(2)

式中:A為ACS的位置集合.式(1)是一組耦合線性方程,故可以將整個(gè)耦合方程組寫成矩陣形式,則所有線圈的校正一致性如下:

x=Gx

(3)

式中:x為從所有線圈中獲取的全部k空間數(shù)據(jù);G為在適當(dāng)位置包含gji的矩陣.

當(dāng)然,除了考慮校正一致性,還應(yīng)滿足數(shù)據(jù)采集一致性.數(shù)據(jù)采集一致性如下:

y=Dx

(4)

式中:y為從所有線圈中獲取的欠采樣數(shù)據(jù);D為從全k空間數(shù)據(jù)進(jìn)行欠采樣的矩陣.

2 fSIDWT-SPIRiT

在重建模型中加入正則項(xiàng)可以有效提升MRI的重建質(zhì)量,因此在SPIRiT模型中引入L1范數(shù)正則項(xiàng),則帶L1范數(shù)正則項(xiàng)的SPIRiT并行MRI重建模型可以表示為以下最優(yōu)化問題:

(5)

式中:Ψ為逐線圈小波變換,用于將圖像稀疏化,是一個(gè)緊框架.本文選用SIDWT作為實(shí)驗(yàn)中的緊框架.F為逐線圈傅里葉變換,F-1為逐線圈傅里葉逆變換.

使用罰函數(shù)技術(shù),可以得到式(5)的無約束版本:

(6)

(7)

針對(duì)式(7)中無約束模型的求解,令α=ΨF-1x,并且有Ψ*Ψ=I,則

(8)

(9)

式中:k為迭代次數(shù);L為f(αk)梯度的Lipschitz常數(shù).于是,易得式(9)的解為

(10)

(11)

式中:β為輸入矩陣;|β|為β的模.由于Ψ*Ψ=I故式(10)可重寫為

(12)

(13)

最后,利用FISTA[18]技術(shù)進(jìn)行加速,則整個(gè)求解過程如下:

(14)

(15)

(16)

(17)

綜上所述,得到基于SIDWT和SPIRiT模型的快速并行磁共振成像重建方法fSIDWT-SPIRiT,實(shí)現(xiàn)步驟如下.

算法一: fSIDWT-SPIRiT

1: Setx0=0,z0=0,θ0=0,k=0

2: Repeat

3:k=k+1

7: 直到達(dá)到最大迭代次數(shù),否則返回步驟3

3 實(shí)驗(yàn)仿真及分析

3.1 實(shí)驗(yàn)設(shè)置

所有實(shí)驗(yàn)均在計(jì)算機(jī)上進(jìn)行,計(jì)算機(jī)的配置為:Intel(R) Core(TM) i5-7200U CPU@2.50 GHz處理器、16 GB內(nèi)存和Windows 10操作系統(tǒng)(64位),所有方法均使用MATLAB實(shí)現(xiàn).

在實(shí)驗(yàn)中,將fSIDWT-SPIRiT的重建性能與pFISTA-SPIRiT和SIDWT-SPIRiT進(jìn)行對(duì)比.后兩種方法都基于SIDWT和SPIRiT模型,直接使用pFISTA技術(shù)進(jìn)行求解得到,區(qū)別在于pFISTA-SPIRiT在圖像域SPIRiT進(jìn)行重建,而SIDWT-SPIRiT在k空間域SPIRiT進(jìn)行重建.使用信噪比(Signal Noise Ratio, SNR)、結(jié)構(gòu)相似性(Structural SIMilarity, SSIM)指標(biāo)和高頻誤差范數(shù)(High-Frequency Error Norm, HFEN)3個(gè)評(píng)價(jià)指標(biāo)來衡量圖像的重建質(zhì)量.3個(gè)評(píng)價(jià)指標(biāo)的定義分別如下:

(18)

(19)

(20)

式中:VHFEN為HFEN值;filter(·)是一個(gè)拉普拉斯高斯濾波器,用于捕捉圖像邊緣,濾波核的大小為15×15像素,標(biāo)準(zhǔn)差為1.5像素.

VSNR和VSSIM的數(shù)值越高,VHFEN的數(shù)值越低,說明圖像的重建質(zhì)量越好.這3個(gè)評(píng)價(jià)指標(biāo)均在感興趣區(qū)域內(nèi)進(jìn)行計(jì)算,實(shí)驗(yàn)中手動(dòng)調(diào)整所有方法的參數(shù)使得VSNR值達(dá)到最優(yōu).

為驗(yàn)證fSIDWT-SPIRiT的有效性,在不同人體器官數(shù)據(jù)集上對(duì)各種方法的重建性能進(jìn)行比較,分別選取活體人類腦部切片圖GE_human_brain[21]、train_brain_AXT1POST_200_6001959[22-23]、心臟切片圖data_v1_k1[24]以及人類膝蓋切片圖train_knee_file1000005[22-23],并依次命名為數(shù)據(jù)集1,數(shù)據(jù)集2,數(shù)據(jù)集3和數(shù)據(jù)集4.其中,數(shù)據(jù)集1是尺寸為256×256像素的8通道腦部切片數(shù)據(jù);數(shù)據(jù)集2是使用20通道線圈獲取的腦部數(shù)據(jù)集的第1幀,然后使用線圈壓縮技術(shù)將其壓縮為8個(gè)虛擬線圈,尺寸為320像素×320像素;數(shù)據(jù)集3是使用28通道線圈獲取的心臟數(shù)據(jù)集,然后使用線圈壓縮技術(shù)將其壓縮為12個(gè)虛擬線圈,尺寸為192像素×192像素;數(shù)據(jù)集4是在15通道線圈獲取的膝蓋數(shù)據(jù)集的第20個(gè)切片,然后使用線圈壓縮技術(shù)將其壓縮為8個(gè)虛擬線圈,尺寸為320像素×320像素.

實(shí)驗(yàn)中的所有數(shù)據(jù)集都使用具有不同加速因子R(不包括ACS)的二維泊松圓盤采樣模式進(jìn)行欠采樣,并且所有方法都使用大小為24像素×24像素的校準(zhǔn)區(qū)域和5像素×5像素的SPIRiT核.本文主要從重建質(zhì)量和重建速度兩個(gè)方面來對(duì)比各種方法的優(yōu)劣.在重建質(zhì)量方面,主要通過重建圖和誤差圖從主觀上比較各種方法的重建性能,并采用3個(gè)評(píng)價(jià)指標(biāo)從客觀上比較各種方法的重建性能;在重建速度方面,主要通過時(shí)間來比較各種方法的重建性能.

3.2 不同方法的重建質(zhì)量比較

3.2.1不同重建方法的視覺比較 首先,從視覺上比較各種方法的重建性能.當(dāng)加速因子R=5(即欠采樣率為20%)時(shí),用3種方法對(duì)4個(gè)數(shù)據(jù)集進(jìn)行重建.圖1～4分別給出4個(gè)數(shù)據(jù)集對(duì)應(yīng)的全采樣圖像、二維泊松圓盤欠采樣掩膜、重建圖及誤差圖.

圖1 在5倍加速時(shí)不同方法對(duì)數(shù)據(jù)集1進(jìn)行重建的視覺比較Fig.1 Comparison of visual reconstruction for dataset 1 by different methods at 5 times acceleration

圖2 在5倍加速時(shí)不同方法對(duì)數(shù)據(jù)集2進(jìn)行重建的視覺比較Fig.2 Comparison of visual reconstruction for dataset 2 by different methods at 5 times acceleration

由圖1～4可以看出,對(duì)于4個(gè)不同人體器官的數(shù)據(jù)集,fSIDWT-SPIRiT與其余兩種方法相比,重建圖和誤差圖都無明顯差異.其中,fSIDWT-SPIRiT在不同數(shù)據(jù)集上的重建圖都有較清楚的紋理細(xì)節(jié),并保留了較好的邊緣輪廓信息,表明所提方法能夠?qū)崿F(xiàn)較好的重建,與其余方法的重建質(zhì)量相當(dāng).

3.2.2不同重建方法的評(píng)價(jià)指標(biāo)比較 從3個(gè)評(píng)價(jià)指標(biāo)上比較各種方法的重建性能.在加速因子R為3～7的情況下,即欠采樣率分別為33.3%、25%、20%、16.7%和14.3%時(shí),用3種方法對(duì)4個(gè)不同數(shù)據(jù)集進(jìn)行重建.表1～4分別給出4個(gè)不同數(shù)據(jù)集下不同方法對(duì)應(yīng)的VSNR值、VSSIM值和VHFEN值.

表1 在3～7倍加速時(shí)不同方法對(duì)數(shù)據(jù)集1進(jìn)行重建的數(shù)值比較Tab.1 Comparison of the values of the reconstruction for dataset 1 by different methods at 3 to 7 times acceleration

表2 在3～7倍加速時(shí)不同方法對(duì)數(shù)據(jù)集2進(jìn)行重建的數(shù)值比較Tab.2 Comparison of the values of the reconstruction for dataset 2 by different methods at 3 to 7 times acceleration

表3 在3～7倍加速時(shí)不同方法對(duì)數(shù)據(jù)集3進(jìn)行重建的數(shù)值比較Tab.3 Comparison of the values of the reconstruction for dataset 3 by different methods at 3 to 7 times acceleration

表4 在3～7倍加速時(shí)不同方法對(duì)數(shù)據(jù)集4進(jìn)行重建的數(shù)值比較Tab.4 Comparison of the values of the reconstruction for dataset 4 by different methods at 3 to 7 times acceleration

從表1～4可以看出,對(duì)于數(shù)據(jù)集1和數(shù)據(jù)集3,與pFISTA-SPIRiT相比,fSIDWT-SPIRiT的重建圖像VSNR值和VSSIM值更高,VHFEN值更低,說明所提方法的重建質(zhì)量更好;而對(duì)于數(shù)據(jù)集2和數(shù)據(jù)集4,pFISTA-SPIRiT和fSIDWT-SPIRiT的重建圖像VSNR值、VSSIM值和VHFEN值相差并不大,說明兩種方法的重建質(zhì)量無顯著差異. 對(duì)于4個(gè)不同的數(shù)據(jù)集,fSIDWT-SPIRiT與SIDWT-SPIRiT相比,重建圖像的VSNR值、VSSIM值和VHFEN值都很接近,說明兩種方法的重建質(zhì)量相當(dāng).

3.3 不同方法的重建速度比較

最后,從速度上比較各種方法的重建性能.在加速因子R為3～7的情況下,用3種方法對(duì)4個(gè)數(shù)據(jù)集進(jìn)行重建.圖5～6為R=3和R=5時(shí)3種方法對(duì)4個(gè)不同數(shù)據(jù)集進(jìn)行重建的時(shí)間(t)比較.表5給出R為3～7時(shí)在4個(gè)不同數(shù)據(jù)集下fSIDWT-SPIRiT相對(duì)于其他方法重建時(shí)間提升的倍數(shù)(m).

表5 在3～7倍加速時(shí)不同方法對(duì)4個(gè)數(shù)據(jù)集進(jìn)行重建的時(shí)間比較Tab.5 Comparison of reconstruction time for 4 datasets by different methods at 3 to 7 times acceleration

圖5 在3倍加速時(shí)3種方法對(duì)4個(gè)數(shù)據(jù)集的重建速度比較Fig.5 Comparison of the reconstruction speed of the three methods for the 4 datasets at 3 times acceleration

圖6 在5倍加速時(shí)3種方法對(duì)4個(gè)數(shù)據(jù)集的重建速度比較Fig.6 Comparison of the reconstruction speed of the three methods for the 4 datasets at 5 times acceleration

從圖5～6可以看出,對(duì)于4個(gè)不同的數(shù)據(jù)集,fSIDWT-SPIRiT與pFISTA-SPIRiT和SIDWT-SPIRiT兩種方法相比,無論加速因子R=3還是R=5,均能獲得比較接近的VSNR值.但fSIDWT-SPIRiT的收斂速度明顯快于其他兩種方法,表明所提方法可以實(shí)現(xiàn)更快速度的重建,并保證圖像的重建質(zhì)量,進(jìn)一步驗(yàn)證方法的有效性.

從表5可以看出,從3倍加速增加到7倍加速,對(duì)于4個(gè)不同的數(shù)據(jù)集,與pFISTA-SPIRiT相比,fSIDWT-SPIRiT的重建速度分別提升3.2倍、3.7倍、3.8倍和3.4倍,平均提升3.5倍;與SIDWT-SPIRiT相比,fSIDWT-SPIRiT的重建速度分別提升3.5倍、3.4倍、4.9倍和4倍,平均提升3.9倍.由此可以說明,所提方法能夠大大縮短圖像重建時(shí)間,更進(jìn)一步驗(yàn)證該方法的有效性.

3.4 方法的收斂性分析

探討VSNR值、VHFEN值以及相對(duì)誤差(Relative Error, RE)隨迭代次數(shù)k的變化,分析fSIDWT-SPIRiT的收斂性.RE的值計(jì)算如下:

(21)

如圖7所示,在3倍加速下,對(duì)于數(shù)據(jù)集1,在前10次迭代中VSNR值和VHFEN值變化很快,第10次迭代之后變化較為緩慢,當(dāng)VRE<2×10-3即迭代次數(shù)大于16次時(shí),VSNR值和VHFEN值基本不再變化.此時(shí)認(rèn)為所提方法收斂,利用此特性提前終止算法可以減少計(jì)算量.對(duì)于其他3個(gè)數(shù)據(jù)集可以得出類似的結(jié)論.對(duì)于數(shù)據(jù)集2～4,當(dāng)VRE分別小于5×10-3、3×10-3和4×10-3即迭代次數(shù)分別大于13次、23次和23次時(shí),認(rèn)為所提方法收斂,可以提前終止算法以減少計(jì)算量.

圖7 在3倍加速時(shí)對(duì)4個(gè)數(shù)據(jù)集進(jìn)行重建的收斂性分析Fig.7 Convergence analysis of the reconstruction for 4 datasets at 3 times acceleration

4 結(jié)語

基于SIDWT和SPIRiT模型,利用pFISTA技術(shù)進(jìn)行求解,提出一種新的快速并行磁共振成像重建方法——fSIDWT-SPIRiT.在4個(gè)不同活體數(shù)據(jù)集上的仿真實(shí)驗(yàn)表明:與pFISTA-SPIRiT和SIDWT-SPIRiT兩種方法相比,所提方法能獲得與之相當(dāng)?shù)闹亟ㄙ|(zhì)量,而且收斂速度明顯更快,平均提升3.5倍和3.9倍.因此,所提方法既保證圖像的重建質(zhì)量,還顯著減少圖像的重建時(shí)間.本文只選用緊框架SIDWT進(jìn)行實(shí)驗(yàn),后續(xù)研究中將考慮使用其他緊框架進(jìn)一步提升圖像的重建性能.