加速动态磁共振重建的低秩系数分解模型

Low-Rank and Sparse Decomposition Model for Accelerating Dynamic MRI Reconstruction，本文源自2017年8月8日的Journal of Healthcare Engineering,截图如下
文章可以去这里下载，截图如下

介绍

理论背景

k-t空间的动态磁共振数据采集可以用下式表示

$$\begin{equation} y(k, t) = \int{x(r, t)\exp(-2\pi jk \cdot r)dr} + n(k, t). \tag{1} \end{equation}$$
其中 $y(k, t)$表示观测到的k-t空间信号， $x(r, t)$表示期望的动态图像序列， $n(k, t)$为观测噪声，它可以用一个额外的白噪声的高斯分布[16][17]来进行合理的模拟.

本文中，问题的解决方案是为了从降采样的观测数据$y(k, t)$中找出磁共振图像$x(r, t)$的最近似表示，$(1)$式可以转换为一个逆问题，并且可以重写成向量形式[18]。

$$\begin{equation} \boldsymbol{Y} = \boldsymbol{R} F \boldsymbol{X} + \boldsymbol{n}. \tag{2} \end{equation}$$
其中 $\boldsymbol{Y}=[y_1| \ldots|y_T]$， $\boldsymbol{X} = [x_1|\ldots|x_T]$， $\boldsymbol{n} = [n_1|\ldots|n_T]$， $T$是框架中的元素个数，$F$为傅里叶变换算子，观测矩阵 $R$为降采样掩模，它被用于k空间。

基于压缩感知的磁共振图像重建

压缩感知方法[5,19]被提出用于从部分采样的k空间数据$Y$中重建出磁共振图像 $X$,这一过程是利用稀疏变换和凸优化算法来实习的。如果我们找到满足式$(2)$的稀疏向量，那么问题将得到解决。

$$\begin{equation} \min_{X} \|D\boldsymbol{X}\|_0 \quad s.t. \quad \|\boldsymbol{Y} - \boldsymbol{R} F \boldsymbol{X} \|_2 \le \epsilon \tag{3} \end{equation}$$
其中 $\|\cdot\|$是 $l_0$范数，它用来计算向量中非零元素的个数， $D$是稀疏变换或者字典，$\epsilon$是一个较小的常数。不幸的是， $(3)$是NP难问题，这需要通过蛮力搜索来解决。压缩感知理论[8]提供了凸松弛方法，指的是 $(3)$式中的 $l_0$范数可以用 $l_0$范数最小化来代替，
$$\begin{equation} \min_X\|D\boldsymbol{X}\|_1 \quad s.t. \quad \|\boldsymbol{Y} - \boldsymbol{R} F \boldsymbol{X}\|_2 \le \epsilon, \tag{4} \end{equation}$$
其中 $\|\cdots\|$是 $l_1$范数，表示向量的绝对值之和。

低秩稀疏分解

基于压缩感知的技术已经完全成功地应用于磁共振图像重建，它利用了图像在变换域
的稀疏性。然而， 压缩感知的性能主要依赖于特定的字典或者稀疏算子，这限制了最大可实现的加速度。因此，一些研究者尝试研究一些新方法来重建磁共振图像[20-24]。这些方法中，低秩矩阵恢复是医学图像处理中最流行的一项技术。

低秩的基本假设和[18]一样，也即，图像$X$是同时稀疏（在图像域内）并且低秩的。现在的问题是从给定的少量$k$-空间数据$Y$中恢复$X$, 少量数据是相对矩阵中的元素个数而言的。我们假定矩阵的近似秩为$r$，图像的大小为$M\times N$。当矩阵$X$是低秩时，它只有$r(M+N-r)$个自由度而不是$MN$个， 这就有可能从少量的采样点中通过解决秩最小化问题来恢复矩阵$X$，

$$\begin{equation} \min \space rank(\boldsymbol{x}) \quad s.t. \|\boldsymbol{Y} - M（\boldsymbol{X})\|_2 \le \epsilon \tag{5} \end{equation}$$
然而，秩最小化问题，也即解决式$(5)$， 是组合问题也是NP难问题[25]。因此，凸松弛经常用于使得最小化问题更容易处理。

$$\begin{equation} \min_{\boldsymbol{X}} \|\boldsymbol{X}\|_* \quad s.t. \|\boldsymbol{Y} - M（\boldsymbol{X})\|_2 \le \epsilon \tag{6} \end{equation}$$
其中 $M$表示任意的线性算子，$\| \boldsymbol{X} \|_*$是核范数，可以定义为
$$\begin{equation} \| \boldsymbol{X} \|_* = \sum_{i=1}^{k} \sigma_i, \tag{7} \end{equation}$$
其中 $\sigma_1, \sigma_2, \ldots, \sigma_r$是 $\boldsymbol{X}$的奇异值， $r$是$\boldsymbol{X}$的秩。

为了从给定的$\boldsymbol{Y}$恢复$\boldsymbol{X}$，\boldsymbol{X}可以分解为地址矩阵\boldsymbol{A}和稀疏矩阵\boldsymbol{E}的叠加。

$$\begin{equation} \boldsymbol{X} = \boldsymbol{A}+ \boldsymbol{E} \tag{8} \end{equation}$$
$\boldsymbol{X}$可以有下列优化问题的解决方案中恢复：
$$\begin{equation} \min_{A,E}\|\boldsymbol{A}\|_* + \gamma\|\boldsymbol{E}\|\quad s.t. \quad \|\boldsymbol{Y} - M(\boldsymbol{A}+ \boldsymbol{E})\|_2 \le \epsilon \tag{9} \end{equation}$$
其中低秩矩阵 $\boldsymbol{A}$有非零奇异值并且表示背景元素，稀疏矩阵 $\boldsymbol{E}$有非零元素值，它对应于变化， $\gamma$是可调整参数，用于平衡 $l_1$范数相对核范数的贡献。

提出方法

在主成分追踪模型[26]中，为了解决式$(9)$，它可以用正则化而不是严格的约束[15]来形成一个优化问题。因此，$(9)$可以转换为

$$\begin{equation} \min_{L,S} \|\boldsymbol{Y} - M(\boldsymbol{A} + \boldsymbol{E})\|_F^2 + \lambda_L\|A\|_* + \lambda_S\|T\boldsymbol{E}\|_1, \tag{10} \end{equation}$$
其中参数 $\lambda_S$和 $\lambda_L$用于平衡数据一致性， $T$是稀疏变换基。

方程$(10)$是一个正则化的主成分追踪（RPCA）问题，它包括最小化核范数和 $l_1$范数的组合。Otazo et al. 研究[15]采用了迭代阈值方案去解决 $(10)$；然而，迭代阈值收敛慢。因此，我们提出了一种不严格的增强拉格朗日乘子算法来解决RPCA问题[27]。 根据式 $(6)$中的约束

$$\begin{equation} M_H y = \boldsymbol{A}^{(n)} + \boldsymbol{E}^{(n)} = \boldsymbol{X}^{(n)}, \tag{11} \end{equation}$$
其中 $M_H$是对偶算子， $\boldsymbol{X}^{(n)}$包含观测噪声， $\boldsymbol{A}^{(n)}$和 $\boldsymbol{E}^{(n)}$分别是低秩元素和稀疏元素。我们应用IALM方法来解决下列优化问题
$$\begin{equation} \begin{aligned}[center] \min_{L,S} \|\boldsymbol{A}^{(n)}\|_* + \lambda\|T \boldsymbol{E}^{(n)}\|_1 + <\mathcal{L}, M^H y - \boldsymbol{A}^{(n)} - \boldsymbol{E}^{(n)}> \\ + \frac{\mu}{2} \|M^H y - \boldsymbol{A}^{(n)} - \boldsymbol{E}^{(n)}\|_F^2 \tag{12} \end{aligned} \end{equation}$$

其中$\mathcal{L}$是拉格朗日乘子，用来消除等式约束，$\mu$是一个很小的正标量。条件$\sum_{k=1}^{+ \infty} \mu_{k}^{-1} = +\infty$表明$\mu_k$不会增长太快。IALM方法用于解决RPCA问题的步骤可以参考算法1.

对算法1，如果$\{\mu_k\}$是非递减并且$\sum_{k=1}^{+ \infty} \mu_{k}^{-1} = +\infty$，那么在RPCA问题中$(\boldsymbol{A}_{k}, \boldsymbol{E}_{k})$收敛于一个优化解$(\boldsymbol{A}^{*}, \boldsymbol{E}^{*})$。无约束$\{\mu_k\}$的优点是可行性条件$\boldsymbol{A}_{k} + \boldsymbol{E}_{k} = \boldsymbol{X}$可以快速逼近，因为$\boldsymbol{X} - \boldsymbol{A}_{k} - \boldsymbol{E}_{k} = \mathcal{L}_{k} - \mathcal{L}_{k-1} / \mu_{k-1}$和$\mathcal{L}_{k}$是受约束的。在算法1中，单一阈值算子[28]可以定义为

$$\begin{equation} SVT_{\lambda}(\boldsymbol{D}) = U \Lambda_{\lambda}(\Sigma)V^H, \tag{13} \end{equation}$$
其中 $\boldsymbol{D} = U \Sigma V^H$是 $\boldsymbol{D}$的任意一个奇异值分解。 $\Lambda_\lambda(\Sigma)$是软阈值算子，可以定义为
$$\begin{equation} \Lambda_\lambda (\Sigma) = \frac{x}{|x|} \max{(|x| - \lambda, 0)}. \tag{14} \end{equation}$$

实验结果与讨论

实验是在MATLAB V7.14.0 (R2012a)， Intel Core i7-2640 M CPU, 4.0 GB 内存，64-bit Win7操作系统上运行的。提出的算法采用两个心脏电影数据集， 来进行实验验证。第一个数据集是从Bio Imaging and Signal Processing Lab上获取的，
这里包含了$n_t = 25$个时间框架，大小为$n_x = n_y = 256$，视野为$345 \times 270mm^2$，层厚为$10mm$。

第二个数据集可以从的Dr. Caballero网站上获取
ta是由Caballero et al. [12]引入的，相关的图像参数如下：图像矩阵大小为$256 \times 256(n_x \times n_y)$，时间框架数为30($n_t$), $FOV = 320 \times 320mm^2$, 层厚为$10mm$.

两种广泛使用的采样轨迹，笛卡尔和径向采样策略被用于$k$-空间的MR数据集采样。图1显示了研究中的采样掩模以及他们在一个时间框架上的幅度。

我们从重建精度和重建速度量方面，比较了提出的算法和k-t SLR[10]以及k-t RPCA[16]。定量图像质量评估可以采用峰值信噪比（PSNR）和结构相似性索引（SSIM）这两个准则。PSNR用于估计重建图像和全采样图像之间的差异，可以定义为

$$\begin{equation} PSNR = -10log_{10} \Bigg ( \frac{\|\hat{\boldsymbol{X}} - \boldsymbol{X}\|_F^2}{\|\boldsymbol{X}\|_F^2} \Bigg) \tag{15} \end{equation}$$
其中$\hat{\boldsymbol{X}}$是重建图像，$\boldsymbol{X}$表示全采样图像。

SSIM是测量重建图像和全采样图像之间相似性的一种新方法。对每一个时间框架$\{x_n\}_{n=1}^{n_t}$，我们采用SSIM去测量重建图像和全采样图像之间的差异；同一重建图像$x_{Rec}$和全采样图像$x_F$之间的SSIM索引可以用下式来评估：

$$\begin{equation} SSIM(x, y) = \frac{(2 \mu_{x_{R}} \mu_{x_{F}} + c_1)(2\sigma_{x_{R} x_{F}} + c_2)}{(\mu_{x_{R}}^2 + \mu_{x_{F}}^2) + c_1)(\sigma_{x_{R}}^2 + \sigma_{x_{F}}^2 + c_2)}, \tag{16} \end{equation}$$
其中 $\mu_{x_{R}}$和 $\mu_{x_{F}}$是重建图像 $x_{Rec}$和全采样图像 $x_F$的平均灰度值 $\sigma_F$和 $\sigma_{x_R}$是重建图像 $x_{Rec}$和全采样图像 $x_F$的标准差， $\sigma_{x_{R} x_{F}}$是重建图像 $x_{Rec}$和全采样图像 $x_F$的协方差， $c_1=(K_1 L)^2$， $c_2=(K_2 L)^2$是常数，其中 $L$是动态范围为$[0, 255]$的8位灰度图像。Wang et al. [29]文中建议的参数值为 $K_1=0.01$， $k_2=0.03$/

k-t SLR算法结合了TV和非凸Schatten p范数，其中$p=0.01$; 一些参数选为基于公开的软件包中的推荐值（Schatten惩罚项参数$\beta_1 = 10^{-9}$, TV范数中则为$\beta_2 = 10^-2$， 最大的内部迭代数为50， 最大的外部迭代数为9）。在k-t RPCA算法中，两个正则化参数为$\mu=200$，$\rho = 1.5$，分别对应正则化项和分解项。

类似地，本文提出的这一方法需要三个参数$\lambda$，$\rho$和$\mu$的具体化。我们设定$\mu_0=1.5/\|\boldsymbol{X}\|_2$， $\rho = 1.2$。我们采取$\|\boldsymbol{X}- \boldsymbol{A}_{k}- \boldsymbol{E}_{k}\|_F / \|\boldsymbol{X}\|_F < 10^{-7}$作为算法1的终止条件。在[16]作者的建议下，我们选择了一个固定的加权参数$\lambda = \max(n_x * n_y,n_t)^{{-1}/{2}}$，提出的算法经实验验证有效，采用了全采样的心脏电影（前文提到的两个数据集）以及两种不同的采样轨迹。

为了仿真$k$-空间的加速，全采样$k$-空间数据是通过变密度（采样因子）随机采样方式人为降采样来实现的。为了测试所提出的方法的鲁棒性，两个数据集的$k$-空间数据都被固定标准差$\sigma=15$的高斯白噪声污染了。首先，在第一个心脏数据集上进行测试，采用了不同的采样模型，同时采用了不同的采样率。图2中展示了可视化质量比较，它比较了提出的方法（算法1）和其他算法的重建结果。笛卡尔采样掩模的加速因子近似为4（大约是获取数据的25%)。图3展示了笛卡尔采样和伪径迹采样方式重建结果的PSNR，其中采样方式是关于采样因子的函数。可以看出在伪径迹采样方式情况下，本文提出的算法的性能由于其他两个方法。但是在低采样率笛卡尔采样方式下，本文算法的性能稍微差于k-t SLR方法。此外，图4展示了在笛卡尔和径迹两种采样方式在同一采样因子（加速因子近似为6，大约16.4%的获取数据）下，不同时间框架的SSIM。实验结果表明提出的方法就SSIM而言，实现了较优越的结果。当采用伪径迹采样方式而不是笛卡尔采样方式时， 并且因此我们的方法的优势相对而言更明显。

更多地，我们用同样的方法在第二个心脏数据集上测试我们提出的方法。图5提供了可视化估计

总结

本文中，我们提出了一种快速算法（IALM），用于解决RPCA优化问题，从而从高度降采样 k <script type="math/tex" id="MathJax-Element-141">k</script>-空间数据中恢复出dMRI序列。我们提出的算法具有一般性，可用于将动态磁共振数据分离成低秩部分和稀疏部分。并且这一算法可以从部分观测数据同时重建和分离动态MR数据。在心脏数据集上进行的实验验证了算法的有效性和高效性， 相比最先进的算法而言。