自动化学报  2018, Vol. 44 Issue (3): 490-505   PDF    
联合相邻帧预测的心脏磁共振电影成像方法
董家林1, 洪明坚1, 张海标1, 葛永新1     
1. 重庆大学软件学院 重庆 401331
摘要: 快速动态磁共振成像可以通过减少采样量来缩短信号的采集时间.因此,从下采样的数据中重建出高质量的图像成为研究的热点.目前,常见的重建方法利用动态图像序列的稀疏表示实现高质量的重建.本文提出了一种联合相邻帧预测(Joint adjacent-frame prediction,JAFP)的重建方法,首先根据动态图像序列相邻帧之间高度的相似性,联合预测当前帧图像,获得稀疏的图像差;其次,利用图像差序列在时间域的拟周期特性,通过傅里叶变换进一步提高图像差序列的稀疏度.在此基础上构建动态成像模型,并在压缩感知(Compressed sensing,CS)框架下进行求解.该方法可将前一次的重建结果作为新的输入,从而形成迭代算法.采用两个磁共振心脏电影成像数据集对提出的方法进行了实验验证,并与k-t FOCUSS ME/MC和MASTeR进行了比较.实验结果表明,该方法联合相邻帧改进了预测图像的效果,提升了重建图像的质量,具有广泛的应用价值.
关键词: 动态磁共振成像     联合相邻帧     傅里叶变换     压缩感知    
Joint Adjacent-frame Prediction for Cardiac Cine MR Imaging
DONG Jia-Lin1, HONG Ming-Jian1, ZHANG Hai-Biao1, GE Yong-Xin1     
1. School of Software Engineering, Chongqing University, Chongqing 401331
Manuscript received : May 24, 2016, accepted: January 4, 2017.
Foundation Item: Supported by Fundamental Research Funds for the Central Universities (106112017CDJQJ158834)
Author brief: DONG Jia-Lin Master student at the School of Software Engineering, Chongqing University. He received his bachelor degree from Chongqing University in 2015. His research interest covers MRI image reconstruction and pattern recognition;
ZHANG Hai-Biao Master student at the School of Software Engineering, Chongqing University. He received his bachelor degree from Chongqing University in 2014. His main research interest is dynimaic MRI;
GE Yong-Xin Associate professor at the School of Software Engineering, Chongqing University. He received his bachelor degree in information and computing science, master degree in operational research and cybernetics, and Ph. D. degree in computer science and technology from Chongqing University in 2003, 2006 and 2011, respectively. His research interest covers computer vision, image processing, and pattern recognition
Corresponding author. HONG Ming-Jian Associate professor at the School of Software Engineering, Chongqing University. He received his bachelor degree and master degree in applied mathematics, Ph. D. degree in instrument science and technology from Chongqing University in 1999, 2002 and 2011, respectively. His research interest covers spectroscopy and MRI reconstruction. Corresponding author of this paper
Recommended by Associate Editor ZHANG Dao-Qiang
Abstract: Dynamic magnetic resonance imaging can be accelerated to reduce the signal acquisition time by under-sampling k-space data, so the quality of reconstructed images from under-sampled data has become the focus of research. Currently, the common approaches use sparse representation of dynamic image sequences to improve the reconstruction. This paper proposes a new method named joint adjacent-frame prediction (JAFP) based on the similarity between adjacent frames of dynamic image sequences. The JAFP promotes the quality of the predicted image sequence by jointsing adjacent frames prediction. Meanwhile, observing the quasi-periodicity of the difference of sequences, it further improves the sparsity by applying Fourier transform to the difference sequence along the time direction. Then a dynamic imaging model is setup to incorporate the fidelity constraint and joint sparse promotion, and it is solved in the framework of compressive sensing (CS). Two cardiac cine MR datasets are evaluated to verify the proposed method, and the proposed method is compared with k-t FOCUSS with ME/MC and MASTeR to show the better quality of reconstruction. Experimental results show that JAFP can improve the quality of reconstructed image and has important application value.
Key words: Dynamic magnetic resonance imaging     joint adjacent-frame     Fourier transform     compressed sensing (CS)    

磁共振成像(Magnetic resonance imaging, MRI)作为一种医学成像手段, 以其良好的软组织成像、细致的图像颗粒度、较高的空间分辨率以及无创的特性广泛应用于医学领域.然而, 由于MRI的扫描时间比较长, 患者在成像期间会产生自主或非自主的运动, 导致图像失真, 从而影响了医生诊断, 限制了MRI在临床医学中的应用.动态磁共振成像(Dynamic MRI)主要用于捕捉组织器官的动态变化或运动过程, 对成像速度的要求更高, 因此加快成像速度尤其重要.

受MRI扫描仪采集速度的限制, 一般通过减少采样量以缩短扫描时间, 再从下采样的数据中重建出高分辨率的图像.压缩感知(Compressed sensing, CS) [1-4]理论表明, 只要图像在某个变换域是稀疏的, 通过随机下采样就能从少量的下采样数据中重建出高质量的图像.由于CS理论突破了传统Nyquist采样定理的瓶颈, 目前已在雷达成像[5]、红外成像[6]、磁共振成像[7]等领域得到广泛的应用. CS理论应用条件主要包括信号具有稀疏表示、不相干的(Incoherent)采样以及有效的非线性重建算法三个方面[8-9].常用的稀疏变换有傅里叶变换(Fourier transform, FT)和小波变换等[10-11].近年来, 基于分块的自适应稀疏表示也得到了广泛应用[12-15].

动态磁共振成像序列不仅在空间上具有一定的稀疏特性, 而且在时间上也存在较大冗余[16].因此, 如何有效利用动态磁共振成像在时间域上的冗余特性, 是CS理论运用于动态磁共振成像的关键.常见的基于CS理论的动态磁共振成像方法有k-t SPARSE [17]、k-t FOCUSS [18]、k-t FOCUSS ME/MC [19]和MASTeR [20]等. k-t SPARSE利用动态磁共振成像在空间域和时间域冗余的特性, 采用小波变换作为图像序列在空间域的稀疏约束, 同时加入沿时间轴的傅里叶变换作为图像序列在时间域的稀疏约束, 建立CS成像模型重建图像序列.然而k-t SPARSE求解过程采用非线性共轭梯度算法[21], 求解时计算量大, 且收敛缓慢、求解时间较长. k-t FOCUSS利用了动态磁共振成像在时间域的冗余, 采用沿时间轴的傅里叶变换作为稀疏约束构建成像模型, 并利用沿时间轴的平均图像作为每一帧图像的预测图像, 最后通过FOCUSS [22]方法重建图像差.然而, k-t FOCUSS采用沿时间轴的平均图像对图像序列进行预测, 只能较为准确地估计少部分图像, 无法保证大多数时间节点上的预测效果[23].为了克服k-t FOCUSS的这一不足, k-t FOCUSS ME/MC建议使用一帧或两帧全采样的图像作为参考图像, 再利用块匹配算法[24]计算参考图像与目标图像之间的运动估计, 提升了预测图像的准确性.然而, 在没有全采样参考图像时, 该方法仍然采用沿时间轴的平均图像作为参考图像, 而这样的选择在影响运动估计效果的同时, 也降低了重建图像的质量. MASTeR方法通过计算相邻前后两帧图像之间的运动量, 估计出预测图像, 并利用预测图像与目标图像差的稀疏性, 建立具有两个稀疏约束的成像模型.具体来说, 该方法计算出前一帧与当前帧的运动估计, 得到前向预测图像, 同时计算出后一帧与当前帧的运动估计, 得到后向预测图像, 并分别以前、后向预测图像与目标图像的差作为稀疏约束建立成像模型.

本文提出了一种新的联合相邻帧预测(Joint adjacent-frame prediction, JAFP), 根据动态磁共振成像序列相邻前后帧之间较高的相似性, 分别求解相邻前后帧之间的运动估计, 联合预测中间帧图像.同时, 利用目标图像序列与预测图像序列的图像差在时间上的拟周期性, 增加图像差在时间轴的稀疏约束, 建立基于压缩感知的成像模型. JAFP方法从以下两个方面提升了动态磁共振图像序列的稀疏表示: 1)联合相邻前后两帧运动估计所得图像的预测效果优于前一帧向后运动估计的图像或后一帧向前运动估计的图像, 因此计算所得的图像差更为稀疏; 2)在时间域增加稀疏变换, 进一步稀疏了图像差信号.实验结果表明, 该方法有效提升了动态图像序列的稀疏度, 显著改进了重建图像的质量.

1 成像模型 1.1 稀疏动态成像模型

在动态磁共振成像中, 利用动态磁共振图像在时空域的冗余性, 建立基于压缩感知理论的成像模型.

$ \begin{align} \widehat{{M}}= \arg\min\left \{ \left \| {y} - \phi {M} \right \|_{2} +\alpha \left \| \psi {M} \right \|_{1} \right \} \end{align} $ (1)

其中, ${M}\in {\bf C}^{n\times p\times t}$ $t$ 帧、 $n\times p$ 大小的图像, ${y}$ 为下采样的 $k$ 空间数据, $\phi$ 为下采样傅里叶变换, $\psi$ 为稀疏变换.根据压缩感知理论, 在同等采样量的情况下, 越稀疏的表示, 图像重建质量越高[7].

为了获得更稀疏的表示, k-t SPARSE选择空间域的小波变换和时间域的傅里叶变换作为稀疏变换.而k-t FOCUSS, k-t FOCUSS ME/MC, MASTeR等方法则将目标图像与预测图像的差作为成像序列在空间域的稀疏约束建立成像模型如下:

$ \begin{align} \widehat{{M}}= \arg\min\left \{ \left \| {y} - \phi {M} \right \|_{2} +\alpha \left \| \psi \left ( {M} - \widehat{{M}} \right ) \right \|_{1} \right \} \end{align} $ (2)

其中, $\widehat{{M}}$ 为预测图像.在k-t FOCUSS中, $\widehat{{M}}$ 选择沿时间轴的平均图像, 对心脏成像时, 由于心脏跳动的拟周期性, 取 $\psi$ 为沿时间轴的傅里叶变换, 以获得最稀疏的表示; k-t FOCUSS ME/MC改进了k-t FOCUSS中预测图像 $\widehat{{M}}$ 的计算方法, 通过块匹配算法计算沿时间轴的平均图像与目标图像之间的运动估计, 并求得 $\widehat{{M}}$ ; MASTeR把 $\widehat{{M}}$ 分成前向 $\widehat{{M}}_{{F}}$ 和后向 $\widehat{{M}}_{{B}}$ 两部分, 前向 $\widehat{{M}}_{{F}}$ 通过计算前一帧与当前帧的运动估计获得, 后向 $\widehat{{M}}_{{B}}$ 通过计算后一帧与当前帧的运动估计获得, 最后用两个稀疏约束 $ \| {M}-\widehat{{M}}_{{F}} \|_{1}$ $ \| {M}-\widehat{{M}}_{{B}} \|_{1}$ 构建成像模型.

1.2 联合相邻帧预测方法

相邻帧图像即为当前帧图像的前一帧图像或后一帧图像.由于它们与当前帧的扫描时间间隔较短, 因此相邻帧图像往往差别较小, 如图 1所示.其中, 从左至右为第4 ~ 6、12 ~ 14帧图像的感兴趣区域(Region of interest, ROI), 差别非常小, 肉眼难以辨别, 但从其与第5帧图像ROI的误差, 可以直观地看出相邻帧图像与非相邻帧图像之间的差异.在对当前帧图像进行预测时, 便可以采用当前帧的前一帧图像或者后一帧图像对当前帧图像进行预测[25]. 图 2用结构相似性(Structural similarity, SSIM) [26]比较了各帧之间的相似度.

图 1 相邻帧图像与非相邻帧图像之间的差异 Figure 1 The difference of between adjacent frame and non-adjacent frame
图 2 帧间相似度 Figure 2 The similarity

由于相邻帧之间的差异较小, 因此选择相邻前后两帧运动估计后的图像预测效果更优.记前向运动估计算子为 ${F}$ , 后向运动估计算子为 ${B}$ , 则前向运动估计后的第 $i$ 帧预测图像 $\widehat{m}_{{F}}\left ( i \right )$ 和后向运动估计后的的第 $i$ 帧预测图像 $\widehat{m}_{{B}}\left ( i \right )$ 可表示为

$ \begin{align} \widehat{m}_{{F}}\left ( i \right )={F}m\left ( i-1 \right ) \end{align} $ (3)
$ \widehat{m}_{{B}}\left ( i \right )={B}m\left ( i+1 \right ) $ (4)

其中, $m\left ( i \right )$ 表示成像序列的第 $i$ 帧图像, 由于运动估计后的预测图像和本真图像之间存在差异, 记 $u\left ( i \right )$ 为前向运动估计后的预测图像序列与本真图像序列的差函数, $v\left ( i \right )$ 为后向运动估计后的预测图像序列与本真图像序列的差函数, 则有

$ \begin{align} &m\left ( i \right )=\widehat{m}_{{F}}\left ( i \right )+u\left ( i \right ) \end{align} $ (5)
$ \begin{align} &m\left ( i \right )=\widehat{m}_{{B}}\left ( i \right )+v\left ( i \right ) \end{align} $ (6)

结合式(5)和式(6)可得

$ \begin{align} 2m\left ( i \right )=&\ \widehat{m}_{{F}}\left ( i \right )+\widehat{m}_{{B}}\left ( i \right )+u\left ( i \right )+v\left ( i \right ) \end{align} $ (7)
$ \begin{align} e\left ( i \right )=&\ m\left ( i \right )-\frac{1}{2}\left (\widehat{m}_{{F}}\left ( i \right )+\widehat{m}_{{B}}\left ( i \right ) \right )= \nonumber\\ &\ \frac{1}{2}\left (u\left ( i \right )+v\left ( i \right ) \right ) \end{align} $ (8)

其中, $e\left ( i \right )$ 为联合前后帧运动估计后的预测图与本真图像的差.然而, 第一帧图像和最后一帧图像只有一个相邻帧.由于心脏跳动是拟周期的, 如图 2(b)所示, 因此可以通过周期延拓进行边界处理.因此, 有

$ \begin{align} m\left ( i \right )\approx m\left ( i+T \right ) \end{align} $ (9)

同时, $\widehat{m}_{{F}}\left ( i \right )$ $\widehat{m}_{{B}}\left ( i \right )$ 仍然保留了 $m(i)$ 的拟周期性, 即

$ \begin{align} \widehat{m}_{{F}}\left ( i \right )\approx \widehat{m}_{{F}}\left ( i+T \right ), \widehat{m}_{{B}}\left ( i \right )\approx \widehat{m}_{{B}}\left ( i+T \right ) \end{align} $ (10)

结合式(8)不难看出, 其图像差 $e\left ( i \right )$ 在时间域上仍然是拟周期的.

$ \begin{align} e\left ( i \right )=&\ m\left ( i \right )-\frac{1}{2}\left (\widehat{m}_{{F}}\left ( i \right )+\widehat{m}_{{B}}\left ( i \right ) \right )\approx \nonumber\\ &\ m\left ( i+T \right )-\frac{1}{2}\left(\widehat{m}_{{F}}\left ( i+T \right)+\right.\notag\\ & \left.\widehat{m}_{{B}}\left ( i+T \right) \right)= e\left ( i+T \right ) \end{align} $ (11)

由此可知, 图像差序列在傅里叶变换域上是稀疏的.因此, 增加沿时间轴的傅里叶变换, 更有利于重建质量的提升.

1.3 成像模型

本文提出的成像方法JAFP在联合相邻帧预测的同时, 增加时间域上的稀疏变换, 获得更稀疏的图像差表示, 以提升重建质量.根据式(2)和式(8)可得JAFP的成像模型如下:

$ \begin{align} \widehat{{M}} =&\ \arg\min \Bigg\{ \left \| {y} - \phi {M} \right \|_{2} + \nonumber\\ &\ \alpha \left \| \psi \left ( {M} -\frac{1}{2} \left (\widehat{{M}}_{{F}}+\widehat{{M}}_{{B}} \right ) \right ) \right \|_{1} \Bigg\} \end{align} $ (12)

其中, $\psi$ 为沿时间轴的傅里叶变换, $\widehat{{M}}_{{F}}$ 表示前向运动估计后的预测图像序列, $\widehat{{M}}_{{B}}$ 表示后向运动估计后的预测图像序列.运动估计算子 ${F}$ ${B}$ 是通过计算初始图像序列相邻前后帧图像的运动估计所得.

k-t FOCUSS ME/MC和MASTeR在模型上都利用了图像差作为成像序列在空间域上的稀疏约束, 然而, 在没有全采样参考图像时, k-t FOCUSS ME/MC采用的沿时间轴平均图像作为参考图像, 尽管其利用了图像差序列在时间上的冗余进一步稀疏了图像差序列, 但其质量较差的参考图像影响了运动估计的结果, 限制了重建图像质量的进一步提升.而MASTeR选用相邻帧图像的运动估计预测当前帧图像, 在预测准确性上优于k-t FOCUSS ME/MC, 然而MASTeR忽略了图像差序列在时间上的冗余.而本文提出的成像模型在综合二者优势的同时, 采用联合相邻帧预测的方式进一步提升了预测效果.

1) JAFP与MASTeR类似地都利用了相邻帧图像的运动估计, 但JAFP联合相邻前后两帧图像计算所得的运动估计, 替代MASTeR所用的单帧运动估计.从模型上看, MASTeR的稀疏约束为

$ \begin{align} \left \| {M}-\widehat{{M}}_{{F}}\right \|_{1}+\left \| {M}-\widehat{{M}}_{{B}}\right \|_{1} \end{align} $ (13)

而JAFP的稀疏约束为

$ \begin{align} \left \| \left ({M}-\widehat{{M}}_{{F}} \right )+ \left ({M}-\widehat{{M}}_{{B}} \right )\right \|_{1} \end{align} $ (14)

显然, 根据1范数的性质[27]

$ \begin{align} &\left \| \left ({M}-\widehat{{M}}_{{F}} \right )+ \left ({M}-\widehat{{M}}_{{B}} \right )\right \|_{1} \leq \nonumber\\ &\qquad\left \| {M}-\widehat{{M}}_{{F}}\right \|_{1}+\left \| {M}-\widehat{{M}}_{{B}}\right \|_{1} \end{align} $ (15)

JAFP具有更优的稀疏表示.

2) JAFP利用目标图像与预测图像的图像差在时间轴上的拟周期性(由式(11)可知)增加时间轴的稀疏约束 $\psi$ , 进一步提升了稀疏性, 弥补了MASTeR的不足.

2 重建算法

JAFP成像方法通过初始图像序列计算相邻前后帧图像的运动估计算子求得前向预测图像序列和后向预测图像序列, 进而计算预测图像序列与目标图像序列的图像差, 建立CS成像模型完成重建.重建完成后, 可以用重建所得结果更替初始图像序列, 再次运用JAFP方法完成重建.该过程可迭代多次, 随着迭代次数的增加, 被更替的预测图像序列愈加接近目标图像, 使得求解目标(图像差)愈稀疏, 有利于提升下一次迭代的图像质量, 当图像质量的改进小于阈值时, 迭代停止.

算法流程如下:

步骤1. 初始化:根据CS理论, 从下采样 $k$ 空间中重建获得初始图像序列, 即根据式(1)求得 $\widetilde{{M}}$ 后, 令 $k=1$ , ${M}^{0}= \widetilde{{M}}$ ;

步骤2. 运动估计:使用重建的图像序列 ${M}^{k-1}$ 估计相邻帧之间的运动, 利用复小波变换(Complex wavelet transform, CWT)运动估计方法, 求得运动估计算子 ${F}$ ${B}$ ;

步骤3. 运动补偿:利用图像差稀疏这一特性建立新CS成像模型

$ \begin{align} \widehat{{M}} =&\ \arg\min \Bigg\{ \left \| {y} - \phi {M}^{k} \right \|_{2} + \nonumber\\ &\ \alpha \left \| \psi \left ( {M}^{k} -\frac{1}{2} \left (\widehat{{M}}_{{F}}^{k}+\widehat{{M}}_{{B}}^{k} \right ) \right ) \right \|_{1} \Bigg\} \end{align} $ (16)

并求解得到 $\widehat{{M}}$ ;

步骤4. 如果 $ \| \widehat{{M}}-{M}^{k} \|/ \| {M}^{k} \|\leq \varepsilon $ , 则停止迭代, 输出 $\widehat{{M}}$ .否则, 令 ${M}^{k}=\widehat{{M}}, k=k+1$ , 重复步骤2.

算法流程图如图 3所示.

图 3 算法流程图 Figure 3 Algorithm flow chart
3 实验方法 3.1 数据集

数据集1 [18]:该数据集为k-t FOCUSS实验过程所用数据集, 系利用飞利浦1.5 T磁共振扫描仪采集得到一名病人的心脏电影成像数据, 包括25帧全采样 $k$ 空间数据, 视场(Field of view, FOV)大小为306.0 mm $\times $ 260.00 mm, 扫描矩阵为220像素× $256$ 像素, 采集顺序为平衡稳态自由进动(Balance steady state free precession, bSSFP), 倾角50°, TR = $3.45$ ms, 心脏跳动频率75 bpm.数据集包括心脏舒张和收缩两个阶段.

数据集2:该心脏电影数据共25帧, 利用飞利浦1.5 T磁共振扫描仪采集获得一名自愿者的心脏电影成像数据, FOV为345.00 mm × 270.00 mm, 扫描矩阵为 $256$ 像素× $256$ 像素.数据是使用稳态自由进动(SSFP)顺序的方式获得, 稳态自由进动过程倾角为50°, TR = 3.45 ms, 心脏跳动频率66 bpm.数据集包括心脏舒张和收缩两个阶段.

3.2 采样模式

由于 $k$ 空间中低频区域即中心部分能量较大, 高频区域能量较低, 低频区域数据包含图像的对比度和信噪比, 高频区域数据包含图像的边缘和细节等特征[28].因此采样时尽可能多地采集低频区域数据, 其他区域在采样时尽可能保证随机.在保证图像一定信噪比的同时, 也满足CS理论中不相干(Incoherent)采样的要求.

实验主要采用笛卡尔下采样模式, $k$ 空间中心的8行数据全采样以保证重建图像的信噪比, 余下 $k$ 空间数据通过标准高斯分布随机采集获得.尽管在实验过程中选择了笛卡尔采样, 但该算法并不限定采样模式.

3.3 求解方法

利用NESTA [29]算法求解式(1)获得初始图像, 算子 ${F}$ ${B}$ 可通过运动估计方法求解.常用的运动估计方法包括块匹配算法(Overlapped block motion compensation, OBMC)、复小波变换(CWT) [30]和光流法(Optical flow, OP) [31]等, 其中OBMC和OP一般不适宜于处理包含噪声或低信噪比的动态图像序列[32], 而JAFP算法是一个迭代收敛的过程, 初始图像序列往往有较大的噪声, 因此选择能抑制噪声影响的CWT运动估计方法.

由于心脏图像差在时间域是拟周期的, 因此稀疏变换Ψ选择沿时间轴的傅里叶逆变换.然后, 再次利用NESTA算法求解式(16), 求解出 $\widehat{{M}}$ 后, 更新 ${F}$ ${B}$ , 并重复该过程, 直至收敛.

3.4 实验结果的评估方法

为了更好地对比实验的重建效果, 通过计算原动态磁共振成像数据集与下采样情况下的重建图像数据的差别定量分析重建的效果.结构相似性(SSIM)的值可以更好地反映人眼的主观感受[33], 因此计算每一帧重建图像与本真图像的结构相似性评估重建结果的优劣.

SSIM的计算公式如下:

$ \begin{align} {\rm SSIM}\left ( x, y \right )=\frac{\left ( 2\mu_{x}\mu_{y} + c_{1} \right )\left ( 2\sigma_{xy} + c_{2} \right )}{\left ( \mu_{x}^{2}+\mu_{y}^{2} + c_{1} \right )\left (\sigma_{x}^{2}+ \sigma_{y}^{2} + c_{2} \right )} \end{align} $ (17)

其中, $x$ , $y$ 分别为两幅图像, $\mu_{x}$ 表示 $x$ 的平均值, $\sigma_{x}^{2}$ 表示 $x$ 的方差, $\sigma_{xy}$ 表示 $x$ $y$ 的协方差, $c_{1}$ = $(k_{1}L)$ , $c_{2}=(k_{2}L)$ , $L$ 为像素值的动态范围, $k_{1}$ 取0.01, $k_{2}$ 取0.03.

同时, 高频误差范数(High-frequency error norm, HFEN) [34]能够评估重建图像边缘和细节特征的质量.因此, 计算原动态磁共振成像数据集与下采样重建图像数据集的高频误差范数(HFEN)定量分析重建结果的质量.

HFEN的计算公式如下:

$ \begin{align} {\rm HFEN}\left ( x, \widehat{x} \right )=\frac{1}{N}\sum _{i=1}^{N}\frac{\left \| LoG(x_{i})-LoG(\widehat{x_{i}}) \right \|_{F}^{2}}{\left \| LoG(x_{i}) \right \|_{F}^{2}} \end{align} $ (18)

其中, $LoG(x)$ 为拉普拉斯高斯滤波函数, 模版尺寸15, 滤波器标准差为1.5.

4 实验结果及分析 4.1 数据集1的实验结果及分析

图 4(a)图 5(a)分别展示了数据集1在4倍下采样和8倍下采样时, 三种重建方法重建后每一帧图像感兴趣区域(ROI)的SSIM对比图. 图 4(b)图 5(b)分别展示了ROI的HFEN对比图.由于JAFP方法进一步提升了预测心脏运动信息的准确性, 其重建结果在SSIM和HFEN两项评价指标上均超越MASTeR和k-t FOCUSS ME/MC.同时, 实验选取了两个心脏壁区域(图 6(a)图 7(a)方框标注所示), 进而比较了各重建方法重建该区域的平均灰度变化, 具体如图 6图 7所示.根据图中所示重建结果与真实结果的误差绝对值可知, 基于JAFP方法的重建结果在8倍下采样时收缩过程的灰度变化总体上也优于MASTeR和k-t FOCUSS ME/MC.

图 4 数据集1在4倍下采样时重建质量比较 Figure 4 The comparison of reconstruction quality with acceleration factor of 4 for dataset 1
图 5 数据集1在8倍下采样时重建质量比较 Figure 5 The comparison of reconstruction quality with acceleration factor of 8 for dataset 1
图 6 数据集1在8倍下采样时左心脏壁灰度均值的比较 Figure 6 The comparison of reconstructed gray average values and errors with acceleration factor of 8 for left myocardial wall of dataset 1
图 7 数据集1在8倍下采样时右心脏壁灰度均值的比较 Figure 7 The comparison of reconstructed gray average values and errors with acceleration factor of 8 for right myocardial wall of dataset 1

图 8呈现了数据集1在4倍下采样和8倍下采样时三种重建方法重建的ROI及其误差图, 图中选取了三帧图像, 分别为心脏收缩阶段的第4帧、收缩完全的第9帧以及舒张阶段的第18帧. 图 9比较了各重建方法重建的变化较明显且靠近心脏壁位置(原成像序列的第130列和第160列)的截面及其误差.由于JAFP在综合MASTeR和k-t FOCUSS ME/MC优势的同时, 还联合相邻帧预测进一步提升了预测的准确性, 其重建的ROI误差以及截面误差都明显少于MASTeR和k-t FOCUSS ME/MC.

图 8 数据集1在4倍和8倍下采样时, 重建的ROI及其误差图 Figure 8 The reconstructed ROI and errors with acceleration factors of 4 and 8 for dataset 1
图 9 在4倍和8倍下采样时, 重建的时间域截面及其误差图 Figure 9 The reconstructed temporal slices and errors with acceleration factors of 4 and 8

与此同时, 定量分析在不同加速倍数下k-t FOCUSS ME/MC、MASTeR和JAFP三种方法重建的图像序列的重建效果, 计算了在2, 4, 8, 10, 16倍加速倍数下共25帧图像序列在SSIM和HFEN两项评价指标上的均值和方差, 如图 10所示.从图中不难看出, JAFP重建方法在各加速倍数下, 重建结果的SSIM以及HFEN两项评价指标均优于k-t FOCUSS ME/MC和MASTeR.由于JAFP利用相邻帧之间的运动估计, 联合预测中间帧图像, 提升了预测效果, 同时, 时间域稀疏约束的增加, 使得其重建结果在更高倍的下采样时重建效果提升明显.

图 10 数据集1在不同下采样倍数(2, 4, 8, 10, 16)下, 不同方法重建综合质量对比图 Figure 10 The comparison of reconstructed quality with acceleration factors of 2, 4, 8, 10 and 16 for dataset 1
4.2 数据集2的实验结果及分析

图 11(a)图 12(a)展示了数据集2分别在4倍下采样和8倍下采样时, 三种重建方法重建后每一帧图像ROI区域的SSIM对比图. 图 11(b)图 12(b)展示了ROI区域的HFEN对比图.显然, 具有更准确预测的JAFP方法重建结果在SSIM和HFEN两项评价指标上均为最优. 图 13图 14分别比较了在8倍下采样时各方法重建结果在两个心脏壁区域(方框标注区域)收缩过程的灰度均值变化情况, 由于JAFP方法较准确地预测了心脏的运动过程, 使其在收缩过程重建结果与真实结果灰度均值的误差绝对值明显小于MASTeR和k-t FOCUSS ME/MC.

图 11 数据集2在4倍下采样时重建质量比较 Figure 11 The comparison of reconstruction quality with acceleration factor of 4 for dataset 2
图 12 数据集2在8倍下采样时重建质量比较 Figure 12 The comparison of reconstruction quality with acceleration factor of 8 for dataset 2
图 13 数据集2在8倍下采样时左心脏壁灰度均值的比较 Figure 13 The comparison of reconstructed gray average values and errors with acceleration factor of 8 for left myocardial wall of dataset 2
图 14 数据集2在8倍下采样时右心脏壁灰度均值的比较 Figure 14 The comparison of reconstructed gray average values and errors with acceleration factor of 8 for right myocardial wall of dataset 2

在重建误差上, 图 15分别呈现了数据集2在4倍下采样和8倍下采样时三种重建方法重建的ROI及其误差图.图 16呈现了数据集2在4倍下采样和8倍下采样时三种重建方法重建的时间域截面及其误差图. 图 15共选取三帧图像, 分别为心脏收缩阶段的第4帧、心脏收缩完全的第9帧以及舒张阶段的第16帧. 图 16选取了重建截面为原成像序列的第130列和第160列.尽管与数据集1完全不同, 但从图像上看, 选取的这两个位置仍然较为靠近心脏外壁, 变化幅度较大, 便于对比分析.由于预测效果的提升以及时间域稀疏约束的增加, 无论是在重建的ROI误差上, 还是在重建截面上, JAFP的重建误差明显少于MASTeR和k-t FOCUSS ME/MC.

图 15 数据集2在4倍和8倍下采样时, 重建的ROI及其误差图 Figure 15 The reconstructed ROI and errors with acceleration factors of 4 and 8 for dataset 2
图 16 数据集2在4倍和8倍下采样时, 重建的时间域截面及其误差图 Figure 16 The reconstructed temporal slices and errors with acceleration factors of 4 and 8 for dataset 2
5 讨论

本文提出的JAFP方法与k-t FOCUSSME/MC、MASTeR在模型上都利用了图像差作为成像序列在空间域上的稀疏约束, 但k-t FOCUSS ME/MC运动预测效果较差, 而MASTeR又忽略了图像差序列在时间上的冗余, 限制了重建质量的进一步提升. JAFP在弥补二者不足的同时, 还联合相邻帧预测进一步提升了预测的准确性.此外, 由于JAFP方法在求解时仅需计算带有一个约束的成像模型, 使得其在计算效率上优于MASTeR.

5.1 联合相邻帧预测与其他预测方式的比较

JAFP在成像模型上与k-t FOCUSS ME/MC、MASTeR都用到了预测图像, 但JAFP的预测效果明显优于二者. JAFP利用相邻前后帧计算运动估计, 在不需要引入额外的参考图像的情况下, 弥补了k-t FOCUSS ME/MC预测图像质量较差的缺点.同时, JAFP联合了相邻前后帧的运动估计, 其预测效果也优于MASTeR的单帧运动估计.

图 17展示了k-t FOCUSS ME/MC利用参考帧的运动估计的预测效果(ME based on ref Img), 同时也展示了MASTeR利用相邻的前一帧的运动估计的预测效果(F-Img)或利用相邻的后一帧的运动估计的预测效果(B-Img), 以及联合相邻帧运动估计的预测效果 $(0.5F+0.5B)$ .从图中可以看出, 联合相邻帧运动估计的预测效果为最优.

图 17 运动估计效果 Figure 17 The quality of motion estimation
5.2 不同运动估计方案比较

JAFP依赖运动估计方法计算预测图像序列 $\widehat{{M}}_{{F}}$ $\widehat{{M}}_{{B}}$ , 为了比较不同的运动估计算法对重建结果的影响, 通过实验比较了在不同运动估计方案下JAFP的重建误差.实验结果如图 18, 图 18(b)图 18(c)图 18(d)分别为基于CWT、OBMC和OP三种运动估计方案的JAFP重建ROI及其误差图, 从图中可以看出, 基于CWT运动估计的重建方案重建的ROI误差相比另外两种方法要少.因此, 实验中选用能抑制噪声影响的CWT运动估计方案.

图 18 在8倍下采样时, 不同运动估计方法重建的ROI及其误差图(从左至右依次为心脏电影成像序列的第7帧、第12帧和第24帧ROI图像 Figure 18 The reconstructed ROI and errors with different motion estimation plan (the reconstructed ROI of the 7th, 12th, 24th frames with acceleration factor of 8)
5.3 算法收敛性分析

为了进一步说明该算法的收敛性, 对数据集2计算过程的收敛速率进行了初步定量分析.实验分别截取了该算法前7次迭代的重建图像, 如图 19所示.如图 19选取了心脏收缩阶段的第5帧, 并分别计算了前7次迭代重建的第5帧图像ROI与原图像的误差图, 从图中可以看出, 其在第5次迭代时图像误差趋于稳定.

图 19 在8倍下采样时, 各迭代次数(0 ~ 7)下重建的第5帧图像ROI及其误差图 Figure 19 The reconstructed the 5th frame ROIs and errors by 0 ~ 7 times iteration

同时, 实验对重建时前7次迭代下共25帧图像序列的SSIM和HFEN两项评价指标的均值和方差进行了对比, 如图 20所示, JAFP重建方法在第5次迭代时, 其重建结果的SSIM以及HFEN两项评价指标均收敛.不难看出, JAFP初次迭代收敛迅速, 且迭代的初始图像越精确, 收敛越快.对于阈值 $\varepsilon$ , 实验中的两个数据集都设置 $\varepsilon=0.01$ , 即前后两次迭代的相对误差小于1 %时停止迭代.由于JAFP每一次迭代都将在前次迭代结果的基础上前后向各做一次运动估计, 根据式(11)可知, 当迭代次数大于心脏运动周期时, 运动估计的结果将与前一周期相近, 不会持续提升重建质量.因此, JAFP算法的迭代次数受限于心脏运动周期, 对于阈值的选取并不敏感.

图 20 在不同迭代次数(0 ~ 7)下, JAFP方法重建综合质量对比图 Figure 20 The comparison of reconstructed quality with 0 ~ 7 iterations
5.4 算法效率比较

k-t FOCUSS ME/MC、MASTeR和JAFP三种算法均通过MATLAB实现仿真.为了比较其计算效率, 三种算法均在相同的实验环境下运行(双核3.20 GHz, 4 GB RAM, Windows 10), 其计算时间及其迭代次数如表 1所示.

表 1 运行时间对比 Table 1 The comparison of runtime

由于MASTeR在求解时需计算两个约束模型(前向图像差稀疏约束和后向图像差稀疏约束), 而JAFP联合前后帧预测中间帧图像, 仅需计算一个约束模型, 因此其计算效率较MASTeR高.

6 结束语

本文利用相邻帧之间的相似性, 联合相邻前后两帧为动态磁共振成像方法提出一种新的重建方法联合相邻帧预测(JAFP). JAFP分别计算前后两帧与中间帧之间的运动估计, 并联合预测中间帧图像.同时, 利用心脏电影成像图像序列与预测图像序列的图像差在时间域仍然具备拟周期的特性, 增加时间域的稀疏约束, 建立成像模型, 在压缩感知框架下完成重建.实验结果表明, JAFP方法重建的图像质量在SSIM和HFEN两项评价指标上均获得提升.虽然在实验中只采用了动态磁共振成像的数据, 但任意的动态成像问题均能根据相邻帧之间的相似性, 应用JAFP算法进一步提高成像质量, 从而扩展了JAFP算法的应用范围.尽管如此, 由于JAFP方法在重建时, 需要计算图像序列在时间域的稀疏变换, 因此JAFP方法难以扩展成为一种在线的重建方法.此外, 在重建时更高质量的运动估计以及多样的联合预测方式, 仍是JAFP重建方法在未来探究的重点.

参考文献
1
Donoho D L. Compressed sensing. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
2
Li Shu-Tao, Wei Dan. A survey on compressive sensing. Acta Automatica Sinica, 2009, 35(11): 1369-1377.
( 李树涛, 魏丹. 压缩传感综述. 自动化学报, 2009, 35(11): 1369-1377.)
3
Jing Nan, Bi Wei-Hong, Hu Zheng-Ping, Wang Lin. A survey on dynamic compressed sensing. Acta Automatica Sinica, 2015, 41(1): 22-37.
( 荆楠, 毕卫红, 胡正平, 王林. 动态压缩感知综述. 自动化学报, 2015, 41(1): 22-37.)
4
Peng Yi-Gang, Suo Jin-Li, Dai Qiong-Hai, Xu Wen-Li. From compressed sensing to low-rank matrix recovery:theory and applications. Acta Automatica Sinica, 2013, 39(7): 981-994.
( 彭义刚, 索津莉, 戴琼海, 徐文立. 从压缩传感到低秩矩阵恢复:理论与应用. 自动化学报, 2013, 39(7): 981-994.)
5
Qi Cong-Hui, Zhao Zhi-Qin, Xu Jing, Zhang Hai. Electromagnetic scattering and image processing of targets under complex environment based on compressive sensing method. High Power Laser and Particle Beams, 2014, 26(7): Article No. 073206.
( 齐聪慧, 赵志钦, 徐晶, 张海. 复杂场景下基于压缩感知的目标电磁散射与成像. 强激光与粒子束, 2014, 26(7): Article No. 073206.)
6
Qin Han-Lin, Han Jiao-Jiao, Yan Xiang, Zhou Hui-Xin, Li Jia, Zeng Qing-Jie. Infrared image reconstruction based on a modified block compressed sensing. High Power Laser and Particle Beams, 2014, 26(12): Article No. 121011.
( 秦翰林, 韩姣姣, 延翔, 周慧鑫, 李佳, 曾庆杰. 基于改进的分块压缩感知红外图像重建. 强激光与粒子束, 2014, 26(12): Article No. 121011.)
7
Otazo R, Candés E, Sodickson D K. Low-rank plus sparse matrix decomposition for accelerated dynamic MRI with separation of background and dynamic components. Magnetic Resonance in Medicine, 2015, 73(3): 1125-1136. DOI:10.1002/mrm.25240
8
Ren Yue-Mei, Zhang Yan-Ning, Li Ying. Advances and perspective on compressed sensing and application on image processing. Acta Automatica Sinica, 2014, 40(8): 1563-1575.
( 任越美, 张艳宁, 李映. 压缩感知及其图像处理应用研究进展与展望. 自动化学报, 2014, 40(8): 1563-1575.)
9
Lai Z Y, Qu X B, Liu Y S, Guo D, Ye J, Zhan Z F, Chen Z. Image reconstruction of compressed sensing MRI using graph-based redundant wavelet transform. Medical Image Analysis, 2016, 27: 93-104. DOI:10.1016/j.media.2015.05.012
10
Song L X, Zhang J G, Wang Q. MRI reconstruction based on three regularizations:total variation and two wavelets. Biomedical Signal Processing and Control, 2016, 30: 64-69. DOI:10.1016/j.bspc.2016.06.003
11
Ye J, Qu X B, Guo H, Liu Y S, Guo D, Chen Z. Patch-based directional redundant wavelets in compressed sensing parallel magnetic resonance imaging with radial sampling trajectory. Journal of Medical Imaging and Health Informatics, 2016, 6(2): 387-398. DOI:10.1166/jmihi.2016.1715
12
Akçakaya M, Basha T A, Pflugi S, Foppa M, Kissinger K V, Hauser T H, Nezafat R. Localized spatio-temporal constraints for accelerated CMR perfusion. Magnetic Resonance in Medicine, 2014, 72(3): 629-639. DOI:10.1002/mrm.v72.3
13
Li Q Y, Qu X B, Liu Y S, Guo D, Lai Z Y, Ye J, Chen Z. Accelerating patch-based directional wavelets with multicore parallel computing in compressed sensing MRI. Magnetic Resonance Imaging, 2015, 33(5): 649-658. DOI:10.1016/j.mri.2015.01.014
14
Majumdar A, Ward R. Learning space-time dictionaries for blind compressed sensing dynamic MRI reconstruction. In: Proceedings of the 2015 IEEE International Conference on Image Processing. Quebec, Canada: IEEE, 2015. 4550-4554
15
Li J S, Sun J Q, Song Y, Zhao J. Accelerating MRI reconstruction via three-dimensional dual-dictionary learning using CUDA. The Journal of Supercomputing, 2015, 71(7): 2381-2396. DOI:10.1007/s11227-015-1386-z
16
Tsao J, Kozerke S. MRI temporal acceleration techniques. Journal of Magnetic Resonance Imaging, 2012, 36(3): 543-560. DOI:10.1002/jmri.v36.3
17
Kim D, Dyvorne H A, Otazo R, Feng L, Sodickson D K, Lee V S. Accelerated phase-contrast cine MRI using k-t SPARSE-SENSE. Magnetic Resonance in Medicine, 2012, 67(4): 1054-1064. DOI:10.1002/mrm.23088
18
Jung H, Ye J C, Kim E Y. Improved k-t BLAST and k-t SENSE using FOCUSS. Physics in Medicine and Biology, 2007, 52(11): 3201-3226. DOI:10.1088/0031-9155/52/11/018
19
Jung H, Sung K, Nayak K S, Kim E Y, Ye J C. k-t FOCUSS:a general compressed sensing framework for high resolution dynamic MRI. Magnetic Resonance in Medicine, 2009, 61(1): 103-116. DOI:10.1002/mrm.v61:1
20
Asif M S, Hamilton L, Brummer M, Romberg J. Motion-adaptive spatio-temporal regularization for accelerated dynamic MRI. Magnetic Resonance in Medicine, 2013, 70(3): 800-812. DOI:10.1002/mrm.24524
21
Dong X L, Liu H W, Xu Y L, Yang X M. Some nonlinear conjugate gradient methods with sufficient descent condition and global convergence. Optimization Letters, 2015, 9(7): 1421-1432. DOI:10.1007/s11590-014-0836-5
22
Gorodnitsky I F, Rao B D. Sparse signal reconstruction from limited data using FOCUSS:a re-weighted minimum norm algorithm. IEEE Transactions on Signal Processing, 1997, 45(3): 600-616. DOI:10.1109/78.558475
23
Ruan Qin-Jie. The Study of Fast Dynamic Magnetic Resonance Imaging Based on Compressed Sensing Technique[Master thesis], Hangzhou Dianzi University, China, 2014.
(阮钦杰. 基于压缩感知技术的快速动态磁共振成像技术研究[硕士学位论文], 杭州电子科技大学, 中国, 2014.) http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D485080
24
Orchard M T, Sullivan G J. Overlapped block motion compensation:an estimation-theoretic approach. IEEE Transactions on Image Processing, 1994, 3(5): 693-699. DOI:10.1109/83.334974
25
Ji Li-Xia, Li Xue-Xiang. Algorithm and simulation of image stabilization for high speed moving target images based on adjacent frames compensation. Computer Science, 2014, 41(7): 310-312, 317.
( 姬莉霞, 李学相. 基于相邻帧补偿的高速运动目标图像稳像算法及仿真. 计算机科学, 2014, 41(7): 310-312, 317. DOI:10.11896/j.issn.1002-137X.2014.07.064)
26
Zhao T S, Wang J H, Wang Z, Chen C W. SSIM-based coarse-grain scalable video coding. IEEE Transactions on Broadcasting, 2015, 61(2): 210-221. DOI:10.1109/TBC.2015.2424012
27
Maligranda L. Some remarks on the triangle inequality for norms. Banach Journal of Mathematical Analysis, 2008, 2(2): 31-41. DOI:10.15352/bjma/1240336290
28
Fang Sheng. Regularization-based Parallel Magnetic Resonance Imaging Technique for Highly Accelerated Data Acquisition[Ph. D. dissertation], Tsinghua University, China, 2010.
(方晟. 基于正则化的高倍加速并行磁共振成像技术[博士学位论文], 清华大学, 中国, 2010.) http://cdmd.cnki.com.cn/Article/CDMD-10003-1011280424.htm
29
Sun S J, Gill M, Li Y F, Huang M, Byrd R A. Efficient and generalized processing of multidimensional NUS NMR data:the NESTA algorithm and comparison of regularization terms. Journal of Biomolecular NMR, 2015, 62(1): 105-117. DOI:10.1007/s10858-015-9923-x
30
Magarey J, Kingsbury N. Motion estimation using a complex-valued wavelet transform. IEEE Transactions on Signal Processing, 1998, 46(4): 1069-1084. DOI:10.1109/78.668557
31
Brox T, Bruhn A, Papenberg N, Weickert J. High accuracy optical flow estimation based on a theory for warping. In: Proceedings of the 8th European Conference on Computer Vision. Berlin Heidelberg, Germany: Springer, 2004. 25-36
32
Wang Bing, Zhao Rong-Chun, Jiang Xiao-Yue, Yu Hong-Bo. Motion estimation based on wavelet transform. Journal of Northwestern Polytechnical University, 2004, 22(4): 417-421.
( 王兵, 赵荣椿, 蒋晓悦, 俞鸿波. 一种基于小波变换的运动估计方法. 西北工业大学学报, 2004, 22(4): 417-421.)
33
Li Zhi-Qing, Shi Zhi-Ping, Li Zhi-Xin, Shi Zhong-Zhi. Sparse coding model based on structural similarity. Journal of Software, 2010, 21(10): 2410-2419.
( 李志清, 施智平, 李志欣, 史忠植. 基于结构相似度的稀疏编码模型. 软件学报, 2010, 21(10): 2410-2419.)
34
Ravishankar S, Bresler Y. MR image reconstruction from highly undersampled k-space data by dictionary learning. IEEE Transactions on Medical Imaging, 2011, 30(5): 1028-1041. DOI:10.1109/TMI.2010.2090538