CN109389651A - 磁共振化学位移编码成像方法与装置 - Google Patents

Abstract

本发明涉及磁共振化学位移编码成像方法与装置以及执行所述方法的各步骤的计算机可读介质。所述方法包括：(1)求得二维原始磁共振图像的两组场图候选解；(2)在不同的尺度下，将前述两组场图候选解分别分割为不交叠的子块；(3)在每一个尺度下分别采用局部迭代场图提取算法，得到相应的原始图像空间下的场图；(4)将各个尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取，得到最终场图；(6)对最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图。本发明所述方法能实现正确的场图估计和水脂分离。

Description

磁共振化学位移编码成像方法与装置

技术领域

本发明涉及磁共振成像技术领域，尤其涉及基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法与装置。

背景技术

磁共振化学位移编码成像是一种基于组织中各成分之间化学位移差异的成像方法。化学位移编码成像在不同的回波时间采集信号，使得每一组信号中各成分之间存在相位差。临床上最常用的化学位移编码成像是水脂分离成像，它主要用于脂肪抑制和脂肪定量等应用中。

多回波化学位移编码成像假设组织信号由多种不同成分(本发明中以水和脂肪分离成像为例)被同时激发得到，并已知各成分相对于水的化学位移，然后拟合信号与各成分、B₀场图和回波时间(echo time,TE)之间的数学模型。实际应用最多的化学位移编码成像是水脂分离成像，即通过多回波信号分离水和脂肪两种成份，以下均以水脂分离成像为例说明两点化学位移编码成像方法。

1984年，Dixon^[1]最早提出了两点水脂分离技术，该技术假设场图分布均匀，选取两个特定的TE分别采集水和脂肪质子相位同向(In-phase，IP) 和反向(opposed-phase，OP)的图像，通过简单的相加/相减操作得到水图和脂肪图。为了解决B0场不均匀造成的影响，2006年，Xiang^[2]提出了一种新的两点水脂分离技术，该技术不再限制第二幅采集图像中水和脂肪相位反向，而是采集水脂相位部分反向(partially-opposed-phase，POP)图像，利用余弦定理得到两组水和脂肪的幅值图像，在此基础上解到两个场图候选值然后将原始图像分割为多个不交叠的子块，根据子块内场图分布的光滑特性，选取出每一个子块内具有最大场图相似性的场图集合作为初始场图，利用区域迭代场图提取(RegionalIterative Phasor Extraction，RIPE)算法从场图候选解中将整个原始图像空间内的场图估计出来，最后利用最小二乘法得到分离的水和脂肪图像。2011年，Eggers等^[3]在Xiang方法的基础上提出了新的两点分离技术，该技术也不再限制第一幅采集图像中水和脂肪相位同向，只是限制第一幅采集图像中的水脂相位差θ₁和第二幅图像中的水脂相位差θ₂之间的差或和不能等于2π或其整数倍。该技术借用Xiang的方法对场图进行估计。

多回波化学位移编码成像面临的最大挑战是由主磁场B₀不均匀性 (本发明称之为场图，B₀均匀时，场图连续缓慢变化)造成的局部分反，因此该技术的分离效果依赖于B₀场图的估计。

前述现有算法存在其缺点，文献[2-3]在不同的子块得到各自的场图集合，将其直接拼接后作为初始场图进行区域迭代场图提取，这样做的缺陷是在场图提取过程中容易出现选取矛盾，即从不同方向进行场图选取时得到不同的选取结果，造成最后的场图中存在跳变。

发明内容

针对上述缺点，本发明提出一种基于多尺度局部迭代场图提取 (multi-scaleRIPE)算法的磁共振化学位移编码成像方法，以实现正确的场图估计和待分离成分的分离(例如，水脂分离)。

本发明所述的磁共振化学位移编码成像方法，所述方法包括：

(1)求得二维原始磁共振图像的两组场图候选解(又称场图候选值)；

(2)在不同的分块尺度下，将前述两组场图候选解分别分割为不交叠的子块；

(3)在每一个分块尺度下分别采用局部迭代场图提取算法，得到相应的原始图像空间下的场图；

(4)将各个分块尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；

(5)利用局部迭代场图提取算法对置0像素点的场图进行选取，得到最终场图；

(6)对最终场图进行平滑处理，分离原始磁共振图像(即，利用平滑处理后的场图对原始磁共振图像进行分离)，得到准确的各待分离成分(例如水和脂肪，或者水与另一待分离成分)的图。

在一实施方式中，所述待分离成分为水和脂肪。

在一实施方式中，所述待分离成分为水和另一成分，所述成分可以是硅胶或其它任意一种和水具有不同化学位移的成分；也可以分离任意两种具有不同化学位移的成分。

在一实施方式中，所述的步骤(1)根据以下执行：

描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线，并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值，构成该像素点的两组场图候选解。

在一子实施方式中，所述的步骤(1)可根据以下执行：根据式(5) 描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线(phasor-error 谱)，

p＝argminerr(p)＝argmin||(I-AA⁺)S||² (5)

其中p是相位矢量(或相位因子)，也称为场图值；I是4×4单位矩阵； A⁺＝(A^TA)^-1A^T，上标“+”表示取矩阵伪逆，上标“T”表示矩阵的转置； err(p)描述了相位矢量p和拟合误差之间的关系，表示为phasor-error谱，其定义为不同相位矢量值下对应的模型拟合误差；S如下面式(4)中所定义；

根据该phasor-error谱寻找使式(5)表征的二范数达到局部极小值的所有场图值p，构成该像素点的场图候选解，个数为1个或2个，以符号p₁和 p₂表示，场图值数目为1时，p₁＝p₂。

在一实施方式中，所述式(5)通过以下获得：

采用以下式(1)的含有水和另一成分(例如脂肪)的两种成分的MR信号模型：

其中TE为回波时间，S_n是在回波时间TE_n下的信号强度；ρ_w和ρ_f是水和另一成分(例如脂肪)的信号强度值；水的化学位移设定为0，f_F是另一成分相对于水的化学位移，为已知量，例如脂肪，为-3.5ppm；f_B是静磁场B₀的场图值，i表示虚数；

利用式(2)定义以下参数c₁、c₂、q和p：

其中ΔTE＝TE₁-TE₂，q代表TE₁时间内由于f_B造成的水的相位，p代表ΔTE时间内由于f_B造成的水的相位累积，即如上所定义的相位矢量或场图值，其角度与f_B成正比，(角度)范围是[-ππ]；

其中利用式(2)的p，由式(1)和(2)解得以下式(3)表示的q：

通过利用式(2)、(3)，将式(1)改写成以下式(4)的矩阵形式：

其中是4×1矩阵， S₁、S₂分别独立地为在回波时间TE₁、TE₂下一个像素点的信号强度，Re代表取复数实部，Im代表取复数虚部；

通过式(4)，根据变量投影方法，使用最小二乘法解得式(5)表示的 p：

p＝argminerr(p)＝argmin||(I-AA⁺)S||² (5)。

在一实施方式中，所述的步骤(2)根据以下执行：

根据预先定义的分块尺度值将待分割的二维原始磁共振图像进行分块，得到的每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度。

在一实施方式中，所述分块尺度值为整数，选自1至16中的一个或多个数值。例如，所述分块尺度值可以是任选的4个数值：6、8、10和 12。

在一实施方式中，所述步骤(3)根据以下执行：

(3a)场图的初始化

对每一个分块尺度下的每一个子块，求得子块场图候选解，利用每一个子块场图候选解为初始值，对子块内所有像素点的两组场图候选解进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小，符合所述规则的像素点场图候选解组成一个场图值集合，未被选取的像素点场图候选解组成另一个场图值集合；计算每一个场图候选解对应场图值集合的矢量和；比较所述矢量和的幅值，选取最大幅值对应的场图值集合；

对每一个分块尺度下每一个子块均进行上述场图初始化，将每一个分块尺度下所有子块的场图值集合合并在一起，即得到该分块尺度下的初始场图值集合；

(3b)场图的更新

对上述初始场图值集合进行低通滤波，且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理；然后以每一个像素点的场图值为初始值，对图像内所有像素点的两组场图候选解再次进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小；

并记录场图值发生改变的像素点的个数N；

(3c)迭代更新场图

重复步骤(3b)，直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止，得到一组场图；

对应于多组分块尺度，经过步骤(3a)、(3b)和(3c)得到多组的场图。

在一子实施方式中，所述步骤(3)可根据以下执行：

(3a)场图的初始化

每一个分块尺度下的每一子块内所有像素点作为整体根据式(6)和(7)描述场图值拟合误差曲线，根据该场图值拟合误差曲线寻找使式(6)达到局部极小值的场图值，所述场图值为该子块的场图候选解，数目最多为2个，该子块的场图候选解以符号p_u和p_v表示，当场图值数目为1时，p_u＝p_v，

式(6)由式(5)改写得到，其中，I是4×4单位矩阵；A⁺＝(A^TA)^-1A^T，上标“+”表示取矩阵伪逆，上标“T”代表矩阵的转置，表示点乘；err(p)描述了场图值p和拟合误差之间的关系，表示为场图值拟合误差曲线，其定义为不同场图值p下对应的模型拟合误差，

式(7)中，和分别代表一组像素点的信号强度，定义如式(4)中表述，和都是行向量，行向量中每一个元素代表一个像素点的信号强度；

以上述子块的一个场图候选解为该子块的初始值，对该子块内所有像素点的两个场图候选解(p₁，p₂)进行选取，选取的规则是：计算该子块内每个像素点的两个场图候选解与子块初始值之间的差值，并对差值取绝对值，较小绝对值对应的场图候选解即为该像素点的场图值，以此方法对子块内所有像素点的场图值进行选取，并计算这些场图值的矢量和，同时计算所有像素点未被选取的场图候选解的矢量和，比较上述两组矢量和的幅值，保留较大值对应的幅值，舍弃另一个幅值；

同样地，以该子块另一个场图候选解为初始值，对该子块内所有像素点的场图值进行选取，也得到两组场图值矢量和的幅值，比较这两组幅值，保留较大值对应的幅值，舍弃另一个幅值；

比较该子块2个场图候选解所得到的保留幅值，并保留较大值的保留幅值对应的子块场图候选解，舍弃另一组较小值对应的场图候选解；根据保留的子块场图候选解对应的两组场图值的矢量和得到加权因子，进而得到该子块内所有像素点的场图值集合P；

或者，根据下文描述的具体方法对该子块内所有像素点的两个场图候选解(p₁，p₂)进行选取；

对每一个分块尺度下每一个子块均进行上述场图初始化，将每一个分块尺度下所有子块的场图值集合P合并在一起，即得到该分块尺度下的初始场图值集合P₀；

(3b)场图的更新

对上述初始场图值集合P₀进行低通滤波，且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理；

假定迭代后得到的场图值集合为P_n，则滤波后得到P_n ^s，然后以每一个像素点的场图值p_n ^s为初始值，对图像内所有像素点的两个场图候选解(p₁， p₂)再次进行选取，选取方法类似于步骤(3a)场图的初始化中对子块内所有像素点的两个场图候选解(p₁，p₂)的选取方法，

并记录场图值发生改变的像素点的个数N；

(3c)迭代更新场图

重复步骤(3b)，直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止，得到一组场图；

对应于多组分块尺度，经过步骤(3a)、(3b)和(3c)得到多组的场图。

在一实施方式中，在前述步骤(3a)中，在根据式(6)和(7)描述每一个分块尺度下的每一子块内所有像素点整体的场图值拟合误差曲线时，根据该场图值拟合误差曲线寻找使式(6)达到局部极小值的场图值时，设置以下阈值排除一定数目的场图值，使之数目为1个或2个，阈值设定为：

A_min+α*(A_max-A_min)，

其中A_max和A_min分别表示场图值拟合误差曲线的最大值和最小值，α的选取要排除拟合误差大的场图值，α取0.3-0.8。在一实施方式中，α取0.4。

在一实施方式中，在前述步骤(3a)中，对每一子块内所有像素点的两个场图候选解(p₁，p₂)进行选取的方法具体为：

针对一个子块，以该子块的场图值p_u和p_v为初始值，分别对该子块内所有像素点的两个场图候选解(p₁，p₂)进行选取，以p_u为例，选取方法如下：

这样对应p_u，子块内所有像素点都选取了相应的场图值，组成一个集合 P_u，余下的场图值也组成一个集合P_uc；

同样地，对应p_v，存在集合P_v和P_vc；

分别计算上述四个集合的矢量和的幅值，记为M_u、M_uc、M_v和M_vc；

取{M_u，M_uc}最大值，记为U_max，取{M_v，M_vc}最大值，记为V_max，将 U_max和V_max进行比较，舍弃较小数值对应的数据集合；假设保留的集合为P_u和P_uc，且M_u大于M_uc；根据保留的子块场图候选解对应的两组矢量和幅值计算得到一个加权因子C，定义如下：

则该子块所有像素点对应的场图值集合为P＝P_u*C；如果M_uc大于M_u，则P＝P_uc*C。

在一实施方式中，在步骤(3b)中，以每一个像素点的场图值p_n ^s为初始值，对图像内所有像素点的两个场图候选解(p₁，p₂)再次进行选取的方法类似于步骤(3a)中的式(8)，即采用如下式(10)进行选取：

在一实施方式中，在步骤(3c)中，N为像素点个数的8％- 10％。

在一实施方式中，所述步骤(4)根据以下执行：

(4a)场图合并

将上述多组场图进行合并得到一幅场图，具体实现方式是：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

所述重新置0是指重新设置为未知或未选取状态。

在一实施方式中，所述步骤(5)根据以下执行：

对合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

在一实施方式中，所述步骤(6)根据以下执行：

利用得到的最终场图，进行平滑处理，然后计算得到水图和另一成分(例如脂肪)的图。

在一子实施方式中，所述步骤(6)可根据以下执行：利用得到的最终场图，进行平滑处理，然后利用下式计算得到水图和另一成分(例如脂肪)的图：

其中ρ_w、ρ_f、A⁺和S如上所定义，根据ρ_w和ρ_f得到的水图和另一成分(例如脂肪)的图。

术语“场图值拟合误差曲线”或“phasor-error谱”，是指由两点水脂分离信号模型推导出了单个像素点和子块内所有像素点为一个整体的 phasor-error谱。

术语“场图的初始化”，是指根据子块phasor-error谱得到的场图值对该子块内所有像素的场图值进行了初始化。

在一实施方式中，当待分离成分为水和其他成分时，只需将式 (1)中的脂肪相对于水的化学位移f_F替换该其他成分对应的化学位移值，即可分离水图和该其他成分图。

此外，本发明还提供一种计算机可读介质，该计算机可读介质具有存储在其中的程序，该程序是计算机可执行的，以使计算机执行上述各种步骤的处理。

本发明还提供一种利用上述方法的磁共振化学位移编码成像装置。该装置包含了不同模块用于实现上述方法中提及的各个步骤。

本发明的振化学位移编码成像装置可以包括：场图候选解求解模块，其从二维原始磁共振图像求得该图像的两组场图候选解；场图候选解分割模块，其在预先定义的不同的分块尺度下，将在场图候选解求解模块中求得的前述两组场图候选解分别分割为不交叠的子块；局部迭代场图提取模块，其在场图候选解分割模块定义的每一个分块尺度下分别采用局部迭代场图提取算法处理在场图候选解分割模块中获得的场图候选解子块，得到相应的原始图像空间下的场图；场图合并模块，其将在局部迭代场图提取模块中各个分块尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；最终场图获取模块，其利用局部迭代场图提取算法对在场图合并模块中置0像素点的场图值进行选取，得到最终场图；图像分离模块，其对在最终场图获取模块中得到的最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图。在一子实施方式中，所述待分离成分为水和另一成分。

在一实施方式中，所述场图候选解求解模块还包括：描述各像素点场图值拟合误差曲线模块，其描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线；以及获取各像素点场图候选解模块，其根据所述场图值拟合误差曲线寻找使之达到局部极小值的所有场图值。

在一子实施方式中，所述场图候选解求解模块可包括：描述各像素点场图值拟合误差曲线模块，其根据权利要求3中所采用的式(5)描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线；以及获取各像素点场图候选解模块，其根据在描述各像素点场图值拟合误差曲线模块中获得的所述场图值拟合误差曲线寻找使式(5)达到局部极小值的所有场图值p，构成该像素点的场图候选解，个数为1个或2个，以符号p₁和p₂表示，其中场图值数目为1时，p₁＝p₂。

在一实施方式中，所述场图候选解分割模块根据预先定义的分块尺度值将二维原始磁共振图像进行分块，得到的每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度。在一实施方式中，所述分块尺度值为整数，选自1至16中的一个或多个数值。在另一实施方式中，所述分块尺度值为6、8、10和12。

在一实施方式中，所述局部迭代场图提取模块包括：场图初始化模块，其实现前面提及的本发明所述的磁共振化学位移编码成像方法中所述的(3a)场图的初始化步骤；场图更新模块，其实现前面提及的(3b)场图的更新步骤；以及迭代更新场图模块，其实现前面提及的(3c)迭代更新场图步骤。

在一实施方式中，所述场图合并模块通过以下实现方式将在局部迭代场图提取模块中获得的多组场图进行合并得到一幅场图：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

在一实施方式中，所述最终场图获取模块对在场图合并模块中合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

在一子实施方式中，所述图像分离模块对在最终场图获取模块中得到的最终场图进行平滑处理，然后(可选地根据前面提及的式(11))计算得到水图和另一成分的图。

前述任一实施方式或子实施方式中，所述另一成分为脂肪。

本发明与现有技术相比具有以下优点：

本发明提出的两点磁共振化学位移编码成像技术不对TE(回波时间)施加任何的限制，达到更快的图像采集速度，使其适用于一些快速成像的医学应用中，如动态成像和腹部屏气成像等；应用多尺度的局部迭代场图提取算法可独立得到多组场图，然后利用自检验机制将其进行合并，保证了场图估计的正确性。

附图说明

图1为本发明提出的基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法应用于两点水脂分离算法的流程图。

图2显示根据本发明所述方法对二维原始磁共振图像进行处理，最终分离得到水图和脂肪图；其中图中K＝6、8、10和12分别表示一个尺度，在每一个尺度下独立地完成区域迭代场图提取算法(RIPE)。

图3显示了本发明所述方法在体组织实验的颈部的分离结果图； (a)颈部冠状位分离结果；(b)颈部矢状位分离结果，图像从左向右依次为分离脂肪图和分离水图。

具体实施方式

为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚，下面结合附图和实施例对本发明作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅仅用于解释本发明，而非对本发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与本发明相关的部分而非全部内容。

通过以下获得本发明所述方法的式(5)：

含有水和脂肪的两种成分的MR信号模型为：

其中TE为回波时间，S_n是在回波时间TE_n下的信号强度；ρ_w和ρ_f是水和脂肪的信号强度值；水的化学位移设定为0，f_F是脂肪相对于水的化学位移，为-3.5ppm；f_B是静磁场B₀的场图值，i表示虚数。

利用式(2)定义参数c₁、c₂、q和p如下：

其中ΔTE＝TE₁-TE₂，q代表TE₁时间内由于f_B造成的水的相位，p代表ΔTE时间内由于f_B造成的水的相位累积，本发明将p称之为相位矢量 (phasor)，其角度与f_B成正比，范围是[-ππ]。在本发明中，在实际应用中，只要估计出p，无需得到实际的f_B，就可得到分离的图像，因此本发明将p也称之为场图。

利用式(2)的p，由式(1)和(2)可解得q：

通过利用式(2)、(3)，将式(1)写作矩阵的形式：

其中是4×1矩阵， 其中S₁、S₂分别独立地为在回波时间TE₁、TE₂下一个像素点的信号强度，Re代表取复数实部，Im代表取复数虚部。

通过式(4)，根据变量投影方法(variable projection， VARPRO)，使用最小二乘法可解得p：

p＝argminerr(p)＝argmin||(I-AA⁺)S||² (5)。

上式(5)可以改写为以下形式：

I是4×4单位矩阵；A⁺＝(A^TA)-¹A^T，上标“+”表示取矩阵伪逆，上标“T”代表矩阵的转置，表示点乘。err(p)描述了场图p和拟合误差之间的关系，表示为场图值拟合误差曲线(phasor-error谱)，其定义为不同场图 (phasor)值p下对应的模型拟合误差。

公式(6)用于描述感兴趣区域内所有像素点作为一个整体的 phasor-error谱。公式(6)达到局部极小值时的p值，数目为1个或2个。

对于描述感兴趣区域(volume of interest，VOI)内所有像素点作为一个整体的phasor-error谱，式(6)依然成立，此时S为：

其中和分别代表一组像素点的信号强度，和都是行向量，定义如式(4)中表述，行向量中每一个元素代表一个像素点的信号强度。

本发明的基于多尺度局部迭代场图提取算法的磁共振化学位移编码成像方法，实现了对场图值的正确选取和待分离成分的分离。流程图见附图1。

参考图1和图2，本发明的所述方法具体实施如下：

(1)求得二维原始磁共振图像的两组场图候选解，根据以下执行：

1)场图候选解的求解：

根据式(5)描述二维原始磁共振图像中每个像素点的phasor-error谱，根据该phasor-error谱寻找使式(5)达到局部极小值的所有场图值，构成该像素点的场图候选解，个数为1个或2个，以符号p₁和p₂表示，在场图值数目为 1时，p₁＝p₂。

(2)在不同的尺度下，将前述两组场图候选解分别分割为不交叠的子块：

1)图像的分块：

本发明将二维原始图像根据预先定义的分块尺度值将图像进行分块，比如选取分块尺度6、8、10和12可得到四组分块的图像，每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度，如当分块尺度是8 时，图像中每一个子块的大小均为8×8。

所述分块尺度值可以选自1至16中的一个或多个数值。例如，所述分块尺度值是6、8、10和12。

(3)在每一个尺度下分别采用局部迭代场图提取算法，得到相应的原始图像空间下的场图：

(3a)场图的初始化：

每一个尺度下的每一子块内所有像素点作为整体根据式(6)和(7)描述 phasor-error谱，根据该phasor-error谱寻找使式(6)达到局部极小值的场图值，所述场图值为该子块的场图候选解，并设置以下阈值排除一定数目的场图值，使之数目为1个或2个，阈值设定为A_min+α*(A_max-A_min)，其中A_max和A_min分别表示phasor-error谱的最大值和最小值，α的选取要排除拟合误差大的场图值，α可以取0.4。该子块的场图候选解以符号p_u和p_v表示，当场图值数目为1时，p_u＝p_v。

以上述子块的场图值p_u和p_v为初始值，分别对该子块内所有像素点的两个场图候选解(p₁，p₂)进行选取，以p_u为例，选取方法如下：

这样对应p_u，子块内所有像素点都选取了相应的场图值，组成一个集合P_u，余下的场图值也组成一个集合P_uc；同样地，对应p_v，存在集合P_v和P_vc。

分别计算上述四个集合的矢量和的幅值，记为M_u、M_uc、M_v和 M_vc。取{M_u，M_uc}最大值，记为U_max，取{M_v，M_vc}最大值，记为V_max，将 U_max和V_max进行比较，舍弃较小数值对应的数据，比如U_max大于V_max，则舍弃集合P_v和P_vc，保留P_u和P_uc，反之，保留P_v和P_vc。为方便表述，假设保留的集合为P_u和P_uc，且M_u大于M_uc。根据保留的子块场图候选解对应的两组矢量和幅值计算得到一个加权因子C，定义如下：

则该子块所有像素点对应的场图值集合为P＝P_u*C；如果M_uc大于M_u，则P＝P_uc*C。

对每一个尺度下每一个子块均进行上述场图初始化，将每一个尺度下所有子块的场图值集合P合并在一起，即得到该尺度下的初始场图值集合P₀。

(3b)场图的更新：

对上述初始场图值集合P₀进行低通滤波，且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理。假定迭代后得到的场图值集合为P_n，则滤波后得到P_n ^s，然后以每一个像素点的场图值p_n ^s为初始值，对图像内所有像素点的两个场图候选解(p₁，p₂)再次进行选取，选取方法类似于式 (8)，如下式(10)所示：

并记录场图值发生改变的像素点的个数N。

(3c)迭代更新场图：

重复步骤(3b)，直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止，得到一组场图。N可以为像素点个数的8％-10％。

对应于多组分块尺度，经过步骤4)和5)可得到多组的场图。

(4)将各个尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0：

(4a)场图合并：

将上述多组场图进行合并得到一幅场图，具体实现方式是：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取，得到最终场图：

对合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

(6)对最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图：

水、脂肪图重建：

利用得到的最终场图，进行平滑处理，然后利用下式计算得到水图和脂肪图：

其中ρ_w、ρ_f、A⁺和S如上所定义，根据ρ_w和ρ_f得到的水图和脂肪图。

本发明的多尺度的局部迭代场图提取算法(multi-scale RIPE)，是在多个尺度下分别进行场图的估计，各尺度之间是相互独立的，提供了一种自校验机制，保证场图估计的正确性。

体组织实验

采用本发明方法在体组织实验进行了测试，其中磁共振扫描操作采用本领域已知的磁共振扫描标准流程，对象为人体颈部，磁共振扫描采用仪器为磁共振成像系统(MR)(联影，上海)。采集序列为TSE，B₀＝1.5T，成像参数为：矩阵大小＝304×320，TE＝[0 1.7]ms。实验数据的处理环境是具有Intel E5-2650v2 CPU和64GB RAM的工作站，采用的数据处理软件是 MATLAB。

参考图3。图3为颈部的水脂分离结果。(a)颈部冠状位分离结果；(b)颈部矢状位分离结果，图像从左向右依次为分离脂肪图和分离水图。可以看出在体组织实验中，得到的水图和脂肪图中无明显分离错误的组织。

本发明实施例所述的方法不仅可以用于水脂分离成像，也可以用于分离其他化学位移成份的应用中，区别仅仅在于在式(1)中将脂肪相对于水的化学位移值f_F替换为待分离的其他化学位移成份对应的化学位移值。

本发明不受上面所显示和描述的实施例的限制，而是在权利要求的范围内可以改变、置换以及等同替代。

参考文献

[1]Dixon WT.Simple proton spectroscopic imaging.Radiology 1984；153(1):189-194.

[2]Xiang QS.Two-point water-fat imaging with partially-opposed phase(POP)acquisition:An asymmetric Dixon method.Magnetic resonance in medicine2006；56(3):572-584.

[3]Eggers H,Brendel B,Duijndam A,Herigault G.Dual-echo Dixon imaging withflexible choice of echo times.Magnetic resonance in medicine 2011；65(1):96-107。

Claims (35)

1.磁共振化学位移编码成像方法，所述方法包括：

(1)求得二维原始磁共振图像的两组场图候选解；

(2)在不同的分块尺度下，将前述两组场图候选解分别分割为不交叠的子块；

(3)在每一个分块尺度下分别采用局部迭代场图提取算法处理场图候选解子块，得到相应的原始图像空间下的场图；

(4)将各个分块尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；

(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取，得到最终场图；

(6)对最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图。

2.根据权利要求1所述的磁共振化学位移编码成像方法，其特征在于，所述待分离成分为水和另一成分。

3.根据权利要求2所述的磁共振化学位移编码成像方法，其特征在于，步骤(1)根据以下执行：

描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线，并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值，构成该像素点的两组场图候选解。

4.根据权利要求3所述的磁共振化学位移编码成像方法，其特征在于，步骤(2)根据以下执行：

根据预先定义的分块尺度值将二维原始磁共振图像进行分块，得到的每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度。

5.根据权利要求4所述的磁共振化学位移编码成像方法，其特征在于，所述分块尺度值为整数，选自1至16中的一个或多个数值。

6.根据权利要求5所述的磁共振化学位移编码成像方法，其特征在于，所述分块尺度值为6、8、10和12。

7.根据权利要求4所述的磁共振化学位移编码成像方法，其特征在于，步骤(3)根据以下执行：

(3a)场图的初始化

对每一个分块尺度下的每一个子块，求得子块场图候选解，利用每一个子块场图候选解为初始值，对子块内所有像素点的所述两组场图候选解进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小，符合所述规则的像素点场图候选解组成一个场图值集合，未被选取的像素点场图候选解组成另一个场图值集合；计算每一个场图候选解对应场图值集合的矢量和；比较所述矢量和的幅值，选取最大幅值对应的场图值集合；

对每一个分块尺度下每一个子块均进行上述场图初始化，将每一个分块尺度下所有子块的场图值集合合并在一起，即得到该分块尺度下的初始场图值集合；

(3b)场图的更新

对上述初始场图值集合进行低通滤波，且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理；然后以每一个像素点的场图值为初始值，对图像内所有像素点的两个场图候选解再次进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小；

并记录场图值发生改变的像素点的个数N；

(3c)迭代更新场图

重复步骤(3b)，直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止，得到一组场图；

对应于多组分块尺度，经过步骤(3a)、(3b)和(3c)得到多组的场图。

8.根据权利要求7所述的磁共振化学位移编码成像方法，其特征在于，在步骤(3c)中，N为像素点个数的8％-10％。

9.根据权利要求7所述的磁共振化学位移编码成像方法，其特征在于，步骤(4)根据以下执行：

(4a)场图合并

将所述多组场图进行合并得到一幅场图，具体实现方式是：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

10.根据权利要求9所述的磁共振化学位移编码成像方法，其特征在于，步骤(5)根据以下执行：

对合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

11.根据权利要求10所述的磁共振化学位移编码成像方法，其特征在于，步骤(6)根据以下执行：

利用得到的最终场图，进行平滑处理，然后计算得到水图和另一成分的图。

12.根据权利要求2-11中任一项所述的磁共振化学位移编码成像方法，其特征在于，所述另一成分为脂肪。

13.一种磁共振化学位移编码成像装置，其特征在于，其包括：

场图候选解求解模块，其从二维原始磁共振图像求得该图像的两组场图候选解；

场图候选解分割模块，其在预先定义的不同的分块尺度下，将在场图候选解求解模块中求得的前述两组场图候选解分别分割为不交叠的子块；

局部迭代场图提取模块，其在场图候选解分割模块定义的每一个分块尺度下分别采用局部迭代场图提取算法处理在场图候选解分割模块中获得的场图候选解子块，得到相应的原始图像空间下的场图；

场图合并模块，其将在局部迭代场图提取模块中各个分块尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；

最终场图获取模块，其利用局部迭代场图提取算法对在场图合并模块中置0像素点的场图值进行选取，得到最终场图；

图像分离模块，其对在最终场图获取模块中得到的最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图。

14.根据权利要求13所述的磁共振化学位移编码成像装置，其特征在于，所述待分离成分为水和另一成分。

15.根据权利要求14所述的磁共振化学位移编码成像装置，其特征在于，所述场图候选解求解模块包括：

描述各像素点场图值拟合误差曲线模块，其描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线；以及

获取各像素点场图候选解模块，其根据所述场图值拟合误差曲线寻找使之达到局部极小值的所有场图值。

16.根据权利要求15所述的磁共振化学位移编码成像装置，其特征在于，所述场图候选解分割模块根据预先定义的分块尺度值将二维原始磁共振图像进行分块，得到的每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度。

17.根据权利要求16所述的磁共振化学位移编码成像装置，其特征在于，所述分块尺度值为整数，选自1至16中的一个或多个数值。

18.根据权利要求17所述的磁共振化学位移编码成像装置，其特征在于，所述分块尺度值为6、8、10和12。

19.根据权利要求16所述的磁共振化学位移编码成像装置，其特征在于，所述局部迭代场图提取模块包括：

场图初始化模块，其实现权利要求7中所述的(3a)场图的初始化步骤；

场图更新模块，其实现权利要求7中所述的(3b)场图的更新步骤；以及

迭代更新场图模块，其实现权利要求7中所述的(3c)迭代更新场图步骤。

20.根据权利要求19所述的磁共振化学位移编码成像装置，其特征在于，所述场图合并模块通过以下实现方式将在局部迭代场图提取模块中获得的多组场图进行合并得到一幅场图：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

21.根据权利要求20所述的磁共振化学位移编码成像装置，其特征在于，所述最终场图获取模块对在场图合并模块中合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

22.根据权利要求21所述的磁共振化学位移编码成像装置，其特征在于，所述图像分离模块对在最终场图获取模块中得到的最终场图进行平滑处理，然后计算得到水图和另一成分的图。

23.根据权利要求14-22中任一项所述的磁共振化学位移编码成像装置，其特征在于，所述另一成分为脂肪。

24.一种计算机可读介质，该计算机可读介质具有存储在其中的程序，该程序是计算机可执行的以使计算机执行包括以下步骤的处理：

(1)求得二维原始磁共振图像的两组场图候选解；

(2)在不同的分块尺度下，将前述两组场图候选解分别分割为不交叠的子块；

(3)在每一个分块尺度下分别采用局部迭代场图提取算法处理场图候选解子块，得到相应的原始图像空间下的场图；

(4)将各个分块尺度下得到的场图合并，保留具有相同场图值的像素点的场图值，其余像素点的场图值置0；

(5)利用局部迭代场图提取算法对置0像素点的场图值进行选取，得到最终场图；

(6)对最终场图进行平滑处理，分离原始磁共振图像，得到准确的各待分离成分的图。

25.根据权利要求24所述的计算机可读介质，其特征在于，所述待分离成分为水和另一成分。

26.根据权利要求25所述的计算机可读介质，其特征在于，步骤(1)根据以下执行：

描述二维原始磁共振图像中每个像素点的场图值拟合误差曲线，并根据该场图值拟合误差曲线寻找使之达到局部极小值的所有场图值，构成该像素点的两组场图候选解。

27.根据权利要求26所述的计算机可读介质，其特征在于，步骤(2)根据以下执行：

根据预先定义的分块尺度值将二维原始磁共振图像进行分块，得到的每一组图像中均包含数目不等的子块，每一个子块的大小对应其分块尺度。

28.根据权利要求27所述的计算机可读介质，其特征在于，所述分块尺度值为整数，选自1至16中的一个或多个数值。

29.根据权利要求28所述的计算机可读介质，其特征在于，所述分块尺度值为6、8、10和12。

30.根据权利要求27所述的计算机可读介质，其特征在于，步骤(3)根据以下执行：

(3a)场图的初始化

对每一个分块尺度下的每一个子块，求得子块场图候选解，利用每一个子块场图候选解为初始值，对子块内所有像素点的所述两组场图候选解进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小，符合所述规则的像素点场图候选解组成一个场图值集合，未被选取的像素点场图候选解组成另一个场图值集合；计算每一个场图候选解对应场图值集合的矢量和；比较所述矢量和的幅值，选取最大幅值对应的场图值集合；

对每一个分块尺度下每一个子块均进行上述场图初始化，将每一个分块尺度下所有子块的场图值集合合并在一起，即得到该分块尺度下的初始场图值集合；

(3b)场图的更新

对上述初始场图值集合进行低通滤波，且对于下面每一步迭代后得到的场图值集合均进行低通滤波处理；然后以每一个像素点的场图值为初始值，对图像内所有像素点的两组场图候选解再次进行选取，选取规则为像素点场图候选解与初始值差值的幅值距离最小；

并记录场图值发生改变的像素点的个数N；

(3c)迭代更新场图

重复步骤(3b)，直至N小于一定的数值或在数次迭代中均固定在某一数值则迭代停止，得到一组场图；

对应于多组分块尺度，经过步骤(3a)、(3b)和(3c)得到多组的场图。

31.根据权利要求30所述的计算机可读介质，其特征在于，其特征在于，在步骤(3c)中，N为像素点个数的8％-10％。

32.根据权利要求30所述的计算机可读介质，其特征在于，步骤(4)根据以下执行：

(4a)场图合并

将所述多组场图进行合并得到一幅场图，具体实现方式是：如果像素点在所有的场图中具有相同的数值，则该保留该像素点的场图值；如果像素点在所有的场图中具有不同的数值，则该像素点的场图值重新置0。

33.根据权利要求32所述的计算机可读介质，其特征在于，步骤(5)根据以下执行：

对合并的场图中置0像素点的场图值重复步骤(3b)和(3c)，得到最终的场图。

34.根据权利要求33所述的计算机可读介质，其特征在于，步骤(6)根据以下执行：

利用得到的最终场图，进行平滑处理，然后计算得到水图和另一成分的图。

35.根据权利要求25-34中任一项所述的计算机可读介质，其特征在于，所述另一成分为脂肪。