[0001] 本发明受美国政府支持及美国国立卫生研究院（NIH）颁发的R43NS076092号研究经费支持。美国政府对本发明享有特定权利。本声明仅为符合美国37 C.P.R.§401.14(f)(4)条款，且不应视为声明或承认本申请对特定发明的公开和/或权利要求。

[0002] 本发明主要涉及提供磁化率定量成像（磁量图）的系统和方法。特别地，本发明是系统、流程及计算机可访问介质的示例性实施例，旨在从磁共振成像（MRI）信号中提取以磁化率为例的物质的固有物理属性。

[0003] 磁化率是物质固有的物理性质，决定着物质置于外部磁场（如MRI扫描仪中的主磁场B0）中时物质内部经诱导产生的磁化强度。通过磁共振成像（MRI）得到磁化率的量化图像，为测量白质束中的髓鞘（见Li et al. Neuroimage 2011;55,pp.1645-1656;Liu.Magn Reson Med 2010;63,pp.1471-1477;Wu et al. Magn Reson Med 2012;67,pp.137-147)）、大脑皮层灰质（见Shmueli et al. Magn Reson Med 2009;62,pp.1510-1522)）、脱氧或退化血液中的铁及非血红素铁的沉积（见Haacke et al. J Magn ResonImaging2010;32,pp.663-676;Liu et al. Magn Reson Med 2011;66,pp.777-783;Liu
et al.Magn Reson Med 2009;61,pp.196-204;Schweser et al. Neuroimage 20
11;54,pp.2789-2807;Wharton et al. Magn Reson Med 2010;63,pp.1292-1304;Wharton and Bowtell.Neuroimage 2010;53,pp.515-525）以及临床和预临床检验中造影剂的生物分布（见Liu.Magn Reson Med 2010;63,pp.1471-1477;de Rochefort et al. Magn Reson Med 2008;60,pp.1003-1009;de Rochefort et al. Med Phys 2008;35,pp.5328-5339;Liu et al. Magn ResonImaging 2010;28,pp.1383-1389）提供了一种非介入式的手段。正因为它有如此广阔的潜力，磁量图（QSM）已经受到越来越多的科研和临床的关注。
组织磁化率的MRI信号模型

[0004] 磁场通过控制自旋进动过程来决定MRI信号的相位，这使得对组织局部磁场进行估计成为可能（见Li.Magn Reson Med 2001;46,pp.907-916;Li and Leigh.Magnetic Resonance in Medicine 2004;51,pp.1 077-1 082;Li.American Journal of Physics2002;70,pp.1029-1033;Li and Leigh.J Magn Reson 2001;148,pp.442-448）。组织局部磁场强度平行于主磁场强度B0，它相对于B0在图像域（称为r空间）的比值bL可以看作是
2 3
偶极子核函数d(r)=(1/4π)(3cosθ-1)/|r| 与组织磁化率的三维分布χ(r)的卷积，再加上观察噪声n(r):

[0005] 这种组织磁化率的信号模型也可以在傅里叶双重k空间中用矩阵的形式进行表示：Bl(k)=D(k)X(k)+N(k) [0.2]
2 2
这里Bl＝FbL其中F为傅里叶变换算子，X=Fχ，对角矩阵D＝Fd=1/3-kz/k,以
及 N ＝ Fn（见 Koch et al. Phys Med Biol 2006;51,pp.638 1-6402;Marques and Bowtell.Concepts in Magnetic Resonance Part B-Magnetic Resonance Engineering
2005;25B,pp.65-78;Salomir et al. Concepts in Magnetic Resonance Part B-Magnetic Resonance Engineering 2003;19B,pp.26-34）。

[0006] 偶极子核具D=0的一个非平凡零空间，它由一对与B0方向成54.7°度角的对立的圆锥面组成。此零空间使得偶极子核对磁化率进行降采样或部分采样，并且使得磁化率发生卷折。因此，从局部磁场定量地推导磁化率是一个病态的逆问题，从而需要额外信息来唯一确定磁化率。磁化率定量成像（QSM），正如方程0.1或方程0.2所示，其数学本质可以等同于在部分采样的、卷绕的以及模糊的、含有噪声的图像数据基础上进行超分辨率图像重建（见Park et al. IEEE Signal Processing Magazine 2003;20,pp.16;Farsiu et al. IEEE transactions on image processing:a publication of the IEEE Signal Processing Society 2004;13,pp.1327-1344）。利用先验知识的正则化方法，尤其是贝叶斯随机的方法，可以用来对磁化率进行唯一地估计。然而，任何正则化方法的使用都需要严谨的证实和评价。

[0007] 总的来说，从MRI信号相位估计得到的局部磁场中的噪声具有复杂的概率分布，只有在 MRI 信号信噪比（SNR）远大于1的时候可以近似为高斯分布（见Gudbjartsson and Patz.Magn Reson Med 1995;34,pp.910-914）。高斯分布近似由于其简单性而在文献中被广泛使用。然而，在某些高磁化率区域，信噪比较低使得高斯分布近似不再成立。高斯分布近似中相位噪声的方差等于信噪比倒数的平方，而且在MRI中随空间变化。当使用贝叶斯方法对数据保真项进行处理时需要考虑这种随空间变化的噪声方差。磁化率定量成像（QSM）的先前方法

[0008] 从磁场测量中确定磁化率的逆问题的数学本质，是经典的超分辨率（SR）图像重建问题 （见Park et al. IEEE Signal Processing Magazine 2003;20,pp.16）。已知的超分辨率问题方法可用于解决QSM重建问题，其中包括通过最大化后验项，即最大化后验概率分布的贝叶斯统计方法（见Park et al. IEEE Signal Processing Magazine 2003;20,pp.16;Elad and Feuer.IEEE transactions on image processing:a publication of the IEEE Signal Processing Society 1997;6,pp.1646-1658）。这就给出了求解磁化率的一般方法：最小化代价方程。这个代价方程由数据保真项L[χ(r)]（与似然函数有关）和正则化项R[χ(r)]（与先验概率分布有关）组成：
χ(r)=argminχ(r)L[χ(r)]+αR[χ(r)], [0.3]
这里 其中w=SNR为了考虑MRI信号相位高斯分布
中的噪声方差随空间的变化；正则化项通常表示成某个惩罚函数C的Lp的p次幂形式，此惩罚函数C需指数地满足 解剖结构先验信息可以以掩膜的形式（见
Gindi et al. Ieee Transactions on Medical Imaging 1993;12,pp.670-680）或梯度权重的形式（见Elad and Feuer.IEEE transactions on image processing:a publication of the IEEE Signal Processing Society 1997;6,pp.1646-1658;Baillet and Garnero.IEEE transactions on bio-medical engineering 1997;44,pp.374-385）引入到正则化项中。先前的QSM工作致力于正则化项R的构建。人们尝试了各种正则化方法，包括L1范数及梯度方程G的L2范数的平方（见Kressler et al. Ieee Transactions on Medical Imaging 2010;29,pp.273-281）以及小波分解的L1范数（见Wu et al. Magn Reson Med
2012;67,pp.137-147）。引入形态学的偶极子求逆（MEDI）过程是一种示例性过程（见Liuet al. Neuroimage 2012;59,pp.2560-2568），它在正则化项中引入一个从解剖结构先验a(r)（如幅度图像）的梯度中获得的边缘掩膜M：当|Ga(r)|小于阈值时M(r)=1，否则M＝0。经过向量化后的MEDI正则化项为：
R=||MGχ||1. [0.4]

[0009] 另外还有非贝叶斯的方法，包括截断k空间分解方法，其中当偶极子核绝对值小于某阈值时可以用该阈值代替偶极子核（见Shmueli et al. Magn Reson Med2009;62,pp.1510-1522）以及加权k空间导数（WKD）方法，其中当偶极子核绝对值小于某阈值时使用偶极子核的线性内插代替偶极子核(见Li et al. Neuroimage
2011;55,pp.1645-1656)。这些分解方法在偶极子核值较小的区域，噪声会被严重地放大。
先前方法的示例性局限

[0010] 本发明介绍了先前方法中一些主要的示例性局限。一个局限是在方程0.4中，解剖结构先验信息以原始边缘掩膜的形式引入，这使得重建后的磁化率图在L1范数上体现为块状或像素状。这是由于细小结构因对比度噪声比（CNR）过低而没有包含在方程0.4中。

[0011] 另一个示例性局限是在高磁化率区域，由于磁化率非均一性引起严重失相，导致信号丢失，从而使得信噪比过低。此时，MRI信号相位的高斯模型变得十分不准确，使得方程0.3中数据保真项导致在重建的磁化率图中产生严重的伪影。

[0012] 因此，有必要根据本发明中的示例性实施例，使用系统、流程以及计算机可访问介质的示例性实施例，对其中一些示例性缺陷和/或局限进行阐述。

[0013] 根据本发明中使用的示例性实施例提供的针对MRI信号采集和处理的方法，系统和计算机可访问介质，以磁化率图像为例，目的在于将物质固有物理属性进行重建显示。

[0014] 本发明示例性实施例的目标之一是克服解剖结构先验信息使用中的示例性局限。这提供了示例性的正则化方法，系统及计算机可访问介质，用来选择合适的方法求解符合已知解剖图像结构信息的磁化率图，包括在磁化率χ(r)和已知解剖信息间极大化相关度，极大化互信息以及极大化联合熵的方法。

[0015] 本发明示例性实施例的另一个目标是提供系统和计算机可访问介质来克服传统磁化率定量成像（QSM）技术中，局部磁场噪声的高斯模型的误差。某些实施例通过如在复数MRI信号域改写数据保真项的方式实现了这一目标，原因是复数域中噪声的确满足高斯分布。一种方式是将线性信号方程0.1改写为复数信号方程，并且将方程0.3中二次数据保真项替换为非线性复数项。

[0016] 本发明示例性实施例进一步的目标在于，将磁化率定量成像用于解决临床问题。示例性的实施例包括改善脑结构显示从而研究脑深部电刺激，改善多发性硬化、出血及钙化中病变的评估以及改善造影剂增强MRI中，造影剂生物分布的评价。

[0017] 本发明示例性实施例包括系统，流程及计算机可访问介质，从而确定与至少一个结构相关的信息。这些示例性实施例包括至少一个的计算环境。这种计算环境可以接收与至少一个结构相关的信号，而这种信号包括复数域中的复数数据。这些示例性实施例也可以在真实复数信息的基础上，产生针对至少一个结构的磁化率图像，同时至少一部分真实信息和至少一个结构的噪声属性有联系。

[0018] 某些示例性实施例通过贝叶斯方法产生磁化率图像。这种示例性方法包括通过利用一个含有与局部磁场有关的相位因子的指数方程和一个含有与偶极子核卷积后磁化率有关的相位因子的指数方程，来构建高斯似然项。这种示例性方法还包括针对至少一个结构数据产生先验信息。再有，这种示例性方法也包括对由似然函数和先验信息决定的后验概率项的极大化。

[0019] 在某些示例性实施例中，通过引入复数数据的非线性最小二乘拟合，可以用复数数据确定局部磁场。某些示例性实施例利用拉普拉斯算子，隐式消除因目标区域外磁场源引起的背景场，从而确定局部磁场。在某些示例性实施例中，拉普拉斯算子可以通过狄拉克δ（Dirac delta）函数与球形均值（spherical mean value）函数求差来实现。在某些示例性实施例中，后验概率可以通过对与其相关的代价方程迭代线性化来进行极大化。

[0020] 某些示例性实施例中，后验概率可以通过似然函数中变化的加权过程进行迭代极大化。某些示例性实施例中，初始权重由与邻域相关的校正噪声方差提供。某些示例性实施例中，先验信息可以整理为L1范数，与一种结构相关的结构信息从该结构图像中估计获得。某些示例性实施例中，结构信息可以是图像的梯度。某些示例性实施例可以通过引入空间解卷折过程来确定局部磁场。

[0021] 某些示例性实施例可以将磁化率图像同解剖结构图融合显示。某些示例性实施例可以利用生成的磁化率图显示丘脑底核，或观察多发性硬化病变中一种或多种、癌化病变、钙化病变或者出血。某些示例性实施例可以利用生成的磁化率图对造影剂增强磁共振成像中造影剂生物分布进行量化。

[0022] 另外一种示例性实施例，为了从磁共振信号中确定至少与一个结构有关的信息，提供系统、流程与计算机可访问介质。它们可以接收与至少一个结构有关的信号，并且在与至少一种结构有关数据的基础上，利用信息论过程，生成至少一个结构的物理属性的图像。

[0023] 某些示例性实施例中，这种物理性质可以包括与至少一个结构有关的磁化率，并且这些示例性实施例可以用贝叶斯方法生成磁化率图像。一种示例性贝叶斯方法包括利用局部磁场和偶极子核卷积的磁化率构建似然函数，利用至少一个结构的信息论量度及与至少一个结构有关的磁化率图生成先验概率，和/或极大化似然函数和先验概率决定的后验概率。

[0024] 某些示例性实施例通过引入非线性最小二乘的拟合过程，从与至少一个结构有关的信号确定磁场。某些示例性实施例利用拉普拉斯算子，隐式消除因目标区域外磁场源引起的背景场，从而确定局部磁场。在某些示例性实施例中，拉普拉斯算子可以通过Dirac delta函数与spherical mean value函数求差来实现。在某些示例性实施例中，似然函数可由局部磁场与偶极子核卷积的磁化率的差值的加权L2范数的平方构成。

[0025] 在某些示例性实施例中，似然函数可由一个含有与局部磁场相关的相位因子的指数方程与一个含有与偶极子核卷积的磁化率相关的相位因子的指数方程的差值的L2范数的平方构成。在某些示例性实施例中，信息论量度可以是与至少一个结构有关的至少两幅图像的联合熵。在某些示例性实施例中，信息论度量可以是与至少一个结构有关的至少两幅图像的互信息。

[0026] 某些示例性实施例可以通过线性化同后验概率有关的代价方程，迭代地极大化后验概率。某些示例性实施例可以利用似然函数中随迭代过程变化的加权过程，迭代地极大化后验概率。在某些示例性实施例中，初始权重由与邻域相关的校正噪声方差提供。某些示例性实施例可以将磁化率分布同解剖结构图融合显示。某些示例性实施例可以利用生成的磁化率图显示丘脑底核，或观察多发性硬化病变中一种或多种、癌化病变、钙化病变或者出血。某些示例性实施例可以利用生成的磁化率图对造影剂增强磁共振成像中造影剂生物分布进行量化。某些示例性实施例可以通过引入空间解卷折过程来确定局部磁场。

[0027] 本发明中示例性实施例的以上及其它目标、特性以及优点，结合所附权利要求，将在阅读关于本发明中示例性实施例的如下详细说明之后变得显而易见。附图说明

[0028] 本发明的进一步目标、特性和优点，从如下详细描述并结合阐述本发明实施例的附图中，将变得显而易见，其中：

[0029] 图1为磁化率定量成像生成过程的功能性框图。此过程可以用复数三维磁共振图像作为输出，输出三维磁化率图像，这将由本发明中示例性实施例系统来完成；

[0030] 图2为描述图1中场图估计过程中示例性操作的流程图；

[0031] 图3为描述图1中场图求逆过程中示例性操作的流程图；

[0032] 图4a-4c为一幅3x3图像的相位图的示例，其中被测量噪声所影响。这包括箭头所示的旋涡状噪声，而灰色框中显示对应的相位解卷折误差；

[0033] 图5a-5c横轴方向示例：解剖结构先验信息，未正则化的磁化率的解以及利用示例性联合熵正则化过程后显著改善的磁化率的解；

[0034] 图5d-5f矢状方向图像：解剖结构先验信息，未正则化的磁化率的解以及利用图5a-5c中示例性联合熵正则化过程后显著改善的磁化率的解；

[0035] 图6a和6b为线性相位拟合中，测量相位中噪声（如圆圈）引起拟合相位（如十字）对真实相位（实线）的显著偏离的示例性示例图。而这种噪声干扰现象会在非线性拟合中的复平面上显著减弱；

[0036] 图6c和6d为1000次独立生成噪声的重复实验中，估计场与真实场的差值示意图。图中显示线性拟合方法相对非线性拟合有着更大的差值；

[0037] 图7a-7e为示例性幅度图像，真实磁化率分布，局部场，线性QSM重建与金标准的差值图（出血区域误差较大）以及非线性QSM重建与金标准的差值图（出血区域误差显著减小）；

[0038] 图8a-8e为矢状方向示意图：示例性幅度图像，开放边缘的相位图（黑色箭头），包含对应解卷折误差的频率图（黑色箭头），线性方法重建的QSM以及非线性方法重建的QSM。其中线性方法结果中有严重的条状伪影，而非线性结果中没有条状伪影；

[0039] 图8f-8j为图8a-8e对应图像的轴向视角显示（图8d中用白色虚线标明的轴向区域），分别为：示例性幅度图像，开放边缘的相位图（黑色箭头），包含对应解卷折误差的频率图（黑色箭头），线性方法重建的QSM以及非线性方法重建的QSM。其中线性方法结果中有严重的条状伪影，而非线性结果中没有条状伪影；

[0040] 图9a-9e为示例性幅度图像，线性QSM，非线性QSM，放大后的线性QSM（虚线框）以及放大后的非线性QSM。箭头标示处显示非线性QSM对线性QSM的细微改善；

[0041] 图10a-10e通过第一个颅内出血的例子，显示了：出血幅度图，含有影状及噪声状伪影的线性QSM结果，含有显著减少的影状伪影的MERIT线性QSM结果，含有影状伪影但不含噪声伪影的非线性非MERIT的QSM结果以及同时含有显著减少的噪声及影状伪影的非线性MERIT的QSM结果；

[0042] 图10f-10j通过第二个颅内出血的例子，显示了：出血幅度图，含有影状及噪声状伪影的线性QSM结果，含有显著减少的影状伪影的MERIT线性QSM结果，含有影状伪影但不含噪声伪影的非线性非MERIT的QSM结果以及同时含有显著减少的噪声及影状伪影的非线性MERIT的QSM结果；

[0043] 图10k-10o通过第三个颅内出血的例子，显示了：出血幅度图，含有影状及噪声状伪影的线性QSM结果，含有显著减少的影状伪影的MERIT线性QSM结果，含有影状伪影但不含噪声伪影的非线性非MERIT的QSM结果以及同时含有显著减少的噪声及影状伪影的非线性MERIT的QSM结果；

[0044] 图11a-11e为脑深部核团（黑箭头）与白质纤维束（白箭头）的轴位与矢状位示意图；

[0045] 图12a-12e为丘脑底核在QSM图，磁化率加权图，幅度图，R2*图以及T2加权图的冠状位示意图。其中QSM提供了丘脑底核（箭头）和周围组织之间的最佳对比度；

[0046] 图13a-13d为白质病变在QSM图，T1加权图，T2 FLAIR图以及T2加权图的示意图。这为病变分级与疾病分阶段提供了机会；

[0047] 图14a和14b为脑部肿瘤在幅度图和QSM图上的示例图；

[0048] 图15a-15d为脑部肿瘤中造影剂在QSM（顶部）和T1加权幅度图（底部）上的标示示意图；

[0049] 图16为钙化病变的QSM，CT显示图以及QSM与CT相关度的示意图；

[0050] 图17为将QSM与解剖结构图融合后，利用灰度及彩色方案显示的示例图；

[0051] 图18为本发明中示例性实施例系系的框图；

[0052] 在整个附图中，相同的参考数字和字符，除另有注明外，都用于表示一致的特征、元件、组件以及图示实施例的部分。此外，尽管本发明将参照附图来进行详细的介绍，但这仅供说明性地阐述实施例且不限于附图中或所附权利要求中阐述的特定的实施例；

[0053] 应该指出，这里的示例性实施例的附图与说明已被简化，着重体现了主要元素以更好地理解整个示例性实施例，同时为了清楚起见，删除了其他众所周知的元素。接下来将参照附图进行详细的介绍。

[0054] 本发明中的示例性实施例可以克服同多种可能存在的误差源有关的缺陷。为了克服先前方法中出现的某些示例性局限，可以施加多种改善正则化项的设置，而且人们也提出了其他改善数据保真项的示例性设置方式。于是，以下将介绍这两种示例性方法，系统及计算机可访问介质的具体实施例，并附有实验结果来阐述其改进点。1.MRI信号方程在磁化率测定中的误差源分析
MRI相位数据中的误差源

[0055] 在QSM中，描述磁场bL变化的前向问题可以定义为磁化率分布χ与单位偶极子d磁场（以下称为偶极子核）的卷积。这个卷积（如方程0.1）可以重新表示为偶极子核矩阵D作用在磁化率分布向量上：bL=Dχ+n. [1.1]

[0056] 在MRI中，磁场bL等于失谐频率图f与旋磁比 倒数的乘积。频率图f可从涉及时间和/或空间解卷折的含有噪声的磁共振（MR）相位图得到。这种频率估计可能含有两种误差源：噪声与解卷折误差。这些误差可能传播至计算出的磁化率分布之中。

[0057] 考虑到MR相位数据中的噪声，观察到的MR信号可被建模为2
这里A和 分别表示没有噪声时信号幅度和相位，i=-1，ε表示0均值的复高斯噪声，其实部与虚部在图像域中以常方差σ2独立同分布。相位图可以逐像素地求解计算S的实部与虚部的反正切来计算。相位噪声的分布， 可能如(见Gudbjartsson and Patz.Magn Reson Med 1995;34,pp.910-914)给出的那样具有复杂的表达式，但在A远大于σ时
2 2
可以近似为均值为0方差为(σ/A)=SNR 的高斯分布，或当A=0可以近似为(-π，π)间均匀分布的随机变量。但A≈σ时没有简单的近似模式。

[0058] 相位噪声会进一步引起解卷折误差。在解卷折中，当邻近体素（空间或时间上）相位差值小于π时，每个体素中的测量相位值被加上或减去2π的整数倍。然而，噪声会导致解卷折单个体素时错误地加上或减去2π的整数倍，进而导致顺序相位解卷折算法中后方体素出现错误(见Goldstein et al. Radio Science 1988;23,pp.713-720).再有，噪声会在三维相位图中产生“漩涡”，漩涡出现处每四个相邻点的卷绕后的相位差值顺时针求和不等于0。这种情况下，每一对相邻体素之间的空间连续性无法得到保证，使得出现接近2πk的相位跳变，k为整数。图4说明了这个问题，其中数值用π进行了归一化。图4(a)表示真实相位，图4(b)表示含噪声的测量相位，图4(c)表示解卷折后的相位。由于所示漩涡的出现（如箭头410），可能无法消除解卷折相位中大于π的间断点，如420框所示。
从相位估计f及偶极子卷积中的模型误差

[0059] MRI中测量到体素中相位是该体素复数信号的相位∠∑sμsexp(-iγbsTE+φs)，这里∠表示复数的相位，∑s为该体素中所有自旋等色原子求和，μs表示磁化强度，bs为自旋等色原子s受到的局部磁场（如沿B0方向），而φs为该自旋等色原子的初始相位。体素中单个自旋等色原子的信号相位可能不同于该体素中自旋等色原子信号的平均相位，∑s(γbsTE+φs)/∑s1。当bs与φs在体素内相对均一时，前者是后者的一个良好近似。当此条件不满足时，这种近似将不再成立。后者反映了体素内平均磁场，是方程1.1中求解磁化率图所关心的量。其他来源，如运动，也会引发相位的变化。这些未考虑的相位影响也会引起模型误差。2.解剖结构先验信息的组成
I.对磁化率图与先验图像间边缘匹配整合方式的改进

[0060] 解剖结构的边缘匹配可通过多种边缘方程形式，整合到正则化项之中。它们包括：
a.梯度(G)的L2范数中体现的空间平滑性(GL2)(Tikhonov正则化)(见de Rochefort et al. Magn Reson Med 2010;63,pp.194-206;Kressler et al. IEEE Trans Med Imaging
2010;29,pp.273-281)
b.梯度L1范数中体现的稀疏性(GL1)(见Kressler et al. IEEE Trans Med Imaging
2010;29,pp.273-281;Liu et al. Neuroimage 2012)
c.分数阶Lp
其中向量v的Lp范数定义为
d.总变差（Total variation）范数中体现的稀疏性(TV)(见Rudin et al.Physica D
1 992;60,pp.259-268)
e.小波域（如德比契斯（Daubechies）小波）中体现的稀疏性(Φ)(见Wu et al. Magn Reson Med 2012;67,pp.137-147)
R=||Φχ||1;
f.两种稀疏项的组合，如total variation和小波(见Wu et al. Magn Reson Med
2012;67,pp.137-147)
αR→α|Φχ|1+βTV(χ)
αR→α|Gχ|1+β||Gχ|1

[0061] 由于在一种应用中，这些通用数学方程不总含有特定解剖结构信息，所以这些方程需要依据解剖先验信息进行调整，例如边缘掩膜。边缘掩膜可以对任意反映解剖结构的图像作canny边缘检测来生成。这些图像包括幅度图，相位图，T1加权图及T2加权图。不同图像得到的边缘掩膜可以利用联合处理整合为一个复合边缘掩膜。

[0062] 小波变化可以看做是一种边缘检测处理，例如Daubechies小波和哈尔（Haar）小波。对每个小波组分，解剖先验a(r)的二值掩膜可以从先验图像的对应小波变换生成。这种方式同将解剖先验掩膜引入梯度项的方法类似，且惩罚函数可以依据此先验掩膜进行加权：M(r)=1如果|Φa(r)|<阈值;否则=0.
R=||MΦχ||1;
II.利用信息论量度的解剖先验整合

[0063] 先前研究解剖先验信息使用的工作显式地从复数MR数据中推导边缘图像，通常涉及对图像求梯度并确定强梯度所在位置。于是，在所有候选解中，边界与解剖先验边缘的差别最稀疏的解被选出来。

[0064] 本发明的示例性实施例利用信息论知识，为QSM提供了一个引入解剖先验信息的过程。某些示例性方法不需要显式地计算边缘或者梯度，而是在信息论量度下，让最终的QSM解包括候选解中与解剖先验相似度最高的一个解：*
χ=argminχIT(X，Y),subject toL(χ)＝μ [2.1]

[0065] 信息论量度IT(X，Y)可以表示候选解X与先验Y的联合熵或者负的互信息。根据本发明中某些示例性实施例，可以用pX，Y(x，y)表示X和Y的联合概率密度函数，而pX(x)与pY(y)分别表示X和Y的边缘概率密度函数。因此联合熵定义如下：，而负的互信息表示为：

[0066] 约束条件L(χ)＝μ，或者等效地,||W(Dχ-bL)||22＝μ,可以保证所选的磁化率分布满足数据保真性，即该解生成的场与MR数据中测量的场相同，且二者差别仅处于噪声水平。

[0067] 由于对利用任何基于梯度下降的求解器对极小化问题进行数值求解时，正则化方程的导数信息十分重要，所以我们在此指出，联合熵的导数可以有效地近似计算。由于二者的相似性，同样的近似可引申至互信息情况。

[0068] 子过程1.联合概率密度函数的估计

[0069] 尽管图像X与Y的联合直方图可以直接从二者信号强度χ及a直接计算，作为联合概率密度函数，可以是χ和a的不连续方程，这为下面的求导计算提供了困难。所以，我们用Parzen窗的方法对联合概率密度方程进行估计：这里N为体素数，φ(x/σ)表示均值为0标准差为σ的高斯窗，这里σ为设计参数。

[0070] 子过程2.函数导数计算

[0071] 正则化函数对χ的导数等于：

[0072] 这里有两种计算 的示例性方法，称为A)内插近似和B)单高斯近似A.内插近似

[0073] 梯度向量 第k个分量可以从方程[2.2&2.5]利用数值近似进行计算：这里B是概率密度函数总计算点的数目，Δx和Δy分别等于X和Y的区间大小除以B，φ′[(xi-χk/σx]＝φ[(xi-χk/σx][(xi-χk/σx2].

[0074] 某些示例性实施例中，如果χ和a中所有体素强度均分布在B个等间距网格上，联合熵及联合概率密度函数的计算实际上将以卷积的形式出现，从而可以用快速傅里叶变换（FFT）有效地计算。所以，计算 可以先将χ和a的强度离散化为 和 接着估计 最后通过双线性内插得到
B.单高斯近似

[0075] 单个二维高斯近似可用来近似联合概率密度函数。特别的，方程2.5中N个二元高斯项可以近似为具有相同的一阶矩及二阶矩的一个二元高斯项：这里 而 拟合后单个高斯项的熵为：

[0076] 于是，为了极小化惩罚函数H(X，Y)，根据本发明的示例性实施例，可以简单地极小* *化|∑ |。利用变量代换，|∑ |简化为：
* 2 2 2 2
R（χ)=| ∑ |=[σ+E(a·a)-E(a)][E(χ·χ)-E(χ)]-[E(χ·a)-E(χ)E(a)],
[2.10]这里E(a)表示图a中体素亮度均值，符号·表示逐点相乘，且：
这里u表示元素全为1的向量。
3.精确磁化率定量成像（QSM）中数据保真的非线性示例性过程
非线性频率图估计

[0077] 在梯度回波MRI中，由于场不均匀性的原因，相位可能随回波时间而改变。尽管相位噪声有着复杂的分布，测量的复数MR信号中噪声实部与虚部实际上满足正态分布或高斯分布。于是，从多TE测量的复数MR信号中估计失谐频率f的问题可被整理为一个非线性最小二乘拟合问题，* * 2
fr,φ0r=argminfr,φ0rΣj||Sr(TEj)-Ar(TEj)exp[i(fr×TEj+φ0r)]||2, [3.1]这里下标r表示空间位置，下标j表示第j个回波。真实复数信号的相位建模为
fr×TE+φ0r其中f为场不均匀性引起的失相频率，φ0r为初始相位，它可能由激发频率与局部拉莫（Larmor）频率不匹配产生。
针对偶极子求逆的非线性数据保真项

[0078] 从估计磁场计算磁化率是病态的逆问题，因为多个候选磁化率解可以生成同样的磁场。先前数据保真项表示成线性最小二乘拟合的形式。本发明的示例性实施例将数据保真约束表示为关于场和磁化率的复指数方程形式，* 2
χ=argminχR(χ),subject to||W(exp(iDχ)-exp(ibL))||2 ＝μ. [3.2]

[0079] 这里，W可以是反映每个体素中估计频率可靠性的对角权重矩阵。bL可从b中减去背景场bB估计得到。μ为期待的噪声水平。复指数方程对磁化率及场来讲可以是高度非线性的，进而准确的分析MR信号中的噪声并且不受错误的2π相位跳变的影响。

[0080] 此非线性数据保真函数的导数可以通过如下方式计算。首先，数据保真项2 (n)
L(χ)=||Wexp(iDχ)-Wexp(ibL)||2 可以在某个解χ 处用一阶Taylor展开近似为二次函数，
(n) 2
L(χ)=||Wexp(iDχ )(1+iDχ)-Wexp(ibL)||2. [3.3]
(n)
使用变量代换w=Wexp(iDχ ),b=Wexp(ibL),梯度
4.示例性QSM过程的实现

[0081] 图1为从MR相位数据计算磁化率一个示例性过程的流程图。首先，将一个或多个复数图像输入至100，接着进行如下示例性步骤。A.多通道数据整合102

[0082] 子过程102只对多通道数据有效。当在第j个回波时对自旋分布ρj进行成像时，第i个通道采集的信号si，j会被该通道的线圈敏感度ci调制，这里ρ，si，j及ci均由复数组成。线圈敏感度已知时，自选分布可以通过ρj(r)=(cHc)-1cHsj估计得到，这里c和sj为所有通道的ci(r)及si，j(r)的向量表示。当线圈敏感度ci未知时，第k个回波的复数图像si,k可以充当线圈敏感度。这种方法会导致ρk仅具有0相位，但不会改变ρ1,ρ2，…,ρn中相位的改变速率，其中包括了接下来将要估计的目标量。B场图估计110

[0083] 子过程110只对多回波数据有效，否则场图被认为与相位图像成正比，比例因子为回波时间的倒数。图2阐释一个示例性子过程，包括步骤110的示例性实现步骤。

[0084] 例如，在200中，示例性过程接受单通道复数图像。210中f与θ0的非线性估计可以用迭代地方式得到。再例如，此过程设第n次迭代中场及初始相位分别为fn和φ0n(简单起见，空间下标r可省略).方程3.1可重写为：Σj|S(TEj)-A(TEj)exp[i(fn+f-fn)×TEj+(φ0n+φ0-φ0n))]|2,或
Σj|S(TEj)-A(TEj)exp[i(fn×TEj+φ0n)]exp[i(dfn×TEj+dφ0n)]|2,
这里dfn=f-fn，dφ0n=φ0-φ0n两处更新是为了极小化以上方程。可以使用小幅更新，并且以上方程可以用一阶泰勒（Taylor）展开进行线性化：
Σj|S(TEj)-A(TEj)exp[i(fn×TEj+φ0n)]-A(TEj)exp[i(fn×TEj+φ0n)]
i(dfn×TEj+dφ0n)|2,
得到关于dfn和dφ0n的二次方程，可用线性最小二乘拟合有效地求解。

[0085] 鉴于真实信号幅度A未知，s的幅度可以作为A的近似。示例性的初始f0及φ00可以从前2个回波使用线性拟合进行估计，其中TE=0时的相位设为φ00。当[dfn dφ0n]T的模小于[f0 φ00]T模的10-4，或者例如n大于10时，迭代结束。[dfn dφ0n]的噪声相关矩阵可以在线性最小二乘拟合中估算，且用dfn方差的平方根的逆构成例如220中权重矩阵W中的权重。感兴趣区域（ROI）之外体素的权重可设为0。余下非零权重可进一步用他们的均值进行归一化。

[0086] 为了解决频率图中的频率混叠问题，两种独立的方法可用来进行空间解卷折（如230中）。一种示例性实施例使用了一种信号质量图引导的空间解卷折算法，其中解卷折路径由频率图的估计噪声水平及空间连接性决定，所以噪声较少的区域优先被解卷折。在这种方法中，从预定义的ROI质量中心开始使用一种区域生长方法。在解卷折后区域边缘的体素会被扩张，而其中具有最低噪声水平的体素会被解卷折。这种示例性过程会被重复直到ROI中所有体素均被访问。所有边缘体素的噪声水平可以存储在一个堆数据结构中，以便对最低噪声体素的快速定位。

[0087] 外一种示例性方法应用了基于拉普拉斯算子的方法。它首先计算了解卷折的频率图的拉普拉斯：这里 表示未卷绕的频率图，而 表示卷绕的频率图，接着对 使用逆拉普拉斯
算符来恢复频率图

[0088] 最后，解卷折后的频率图可以通过除以如240中的中心频率来得到相对差值场图。C.偶极子反演

[0089] 图3中显示了在计算局部磁化率分布时，为了消除背景场可以使用的两个示例性、独立的的方案。这两个方案是a)执行明确的背景场移除和b)在偶极子核中隐含嵌入背景场移除，其实现参见300和310子程序。

[0090] 明确的背景场移除，偶极子域投影（projection onto dipole fields）(PDF)算法(见Liu et al. NMR Biomed 2011;24,pp.1 129-1136)或者相位数据精密谐波伪影缩减（sophisticated harmonic artifact reduction for phase data）(SHARP)算法可以使用在320中。PDF把背景场和由ROI中的磁源产生的局部磁场分离。在人脑成像中ROI可以是脑实质，并且可以使用自动脑区域剪裁算法来得到(见Smith.Hum Brain Mapp2002;17,pp.143-155)。在PDF中，所有脑区域以外且在视场中的体素都可以被假定对脑区域中的背景场有贡献。体素的背景场的强度可以通过对脑区域中的磁场的加权最小二乘法拟合来得到。校正以后的磁场作为场到源逆问题的输入。SHARP算法利用了背景场在大脑中是一个谐波方程的事实。SHARP算法首先对由背景场和局部场构成的总场b进行一个拉普拉斯卷积操作来移除背景场；然后在重置ROI的边界条件为零后，使用反拉普拉斯操作来恢复改动后的局部场。

[0091] 在322中的隐含背景场移除中，场图b、偶极子核d和数据加权W都可以被修正。首先在b上使用拉普拉斯操作来移除所有背景场的贡献：
可以进一步表示为：
这种方法的优点在于对b使用了三角函数，从而对相位卷折不敏感，进而防止了相位反卷折不完全造成的错误扩散。

[0092] 为了在数据保真项中保持一致，需要对可能的解进行同样的卷积运算，即：在本发明中某些示例性实施例，使用标记 可以使用FFT快速有效的
来实现。

[0093] 在方程[4.1&3]中实现拉普拉斯操作 的一种方法可以写为( 其中3
δ(r)是Dirac delta函数，s(r)在|r|＜R0时等于1/(4/3πR0)，其他情况等于0，y表示标量场。这个方法也就是求特定的体素y与它周围特定半径内体素的spherical mean value(SMV)的差。卷积可以快速有效的使用FFT来实现。

[0094] 当使用隐含背景场移除时，数据加权W需要随bL噪声特性的改变而改变。因为2
在bL中的噪声方差变为 其中σb 是b中的噪声方
差(或者等价的，W的对角线元素)，σbL组成了W’的对角线元素。在选择隐含背景场移除时，后续文本中的d和W将被修正后的d’和W’代替。

[0095] 在选择了数据保真项和正则化项330以后,方程[2.1&3.2]中的限定性最小化问题可以改写为非限定性拉格朗日（Lagrangian）最小化问题，*
χ=argminχL(χ)+αR(χ) [4.4]

[0096] 340中的初始化Lagrangian乘子α可以按照经验选择。在本发明的一个实施例-3中，α＝10 。当使用非线性正则化项或者数据保真项时，方程4.4中的非线性最小二乘拟合问题可以使用迭代方法来求解，迭代的初始估计可以采用非正则化的解，例如350。

[0097] 在每次迭代中，代价函数Q(χ)=L(χ)+αR(χ)的导数可以使用当前的解来计算360。这个导数进行长度归一化以后的值将确定线性变换解的更新方向。步长将确定沿此方向前进的距离，通常应该满足阿米霍（Armijo）、沃尔夫（Wolfe）或者金斯坦（Goldstein）条件。

[0098] 两种不同的方法可以用来计算解更新362。如果可以有效构造代价函数的海森（Hessian）矩阵，例如当 是χ的线性函数或者使用非线性数据保真项和MEDI正则化项时，可以在导数上使用Hessian矩阵的逆作为方向的更新，而且准牛顿（quasi-Newton）方法可以用来计算这个更新。a.示例的准牛顿固定点（Quasi-Newton Fixed Point）方法。
因为 是一个矩阵， 操作使用共轭梯度法来实现。

[0099] 如果 不是χ的线性函数或者不能有效构造代价函数的Hessian矩阵，例如使用差值逼近的联合熵时，这个更新可以使用带有Armijo线搜索方法的共轭梯度来计算。b.共轭梯度和Armijo线搜索法

[0100] 示例程序364中可以包括迭代调整模型误差缩减（Model Error Reduction Through Iterative Tuning）(MERIT)方法。模型和观测值的一致性可以由模型错误在体素中导致的误差来折中得到，所以在这些体素上的权重在反解过程中应该递减，以减少误差向临近点扩散。因此，可以通过减少权重来满足磁场和磁化率解的全局一致性。具体来说，可以用一个不断更新的权重项WN代替W，从而抑制这些体素。χ*=argminχR(χ),subject to ||WN(exp(iDχ)-exp(ibL))||22＝μ. [4.5][0101] 公式2.1和3.2中的数据权重WN在每次迭代结束时可以进行调整。具体来说，公式3.2中的权重W可以基于逐体素残差ρn=Wn|exp(iDχn)-Wexp(ibL)|来动态减少，其中下标n表示第n次迭代。所有体素中的噪声标准差σn可以估计得到。残差大于6σn的体素v可以认为是损坏的，它的权重重新计算为Wn+1，v=Wn,v/(ρn，v/σn)t，这里的指数项t决定了调谐权重， 例如一个示例性实施例中使用t=2。这个随迭代减少的权重Wn+1接着被用在下一次迭代中。

[0102] 当更新值的L2范数小于当前解L2范数的一定比例，例如1%；或者达到迭代次数时在366中迭代退出。退出时解的残差可以同期望残差μ做比较。如果残差明显大于μ，可以减小正则化参数α来加大数据保真项，反之亦然，例如372。 这个步骤可以一直继续直到残差跟期望噪声水平大致相差50%μ以内。5.实验

[0103] 一些实验被用来验证在某些示例性实施例中的正则化项和非线性数据保真项的有效性。信息论度量项评估

[0104] 在一个示例性实施例中，五个健康志愿者(3男2女，年龄25±8)被召集来在一个3T扫描仪(HDx,GE Healthcare,Waukesha,WI,USA)中进行脑成像。扫描使用一个在所有回波的读出方向（前/后）上带流动补偿的单极多回波梯度回波序列。其他参数包括：重复时间(TR)=35毫秒（ms），从4ms到32ms均匀分布的8个回波时间（TE），带宽300赫兹（Hz）/体素，视场24厘米（cm），层厚1毫米（mm），采集矩阵240×240×150，重建到256×256×150。采集中使用了加速系数为2的并行采集技术(Array Spatial Sensitivity Encoding Technique)。实部和虚部的医学数字影像和通讯格式（DICOM）图像被保存来进行后处理。不使用正则化的重建被用来与使用正则化项的代价函数Q的重建做比较，Q由线性数据保真项和使用高斯近似的联合熵的正则化项组成。线性数据保真项的导数为H
其中A 表示A的共轭转置。
非线性数据保真项评估

[0105] 在这一节中，代价函数Q由非线性数据保真项和MEDI正则化项R=||MGχ|1组成，其中M是一个从幅度图或者局部场bL得到的二值图，M对应明显梯度的地方设置为0，否则为1。当使用弱导数的绝对值时，其中diag(y)表示一个方阵，其对角线元素为y，其
余为0。由于 和 都是Aχ形式，所以L和R的Hessian矩阵可以高效地计算得
到，并且quasi-Newton fixed point方法可以用来计算更新值。

[0106] 在以下实验中，QSM的重建同时使用本发明中提出的非线性QSM方法与先前开发的被称为线性QSM的方法。所有下列回顾性人类研究都经过了我们机构审查委员会的批准。所有算法在一台包含Intel Core i7处理器和6GB内存的Linux工作站上用C实现。所有人类数据的重建时间都被记录下来。
A.相位噪声扩散的减少

[0107] 数值模拟被用来检测本发明提出的方法减少由相位噪声带来误差的有效性。在频率图估计中，实验仿真了在6ms到96ms之间均匀分布的16个不同TE上的复MR信号。信号幅度在TE=0ms的初始值是25任意单位（a.u），并且以T2*=20ms进行指数T2*递减。真实频率被设置为50HZ。方差为1a.u的随机零均值高斯噪声被独立的加在实部和虚部上。估计频率采用了线性和非线性场图估计法。在线性场图估计中，每个体素首先进行以弧度表示的相位上的空域反卷折，然后使用加权最小二乘拟合估计该体素上的频率。以弧度表示的相位上的空域反卷折法为在相邻两个回波的绝对跳变大于或等于π时，对其进行2π补足。估计频率和真实频率的差别被记录下来。这个实验重复了1000次来研究频率估计错误的均值和标准差。在每个估计方法中，学生t检验（student’s t-test）被用来检验每种方法是否有偏差。

[0108] 在偶极子反演中，使用了一个祖巴（Zubal）大脑模型(见Zubal et al. Med Phys1994;21,pp.299-302)。不同脑区的幅度和磁化率值被赋予如表1所示的值。利用基于傅里叶变换的方法，相位可以从假设的磁化率分布直接计算得到。方差为0.01a.u.的随机零均值高斯噪声被独立的加在复图像的实部和虚部上。磁量图使用线性和非线性QSM方法来得到。它们相对于真实值的差别被绘制出来。
表1.不同脑区的仿真幅度和磁化率值
区域 幅度(a.u) 磁化率值(ppm)
白质 1 -0.05
静脉血 0.7 0.3
尾状核 0.6 0.09
壳核 0.6 0.09
丘脑 0.7 0.07
苍白球 0.5 0.19
出血 0 1
其他脑区 1 0
B.反卷折误差扩散的消除

[0109] 这里对一个颅内出血的病人使用不带流动补偿的单极多回波梯度回波序列在一台3T MR扫描仪(西门子医疗集团，埃朗根，德国)(Trio Tim Siemens Medical Solutions,Erlangen,Germany)中进行了扫描。其他参数包括：TR＝45ms，TEs=5,10,15,20ms，带宽200Hz/体素，视场24cm，层厚2mm，矩阵大小256×256×64。幅度和相位DICOM图像被同时保存以进行后处理。值得注意的是相位图上存在开放性条纹线，这对空域反卷折带来潜在困难。
C.在志愿者中比较线性和非线性QSM

[0110] 在一个示例性实施例中，五个健康志愿者(3男2女，年龄37±16)被召集在一个3T扫描仪(通用电气医疗集团，瓦克夏，威斯康星州，美国)(HDx,GE Healthcare,Waukesha,WI,USA)中进行脑成像。扫描使用了在所有回波的读出方向（前/后）上带流动补偿的单极多回波梯度回波序列。其他参数包括：TR＝54ms，从6ms到48ms均匀分布的8个TE，带宽300Hz/体素，视场24cm，层厚2mm，采集矩阵416×320×60，重建到512×512×60。采集中使用了加速系数为2的并行采集技术(阵列空间敏感性编码技术，Array Spatial Sensitivity Encoding Technique)。采集中同时使用了系数为0.75部分层编码来进一步缩短采集时间，总的采集时间少于7分钟。实部和虚部DICOM图像被保存来进行后处理。本发明提出的非线性QSM方法被用来和线性QSM方法进行比较。
D.颅内出血成像

[0111] 十三个去除了身份信息的颅内出血病人被用来作分析。这些图像来自于常规病人扫描，使用跟志愿者相同的脉冲序列在3T扫描仪(HDx,GE Healthcare)中进行脑成像。在临床扫描时使用以下参数：TR＝45ms，从5ms到40ms均匀分布的8个TE，带宽244Hz/体素，视场24cm，层厚2mm，采集矩阵240×240×28~54，重建到256×256×28~60。扫描中使用了系数为0.70的部分层编码。总扫描时间正比于扫描层数（~10层/分钟）。幅度和相位DICOM图像被保存来进行后处理。原始的线性QSM、使用MERIT的线性QSM、不使用MERIT的非线性QSM和本发明提出的非线性QSM方法被同时用来进行重建以理解MERIT的作用。一个经验丰富的放射科医生在未知重建方法的情况下对使用线性和非线性QSM方法得到的重建图像进行了浏览，并对图像伪影进行了评分。因为椒盐噪声和棋盘状图案、显现在轴状平面是阴影的条状伪影构成QSM中的主要质量问题。图像评分基于噪声和伪影的严重程度，每种评分分为三级：0=没有，1=适中，2=严重。两种评分之和（从0到4）被记录下载。线性方法和新方法的重建图像的平均评分被记录下来。威氏符号等级（Wilcox signed-rank）测试被用来检测评分间是否有显著性不同。
6.结果
信息论度量项评估

[0112] 在一个示例性例子中，非正则化重建图像有严重的条状伪影，在矢状面最明显。在联合熵正则化重建中（图5），这些条状伪影被成功的消除。使用解剖先验条件的解（图
5a&d）、非正则化的解（图5b&e）和使用联合熵正则化的解（图5c&d）的矢状面和轴状面如图
5所示。如图5b&e所示，使用正则化项的解成功消除了条状伪影。
非线性数据保真项评估
A.相位噪声扩散的减少

[0113] 线性和非线性拟合的频率图估计的比较如图6所示。图6a显示了一个线性拟合方法失败的示例性案例，噪声干扰了较晚回波的相位值。使用非线性拟合方法减少了复平面中的噪声干扰问题（图6b）。图6b中的最右边的点610对应的是第一个TE对应的信号测量值，其信号幅度最强。随着TE的增加，信号在复平面上以逆时针方向向（0,0）点逼近。两种方法的频率图估计都没有显著的偏差（线性和非线性拟合对应值分别为：均值=-0.39Hz和-0.12Hz,P=0.44和0.42）（图6c&d）。但是，非线性拟合显著地减少了误差的标准差（线性和非线性拟合对应的方差为16.1Hz和4.8Hz）。

[0114] 数值模型的偶极子反演中使用线性和非线性数据保真项的比较如图7所示。因为出血710（图7a）区域的信号强度很低，在出血区域的场信号被噪声干扰（图7c）。当使用线性偶极子反演时，干扰的场信号导致出血区域的磁化率值遭到干扰（图7d）。非线性拟合消除了出血区域中的严重噪声干扰，从而使得磁化率值相对于真值的准确性大大提高。在使用线性QSM和非线性QSM重建时，出血区域的均方误差分别为0.39ppm和0.04ppm，相对误差分别为39%和4%。B.反卷折误差扩散的消除

[0115] 直接从扫描仪上采集到的幅度图的矢状面和轴状面如图8a和图8f所示。相位图中的开放性条纹线（图8b&g）导致频率图（图8c&h）出现不完全空间解卷折，导致了在使用线性QSM时的严重条状伪影（图8d&i），其原因是线性数据保真项无法对此突变不连续进行建模。如图8e&j，非线性QSM重建方法消除了该条状伪影。线性和非线性QSM的重建时间都为4分钟。C.在志愿者中比较线性和非线性QSM

[0116] 在5个示例性测试例的其中3个中，非线性QSM在图像质量上对于线性QSM在以下区域中出现的伪影有改善：1）血管内出现的椒盐噪声；2）血管外出现的棋盘状图案和黑色伪影；3）基底节区的高磁化率区域出现的棋盘状图案。其中一个例子如图9所示。图9a所示为该例的幅度图。线性和非线性QSM重建图像目测时有高度相似性（图9b&c）。虚线框910A和915A的放大图如图9d&e，即915A和915B所示。如箭头所示，非线性QSM（图9e）比线性QSM（图9d）有细微提高。其他区域的线性和非线性QSM重建图像目测时有高度相似性。线性和非线性的重建时间都为17±2分钟。
D.颅内出血成像

[0117] 三个代表性例子的幅度图如图10a,f&k所示。 对比来看，MERT有效的降低了阴影状伪影（图10c,h,m），非线性数据保真项有效的降低了出血区域的噪声，这些发现跟仿真结果一致（图7e）。同时使用MERIT和非线性数据保真项的非线性QSM（图10e，j&o）对比线性QSM（图10b,g&l）在图像质量上有显著提高（P<0.001）。线性QSM的平均评分是2.6±1.1（均值±标准差），非线性QSM是0.5±0.7。非线性QSM的伪影评分一致低于线性QSM。线性和非线性QSM的重建时间都为3±1分钟。进一步示例性实施例

[0118] 在上述某些示例性实施例中，稳健性可能是阻碍常规QSM实践的一个瓶颈。本发明所示的某些示例性实施例中，可以使用非线性QSM重建方法。这种方法正确考虑了相位噪声和反卷折错误，降低了未期望错误的权重。非线性QSM在保持与线性QSM在相同重建时间的情况下，显著提高了低信噪比区域的解的精度。图像质量的提高有助于对颅内出血进行成像。

[0119] 颅内出血（ICH）是一种毁灭性疾病，在美国首次中风中比例10-15%，在中国的中风中比例30%(见Broderick et al. Circulation 2007;116,pp.e391-413)，并且在所有中风类型中有最高的死亡率，一年存活率小于50%(见Flaherty et al. Neurology2006;66,pp.1182-1186)。由于MRI能提供丰富的组织对比、对ICH存在具有高敏感性(见Kidwell et al. JAMA 2004;292,pp.1823-1830)，MRI在急性ICH临床实践和治疗研究中扮演着重要角色。传统的MRI技术对ICH的探测有局限性，并且通常是非量化的。ICH治疗需要对其进行量化，其中包括对血肿体积的测量(见Christoforidis et al. Stroke *
2007;38,pp.1799-1804)。QSM方法通过反卷积可以消除T2 加权幅度图像中的开花状伪影，可以得到可靠的血肿体积的量化值。血液随在血肿时间变化过程中从氧合血红蛋白(超急性ICH)、脱氧血红蛋白(急性),高铁血红蛋白(亚急性)含铁血黄素和铁蛋白(慢性)逐渐退化，QSM可以用来精确测量这种随血肿时间变化的磁化率值的变化。因此QSM可以对ICH病人管理提供有用的量化信息。

[0120] 非线性QSM的一个潜在示例性用途将是定量阴性对比剂，例如超顺磁性氧化铁（SPIO），该类对比剂拥有很强的磁化率值，在梯度回波上会造成类似出血的信号盲区。临床前试剂的SPIO定量将能有效监测癌症治疗或者发炎中的药物输送。类似的，QSM可以用来量化钆基阳性对比剂。该类对比剂已经被广泛使用在临床实践中，但是在传统的T1/T2增强中由于容易受到包括信号饱和、偏转角不精确、运动等很多错误影响，进而很难进行在体定量(见de Rochefort et al. Med Phys 2008;35,pp.5328-5339;Schabel and Parker.Phys Med Biol 2008;53,pp.2345-2373)。在体精确定量可以提供精确的诸如ktr等跟踪动力学参数，从而提高诸如脑肿瘤等疾病的诊断特异性。

[0121] 在技术上，非线性QSM的改进就是更好的管理MR信号中的噪声。MEDI中的线性最小二乘拟合仅在相位中的噪声是高斯的情况下，才能得到统计最优解。健康志愿者的SNR比较高，相位噪声可以假设为高斯的，所以使用线性QSM是可以的。但有出血病变时，这个假设将不再成立，例如图4。在这种情况下，非线性QSM方法正确考虑了实部和虚部中的零均值高斯噪声，利用偶极子场的长程本质，出血区域的磁化率值可以通过周围磁场来估计。因此，非线性QSM显著减少了病变区域的噪声。

[0122] 反卷折错误可能是阻碍稳健QSM重建的一个主要问题。虽然当前有很多反卷折方法，但是传统的在QSM中集成相位反卷折的方法不总是可靠的。由于偶极子场的长程本质，局部的相位跳变会在频率图中波及到很远的位置，例如图8。在本发明的某些示例性实施例中，非线性频率图估计中的空域反卷折被完全绕过，从而消除了空域反卷折错误。有几个示例性方案可以进一步减少空域反卷折错误。如果其他方法(见Chen et al. Neuroimage2010;49,pp.1289-1300;Robinson et al. Magn Reson Med 2011;65,pp.1638-1648)用来将多通道原始数据合成一个单一的相位图，则可以完全避免相位图中的开放性条纹线（图
8b&g）。如图8c&h中的空域不连续可以使用基于最小范数的空域反卷折方法来消除(见Li et al. Neuroimage 20 11;55,pp.1 645-1 656;Ghiglia and Pritt.1998;xiv,493 p.178-277;Schofield and Zhu.Optics Letters 2003;28,pp.1194-1196)。但是，因为反卷折和卷折的相位的差几乎从来不是2π的倍数，所以基于最小范数的方法只是定性的，破坏了磁化率图的量化本质。作为一个完全绕过空域反卷折的尝试，某些示例性实施例用复平面的非线性拟合重写了PDF中的最小二乘拟合问题的同时，使用了与非线性PDF相同的迭代求解算法。为了使非线性PDF快速收敛，可以使用线性PDF来得到初始背景场偶极子分布。图8中的某些示例性非线性PDF结果可能相比线性PDF的结果没有显著性提高，表明对比偶极子反演而言背景场移除中的反卷折错误是一个较小的问题。所以根据某些示例性实施例，对比线性QSM而言，完全非线性的QSM步骤将是可实行的且对空域反卷折错误更稳健；这是因为相位被映射到复平面上进而对错误的2π跳变免疫。

[0123] 示例性MERIT方法类似于迭代冲加权最小二乘（iteratively reweighted least squares）(IRLS)方法。IRLS通常作为减少非正态分布数据集中的异常值的影响而用在稳健的回归中。IRLS可以使用Lp范数来最小化数据保真项，其中p=2-2t(见1996;xvii,408 p.173-405)，t定义在实现部分。某些示例性实施例中的方法跟传统的IRLS有区别。这些示例性实施例中选择性地对一些体素使用了递减的权重，而传统的IRLS对所有的体素使用重加权（减少或增加权重）。某些示例性实施例的初步研究表明MERIT在t=2时能有效地减少与不能很好建模的被干扰数据点有关的伪影。场模型错误可能在高磁化率区域（～1ppm）很明显,这将导致在TE=20ms体素内场差别为 这种情况
也可能出现在有流动等运动的位置上。

[0124] 虽然本发明提出的非线性最小二乘拟合使用的L1范数正则化项倾向于与解剖结构相同，但是非线性数据保真项可以使用各种正则化方法，不一定与这个特殊的正则化方法绑定。非线性拟合和MEDI可以很好的匹配是因为MEDI中的L1正则化项已经需要非线性求解，而增加的非线性数据保真项不需要额外的求解和计算时间。示例性方法小节中概述的线性化过程工作良好，使得非线性收敛性能与线性QSM相同。

[0125] 非线性频率图估计方法可能类似于基于回波不对称性与最小二乘估计的迭代水脂分解(the Iterative Decomposition of water and fat with Echo Asymmetry and Least-squares estimation)(IDEAL)方 法( 见 Reeder et al. Magn Reson Med2005;54,pp.636-644)。因为未知数数量比IDEAL少且有更多的回波可用，某些示例性实施例中的频率图估计可能比IDEAL更稳健。水脂分离中一个已知问题是解常取决于初始值，这将会得到局部而不是全局最小解。在某些示例性实施例中，这个问题可以被避免，通过使用前两个回波来计算得到已经靠近最终解的初始频率。但是偶极子反演中的数据保真项的收敛性依旧受初始解的影响。这个问题可以通过在某些示例性实施例中使用更靠近最终解的初始解来解决，例如使用线性QSM的估计结果或者使用更好的求解法，例如基于图割的方法(见Hernando et al. Magn Reson Med 2010;63,pp.79-90)。

[0126] 本发明中的示例性实施例包括非线性QSM重建方法，其中适当考虑了相位噪声和反卷折错误，以及为了更稳健的拟合而使用的重加权方法。通过使用此非线性方法克服了相位反卷折问题和非高斯相位噪声问题。数值仿真和病人数据都表明此方法对低信噪比区域的重建有显著提高。图像质量的提高特别对颅内出血成像有帮助。示例性应用

[0127] 在以下示例性应用中，复MR信号使用非线性场图估计、隐含背景场移除、MEDI正则化项和非线性数据保真项进行处理。可视化脑深部核团和白质纤维束

[0128] 因为铁具有强顺磁性且深部脑核团通常相对于周围组织具有丰富的铁沉积，磁量图可以提供清晰的深部脑核团和白质纤维束的划分，如图11所示。在这个示例性例证中，复MR信号使用以下示例性参数的多回波梯度回波序列在一个3T扫描仪中进行脑成像。扫描使用在读出方向（前/后）带流动补偿、TR＝54ms、从3.8ms到48.9ms均匀分布的
12个TE、带宽365Hz/体素、轴向视场24cm、层厚0.7mm、采集矩阵342×342×150、重建到
512×512×150。采集中使用了加速系数为2的并行采集技术和系数为0.7部分层编码来缩短采集时间。深部脑核团包括齿状核、红核、黑核、丘脑、苍白球、壳核和尾状核可以清晰地从周围脑实质中分辨出来，如图11中的箭头1110所示。此外白质纤维束包括视放射、胼胝体膝部、压部、内囊和外囊显示出相比周围脑组织具有负磁化率值，如图11的箭头11 15所示。这些脑结构的清晰可视化提供了有用的脑结构信息。
在脑深部刺激中可视化丘脑底核

[0129] 在众多脑核中，丘脑底核由于处在深部脑刺激的位置而特别受重视。对该区域的更好的可视化对治疗病人时的手术方案特别有用，例如帕金森病人。在QSM（图12a）、磁化率加权成像（图12b）、幅度图（图12c）、R2*图（图12d）和T2加权图（图12e）中，QSM能提供丘脑底核和周围组织最好的对比。这个案例的复MR信号使用以下参数的多回波梯度回波序列在一个3T扫描仪中进行脑成像。扫描使用在读出方向（上/下）上带流动补偿、TR＝45ms、从4ms到40ms均匀分布的11个TE、带宽390Hz/体素、轴向视场24cm、层厚2mm、采集矩阵320×320×40、重建到512×512×40。采集中使用了加速系数为2的并行采集技术和系数为0.7部分层编码来缩短采集时间。为了对比，也采集了T2加权图像（T2WI），其采集参数为TE=86ms，TR＝8s，带宽=390Hz/体素，平均次数=4。两种方法采集时间都小于5分钟。在T2WI和从同一个多回波图像中得到的QSM、磁化率加权图像、幅度图、R2*图中，QSM提供了最佳对比噪声比的丘脑底核的最佳可视化图像。对比噪声比定义为丘脑底核区的峰值信号与相邻组织的信号差除以相邻丘脑感兴趣区域中信号的标准差。
在多发性硬化症患者中可视化白质病变

[0130] 磁量图也可以提供多发性脑硬化病人的白质病变的额外信息，如图13所示。病变在QSM（图13a）、T1加权图像（图13b）、T2液体衰减反转恢复序列（FLAIR）图像（图13c）和T2加权图像（图13d）中可能会有不同的表现，从而提供了对病变分级、疾病分期的可能性。在这个例证中，复MR信号使用以下参数的多回波梯度回波序列在一个3T扫描仪中进行脑成像。扫描使用读出方向（前/后）上的带流动补偿、TR＝57ms、从4.3ms到52.4ms均匀分布的11个TE、带宽300Hz/体素、轴向视场24cm、层厚2mm、采集矩阵416×320×60、重建到512×512×60。采集中使用了加速系数为2的并行采集技术和系数为0.7部分层编码来缩短采集时间。在重建的QSM图像中，道森指（Dawson Finger）是高信号病变，可以通过其形态和相关的髓静脉1310来识别。图中也可以看到其他椭圆形高信号病变，表示铁沉积或者脱髓鞘。这些跟病变物理特征相对应的额外信息可以给病变分期和疾病分级提供相对于T1加权、T2加权和T2FLAIR图像的补充信息。
示例性癌症成像

[0131] 一个示例性磁量图可以用在脑肿瘤中，如图14所示。在这个示例性的复MR信号使用以下示例性参数的多回波梯度回波序列在一个3T扫描仪中进行脑成像。扫描在读出方向（前/后）上带流动补偿、TR＝57ms、从4.3ms到52.4ms均匀分布的11个TE、带宽300Hz/体素、轴向视场24cm、层厚2mm、采集矩阵416×320×60、重建到512×512×46。采集中使用了加速系数为3的并行采集技术和系数为0.7部分层编码来缩短采集时间。脑肿瘤中的磁化率值从负磁性到正磁性都有。对比剂成像

[0132] 例如钆（Gd）和SPIO粒子等的对比剂已经在体模和体外实验中证明是顺磁性的，而且其磁化率值与浓度和摩尔磁化率成线性关系。因此磁量图提供了在体对比剂量化的方法。

[0133] 在一个示例中（如图15所示），脑肿瘤分别在对比剂注射之前和之后成像。这个示例性例证使用以下示例性参数的多回波梯度回波序列进行脑成像。扫描在读出方向（前/后）上带流动补偿、TR＝45ms、从3ms到36ms均匀分布的12个TE、带宽521Hz/体素、轴向视场24cm、层厚3mm、采集矩阵240×240×52、重建到256×256×52。采集中使用了加速系数为2的并行采集技术来缩短采集时间。对比注射前和注射后的图像，外围组织增强更明显（图15a&b），表示肿瘤的外围高血管系统特性。而同时幅度图上显示几乎是均匀增强，可能是由于信号强度和钆浓度的非线性关系而达到饱和。QSM测量的钆浓度在外围是0.22mM，中央是0.017mM。明确的钙化鉴别

[0134] 磁量图可以用来进行明确的钙化鉴别，因为钙化在人脑常见已知物质中具有独特的逆磁性特征。钙化病变在QSM中是低信号的，可以很容易的从其他病变中区分和鉴别出来。一个示例性例证如图16所示。这个示例性例证的复MR信号使用以下示例性参数的多回波梯度回波序列在一个3T扫描仪中进行脑成像。扫描在读出方向（前/后）上带流动补偿、TR＝53ms、从5ms到40ms均匀分布的8个TE、带宽365Hz/体素、轴向视场24cm、层厚2mm、采集矩阵416×320×70、重建到512×512×70。采集中使用了加速系数为2的并行采集技术来缩短采集时间。为了确认钙化存在，病人额外进行了计算机断层成像（CT）扫描，如图16b所示。如果钙化病变的总磁化率值定义为手工分割的感兴趣区域的磁化率值的和与体素大小之积，且总亨斯菲尔德（Hounsfield）单位定义使用CT中的定义，则钙化病变的总磁化率值和总Hounsfield单位具有很好的相关性，如图16c所示。
7.可视化

[0135] 因为QSM是后处理技术，生成的图像自动会对准到输入图像上，从而可以被融合到源图像中。在某个示例性实施例中，QSM图像作为前景图像被融合进作为背景图像的幅度图像中。前景图像的透明度可以从0%到100%变化，而且可以使用各种伪彩色进行彩色显示。显示窗口和背景图相对于前景图像可以独立调节，如图17中的一系列示例性图像所示。体积QSM图像可以进行三维体渲染显示和基于数层的平均密度投影、最小密度投影、最大密度投影显示。附加的示例性实施例

[0136] 图18显示了按照本发明中的一个系统的示例性实施例所绘的示例性方框图。例如本发明中的示例性实施例可以由所述处理装置和/或一个计算装置1810来实现。该处理/计算装置1810可以是，例如全部或者一部分、包括但不仅限于计算/处理1820。计算/处理1820可以包括，例如一个或多个处理器，使用存储在计算机访问介质（如随机存储器（RAM）、只读存储器（ROM）、硬盘或其他存储设备）中的指令。

[0137] 如图18所示，例如1830的计算机访问介质（正如上面所述,存储设备如硬盘、软盘、记忆棒、光盘、RAM、ROM等，或其中一个集合）可以提供来于例如处理装置1810进行通信。计算机访问介质1830中可以包含可执行指令1840。额外或者备选的，可以提供独立于计算机访问介质1830的存储装置1850，为处理装置1810提供指令配置用于执行例如如上所述的某些示例性实施例、过程和方法。

[0138] 此外，示例性处理装置1810可以提供带有或包括一个输入/输出装置1870；1870可以包括例如有线网络、无线网络、互联网、内部网、数据采集探头和传感器等。如图18所示，示例性处理装置18可以与示例性显示装置1860进行通信；1860按照本发明的示例性实施例可以配置为给处理装置提供输入信息和输出信息而使用的设备，例如触摸屏。此外，示例性显示装置1860和/或存储装置1850可以以用户可访问和/或用户可读的格式来显示和/或存储数据。事实上，示例性处理装置可以按照本发明中所述示例性实施例来执行上述的示例性过程。

[0139] 上述仅仅说明了本发明的原理。本领域的技术人员很明显可以按照上述指导对所述实施例进行修正和改动。本领域的技术人员设计出众多的以上没有明确显示或描述的、体现本发明的原理的系统、装置和程序都属于本发明的精神和范围。此外所有出版和引用以上所述的，在其所有权中应该将本发明作为引用全部纳入。应当理解上述示例性程序可以存储在任何计算机访问介质中，包括硬盘、RAM、ROM、移动硬盘、光盘、记忆棒等中，且可以被执行在处理装置和/或计算机装置中，可以是和/或包括硬件处理器、微处理器、小型、宏、大型机等，以及其多元和/或组合。此外，本发明中使用的某些条款包括其中的描述、图表和要求在某些情况下可以同义地使用，包括但不限于数据和信息。应当理解，虽然这些词和/或其他词彼此同义，可以作为同义词使用，有可能出现故意不使用同义词的情况。此外，作为先验知识而没有显式引用以上所述，应当在其所属权中明确的纳入本发明。所有引用本发明的出版物在其所属权中应列入本发明。参考文献
1.Li，W.，Wu，B.，Liu，C.，201 1.Quantitative susceptibility mapping of human brain reflects spatial variation in tissue composition.Neuroimage 55，
1645-1656.
2.Liu，C.，2010.Susceptibility tensor imaging.Magn Reson Med 63，1471-1477.
3.Wu，B.，Li，W.，Guidon，A.，Liu，C.，201 2.Whole brain susceptibility mapping using compressed sensing.Magn Reson Med 67，137-147.
4.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 usingMRI phase data.Magn Reson Med 62，1510-1522.
5.Haacke，E.M.，Tang，J.，Neelavalli，J.，Cheng，Y.C.，201 0.Susceptibility mapping as a means to visualize veins and quantify oxygen saturation.J Magn Reson Imaging 32，663-676.
6.Liu，T.，Liu，J.，de Rochefort，L.，Spincemaille，P.，Khalidov，I.，Ledoux，J.R.，Wang，Y.，201 1.Morphology enabled dipole inversion(MEDI)from a single-angle acquisition：comparison with COSMOS in human brain imaging.Magn Reson Med 66，
777-783.
7.Liu，T.，Spincemaille，P.，de Rochefort，L.，Kressler，B.，Wang，Y.，2009.Calculation of susceptibility through multiple orientation sampling(COSMOS)：a method for conditioning the inverse problem from measured magneticfield map to susceptibility source image in MRI.Magn Reson Med 61，1 96-204.
8.Schweser，F.，Deistung，A.，Lehr，B.W.，Reichenbach，J.R.，2011.Quantitative imaging of intrinsic magnetic tissue properties usingMRI signal phase：an approach to in vivo brain iron metabolism?Neuroimage 54，2789-2807.
9.Wharton，S.，Schafer，A.，Bowtell，R.，201 0.Susceptibility mapping in the human brain using threshold-based k-space division.Magn Reson Med 63，
1292-1304.
1 0.Wharton，S.，Bowtell，R.，201 0.Whole-brain susceptibility mapping at highfield：acomparison of multiple-and single-orientation methods.Neuroimage
53，515-525.
11.de Rochefort，L.，Brown，R.，Prince，M.R.，Wang，Y.，2008.Quantitative MR susceptibility mapping using piece-wise constant regularized inversion of the magneticfield.Magn Reson Med 60，1 003-1 009.
12.de Rochefort，L.，Nguyen，T.，Brown，R.，Spincemaille，P.，Choi，G.，Weinsaft，J.，Prince，M.R.，Wang，Y.，2008.In viyo quantification of contrast agent concentration usingthe induced magnetic field for time-resolved arterialinput function measurement with MRI.Med Phys 35，5328-5339.
13.Liu，T.，Spincemaille，P.，de Rochefort，L.，Wong，R.，Prince，M.，Wang，Y.，
201 0.Unambiguous identification of superparamagnetic iron oxide particles through quantitative susceptibility mapping of the nonlinear response to magneticfields.Magn Reson Imaging 28，1383-1389.
14.Li，L.，2001.Magnetic susceptibility quantification for arbitrarily shaped objects in inhomogeneous fields.Magn Reson Med 46，907-916.
15. Li，L.，Leigh，J.S.，2004.Quantifying arbitrary magnetic susceptibility distributions with MR.Magnetic Resonance in Medicine 51，1077-1082.
16.Li，L.，2002.Averaging of harmonic physicalfields over an annular region enclosing field sources.American Journal of Physics 70，1029-1033.
17.Li，L.，Leigh，J.S.，2001.High-precision mapping of the magnetic field utilizing the harmonic function mean value property.J Magn Reson 148，442-448.
18.Koch，K.M.，Papademetris，X.，Rothman，D.L.，de Graaf，R.A.，2006.Rapid calculations of susceptibility-induced magnetostaticfield perturbations for in viyo magnetic resonance.Phys Med Biol 51，6381-6402.
19.Marques，J.P.，Bowtell，R.，2005.Application of a fourier-based method for rapid calculation of field inhomogeneity due to spatial variation of magnetic susceptibility.
Concepts in Magnetic Resonance Part B-Magnetic Resonance Engineering 25B，
65-78.
20.Salomir，R.，De Senneville，B.D.，Moonen，C.T.W.，2003.A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility.
Concepts in Magnetic Resonance Part B-Magnetic Resonance Engineering 19B，
26-34.
21.Park，S.C.，Park，M.K.，Kang，M.G.，2003.Super-resolution image
reconstruction：a technical overview.IEEE Signal Processing Magazine 20，16.
22.Farsiu，S.，Robinson，M.D.，Elad，M.，Milanfar，P.，2004.Fast and robust multiframe super resolution.IEEE transactions on image processing：a publication of the IEEE Signal Processing Society 13，1327-1344.
23.Gudbjartsson，H.，Patz，S.，1995.The Rician distribution of noisy MRI data.Magn Reson Med 34，910-914.
24.Elad，M.，Feuer，A.，1 997.Restoration of a single superresolution image from several blurred，noisy，and undersampled measured images.IEEE transactions on image processing：a publication of the IEEE Signal Processing Society 6，
1646-1658.
25.Gindi，G.，Lee，M.，Rangarajan，A.，Zubal，I.G.，1993.Bayesian reconstruction of functionalimages using a natomicalinformation as priors.Ieee Transactions on Medical Imaging 12，670-680.
26.Baillet，S.，Garnero，L.，1997.A Bayesian approach to introducing anatomo-functional priors in the EEG/MEG inverse problem.IEEE transactions on bio-medical engineering 44，374-385.
27.Kressler，B.，de Rochefort，L.，Liu，T.，Spincemaille，P.，Jiang，Q.，Wang，Y.，2010.Nonlinear regularization for per voxel estimation of magnetic susceptibility distributions from MRI field maps.Ieee Transactions on Medical Imaging 29，273-281.
28.Liu，J.，Liu，T.，de Rochefort，L.，Ledoux，J.，Khalidov，I.，Chen，W.，Tsiouris，A.J.，Wisnieff，C.，Spincemaille，P.，Prince，M.R.，Wang，Y.，2012.Morphology enabled dipole inversion for guantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map.Neuroimage
59，2560-2568.
29.Goldstein，R.M.，Zebker，H.A.，Werner，C.L.，1988.Satellite Radar Interferometry-Two-Dimensional Phase Unwrapping.Radio Science 23，713-720.
30.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 bayesian regularization：validation and application to brain imaging.Magn Reson Med 63，1 94-206.
31.Kressler，B.，de Rochefort，L.，Liu，T.，Spincemaille，P.，Jiang，Q.，Wang，Y.，
2010.
Nonlinear regularization for per voxel estimation of magnetic
susceptibility distributions from MRI field maps.IEEE Trans Med Imaging 29，
273-281.
32.Liu，J.，Liu，T.，de Rochefort，L.，Ledoux，J.，Khalidov，I.，Chen，W.，Tsiouris，A.J.，Wisnieff，C.，Spincemaille，P.，Prince，M.R.，Wang，Y.，2012.Morphology enabled dipole inversion for guantitative susceptibility mapping using structural consistency between the magnitude image and the susceptibility map.Neuroimage.
33.Rudin，L.I.，Osher，S.，Fatemi，E.，1992.Nonlinear Total Variation Based Noise Removal Algorithms.Physica D 60，259-268.
34.Liu，T.，Khalldov，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.
35.Smith，S.M.，2002.Fast robust automated brain extraction.Hum Brain Mapp
17，143-155.
36.Zubal，I.G.，Harrell，C.R.，Smith，E.O.，Rattner，Z.，Gindi，G.，Hoffer，P.B.，
1994.
Computerized three-dimensional segmented human anatomy.Med Phys 21，
299-302.
37.Broderick，J.，Connolly，S.，Feldmann，E.，Hanley，D.，Kase，C.，Krieger，D.，Mayberg，M.，Morgenstern，L.，Ogilvy，C.S.，Vespa，P.，Zuccarello，M.，2007.Guidellnes for the management of spontaneous intracerebral hemorrhage in adults：2007 update：a guidellne from the American Heart Association/American Stroke Association Stroke Council，High Blood Pressure Research Council，and the Quality of Care and Outcomes in Research Interdisciplinary Working Group.Circulation
116，e391-413.
38.Flaherty，M.L.，Haverbusch，M.，Sekar，P.，Kissela，B.，Kleindorfer，D.，Moomaw，C.J.，Sauerbeck，L.，Schneider，A.，Broderick，J.P.，Woo，D.，2006.Long-term mortality after intracerebral hemorrhage.Neurology 66，1182-1186.
39.Kidwell，C.S.，Chalela，J.A.，Saver，J.L.，Starkman，S.，Hill，M.D.，Demchuk，A.M.，Butman，J.A.，Patronas，N.，Alger，J.R.，Latour，L.L.，Luby，M.L.，Baird，A.E.，Leary，M.C.，Tremwel，M.，Ovbiagele，B.，Fredieu，A.，Suzuki，S.，Villablanca，J.P.，Davis，S.，Dunn，B.，Todd，J.W.，Ezzeddine，M.A.，Haymore，J.，Lynch，J.K.，Davis，L.，Warach，S.，2004.
Comparison of MRI and CT for detection of acute intracerebral hemorrhage.JAMA 292，1823-1830.
40.Christoforidis，G.A.，Slivka，A.，Mohammad，Y.，Karakasis，C.，Avutu，B.，Yang，M.，2007.Size matters：hemorrhage volume as an objective measure to define significant intracranial hemorrhage associated with thrombolysis.Stroke 38，
1799-1804.
41.Schabel，M.C.，Parker，D.L.，2008.Uncertainty and bias in contrast concentration measurements using spoiled gradient echo pulse sequences.Phys Med Biol 53，2345-2373.
42.Chen，Z.，Johnston，L.A.，Kwon，D.H.，Oh，S.H.，Cho，Z.H.，Egan，G.F.，2010.An optimised framework for reconstructing and processing MR phase images.Neuroimage 49，1289-1300.
43.Robinson，S.，Grabner，G.，Witoszynskyj，S.，Trattnig，S.，2011.Combining phase images from multi-channel RF coils using 3D phase offset maps derived from a dual-echoscan.Magn Reson Med 65，1638-1648.
44.Ghiglia，D.C.，Pritt，M.D.，1998.Two-dimensional phase unwrapping：theory，algorithms，and software.Wiley，New York.
45.Schofield，M.A.，Zhu，Y.M.，2003.Fast phase unwrapping algorithm for interferometric applications.Optics Letters 28，1194-1196.
46. 1996.Numerical methods for least squares problems.SIAM，
Philadelphia.
47.Reeder，S.B.，Pineda，A.R.，Wen，Z.，Shimakawa，A.，Yu，H.，Brittain，J.H.，Gold，G.E.，Beaulieu，C.H.，Pelc，N.J.，2005.Iterative decomposition of water and fat with echo asymmetry and least-squares estimation(IDEAL)：application with fast spin-echo imaging.Magn Reson Med 54，636-644.
48.Hernando，D.，Kellman，P.，Haldar，J.P.，Liang，Z.P.，2010.Robust water/fat separation in the presence of large field inhomogeneities using a graph cut algorithm.MagnReson Med 63，79-90.