CN104267361A - 基于结构特征的自适应定量磁化率分布图复合重建的方法 - Google Patents

Abstract

基于结构特征的自适应定量磁化率分布图复合重建的方法，涉及定量磁化率成像。基于幅值图像结构特征先验的磁化率分布图基础重建，重建模型包含具有压缩感知特性的保真项和具有稀疏特性的幅值先验正则化约束项，并且通过添加感兴趣区域从幅值图像中提取二值加权矩阵，对原始磁化率分布进行二值加权；基于磁化率分布结构特征的磁化率分布图复合重建，重建模型包含最小二乘保真项、由基础重建获得的磁化率分布图构造的结构特征正则化约束项以及改善重建磁化率分布效果的平滑项，磁化率结构特征先验定义为3D图像数据在三个方向上的一阶梯度信息；对于l1范数优化问题，采用迭代阈值法进行处理；之后基于l2范数的凸函数性质，用共轭梯度法进行求解。

Description

基于结构特征的自适应定量磁化率分布图复合重建的方法

技术领域

本发明涉及定量磁化率成像，尤其是涉及一种基于结构特征的自适应定量磁化率分布图复合重建的方法。

背景技术

定量磁化率成像(Quantitative Susceptibility Mapping,QSM)利用梯度回波数据的相位信息产生组织的磁场特性图[1]。磁化率和磁场分布关系在本质上是非局部的，它依赖于磁化率的空间分布和相对于主磁场的方向。磁化率信息以前常被认为是无用的图像对比信息，甚至会导致信号丢失，扭曲变形和成像伪影。早期研究者致力于有效计算出因体磁化率的任意分布产生的感应磁场的不均匀性而导致的主磁场扰动，这些研究为QSM的提出奠定了基础[2]。QSM通过求解关于感应磁场分布的病态逆问题，由测得的相位图像导出潜在的组织体磁化率的定量图，使得对组织本质上的磁化率对比成像成为可能。QSM可视为磁敏感加权成像(Susceptibility Weighted Imaging,SWI)的进一步发展，它克服了SWI的非量化特性以及缺乏特异性的问题，例如SWI无法区分钙沉积还是铁沉积。

QSM是一种定量成像技术，是近年来磁共振研究的热点之一，在医学磁共振成像领域具有诱人的应用前景[3]。随着磁共振成像设备硬件的进步，如超高场MRI(7T)、多通道线圈阵列和高级图像处理技术的发展，进一步推动了QSM的研究。该技术能够对局部组织的磁化率特性进行量化研究，有助于分辨并量化特定生物标记，如铁、钙、钆和超顺磁性氧化铁纳米颗粒。磁化率图凭借其能够检测出铁沉积，神经组织脱髓鞘作用和神经退行性疾病中的连通性中断，创伤性脑损伤中的微小出血，脑卒中的出血和缺血性中风的出血转化风险，骨骼矿化，以及动脉粥样硬化斑块成分和脆弱性的优势，有望在手术风险评估和指导中发挥重要作用。QSM还能够对活体组织神经束和脑白质纤维进行非入侵式检测，为神经纤维成像研究提供新的对比机制，这对于神经影像学定量的连接性研究和生物物理研究非常重要。

QSM是近年来兴起并得到快速发展的一个研究方向，目前国际上主要有德国Jena大学的Reichenbach J.，英国Nottingham大学的Bowtell R.，美国Cornell大学的Wang Y.，美国JohnHopkins医学院的Craig J.，美国国家健康研究院的Duyn J.和美国Duke大学的Liu C.L.等为首的研究小组开展相关研究。首先，QSM中需要解决的一个重要问题是如何更好地对相位图进行校正，包括两个主要步骤：相位解缠绕及不均匀背景场的去除。目前已有的相位解缠绕方法及相关论文中，Jenkinson 2003提出的相位区域扩展方法广为应用[4]，Witoszynskyj2009是最快的相位解缠绕方法[5]。受物理原理启发，近来研究人员提出两个背景场消除技术：相位数据复杂谐波伪影抑制算法(SHARP)[6]和偶极场投影方法(PDF)[7]，这两个方法都依赖于对内部和外部磁化率源的正确分辨。其次，定量磁化率图重建是一个受到伪影影响的不适定逆问题，奇异角的存在使得由校正后的相位图像和单位偶极核矩阵直接反演求取磁化率图存在困难。为此，研究者们提出了几种不同的方法，其中有多方向测量的过采样法(COSMOS)[8]，基于分片常数先验[9]或是从幅值图像提取边缘信息先验[10,11]的正则化方法。以上方法可以在不同程度上克服重建伪影，但在重建对比度和准确性方面还有待进一步的提高。

目前，国内磁化率对比成像的研究主要集中于SWI，QSM的研究还有待全面开展。华东师范大学李建奇与美国Cornell大学的Wang Y.合作，采用Siemens3T成像仪研究了定量磁化率成像中磁化率值的方向依赖性。浙江大学生物医学工程金朝阳、杜一平等对SWI数据采集、重建和静脉增强等问题做了较为系统和深入的研究。山东医学影像学研究所、上海交通大学附属仁济医院、浙江大学医学院、天津医科大学和厦门市第二医院等对SWI在脑部恶性肿瘤、脑内血管畸形、脑铁代谢、脑铁沉积、脑梗死等的诊断应用开展了相关临床研究实验。《磁共振成像》杂志近两年来也报道了SWI的原理综述和临床应用研究进展。作为一项新兴的成像技术，QSM在背景场的去除、定量磁化率图重建、磁化率张量成像以及临床应用的探索等方面仍有大量工作有待开展和改进，但是QSM已经展现出广阔的研究前景和重要的临床诊断价值。

参考文献：

[1]Shmueli,K.,de Zwart,J.A.,van Gelderen,P.,Li,T.Q.,Dodd,S.J.,Duyn,J.H.,2009.Magnetic susceptibility mapping of brain tissue in vivo using MRI phase data.Magn.Reson.Med.62,1510-1522.

[2]Marques,J.P.,Bowtell,R.,2005.Application of a Fourier-based method for rapidcalculation of field inhomogeneity due to spatial variation of magnetic susceptibility.ConceptsMagn.Reson.B Magn.Reson.Eng.25B(1),65-78.

[3]Reichenbach,J.R.,2012.The future of susceptibility contrast for assessment of anatomyand function.NeuroImage 62,1311-1315.

[4]Jenkinson,M.,2003.Fast,automated,N-dimensional phase-unwrapping algorithm.MagnReson Med 49,193-197.

[5]Witoszynskyj,S.,Rauscher,A.,Reichenbach,J.R.,Barth,M.,2009.Phase unwrapping ofMR images using Phi UN--a fast and robust region growing algorithm.Med Image Anal13,257-268.

[6]Schweser,F.,Deistung,A.,Lehr,B.W.,Reichenbach,J.R.,2011.Quantitative imaging ofintrinsic magnetic tissue properties using MRI signal phase:an approach to in vivo brain ironmetabolism？NeuroImage 54,2789-2807.

[7]Liu,T.,Khalidov,I.,de Rochefort,L.,Spincemaille,P.,Liu,J.,Tsiouris,A.J.,Wang,Y.,2011.A novel background field removal method for MRI using projection onto dipole fields(PDF).NMR Biomed.24,1129-1136.

[8]Liu,T.,Spincemaille,P.,de Rochefort,L.,Kressler,B.,Wang,Y.,2009.Calculation ofsusceptibility through multiple orientation sampling(COSMOS):a method for conditioning theinverse problem from measured magnetic field map to susceptibility source image in MRI.Magn.Reson.Med.61,196-204.

[9]de Rochefort,L.,Brown,R.,Prince,M.R.,Wang,Y.,2008.Quantitative MR susceptibilitymapping using piece-wise constant regularized inversion of the magnetic field.Magn.Reson.Med.60,1003-1009.

[10]Liu,T.,Liu,J.,de Rochefort,L.,Spincemaille,P.,Khalidov,I.,Ledoux,J.R.,Wang,Y.,2011.Morphology enabled dipole inversion(MEDI)from a single-angle acquisition:comparisonwith COSMOS in human brain imaging.Magn.Reson.Med.66,777-783.

[11]de Rochefort,L.,Liu,T.,Kressler,B.,Liu,J.,Spincemaille,P.,Lebon,V.,Wu,J.,Wang,Y.,2010.Quantitative susceptibility map reconstruction from MR phase data using bayesianregularization:validation and application to brain imaging.Magn.Reson.Med.63,194-206.

发明内容

本发明的目的在于提供一种基于结构特征的自适应定量磁化率分布图复合重建方法。

本发明包括如下步骤：

1)基于幅值图像结构特征先验的磁化率分布图基础重建(M-Step)，重建模型为

${\mathrm{\χ}}_{\mathrm{mag}}=\underset{{\mathrm{\χ}}_{\mathrm{mag}}}{\arg \min }{\left|\right|\mathrm{\χ}\left(k\right)\mathrm{gH}-\mathrm{diag}\left(H\right)\mathrm{gF\χ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|P_{\mathrm{mag}}g\▿\mathrm{\χ}\left|\right|}_1+\mathrm{\δ}{\left|\right|W_{\mathrm{mag}}\mathrm{g\χ}\left|\right|}_2$

包含具有压缩感知特性的保真项和具有稀疏特性的幅值先验正则化约束项并且通过添加感兴趣区域从幅值图像中提取二值加权矩阵W_mag，对原始磁化率分布进行二值加权，减小计算量的同时可以改善磁化率分布重建效果；

2)基于磁化率分布结构特征的磁化率分布图复合重建(S-Step)，重建模型为

包含最小二乘保真项由基础重建(M-Step)获得的磁化率分布图构造的结构特征正则化约束项以及改善重建磁化率分布效果的平滑项磁化率结构特征先验定义为3D图像数据在三个方向上的一阶梯度信息；

所述重建模型中，χ表示待求解的空间离散的磁化率分布，λ和δ为可调节的系数，为梯度算子，C为卷积核矩阵，表示相位图，F表示离散傅里叶变换；对于l1范数优化问题，采用迭代阈值法(Iterative threshold algorithm，ITA)来进行处理；之后基于l2范数的凸函数性质，采用共轭梯度法进行求解，以保证收敛速度。

在步骤1)中，所述保真项采用区域自适应截值法算出的初始磁化率值的傅里叶空间数据χ(k)，并对其进行三维二值加权处理，舍弃那些会带来重建伪影的较小值，类似于压缩感知理论中的欠采样过程；所述幅值先验正则化约束项代表3D图像数据分别在三个方向的一阶梯度信息；由于图像在全变分(total variation,TV)变换域内是稀疏的，因此设计基于l2范数的全变分平滑约束；约束项并不对所有的值都进行处理，而是保留在幅值图像内边缘梯度较大的值，将其权重设置为0，对于那些幅值梯度较小的区域权重设置为1，进行稀疏域内的l1范数最小化处理，以实现幅值结构特征先验的约束；加权矩阵P_mag的提取目的是使重建磁化率分布图与幅值图在结构上趋于一致，消除磁化率分布图在平滑区域的细小结构伪影，重建出磁化率分布图中的梯度较大的边界信息；二值加权矩阵W_mag将感兴趣区域权重设置为1，不需要计算的区域权重设置为0。

在步骤2)中，所述保真项在变换域区分开，选取时域内的保真项，以使整个求解过程始终在可行域内；所述约束项采用三维一阶梯度差分算子对原始磁化率进行处理，代表着图像三个方向的一阶梯度信息，相比于基础重建中所使用的幅值先验加权，复合重建中所使用的新的从初步重建磁化率分布图中提取的先验加权矩阵能够为复合重建提供更加准确的先验结构信息；同时为改善重建磁化率分布效果，引入从初步重建的磁化率分布图中提取的新的二值加权矩阵其获得方式与初步重建时一致，只是对于获取先验信息的对象由幅值图变为基于幅值先验重建出的磁化率分布图。

本发明提供一种基于结构特征的自适应定量磁化率分布图的复合稀疏重建方法。已知在磁共振成像中，幅值图与磁化率分布图的结构信息具有整体相似性的同时也存在着局部差异。以此为依据，该方法首先利用幅值结构特征先验作为正则化约束项，实现磁化率病态逆问题的初步求解；然后在此基础上，从初步重建出的磁化率分布图中提取磁化率的结构特征构造更准确的先验加权矩阵，从而建立基于磁化率结构特征分布的自适应定量磁化率分布图复合重建模型。与基于幅值特征先验的磁化率重建方法相比，复合重建模型能够有效抑制由于幅值结构特征与磁化率结构特征的局部不一致性在初步重建过程中可能引入的先验加权误差，以确保更准确的计算出定量磁化率分布图。

附图说明

图1为数值仿真图。a)为真实磁化率分布图，b)为幅值图像，c)为相位图像。

图2为仿真图实验结果。a)为幅值图像，b)为基于幅值结构特征先验加权重建磁化率分布图，c)为从幅值图像提取的位置加权矩阵图像，其中黑色的地方值为0，白色的地方值为1，d)为从重建磁化率图像中提取的二值加权矩阵图像，其中黑色的地方值为0，白色的地方值为1。

图3为仿真图实验结果对比。a)为FDRI方法重建磁化率分布图像，b)为MEDI方法重建磁化率分布图像，c)为RASF方法重建磁化率分布图像。

图4为仿真图实验定量分析区域选取。标记A白线所在位置为采样分析区域。

图5为仿真图实验定量分析曲线对比。

图6为人脑磁化率分布图37层、39层和41层重建结果差值图对比。a)为幅值图，b)为RASF重建的磁化率分布图，c)为RASF重建结果与FDRI重建结果差值图，d)为RASF重建结果与MEDI重建结果差值图。

具体实施方式

本发明所提出的基于结构特征的自适应定量磁化率分布图复合重建的方法，能重建出更多的结构信息，提供更为精确的重建磁化率分布定量分析。方法具体实施过程如下：

首先，进行数值仿真实验。在一个模拟的三维体数据中放置一个圆盘，如图1.a)所示，圆盘中均匀放置的8个磁化率球体，每个磁化率球体的磁化率值服从0.5～4ppm间线性分布，在球体中间放置了3个彼此相互垂直的柱体，其磁化率值为0.5ppm。图1.b)所示为仿真数据对应的幅值图像，其中黑色箭头所指向的球体为模拟大脑中出血点，其幅值为0，白色箭头所指向的球体幅值比周边背景高30％，3个柱体的幅值比周边背景高10％，其余6个球的幅值设置为均匀统一大小，并且为周围背景幅值的二倍。图1.c)为对应的相位图像。

从图2.a)和2.b)中可以看出，在白色箭头所指向的位置处，幅值图像中的垂直柱体并没有很清晰的轮廓，而在基于幅值结构特征先验加权重建的磁化率分布图中，该柱体较为清晰的显示出来。图2.c)和2.d)分别为幅值图像与基于幅值结构特征先验加权重建的磁化率分布图所提取的先验加权矩阵图像。在图2.c)和2.d)中，黑色箭头所显示的位置充分说明了由基于幅值结构特征先验加权重建的磁化率分布图提取的结构先验信息要明显优于由幅值图像提取的结构先验信息。由基于幅值结构特征先验加权重建的磁化率分布图提取的结构先验信息具有更加清晰的边界信息，以及在磁化率变化较平缓的区域的确表现出更加良好的平滑性，这样就可以从中提取更为准确的结构先验信息进行复合磁化率分布图重建。

选取Berkin Bilgic在2012年提出的FDRI(Field-Dependent Relaxation Rate Increase)算法的磁化率分布图重建方法以及MEDI(Morphology Enabled Dipole Inversion)方法作为对比，以说明本发明提出的基于幅值及磁化率分布结构特征先验加权的复合磁化率分布稀疏重建方法(RASF,Regional Adaptivity and Structural Features)在磁化率分布图重建中的良好效果。

从图3中白色箭头所指向的位置处可以看出，在三个垂直柱体位置处的磁化率分布重建中，RASF方法重建磁化率分布图效果要明显优于FDRI方法重建磁化率分布图与MEDI方法重建磁化率分布图，RASF方法重建磁化率分布图的边缘结构更为清晰，对比度更好。从白色矩形所圈位置处可以看出，RASF方法重建磁化率分布图也要明显优于FDRI方法及MEDI方法重建磁化率分布重建图，RASF方法磁化率分布重建图在对背景噪声的抑制效果上有着明显的优势。

为了验证本方法在磁化率分布图重建定量分析中的效果，选取磁化率分布图中一些体素点来进行验证分析，图4中所示白线位置为选取的重建磁化率图的体素点，将它们的磁化率值分布曲线绘制出并与真实磁化率体素进行对比。如图5所示，从图中黑色箭头所指向的位置可以看出，RASF方法重建的磁化率分布图要更加接近于真实磁化率分布，表现出更为良好的重建定量精确性，基于FDRI方法重建的磁化率分布图中出现了较大的磁化率值体素点，而MEDI方法重建磁化率分布图在某些特定的位置上会出现磁化率值浮动较大的现象，所以RASF方法重建的磁化率分布图在定量性上更加接近真实的磁化率值，表现出了良好的定量磁化率分布重建性质。

同时采用信噪比(SNR)、对比噪声比(CNR)、峰值信噪比(PSNR)和平均结构相似度(MSSIM)作为评价参数对图像质量做出定量评价，其中SNR定义为：

$\mathrm{SNR}=10\log (\mathrm{Var}\left(x\right)/\mathrm{Var}(\hat{x}-x))$

其中x表示原始信号，表示加入噪声后的信号，Var表示数据方差。

CNR是一种常用的医学图像质量评定标准，它能够有效地反应图像的视觉对比度，定义为：

$\mathrm{CNR}=\frac{\stackrel{\‾}{X_1}-\stackrel{\‾}{X_2}}{{({{\mathrm{\σ}}^2}_{X_1}-{{\mathrm{\σ}}^2}_{X2})}^{1/2}}$

其中分别表示载波信号和噪声的平均值，σ表示数据的标准差。

PSNR是原图像与被处理图像之间的均方误差相对于(2^n-1)^2的对数值(信号最大值的平方，n是每个采样值的比特数)，定义为：

$\mathrm{PSNR}=10\log \left\{\frac{\max {\left(X\right)}^2}{\frac{1}{N}\underset{i=1}{\overset{N^2}{\mathrm{\Σ}}}{(X\left(i\right)-\hat{X}\left(i\right))}^2}\right\}$

其中X表示原图像，表示处理后的图像。

MSSIM是两幅图像平均结构相似度的一个客观评价，它的取值取决于三个因素：图像亮度的相似程度、图像对比度的相似程度和图像结构的相似程度，在已知真实图像时，可以通过计算两幅图像的MSSIM评价图像重建质量。从表1中可以看出，SNR、PSNR、CNR、MSSIM这些评价指标中RASF方法磁化率分布重建图都要优于FDRI方法及MEDI方法重建磁化率分布图。

表1：仿真图实验结果图像质量分析对比(单位：dB)

随后，进行人脑数据实验，该数据是从通用电气3T成像仪中获得。采样序列为多回波，成像参数为TE＝2.8～44.2ms，TR＝50～60ms，成像视野为24×24cm，层厚为2mm，获得的脑部扫描数据维度为256×256×74。在对原始磁化率分布数据进行复合重建后，我们给出RASF方法重建结果与FDRI、MEDI方法结果差值图对比，如图6所示。从图6.c)这一列中可以看出，对于没有结构先验约束的FDRI方法重建的磁化率分布图，RASF在重建结构上明显对结构进行了加强，突出了更多的结构信息，具有更高的对比度，从图6.d)这一列中可以看出，相比于MEDI方法，RASF方法重建的结构在苍白球及额叶位置处提高了重建的准确性，这说明利用磁化率分布图提取的结构先验要优于从幅值图提取的结构先验信息，达到了更好的重建效果。同时可以看到，在图6.d)这一列大脑边缘脑沟回处也有着高亮的区域，这说明利用磁化率分布图提取结构先验信息对大脑边缘沟回处的磁化率分布重建也有着明显的增强效果。

Claims (4)

1.基于结构特征的自适应定量磁化率分布图复合重建的方法，其特征在于包括如下步骤：

1)基于幅值图像结构特征先验的磁化率分布图基础重建(M-Step)，重建模型为

${\mathrm{\χ}}_{\mathrm{mag}}=\underset{{\mathrm{\χ}}_{\mathrm{mag}}}{\arg \min }{\left|\right|\mathrm{\χ}\left(k\right)\mathrm{gH}-\mathrm{diag}\left(H\right)\mathrm{gF\χ}\left|\right|}_2^2+\mathrm{\λ}{\left|\right|P_{\mathrm{mag}}g\▿\mathrm{\χ}\left|\right|}_1+\mathrm{\δ}{\left|\right|W_{\mathrm{mag}}\mathrm{g\χ}\left|\right|}_2$

包含具有压缩感知特性的保真项和具有稀疏特性的幅值先验正则化约束项||P_magg▽χ||₁，并且通过添加感兴趣区域从幅值图像中提取二值加权矩阵W_mag，对原始磁化率分布进行二值加权，减小计算量的同时可以改善磁化率分布重建效果；

2)基于磁化率分布结构特征的磁化率分布图复合重建(S-Step)，重建模型为

包含最小二乘保真项由基础重建(M-Step)获得的磁化率分布图构造的结构特征正则化约束项以及改善重建磁化率分布效果的平滑项磁化率结构特征先验定义为3D图像数据在三个方向上的一阶梯度信息；

所述重建模型中，χ表示待求解的空间离散的磁化率分布，λ和δ为可调节的系数，▽为梯度算子，C为卷积核矩阵，表示相位图，F表示离散傅里叶变换；对于l1范数优化问题，采用迭代阈值法(Iterative threshold algorithm，ITA)来进行处理；之后基于l2范数的凸函数性质，采用共轭梯度法进行求解，以保证收敛速度。

2.如权利要求1所述基于结构特征的自适应定量磁化率分布图复合重建的方法，其特征在于在步骤1)中，所述保真项采用区域自适应截值法算出的初始磁化率值的傅里叶空间数据χ(k)，并对其进行三维二值加权处理，舍弃那些会带来重建伪影的较小值，类似于压缩感知理论中的欠采样过程；所述幅值先验正则化约束项代表3D图像数据分别在三个方向的一阶梯度信息；由于图像在全变分(total variation,TV)变换域内是稀疏的，因此设计基于l2范数的全变分平滑约束；约束项并不对所有的值都进行处理，而是保留在幅值图像内边缘梯度较大的值，将其权重设置为0，对于那些幅值梯度较小的区域权重设置为1，进行稀疏域内的l1范数最小化处理，以实现幅值结构特征先验的约束；加权矩阵P_mag的提取目的是使重建磁化率分布图与幅值图在结构上趋于一致，消除磁化率分布图在平滑区域的细小结构伪影，重建出磁化率分布图中的梯度较大的边界信息；二值加权矩阵W_mag将感兴趣区域权重设置为1，不需要计算的区域权重设置为0。

3.如权利要求1所述基于结构特征的自适应定量磁化率分布图复合重建的方法，其特征在于在步骤2)中，所述保真项在变换域区分开，选取时域内的保真项，以使整个求解过程始终在可行域内。

4.如权利要求1所述基于结构特征的自适应定量磁化率分布图复合重建的方法，其特征在于在步骤2)中，所述约束项采用三维一阶梯度差分算子对原始磁化率进行处理，代表着图像三个方向的一阶梯度信息，相比于基础重建中所使用的幅值先验加权，复合重建中所使用的新的从初步重建磁化率分布图中提取的先验加权矩阵能够为复合重建提供更加准确的先验结构信息；同时为改善重建磁化率分布效果，引入从初步重建的磁化率分布图中提取的新的二值加权矩阵其获得方式与初步重建时一致，只是对于获取先验信息的对象由幅值图变为基于幅值先验重建出的磁化率分布图。