基于张量积复小波紧框架的磁共振图像快速重构方法

摘要

本发明公开了基于张量积复小波紧框架的磁共振图像快速重构方法，为磁共振快速成像提供了新的研究方法。本发明利用笛卡尔采样轨迹模式对K空间数据进行欠采样，极大地提高设备扫描速度；基于张量积复小波紧框架的稀疏变换，对图像从多方向进行分解，提高了磁共振成像的精度；利用投影快速迭代软阈值算法求解图像重构中的凸优化问题，加快了磁共振图像重构的速度；且在投影快速迭代软阈值算法中，采用基于双变量收缩方法自适应地计算正则化参数，不仅舍去了对参数盲目选取的过程，更是提高了图像的重构质量。

说明书

技术领域

本发明属于磁共振成像技术领域，具体涉及一种基于张量积复小波紧框架的磁共振图像快速重构方法。

背景技术

磁共振成像(Magnetic Resonance Imaging,MRI)具有无电离辐射、多角度成像、对人体组织无损伤等特点，因此成为医学临床和医学科研中非常重要的检测手段；但是MRI存在的不足之处就是成像速慢，动态MRI时间分辨力较低；同时，单个患者扫描时间较长、成本较高，限制了MRI的进一步推广。为了加快扫描速度，提高扫描效率，近些年并行成像(Parallel Imaging,PI)技术和压缩感知理论成为研究热点。

并行磁共振成像技术，采样多通道线圈同时接受磁共振信号，减少相比使用单个线圈时所需的相位编码的次数，从而降低信号采集时间，加快成像速度。目前，有两种方法被广泛应用于从获取的部分并行数据中重构图像，一种方法是利用线圈灵敏度信息，如SENSE(Sensitivity Encoding)；另一种则利用多通道K空间数据之间的相关性，如GRAPPA(Generalized Auto-calibrating Partially Parallel Acquisitions)。然而，由于严重病态的系统矩阵和重构算法的性能局限，并行磁共振的加速因子仍然受到限制。

压缩感知理论的兴起，为磁共振快速成像提供了理论支持。压缩感知理论认为，只要信号是稀疏的，或者在某个变换域内是稀疏的，我们就可以突破奈奎斯特采样率的限制，使用较少的数据，来恢复大量的信号。这种“压缩”数据的特性与MRI相结合，能够极大地减少K空间的采样数据，从而降低扫描时间。目前，基于稀疏表示的磁共振重构算法被提出，这些算法通常在频域或者时域上来稀疏表示信号，例如小波变换，Total Variation等，最后通过反解凸优化问题来重构图像。为了得到自适应的字典来稀疏表示信号，字典学习被用来解决磁共振成像问题，如DLMRI。然而，这些算法要么是无法提供精准的基来稀疏表示信号，要么是需要花费大量的时间来求解凸优化问题，因此限制了基于压缩感知的磁共振快速成像的应用。

发明内容

本发明的目的是为克服上述方法存在的缺陷，提出一种基于张量积复小波紧框架的磁共振图像快速重构方法，该方法基于张量积复小波紧框架作为稀疏变换，从多方向上对信号进行分解，提高了磁共振成像的精度；利用投影快速迭代软阈值算法求解图像重构问题，提高了磁共振成像重构的速度；基于双变量收缩方法自适应地计算正则化参数，进一步提高了图像的重构质量。

本发明具体步骤如下：

1)利用笛卡尔采样轨迹对磁共振K空间数据进行欠采样，得到K空间欠采样数据。

2)对K空间欠采样数据进行零填充，再通过傅里叶反变换求得初始磁共振图像。

3)将初始磁共振图像，作为输入图像，利用张量积复小波紧框架进行稀疏分解，分解后得到的分解系数表达式如下：

其中α为分解系数，x为输入图像，T为转置符号；为张量积复小波紧框架；

4)采用投影快速迭代软阈值算法对磁共振图像的目标函数进行最小化求解，优化问题表述如下：

其中，N维复数向量x和M维复数向量y分别代表要重构的磁共振图像和K空间欠采样数据；N×N复数矩阵F表示离散傅里叶变换矩阵；M×N实数矩阵U表示对K空间数据进行欠采样时的欠采样矩阵，M<N；λ为正则化参数，||·||₂为二范数运算，||·||₁为一范数运算。

迭代方程如下：

μ＝x_k+F^*U^T(y-UFx_k)

其中，x_k表示迭代第k次的重构图像；c、μ和s_k+1均为中间过渡向量，F^*和U^T分别表示离散傅里叶反变换和采样矩阵的转置；为软阈值函数，用于收缩小波系数；t_k为迭代第k次的步长。

当迭代次数达到设定值时，则输出重构出的磁共振图像。

进一步，分解系数α中每一层α_q对应的λ采用双变量收缩函数来评估，具体如下：

其中，α_q是当前尺度系数，α_p是α_q的父级尺度系数，σ_α是当前子带(当前层)q内以α_q为中心的3×3窗口预测方差，噪声方差用一个中值估计器来估计，即：

其中α_h是第一层分解中高频子带系数，median为取中值操作符。

本发明是基于张量积复小波紧框架的磁共振图像快速重构方法，为磁共振快速成像提供了新的研究方法。本发明利用笛卡尔采样轨迹对K空间数据进行欠采样，极大地提高设备扫描速度；基于张量积复小波紧框架的稀疏变换，对图像从多方向进行分解，提高了磁共振成像的精度；利用投影快速迭代软阈值算法求解图像重构中的凸优化问题，加快了磁共振图像重构的速度；且在投影快速迭代软阈值算法中，采用基于双变量收缩方法自适应地计算正则化参数，进一步提高了图像的重构质量。

附图说明

图1为本发明的流程示意图。

图2为笛卡尔采样轨迹模式图。

图3(a)为肩膀磁共振图像。

图3(b)为胸部磁共振图像。

图3(c)为大脑磁共振图像。

图4(a)为在不同加速因子下，本发明与目前流行的DLMRI，PBDWS，PANO和GBRWT算法对大脑磁共振图像关于峰值信噪比PSNR的对比示意图。

图4(b)为在不同加速因子下，本发明与目前流行的DLMRI，PBDWS，PANO和GBRWT算法对大脑磁共振图像关于结构相似性测量SSIM的对比示意图。

具体实施方式

为了更加具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，本发明基于张量积复小波紧框架的磁共振图像快速重构方法，具体步骤如下：

1、基于如图2所示的笛卡尔采样轨迹对磁共振K空间数据进行欠采样，得到K空间欠采样数据，其中白色像素点值为1，表示采样点；黑色像素点值为0，表示非采样点。

2、对K空间欠采样数据进行零填充(将非采样点值设为0)，再通过傅里叶反变换求得初始磁共振图像。

3、将初始磁共振图像，作为输入图像，利用张量积复小波紧框架进行稀疏分解，分解后得到的分解系数表达式如下：

其中α为分解系数(也称为小波系数)，x为输入图像，T为转置符号；为张量积复小波紧框架，构建方法采用现有技术，如2014年7月发表于期刊SIAM Journal on ImagingSciences上的Han B,Zhao Z.Tensor product complex tight framelets withincreasing directionality。

4、采用投影快速迭代软阈值算法(pFISTA)对磁共振图像的目标函数进行最小化求解，优化问题表述如下：

其中N维复数向量x和M维复数向量y分别代表要重构的磁共振图像和K空间欠采样的数据；N×N复数矩阵F表示离散傅里叶变换矩阵；M×N(M<N)实数矩阵U表示对K空间数据进行欠采样时的欠采样矩阵；λ为正则化参数，||·||₂为二范数运算，||·||₁为一范数运算。

迭代方程如下：

μ＝x_k+F^*U^T(y-UFx_k)

其中，x_k表示迭代第k次的重构图像；c、μ和s_k+1均为中间过渡向量，F^*和U^T分别表示离散傅里叶反变换和采样矩阵的转置；为软阈值函数，用于收缩小波系数；t_k为迭代第k次的步长。

其中，利用双变量收缩函数来评估分解系数α中每一层α_q对应的λ，具体如下：

其中，α_q是当前尺度系数，α_p是α_q的父级尺度系数，σ_α是当前子带(当前层)q内以α_q为中心的3×3窗口预测方差，噪声方差用一个中值估计器来估计，即：

其中α_h是第一层分解中高频子带系数，median为取中值操作符。

当迭代次数达到设定值时，则输出重构出的磁共振图像。

为了验证基于张量积复小波紧框架的磁共振图像快速重构方法(TPCTF-BS)的重构性能，采用三张全采样的磁共振图像进行模拟实验，如图3(a)、3(b)、3(c)所示，它们分别是肩膀，胸部和大脑的磁共振图像数据。图2是加速因子为4(相当于采样率为25％)的笛卡尔采样矩阵，它们被用来欠采样K空间上的数据。将本发明方法与目前流行的磁共振重构算法进行比较，主要包括DLMRI，PBDWS，PANO和GBRWT，用峰值信噪比(PSNR)、结构相似度(SSIM)作为客观评价指标来衡量不同算法的重构效果，以验证本发明方法的重构性能。两种客观评价指标的定义分别为：

其中，x为重构后的图像，为全采样下的图像，W表示图像x及的宽度，H表示图像x及的高度，μ_x与分别为x和的均值，与分别为x和的方差，为两幅图像的协方差，C₁与C₂是用来维持稳定的常数。

表1和表2记录了在笛卡尔欠采样矩阵下，欠采样25％的K空间数据，利用不同算法重构磁共振图像的PSNR和SSIM，可以发现，本发明TPCTF-BS具有最好的重建质量。

表3是磁共振图像在不同算法下的平均重构时间，可以发现，本发明TPCTF-BS拥有最快的重构速度。

图4(a)和图4(b)分别是大脑磁共振图像在不同加速因子下重构出的磁共振图像的PSNR和SSIM，可以看出，随着加速因子的增加，即采样K空间的数据的减少，所有算法的重建质量下降，但是在同等的加速因子下，本发明TPCTF-BS仍然具有最高的PSNR和SSIM值，表明该方法在任何重构条件下具有适用性。

表1

图像算法DLMRIPBDWSPANOGBRWTTPCTF-BS肩膀图像35.6937.4937.3736.6437.69胸部图像26.3626.3326.6926.0328.50大脑图像28.5530.8429.3528.9231.33

表2

图像算法DLMRIPBDWSPANOGBRWTTPCTF-BS肩膀图像0.92910.94330.94360.93760.9484胸部图像0.73900.82850.78290.79980.8620大脑图像0.84880.90790.87040.88900.9145

表3

算法DLMRIPBDWSPANOGBRWTTPCTF-BS重构时间(s)208.792.654.695.645.8

以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步的详细说明，应当理解，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围。特别指出，对于本领域技术人员来说，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。