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

摘要

本发明涉及一种基于稀疏表示和非局部相似的核磁共振图像重建方法，主要目的是提高核磁共振图像的重建质量。本发明的具体步骤是：首先，对磁共振图像对应的傅里叶变换系数，采用变密度随机下采样的方式采样，对采样数据进行傅里叶反变换得到用于重建的初始参考图像；其次，对参考图像进行分块，得到的每类图像子块具有相似的结构特征，并得到每类图像子块对应的字典和图像子块的稀疏表示系数；最后，利用图像子块的非局部相似性估计原始图像，对图像子块的稀疏系数进行约束，结合图像在小波域的稀疏性，通过混合正则项求解模型进行迭代重建。本发明充分利用图像的非局部相似性，能有效重建图像中的复杂纹理，提高重建图像质量。

主权项

一种基于稀疏表示和非局部相似的核磁共振图像重建方法，其特征在于该方法包括如下步骤：步骤(1)获取用于重建的初始参考图像，具体是：对磁共振图像对应的傅里叶变换系数，采用变密度随机下采样的方式采样，即对傅里叶系数对应图像低频信息的部分较多的采样，对傅里叶系数对应图像高频信息的部分较少采样；对得到的采样数据矩阵缺失的部分补零值，然后用二维傅里叶反变换得到用于重建的初始参考图像x(0)；步骤(2)对参考图像分块并对图像子块分类，求得每类对应的字典和图像子块对应类字典的稀疏系数，具体是：将输入图像x＝x(0)进行分块，即利用提取图像子块的矩阵Pi通过式xi＝Pix从x中得到大小为n×n的图像子块xi，总的分块数为Π，其中对图像子块进行高通滤波，然后使用K均值算法对图像子块进行分类，得到具有相似结构特征的K类图像子块；利用主成分分析法求得每一类的字典，由K类图像子块对应的K个子字典构成整个图像的字典；对于每一图像子块xi选择该块所属类k对应的子字典则图像子块xi对应的稀疏系数通过求解l1范数最小问题得到，如式(1)所示：其中常数λ表示正则化参数，由此，得到对应的稀疏表示为步骤(3)对图像子块稀疏系数进行约束，迭代重建磁共振图像，具体是：对于每个图像子块xi，以xi的中心为中心、大小为S×S的范围内搜索其非局部相似的块，得到xi所有非局部相似图像子块的集合表示为Ci，Q(Ci)为该集合中的元素个数；利用图像的非局部相似性对重建图像进行估计，令xi,j为xi的非局部相似的图像子块，其稀疏系数表示为则重建图像子块xi对应的稀疏系数的无偏估计表示为通过式(2)计算：其中，ωi,j是权重，通过式(3)计算，其中h为权重控制常数；ωi,j=exp(-||xi-xi,j||22/h)/W,W=Σj=1Q(Ci)exp(-||xi-xi,j||22/h)---(3)]]>非局部中心化的稀疏表示约束项为其中表示的转置；利用该约束项使重建图像的稀疏系数逼近由估计得到的稀疏系数，从而得到重建模型为式(4)：x^=argminx{12||Rx-b||22+β||Φ~x||1+γΣi=1Π||ΦkiTxi-E[ΦkiTxi,j]||1}---(4)]]>其中，R表示局部傅里叶变换矩阵，K空间下采样数据为b＝Rx+noise，noise为高斯白噪声，正则项中的稀疏变换矩阵采用小波变换矩阵，常数β、γ为正则化参数，为图像子块xi所属类的子字典；模型求解之前给出如下定义：给定连续凸函数g(s)和参数ρ＞0，表示函数f在点u处的梯度，则s到u的近似映射过程定义为：proxρ(g)(u):=argmins{g(s)+12ρ||s-(u-ρ▿f(u))||2}---(5)]]>模型求解的具体步骤如下：1)输入：步骤(1)中重建的初始参考图像x(0)，迭代次数count＝1，最大迭代次数MaxIter，重建误差ε，初始化参数ρ，β，γ，t(count)＝1，r(count)＝x(0)；2)更新重建图像：其中3)通过迭代阈值收缩算法求解，即分别求解：xv1=proxρ(2β||Φ~x||1)(xg)---(6)]]>xgi＝Pixgxiv2=proxρ(2γΣi=1Π||ΦkiTxi-E[ΦkiTxi,j]||1)(xgi)---(7)]]>xv2=(Σi=1ΠPiTPi)-1Σi=1Π(PiTxiv2)]]>4)求解的算术平均：5)求解x(count)在图像数据范围的投影，即求解x(count)＝project(x(count),[pl,ph])：投影函数vp＝project(vp,[pl,ph])定义为：pl，ph为常数，变量vp满足①当pl≤vp≤ph时，vp＝vp；②当vp＜pl时，vp＝pl；③当vp＞ph时，vp＝ph；6)更新迭代过程中的参数t(count)，r(count)：count＝count+1                                       (8)t(count)=(1+1+4(t(count-1))2)/2---(9)]]>r(count)＝x(count‑1)+((t(count‑1)‑1)/t(count))(x(count‑1)‑x(count‑2))      (10)7)判断迭代终止条件：当满足count＞MaxIter或者满足式(11)时停止迭代，执行步骤8)，否则返回步骤2)继续迭代；||r(count)-r(count-1)||22/||r(count-1)||22≤ϵ---(11)]]>8)输出重构图像x＝x(count‑1)。