基于自适应正交字典学习的动态磁共振并行重建方法与流程

本发明属于医学图像重建技术领域，具体是涉及一种基于自适应正交字典学习的动态磁共振并行重建方法。

背景技术：

磁共振成像(magneticresonanceimaging，mri)技术具有无伤口、无辐射、分辨率高和可多维成像等优点，不仅可显示人体组织的解剖信息，而且可显示其功能信息。mri被广泛的应用于临床医学各个系统，是继ct以后的又一重要临床检测方法。但mr成像速度慢是其一大缺点，尤其是动态磁共振成像(dynamicmagneticresonanceimaging，dmri)，需要在较短时间内获得高时空分辨率的mri图像序列，目前是医学界的一个难题。过长的扫描时间加上病人的器官运动(如呼吸，吞咽等)，会导致成像模糊，同时也无法满足动态实时成像和功能成像的需求。在k空间进行降采样是提高成像速度的一种方法，但如果直接从k空间逆里叶逆变重建图像，根据奈奎斯特采样定理，会导致重建图像产生混叠效应。动态磁共振成像数据在时空域的具有很强的稀疏特性，使得压缩感知(compressivesensing，cs)技术被广泛应用到mr图像重建当中来。cs理论指出，如果一个信号在变换域是稀疏的，且变换基和测量矩阵是不相关的(又称有限等距性质，rip)，则采用cs方法可以从降采样(远小于奈奎斯特采样率条件下)的数据样本当中，通过非线性重建算法完美重建该信号。

dmri重建方法可以分为在线模式和离线模式。采用离线模式重建时，需要在重建之前获得整个dmri序列的采样数据。常见的离线重建方法有运动校正，字典学习，和低秩近似等。这些方法充分利用整个dmri序列稀疏特性进行高精度重建，其缺点是重建速度较慢，并且需要等待较长的扫描时间。采用在线方法重建时，每一帧的重建仅仅跟之前的帧有关，可以边扫描边重建，节省了等待时间，但由于缺乏整个序列的完整信息，重建精度无法保障，同时对重建速度也提出更高的要求。

在线重建有两种常用的方案：串行和并行。串行方案通常利用图像或变换域中相邻帧之间的稀疏性，这也是大多数现有在线方法中常用的策略。然而，这些方法往往会导致累积误差。chenc等人采用了一种新的动态全变差(dtv)的并行重建方案来解决误差累积问题。该方法采用第一帧高精度采样作为参考帧，其余的帧与其逐一比较并行重建。该方法每次只能利用两帧之间的稀疏作为重建的先验知识，导致比离线方法精度低。

技术实现要素：

为解决现有技术中存在的问题，本发明提出了一种基于自适应正交字典学习的实时动态磁共振并行重建方法。

为达到上述目的，本发明的技术方案如下：

一种基于自适应正交字典学习的动态磁共振并行重建方法，包括如下步骤：

s1：输入原始dmri序列x，输入测量值y，所述测量值y为k-t空间的欠采样数据，采用伪随机射线欠采样模式，输入算法的第一循环迭代次数outloop，输入算法的第二循环迭代次数innerloop，输入字典学习参数；

s2：初始化，将重建图像初始值设为这里为第k次迭代后重建的mr子序列图像，xzf为k-t空间欠采样后零填充数据，初始化字典d，所述字典d为dct字典；

s3：迭代更新，

fori＝1:outloop

forj＝innerloop

更新自适应字典d；

更新图像块稀疏表示系数αi；

更新重建图像的频域值逆傅里叶变换得到第j个重构子序列

end

输出重建子序列图像xs(j)，等待下一个子序列xs(j+1)；

end

s4：将各子序列图像xs(j)重新组合成重建后的dmri序列。

作为上述技术方案的优选，所述步骤s3中：

自适应字典d和稀疏表示系数αi的求解包括如下步骤：

s310：给定一个dmri序列表示为x，及其k-t空间的欠采数据为y，压缩感知dmri重建问题可以归结为以下l0范数最小化问题：

其中，fu＝diag[fu(1),fu(2),...,fu(nt)]为k-t空间的采样矩阵，fu(t)＝f2dpt，其中f2d为二维傅里叶变换算子，pt是第t帧欠采样矩阵，y表示欠采样k域，||x||0是x的l0范数，ψ是稀疏变换矩阵，λ是与采样噪声相关的常数；

s311：利用压缩感知系统采集到的测量值作为图像观测样本数据来训练字典，将待处理的图像分为重叠的小块，来代替整幅图像进行稀疏表示，字典学习问题可描述为：

其中，x为待重建图像序列，d为过完备字典，ri为重叠取块的算子，αi为x的稀疏表示系数，γ＝{α1,...,αi}为稀疏表示系数αi的集合，t0表示稀疏度阈值常数，s.t.为满足约束条件的意思，||αi||0≤t0为约束条件，为对任意的i，其中第一项和稀疏约束条件||αi||0≤t0确保了过完备字典对每个图像块的最优稀疏逼近，第二项为数据保真项，变量ν为常数，与k-t空间采样时叠加的高斯白噪声的标准差σ有关；

s312：字典学习过程中加入正交限制约束d^td＝i，步骤s311中的目标函数变为：

其中，这里ν和β为正则化参数，目的是减少后面两项的贡献值，防止方程产生过拟合，i为单位对角阵；

s313：各子序列图像xs数据保持不变，将问题转化为求解步骤s312公式中d和αi最优解的子问题：

s314：进行第一次迭代时，为对应子序列图像xs在k空间欠采样后直接进行零填充得到的图像数据，首先将采样子序列图像xs重叠采样间隔为1的三维重叠分块，随机选取部分量图像块，用dct字典作为初始字典，固定d，采用下述公式算法更新稀疏表示系数αi，

该问题的求解可采用硬阈值函数，具体体现为：

其中t(g)为硬阈值函数，

s315：更新完稀疏系数αi以后，固定αi，采用奇异值分解的方法更新字典d，字典更新问题可以转化为：

使得d^td＝i

这里，

其中，x＝{x1,x2,lxm}∈r^n×m为图像块矩阵，v＝{v1,v2,lvm}∈r^k×m为稀疏系数矩阵，tr(g)为矩阵求迹运算，则字典更新问题转变为：

该问题由奇异值分解算法来实现：

xv^t＝p∑q^t,d^k+1＝pq^t

这个是典型的矩阵的svd分解，σ是一个m×n的矩阵，除了主对角线上的元素以外全为0，主对角线上的每个元素都称为奇异值，p和q都是酉矩阵，即满足p^tp＝i，q^tq＝i。

作为上述技术方案的优选，所述步骤s3中：

更新重建图像的频域值逆傅里叶变换得到第j个重构子序列具体包括如下步骤：

s321：由压缩感知字典学习重建模型可得：

s322：步骤s321的公式中d与αi固定不变，图像重建子问题变成一个普通的最小二乘法问题，对唯一的变量xs求导并令其等于0得到：

其中：为fu的共轭转置矩阵；

s323：对步骤s322中的公式两边进行傅里叶变换得到：

i代表单位对角阵，n为任意一个像素被不同三维图像小块包含的个数，当分块间隔取最小值1时，n即为图像分块的向量维数，为降采样零填充的k空间数据，可以表示为是一个p×p的对角矩阵，p是整个子序列图像排成向量的维数；

上式可简化为：

其中，为待重建序列在k空间对应位置(kx,ky)处的取值，ω是采样矩阵中取值为1的位置的集合，λ＝q/σ是由k空间采样噪声标准差决定，σ为噪声方差，q为含噪采样条件下的可调参数，在无噪声采样条件下λ为无穷大，采样点的重建信号可以直接令

本发明的有益效果在于：本发明提出的一种基于自适应正交字典学习的实时动态磁共振并行重建方法，将原本在离线模式下运行的相对较慢的的字典学习算法应用到在线模式当中来，以高精度采样的第一帧为参考，实现对任意n个相邻帧mr图像的实时在线重建，以三维图像小块作为重建对象，采用正交字典作为稀疏约束条件和奇异值分解(singularvaluedecomposition,svd)算法提高重建速度和精度。

附图说明：

以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范围。其中：

图1为本发明一个实施例的一种并行自适应字典学习dmri重建原理图；

图2为本发明一个实施例的第一帧采样矩阵和其他帧采用矩阵示意图；

图3为本发明一个实施例的迭代次数与psnr关系示意图；

图4为本发明一个实施例的重建前后图像对比示意图；

图5为本发明一个实施例的各算法重建mr图像和误差图像比较示意图；

图6为本发明一个实施例的重建均方根误差与q的关系示意图。

具体实施方式：

本发明的一种基于自适应正交字典学习的动态磁共振并行重建方法，包括如下步骤：

s1：输入原始dmri序列x，输入测量值y为k-t空间的欠采样数据，采用伪随机射线欠采样模式(见附图2)，输入算法的第一循环迭代次数outloop，输入算法的第二循环迭代次数innerloop，输入字典学习参数；

s2：初始化，将重建图像初始值设为这里为第k次迭代后重建的mr子序列图像，xzf为k-t空间欠采样后零填充数据，初始化字典d为dct字典。

s3：迭代更新，

fori＝1:outloop

forj＝1：innerloop

更新自适应字典d；

更新图像块稀疏表示系数αi；

更新重建图像的频域值逆傅里叶变换得到第j个重构子序列

end

输出重建子序列图像xs(j)，等待下一个子序列xs(j+1)；

end

s4：将各子序列图像xs(j)重新组合成重建后的dmri序列。

自适应字典d和稀疏表示系数αi的求解包括如下步骤：

s310：给定一个dmri序列表示为x，及其k-t空间的欠采数据为y，压缩感知dmri重建问题可以归结为以下l0范数最小化问题：

其中，fu＝diag[fu(1),fu(2),...,fu(nt)]为k-t空间的采样矩阵，fu(t)＝f2dpt，其中f2d为二维傅里叶变换算子，pt是第t帧欠采样矩阵，||x||0是x的l0范数，ψ是稀疏变换矩阵，λ是与采样噪声相关的常数。

mri采集k域(傅里叶变换域)信号，再通过逆傅立叶变换重建空域信号。傅立叶变换是一种线性变换，要求k域信号数必须等于图像域的像素数。2007年lustig等人提出了压缩感知磁共振成像(csmri)的概念，以远低于奈奎斯特频率的采样频率对k空间信号进行采集，从而极大降低成像时间，并通过非线性重建算法完美重建原始图像。

csmri需要满足3个条件：

(1)稀疏矩阵ψ，使得原始mr图像在稀疏域中仅有少量非零系数；

(2)矩阵fuψ-1满足rip条件，一般采用随机或伪随机采样矩阵，降低采样数据的相关性。

(3)非线性重建算法，能快速准确重建原始图像。

s311：过完备字典在信号的稀疏表示当中有着广泛的应用。各种基于图像库训练的非自适应图像重建方法因其适应性较差、重建效果一般等缺点逐渐被基于自适应字典学习的图像重建算法所取代。这里的自适应是指由压缩感知系统采集到的测量值作为图像观测样本数据来训练字典。由于字典学习的计算量会随着字典尺寸的增加而快速增加，我们通常将待处理的图像分为重叠的小块，来代替整幅图像进行稀疏表示。字典学习问题可描述为：

其中，x为待重建图像序列，d为过完备字典，ri为重叠取块的算子，αi为x的稀疏表示系数，γ＝{α1,...,αi}为稀疏表示系数αi的集合(，t0表示稀疏度阈值常数。s.t.为满足约束条件的意思(下同)，||αi||0≤t0为约束条件，为对任意的i。其中第一项和稀疏约束条件||αi||0≤t0确保了过完备字典对每个图像块的最优稀疏逼近。第二项为数据保真项，变量ν为常数，与k空间采样时叠加的高斯白噪声的标准差σ有关；

由于dmri动态扫描时，每一帧的主要器官相对位置近似不变，我们选取序列的第一帧图像为高精度采样(一般取大于50％的采样率)，为后续帧的重建提供高精度参考，这可由磁共振机的一次预扫描完成。第一帧与其余任意相邻的n-1帧图像组成一系列n帧的子序列图像xs(j)，以3帧为例，如图1所示，xs(j)＝[x1,x2j,x2j+1]，j＝1,2,3…，x2j表示原dmri序列x中的第2j帧图像。只要相邻两帧采样完毕，即可进行实时并行重建，无需等待后续帧采样完成。每个子序列的重建与其他子序列不相关。

正交字典能够最大限度的保证字典原子之间的不相关性，从而保证了稀疏编码过程中的效率最大化。本文提出一种基于自适应正交字典的dmri并行重建算法，旨在缩短扫描等待时间，提升重建速度，提高图像重建的质量。

s312：字典学习过程中加入正交限制约束d^td＝i，步骤s311中的目标函数变为：

其中，这里ν和β为正则化参数，目的是减少后面两项的贡献值，防止方程产生过拟合，i为单位对角阵(下同)。

上式中提出的重建模型可分解为求解自适应正交字典学习、图像稀疏系数的稀疏优化子问题和图像重建子问题，实现对目标函数最优解的求解。

s313：各子序列图像xs数据保持不变，将问题转化为求解步骤s312公式中d和αi最优解的子问题：

s314：进行第一次迭代时，为对应子序列图像xs在k空间欠采样后直接进行零填充得到的图像数据，首先将采样子序列图像xs重叠采样间隔为1的三维重叠分块，随机选取部分量图像块，用dct字典作为初始字典，固定d，采用下述公式算法更新稀疏表示系数αi，

该问题的求解可采用硬阈值函数，具体体现为：

其中t(g)为硬阈值函数。

s315：更新完稀疏系数αi以后(上标k+1就是代表更新后的αi)，固定稀疏系数不变，采用奇异值分解的方法更新字典d，字典更新问题可以转化为：

这里，

其中，x＝{x1,x2,lxm}∈r^n×m为图像块矩阵，v＝{v1,v2,lvm}∈r^k×m为稀疏系数矩阵，tr(g)为矩阵求迹运算，则字典更新问题转变为：

该问题由奇异值分解算法来实现：

xv^t＝p∑q^t,d^k+1＝pq^t。

这个是典型的矩阵的svd分解，σ是一个m×n的矩阵，除了主对角线上的元素以外全为0，主对角线上的每个元素都称为奇异值，p和q都是酉矩阵，即满足p^tp＝i，q^tq＝i。

图像重建子问题：更新重建图像的频域值逆傅里叶变换得到第j个重构子序列具体包括如下步骤：

s321：

这个公式是经典的压缩感知字典学习重建模型，上标代表第k+1次迭代更新的意思，参数指标参考s311。

s322：步骤s321的公式中d与αi固定不变，图像重建子问题变成一个普通的最小二乘法问题，对唯一的变量xs求导并令其等于0得到：

其中：为fu的共轭转置矩阵。

s323：直接对上式进行求解，计算量是相当大的，因为要对一个p×p的矩阵进行求逆才可得到重建结果，这使得求解此问题的时间复杂度达到o(p³)。我们可将运算由图像空间变换到傅里叶空间进行，对步骤s322中的公式两边进行傅里叶变换得到：

n为任意一个像素被不同三维图像小块包含的个数，当分块间隔取最小值1时，n即为图像分块的向量维数，为fu的共轭转置矩阵，为降采样零填充的k空间数据，可以表示为是一个p×p的对角矩阵，p是整个子序列图像排成向量的维数；

上式可简化为：

其中，为待重建序列在k空间对应位置(kx,ky)处的取值，ω是采样矩阵中取值为1的位置的集合。λ＝q/σ是由k空间采样噪声标准差决定σ为噪声方差，q为含噪采样条件下的可调参数。在无噪声采样条件下λ为无穷大，采样点的重建信号可以直接令重建序列是在k空间更新数据后逆傅里叶变换得到的，在含噪声采样条件下，参数q对于重建性能至关重要，无噪声条件下，重建结果对q的取值并不敏感。

如图1所示，本实施例以高精度采样的第一帧为参考，实现对任意n个相邻帧mr图像的实时在线重建。以三维图像小块作为重建对象，采用正交字典作为稀疏约束条件和奇异值分解(singularvaluedecomposition,svd)算法提高重建速度和精度。

本实施例还提供一种基于自适应正交字典学习的动态磁共振并行重建方法的实验如下：

所有的实验都是在一台带有3.0ghzinteli7cpu的笔记本电脑上用matlab2013a进行仿真的。我们选取一组心肌灌注dmri序列(分辨率为192×192×30帧)用于评估算法的性能。伪射线采样矩阵用来对测试图像在k-t空间进行欠采样，模拟磁共振机加速成像。如图2所示，第一帧采样率为50％，其他帧采样率为15％。图像质量的客观评价方法采用峰值信噪比(peaksignaltonoiserate,psnr)，单位为db。psnr计算公式如下：

psnr(x,y)＝20lg(imax/rmse(x,y))

这里x和y分别为原始图像和重建图像，imax为图像像素最大值，一般为了减少矩阵运算量，要对图像进行归一化操作，则imax＝1。rmse为两幅图像的均方根误差。

实验参数设置：为了快速并行重建dmri序列，本算法将三维图像块大小设置为最小的3*3*3，经过实验测试，重建图像质量与图像块尺寸关系不大，但是重建速度会随着图像块尺寸的增加而迅速下降。字典大小设为27*108，每次在线学习的原子个数为5000。对比的算法有离线模式的k-tslr和在线模式的tv，dtv这三种，三个对比实验算法的参数均按出处论文的最优参数设置。

图3给出了本文算法重建一个dmri子序列(以3帧为例)的过程。经过14次迭代，耗时2.8秒，就将k空间欠采样后零填充的低信噪比图像(psnr＝31.1db)，重建为高信噪比的图像(psnr＝36.3db)。图4给出了重建前后的图像视觉对比，可见，重建后的图像已经完全去除了欠采样造成的伪影。

图5对比了各种算法的重建效果(无采样噪声条件下)，以列为单位，第一行为各自的重建图像，第二行为重建图像与原始图像的误差图，为凸显效果，我们将误差放大5倍进行比较。显然我们的算法误差的最小的，尤其是关键的心脏部位，视觉效果最为逼真。

在有采样噪声情况下，假设k空间叠加的高斯白噪声标准差为σ，公式(9)中λ＝q/σ的取值就与重建效果息息相关。图6给出了参数q与重建图像的rmse之间的关系。采样信噪比分别取23db(强噪声环境)、37db(中等噪声环境)和57db(低噪声环境)。可见，q＝0.005是考虑了在不同的采样噪声环境下的一个最优选择。尤其是在强噪声环境下，q＝0.005时重建效果最佳。

本文算法可适用于任意相邻的n帧子序列构成的并行重建，表格1给出了子序列帧数n与重建psnr值和算法重建时间之间的关系(20％的采样率条件下)。可见，以3帧作为一个子序列重建单元，重建的精度(psnr值)是最高的，而重建的时间又是最少的。并且各种参数下，重建精度都在可控的误差范围内，体现了本算法的稳定性。

表1重建子序列帧数与重建效果之间关系

本文算法与三个对比算法的重建psnr值对比结果如表1所示。在不同的采样率下，本文的算法均是最佳的。尤其在15％的低采样率条件下，相比其他算法有0.4db以上的性能提升。以3帧为一个dmri子序列为例，本文算法平均每1.4秒即可重建一帧图像，重建速度远远快于离线模式(k-tslr)，但慢于两种采用tv方法的算法。

表2不同采样率下重建图像的psnr值(db)

本文提出的一种基于自适应正交字典学习的动态磁共振并行重建的方法。采用在线并行重建模式，实现了对dmri图像的实时高精度重建，并对有噪声环境下的参数q作了优化。实验结果验证了该方法在重建精度和速度上的优越性。该方法为实时在线重建dmri提供了一种解决方案。

显然，上述实施例仅仅是为清楚地说明所作的举例，而并非对实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。而由此所引伸出的显而易见的变化或变动仍处于本发明创造的保护范围之中。