基于圖像序列差分正則項(xiàng)的實(shí)時心臟磁共振重構(gòu)

衡陽，徐劍峰，陳峰，湯敏*

1.南通大學(xué)電子信息學(xué)院，江蘇南通 226007；2.通科微電子學(xué)院，江蘇南通 226007；3.南通大學(xué)-南通智能信息技術(shù)聯(lián)合研究中心，江蘇南通 226007；4.南通大學(xué)附屬醫(yī)院醫(yī)學(xué)影像科，江蘇南通 226007；5.南通大學(xué)電氣工程學(xué)院，江蘇南通 226007；

心臟磁共振（cardiac magnetic resonance，CMR）具有多平面成像、多序列成像、無電離輻射、圖像質(zhì)量好、可重復(fù)性強(qiáng)等優(yōu)點(diǎn)，能夠全面提供心臟及其鄰近結(jié)構(gòu)的信息，可對心臟形態(tài)結(jié)構(gòu)、心室功能、心肌活性等生理方面做出準(zhǔn)確評價(jià)[1-4]。然而，CMR成像時間有時長達(dá)1 h以上，制約了其臨床應(yīng)用。隨著場強(qiáng)、梯度硬件、脈沖序列等性能上的改進(jìn)，成像速度有所提高，但是高速切換的梯度場制造成本昂貴，還會刺激受檢者的神經(jīng)與肌肉，引起渦流擾動。因此，從硬件設(shè)備角度提高成像速度已基本達(dá)到極限，從算法角度入手是目前解決這一瓶頸問題的突破口[5-6]。

壓縮感知（compressed sensing，CS）理論是指在某一變換域內(nèi)具有稀疏表達(dá)的信號，以遠(yuǎn)低于奈奎斯特采樣定理標(biāo)準(zhǔn)的方式采集數(shù)據(jù)，通過與變換基非相干的隨機(jī)投影采樣，運(yùn)用合適的優(yōu)化算法高概率精確重構(gòu)原始信號[7]。由于 CMR成像屬于動態(tài) MRI范疇，時空相關(guān)度較高，即序列圖像上保持不變的背景部分較多，變化的動態(tài)部分較少，相鄰像素的相似程度較高。因此，可將目前常見的CMR圖像重構(gòu)算法大致分為串行和并行兩大類（圖1）。串行重構(gòu)是指在重構(gòu)算法運(yùn)行之前必須搜集所有幀的K空間數(shù)據(jù)，利用整個數(shù)據(jù)的冗余性重構(gòu)圖像序列，這類算法通常重構(gòu)精度較高，但是需要調(diào)整的參數(shù)較多，因此速度較慢，難以達(dá)到實(shí)時性的要求；并行重構(gòu)則是假設(shè)相鄰兩幀的圖像之間差異很小，通常在密度域或小波域進(jìn)行，通過稀疏正則化重構(gòu)并更新圖像，這類算法的特點(diǎn)是計(jì)算速度較快，但是由于存在誤差累積，易導(dǎo)致重構(gòu)結(jié)果不夠精確，因此不適用于相對較長的序列圖像的重構(gòu)。

圖1 CMR圖像重構(gòu)算法分類及流程。紅色實(shí)線表示并行重構(gòu)，藍(lán)色虛線表示串行重構(gòu)

鑒于上述兩類算法的差異，結(jié)合串行和并行算法的各自優(yōu)點(diǎn)，本研究提出基于圖像序列差分正則項(xiàng)的實(shí)時CMR圖像重構(gòu)算法，力爭在重構(gòu)結(jié)果精確性和實(shí)時性之間取得較好的平衡。該算法的基本思路是：以第1幀圖像作為參考圖像，從K空間數(shù)據(jù)重構(gòu)出第1幀圖像后，利用圖像序列差分正則化使得后續(xù)幀的重構(gòu)均以第1幀圖像為參考，從而同時利用時間域和空間域的稀疏性。采用改進(jìn)的NESTA算法對模型進(jìn)行求解，提高算法的實(shí)時性。以峰值信噪比（peak signal to noise ratio，PSNR）、相對誤差（relativel2norm error，RLNE）和均方根誤差（root mean square error，RMSE）作為定量評價(jià)指標(biāo)，結(jié)合人眼視覺感受以及局部區(qū)域放大，對比分析本算法與 kt FOCUSS[8]、基于運(yùn)動補(bǔ)償?shù)?kt FOCUSS[9]、kt SLR[10]、DTV[11]等主流算法的性能優(yōu)劣。

1 算法關(guān)鍵技術(shù)

1.1 CS-MRI一般模型 CS重構(gòu)算法的一般形式見公式（1），亦可寫作拉格朗日形式，見公式（2）。

一般而言，CS重構(gòu)算法主要分為貪婪算法、約束優(yōu)化算法、迫近算法和同倫算法4大類。其中，貪婪算法采用l0范數(shù)最小化實(shí)現(xiàn)，因其非多項(xiàng)式難解問題，一般只能尋找局部最優(yōu)解，算法簡單快捷，但對于具體應(yīng)用條件有嚴(yán)格規(guī)定。約束優(yōu)化算法采用l1范數(shù)最小化來替代原始的無約束問題，盡管重構(gòu)效果良好，但是計(jì)算量劇增。迫近算法將原始問題轉(zhuǎn)換為相應(yīng)的迫近算子，計(jì)算復(fù)雜度低，易于實(shí)現(xiàn)，被視為解決非光滑、有約束的大規(guī)模數(shù)據(jù)優(yōu)化問題的有力工具，但重構(gòu)精度有待提高。同倫算法通過調(diào)整同倫參數(shù)，用于l1范數(shù)最小化不斷迭代尋找最優(yōu)解。對于算法特點(diǎn)和數(shù)據(jù)量而言，壓縮感知重構(gòu)算法可以分成一階和二階算法兩大類，其中一階算法最多只求一階導(dǎo)數(shù)，計(jì)算速度較快，但是不夠精確；二階算法在求解時較為準(zhǔn)確，但是計(jì)算復(fù)雜度很高，難以適應(yīng)和滿足大規(guī)模計(jì)算的應(yīng)用場合。

1.2 圖像序列差分正則項(xiàng) TV模型利用相鄰像素值差別小的特征對圖像進(jìn)行稀疏表示，在去除噪聲和不需要的小尺度細(xì)節(jié)信息外，保留圖像不連續(xù)的邊緣信息，提高二維圖像重構(gòu)精度，一般模型見公式（3）[12-13]。

其中，下標(biāo)i和j表示像素x的空間位置。

盡管傳統(tǒng)TV在動態(tài)MRI重構(gòu)算法中已取得一定的積極進(jìn)展[14-16]，但是傳統(tǒng)TV是基于圖像的梯度模稀疏這一假設(shè)前提的，對于細(xì)節(jié)和紋理信息較少的圖像具有較好的重構(gòu)結(jié)果；但對于細(xì)節(jié)和紋理信息豐富的CMR圖像而言，傳統(tǒng)TV時常導(dǎo)致過度平滑、邊緣模糊、細(xì)節(jié)失真、階梯效應(yīng)等缺點(diǎn)。本研究提出圖像序列差分（total variation of image sequence，TVIS）正則項(xiàng)的思路，以圖像序列的第1幀作為參考圖像，其余后續(xù)幀的重構(gòu)均以此作為比較，其定義見公式（4）。

其中，r表示參考圖像，即整個圖像序列的第 1幀圖像；?x和?y分別表示沿x和y方向的梯度。比較公式（3）和公式（4），與傳統(tǒng)TV方法不同的是，本研究提出的TVIS可以同時利用時間域和空間域的稀疏性，在CMR圖像的實(shí)時重構(gòu)中應(yīng)該能發(fā)揮更好的作用。因此，基于TVIS模型的壓縮感知重構(gòu)算法的一般形式可寫作公式（5）。

1.3 改進(jìn)的 NESTA算法 NESTA算法是在 Yuri Nesterov思想的基礎(chǔ)上拓展而來的，是屬于迫近算法的一階重構(gòu)算法，可應(yīng)用于求解有約束或無約束的優(yōu)化問題，理論收斂速度快，重構(gòu)結(jié)果精確，魯棒性好，而且算法所需調(diào)節(jié)參數(shù)少[17-18]。NESTA算法用于求解如下問題：

以下重點(diǎn)研究yk和zk的計(jì)算。因?yàn)橹灰?jì)算出它們的值，xk便是兩者的加權(quán)平均。首先，yk的求解因其目標(biāo)函數(shù)是凸函數(shù)，可行域是凸集，因此可采用標(biāo)準(zhǔn)的拉格朗日乘法實(shí)現(xiàn)。構(gòu)造拉格朗日函數(shù)是拉格朗日問題的對偶解，則滿足 Karush-Kuhn-Tucker（KKT）條件：

根據(jù)公式（8）的第4個條件和Sherman-Mprrison-Woodbury求逆公式，可得yk解如公式（9）所示。同理，zk解如公式（10）所示。將公式（9）、（10）代入改進(jìn)的NESTA算法的計(jì)算步驟，便可完成NESTA算法的迭代過程。

綜上所述，改進(jìn)后的NESTA算法迭代步驟主要分成3步：求出的基礎(chǔ)上，加權(quán)平均得到其核心思路是采用平滑技術(shù)，把不光滑的函數(shù)平滑化，然后再用高速便捷的迭代算法來求解。改進(jìn)后的NESTA算法收斂速度可達(dá)級，基本逼近一階算法的理論極限。

1.4 性能評價(jià)指標(biāo) 采用峰值信噪比（peak signal to noise ratio，PSNR）、相對誤差（RLNE）、均方根誤差（root mean square error，RMSE）3個指標(biāo)定量分析和評價(jià)算法的性能特點(diǎn)，其中PSNR是一種評價(jià)圖像品質(zhì)的客觀標(biāo)準(zhǔn)，其值越大表示圖像失真越??；RLNE用于比較完整采樣圖像和重構(gòu)圖像的相對l2范數(shù)誤差，其值越小越好；RMSE說明完整采樣圖像和重構(gòu)圖像之間樣本的離散程度，能夠較好地反映重構(gòu)圖像的精密程度，其值越小越好。上述3個指標(biāo)的計(jì)算分別見公式（11）～（13）。

2 結(jié)果

本實(shí)驗(yàn)在HP工作站Z600上采用MATLAB 2016編程實(shí)現(xiàn)，硬件配置為：Intel Xeon CPU E5606@2.13GHz，24 GB內(nèi)存，64位Windows 7操作系統(tǒng)。使用的完整采樣數(shù)據(jù)來自 kt FOCUSS程序中的full_rad數(shù)據(jù)，大小為384 384 32× ×。該數(shù)據(jù)集采用徑向采樣模式，沿通過K空間中心的徑向直線采集數(shù)據(jù)，具有以下優(yōu)點(diǎn)：①徑向采樣數(shù)據(jù)的每條線含有等量的低頻到高頻信息，有利于圖像的欠采樣重構(gòu)；②徑向采樣模式對決定圖像主要信息的中心數(shù)據(jù)過采樣，因此徑向釆樣對物體運(yùn)動沒有笛卡爾采樣那樣敏感，更適用于CMR圖像重構(gòu)。

2.1 對比算法及參數(shù)設(shè)置 對比算法包括：直接填0的DFFT算法、kt FOCUSS算法[8]、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS 算法[9]、kt SLR 算法[10]、DTV 算法[11]。各個算法的相關(guān)參數(shù)設(shè)置見表1。

表1 本實(shí)驗(yàn)算法與5種對比算法的相關(guān)參數(shù)設(shè)置

2.2 算法精確性分析 本研究算法中，由于第1幀圖像將作為參考圖像，后續(xù)幀圖像的壓縮感知重構(gòu)均以此作為參照，因此第1幀圖像的欠采樣率不能過低，可選定為40%左右，后續(xù)各幀的欠采樣率設(shè)定為 12%。完整采樣圖像序列的第 16幀圖像見圖2A，其中紅色正方形方框放大見圖2B。本研究算法和上述對比算法同時選取第 16幀及其局部區(qū)域放大分別見圖3、4，其中A～F依次為直接填0的DFFT算法、kt FOCUSS算法、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法、kt SLR算法、DTV算法和本研究算法結(jié)果圖。

將本研究算法和上述對比算法所得結(jié)果，沿圖2A中黃色虛線展開可得圖5。

圖2 完整采樣圖像第 16幀及其局部區(qū)域放大顯示。A為完整采樣圖像第 16幀，其中紅色虛線框放大為B圖

圖3 本研究算法和其他對比算法第16幀重構(gòu)結(jié)果。A～F依次為直接填0的DFFT算法、kt FOCUSS算法、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法、kt SLR算法、DTV算法、本研究算法

圖4 本研究算法和對比算法第16幀重構(gòu)結(jié)果局部區(qū)域放大圖。A～F依次為直接填0的DFFT算法、kt FOCUSS算法、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法、kt SLR算法、DTV算法、本研究算法

圖5 本研究算法和對比算法沿圖2黃色虛線展開圖。A～G依次為完整采樣圖像、直接填0的DFFT算法、kt FOCUSS算法、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法、kt SLR算法、DTV算法、本研究算法。與完整采樣圖像相比，直接填0的DFFT方法展開圖中存在不少偽影；箭所指位置，TV方法的展開圖中缺失了上方的灰色條帶，kt FOCUSS和基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法中上方的灰色條帶不夠連續(xù)，kt SLR算法中下方箭所指位置的灰色條帶明顯偏暗

由圖2～4可見，本研究算法得到的重構(gòu)結(jié)果人眼視覺效果最佳，尤其在圖4中標(biāo)出的圓形和矩形框中，黑色像素點(diǎn)的位置及大小與圖2B中所示的完整采樣圖像放大結(jié)果最接近。相比之下，直接填 0的DFFT方法、kt SLR算法和TV算法放大區(qū)域明顯缺少黑色像素點(diǎn)聚集，而kt FOCUSS和基于運(yùn)動補(bǔ)償?shù)?kt FOCUSS算法圓形框中的明亮像素明顯缺失。本研究算法所得重構(gòu)圖像在人眼視覺系統(tǒng)的定性評價(jià)上性能最優(yōu)，最接近完整采樣圖像。

從定量指標(biāo)上進(jìn)行分析，上述算法32幀圖像的性能指標(biāo)對比見圖6，性能指標(biāo)的平均值見表2。由于直接填0的DFFT方法效果最差，PSNR最低，因此在RLNE和RMSE的對比圖上未參與顯示，重點(diǎn)比較其他幾種算法的RLNE和RMSE指標(biāo)。本研究算法幾乎在各幀圖像重構(gòu)中均能得到最高的PSNR值和最小的RMSE值，32幀的性能指標(biāo)平均值也是如此；僅RLNE指標(biāo)，是基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法性能最優(yōu)，可能因?yàn)樵谶\(yùn)動補(bǔ)償環(huán)節(jié)對重構(gòu)圖像進(jìn)行了必要的修正所致。因此，從性能指標(biāo)定量分析來看，本研究算法總體性能仍為最佳（圖6、表2）。

圖6 本研究算法和對比算法的性能指標(biāo)對比。A為PSNR指標(biāo)，B為RLNE指標(biāo)，C為RMSE指標(biāo)

表2 本研究算法和對比算法的性能指標(biāo)平均值

2.3 首幀圖像欠采樣率對算法性能的影響 如前所述，第 1幀圖像作為整個圖像序列的參考圖像，其重構(gòu)精度將顯著影響后續(xù)幀的重構(gòu)精度，因此有必要考查首幀圖像欠采樣率對整個算法性能的影響。首先，固定第1幀圖像的欠采樣率為40%，后續(xù)幀的欠采樣率分別設(shè)置為10%、15%、20%、25%、40%，算法性能指標(biāo)見圖7。然后，將第1幀圖像的欠采樣率分別設(shè)置為100%、75%、50%、25%、15%，而將后續(xù)幀的欠采樣率固定為 25%，算法的性能指標(biāo)見圖8。

圖7 首幀圖像欠采樣率對算法性能的影響。固定第1幀圖像欠采樣率為40%，后續(xù)幀欠采樣率可調(diào)。A為PSNR指標(biāo)，B為RLNE指標(biāo)，C為RMSE指標(biāo)

圖8 首幀圖像欠采樣率對算法性能的影響。第1幀圖像欠采樣率可調(diào)，后續(xù)幀欠采樣率固定為25%。A為PSNR指標(biāo)，B為RLNE指標(biāo)，C為RMSE指標(biāo)

由圖7可見，首幀圖像欠采樣率固定為40%，即使后續(xù)幀欠采樣率僅10%時，本研究算法依然能取得較好的性能指標(biāo)，其PSNR均值為44.6268 dB，RLNE均值0.0374，RMSE均值0.0059，完全能夠滿足臨床輔助診療的需要。隨著后續(xù)幀欠采樣率上升，各項(xiàng)性能指標(biāo)逐步趨優(yōu)，當(dāng)后續(xù)幀欠采樣率也達(dá)40%時，性能指標(biāo)最高可達(dá)PSNR均值為49.3388 dB，RLNE均值0.0185，RMSE均值0.0034。由圖8可見，當(dāng)固定后續(xù)幀欠采樣率為25%時，首幀圖像的欠采樣率對圖像重構(gòu)質(zhì)量的影響不是特別明顯。若首幀圖像完整采樣，則此時的性能指標(biāo)最高可達(dá) PSNR均值為 49.1890 dB，RLNE均值0.0149，RMSE均值0.0035；若首幀圖像欠采樣率僅為15%，此時的性能指標(biāo)也能達(dá)到PSNR均值為45.9503 dB，RLNE均值0.0303，RMSE均值0.0051。綜上所述，為使本研究算法的精確度能夠達(dá)到臨床輔助診療的作用，可以將首幀圖像欠采樣率設(shè)置為25%～40%，后續(xù)幀的欠采樣率設(shè)置為15%～20%，這樣即能確保RLNE和RMSE分別小于0.05和0.005。2.4 算法實(shí)時性分析 本研究算法和對比算法的運(yùn)行時間結(jié)果見表2。相比而言，盡管直接填0的DFFT算法運(yùn)行時間短到可以忽略不計(jì)，但由于其性能太差，無法應(yīng)用于醫(yī)學(xué)圖像臨床輔助診療，因此排除在外。其余算法中，kt SLR的運(yùn)行時間達(dá)2977 s，約需50 min；kt FOCUSS算法運(yùn)行222 s，基本滿足實(shí)時性的要求；基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法因運(yùn)動補(bǔ)償環(huán)節(jié)，其運(yùn)行時間達(dá)882 s，即15 min左右；本研究算法運(yùn)行時間僅36 s，平均每幀圖像重構(gòu)僅需1.125 s，完全滿足醫(yī)學(xué)圖像重構(gòu)實(shí)時性的要求。

3 討論

本研究將壓縮感知理論引入 CMR，提出基于圖像序列差分正則項(xiàng)的實(shí)時CMR圖像重構(gòu)算法，在取得精確重構(gòu)結(jié)果的同時，提高了程序運(yùn)行的實(shí)時性。

在精確性方面，首幀圖像欠采樣率40%，后續(xù)各幀欠采樣率12%時，相比其他5種對比算法，本研究算法獲得最佳的性能評價(jià)指標(biāo)，同時局部區(qū)域放大后黑色像素點(diǎn)的位置和大小最接近完整采樣圖像，即人眼視覺系統(tǒng)的主觀感受最佳。

首幀圖像欠采樣率對算法性能影響不大。一般來說，首幀圖像欠采樣率為 25%～40%，后續(xù)幀欠采樣率為 15%～20%，即能取得令人滿意的重構(gòu)結(jié)果，完全能夠滿足臨床輔助診斷的需求。

實(shí)時性方面，對于實(shí)驗(yàn)使用的32幀CMR圖像而言，本研究算法重構(gòu)僅需36 s，比DTV算法運(yùn)行時間略短，與kt FOCUSS算法、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS算法以及kt SLR算法相比，本研究算法的運(yùn)行時間僅為其1/10甚至1/100。因此，本研究算法可以在獲得精確重構(gòu)結(jié)果的同時滿足實(shí)時性的要求，在臨床上具有廣闊的應(yīng)用前景。

當(dāng)然，本研究算法尚有進(jìn)一步提升的空間。近年來，基于低秩稀疏分解的CS-MRI圖像重構(gòu)研究方興未艾。有學(xué)者首次將低秩稀疏分解模型應(yīng)用于MRI圖像重構(gòu)，初步應(yīng)用于心臟灌注、心臟電影、血管成像、胸腹部 MRI等[19]。有學(xué)者利用魯棒主成分分析提取動態(tài) MRI中的動態(tài)組織部分，提出基于稀疏和低秩先驗(yàn)分離的重構(gòu)算法[20]。有人提出基于部分可分離模型和壓縮感知相結(jié)合的Stepped-Sparse PS算法[21]。還有研究者初步提出多尺度和多維低秩矩陣分解，并應(yīng)用于人臉圖像、監(jiān)控視頻、動態(tài) MRI等的多尺度建模和協(xié)同濾波[22-23]。

今后的研究方向包括：①將本研究算法與基于低秩稀疏分解算法進(jìn)行性能對比和分析；②將本研究算法應(yīng)用于更多CMR序列的快速重構(gòu)，檢驗(yàn)算法的普適性；③將本研究算法應(yīng)用于兒童心肌病、復(fù)雜先天性心臟病、心臟二尖瓣病變圖像重構(gòu)中，并與臨床醫(yī)師的主觀判斷相對照，從而檢驗(yàn)算法在臨床疾病輔助診療中的實(shí)際作用。

綜上所述，壓縮感知理論引入 CMR，可以加快成像速度，提高成像質(zhì)量，提高檢查舒適度和檢查質(zhì)量，促進(jìn)循環(huán)系統(tǒng)疾病診療新技術(shù)的發(fā)展。本研究結(jié)合串行和并行算法的各自優(yōu)點(diǎn)，提出基于圖像序列差分正則項(xiàng)的實(shí)時CMR圖像重構(gòu)算法，基本思路是以第1幀圖像作為參考圖像，從K空間數(shù)據(jù)重構(gòu)出第1幀圖像，此后利用圖像序列差分正則化使得后續(xù)幀的重構(gòu)均以第1幀圖像為參考，從而同時利用時間域和空間域的稀疏性。采用改進(jìn)的NESTA算法對模型進(jìn)行求解，提高算法的實(shí)時性。以PSNR、RLNE和RMSE作為定量評價(jià)指標(biāo)，結(jié)合人眼視覺感受、局部區(qū)域放大和沿時間軸展開等定性分析方法，充分說明本研究算法較kt FOCUSS、基于運(yùn)動補(bǔ)償?shù)膋t FOCUSS、kt SLR、DTV等主流算法的性能優(yōu)勢，能在獲取較為精確的重構(gòu)結(jié)果的同時，很好地滿足醫(yī)學(xué)圖像處理實(shí)時性的要求。