一种基于稀疏表示和非局部相似的核磁共振图像重建方法

技术领域

本发明属于图像处理技术领域，涉及一种核磁共振图像重建方法，具体是一种基于稀疏表示和非局部相似的核磁共振图像重建方法。>

背景技术

磁共振成像是一种重要的医学成像技术，在临床中有着广泛的应用。近年来的研究表明，压缩感知理论能够利用相比其它方法更少的采样数据较好重建图像。该理论指出在重建磁共振图像的问题中，如果图像在某个确定的变换域中能够稀疏的表示，那么利用核磁共振图像对应的频域下采样数据能够以很高的概率完整重建原始图像。一种重建图像的方法是将图像的重建问题转化为凸优化问题。已有的采用全变差约束和稀疏域l¹范数约束的凸优化重建方法，对纹理结构复杂的图像重建效果不佳，图像的细节很难重建。通过利用非局部相似性先验信息对图像的稀疏表示系数进行约束，可以降低重建图像中的块效应，改善图像中纹理结构的重建效果。进而降低所需的采样数据量，具有实际意义。>

发明内容

本发明的目的是提供一种充分利用图像非局部相似性的核磁共振图像重建方法，使得能够有效重建图像中的各种结构，提高重建质量。>

为了实现上述目的，本发明提供的技术方案如下：>

首先，对磁共振图像对应的傅里叶变换系数，采用变密度随机下采样的方式采样，由采样得到的数据通过傅里叶反变换得到用于重建的初始参考图像；其次，对初始参考图像进行分块，依据图像子块中的边缘和结构对图像子块进行分类，并得到每类图像子块对应的字典和图像子块的稀疏表示系数；最后，利用非局部相似性先验信息引入针对图像稀疏系数的约束项，结合图像在小波域的稀疏性，通过混合正则项求解模型对图像进行迭代重建。>

具体实现方案包括以下步骤：>

步骤(1)获取用于重建的初始参考图像，具体是：>

对磁共振图像对应的傅里叶变换系数，采用变密度随机下采样的方式采样，即对傅里叶系数对应图像低频信息的部分较多的采样，对傅里叶系数对应图像高频信息的部分较少采样；对得到的采样数据矩阵缺失的部分补零值，然后用二维傅里叶反变换得到用于重建的初始参考图像x⁽⁰⁾。>

步骤(2)对参考图像分块并对图像子块分类，求得每类对应的字典和图像子块对应类字典的稀疏系数，具体是：>

将输入图像x＝x⁽⁰⁾进行分块，即利用提取图像子块的矩阵P_i通过式x_i＝P_ix从x中得到大小为n×n的图像子块x_i，总的分块数为Π，其中对图像子块进行高通滤波，然后使用K均值算法对图像子块进行分类，得到具有相似结构特征的K类图像子块；利用主成分分析法(PCA)求得每一类的字典，由K类图像子块对应的K个子字典构成整个图像的字典；对于每一图像子块x_i选择该块所属类k对应的子字典则图像子块x_i对应的稀疏系数可以通过求解l¹范数最小问题得到，如式(1)所示：>

其中常数λ表示正则化参数，由此，得到对应的稀疏表示为

步骤(3)对图像子块稀疏系数进行约束，迭代重建磁共振图像，具体是：>

对于每个图像子块x_i，以x_i的中心为中心、大小为S×S的范围内搜索其非局部相似的块，得到x_i所有非局部相似图像子块的集合表示为C_i，Q(C_i)为该集合中的元素个数；利用图像的非局部相似性对重建图像进行估计，令x_i,j为x_i的非局部相似的图像子块，其稀疏系数表示为则重建图像子块x_i对应的稀疏系数的无偏估计表示为可通过式(2)计算：>

其中，ω_i,j是权重，可以通过式(3)计算，其中h为权重控制常数；>

${\omega }_{i,j}=\exp (-{\left|\right|x_i-x_{i,j}\left|\right|}_2^2/h)/W,W=\underset{j=1}{\overset{Q\left(C_i\right)}{\Sigma }}\exp (-{\left|\right|x_i-x_{i,j}\left|\right|}_2^2/h)---\left(3\right)$

非局部中心化的稀疏表示约束项为其中表示的转置；利用该约束项可使重建图像的稀疏系数逼近由估计得到的稀疏系数，从而得到重建模型为式(4)：>

$\hat{x}=\arg \underset{x}{\min }\{\frac{1}{2}{\left|\right|\mathrm{Rx}-b\left|\right|}_2^2+\beta {\left|\right|\tilde{\Phi }x\left|\right|}_1+\gamma \underset{i=1}{\overset{\Pi }{\Sigma }}{\left|\right|{\Phi }_{k_i}^T{x}_i-E\left[{\Phi }_{k_i}^T{x}_{i,j}\right]\left|\right|}_1\}---\left(4\right)$

其中，R表示局部傅里叶变换(Partial Fourier Transform)矩阵，K空间下采 样数据为b＝Rx+noise，noise为高斯白噪声，正则项中的稀疏变换矩阵采用小波变换矩阵，常数β、γ为正则化参数，为图像子块x_i所属类的子字典。>

模型求解之前给出如下定义：给定连续凸函数g(s)和参数ρ＞0，表示函数f在点u处的梯度，则s到u的近似映射过程定义为：>

${\mathrm{prox}}_{\rho }\left(g\right)\left(u\right):=\arg \underset{\text{s}}{\min }\{g\left(s\right)+\frac{1}{2\rho }{\left|\right|s-(u-\rho ▿f\left(u\right))\left|\right|}^2\}---\left(5\right)$

模型求解的具体步骤如下：>

1)输入：步骤(1)中重建的初始参考图像x⁽⁰⁾，迭代次数count＝1，最大迭代次数MaxIter，重建误差ε，初始化参数ρ，β，γ，t^(count)＝1，r^(count)＝x⁽⁰⁾；>

2)更新重建图像：$x_g=r^{\left(\mathrm{count}\right)}-\rho ▿f\left(r^{\left(\mathrm{count}\right)}\right),$其中$f\left(r^{\left(\mathrm{count}\right)}\right)=\frac{1}{2}{\left|\right|{\mathrm{Rr}}^{\left(\mathrm{count}\right)}-b\left|\right|}^2,$$▿f\left(r^{\left(\mathrm{count}\right)}\right)=R^T({\mathrm{Rr}}^{\left(\mathrm{count}\right)}-b);$

3)通过迭代阈值收缩算法求解，即分别求解：>

$x_{v_1}={\mathrm{prox}}_{\rho }\left(2\beta {\left|\right|\tilde{\Phi }x\left|\right|}_1\right)\left(x_g\right)---\left(6\right)$

x_gi＝P_ix_g

$x_{{\mathrm{iv}}_2}={\mathrm{prox}}_{\rho }\left(2\gamma \underset{i=1}{\overset{\Pi }{\Sigma }}{\left|\right|{\Phi }_{k_i}^T{x}_i-E\left[{\Phi }_{k_i}^T{x}_{i,j}\right]\left|\right|}_1\right)\left(x_{\mathrm{gi}}\right)---\left(7\right)$

$x_{v_2}={\left(\underset{i=1}{\overset{\Pi }{\Sigma }}{P}_i^T{P}_i\right)}^{-1}\underset{i=1}{\overset{\Pi }{\Sigma }}\left(P_i^T{x}_{{\mathrm{iv}}_2}\right)$

4)求解x_v1、x_v2的算术平均：

5)求解x^(count)在图像数据范围的投影，即求解x^(count)＝project(x^(count),[l,h])：投影函数θ＝project(θ,[l,h])定义为：l，h为常数，变量θ满足①当l≤θ≤h时，θ＝θ；②当θ＜l时，θ＝l；③当θ＞h时，θ＝h；>

6)更新迭代过程中的参数t^(count)，r^(count)：>

count＝count+1 (8)>

$t^{\left(\mathrm{count}\right)}=(1+\sqrt{1+4{\left(t^{(\mathrm{count}-1)}\right)}^2})/2---\left(9\right)$

r^(count)＝x^(count-1)+((t^(count-1)-1)/t^(count))(x^(count-1)-x^(count-2)) (10)>

7)判断迭代终止条件：当满足count＞MaxIter或者满足式(11)时停止迭代，执行步骤8)，否则返回步骤2)继续迭代；>

${\left|\right|r^{\left(\mathrm{count}\right)}-r^{(\mathrm{count}-1)}\left|\right|}_2^2/{\left|\right|r^{(\mathrm{count}-1)}\left|\right|}_2^2\le ϵ---\left(11\right)$

8)输出重构图像x＝x^(count-1)。>

本发明在算法的迭代过程中，利用图像中存在的非局部相似性得到估计的图像，在采用传统重建算法中基于稀疏性先验约束的同时，约束重建图像的稀疏系数逼近估计获得的图像的稀疏系数；得到了非局部稀疏表示正则化的磁共振图像重建模型，该模型能够有效降低重建图像中的块效应，对纹理结构复杂的图像能够较好的重建。>

具体实施方式

以下结合实施实例对本发明进行详细说明：>

步骤(1)获取用于重建的初始参考图像，具体是：>

对大小为256×256的核磁共振灰度图像进行傅里叶变换，采用变密度随机下采样的方式对傅里叶变换系数进行采样，即对傅里叶系数对应图像低频信息的部分较多的采样，对傅里叶系数对应图像高频信息的部分较少采样；可以取采样得到的数据量占全部傅里叶变换数据的16%-30%，如取20%；对得到的采样数据矩阵缺失的部分补零值，然后用二维傅里叶反变换得到用于重建的初始参考图像x⁽⁰⁾；>

步骤(2)对参考图像分块并对图像子块分类，求得每类对应的字典和图像子块对应类字典的稀疏系数，具体是：>

将输入图像x＝x⁽⁰⁾进行分块，即利用提取图像子块的矩阵P_i通过式x_i＝P_ix从x中得到大小为n×n的图像子块x_i，总的分块数为Π，其中对图像子块进行高通滤波，然后使用K均值算法对图像子块进行分类，得到具有相似结构特征的K类图像子块，此处K取40；利用PCA求得每一类的字典，由K类图像子块对应的K个子字典构成整个图像的字典；对于每一图像子块x_i选择该块所属类k对应的子字典则图像子块x_i对应的稀疏系数可以通过求解l¹范数最小问题得到，如式(1)所示：>

其中常数λ表示正则化参数，由此，得到对应的稀疏表示为

步骤(3)对图像子块稀疏系数进行约束，迭代重建磁共振图像，具体是：>

对于每个图像子块x_i，以x_i的中心为中心、大小为13×13的范围内搜索其非>i所有非局部相似图像子块的集合表示为C_i，Q(C_i)为该集合中的元素个数；利用图像的非局部相似性对重建图像进行估计，令x_i,j为x_i的非局部相似的图像子块，其稀疏系数表示为则重建图像子块x_i对应的稀疏系数的无偏估计表示为可通过式(2)计算：>

其中，ω_i,j是权重，可以通过式(3)计算，其中h为权重控制常数取75；>

${\omega }_{i,j}=\exp (-{\left|\right|x_i-x_{i,j}\left|\right|}_2^2/h)/W,W=\underset{j=1}{\overset{Q\left(C_i\right)}{\Sigma }}\exp (-{\left|\right|x_i-x_{i,j}\left|\right|}_2^2/h)---\left(3\right)$

非局部中心化的稀疏表示约束项为其中表示的转置；利用该约束项可使重建图像的稀疏系数逼近由估计得到的稀疏系数，从而得到重建模型为式(4)：>

$\hat{x}=\arg \underset{x}{\min }\{\frac{1}{2}{\left|\right|\mathrm{Rx}-b\left|\right|}_2^2+\beta {\left|\right|\tilde{\Phi }x\left|\right|}_1+\gamma \underset{i=1}{\overset{\Pi }{\Sigma }}{\left|\right|{\Phi }_{k_i}^T{x}_i-E\left[{\Phi }_{k_i}^T{x}_{i,j}\right]\left|\right|}_1\}---\left(4\right)$

其中，R表示Partial Fourier transform矩阵，K空间下采样数据为b＝Rx+noise，noise为高斯白噪声，正则项中的稀疏变换矩阵采用小波变换矩阵，常数β、γ为正则化参数，为图像子块x_i所属类的子字典；>

模型求解之前给出如下定义：给定连续凸函数g(s)和参数ρ＞0，表示函数f在点u处的梯度，则s到u的近似映射过程定义为：>

${\mathrm{prox}}_{\rho }\left(g\right)\left(u\right):=\arg \underset{\text{s}}{\min }\{g\left(s\right)+\frac{1}{2\rho }{\left|\right|s-(u-\rho ▿f\left(u\right))\left|\right|}^2\}---\left(5\right)$

模型求解的具体步骤如下：>

1)输入：步骤(1)中重建的初始参考图像x⁽⁰⁾，迭代次数count＝1，最大迭代次数MaxIter，重建误差ε，初始化参数ρ，β，γ，t^(count)＝1，r^(count)＝x⁽⁰⁾；>

2)更新重建图像：$x_g=r^{\left(\mathrm{count}\right)}-\rho ▿f\left(r^{\left(\mathrm{count}\right)}\right),$其中$f\left(r^{\left(\mathrm{count}\right)}\right)=\frac{1}{2}{\left|\right|{\mathrm{Rr}}^{\left(\mathrm{count}\right)}-b\left|\right|}^2,$$▿f\left(r^{\left(\mathrm{count}\right)}\right)=R^T({\mathrm{Rr}}^{\left(\mathrm{count}\right)}-b);$

3)通过迭代阈值收缩算法求解，即分别求解：>

$x_{v_1}={\mathrm{prox}}_{\rho }\left(2\beta {\left|\right|\tilde{\Phi }x\left|\right|}_1\right)\left(x_g\right)---\left(6\right)$

x_gi＝P_ix_g

$x_{{\mathrm{iv}}_2}={\mathrm{prox}}_{\rho }\left(2\gamma \underset{i=1}{\overset{\Pi }{\Sigma }}{\left|\right|{\Phi }_{k_i}^T{x}_i-E\left[{\Phi }_{k_i}^T{x}_{i,j}\right]\left|\right|}_1\right)\left(x_{\mathrm{gi}}\right)---\left(7\right)$

$x_{v_2}={\left(\underset{i=1}{\overset{\Pi }{\Sigma }}{P}_i^T{P}_i\right)}^{-1}\underset{i=1}{\overset{\Pi }{\Sigma }}\left(P_i^T{x}_{{\mathrm{iv}}_2}\right)$

4)求解x_v1、x_v2的算术平均：

5)求解x^(count)在图像数据范围的投影，即求解x^(count)＝project(x^(count),[l,h])：投影函数θ＝project(θ,[l,h])定义为：l，h为常数，变量θ满足①当l≤θ≤h时，θ＝θ；②当θ＜l时，θ＝l；③当θ＞h时，θ＝h；>

6)更新迭代过程中的参数t^(count)，r^(count)：>

count＝count+1(8)>

$t^{\left(\mathrm{count}\right)}=(1+\sqrt{1+4{\left(t^{(\mathrm{count}-1)}\right)}^2})/2---\left(9\right)$

r^(count)＝x^(count-1)+((t^(count-1)-1)/t^(count))(x^(count-1)-x^(count-2))(10)>

7)判断迭代终止条件：当满足count＞MaxIter或者满足式(11)时停止迭代，执行步骤8)，否则返回步骤2)继续迭代；>

${\left|\right|r^{\left(\mathrm{count}\right)}-r^{(\mathrm{count}-1)}\left|\right|}_2^2/{\left|\right|r^{(\mathrm{count}-1)}\left|\right|}_2^2\le ϵ---\left(11\right)$

8)输出重构图像x＝x^(count-1)。>