基于熵和几何方向的分类多字典学习磁共振图像重建方法

摘要

本发明涉及磁共振成像技术领域，具体涉及一种基于熵和几何方向的分类多字典学习磁共振图像重建方法。本发明的目的是为解决现有DLMRI算法对核磁共振图像重建时，细节部分存在不足的问题。本发明利用降采样模型获取部分K空间数据；对获取的部分K空间数据建立磁共振图像重建模型；对部分K空间数据做傅里叶逆变换得到初始图像；将初始图像分成重叠图像块；求出每一图像块的熵，将图像块样本按照熵从小到大分成四类，将后两类图像块根据几何方向进一步分类后，再对图像块样本进行字典训练，并求出对应该字典的稀疏系数，得到重建图像矩阵；再对重建图像矩阵做傅里叶变换，更新图像矩阵，对更新后的图像矩阵做傅里叶逆变换，得到重建磁共振图像。

说明书

技术领域

本发明涉及磁共振成像技术领域，具体涉及一种基于熵和几何方向的分类多字典学习磁共振图像重建方法。

背景技术

磁共振成像(Magnetic Resonance Imaging，MRI)技术是目前应用最广泛的医学成像手段之一，压缩感知(Compressed Sensing，CS)作为一种新的采样理论，为磁共振成像提供关键的理论依据，即在仅获取部分K空间数据的前提下，通过优化重建算法就可以获得满足临床诊断质量需求的磁共振图像，CS-MRI对MRI精确重建的前提是MR图像在某种变换域或字典下可以稀疏表示，以信号的稀疏先验求解图像的反问题，从而利用部分K空间数据实现MR图像的精确重建。

信号的稀疏表示，给定一组训练信号，使用一个包含该训练信号信息的字典，信号可由字典中少量原子线性组合表示，字典可以是预先定义的一组基函数，如DCT基、Gabor基等，也可以是通过某种算法学习得到的学习型字典，学习型字典能够根据训练样本自适应的构建训练字典，因此稀疏误差要比固定基字典小，S.Ravishankar等在2011年提出自字典学习算法(Dictionary Learning Magnetic resonance Imaging，DLMRI)，该算法是一种迭代算法，可以根据不同的参考图像产生自适应字典，对图像进行分块稀疏表示，因此该算法具有自适应性，将参考图像作为训练原子库进行字典训练，K-SVD依次进行字典更新，其中训练字典的每一列原子和相应的稀疏表示同时更新，固定当前训练字典，利用贪婪算法，如正交匹配追踪算法进行稀疏编码求解稀疏系数矩阵，由字典和稀疏系数矩阵重建出图像，使用DLMRI算法对MR图像重建，但重建细节部分还存在不足。

图像纹理作为一种区域特征，是对于图像各像元之间空间分布的一种描述，它能够较好的兼顾图像的宏观性质和细微结构两方面，由于纹理能充分利用图像信息，可以作为图像分类的重要依据，与其他纹理特征相比，熵能够提取图像的纹理先验信息，是图像所具有的信息量的度量，它可以表示图像纹理的非均匀程度和变化幅度的综合信息，当复杂程度高时，熵值较大，反之则较小。

Bandelet变换是一种基于边缘的图像表示方法，能自适应的跟踪图像的几何正则方向，在图像块分类中，充分考虑到图像结构中的几何正则方向信息，利用预先求解图像块的最优几何方向信息，对图像块进行更准确的分类。

发明内容

本发明的目的是为解决现有DLMRI算法对核磁共振图像重建时，细节部分存在不足的问题，提出一种基于熵和几何方向的分类多字典学习的磁共振图像重建方法，实现图像块更准确的分类，提高字典的稀疏能力，同时保留了图像的细节分量消除伪影，提高磁共振图像的重建质量。

本发明为解决上述问题采取的技术方案是：本发明的基于熵和几何方向的分类多字典学习磁共振图像重建方法，具体步骤包括：

步骤a、采用辐射型降采样模型对K空间数据进行降采样，获得部分K空间数据，并将所述部分K空间数据进行傅里叶逆变换，得到初始图像；

步骤b、根据滑动距离s，提取图像块样本，并将图像块样本转换成列向量，组成字典训练矩阵；

步骤c、计算步骤b中获取的每个图像块样本的熵，并按照熵从小到大的顺序将图像块样本分成四类，并获取每一图像块样本在初始图像的位置索引；

步骤d、对于第一类图像块样本，随机抽取若干图像块样本，进行K-SVD字典训练，得到第一冗余字典，采用正交匹配追踪算法对所述第一冗余字典进行运算，得到第一稀疏系数矩阵；

对于第二类图像块样本，随机抽取若干图像块样本，进行K-SVD字典训练，得到第二冗余字典，采用正交匹配追踪算法对所述第二冗余字典进行运算，得到第二稀疏系数矩阵；

对于第三类图像块样本和第四类图像块样本，首先计算每一图像块样本的最优几何方向，再根据最优几何方向对第三类图像块样本和第四类图像样本进一步细分类，然后对细分后的每一类图像块样本，随机抽取若干图像块样本，分别进行K-SVD字典训练，得到与细分后的每一类图像块样本对应的冗余字典，并用正交匹配追踪算法分别对细分后的每一类图像块样本得到的冗余字典进行运算，得到细分后的每一类图像块样本对应的稀疏系数矩阵；

步骤e、根据步骤d所得的每个冗余字典D_c和稀疏系数矩阵{α_ij}，重建图像块样本，并根据分类前图像块样本的位置索引把重建的图像块样本放回原来位置，具体为：令每个点的灰度值都是所有包含它的图像块样本的灰度值的平均值，得到重建图像矩阵X，再对X作傅里叶变换，得到K空间数据S，通过最小二乘算法更新K空间数据S；

步骤f、对更新后的K空间数据S做傅里叶逆变换，得到重建MR图像；

步骤g、判断是否达到设定的迭代次数，如果：

是，则输出重建MR图像，

否，则返回到步骤b。

进一步地，所述步骤b的具体方法是：将提取的图像块样本以列为单位从左至右排列，按照每列首尾依次相接的顺序，将每一个图像块样本转换成一个列向量。

进一步地，所述步骤c中，计算每个图像块样本的熵的具体方法是：将图像块样本中的复数像素点分别取模。

进一步地，所述步骤d中，在每个冗余字典D_c和稀疏系数α_ij已知的条件下，图像重建模型表示为：

其中，R_ij表示将图像分成重叠小块的算子，I表示初始图像中每一类图像块样本集合，图像块样本I_ij＝RI，代表通过K-SVD算法得到的一类冗余字典，n²表示提取的图像块样本大小，K表示冗余字典的原子个数，α_ij表示每个图像块样本在与之对应的冗余字典下的稀疏系数，λ取决于测量噪声的标准差σ，即λ＝ν/σ，其中ν是一个正常数，F_μ表示降采样算子，y表示测量数据，T₀表示稀疏表示系数中非零分量的最大数目。

进一步地，所述步骤d中最优几何方向的计算方法是：根据图像块样本的大小，进行角度预采样，构造一个与图像块样本大小相同的网格点，计算每个网格点在采样角度上的正交投影偏移，网格点对应图像块样本的灰度值按照正交投影偏移从小到大进行排序成一维数组，进行一维Haar小波变换，保留最大的M个系数，重建一维数组，重建误差最小的角度即为最优的几何方向，最优几何方向的选取通过如下模型实现：

其中，θ＝{θ₁,θ₂,...θ_d...θ_D}表示预采样角度，θ_j,d表示第j个图像块样本的第d个预采样角度，Ψ^T表示前向正交一维Haar小波变换，c_i,j(θ_j,d,M)表示M个最大小波系数，P(θ_j,d)b_j表示图像块样本重新排列的平行于θ_j,d的列向量。

进一步地，所述步骤e的具体方法是：

式(1)需要满足下面等式：

其中，H表示复数的共轭转置，β表示像素重叠次数，C为分类个数，λ为加权系数，N²为初始图像的大小，

把更新的图像样本从图像域转换到傅里叶空间，用表示规范化的傅里叶编码矩阵，即则将待重建图像的K空间数据FI带入式(3)可得：

其中，矩阵是一个由0和1组成的对角矩阵，元素1代表的是K空间的采样位置，向量表示零填充的傅里叶采样值，剩余项以及FI_t分别表示对角矩阵和可逆求解的向量，式(5)等号两边同除以β，并用λ'＝λ/β代表权重，将等号右边第一项傅里叶空间用下式表示：

综上所述，式(1)的解可表示为：

其中，表示更新后的图像矩阵，表示对图像矩阵作傅里叶变换得到的矩阵，代表的是零填充的K空间采样值，Ω是由采样值组成的K空间的子集。

有益效果：

本发明的磁共振图像重建方法首先根据图像块的熵的大小排序分类，然后对后两类图像块根据几何方向进一步细分类，对每一类新组成的图像块进行字典训练，并求出对应该字典的稀疏系数，得到重建图像矩阵；再对重建图像矩阵做傅里叶变换，更新图像矩阵，对更新后的图像矩阵做傅里叶逆变换，得到重建磁共振图像，本发明根据熵和几何方向结合的分类算法，对图像块进行多字典训练，提高字典学习能力，在低采样率的情况下，有效消除混叠伪影，提高图像重建质量。

附图说明

图1本发明的方法流程图；

图2选取的初始图像；

图3本发明的图像重建方法的降采样矩阵模型；

图4采用本发明的方法与字典学习分类方法重建图像的效果对比；

图中：(a)和(c)为采用本发明的方法重建的图像，(b)和(d)为采用字典学习分类方法DLMRI重建的图像。

具体实施方式

结合图1～图4说明本实施方式，本实施方式的一种基于熵和几何方向的分类多字典学习磁共振图像重建方法流程图如图1所示，具体步骤包括：

步骤a、采用辐射型降采样模型对K空间数据进行降采样，降采样矩阵模型如图3所示，获得部分K空间数据，并将所述部分K空间数据进行傅里叶逆变换，得到初始图像，如图2所示；

步骤b、根据滑动距离s＝2，提取图像块样本，将提取的图像块样本以列为单位从左至右排列，按照每列首尾依次相接的顺序，将每一个图像块样本转换成一个列向量，组成字典训练矩阵；

步骤c、将步骤b中获取的每个图像块样本中的复数像素点分别取模，计算每个图像块样本的熵，并按照熵从小到大的顺序将图像块样本分成四类，并获取每一图像块样本在初始图像的位置索引；

步骤d、对于第一类图像块样本，随机抽取150个图像块样本进行字典训练，得到第一冗余字典，采用正交匹配追踪算法对所述第一冗余字典进行运算，得到第一稀疏系数矩阵；

对于第二类图像块样本，随机抽取150个图像块样本进行字典训练，得到第二冗余字典，采用正交匹配追踪算法对所述第二冗余字典进行运算，得到第二稀疏系数矩阵；

对于第三类图像块样本和第四类图像块样本，首先按照如下方法计算每一图像块样本的最优几何方向：根据图像块样本的大小，进行角度预采样，构造一个与图像块样本大小相同的网格点，计算每个网格点在采样角度上的正交投影偏移，网格点对应图像块样本的灰度值按照正交投影偏移从小到大进行排序成一维数组，进行一维Haar小波变换，保留最大的M个系数，重建一维数组，重建误差最小的角度即为最优的几何方向，最优几何方向的选取通过如下模型实现：

其中，θ＝{θ₁,θ₂,...θ_d...θ_D}表示预采样角度，θ_j,d表示第j个图像块样本的第d个预采样角度，Ψ^T表示前向正交一维Haar小波变换，c_i,j(θ_j,d,M)表示M个最大小波系数，P(θ_j,d)b_j表示图像块样本重新排列的平行于θ_j,d的列向量。

然后根据最优几何方向对第三类图像块样本和第四类图像块样本进一步细分为六类，再对细分后的每一类图像块样本，随机抽取150个图像块样本进行字典训练，得到与细分后的每一类图像块样本对应的冗余字典，并用正交匹配追踪算法分别对细分后的每一类图像块样本得到的冗余字典进行运算，得到细分后的每一类图像块样本对应的稀疏系数矩阵；

在每个冗余字典D_c和稀疏系数矩阵{α_ij}已知的条件下，图像重建模型表示为：

其中，R_ij表示将图像分成重叠小块的算子，I表示初始图像中每一类图像块样本集合，图像块样本I_ij＝RI，代表通过K-SVD算法得到的一类冗余字典，n²表示提取的图像块样本大小，n＝4，K表示冗余字典的原子个数，K＝100，α_ij表示每个图像块样本在与之对应的冗余字典下的稀疏系数，λ取决于测量噪声的标准差σ，即λ＝ν/σ，其中ν是一个正常数，F_μ表示降采样算子，y表示测量数据，T₀表示稀疏表示系数中非零分量的最大数目。

步骤e、根据步骤d所得每一类的冗余字典D_c和稀疏系数矩阵{α_ij}，重建图像块样本，并根据分类前图像块样本的位置索引把重建的图像块样本放回原来位置，具体为：令每个点的灰度值都是所有包含它的图像块样本的灰度值的平均值，得到重建图像矩阵X，再对X作傅里叶变换，得到K空间数据S，通过最小二乘算法更新K空间数据S，具体计算过程是：

式(1)需要满足下面等式：

其中，H表示复数的共轭转置，β表示像素重叠次数，C为分类个数，λ为加权系数，取值为140，N×N为初始图像的大小，N＝256，把更新的图像样本从图像域转换到傅里叶空间，用表示规范化的傅里叶编码矩阵，即则将待重建图像的K空间数据FI带入式(3)可得：

其中，矩阵是一个由0和1组成的对角矩阵，元素1代表的是K空间的采样位置，向量表示零填充的傅里叶采样值，剩余项以及FI_t分别表示对角矩阵和可逆求解的向量，式(5)等号两边同除以β，并用λ'＝λ/β代表权重，将等号右边第一项傅里叶空间用下式表示：

综上所述，式(1)的解可表示为：

其中，表示更新后的图像矩阵，表示对图像矩阵作傅里叶变换得到的矩阵，代表的是零填充的K空间采样值，Ω是由采样值组成的K空间的子集。

步骤f、对更新后的K空间数据S做傅里叶逆变换，得到重建MR图像；

步骤g、判断是否达到设定的迭代次数10，如果：

是，则输出重建MR图像，

否，则返回到步骤b。

图2为选取的初始图像，图像信息为：3T西门子Trio Tim磁共振扫描仪T2加权项TR/TE＝6100/99毫秒，视野FOV＝220*220mm，切片厚度为3mm。

本发明所选取的衡量图像重建质量的参数是：峰值信噪比(PSNR)、结构相似性(SSIM)和L2范数相对误差(RLNE)；本发明首先按照图像块样本熵的从小到大将图像块样本分成四类，第三和第四图像块计算几何方向，再细分成六类，即：2熵+6几何算法，表1为2熵+6几何算法和DLMEI算法重建质量对比，由表1可知基于熵和几何方向的多分类字典学习方法在峰值信噪比、结构相似性和L2范数相对误差三方面都优于字典学习分类DLMRI方法，且重建图像耗时少。

表1 MR图像重建质量

图4中(a)和(c)为采用本发明的方法重建的图像，(b)和(d)为采用字典学习分类DLMRI方法重建的图像。可以看出，基于熵和几何方向的多分类字典学习方法比字典学习分类DLMRI方法更好的消除伪影，

本发明提供一种基于熵和几何方向的分类多字典学习磁共振图像重建方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方案，应当指出，对于本技术领域的一般技术人员，在不脱离本发明原理的前提下。在具体实施和应用范围上会有所改进之处，在这些改进之处也应视为本发明的保护范围，本实例中未明确的各组成部分均可用现有技术加以实现。