一种基于小波域结构和非局部分组稀疏的MR图像重建方法与流程

本发明属于基于压缩感知原理的图像重建方法，尤其涉及一种基于小波域结构和非局部分组稀疏的mr图像重建方法，结合图像在小波域的系数结构和图像自身的非局部低秩这两种特性。

背景技术：

核磁共振是一种核物理现象，描述的是在外磁场的作用下磁矩不为零的原子发生自旋能级分裂，共振吸收一定频率的射频能量的过程。最初,在1946年bloch、purcell等人发现了磁共振现象。人们利用这一发现用于波谱学来探索物质微观结构，从而诞生了核磁共振这一新兴学科。

磁共振成像(magneticresonanceimaging,mri)技术依靠磁体技术、低温超导技术、电子与计算机技术等相关技术的不断进步已经成为现代医学影像领域中的重要手段，在物理、化学、生物技术、临床诊疗等领域得到了长足发展。特别是技术依靠其成像清晰度高、成像角度多、无电离辐射等优点成为了现代医学影像技术中重要的一员。特别是在医学应用当中，磁共振成像对软组织的检测灵敏度高，组织的成像对比度高，不仅能反映人体的解剖结构信息，还能反映出组织中的生理化学变化过程，为临床医学提供了非常有价值的诊断依据，促进了医学诊断技术的发展。

但磁共振成像在实际的应用中还存在一些问题，如磁共振图像的成像时间由两方面组成，一方面是进行磁共振数据采集的扫描时间，另一方面是对采集到的磁共振数据进行重建的时间，磁共振的成像时间往往较长。尤其是扫描时间，在早期的核磁共振成像中，一次成像过程中的扫描时间往往就需要几分钟，这个过程中要求病人要完全静止，否则会造成运动伪影使成像的质量下降。一些儿童或不能自主的病人往往无法在这个过程中保持不动。因此，过长的扫描时间成为了技术发展中必须要解决的问题。解决这一问题主要有两个途径，一个是通过提高数据采集速度，另一个是对重建算法进行改进，达到利用较少的采样数据量尽可能精确的重建图像。

基于压缩传感理论的图像重建方法。近年来发展起来的压缩感知理论能够突破奈奎斯特数据采样率，很好的解决了数据的欠采样恢复问题。该理论利用信号能够在某个变换域中稀疏表示的先验知识，提出可以通过部分采样得到的数据以很高的概率完整恢复原始信号。2006年candes从理论上证明了利用局部傅里叶变换系数能够精确重构原始信号，提出了压缩感知理论。能够利用远小于传统采样方法所需要的采样数据量重构原始信号，突破了香农采样定理的瓶颈^[1-2]。

mr图像信号可以通过对其进行稀疏变换来满足压缩感知重建的要求，对于有些在图像域就具有稀疏性的图像，可以不经过稀疏变换直接进行采样重建，如血管磁共振成像。因此cs理论在mri重建领域有很好的应用优势。

lusting^[3]等人在2007年将压缩感知理论应用于重建磁共振图像的问题中，基于图像在某个确定的变换域能够稀疏表示的前提下，利用k空间下采样数据重构的图像在图像中出现的混迭伪影是不连贯的，因此能够通过非线性重建的方法以很高的概率完整重建原始图像，提出了经典的重建模型如式(1)所示。通过实验验证了该方法的可行性和有效性，结果表明该法能大幅度减少mr图像重建所需的k空间数据量^[3]。

这样将图像的重建转化为求解一个约束优化问题。在重建模型中，观测到的数据表示为y＝fux+n，其中fu表示局部傅里叶变换，x表示原始图像，n表示高斯白噪声。表示重建的磁共振图像，ψ表示稀疏变换，||x||tv表示图像的全变分约束，α和β是两个非负参数。国内学者屈小波等通过将曲波、方向性小波等变换引入mri重建的过程中，得到了较好的重建效果。在图像通过小波稀疏表示方面取得了不错的成果，提出了基于图像块的方向小波mri重建算法，将图像子块进行训练得到自适应的方向波稀疏表示，通过方向性小波对压缩感知重建模型进行稀疏约束重建^[4]。

由于全变分正则项基于图像是分块光滑的先验，模型对于图像中丰富的纹理重建效果不佳。近几年图像的非局部相似性在图像重建领域起到了很大的作用，非局部分组稀疏性对于自然信号的建模和纹理恢复部分尤为重要。非局部相似块不仅在变换域是稀疏的，而且它们的非零元素共享一个联合稀疏模式。通过建模的稀疏系数之间的相关性，能大幅度减少导致更准确的重建^[5-6]未知信号的不确定性。文献[7]针对结构化稀疏提出了非局部低秩正规化(nlr)的方法，并应用到了压缩感知mri的重建上面。并且取得了很好的重建效果。

因此利用部分k空间数据重建原始图像的过程，本质上是利用了图像存在的一些先验知识，通过先验条件约束求解方程达到精确重建。因此，利用哪些先验条件，如何利用这些先验约束直接关系到重建图像的质量，研究重建模型及先验约束成为了重建算法研究的重点。

参考文献

[1]donohodl.compressedsensing[j].ieeetransactionsoninformationtheory,2006,52(4):1289-1306

[2]candesej,rombergj,taot.robustuncertaintyprinciples:exactsignalreconstructionfromhighlyincompletefrequencyinformation[j].ieeetransactionsoninformationtheory,2006.52(2):489-509

[3]lustingm,donohod,paulyjm.sparsemri:theapplicationofcompressedsensingforrapidmrimaging[j].magnresonmed,2007,58(6):1182-1195.[4]quxb,guod,ningbd.etal.undersampledmrireconstructionwithpatch-baseddirectionalwavelets[j].magneticresonanceimaging,2012,30(7):964-977

[5]a.buades,b.coll,andj.m.morel,“areviewofimagedenoisingalgorithms,withanewone,”multiscalemodel.simul.,vol.4,no.2,pp.490–530,2005.

[6]w.dong,g.shi,x.li,l.zhang,andx.wu,“imagereconstructionwithlocallyadaptivesparsityandnonlocalrobustregularization,”signalprocess.,imagecommun.,vol.27,no.10,pp.1109–1122,2012.

[7]weishengdong,guangmingshi,xinli,yima,fenghuang:compressivesensingvianonlocallow-rankregularization.ieeetrans.imageprocessing23(8):3618-3632(2014)。

技术实现要素：

本发明要解决的技术问题是，提出一种基于小波域结构和非局部分组稀疏的mr图像重建，利用了图像在小波域下系数的结构稀疏特性有利于捕捉原始信号隐含的结构和细节信息，特别是，当与信号的结构特性相匹配时，可以较大程度地改善图像信号的重建质量。同时图像信号中存在着很强的时空相关性，具有非局部相似性，能够很好的保留图像细节。通过将结构和分组稀疏统一到同一个框架中，最终提高了图像的重建质量。本发明主要应用于医学核磁共振图像的重建，这种方法在主客观上都有一定的提高。

为实现上述目的，本发明采用如下的技术方案：

1、一种基于小波域结构和非局部分组稀疏的mr图像重建方法，其特征在于，具体包括以下步骤：

步骤1、首先将图像变换到傅里叶域进行随机采样，得到采样数据y，之后，对其利用基本的压缩感知重建图像原理进行初始化；

步骤2、迭代奇异值阈值法求解低秩矩阵li

对初始化的整副图像提取一系列重叠块，这个过程用表示，对提取出来的每一块，在局部窗口内用k近邻方法寻找与其相似的前m块，则其中，然后再通过奇异值阈值法求解出低秩矩阵li；

步骤3、交替方向乘子法(admm)求解图像x

对式(10)引入两个辅助变量和z∈rⁿ，将(10)式重写成下式：

对(11)式子运用admm算法，其增广拉格朗日函数为：

其中，μ∈rⁿ和γ∈rⁿ是拉格朗日乘子，β1,β2＞0是约束x＝z和的惩罚参数，(12)式可以重写为：

对于(13)式的优化包含下面几部分迭代：

运用admm算法将其分解成3个子问题进行求解图像x，具体过程如下：

1、对于求解z^(l+1)，式(14)本身含有闭合解，可以通过最小二乘法一步求解出来；

2、对于求解x^(l+1)，式(16)同样含有闭合解，可以通过最小二乘法一步求解出来；

3、对于求解式(15)中的小波系数结构为高频和低频(wl2-l1)约束的时候，将(28)式子放在小波域讨论求解，即

当式(15)中的小波系数结构为父子分组(group)约束的时候，同理将(32)式子放在小波域讨论求解，即，

作为优选，对于求解的具体过程如下：

a)、当式(15)可以写成：

在求解的时候将(28)式子放在小波域讨论求解，令和分别表示小波变换后的所有系数、高频系数和低频系数，即，(28)式的优化解为：

根据同理b^(l)也可以分为和所以(29)可以分为2个问题分别求解：

式(30)和(31)可以通过proximal方法很容易解出来，最后，可以通过小波逆变换得到，即，

b)、同理当时，式(15)可以写成：

同a)的处理过程一样，将式(32)放在小波域讨论，则为：

将待求向量b和扩展成和则(33)式可以写成：

其中，扩展成这样存在矩阵g和h，使得和同理和式子(34)可以通过spams算法求得，同理，可以通过小波逆变换得到，即，

附图说明

图1.基于小波域结构稀疏和非局部分组稀疏的mr图像的压缩感知框架；

图2(a)为mri测试的brain图像；

图2(b)为mri测试的chest图像；

图2(c)为mri测试的heart图像；

图2(d)为mri测试的artery图像；

图3.观测矩阵；

图4.本发明的整体流程图。

具体实施方式

压缩感知理论主要研究的是如何通过获得较少的采样观测值，利用信号本身具有的稀疏性先验知识进行信号的重建。不同于传统的信号获取方式，压缩感知理论通过稀疏表示、投影测量和重建算法来获取信号。使得采样率能够远低于奈奎斯特的采样频率，大大降低了信号的获取和传输成本，通过重建算法能够以很高的概率精确重建原始信号。

基于压缩感知理论重建信号本质上是利用了图像可以被稀疏表示的特性和图像本身存在的多种先验知识来对重建模型进行约束求解，这些先验信息以正则项的形式包含进信号的重建模型中。图像信号通常是有结构的，一般在某基底下可以稀疏表示。同时本发明结合图像在小波域的系数结构稀疏和图像自身的非局部相似分组稀疏这两种先验信息，通过交替优化达到最优重建效果，具有以下特征：

1.模型框架

压缩感知理论表明：在某个变换域下稀疏的信号，可以利用优化方法由少量观测数据稀疏重建。这种非自适应的压缩采样将信号中包含的信息凝聚在少量的观测数据上，大幅度地降低了精确重建原始信号所需要的采样数目,然后再利用图像本身存在的多种先验知识来对信号进行重建。本发明提出了在小波变换域下，利用图像在小波域的系数结构特性和图像自身的非局部低秩特性这两种先验信息来重建原图。

1)图像非局部分组稀疏

图像中像素之间或多或少地存在一定的联系，称之为图像的相似性。以像素点为中心的窗口邻域(或图像块)为例，除了中心像素很可能和其周围像素相似以外，该邻域还可能和图像中其他部位的某些邻域像素相似(结构特征相似)，使得处于图像中不同位置处的像素点往往表现出很强的相关性，常称之为图像的非局部相似性。

通过图像块匹配方法在整个图像(为了降低计算量通常是在一个比较大的搜索窗口里进行)内搜索与当前图像块相似的图像块，然后对找到的相似图像块进行加权，从而能很好的保持图像的纹理细节。

2)小波域结构稀疏

小波变换的多分辨特性导致了图像小波域呈现如下特性：(1)小波系数的能量集中在最低频部分：最低子带承载着空域统计相关区域的信息,因此包含了图像的主要能量；(2)高频系数是稀疏的：各个级别的高频子带主要承载了图像中边和物体的边界信息，并且该部分所承载的视觉信息对于整幅图像贡献远远超出其能量对整幅图像的贡献，因此刻画了图像的结构。(3)零树结构：每一个系数节点与四个高一级子带系数节点相关，即父节点与子节点通常同时为非零(或零)。

根据小波域系数的特性，分为低频+高频和低频+高频内父子两种分组情况。本发明利用这2种分组模式，将这2种分组特性作为传统稀疏表示的另一个约束，不仅能进一步提升不同稀疏基的稀疏表示能力，细化不同图像的稀疏约束，同时又在整副图像上做了整体的稀疏约束。

通过将分组稀疏和小波域结构稀疏统一到同一个框架中，提出了基于结构和分组稀疏的mri图像压缩感知重建，如图1所示。

本发明模型的目标函数如如下：

其中，y频域采样数据，φ是观测矩阵，x是需要重建的图像，令xi＝rix，xi表示从图像i位置提取出来的大小为的图像块，xi∈cⁿ，其中ri是一个矩阵，作用是从图像x中的i位置提取图像块xi，令矩阵中的每一列都是xi的相似块，所以具有低秩性质，其中di是过完备字典，是一个稀疏系数矩阵，ω(·)是小波域结构稀疏约束。η，α，θ是正则项系数。

由于具有低秩性质，因为实际中易受噪声污染，脱离低秩约束，所以通过低秩矩阵li∈c^n×m,n≤m逼近低秩矩阵li通过优化下式求解：

代替凸核准则，使用非凸对数logdet(x)作为秩的平滑代理函数。对于一个半正定矩阵x∈r^n×n，秩最小优化问题可以近似为矩阵奇异值的log之和:

对于非半正定的低秩矩阵li∈c^n×m，秩最小优化问题可以通过对其进行奇异值分解来求解：

其中ε是一个很小的常数。σ是一个对角矩阵，矩阵中的对角元素是矩阵lili^t的特征值。σ^1/2也是一个对角矩阵，矩阵中的对角元素是矩阵li的奇异值。

通过选择合适的λ值，(3)式可以重写为：

总的目标函数(2)可以重写为：

对于ω(·)是小波域结构稀疏约束，本发明利用了两种小波域系数结构，在研究内容部分有介绍。第一种是高频+低频结构稀疏，ω(·)＝ω||plψx||2+(1-ω)||phψx||1，其中pl和ph分别表示提取图像x在小波变换域后的低频和高频系数，ω是高低频系数的权重系数。第二种是低频+高频内父子分组结构稀疏，g是变量的不同分组形式的集合，向量表示集合g中某一个分组gi中的系数。ωi表示不同分组的权值。

2.模型的优化

由于模型(6)中含一项保真项和两个正则约束项，而且约束项都是不光滑的，求解起来很困难。本发明提出了一种优化方法，协同优化保真项，小波域稀疏约束项和非局部分组稀疏约束这三项。将这个优化问题分解成了两个子问题，通过交替迭代奇异值阈值法和交替方向乘子法来求解。

1)奇异值阈值法求解低秩矩阵li

根据^[7]的文章，(7)式可以写成下式：

其中，τ＝λ/2η，k表示迭代次数，n0＝min{n,m}，表示带有权重的加权核范数，σj是li的第j个奇异值，可以通过下式求解：

其中，是xi的奇异值分解，(x)+＝max(x,0)，我们初始化w＝[1,1,...,1]^t，这样第一次迭代就等价于求解一个普通核范数的最小化问题。

2)交替方向乘子法求解图像x

交替方向乘子法(admm)是一种非常流行的用于处理大规模优化问题的方法。因此我们首先引入两个辅助变量和z∈rⁿ。将(10)式重写成下式：

对(11)式子运用admm算法，其增广拉格朗日函数为：

其中，μ∈rⁿ和γ∈rⁿ是拉格朗日乘子，β1,β2＞0是约束x＝z和的惩罚参数，(12)式可以重写为：

对于(13)式的优化包含下面几部分迭代：

其中ρ＞1是一个常数，l表示迭代次数，对于这几个变量的求解，可以分成3个子问题来解。

a.求解z^(l+1)

对于式(14)可以用一步就解出来，其闭合形式解为：

其中，是一个对角矩阵，表示表示相似块的平均值。

b.求解x^(l+1)

同理对于式子(16)，其解为：

其中，φ是部分傅里叶变换矩阵，φ＝df，d和f分别表示下采样矩阵和傅里叶变换矩阵。将φ＝df带入(22)得：

解得：

c.求解

对于式(15)，用proximal方法求解，令h(x)作为一个凸函数，h(x)的proximal算子定义为：

这个函数对每个u都有唯一最小值。对于l1范式h(x)＝||x||1的情况，收敛算子proxth为：

对于l2范式h(x)＝||x||2，proxth为：

在本发明中有两种小波域小波域系数结构，分别是

和

a)当式(15)可以写成：

在求解的时候，因为低频和高频系数没有重叠部分，并且小波变换是可逆的，所以，我们将(28)式子放在小波域讨论求解。令和分别表示小波变换后的所有系数、高频系数和低频系数。即，(28)式的优化解为：

根据同理b^(l)也可以分为和所以(29)可以分为2个问题分别求解：

式(30)和(31)可以通过proximal方法很容易解出来。最后，可以通过小波逆变换得到，即，

b)同理当时，式(15)可以写成：

同a)的处理过程一样，将式(32)放在小波域讨论，则为：

由于这个分组稀疏约束项中有交叠分组，可能导致每次迭代求解过程中某一待求系数的值因为处在不同的分组，计算的结果会有差异，因此将待求向量b和扩展成和则(33)式可以写成：

其中，扩展成这样就存在矩阵g和h，使得和同理和式子(34)是一个常见的优化问题，可以通过spams算法求得。同理，可以通过小波逆变换得到，即，

如图4所示，本发明的基于小波域结构和非局部分组稀疏的mr图像重建方法，具体包括以下步骤：

步骤1、首先将图像变换到傅里叶域进行随机采样，得到采样数据y。之后，对其利用基本的压缩感知重建图像的方法(本发明采用dct方法)进行初始化。

步骤2、迭代奇异值阈值法求解低秩矩阵li

对初始化的整副图像提取大小为6×6的一系列重叠块，这个过程用表示，对提取出来的每一块，在局部窗口内用k近邻方法寻找与其相似的前43块，则其中，m＝43。然后再通过奇异值阈值法求解出低秩矩阵li。

步骤3、交替方向乘子法(admm)求解图像x

对式(10)引入两个辅助变量和z∈rⁿ，运用admm算法将其分解成3个子问题进行求解；其中，交替方向乘子法(admm)求解图像x，具体过程如下：

1、对于求解z^(l+1)，式(14)本身含有闭合解，可以通过最小二乘法一步求解出来。

2、对于求解x^(l+1)，式(16)同样含有闭合解，可以通过最小二乘法一步求解出来。

3、对于求解式(15)中的小波系数结构为高频和低频(wl2-l1)约束的时候，因为低频和高频系数没有重叠部分，并且小波变换是可逆的，所以，我们将(28)式子放在小波域讨论求解。和分别表示小波变换后的所有系数、高频系数和低频系数，其中之后将高频l1范式和低频l2范式分别用对应的proximal方法求解出来。最后，可以通过小波逆变换得到，即当式(15)中的小波系数结构为父子分组(group)约束的时候，同理将(32)式子放在小波域讨论求解，由于这个分组稀疏约束项中有交叠分组，可能导致每次迭代求解过程中某一待求系数的值因为处在不同的分组，计算的结果会有差异，因此将待求向量b和扩展成和通过spams算法求得。同理，可以通过小波逆变换得到，即，

本发明的模型的实现，采用了真是的二维核磁共振(mr)图像。4副图像分别是210×210的brain图像，220×220的chest图像，192×192的heart图像和220×220的artery图像，分别如图2(a)，图2(b)，图2(c)，图2(d)所示。

在图像的观测设置上，采用了部分傅里叶观测方法。观测矩阵如图3所示。

假设mr图像x具有n个像素，部分傅里叶变换矩阵r是由n×n的傅里叶变换矩阵中随机选出的m行组成。随机选出的m行对应的观测组成观测值向量y。采样率定义为m/n。如果采样率较小，则mr图像扫描时间较短。

为了验证上述所提到的基于压缩感知的核磁共振图像重建方案的有效性，我们对几幅常用的医学图像进行重建，比较了一下其客观质量，客观质量主要是通过峰值信噪比(peaksignaltonoiseratio，psnr)度量，单位为分贝(db)。其计算公式如下：

两幅大小为m*n的图像的均方误差mse的定义如下：

其中i，j分别表示原始图像和重建图像，而i(x,y),j(x,y)为对应于位置(x,y)处的像素值，则均方误差越小，则psnr越高，则重建图像的质量越高。

为了比较图像重建的质量，表1给出重建的灰度图像的psnr比较结果：

表1.mri图像在不同采样率下的重建结果

表1中，后两列是本发明提出的模型，其中nlr+wl1-l2是小波域系数约束选取高频+低频时候的模型，nlr+group是小波域系数约束选取低频+父子分组时的模型。nlr模型是基于非局部低秩约束的压缩感知重建模型。为了有效说明不同的模型对图像重建质量的作用，在这里具体比较实验分为两个部分：(1)基于非局部低秩约束的nlr模型重建结果和nlr+wl1-l2模型重建的比较；(2)基于非局部低秩约束的nlr模型重建结果和nlr+group模型重建的比较。

如表1所示，本发明模型的优势比较明显，较nlr方法有一定提升。虽然在各个采样率下本文nlr+group模型的优势比模型nlr+wl2-l1的优势小一些，但是也足以证明本文提出的模型要优于nlr模型。