2023年7月29日发(作者：)

基于小波树和互补分解的CS-MRI重建算法

PEI Ying;ZHU Jin-xiu;YANG Yu-chen;WU Wen-xia

【摘 要】针对压缩感知(CS)核磁共振成像(MRI)重建算法中全变分(TV)正则项会导致图像细节丢失的问题,引入互补分解模型,结合小波树结构稀疏(简称小波树),提出一种基于小波树和互补分解的CS-MRI重建算法.利用互补分解将图像分成平滑分量和残差分量两个部分,并将平滑分量用于TV正则项,残差分量用于e1范数,可避免TV正则项在滤除噪声的同时滤除过多的细节信息;利用小波树结构稀疏可进一步补充小波稀疏等先验信息,减少测量值或提高信噪比.针对目标函数中存在平滑和残差两个未知分量,将目标函数分解为相应的两个子问题交替最小化进行求解.实验结果表明,与基于小波树的WaTMRI和基于TV的TVCMRI、FCSA等重建算法相比,其能在滤除噪声的同时有效改善MRI图像的细节信息.

【期刊名称】《计算机技术与发展》

【年(卷),期】2018(028)012

【总页数】5页(P152-156)

【关键词】核磁共振成像;压缩感知;互补分解;小波树结构稀疏(小波树);目标函数;重建算法

【作 者】PEI Ying;ZHU Jin-xiu;YANG Yu-chen;WU Wen-xia

【作者单位】;;;

【正文语种】中 文

【中图分类】TP301.6 0 引 言

核磁共振成像(magnetic resonance imaging，MRI)是医学成像，应用广泛。目前MRI应用的关键在于快速成像。Nyquist采样定理需两倍带宽，不符合实际应用[1]，现常用方法是在压缩感知(compressed sensing，CS)[2-3]框架下，从欠采样k空间中重建MRI数据，能有效减少采样时间，达到快速成像的目的。

关于CS-MRI重建算法的研究有很多。例如，Lusting等[4]利用全变分(total

variation，TV)和小波构建目标函数，提出共轭梯度算法(CG)，但重建时间有待提高；FISTA算法[5]通过计算更合适的起始点以加快收敛速度；TVCMRI[6]、RecPF[7]和FCSA[8]算法利用算子分割、变量分割思想求解联合正则算子，提高重建速度和质量；文献[9]利用图像在频域的特性优化测量矩阵并提出迭代加权算法，提高了重建精度；文献[2,10]利用MR图像低秩特性进行奇异值分解，但过程过于复杂；文献[11-12]提出结构稀疏理论，图像不仅在小波域有稀疏性，其小波稀疏系数也有特定的四叉树结构，在此基础上文献[13]提出YALL1算法，利用小波树结构稀疏代替小波稀疏构建目标函数，以提高图像稀疏度；文献[14]提出WaTMRI算法，联合小波树结构稀疏和小波稀疏分别拥有的结构稀疏和稀疏性，联合改善图像质量；Park等[15]提出互补分解，将完整图像分为平滑和残差两个分量，仅将平滑分量用于TV，残差分量用于1范数，以解决全变分导致的细节过平滑问题；文献[16]利用贪婪算法提高重建速度，但需要确定图像稀疏度，缺乏实际性；文献[17]基于p范数构建重建算法，计算较为复杂。

重建算法常用小波稀疏和TV构建目标函数，但TV会在滤除噪声的同时平滑图像的边缘纹理信息。因此，文中提出一种基于小波树和互补分解的CS-MRI(MRI

based on complementary dual decomposition and wavelet tree sparsity，DualWaTMRI)重建算法,利用互补分解抑制TV造成的细节过平滑问题，小波树结构稀疏用于补充小波稀疏的先验信息，并与同类重建算法进行了对比。

1 基于小波树和互补分解的CS-MRI重建算法

基于小波树和互补分解的CS-MRI重建算法系统框架如图1所示。首先欠采样输入数据得到测量值b，将其傅里叶逆变换预处理生成初始重构图像，再对其互补分解，生成初始残差分量和平滑分量，分别进行TV和1范数处理，同时利用结构化稀疏理论[13]，对整幅图像进行小波树结构稀疏处理，最后与最小二乘拟合共同构建目标函数。同时，提出DualWaTMRI重构算法进行求解，分别得到重构的残差分量和平滑分量，由两分量叠加组成最终的重建图像。

图1 算法系统框架

1.1 基于小波树和互补分解的目标函数

1.1.1 小波树

小波树结构稀疏是图像的小波稀疏系数呈四叉树结构分布，称大尺寸处的系数值是父系数，相邻小尺寸的系数值是子系数，则父子系数间有依赖性，即当某一位置的小波稀疏系数值很大(或接近于0)时，同位置邻近尺度的小波稀疏系数值也很大(或接近于0)。根据该特性将系数分组，生成小波树结构稀疏正则项其中φ表示小波稀疏基，x表示MR图像，ζ表示小波树形分组的集合，g表示其中一组。

1.1.2 互补分解

考虑到对整幅图像进行TV处理会同时平滑噪声和图像的细节信息，引入互补分解模型[13]，将MR图像看成是平滑分量和残差分量的叠加：

x=L+S

(1)

其中，x表示MR图像；L表示灰度值缓慢变化的平滑分量，一般代表图像背景信息；S表示灰度值快速变化的残差分量，一般代表图像中前景信息。

为避免TV过平滑图像，仅对L进行TV处理，S进行小波稀疏处理。文中利用边缘检测算子Sobel[18]对图像进行互补分解，令初始残差分量S0为分解得到的边缘信息，初始平滑分量L0为0。Sobel算子通过核与像素值做卷积和运算，选取合适的阈值提取边缘信息，阈值计算方法为T=μ*255，其中μ是阈值比。

1.1.3 目标函数

综合互补分解和小波树结构稀疏的优势，提出基于小波树和互补分解的CS-MRI重建算法，构建目标函数如下：

(L+S)=

α‖L‖TV+β(‖φS‖1+

(2)

其中,L，S分别是平滑和残差分量；A是测量矩阵；b是测量值；φ是小波稀疏基；ζ是小波树形分组的集合；g表示其中一组；α和β是调优参数，用于平衡所占比重。式2中第一项最小二乘拟合项用于保证重建图像的准确度，第二项平滑分量的TV项用于抑制噪声，避免图像过平滑，第三项残差分量的1范数用于保证小波稀疏，第四项小波树结构稀疏项用于保证图像的结构稀疏性。这四个正则项相互补充，增加图像先验信息，提高算法的鲁棒性。

1.2 DualWaTMRI重构算法

针对式2中的L和S两个未知量，利用交替最小化方法构建DualWaTMRI重构算法进行求解，将目标函数分为L子问题(此时S为固定值)和S子问题(此时L为固定值)，同时将这两个子问题交替迭代，最后重建图像x即是求解得到的L和S值之和。

1.2.1 L子问题求解

求解L子问题时，假设S为固定值，则最小化求解时可省略常数项。L子问题公式表示如下： (3)

为简便求解小波树结构稀疏正则项，令zp=GφL，其中G是ζ的二值矩阵[15]，式3改为：

(4)

其中，λ是调优参数；ζ中共有s组g，gi表示第i组g，即第i个具有父子依赖性的小波系数组，i=1,2,…,s。式4中亦包含zp和L两个未知量，采用交替最小化进一步细化：

(1)L子问题中zp未知量按gi分组求解：

zpgi=

β‖zpgi‖2)

(5)

上式可采用分组软阈值法求解，令rpi=(GφL)gi：

zpgi=

i=1,2,…,s

(6)

由上式求解zpgi值，将其按位置线性组合生成zp。

(2)L子问题中L未知量求解公式如下：

其中，zp值由式6解得。式7可采用FISTA[5]算法求解得到L值。

1.2.2 S子问题求解

与L子问题类似，求解未知量S时假设L为固定值，求解公式如下：

(8)

同样的，令zr=GφS，上式更新为：

(9)

(1)未知量zr可通过分组软阈值法得到，令rri=(GφS)gi，求解公式为：

zrgi=

i=1,2,…,s

(10)

(2)未知量S求解公式如下：

(11)

上式亦可采用FISTA方法求解。

1.2.3 重建算法总结

针对目标函数，利用交替最小化将函数分为两个子问题求解，结合上述重建步骤构建算法DualWaTMRI，该算法重建步骤总结如下：

输入：最大迭代次数N，k空间测量值b，阈值比μ。

步骤1：初始化。

步骤1.1：对b进行傅里叶逆变换预处理，生成初始图像x0；

步骤1.2：对x0取合适的阈值μ互补分解得到初始残差分量S0；

步骤1.3：令初始平滑分量L0=0，n为迭代次数，初始化n=1，ρ=1/Lf。

步骤2：交替最小化求解。

步骤2.1：求解L子问题：给定Sn-1，结合Ln-1求解式5得到zp值，将zp值代入式7求解Ln；

步骤2.2：求解S子问题：给定Ln，结合Sn-1求解式10得到zr值，将zr值代入式11求解Sn；

步骤2.3：令n=n+1，判断n是否等于N，若不等于，回到步骤2.1，否则跳到步骤2.4；

步骤2.4：重建图像是平滑分量和残差分量的叠加：x=L+S(此时L=Ln，S=Sn)。

2 实验结果与分析

对Chest，Heart，Brain和Shoulder等四幅标准MRI图像进行实验。将DualWaTMRI与CG[4]、TVCMRI[6]、FCSA[8]和WaTMRI[15]算法进行比较。根据经验设α=0.02，β=3.5e-2，λ=β/2，Lf=1。同时选小波为稀疏基，伪高斯[6]为采样模板，测量矩阵为采样模板下的部分傅里叶矩阵。结果如表1所示。

表1 四幅MR图像进行重建时的PSNR结果 dB标准图像CGTVCMRIFCSAWaTMRIDualWaTMRIHeart8.827 114.354 917.950518.554

520.2540Chest10.439 813.808 215.539 916.316 716.953 4Shoulder10.017

915.676 520.410 521.779 422.388 3Brain8.706 813.699 515.938 416.950

718.649 0

2.1 实验参数：阈值比μ的选取

对DualWaTMRI算法进行仿真，阈值比μ取值从0到1，结果如图2所示。

从图2可看出，在μ取值为0.05～0.20时，SNR值急剧上升，在μ大于0.20时，SNR值稍有减少，最后趋于平缓。可以得出，在DualWaTMRI算法互补分解模型中μ取值0.16。

图2 仿真结果对比

2.2 重构性能的比较

将MR图像分别通过CG、TVCMRI、FCSA、WaTMRI和DualWaTMRI等算法进行仿真，结果分别如图3、4和表1所示。

图3 四幅MR图像的SNR比较 图4 四幅MR图像的t比较

图3和图4给出了四幅MR图像通过不同算法得到的仿真结果折线图(SNR和t)。表1给出了SNR具体值。从图3和表1中可以看出，四幅MR图像通过DualWaTMRI算法得到的信噪比最高，比WaTMRI算法分别高出约1.69 dB、0.63 dB、0.60 dB和1.70 dB，更优于其他算法。但图4中DualWaTMRI算法上升趋势最慢，表明重建时间长，这是由于目标函数中有两个未知量，需要交替求解，因此文中算法通过牺牲少量计算时间，达到增加信噪比的目的。

综上，DualWaTMRI算法与同类CS-MRI算法相比，虽然由于先验信息丰富导致图像的重建时间有所增加，但重建图像的质量有一定改善。

3 结束语

为改善TV造成的图像细节模糊问题，引入互补分解模型，与小波树结构稀疏相结合，提出一种基于小波树和互补分解的CS-MRI重建算法。利用互补分解模型将MR图像分为平滑和残差两个分量，仅将平滑分量用于TV，残差分量用于1范数，利用两个分量各自的特性，很好地保留了图像的细节信息，同时利用树结构稀疏特性，丰富图像的先验信息，使得在同等数目的扫描数据下，能得到更好的重建图像，或者只需更少的扫描数据，就能获得同等质量的重建图像。实验结果表明，与现有的基于整幅图像全变分、基于小波树的算法相比，该算法以部分运算时间为代价很好地改善了MR图像的重建质量。

参考文献：

【相关文献】

[1] 蒋 沅,苗生伟,罗华柱,等.L_p范数压缩感知图像重建优化算法[J].中国图象图形学报,2017,22(4):435-442. [2] 王郗雨,杨晓梅,胡学姝.基于奇异值分解的压缩感知核磁共振图像重构算法[J].计算机应用研究,2013,30(4):1247-1249.

[3] 赵 扬,王敬时,王 伟,等.基于压缩感知的磁共振血管快速成像的研究进展[J].医疗卫生装备,2016,37(3):104-109.

[4] LUSTIG M,DONOHO D,PAULY J MRI:the application of compressed sensing

for rapid MR imaging[J].Magnetic Resonance in Medicine,2007,58(6):1182-1195.

[5] BECK A,TEBOULLE M.A fast iterative shrinkage-thresholding algorithm for linear

inverse problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202.

[6] MA Shiqian,YIN Wotao,ZHANG Yin,et efficient algorithm for compressed MR

imaging using total variation and wavelets[C]//IEEE conference on computer vision and

pattern age,AK,USA:IEEE,2008:1-8.

[7] YANG Junfeng,ZHANG Yin,YIN Wotao.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.

[8] HUANG Junzhou, ZHANG Shaoting, METAXAS ent MR image reconstruction for

compressed MR imaging[J].Medical Image Analysis,2010,15(5):670-679.

[9] 袁 静.基于压缩感知的核磁共振成像重构算法[J].计算机工程,2015,41(10):270-274.

[10] 陈思吉,杨晓梅,吕雪霜.基于稀疏和低秩先验分离的快速动态MRI重建[J].计算机应用研究,2016,33(10):3196-3200.

[11] HEGDE C,INDYK P,SCHMIDT imation algorithms for model-based

compressive sensing[J].IEEE Transactions on Information Theory,2015,61(9):5129-5147.

[12] 孙 洪,张智林,余 磊.从稀疏到结构化稀疏:贝叶斯方法[J].信号处理,2012,28(6):759-773.

[13] DENG Wei,YIN Wotao,ZHANG sparse optimization by alternating direction

method[C]//SPIE optical engineering+ applications.[s.l.]:International Society for Optics

and Photonics,2013.

[14] CHEN Chen,HUANG ting the wavelet structure in compressed sensing

MRI[J].Magnetic Resonance Imaging,2014,32(10):1377-1389.

[15] PARK S,PARK ssed sensing MRI exploiting complementary dual

decomposition[J].Medical Image Analysis,2014,18(3):472-486.

[16] 李少东,裴文炯,杨 军,等.贝叶斯模型下的OMP重构算法及应用[J].系统工程与电子技术,2015,37(2):246-252.

[17] LI Zechao,TANG rvised feature selection via nonnegative spectral

analysis and redundancy control[J].IEEE Transactions on Image

Processing,2015,24(12):5343-5355.

[18] 沈德海,张龙昌,鄂 旭.基于Sobel算子的改进边缘检测算法[J].信息技术,2015,39(4):15-18.