网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

基于非局部低秩约束的改进灵敏度编码重建算法  PDF

  • 潘婷
  • 段继忠
昆明理工大学信息工程与自动化学院,昆明 650500

中图分类号: TP391

最近更新:2023-01-16

DOI:10.16337/j.1004⁃9037.2023.01.017

  • 全文
  • 图表
  • 参考文献
  • 作者
  • 出版信息
EN
目录contents

摘要

灵敏度编码(Sensitivity encoding, SENSE)是一种应用广泛的并行磁共振成像(Magnetic resonance imaging, MRI)重建模型。目前已有的针对SENSE模型的改进方法的重建图像中依然有较多伪影,尤其在较高加速因子时很难重建出比较清晰的图像。因此,本文基于非局部低秩约束(Nonlocal low‑rank, NLR),提出了一种改进的SENSE模型,称为NLR‑SENSE。该模型使用加权核范数作为秩代理函数,并使用交替方向乘子法(Alternating direction multiplier method, ADMM)进行求解。仿真实验结果表明,与其他几种并行磁共振成像方法相比,NLR‑SENSE方法在视觉比较和3个不同的客观指标上均表现优异,能有效提升重建图像的质量。

引 言

磁共振成像是目前主要的无创医学成像技术之一,在学术研究与临床应用方面都有很大的研究价

1。磁共振的扫描时间较长,这在很大程度上影响了磁共振成像的应用。目前有两种较为常用的技术可以用于加快磁共振成像:一种是压缩感知(Compressed sensing, CS2‑3,该理论表明只需要较少的测量数据就能完成质量较高的数据重建;另一种是并行成像,该技术使用空间灵敏度编码方法同时进行数据采集并完成数据重建。这两种方法都已经广泛应用于医学磁共振成像。

并行成像方法使用一组多通道接收线圈同时采集数据,利用线圈的灵敏度差异完成空间编码,减少了磁共振成像过程中相位编码步骤。在并行成像中对K空间数据进行欠采样,能有效缩短数据采集的时间,从而加快磁共振成像速度。灵敏度编码(Sensitivity encoding, SENSE

4是一种十分有效的并行磁共振重建算法,它显式利用灵敏度信息进行并行磁共振成像的重建,在具有准确的灵敏度信息时有最佳的重建结果。SENSE模型的重建是一个超定问题(线圈数通常大于加速因子)的求解,对这样的病态问题直接求解时,重建的图像中会存在较多伪影,正则化技术能有效提升算法的重建性能。常见的正则项有基于图像稀疏5‑6的L12、基于图像局部特性的全变7‑8等。过去的研究中,已经提出了很多基于正则化的SENSE方法,包含了基于小波域稀疏约束的L1‑SENSE模型以及使用全变分(Total variation, TV)正则化约束的TV‑SENSE模9。文献[10]提出了基于Lp伪范数和全变分约束的LpTV‑SENSE模型,实验结果证明了它能进一步提高图像的重建质量。

事实上,用于正则化约束的图像先验信息不只有图像的稀疏性和局部特性,图像本身还具有非局部特性,如非局部相似性(图像的某些部分存在相似的特性)。在图像去噪以及压缩感知图像重建领域,已经出现了很多对于非局部相似特性的研究。非局部平均(NL‑means)方

11是非局部特性方法的开创者,该方法以图像块为单位,查找得到相似块,通过平均的方式可以有效去除图像中存在的噪声。之后提出的基于补丁的非局部算子(Patch‑based nonlocal operator, PANO12以及块匹配三维去13方法,都是利用了非局部相似特性的有效的去噪方法。在压缩感知磁共振成像方法,Dong14提出了将非局部特性转化为低秩约束正则化,得到一种基于非局部低秩(nonlocal low‑rank, NLR)的压缩感知磁共振成像重建算法,称为NLR‑CS,具有十分优越的图像重建性能。朱俊15在压缩感知成像中将非局部与局部的正则化方法进行结合,得到了更好的重建效果。

然而,现有的结合了SENSE的众多模型的重建图像中依然包含较多伪影,重建性能还有较大的提升空间。使用非局部特性作为先验信息的正则化方法在图像去噪以及压缩感知图像重建领域表现出了优越的性

12‑1416。本文基于SENSE模型,结合NLR,提出NLR‑SENSE模型,并使用交替方向乘子法(Alternating direction multiplier method, ADMM17‑18对所提出的目标优化问题进行有效的求解。此外,对于模型中的秩最小化问题,本文使用加权核范19作为秩的代理函数进行求解,与传统的核范20相比,可以得到更近似低秩的解。

本文将NLR‑SENSE方法与其他几种具有竞争的并行磁共振成像重建算法进行了比较,包括一些SENSE模型的改进方法,如基于L1范数正则化约束的L1‑SENSE方法、基于全变分的TV‑SENSE方

9、基于Lp伪范数的全变分正则化的LpTV‑SENSE方10、结合了TV正则项和局部k空间邻域低秩矩阵模型(Low‑rank modeling of local k‑space neighborhood, LORAKS21的SENSE‑LORAKS‑TV方22。也与在迭代自洽并行成像重建(Iterative self‑consistent parallel imaging reconstruction,SPIRiT)模23中加入同时双向低秩(Simultaneous two‑directional low‑rankness, STDLR)约束的STDLR‑SPIRiT方24进行了比较。在不同数据集上使用不同欠采样方式进行了仿真实验实验,在所有对比的并行磁共振成像方法中NLR‑SENSE重建结果有最佳的视觉重建效果,并且在3个客观图像指标上整体表现优秀。

1 算法原理

1.1 SENSE模型回顾

假设XCN表示待重建的图像,N=Nh×Nv表示单线圈图像的像素点数,NvNh分别是单线圈图像的水平方向和垂直方向的像素点数;Y=[Y1T,,YCT]TCMC表示多线圈欠采样k空间数据,其中YcCM表示第c个线圈的欠采样k空间数据,c=1,,CC表示线圈数,M表示一个线圈的欠采样数据点数,M≪N。在SENSE重建模型中,对角矩阵ScCN×N表示第c个线圈的灵敏度信息,令灵敏度算子S=[S1T,,SCT]TCNC×N,则并行成像模型可表示成

Y=PFSX (1)

式中:P=ICPRMC×NC表示多线圈欠采样算子,IC表示大小为C×C的单位矩阵,表示克罗内克积;F=IC(UhUv)CNC×NC表示傅里叶算子;UhCNh×NhUvCNv×Nv分别表示Nh点和Nv点离散二维傅里叶变换矩阵。过去的很多工作中,SENSE模型的求解常写成以下最小二乘问题

X=arg minXY-PFSX22 (2)

式中:2表示向量的L2范数。

过去已经提出了多种基于SENSE的算法:包含使用图像在投影空间的稀疏性作为约束的L1‑SENSE、使用基于图像局部特性作为先验信息TV‑SENSE

9、结合了Lp伪范数的全变分正则化约束的LpTV‑SENSE10、以及结合了TV正则项和LORAKS方法的SENSE‑LORAKS‑TV22,都能在一定程度上改进SENSE模型的重建性能。

1.2 NLR‑SENSE模型

为了得到一种更加有效的并行磁共振成像重建算法,本文在SENSE中使用NLR约束,得到NLR‑SENSE重建模型。该重建模型充分利用了图像的非局部自相似特性先验信息,得到如下最优化问题

X=arg minX12Y-PFSX22+μi=1NprankVi(X) (3)

式中:Y-PFSX22为数据保真项;i=1NprankVi(X)为正则项;rankVi(X)表示Vi(X)的秩;μ>0表示正则化参数,用于均衡正则项与数据保真项;Vi:XCNVi(X)Cn×m表示块匹配算子。块匹配的步骤如下:将X划分为Np个大小为n×n的重叠块,n表示图像块的像素点个数,i=1,,Np表示块索引;对于每一个参考块,在它的局部邻域窗口中(例如设置局部窗口大小为40×40)搜索出最相似(即搜索块与参考块之间的欧几里得距离最小)的m个块;将m个相似块按相似程度排序并向量化之后作为Vi(X)的列,过程如图1所示,图中上部分为NLR‑SENSE算法过程框图,下部分为包括块匹配步骤的低秩最小化去噪过程原理图。

图1  NLR‑SENSE算法原理图

Fig.1  Schematic diagram of NLR‑SENSE algorithm

对于问题(3)这样的大型非凸最优化问题,可以使用ADMM技

25进行求解。具体来说,引入辅助变量Di=Vi(X)Cn×mZ=SXCNC,以及对应的拉格朗日乘子diCn×mzCNC,问题(3)可以转化为以下步骤交替更新求解

Dik+1=argminDiγ12Di-Vi(Xk)+dikF2+μDiw,* (4)
Zk+1=arg minZY-PFZ22+γ2Z-(SXk+zk)22 (5)
Xk+1=arg minXγ2SX-(Zk+1-zk)22+γ1i=1NpVi(X)-(Dik+1-dik)F2 (6)
zk+1=zk+η2(SXk+1-Zk+1) (7)
dik+1=dik+η1Vi(Xk+1)-Dik+1 (8)

式中:w,*表示加权核范数,其定义为Dw,*=j=1min(m,n)wjσj,其中σj为矩阵D的第j个奇异值,wj表示对应于σj的权重;F表示Frobenius范数;{}表示一系列的值;Di中包含NpDidi中包含Npdiγ1>0γ2>0为惩罚项参数;η1>0η2>0为乘子更新的参数。为了方便叙述,将所有的Di按顺序水平拼接成矩阵D=[D1,D2,,DNp]Cn×mNp,对应的,将所有的di按顺序水平拼接成矩阵d=[d1,d2,,dNp]Cn×mNp

(1)关于Di的秩最小化问题求解:对于问题(4)的求解,常规的方法是使用核范

20作为代理函数进行求解。然而,过去的研究表明使用加权核范数作为秩代理函数求解秩最小化问题可以得到更加近似低秩的19。因此,本文使用加权核范数作为秩的代理函数求解。得到D的最小化问题如下

Dk+1=argminD12D-Vi(Xk)+dikF2+τDw,* (9)

式中τ=μ/γ1

问题(9)是加权核范数最小化问题,可以通过奇异值阈值算法求解:将Vi(Xk)+dik进行奇异值分解,得到ΨΣΓH=Vi(Xk)+dik,其中Σ为奇异值对角矩阵,其第j个对角元素Σjj=σjVi(Xk)+dikσjVi(Xk)+dik表示Vi(Xk)+dik的第j个奇异值,将所有的奇异值组成一个行向量σ=σ1Vi(Xk)+dik,,σmin(m,n)Vi(Xk)+dikΨΓ均为酉矩阵,其列分别是奇异值的左、右奇异向量。然后,可以得到问题(9)的解为

Dik+1=Ψsoft(Σ,τw)ΓH (10)

式中:soft(Σ,τw)=max(Σ-diag(τw),0)表示软阈值算子;w=w1,,wmin(m,n)表示Vi(Xk)的奇异值对应的权重,其第j个值wj计算式为

wj=1σj+ε (11)

式中ε=10-16用于保证分母不为0。

(2)根据式(6),可以得到Z的解为

Zk+1=FHPHP+γ2I-1PHY+γ2F(SXk+zk) (12)

式中:PHPCNC×NC是一个对角矩阵;ICNC×NC表示单位矩阵,问题(13)可以直接求解。

(3)根据式(7),可以得到X的解为

Xk+1=γ2SHS+γ1i=1NpVi*Vi-1γ2SH(Zk+1-zk)+γ1i=1NpVi*Dik+1-dik (13)

式中:Vi*:Cn×mCNVi的伴随算子,表示将相似块组放回图像中对应位置;i=1NpVi*ViCN×N为对角矩阵,它的每个对角元素等于对应的像素点出现在所有{Vi(Xk)}中的次数。由S的定义易知SHSCN×N为对角矩阵。因此,γ2SHS+γ1i=1NpVi*Vi也为对角矩阵,其逆易得,从而问题(14)的求解也十分容易。

通过对问题(4)~(9)交替求解,就能得到最后的重建图像。算法1给出了NLR‑SENSE方法的求解过程。由于块匹配步骤十分复杂,为了减少算法的计算量,NLR‑SENSE算法中每10次迭代更新一次块匹配算子,即算法1中设置T=10

算法1: NLR‑SENSE求解过程

(1) Input: PFSY

(2) Initial:X0=(PFS)HYZ0=0z0=0Di0=0di0=0k=0

(3) Repeat

(4) If modk,T=0, Update the patch grouping Vi by block matching

(5) Construct similar patch group matrix Vi(Xk)

(6) Compute the full SVD ΓΣΨH of Vi(Xk)

(7) Update the weights w via Eq.(12)

(8) Compute Dk+1 via Eq.(11)

(9) Compute Xk+1 via Eq.(14)

(10) Compute Zk+1 via Eq.(13)

(11) Update zk+1 via Eq.(7)

(12) Update dik+1 via Eq.(8)

(13) k=k+1

(14) Until k<K

(15) Output the reconstructed image X=Xk+1

2 实验结果

2.1 实验设置

为了验证本文所提出的NLR‑SENSE方法的性能,将本文算法与几种具有竞争性的并行重建算法进行了比较,包括基于SENSE的L1‑SENSE算法、TV‑SENSE

9、LpTV‑SENSE10和SENSE‑LORAKS‑TV算22,以及基于SPIRiT模23的STDLR‑SPIRiT24算法。所有算法均用MATLAB实现,实验在配置为英特尔酷睿i9 10900X@3.7 GHz处理器、64 GB内存和64位Windows 10操作系统的服务器上进行。

SENSE模型的重建需要依赖线圈灵敏度编码信息,而事实上准确的灵敏度信息在实际应用中是不容易获得的。在过去的研究中已经提出了多种灵敏度估计的方

2627以及一些自适应的灵敏度估计方28。本文工作使用了ESPIRiT模2629所提供的灵敏度估计方法来估计线圈灵敏度,它依赖于中心自动校准数据(Auto‑calibration signal,ACS),使用特征向量分解的方法来估计灵敏度,是一种较为准确的灵敏度估计方法。所有基于SENSE的对比算法也使用了同样的灵敏度估计方法。

实验分别对比了不同加速因子(Acceleration factor, AF)的二维泊松圆盘采样、一维等间隔采样以及一维高斯随机采样数据的重建结果。为了使用ESPIRiT进行灵敏度估计,在二维泊松圆盘采样中使用24×24的自动校准区域,在一维采样中使用20条自动校准线。

为了验证NLR‑SENSE算法的性能,将所有对比算法在不同的公开MR全采样脑部图像上进行了测试。图2(a~g)来自于文献[

30],均为使用12通道线圈采集的大脑数据,大小为218像素×170像素。图2(h)来自数据集brain_8ch23,是一个使用8通道线圈采集的大脑数据,大小为200像素×200像素。图2(i~k)来自于NYU fastMRI31‑32,是大小为320像素×320像素的大脑数据集。使用线圈压缩技33将15通道采集的数据压缩到8通道用于算法性能的测试。本文对其中4个数据集进行了视觉比较,分别命名为Dataset 1(图2(a))、Dataset 2(图2(e))、Dataset 3(图2(h))和Dataset 4(图2(j))。3种欠采样掩模如图2(l~n)所示。

图2  测试图像与欠采样掩模

Fig.2  Test data and undersampling mask

为了评估不同方法的重建性能,本文使用了3个数值评估指标:信噪比(Signal‑to‑noise ratio,SNR

34、高频误差范数(High‑frequency error norm,HFEN35以及结构相似性(structural similarity,SSIM)指36。其中SNR和SSIM的数值越高,或者HFEN数值越低,则表示图像重建质量越好。需要注意的是,为了更加合理地分析重建质量,实验中仅在感兴趣区域(Region of interest,ROI)计算这些指标。3个评估指标的计算方式如下。

对于全采样的参考图像X^与重建图像X,SNR(dB)的定

34

SNR=10log10VarMSE (14)

式中:MSE表示XX^之间的均方误差;Var为参考图像X^的方差。

HEFN定

35

HFEN=filterX^-filterX2filterX^2 (15)

式中filter·为一个拉普拉斯高斯滤波器,用于捕捉图像边缘。

SSIM定

36

SSIM=2uxux^+c12σxx^+c2ux2+ux^2+c1σx2+σx^2+c1 (16)

式中:uxux^XX^的均值;σx2σx^2XX^的方差;σxx^XX^的协方差;c1c2为常数。

为了进行不同算法的重建比较,所有的模型均调整参数到SNR最佳。对于NLR‑SENSE模型,在非局部低秩去噪过程中,将单线圈图像以5为步进分成重叠的参考块,参考块大小为6×6(n=36),对于每个参考块,匹配出与其最相似的43个块用于去噪(m=43),并调整参数μλ1λ2到SNR最优。

2.2 不同算法的重建视觉比较

首先比较二维泊松圆盘欠采样数据的重建结果。图3展现了使用不同算法对加速因子为3的欠采样的Dataset 1~Dataset 4进行重建的结果。图4展现了使用对加速因子为5的欠采样Dataset 1~Dataset 4进行重建的结果。实验还比较了一维欠采样数据的重建结果。图5展示了加速因子为3的一维等间隔欠采样和一维高斯随机欠采样数据的重建结果。为了使视觉比较更加直观,对重建图像的局部区域进行了放大,并绘制了重建图像与参考图像之间的误差图像。

图3  采用不同算法从加速因子为3的二维泊松圆盘欠采样数据的重建图像与误差图像

Fig.3  Reconstruction images and error graphs of different algorithms when reconstructing from the undersampled data based on 2‑D Poisson‑disc undersampling patterns and 3‑fold acceleration factor

图4  采用不同算法从加速因子为5的二维泊松圆盘欠采样数据的重建图像与误差图像

Fig.4  Reconstruction images and error graphs of different algorithms when reconstructing from the undersampled data based on 2‑D Poisson‑disc undersampling patterns and 5‑fold acceleration factor

图5  采用不同算法从加速因子为3的一维欠采样数据的重建图像与误差图像

Fig.5  Reconstruction images and error graphs of different algorithms when reconstructing from the undersampled data based on 1‑D patterns and 3‑fold acceleration factor

对比图3图4不同加速因子情况下二维泊松圆盘欠采样数据的重建结果,可以看出:(1)加速因子为3时的重建图像明显比加速因子为5的重建图像清晰并且更接近于参考图像;(2)加速因子为3时,各个算法均能较好地重建图像,NLR‑SENSE的重建结果含有最少的伪影和误差;(3)加速因子为5的重建图像质量对比加速因子为3时明显变差,图像的误差也明显变多,此时NLR‑SENSE依然能重建出比较清晰的图像,明显优于其他几种算法。

实验对比了不同采样掩模在加速因子为3对Dataset 1的重建结果。分析重建图像的误差图可以看出,二维泊松圆盘欠采样数据的重建结果具有最少的重建误差,其次是一维等间隔采样、一维高斯随机欠采样数据的重建图像重建误差最多。

从图3~5可以看出,对于不同的数据集以及不同的欠采样模式,L1‑SENSE算法重建出的图像比较模糊,具有较多的图像伪影,且它的重建误差最多;TV‑SENSE、LpTV‑SENSE、SENSE‑LORAKS‑TV以及STDLR‑SPIRiT方法能进一步减少重建伪影和重建误差;本文提出的NLR‑SENSE方法重建的图像最清晰,线条最连贯,重建误差最少,重建图像也最接近参考图像。

2.3 不同算法的数值指标对比

首先比较二维泊松圆盘欠采样数据的重建。表1表2列出了对于Dataset 1~Dataset 4,不同算法对于加速因子3~7的二维泊松圆盘欠采样数据的重建结果。NLR‑SENSE方法的SNR、HFEN和SSIM总体上均优于其他算法,这与重建结果的视觉比较结果一致。此外,对于图2中所有的测试数据,均使用2‑D泊松圆盘和加速因子3~7进行欠采样,并使用所有的比较方法进行图像的重建。本文计算了不同加速因子下NLR‑SENSE算法与比较算法之间的SNR、HFEN和SSIM差值,并计算了不同欠采样模式下的均值,如表3所示。从表3可以看出,与其他算法相比NLR‑SENSE算法的SNR和SSIM均有一定提升,HFEN也一定程度减小了。也就是说,NLR‑SENSE在所有测试数据上获得了最佳的SNR、HFEN和SSIM结果。具体地,对于所有数据和所有欠采样模式,NLR‑SENSE的SNR值在L1‑SENSE、TV‑SENSE、LpTV‑SENSE、SENSE‑LORAKS‑TV和STDLR‑SPIRiT的基础上平均提升了2.06、1.31、1.00、1.23和1.06 dB。

表1  从不同加速因子的二维泊松圆盘采样数据重建Dataset 1和Dataset 2时不同算法重建图像的SNR、HFEN以及SSIM对比
Table 1  SNR, HFEN and SSIM of reconstructed images via different algorithms from the undersampled data based on 2‑D Poisson‑disc undersampling patterns and different AFs for Dataset 1 and Dataset 2
算法指标Dataset 1Dataset 2
AF=3AF=4AF=5AF=6AF=7AF=3AF=4AF=5AF=6AF=7
L1‑SENSE SNR 23.40 21.59 20.27 18.85 18.22 22.02 20.41 19.20 17.98 17.38
HFEN 0.103 7 0.138 9 0.175 4 0.221 2 0.240 8 0.129 5 0.166 5 0.204 0.255 4 0.276 6
SSIM 0.962 0 0.946 5 0.932 6 0.914 8 0.905 3 0.941 8 0.922 8 0.904 8 0.884 9 0.873 4
TV‑SENSE SNR 23.77 22.24 20.87 19.79 19.04 22.60 21.18 20.04 18.63 18.15
HFEN 0.105 5 0.134 2 0.167 5 0.198 7 0.220 5 0.126 3 0.157 4 0.189 7 0.241 0 0.252 6
SSIM 0.966 2 0.955 4 0.943 7 0.931 6 0.922 4 0.952 8 0.939 5 0.925 5 0.908 3 0.898 1
LpTV‑SENSE SNR 24.11 22.54 21.27 20.08 19.28 22.70 21.31 20.24 18.97 18.53
HFEN 0.098 5 0.126 1 0.155 0 0.188 8 0.211 1 0.122 1 0.150 3 0.178 4 0.223 7 0.233 9
SSIM 0.969 3 0.958 2 0.947 9 0.933 1 0.922 7 0.952 3 0.938 9 0.925 2 0.909 4 0.900 3
SENSE‑LORAKS‑TV SNR 23.83 22.41 21.14 19.94 19.25 22.37 21.04 19.97 18.55 18.04
HFEN 0.095 1 0.118 7 0.149 1 0.182 3 0.202 5 0.120 8 0.149 0 0.179 5 0.233 8 0.248 0
SSIM 0.966 3 0.956 1 0.945 6 0.933 1 0.924 8 0.948 8 0.935 2 0.922 3 0.904 2 0.893 9
STDLR‑SPIRiT SNR 24.54 22.91 21.45 20.22 19.34 22.88 21.37 19.94 18.59 18.16
HFEN 0.089 8 0.114 1 0.143 8 0.171 5 0.194 3 0.124 7 0.157 8 0.198 7 0.248 5 0.260 5
SSIM 0.968 9 0.957 7 0.945 2 0.932 8 0.922 0 0.952 5 0.938 1 0.921 7 0.903 2 0.893 6
NLR‑SENSE SNR 25.04 23.75 22.71 21.34 20.53 23.49 22.41 21.52 20.49 19.94
HFEN 0.084 6 0.104 6 0.129 1 0.164 2 0.187 6 0.108 7 0.129 6 0.151 9 0.184 6 0.199 7
SSIM 0.975 2 0.967 7 0.960 5 0.946 0 0.936 8 0.958 9 0.950 0 0.940 0 0.927 8 0.919 0
表2  从不同加速因子的二维泊松圆盘采样数据重建Dataset 3和Dataset 4时不同算法重建图像的SNR、HFEN以及SSIM对比
Table 2  SNR, HFEN and SSIM of reconstructed images via different algorithms from the undersampled data based on 2‑D Poisson‑disc undersampling patterns and different AFs for Dataset 3 and Dataset 4
算法指标Dataset 3Dataset 4
AF=3AF=4AF=5AF=6AF=7AF=3AF=4AF=5AF=6AF=7
L1‑SENSE SNR 23.38 21.66 20.36 19.36 18.47 26.34 24.59 23.60 22.48 21.64
HFEN 0.114 6 0.154 1 0.193 4 0.229 4 0.261 4 0.084 4 0.112 2 0.131 3 0.156 9 0.176 4
SSIM 0.965 8 0.953 5 0.940 4 0.929 4 0.917 5 0.978 1 0.969 5 0.964 4 0.957 5 0.951 9
TV‑SENSE SNR 24.06 22.40 20.97 20.04 19.05 26.80 25.10 24.11 23.07 22.20
HFEN 0.113 9 0.1478 0.184 8 0.211 8 0.246 5 0.084 4 0.107 3 0.124 5 0.142 4 0.159 6
SSIM 0.972 9 0.963 5 0.952 2 0.942 6 0.931 3 0.981 8 0.975 4 0.971 4 0.965 9 0.961 2
LpTV‑SENSE SNR 24.19 22.61 21.32 20.43 19.44 26.92 25.36 24.46 23.61 22.80
HFEN 0.111 2 0.140 5 0.173 2 0.199 0.231 8 0.083 5 0.103 8 0.119 3 0.133 4 0.149 6
SSIM 0.973 7 0.964 6 0.954 3 0.944 5 0.933 0 0.981 3 0.974 8 0.970 7 0.965 8 0.961 1

SENSE‑LORAKS‑

TV

SNR 23.67 22.22 21.01 20.19 19.23 26.54 25.00 24.11 23.13 22.23
HFEN 0.110 0 0.140 6 0.173 1 0.197 7 0.231 4 0.079 0 0.099 7 0.114 7 0.132 9 0.152 5
SSIM 0.968 7 0.959 3 0.949 5 0.941 3 0.931 2 0.979 3 0.972 3 0.968 6 0.963 5 0.958 8
STDLR‑SPIRiT SNR 24.33 22.93 21.50 20.58 19.59 27.07 25.37 24.47 23.55 22.84
HFEN 0.101 1 0.131 5 0.169 2 0.195 7 0.232 0 0.076 1 0.099 1 0.114 3 0.132 2 0.147 7
SSIM 0.971 7 0.963 9 0.953 7 0.943 6 0.932 0 0.979 7 0.972 0 0.967 8 0.962 9 0.959 2
NLR‑SENSE SNR 24.76 23.29 22.13 21.40 20.38 27.69 26.35 25.61 24.92 24.10
HFEN 0.101 5 0.128 9 0.153 8 0.177 9 0.210 9 0.073 4 0.091 7 0.104 1 0.116 9 0.131 7
SSIM 0.976 8 0.969 2 0.961 3 0.953 2 0.941 6 0.984 0 0.978 7 0.975 7 0.971 8 0.967 4
表3  从不同加速因子的二维泊松圆盘采样数据重建图像2所有数据时NLR‑SENSE算法与其他算法的性能提升
Table 3  Performance improvement between NLR‑SENSE and other comparing algorithms for reconstructing all datasets in Fig.2 from the undersampled data based on 2‑D Poisson‑disc undersampling patterns and different AFs
对比算法指标AF=3AF=4AF=5AF=6AF=7
L1‑SENSE SNR 1.38 1.87 2.18 2.37 2.49
HFEN -0.015 7 -0.029 6 -0.042 7 -0.054 3 -0.060 9
SSIM 0.009 5 0.015 6 0.020 9 0.025 1 0.027 9
TV‑SENSE SNR 0.85 1.18 1.40 1.53 1.61
HFEN -0.013 1 -0.020 8 -0.027 5 -0.032 1 -0.035 7
SSIM 0.004 0 0.006 5 0.008 9 0.010 1 0.011 1
LpTV‑SENSE SNR 0.66 0.93 1.07 1.15 1.19
HFEN -0.009 7 -0.014 7 -0.018 9 -0.022 0 -0.023 4
SSIM 0.003 5 0.006 0 0.008 1 0.009 7 0.010 7
SENSE‑LORAKS‑TV SNR 0.99 1.16 1.23 1.33 1.42
HFEN -0.005 7 -0.009 8 -0.014 3 -0.017 8 -0.021 6
SSIM 0.005 5 0.007 6 0.009 2 0.009 6 0.010 3
STDLR‑SPIRiT SNR 0.52 0.95 1.20 1.32 1.30
HFEN -0.009 7 -0.018 7 -0.024 8 -0.029 3 -0.029 7
SSIM 0.005 4 0.009 0 0.011 5 0.013 2 0.013 1

对于加速因子为3~5的一维等间隔和一维高斯随机采样方式,在Dataset 1上对所有的比较方法进行了测试,重建结果的SNR、HFEN和SSIM结果如表4所示。从表中可以看出,对不同采样率的两种一维欠采样数据进行重建时,NLR‑SENSE方法也都得到了最优的SNR、HFEN和SSIM。

表4  从3~5倍加速的一维等间隔采样和一维高斯随机采样数据重建Dataset 1时不同算法重建图像的SNR、HFEN以及SSIM对比
Table 4  SNR, HFEN and SSIM values when reconstructing Dataset 1 from the undersampled data based on 1‑D uniform undersampling and 1‑D Gaussian random undersampling patterns with different AFs
算法指标1‑D uniform1‑D Gaussian random
AF=3AF=4AF=5AF=3AF=4AF=5
L1‑SENSE SNR 19.53 16.20 13.93 16.65 13.90 12.58
HFEN 0.202 5 0.298 5 0.425 4 0.297 5 0.429 6 0.521 0
SSIM 0.931 1 0.883 2 0.842 9 0.903 1 0.852 7 0.829 0
TV‑SENSE SNR 20.99 18.23 15.79 18.50 15.55 13.95
HFEN 0.169 5 0.247 4 0.344 2 0.231 3 0.345 9 0.427 0
SSIM 0.949 3 0.922 5 0.892 9 0.930 8 0.893 4 0.872 0
LpTV‑SENSE SNR 20.96 18.90 17.07 18.91 16.55 15.00
HFEN 0.168 1 0.221 7 0.287 3 0.215 3 0.299 9 0.368 9
SSIM 0.946 1 0.923 5 0.897 6 0.929 9 0.895 2 0.875 0
SENSE‑LORAKS‑TV SNR 21.09 18.94 16.44 18.99 16.18 14.62
HFEN 0.154 3 0.209 9 0.311 7 0.204 8 0.307 4 0.383 7
SSIM 0.946 8 0.926 1 0.893 5 0.934 3 0.896 9 0.876 0
STDLR‑SPIRiT SNR 21.11 19.16 17.70 19.21 17.04 15.84
HFEN 0.181 5 0.230 4 0.275 5 0.219 0 0.289 6 0.336 8
SSIM 0.946 8 0.924 7 0.905 9 0.933 8 0.904 7 0.891 5
NLR‑SENSE SNR 22.51 20.10 18.19 20.42 18.11 16.83
HFEN 0.137 8 0.199 2 0.255 2 0.185 6 0.256 4 0.296 2
SSIM 0.960 2 0.936 6 0.914 3 0.944 4 0.911 5 0.895 7

2.4 算法收敛性分析

在结合非局部低秩的SENSE模型中,可以使用核范数作为秩代理函数求解非局部低秩问题。模型命名为NLR‑SENSE‑baseline,使用ADMM算法求解对应的凸优化问题。使用NLR‑SENSE‑baseline和NLR‑SENSE对5倍加速的二维泊松圆盘欠采样数据重建Dataset 1。图6展示了两种算法的SNR和HFEN与迭代次数的关系。可以看出,NLR‑SENSE算法在第22次迭代时SNR和HFEN均达到了最优,NLR‑SENSE‑baseline在第25次迭代时SNR和HFEN达到了最优,并且NLR‑SENSE比NLR‑SENSE‑baseline重建图像的SNR高0.7 dB。

图6  使用NLR-SENSE和NLR-SENSE-baseline方法从5倍加速二维泊松圆盘欠采样数据重建Dataset 1时SNR和HFEN随迭代次数的变化

Fig.6  SNR and HFEN versus the iteration number when reconstructing the Dataset 1 based on the 2‑D Poisson-disc undersampling pattern and 5-fold acceleration factor using NLR-SENSE and NLR-SENSE-baseline algorithms

NLR‑SENSE与NLR‑SENSE‑baseline相比,只需要在低秩去噪过程中格外进行一个简单的权值计算,因此两种算法的计算量相当。也就是说,使用加权核范数作为秩代理进行低秩问题求解的NLR‑SENSE算法能够更快地达到收敛并得到更好的图像重建结果。

2.5 参数敏感性分析

以5倍加速因子泊松圆盘欠采样的Dataset 1重建为例,分析了参数μλ1λ2变化时NLR‑SENSE重建图像的SNR变化情况,如图7所示。从图7可以看出,当μ=0.2γ1=0.025γ2=5e-3时,重建图像的SNR最优,为22.71 dB。当μ[0.18,0.24]γ1[0.021,0.03]以及γ2[3e-3,10e-3]时,重建图像的SNR均能保证在22.6 dB以上,与最优的SNR差异小于0.1 dB。

图7  使用NLR-SENSE方法从5倍加速二维泊松圆盘欠采样数据重建Dataset 1时参数μγ1γ2的变化对重建图像的SNR影响

Fig.7  SNR value for different parameters μ, γ1 and γ2 when reconstructing the Dataset 1 based on the 2‑D Poisson-disc undersampling pattern and 5-fold acceleration factor

3 结束语

为了进一步提高SENSE模型的并行磁共振成像重建质量,本文结合非局部低秩约束,提出了一种有效的并行磁共振重建算法,称为NLR‑SENSE。其对应的模型使用加权核范数作为秩代理函数,并使用ADMM方法求解。与使用核范数进行求解的NLR‑SENSE‑baseline方法相比,两种方法的计算量相当,NLR‑SENSE收敛更快且重建图像的质量更高。本文将NLR‑SENSE与不同的并行磁共振成像方法进行了比较,包括多种不同的SENSE模型改进方法,如L1‑SENSE、TV‑SENSE、LpTV‑SENSE和NLR‑SENSE和基于SPIRiT模型的STDLR‑SPIRiT方法。仿真实验结果表明,对不同欠采样模式和不同数据集的重建,NLR‑SENSE具有最清晰的重建图像,最少的重建误差,并且在总体上达到了最佳的SNR、HFEN以及SSIM。对加速因子较高的欠采样数据,NLR‑SENSE也能重建出质量较好的图像。总之,NLR‑SENSE方法能有效提升并行磁共振成像重建图像的质量,有利于进一步加快磁共振成像数据采集的速度。

本文提出的算法需要对每个参考块匹配与其相似的块,并对每一个相似块矩阵进行奇异值分解,导致了较大的计算量。由于每个参考块匹配和相似块矩阵的奇异值分解之间是独立的,非常适合使用GPU进行并行计算,在后续的研究中将考虑使用GPU优化重构算法从而使其达到实时重构。

参考文献

1

王秋良. 磁共振成像系统的电磁理论与构造方法[M]. 北京: 科学出版社, 2018. [百度学术] 

WANG Qiuliang. Electromagnetic theory and construction method of magnetic resonance imaging system[M]. Beijing: Science Press, 2018. [百度学术] 

2

LUSTIG M, DONOHO D, PAULY J M. Sparse MRI: The application of compressed sensing for rapid MR imaging[J]. Magnetic Resonance in Medicine, 2007, 58(6): 1182-1195. [百度学术] 

3

杜秀丽, 刘晋廷, 吕亚娜,. 基于迭代p阈值投影算法的压缩感知磁共振成像[J]. 数据采集与处理, 2020, 35(6): 1060-1068. [百度学术] 

DU Xiuli, LIU Jinting, LYU Yana, et al. Compressed sensing magnetic resonance imaging based on projected iterative p-thresholding algorithm[J]. Journal of Data Acquisition and Processing, 2020, 35(6): 1060-1068. [百度学术] 

4

PRUESSMANN K P, WEIGER M, SCHEIDEGGER M B, et al. SENSE: Sensitivity encoding for fast MRI[J]. Magnetic Resonance in Medicine, 1999, 42(5): 952-962. [百度学术] 

5

史振威,雷森. 图像超分辨重建算法综述[J]. 数据采集与处理, 2020, 35(1): 1-20. [百度学术] 

SHI Zhenwei, LEI Sen. Review of image super-resolution reconstruction[J]. Journal of Data Acquisition and Processing, 2020, 35(1): 1-20. [百度学术] 

6

李响,蒋敏,彭钰蘅,. 编码曝光图像的L0正则化去模糊重建方法[J]. 数据采集与处理, 2020, 35(1): 65-78. [百度学术] 

LI Xiang, JIANG Min, PENG Yuheng, et al. Deblurring and restoration method of coded exposure images using L0 regularization[J]. Journal of Data Acquisition and Processing, 2020, 35(1): 65-78. [百度学术] 

7

RUDIN L, OSHER S, FATEMI E. Nonlinear total variation based noise removal algorithms[J]. Physica D: Nonlinear Phenomena, 1992, 60(1): 259-268. [百度学术] 

8

许基隆,陈颖频. 一种基于全变分技术和Lp伪范数的椒盐噪声去除方法[J]. 数据采集与处理, 2020, 35(1): 89-99. [百度学术] 

XU Jilong,CHEN Yingpin. Method of removing salt and pepper noise based on total variation technique and Lp pseudo norm[J]. Journal of Data Acquisition and Processing, 2020, 35(1): 89-99. [百度学术] 

9

RAMANI S, FESSLER J A. Parallel MR image reconstruction using augmented lagrangian methods[J]. IEEE Transactions on Medical Imaging, 2011, 30(3): 694-706. [百度学术] 

10

鲍中文. 磁共振成像的稀疏重构算法研究[D]. 昆明: 昆明理工大学, 2020. [百度学术] 

BAO Zhongwen. Research on sparse reconstruction algorithm of magnetic resonance imaging[D]. Kunming: Kunming University of Science and Technology, 2020. [百度学术] 

11

MANJÓN J V, COUPÉ P, MARTÍ-BONMATÍ L, et al. Adaptive non-local means denoising of MR images with spatially varying noise levels[J]. Journal of Magnetic Resonance Imaging, 2010, 31(1): 192-203. [百度学术] 

12

QU X, HOU Y, LAM F, et al. Magnetic resonance image reconstruction from undersampled measurements using a patch-based nonlocal operator[J]. Medical Image Analysis, 2014, 18(6): 843-856. [百度学术] 

13

EKSIOGLU E M. Decoupled algorithm for MRI reconstruction using nonlocal block matching model: BM3D-MRI[J]. Journal of Mathematical Imaging and Vision, 2016, 56(3): 430-440. [百度学术] 

14

DONG W, SHI G, LI X, et al. Compressive sensing via nonlocal low-rank regularization[J]. IEEE Transactions on Image Processing, 2014, 23(8): 3618-3632. [百度学术] 

15

朱俊, 陈长伟, 苏守宝,. 基于局部和非局部正则化的图像压缩感知[J]. 数据采集与处理, 2016, 31(6): 1148-1155. [百度学术] 

ZHU Jun, CHEN Changwei, SU Shoubao, et al. Image compressed sensing based on local and nonlocal regularizations[J]. Journal of Data Acquisition and Processing, 2016, 31(6): 1148-1155. [百度学术] 

16

DABOV K, FOI A, KATKOVNIK V, et al. Image denoising by sparse 3-D transform-domain collaborative filtering[J]. IEEE Transactions on Image Processing, 2007, 16(8): 2080-2095. [百度学术] 

17

YANG J, ZHANG Y, YIN W. A fast alternating direction method for TVL1-L2 signal reconstruction from partial fourier data[J]. IEEE Journal of Selected Topics in Signal Processing, 2010, 4(2): 288-297. [百度学术] 

18

侯榆青,金明阳,贺小伟,. 基于随机变量交替方向乘子法的荧光分子断层成像[J]. 光学学报, 2017, 37(7): 717001. [百度学术] 

HOU Yuqing, JIN Mingyang, HE Xiaowei , et al. Fluorescence molecular tomography using a stochastic variant of alternating direction method of multipliers[J]. Acta Optica Sinica, 2017, 37(7): 0717001. [百度学术] 

19

GU S, ZHANG L, ZUO W, et al. Weighted nuclear norm minimization with application to image denoising[C] //Proceedings of 2014 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). Columbus, OH, USA: IEEE, 2014: 2862-2869. [百度学术] 

20

CAI J, CANDÈS E J, SHEN Z. A singular value thresholding algorithm for matrix completion[J]. SIAM Journal on Optimization, 2010, 20(4): 1956-1982. [百度学术] 

21

HALDAR J P. Low-rank modeling of local k-space neighborhoods (LORAKS) for constrained MRI[J]. IEEE Transactions on Medical Imaging, 2014, 33(3): 668-681. [百度学术] 

22

KIM T H, SETSOMPOP K, HALDAR J P. LORAKS makes better SENSE: Phase-constrained partial fourier SENSE reconstruction without phase calibration[J]. Magnetic Resonance in Medicine, 2017, 77(3): 1021-1035. [百度学术] 

23

LUSTIG M, PAULY J M. SPIRiT: Iterative self-consistent parallel imaging reconstruction from arbitrary k-space[J]. Magnetic Resonance in Medicine, 2010, 64(2): 457-471. [百度学术] 

24

ZHANG X, GUO D, HUANG Y, et al. Image reconstruction with low-rankness and self-consistency of k-space data in parallel MRI[J]. Medical Image Analysis, 2020, 63: 101687. [百度学术] 

25

AFONSO M V, BIOUCAS-DIAS J M, FIGUEIREDO M A T. Fast image recovery using variable splitting and constrained optimization[J]. IEEE Transactions on Image Processing, 2010, 19(9): 2345-2356. [百度学术] 

26

UECKER M, LUSTIG M. Estimating absolute-phase maps using ESPIRiT and virtual conjugate coils[J]. Magnetic Resonance in Medicine, 2017, 77(3): 1201-1207. [百度学术] 

27

ALLISON M J, RAMANI S, FESSLER J A. Regularized MR coil sensitivity estimation using augmented Lagrangian methods[C]//Proceedings of 2012 9th IEEE International Symposium on Biomedical Imaging (ISBI). Barcelona, Spain: IEEE, 2012: 394-397. [百度学术] 

28

KE L, JINGXIN Z. Dynamic-parallel MR image reconstruction based on adaptive coil sensitivity estimation[C] // Proceedings of 2008 5th IEEE International Symposium on Biomedical Imaging: From Nano to Macro. Paris, France: IEEE, 2008: 1015-1018. [百度学术] 

29

UECKER M, LAI P, MURPHY M J, et al. ESPIRiT—An eigenvalue approach to autocalibrating parallel MRI: Where SENSE meets GRAPPA[J]. Magnetic Resonance in Medicine, 2014, 71(3): 990-1001. [百度学术] 

30

SOUZA R, LUCENA O, GARRAFA J, et al. An open, multi-vendor, multi-field-strength brain MR dataset and analysis of publicly available skull stripping methods agreement[J]. NeuroImage, 2018, 170: 482-494. [百度学术] 

31

ZBONTAR J, KNOLL F, SRIRAM A, et al. fastMRI: An open dataset and benchmarks for accelerated MRI[EB/OL]. (2019-12-11)[2021-06-30]. https://arxiv.org/abs/1811.08839. [百度学术] 

32

KNOLL F, ZBONTAR J, SRIRAM A, et al. fastMRI: A publicly available raw k-space and DICOM dataset of knee images for accelerated MR image reconstruction using machine learning[J]. Radiol Artif Intell, 2020, 2(1): e190007. [百度学术] 

33

BAHRI D, UECKER M, LUSTIG M. ESPIRIT-based coil compression for cartesian sampling[C]// Proceedings of International Society for Magnetic Resonance in Medicine(ISMRM). Salt Lake City: ISMRM, 2013: 2657. [百度学术] 

34

CHARTRAND R. Fast algorithms for nonconvex compressive sensing: MRI reconstruction from very few data[C] // Proceedings of 2009 IEEE International Symposium on Biomedical Imaging: From Nano to Macro. Boston, MA, USA: IEEE, 2009: 262-265. [百度学术] 

35

RAVISHANKAR S, BRESLER Y. MR image reconstruction from highly undersampled k-space data by dictionary learning[J]. IEEE Transactions on Medical Imaging, 2011, 30(5): 1028-1041. [百度学术] 

36

WANG Z, BOVIK A C, SHEIKH H R, et al. Image quality assessment: from error visibility to structural similarity[J]. IEEE Transactions on Image Processing, 2004, 13(4): 600-612. [百度学术]