本发明涉及磁共振成像中的水脂技术，特别是涉及一种磁共振成像的图像水脂分离方法和系统。
背景技术：
在常规的临床磁共振成像(Magnetic Resonance Imaging，MRI)中，脂肪通常呈现高信号，将会掩盖水激发所产生的信号，而一般认为组织中水激发所产生的信号包含了大部分的病理学信息，因此脂肪的高信号会使得某些病变无法被显示，甚至造成临床医生对病变的误判，因此脂肪信号的存在严重降低了MR信号的临床病理学研究价值。现在常用的方法是多回波Dixon(或其变种、扩展)水脂分离方法，它利用水和脂肪之间存在的频率化学位移特性，将水和脂肪信号从MR图像中分离出来。但是Dixon水脂分离算法在处理水脂含量差异较大的组织时，面临一个挑战：在这类组织中，Dixon算法估计的主磁场B0不均匀性，即场图值(单位Hz)fB具有两个解值，其中一个是正确值fB,true(通常叫真实解)，而另外一个是水脂分反的值fB,swap(通常叫反转解)，这两个解不能够在单个像素的基础上进行正确的区分，需用到其他一些额外的辅助信息。
最早的Dixon方法假设静磁场B0均匀，分别采集水脂同向(相位差为0)和异向(相位差为π)两个时刻的图像，通过简单的加/减计算得到水和脂肪的分离图像。后来，出现了一种迭代水脂分离的方法(Iteration Decomposition of water and fat with Echo Asymmetry and Least-squares,IDEAL)，IDEAL水脂分离方法最大的缺点是该方法利用迭代算法估计静磁场B0的不均匀性，将不可避免的会收敛到一个错误的场图值fB,swap，造成水脂分反的后果。为了克服这个缺点，不同的科研工作者提出了多种策略进行场图的估计，比如首先对图像进行采样，得到低分辨率的原始图像，然后利用IDEAL算法估计得到场图值，从中选出场图值与中值最接近的15个像素，再在这一组像素中选取与图像重心几何位置最接近的一个像素作为种子点，以此种子点作为其所对应的一簇高分辨率像素的初始值，利用局部增长算法和IDEAL算法对高分辨率图像中未估计像素进行场图值的估计；或者是求得每一个像素的phasor(即相位复矢量)，即b＝exp(i2πγfBΔt)，其中Δt为回波时间间隔，如前表述，由于fB存在两个解，因此b也是两个值btrue和bswap，然后选取信噪比高并表现为水脂比例接近的像素作为合理的种子点，然后利用局部增长算法对未估计像素正确的b值进行选取。此外还有一种使用多分辨率的IDEAL算法，首先将图像作为一个整体，估计其phasor，然后将图像分成多个具有交叠区域的子图像，利用前述phasor作为初始值，利用IDEAL算法对子图像进行场图phasor估计，以此类推，将子图像再次进行细分、估计phasor，直至得到满意的水脂分离结果。
基于前述各种算法均存在其缺点，比如，只在低分辨率的场图图像上选取了一个像素作为种子点，这单个种子点对应于高分辨率图像中的一簇像素，如果原始MR图像中的组织是空间分离的，则以这一簇像素作为种子点是不可能增长到相分离组织之中的，即这些组织中的场图值不能够得到正确的估计；又如后两种方法中使用了不同的策略去完成水脂分离，但这两种方法只是得到了场图的phasor，而没有得到真正的场图值fB。
基于上述现有技术的缺点，有待进一步地提高。
技术实现要素：
基于此，有必要针对现有技术的问题，提供一种磁共振成像的图像水脂分离方法和系统，通过本发明的方法和系统可找到适当数目的种子点，以局部增长的方法完成对每个像素的场图值两个解(即真实解fB,true和反转解fB,swap)的正确选取，实现水脂分离。
本发明提供了一种磁共振成像的图像水脂分离方法，其包括：
基于磁共振成像获得的原始图像，估计所述原始图像中每个像素的场图值，获得所述每个像素的初始场图值，生成相应的场图；
将所述场图进行分块处理，获得多个子块，生成每一个子块中初始场图值的分布曲线；
基于所述每一个子块中初始场图值的分布曲线，选择所述初始场图值的分布曲线为三段式结构对应的子块，作为种子块；
取所述种子块内初始场图值的分布曲线中的中间段数值作为所述种子块内相应像素的场图值最终解，以所述种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计所述每个像素的场图值最终解；
根据所述每个像素的场图值最终解，获得水的图像和脂肪的图像。
在其中一个实施例中，所述以所述种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计所述每个像素的场图值最终解的过程包括低分辨率的增长步骤，该步骤为：
将所述初始场图值的分布曲线为两段式结构的子块作为待估子块，以所述种子块内每一个像素的场图值最终解为初始值，对所述待估子块进行局部增长处理，估计获得所述待估子块内像素的场图值最终解。
在其中一个实施例中，所述对所述待估子块进行局部增长处理、估计获得所述待估子块内像素的场图值最终解的步骤为：
步骤a、从所述场图中与所述种子块相邻的子块中选择所述待估子块；
步骤b、分别计算所述待估子块内初始场图值的分布曲线中两段数值与所述种子块内相应像素的场图值最终解之间的欧氏距离，选取所述欧氏距离相对较小的一段数值作为待估子块内相应像素的场图值最终解；
步骤c、将所述待估子块和所述种子块进行合并，将合并后获得的子块作为所述步骤a中的种子块，重复执行步骤a至步骤c，直至对全部待估子块完成局部增长处理，获得相应像素的场图值最终解。
在其中一个实施例中，所述以所述种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计所述每个像素的场图值最终解的过程还包括高分辨率的增长步骤，该步骤为：
将所述种子块和待估子块中的像素作为初始起点，对所述原始图像中除所述初始起点之外的剩余像素进行局部增长处理，估计获得所述原始图像中剩余像素对应的场图值最终解。
在其中一个实施例中，所述对所述原始图像中剩余像素进行局部增长处理、估计获得所述原始图像中剩余像素对应的场图值最终解的步骤包括：
分别计算每个剩余像素的所述初始场图值和与该剩余像素相邻的所述初始起点的场图值最终解之间的欧氏距离，选取所述欧氏距离最小的对应的所述初始场图值作为所述原始图像中该剩余像素的场图值最终解。
在其中一个实施例中，所述获得所述每个像素的初始场图值生成相应的场图的步骤包括：
针对所述每个像素判断所述初始场图值的个数，当所述初始场图值为两个时，该像素的初始场图值包括第一解和第二解；当所述初始场图值为一个时，将该像素对应的第一解和第二解均赋值为所述初始场图值；
汇总所述原始图像中每个像素对应的第一解和第二解，分别生成所述第一解对应的第一场图和所述第二解对应的第二场图。
在其中一个实施例中，所述将所述场图进行分块处理、获得多个子块及生成每一个子块中初始场图值的分布曲线的步骤包括：
分别将所述第一场图和第二场图进行分块处理，获得所述多个子块；
针对每一个子块，合并统计该子块中相应像素的所述第一解和第二解，获得所述每个子块内所有的初始场图值；
将所述每一个子块内所有的初始场图值进行排序，获得场图值排序位置；
绘制表征初始场图值与场图值排序位置的对应关系曲线，生成所述每一个子块中初始场图值的分布曲线。
在其中一个实施例中，所述对所述原始图像中剩余像素进行局部增长处理、估计获得所述原始图像中剩余像素对应的场图值最终解的步骤中，分别计算每一个剩余像素对应的场图值的第一解和第二解与所述初始起点的场图值最终解之间的欧氏距离，选取所述场图值的第一解和第二解之中与所述初始起点的场图值最终解之间的欧氏距离最小的解，作为所述原始图像中该剩余像素的场图值最终解。
在其中一个实施例中，所述原始图像采用至少四个回波时间对应的回波序列信号进行图像重建。
基于上述方法，本发明还提供了一种磁共振成像的图像水脂分离系统，其包括：
场图获取模块，用于基于磁共振成像获得的原始图像，估计所述原始图像中每个像素的场图值，获得所述每个像素的初始场图值，生成场图；
曲线绘制模块，用于将所述场图进行分块处理，获得多个子块，生成每一个子块中初始场图值的分布曲线；
种子选择模块，用于基于所述每一个子块中初始场图值的分布曲线，选择所述初始场图值的分布曲线为三段式结构对应的子块，作为种子块；
局部增长处理模块，用于取所述种子块内初始场图值的分布曲线中的中间段数值作为所述种子块内像素的场图值最终解，以所述种子块内像素的场图值最终解为初始值、利用局部增长方法估计所述所述每个像素的场图值最终解；及
图像分离模块，用于根据所述每个像素的场图值最终解，获得水的图像和脂肪的图像。
本发明提出了一种基于多分辨率的场图值局部增长策略，实现对场图值的正确选取，达到正确水脂分离的目的，其首先完成对图像中每一个像素的场图值估计，基于初始估计的场图值进行分块处理，绘制每一个子块中场图值的分布曲线，利用其分布规律找到可作为种子点的子图像，利用该组子图像采用局部增长策略完成对其他子图像场图值的选取以及完成对余下像素的场图选取。
附图说明
图1为本发明的方法流程示意图；
图2(A)和图2(B)为本发明初始场图值的分布曲线图；
图3为本发明的方法的优选实施例流程示意图；
图4为本发明的系统结构示意图；
图5为本发明的系统的优选实施例结构示意图；
图6为利用本发明方法的实验效果示意图。
具体实施方式
基于磁共振成像中的水脂分离技术，本发明为了提供一种磁共振成像的图像水脂分离方法和系统，其提出一种基于多分辨率的场图值局部增长策略，实现对场图值的正确选取，达到正确水脂分离的目的，其首先由离散化场图值对应的模型拟合误差曲线估计得到图像中每一个像素的场图值估计，基于初始估计的场图值进行分块处理，绘制每一个子块中场图值的分布曲线，利用其分布规律找到可作为种子点的子块，利用该组子块采用局部增长策略完成对其他子块场图值的选取以及完成对余下像素的场图选取。本发明提到的场图概念基于IDEAL算法，其是指假设组织信号由多种不同成分被同时激发得到，并已知各成分相对于水的化学位移，然后构造组织信号与各成分和静磁场B0不均匀性之间关于时间的数学模型，本发明基于这一模型量化构造组织信号与各成分和静磁场B0不均匀性形成的场图，静磁场B0均匀时，场图值为0；静磁场B0不均匀效应越大，场图值越大。以下结合附图详细说明本发明的具体实施例。
如图1所示，本发明提供的一种磁共振成像的图像水脂分离方法，其包括以下步骤：
步骤100，基于磁共振成像获得的多个不同回波时间下的原始图像，估计原始图像中每个像素的场图值，获得每个像素的初始场图值，生成相应的场图；这里的场图可通过矩阵数据表示，对应于下述提到的真实解F1和反转解F2，可以简单理解为按照像素的位置对应于初始场图值而构成的矩阵数据。
步骤200，基于每个子块中的待估场图值是平滑连续的假设，将场图进行分块处理，获得多个子块，生成每一个子块中初始场图值的分布曲线。
步骤300，基于每一个子块中初始场图值的分布曲线，选择初始场图值的分布曲线为三段式结构对应的子块，作为种子块；这里提到每个子块对应的初始场图值的分布曲线是指表征场图值与场图值排序位置的对应关系曲线，即将子块中的初始场图值由小至大排序，绘制的纵坐标代表场图值、横坐标代表相应场图值对应的排序位置的分布曲线，或者绘制的横坐标代表场图值、纵坐标代表相应场图值对应的排序位置的分布曲线。
步骤400，取种子块内初始场图值的分布曲线中的中间段数值作为种子块内相应像素的场图值最终解，以种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计上述原始图像中每个像素的场图值最终解。
步骤500，根据上述原始图像中每个像素的场图值最终解，获得水的图像和脂肪的图像。
上述步骤100基于水和脂肪化学位移的MR信号模型来获得离散化场图值对应的模型拟合误差曲线，基于此场图-模型拟合误差曲线通过数值最优化问题来估计原始图像中每一个像素的场图值，获得每个像素的初始场图值。以下将具体说明上述步骤100的一种具体实现方法。
通常将含有两种成分(常见为水和脂肪)的简化MR信号模型为：
其中，N为回波个数；ρw和ρf是水和脂肪的强度值，其为复数；fF是脂肪相对于水的化学位移，为-3.5ppm，在3T系统中，fF＝-431Hz；fB是静磁场B0的场图值；tn表示回波时间。
如果组织中水和脂肪的含量差别比较大，则上式简化为下述公式(2)(以纯水为例，ρf＝0)：
从上式可以看出，在纯水的情况下场图值有两个解：真实解fB(即fB,true)和水值ρw；反转解(即fB,swap)和脂肪值同样的，对于纯脂肪组织也会有两个解：真实解fB(即fB,true)和脂肪值ρf；反转解(即fB,swap)和水值而对于水脂比例接近的组织，只有一组解：场图值fB、水值ρw和脂肪值ρf。在纯水或纯脂肪的情况下，真实解和反转解相差fF，又由于fF是负值，因此在纯水的组织中，反转解比真实解大fF，而在纯脂肪组织中，反转解比真实解小fF。
寻找场图值使下述公式(3)E(fb,q)达到最小，即为场图的估计值：
公式(3)
其中，E(fB)表示模型拟合误差，I表示N×N的单位矩阵；上标“+”代表M-P广义逆，diag(X)表示利用向量X的元素构造对角矩阵，S＝[S(t1),S(t2),…,S(tN)]H，A为N×2的矩阵(见下述公式(4))，||||2表示取L2范数。
基于上述各个公式，在本发明的一个实施例中，上述步骤100可以包括以下步骤：
首先，针对至少四个回波时间对应的回波序列进行图像重建，获得上述原始图像；在本发明的一个实施例中优选采用4个回波时间对应的回波序列，比如回波时间表示为{TE1TE2TE3TE4}；
然后，选定一个合理的场图值分布区间进行离散化处理，获得一组离散化的离散场图值；
其次，针对至少四个回波时间进行排列组合获得至少两组回波时间组合，依照上述公式(3)可以获得每个像素在每组回波时间组合对应的回波序列中离散场图值与模型拟合误差之间的对应关系，生成离散场图值对应的模型拟合误差曲线，每个像素在至少两组回波时间组合对应的回波序列中相应生成至少两条模型拟合误差曲线；
最后，寻找使模型拟合误差曲线达到极小的离散场图值，即为每个像素对应的初始场图值。
本实施例在进行离散化处理时，优选利用{TE1 TE2 TE4}(等间距)和{TE1TE2TE3TE4}组合分别依照上述公式(3)得到两条误差曲线，其中TE3为TE2和TE4的黄金分割点，即基于上述两组回波时间组合得到的场图-模型拟合误差曲线，在纯水或纯脂肪情况下，像素对应可以获得最佳估计值为2个，分别为表示为场图值的真实解F1和反转解F2，即上述fB,true和fB,swap；在水脂含量接近的情况下，像素对应的最佳场图值为1个，可认为真实解F1和反转解F2是同一个数值。
基于上述实施例，如果基于步骤100获得的部分像素对应的场图值分为真实解F1和反转解F2，分别基于真实解F1和反转解F2会对应生成第一场图和第二场图，基于每个子块中所包含像素的待估场图值是平滑连续的这一假设，将第一场图和第二场图分别进行分块处理，每个场图对应获得多个子块。于是，在本发明的一个实施例中，上述步骤100中获得上述每个像素的初始场图值生成相应的场图的步骤包括：
首先，针对上述每个像素判断上述初始场图值的个数，当上述初始场图值为两个时，该像素的初始场图值包括第一解和第二解；当上述初始场图值为一个时，将该像素对应的第一解和第二解均赋值为上述初始场图值。
例如，当上述初始场图值包括真实解和反转解时，则第一解和第二解分别为真实解和反转解；如果当上述初始场图值只包括一个真实解时，则第一解和第二解均为该真实解。
然后，汇总上述原始图像中每个像素对应的第一解和第二解，分别生成上述第一解对应的第一场图和上述第二解对应的第二场图，于是针对原始图像中的所有像素则会出现两幅场图。
当原始图像中其中一个或多个像素对应的初始场图值包括场图值的真实解和反转解时，将上述场图值的真实解和反转解对应的场图中相应的子块进行合并后，获得每一个子块中初始场图值的分布曲线，从而提高计算的精确度。因此，在本发明的一个优选实施例中，上述步骤200中将上述场图进行分块处理、获得多个子块及生成每一个子块中初始场图值的分布曲线的步骤包括：
首先，分别将上述第一场图和第二场图进行分块处理，获得上述多个子块；
然后，针对每一个子块，合并统计该子块中相应像素的上述第一解和第二解，获得上述每个子块内所有的初始场图值；
其次，将上述每一个子块内所有的初始场图值进行排序，获得场图值排序位置；
最后，绘制表征初始场图值与场图值排序位置的对应关系曲线，生成上述每一个子块中初始场图值的分布曲线。
基于上述各个实施例，获得每一个子块中初始场图值的分布曲线时，对每一个子块单独处理，例如，上述第一场图和第二场图分别对应N个子块，N取从1到n的自然数，那么针对第N个子块，则将真实解F1和反转解F2分别对应的第一场图和第二场图中各自的第N个子块取出来，将真实解F1和反转解F2合并统计统一归入该子块中相应像素的初始场图值，并由小至大排序，绘制初始场图值的分布曲线，纵坐标代表场图值，横坐标代表将子块中的初始场图值由小至大排序后相应场图值对应的排序位置，即为1,2,3，…，K，其中K代表初始场图值的个数。具体如图2所示，绘制的初始场图值的分布曲线，结构形式大概分为二种样式，第一种样式见图2(A)所示的两段式结构，表示该子块为纯水或纯脂肪，曲线分为两段，有近似平行的两段曲线；第二种样式见图2(B)所示的三段式结构，表示该子块为水脂混合子块，曲线分为三段，存在近似平行的三段曲线。根据前述先验知识，三段曲线中的中间段的数值即为真实场图值，第一段数值代表纯脂肪组织的反转值，第三段数值代表纯水组织的反转值。基于此初始场图值的分布曲线的结构形式不同，本实施例选择初始场图值的分布曲线为三段式结构(即图2(B)所示分布曲线形式)对应的子块，作为种子块，取种子块对应的初始场图值的分布曲线中的中间段数值作为种子块内相应像素的场图值最终解。当然也可能存在其他结构的分布曲线，但是对于其他结构的分布曲线，采用后文提到的高分辨率的增长步骤进行处理。
上述步骤400中，以种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计每个像素的场图值最终解，这一过程可以采用低分辨率的局部增长方法。而在本发明的一个实施例中，如图3所示，上述步骤400中执行完步骤401取种子块内初始场图值的分布曲线中的中间段数值作为种子块内相应像素的场图值最终解之后，以种子块的场图值内每一个像素的最终解为初始值、利用局部增长方法估计上述原始图像中每个像素的场图值最终解的过程可以包括如下低分辨率的增长步骤402，该步骤为：
将初始场图值的分布曲线为两段式结构的子块作为待估子块，以种子块内每一个像素的场图值最终解为初始值，对待估子块进行局部增长处理，估计获得待估子块内像素的场图值最终解。这里的对待估子块进行局部增长处理、估计获得待估子块内像素的场图值最终解的步骤优选为：分别计算待估子块初始场图值的分布曲线中两段数值(是指图2(A)中近似平行的两段曲线对应的初始场图值)与种子块内相应像素的场图值最终解之间的欧氏距离，选取欧氏距离相对较小的一段数值作为待估子块内相应像素的场图值最终解。局部增长处理法利用待估子块周围的已估子块的信息，由于待估子块已知其只有两段数值分布，因此分别计算这两段数值与周围已估子块的数值之间的欧氏距离d，则具有较小d值的数值段即认为是真实场图值。当然本发明不限于只采用欧氏距离的计算来估计获得待估子块的场图值最终解，还可以采用其他局部增长处理的方式。
为提高本发明方法的处理精确度，使场图值的最终解更接近真实解，则本发明的一个实施例中，在上述低分辨率的增长步骤中优先对与种子块相邻的待估子块进行局部增长处理，估计获得相应像素的场图值最终解，其具体步骤如下所述：
步骤a、从上述场图中与上述种子块相邻的子块中选择上述待估子块；
步骤b、分别计算上述待估子块内初始场图值的分布曲线中两段数值与上述种子块内相应像素的场图值最终解之间的欧氏距离，选取上述欧氏距离相对较小的一段数值作为待估子块内相应像素的场图值最终解；
步骤c、将上述待估子块和上述种子块合并构成所述步骤a中的种子块，重复执行步骤a至步骤c，直至对全部待估子块完成局部增长处理，获得相应像素的场图值最终解。
上述步骤c中将上述局部增长处理后的待估子块和原种子块合并得到新的种子块，重复上述步骤，直至对全部待估子块完成局部增长处理，获得相应像素的场图值最终解。而对于除这些已估计的像素(即上述已处理子块中的像素)之外的剩余像素采用另一种高分辨率的局部增长处理方法进行处理。即在本发明的一个实施例中，如图3所示，上述步骤400中，以种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计每个像素的场图值最终解的过程还包括高分辨率的增长步骤403，该步骤为：
将种子块和待估子块中的像素作为初始起点，对原始图像中除初始起点之外的剩余像素进行局部增长处理，估计获得原始图像中剩余像素对应的场图值最终解。这里的对原始图像中剩余像素进行局部增长处理、估计获得原始图像中剩余像素对应的场图值最终解的步骤优选包括：分别计算每个剩余像素的初始场图值和与该剩余像素相邻的初始起点的场图值最终解之间的欧氏距离，选取欧氏距离最小的对应的初始场图值作为原始图像中该剩余像素的场图值最终解。在上述过程中，针对每一个剩余像素，计算其周围的初始起点像素的个数，对具有最多周围起始点的剩余像素进行初始场图值的选取，计算所有周围初始点场图值的均值，以此均值对当前像素的场图值进行判断选取，重复以上过程，直至所有剩余像素被处理。当然本发明不限于只采用欧氏距离的计算来估计获得剩余像素的场图值最终解，还可以采用其他局部增长处理的常用方式。
当原始图像中其中一个或多个像素对应的初始场图值包括场图值的真实解和反转解时，上述基于计算欧氏距离的高分辨率的增长步骤403可以为：分别计算每个剩余像素对应的场图值的第一解和第二解与初始起点的场图值最终解之间的欧氏距离，选取场图值的第一解和第二解之中与初始起点的场图值最终解之间的欧氏距离最小的解，作为原始图像中该剩余像素的场图值最终解。上述步骤403主要是利用高分辨率的局部增长法进行处理，即利用已知场图值的像素去估计未知场图值的像素。若未知像素具有两个数值，分别计算这两个数值与周围已估像素的数值之间的欧氏距离d，则具有较小d值的数值段即认为是真实场图值。
上述步骤500中根据每个像素的场图值最终解获得水的图像和脂肪的图像的的过程可以参考以下方法。
选取至少4个回波时间对应的MR图像作为原始图像，利用估计得到的场图值最终解(即fB,ture)和下述公式(5)计算得到水图和脂肪图：
ρ＝A+diag(ψ(fB))S 公式(5)
其中A+,diag(ψ(fB))的定义如上述公式(3)，ρ＝[ρw,ρf]即为计算得到的水图和脂肪图，其结果为复数。
图3是本发明的最优实施例，其利用纯水、纯脂肪组织的正确场图值和反转场图值之间的规律及子块场图值平滑连续的假设，选取场图值分布曲线为三段式的子块作为种子块，而非一个像素作为种子点，其中间曲线段数值即为该子块中真实场图值的分布，这些该子块即选为种子块。然后同时基于低分辨率的增长和高分辨率的增长来对其余像素进行处理，即低分辨率的增长以上述种子块为起点，对场图分布曲线为两段式的待估子块进行局部增长处理，分别计算两段数值与种子块数值之间的欧氏距离，选取欧氏距离小的数值作为待估子块的真实场图值；而高分辨率的增长以上述得到真实场图值的像素为起点，对未估计像素进行增长处理，选取两个场图值之中与已估像素数值欧氏距离最小的值作为fB,ture。
图1或图3为本发明一个实施例方法的流程示意图。应该理解的是，虽然图1或图3的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，其可以以其他的顺序执行。而且，图1或图3中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，其执行顺序也不必然是依次进行，而是可以与其他步骤或者其他步骤的子步骤或者阶段的至少一部分轮流或者交替或者组合地执行。
基于上述磁共振成像的图像水脂分离方法，本发明还提供了一种磁共振成像的图像水脂分离系统，如图4所示，其包括：
场图获取模块601，用于基于磁共振成像获得的原始图像，估计上述原始图像中每个像素的场图值，获得上述每个像素的初始场图值，生成相应的场图；这里的原始图像优选采用至少四个回波时间对应的回波序列信号进行图像重建；
曲线绘制模块602，用于将上述场图进行分块处理，获得多个子块，生成每一个子块中初始场图值的分布曲线；
种子选择模块603，用于基于上述每一个子块中初始场图值的分布曲线，选择上述初始场图值的分布曲线为三段式结构对应的子块，作为种子块；
局部增长处理模块604，用于取上述种子块内初始场图值的分布曲线中的中间段数值作为上述种子块内相应像素的场图值最终解，以上述种子块内每一个像素的场图值最终解为初始值、利用局部增长方法估计上述每个像素的场图值最终解；及
图像分离模块605，用于根据上述每个像素的场图值最终解，获得水的图像和脂肪的图像。
基于上述系统，在本发明的一个实施例中，在上述场图获取模块601中包括：
第一单元，用于针对上述每个像素判断上述初始场图值的个数，当上述初始场图值为两个时，该像素的初始场图值包括第一解和第二解；当上述初始场图值为一个时，将该像素对应的第一解和第二解均赋值为上述初始场图值；
第二单元，用于汇总上述原始图像中每个像素对应的第一解和第二解，分别生成上述第一解对应的第一场图和上述第二解对应的第二场图。
基于上述系统，在本发明的一个实施例中，上述曲线绘制模块602包括：
第三单元，用于分别将上述第一场图和第二场图进行分块处理，获得上述多个子块；
第四单元，用于针对每一个子块，合并统计该子块中相应像素的上述第一解和第二解，获得上述每个子块内所有的初始场图值；
第五单元，用于将上述每一个子块内所有的初始场图值进行排序，获得场图值排序位置；
第六单元，用于绘制表征初始场图值与场图值排序位置的对应关系曲线，生成上述每一个子块中初始场图值的分布曲线。
基于上述系统，在本发明的一个实施例中，如图5所示，上述局部增长处理模块604包括低分辨率的增长单元614，用于将初始场图值的分布曲线为两段式结构的子块作为待估子块，以种子块内每一个像素的场图值最终解为初始值，对待估子块进行局部增长处理，估计获得待估子块内像素的场图值最终解。其中优选，分别计算待估子块初始场图值的分布曲线中两段数值与种子块对应的初始场图值分布曲线中的中间段数值之间的欧氏距离，选取欧氏距离相对较小的一段数值作为待估子块内相应像素的场图值真实解。
如果上述待估子块为在场图中与种子块相邻的子块，则上述局部增长处理模块604还包括高分辨率的增长单元624，用于将种子块和待估子块中的像素作为初始起点，对原始图像中剩余像素进行局部增长处理，估计获得原始图像中剩余像素对应的场图值最终解。其中优选，分别计算每个剩余像素的初始场图值和与该剩余像素相邻的初始起点的场图值最终解之间的欧氏距离，选取欧氏距离最小的对应的初始场图值作为原始图像中该剩余像素的场图值最终解。在本发明的一个实施例中，当原始图像中其中一个或多个像素对应的初始场图值包括场图值的真实解和反转解，则高分辨率的增长单元624中，分别计算每一个剩余像素对应的场图值的第一解和第二解与上述初始起点的场图值最终解之间的欧氏距离，选取场图值的第一解和第二解之中与初始起点的场图值最终解之间的欧氏距离最小的解，作为原始图像中该剩余像素的场图值最终解。
上述系统基于方法建立，其各个功能模块具体实现方式可参见上述方法各个步骤的相关详细说明，在此不累述。以上各个实施例在具体说明中仅只针对相应步骤的实现方式进行了阐述，然后在逻辑不相矛盾的情况下，上述各个实施例是可以相互组合的而形成新的技术方案的，而该新的技术方案依然在本具体实施方式的公开范围内。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个非易失性计算机可读存储介质(如ROM、磁碟、光盘)中，包括若干指令用以使得一台终端设备(可以是手机，计算机，服务器，或者网络设备等)执行本发明各个实施例所述的系统结构和方法。
上述方法和系统中对估计场图进行分块处理，利用纯水和纯脂肪组织的估计场图值分布规律选取种子块，并以这些种子块为起点依次按照低分辨率和高分辨率的局部增长策略，完成对其它种子块和像素的场图值选取，得到准确的待估场图值。本发明相对于现有技术，可选取足够数目的种子点，利用多分辨率的局部增长策略完成对成像区域中各组织场图值的估计和选取，避免对孤立组织中场图值的误判，得到正确的场图值和精确的水脂分离结果；同时本算法最终得到的是场图值，而不是所谓的phasor，该值真实、客观地描述了静磁场B0的不均匀性。为了验证本发明的可行性，利用在体组织实验分别进行了测试。在体组织实验在3T MR系统中进行扫描，在体实验对象为颈部。用GRE多回波序列扫描对象，扫描的序列参数为：TR＝13ms,TE＝5.62，6.88，7.66，8.15ms,翻转角flip angle＝15°,视野FOV＝28.2cm,矩阵大小matrix size＝320×320,层厚slice thickness＝3mm,采集层数slice＝5,带宽bandwidth＝391Hz/pixel。实验数据的处理环境是具有Intel E5-2650v2CPU和64GB RAM的工作站，采用的数据处理软件是MATLAB。附图6在体颈部实验，其中(A)TE＝5.62ms时的原图；(B)选取的种子块；(C)低分率局部增长结果；(D)高分辨率局部增长结果；(E)分离得到的水图；(F)分离得到的脂肪图。附图6中为其中一层的分离结果，从其中可以看出在体实验中，得到的水图和脂肪图中无明显分离错误的组织。
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。