CN105844626B - 一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法 - Google Patents

Abstract

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括：(1)对输入磁共振数据进行非局部均值滤波，计算滤波前后复向量；(2)计算滤波前后缠绕相位和滤波前后缠绕相位差；(3)计算每个像素的向量相位局部偏差和缠绕相位局部偏差；(4)计算缠绕附近像素个数和极点个数；(5)将掩模内像素按其缠绕相位局部偏差值降序排列，计算像素分类阈值；(6)对掩膜内像素分类：平滑无缠绕不连通块和残余像素；(7)解缠绕并合并所有不连通块，根据已解缠绕区域相位对残余像素解缠绕；(8)对解缠绕结果从新解缠绕，依次重复(3)——(7)；(9)把与差添加到解缠绕结果上，重复(8)迭代解缠绕两次。

Description

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法

技术领域

本发明属于磁共振成像技术领域，具体涉及一种在图像域内基于缠绕识别和局部多项式曲面拟合的二维磁共振相位解缠绕方法。

背景技术

磁共振成像(Magnetic resonance imaging,MRI)由于其无电离辐射、能够获得丰富的组织对比度信息以及具有非入侵式检测等优点，已广泛应用于临床医学影像检查。在磁共振成像中，在每个像素上获得的信号数据都是复数，同时包含幅值和相位两部分信息。在大多数临床应用中，仅幅值图像被用于疾病的诊断与研究，而忽视了相位信息。也存在一些磁共振成像技术，如化学偏移成像、磁敏感度成像、磁共振相位对比血流成像、血流测量、人脑相位成像等，利用采集到相位数据去获得对疾病诊断有价值的信息。这些基于相位的磁共振成像技术通常需要获得准确的相位信息。然而，从复数数据中获得的相位被定义在(-π,π]这一区间内。如果真实相位值不在该区间内，那么它将被缠绕到这个区间内，引起相位混叠，这种结果性模糊通常被称为相位缠绕。理论上，真实相位与缠绕相位的差为2kπ，k是整数。可以表示为如下形式：

上式中，φ(x,y)和分别表示位于空间位置(x,y)处的真实相位和缠绕相位。从缠绕相位中准确可靠地恢复出真实像位，对于和相位相关的磁共振成像是非常重要的。

相位解缠绕通常假设相邻像素的真实相位差小于π，如果该条件在所有像素上都得到满足，获得准确的解缠绕结果是比较容易的。然而，当待处理的相位图像中含有强烈的噪声、相位快速变化、或不连通区域时，相邻像素的真实相位差有可能大于π。这种情况下，从缠绕相位图像中准确地恢复出真实相位是非常困难的，尤其对二维及以上的数据来说。在过去很多年里，提出了许多相位解缠绕方法。这些方法主要分成两类：(1)路径跟踪法，(2)最小范数法。

路径追踪法是利用相邻像素的相位梯度进行路径积分来实现相位解缠绕的。如果存在极点，则解缠绕结果依赖于积分路径。大多数路径追踪法尝试通过优化积分路径来处理这种由于极点导致的不连续。例如，Goldstein’s枝切算法先识别出缠绕相位图像中的极点并使用枝切连接它们，然后沿着不穿过枝切线的任何路径进行积分以获得准确的解缠绕相位。基于区域增长的相位解缠绕方法是另外一类著名的路径追踪方法。区域增长解缠绕方法选取相位相对均匀区域的某个像素作为起始的种子点，然后根据已解缠绕相邻像素的相位信息来预测将要被解缠绕像素的相位值。该方法相位解缠绕的结果依赖于区域增长的路径，因此需要鲁棒的区域增长准则来确保相位解缠绕沿着一个可靠的路径。Xu和Cumming提出通过检测相位预测的连贯性来引导解缠绕路径沿着最可靠方向。质量图被提出用来指导区域增长路径：按照从质量由好到坏的顺序依次进行逐像素的相位解缠绕。局部一阶相位偏差和二阶相位偏差曾作为质量来指导区域增长的路径和种子点选取。最近，Witoszynskyj等利用幅值信息来指导解缠绕路径和相位局部相关信息来选取多个种子点。然而上述路径追踪方法中，整个路径积分或区域增长中的任何错误都可能累计至后续像素的解缠绕过程，导致最终解缠绕结果出现严重的缠绕残留。

最小范数法通过最小化未知的真实相位局部偏差与已知的缠绕相位局部偏差的差值来实现相位解缠绕。最简单易用的最小范数相位解缠绕方法可能是最小二乘法，最小二乘法最小化缠绕图像梯度与估计图像梯度的差平方和。更常用的最小范数法是假设真实相位是平滑的，并可以用经验的数学模型表示，如多项式模型、truncated泰勒级、Markov模型等，从而把把相位解缠绕问题转化为模型的参数估计问题。而以模型为基础最小范数方法的优点是相位解缠绕结果一般不受局部区域噪声或相位不连贯的影响，甚至是对某些被较大无信号(signal voids)区域分割开的区域，该方法仍能获得准确的解缠绕结果。然而该方法主要局限是经验的模型很难准确的描述快速或复杂的相位变化。

把缠绕相位图像分成若干个块，先进行块内解缠绕，再利用以块为基本单元的区域增长策略并结合最小二乘方法进行块间解缠绕，直至合并所有的块，即“split andmerge”，被表明是一种有效准确鲁棒的相位解缠绕方案。其中，Jenkinson提出的PRELUDE方法首先把(-π,π]这一区间划分成偶数个等间隔子区间，然后利用这些子区间把相位图像分成若干块，而不需块内解缠绕，因此相位解缠绕只需要进行块与块之间的解缠绕并合并。当信噪比较低时，相邻像素的原始相位差可能会大于2π，PRELUDE方法获得的块将存在块内缠绕，导致该方法不能成功地估计出真实相位。此外，该方法利用块间相邻的像素进行块合并，如果这些相邻的像素位于真实相位存在快速变化的区域，如组织间隙等，基于这些低质量的像素进行块合并可能会降低相位解缠绕的准确性与可靠性。

因此，针对现有技术不足，提出方法利用修改的复数非局部均值滤波器特点，根据缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。本技术通过仿真实验与真实磁共振水脂分离实验对所提方法的性能进行了可靠性验证。

发明内容

本发明的目的在于针对现有技术不足，提供一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，该方法能够较好的在图像域中对缠绕图像实现解缠绕，获得平滑解缠绕结果，应用在三点Dixon技术中能获得较好的水脂分离结果。

本发明的上述目的通过如下技术方案实现：

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：依次包括如下步骤：

步骤(1)根据采集的复数磁共振数据计算其复向量并对采集的复数磁共振数据进行非局部均值滤波以获得滤波后的复向量

步骤(2)根据滤波后的复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤波后缠绕相位的差；

步骤(3)根据滤波后的缠绕相位图，计算每个像素的向量相位局部偏差和缠绕相位局部偏差；

步骤(4)根据向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据极点定义计算感兴趣区A内掩模内(用户兴趣的区域，通过人为设定)极点个数；

步骤(5)根据每个像素的缠绕相位局部偏差值，将掩模感兴趣区内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模感兴趣区内像素分类阈值；

步骤(6)根据计算的像素分类阈值对滤波后的掩模感兴趣区内的像素分类：一类为相位平滑无缠绕的不连通块和，一类为残余像素；

步骤(7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕区域的相位值，按照结合质量引导的局域增长方案法和多项式曲面拟合方法对残余像素解缠绕；

步骤(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)；

步骤(9)把步骤(2)中得到的滤波前缠绕相位与滤波后缠绕相位的差滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)迭代解缠绕两次。解缠绕两次。

上述步骤(1)中利用非局部均值滤波器以相同参数同时对输入的磁共振数据实部和虚部滤波，获得滤波后的复向量上式中表示滤波后缠绕相位，其滤波参数分别是：t＝5；f＝3；k＝1.2,噪声标准差为0.3，i×i＝-1。

上述步骤(2)中计算出滤波后的缠绕相位，并计算滤波前缠绕相位与滤波后缠绕相位的相位差具体是：

根据滤波后的复向量，对其做缠绕操作使得每个像素相位值都在区间(-π,π]内，从而获得滤波后的缠绕相位图像，同理获得滤波前的缠绕相位图像

根据前面计算的滤波前后缠绕相位图像，计算由滤波导致的缠绕相位值的变化；

上述步骤(3)中计算向量相位局部偏差和缠绕相位局部偏差具体是：

根据复向量计算出的向量相位局部偏差(LD-PP)能把掩模内具有真实相位快速变换特点的像素更直观表现出来，即在LD-PP图中，具有较大值得像素代表其所在的区域真实相位变化较快，因此，LD-PP可由下式获得：

上式中，N(x₀,y₀)表示以空间位置(x₀,y₀)为中心的八邻域像素空间坐标集合，为与空间位置(x₀,y₀)相邻的位于(x,y)处像素的复向量，为位于(x₀,y₀)处的复向量。为的复共轭。

由缠绕相位计算出的缠绕相位局部偏差(LD-WP)不仅受到真实相位快速变化的影响，同时也受到缠绕影响，即在LD-WP图中，某像素的值较大可能由真实相位变化较快或有缠绕引起的。LD-WP图可由下式获得：

上式中，为以空间位置(x₀,y₀)为中心的位于(x,y)处的缠绕相位，为位于(x₀,y₀)的缠绕相位。

上述步骤(4)中根据计算向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数和极点个数具体是：

M_W＝|(LD-WP)-(LD-PP)|>0 [3]

上式中，M_w表示缠绕附近像素个数，LD-WP和LD-PP分别是缠绕相位偏差和向量相位局部偏差。

N_P＝|Residues|>0 [4]

上式中，N_p代表极点个数，Residues表示根据极点定义计算出的所有极点。

上述步骤(5)中根据获得的缠绕相位局部偏差值，求得掩模区内像素分类阈值具体是：

对掩模内缠绕像素按照其缠绕相位局部偏差值进行降序排列，选择靠近序号(M_W+A*N_P)处的某个像素缠绕相位局部偏差值作为像素分类阈值Thre,如下式：

Thre_se＝round(M_W+A*N_P) [5]

上式中，A的取值为1，round(Z)表示一个计算距离Z最近整数的函数，Thre_se表示相位分类阈值的位置。

上述步骤(6)根据获得的像素分类阈值Thre对掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素；

Blocks＝(LD-WP)>Thre [6]

上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连通块被重新分类为残余像素，同时，如果不是面积最大的某不连通块距离其他块的边沿欧式距离都大于20个像素，那么这个块内像素也将被归类残余像素。

上述步骤(7)中利用局部多项式曲面拟合方法进行块间相位解缠绕并合并的方法具体为：

步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始块内相位已解缠绕。

步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素具有相同的相位补偿2kπ，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是0.9和0.1。距离最小的块被选为生长块。

步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块，并作为新的已解缠绕区域，具体为：我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近生长块的解缠绕区域中选取n个像素参加曲面拟合,从靠近解缠绕区域的生长块中选取n个像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的整数区间，如[-3,3]。

步骤d，根据步骤c选取的已解缠绕区域的n个像素相位和生长块中的n个像素相位以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优整数补偿。

上述步骤(7)中根据已解缠绕区域相位，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠绕的方法具体为：

步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个，那么从这些生长点中选取LD-PP值最小的像素为生长点。

步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19×19。

步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合并为新的已解缠绕区域，具体为：在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合窗内，假设有P个已解缠绕像素：

上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个P×(M+1)(N+1)大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x,y)的行向量为[1,x,y,xy,...,x^M,y^N,x^My^N]，M和N分别是沿着x和y方向的多项式指数。)。是含有多项式系数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵然后生长点的解缠绕相位可以估计为：

上式中，表示生长点相位估计值，是由生长点的多项式基构成。因此生长点的整数补偿可以计算为：

上述步骤(8)中对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)从新解缠绕。在重复步骤3中，计算向量相位偏差和缠绕相位偏差时，我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

上述步骤(9)中把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这可以表示为：

上式中，表示了添加滤波前后缠绕相位差的缠绕相位，代表滤波后缠绕相位，φ_NLM表示对滤波后缠绕相位解缠绕后的结果。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，第一次我们根据的第二次我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

本发明的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括如下步骤：(1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行非局部均值滤波以获得滤波后的复向量；(2)根据滤波后复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤波后缠绕相位的差；(3)根据滤波后缠绕相位图，计算向量相位局部偏差和缠绕相位局部偏差；(4)根据计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据极点定义计算掩模区内极点个数；(5)根据计算的缠绕相位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模区内像素分类阈值；(6)根据计算的像素分类阈值对滤波后的掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素；(7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕区域的相位值，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠绕；(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)；(9)把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)解缠绕两次。因此，针对现有技术不足，提出方法利用修改的复数非局部均值滤波器特点，缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。本技术通过仿真实验与真实磁共振水脂分离实验对所提方法的性能进行了可靠性验证。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法的流程图；

图2为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中计算得到的幅值图像，滤波前后缠绕相位图像，滤波前后缠绕相位差，LD-PP，LD-WP，LD-WP与LD-PP之差和像素分类结果；

图3为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用具有不同信噪比和不同快速相位变化的仿真数据集1来比较验证提出方法；

图4为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用具有相同信噪比但不同快速相位变化的仿真数据集2来判断提出方法对具有相同信噪比但不同快速相位变换的数据解缠绕结果示意图；

图5为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真实人体膝关节横断面磁共振数据判断提出方法应用在三点Dixon水脂分离技术中的效果示意图；

图6为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真实人体膝关节矢状位磁共振数据判断提出方法应用在三点Dixon水脂分离技术的效果示意图；

图7为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真实人体踝关节矢状位磁共振数据判断提出方法应用在三点Dixon水脂分离技术中的效果示意图；

具体实施方式

下面结合具体的实施例对本发明进行详细描述。

实施例1

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括如下步骤：

(1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行非局部均值滤波以获得滤波后的复向量；

滤波前后复向量采用如下方式计算获得：

滤波前复向量可由获取的复数磁共振数据除以其幅值获得。滤波后复向量由滤波后的复数磁共振数据除以相应幅值获得。其中非局部均值滤波参数分别是：t＝3；f＝1；k＝1.0；噪声标准差为0.3,以相同滤波参数同时对磁共振数据的实部和虚部进行滤波。

(2)根据滤波后复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤波后缠绕相位的差；

计算出滤波后的缠绕相位滤波前缠绕相位与滤波后缠绕相位的相位差具体是：

根据计算出的滤波前后复向量和分别对其做缠绕操作使得每个像素相位值都在区间(-π,π]内，从而获得滤波前后的缠绕相位图像和

根据前面计算的滤波前后缠绕相位图像和计算由滤波导致的缠绕相位变化

(3)根据滤波后缠绕相位图，计算向量相位的局部偏差和缠绕相位的局部偏差；

计算向量相位的局部偏差(LD-PP)和缠绕相位的局部偏差(LD-WP)具体是：

根据复向量计算出的向量相位局部偏差LD-PP能把掩模内具有真实相位快速变换特点的像素更直观表现出来，即在LD-PP图中，具有较大值的像素代表其所在的区域真实相位变化较快，LD-PP可由下式获得：

上式中，N(x₀,y₀)表示以空间位置(x₀,y₀)为中心的八邻域像素空间坐标集合，为与空间位置(x₀,y₀)相邻的位于(x,y)处像素的复向量，为位于(x₀,y₀)处的复向量。为的复共轭。

由缠绕相位计算出的缠绕相位局部偏差LD-WP不仅受到真实相位快速变化的影响，同时也受到缠绕影响，即在LD-WP图中，如果某像素的LD-WP值较大，这可能是由真实相位变化较快或缠绕引起的。LD-WP图可由下式获得：

上式中，为与空间位置(x₀,y₀)相邻的位于(x,y)处的缠绕相位，为位于(x₀,y₀)的缠绕相位。

(4)根据计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据极点定义计算掩模区内极点个数；

根据计算向量相位局部偏差LD-PP和缠绕相位局部偏差LD-WP估计缠绕附近像素个数和极点个数具体是：

M_W＝|(LD-WP)-(LD-PP)|>0 [3]

上式中，M_w表示缠绕附近像素个数。

N_P＝|Residues|>0 [4]

上式中，N_p代表极点个数，Residues表示根据极点定义计算出的所有极点。

(5)根据计算的缠绕相位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模区内像素分类阈值；

根据缠绕相位局部偏差值的大小，求得掩模区内像素分类阈值具体是：

对掩模内缠绕像素按照其LD-WP值的大小进行降序排列，选择靠近(M_W+A*N_P)处的某个像素的LD-WP值作为像素分类阈值Thre,如下式：

Thre_Se＝round(M_W+A*N_P) [5]

上式中，A的取值为1，round(Z)表示一个计算距离Z最近整数的函数，Thre_se表示相位分类阈值的位置。

(6)根据计算的像素分类阈值对滤波后的掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素；

根据像素分类阈值Thre对掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素；

Blocks＝(LD-WP)<Thre [6]

上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连通块被重新分类为残余像素。如果某个不是面积最大的不连通块距离其他块的最小边沿欧式距离都大于20个像素，那么这个块被重新归类为残余像素。

(7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕区域的相位值，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠绕；

其中利用局部多项式曲面拟合方法进行块间相位解缠绕的方法具体为：

步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始块内相位已解缠绕。

步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素具有相同的相位补偿2kπ，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是0.9和0.1。距离最小的块被选为生长块。

步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块为新的已解缠绕区域，具体为：我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近生长块的已解缠绕区域中选取n个像素参加曲面拟合,从靠近已解缠绕区域的生长块中选取n个像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的整数区间，如[-3,3]。

步骤d，根据步骤c选取的已解缠绕区域的n个像素相位和生长块中的n个像素相位以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优整数补偿。

根据合并后的块，按照结合质量引导的区域增长方案和多项式曲面拟合方法对残余像素解缠绕的方法具体为：

步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个，那么从这些生长点中选取LD-PP值最小的像素为生长点。

步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19×19。

步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合并为新的已解缠绕区域，具体为：在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合窗内，假设有P个已解缠绕像素：

上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个P×(M+1)(N+1)大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x,y)的行向量为[1,x,y,xy,...,x^M,y^N,x^My^N]，M和N分别是沿着x和y方向的多项式指数。)。是含有多项式系数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵然后生长点的解缠绕相位可以估计为：

上式中，表示生长点相位估计值，是由生长点的多项式基构成。因此生长点的整数补偿可以计算为：

(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)；

对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)从新解缠绕。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

(9)把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)解缠绕两次。

把步骤(2)中得到的滤波前与滤波后缠绕相位差添加到步骤(8)得到的解缠绕相位图像上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这可以表示为：

上式中，表示添加滤波前与滤波后缠绕相位差的缠绕相位，代表非局部均值滤波后的缠绕相位，φ_NLM表示对滤波后的缠绕相位解缠绕后的结果。在重复步骤3中，计算向量相位偏差和缠绕相位偏差时，第一次我们根据的第二次我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

与现有技术相比，提出方法利用修改的复数非局部均值滤波器特点，根据缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。

本技术通过仿真实验与真实磁共振水脂分离实验对所提方法的性能进行了实验验证，实验结果表明提出的方法能很好实现对缠绕相位图像解缠绕，水脂分离结果中无明显水脂互换现象存在。

实施例2

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括以下步骤：

(1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行非局部均值滤波以获得滤波后的复向量。在图1中，第一行分别显示了一个膝关节横断面磁共振数据的幅值(Magnitude)图像,缠绕相位(Wrapped phase)图像，非局部均值滤波后获得缠绕相位(Filtered wrapped phase)和滤波前后缠绕相位差。第二行显示了由滤波前磁共振数据获得的缠绕相位局部偏差(LD-WP)、向量相位局部偏差(LD-PP)、缠绕附近像素(Wraps locations)和像素分类后的不连通块(Blocks)和残余像素(Residual pixels)；第三行显示了由滤波后磁共振数据获得的缠绕相位局部偏差、向量相位局部偏差、缠绕附近像素和像素分类后的不连通块和残余像素，

在对掩模区像素分类中，由滤波前磁共振数据获得的缠绕附近像素和极点个数分别是3374和268，在缠绕相位局部偏差获得的分块阈值为1.997，像素分类结果如图1所示。由滤波后磁共振数据获得的缠绕附近像素和极点个数分别是2185和0，在缠绕相位局部偏差图获得的分块阈值为0.5679，分块结果如图1所示。

因此，利用非局部均值滤波能有效降低极点对像素分类的影响，获得较好的像素分类结果。

(2)根据滤波后复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤波后缠绕相位的差。利用步骤(1)中计算的滤波前后复向量，再利用Matlab函数库中angle()函数，可以分别获得滤波前后缠绕相位及滤波前后缠绕相位差(请参阅图1)。

(3)根据滤波缠绕相位图，计算向量相位局部偏差和缠绕相位局部偏差。

分别利用步骤(1)和(2)中计算的滤波后复向量和缠绕相位，可以分别获得向量相位局部偏差(LD-PP)和缠绕相位局部偏差(LD-WP)(请参阅图1)。其中计算向量相位局部偏差和缠绕相位局部偏差具体是：

根据复向量计算出的向量相位局部偏差LD-PP能把掩模内具有真实相位快速变换特点的像素更直观表现出来，即在LD-PP图中，具有较大值的像素代表其所在的区域真实相位变化较快，LD-PP可由下式获得：

上式中，N(x₀,y₀)表示以空间位置(x₀,y₀)为中心的八邻域像素空间坐标集合，为与空间位置(x₀,y₀)相邻的位于(x,y)处像素的复向量，为位于(x₀,y₀)处的复向量。为的复共轭。如图1所示，LD-PP的较大值发生在真实相位快速变换区域。

由缠绕相位计算出的缠绕相位局部偏差LD-WP不仅受到真实相位快速变化的影响，同时也受到缠绕影响，即在LD-WP图中，如果某像素的LD-WP值较大，这可能是由真实相位变化较快或缠绕引起的。LD-WP图可由下式获得：

上式中，为与空间位置(x₀,y₀)相邻的位于(x,y)处的缠绕相位，为位于(x₀,y₀)的缠绕相位。如图1所示，LD-WP的较大值出现在缠绕发生区域和/或相位快速变换区域。

因此，在本实施方案中，根据缠绕在计算LD-WP和LD-PP中的不同表现，可以有效的估计缠绕附近像素的位置。

(4)根据计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据极点定义计算掩模区内极点个数。

根据步骤(3)计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数和极点个数具体是：

M_W＝|(LD-WP)-(LD-PP)|>0 [3]

上式中，M_w表示缠绕附近像素个数。在本例中M_w为2185(请参阅图1)。

N_P＝|Residues|>0 [4]

上式中，N_p代表极点个数，Residues表示根据极点定义计算出的所有极点。在本例中N_p为0(请参阅图1)。

(5)根据计算的缠绕相位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模区内像素分类阈值。

根据步骤(3)计算的缠绕相位局部偏差值的大小，求得掩模区内像素分类阈值具体是：

对掩模内缠绕像素按照其缠绕相位局部偏差值的大小进行降序排列，选择最接近序号(Mw+A*N_P)处的某个像素的LD-WP值作为像素分类阈值Thre,如下式：

Thre_Se＝M_w+N_P [5]

上式中，A的取值为1，round(Z)表示一个计算距离Z最近整数的函数，Thre_Se表示相位分类阈值的位置。在本例中Thre_Se和Thre分别等于2185和0.5679(请参阅图1)。

(6)根据计算的像素分类阈值对滤波后的掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素。

根据步骤(5)计算的像素分类阈值Thre对掩模区内像素分类，具体如下式：

Blocks＝(LD-WP)<Thre [6]

上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连通块被重新分类为残余像素。如果某个不是面积最大的不连通块距离其他块的最小边沿欧式距离都大于20个像素，那么这个块被重新归类为残余像素。

(7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕区域的相位值，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠绕。

其中利用局部多项式曲面拟合方法进行块间相位解缠绕的方法具体为：

步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始块内相位已解缠绕。

步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素具有相同的相位补偿2kπ，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是0.9和0.1。距离最小的块被选为生长块。

步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块为新的已解缠绕区域，具体为：我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近生长块的已接缠绕区域中选取n个像素参加曲面拟合,从靠近已接缠绕区域的生长块中选取n个像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的整数区间，如[-3,3]。

步骤d，根据步骤c选取的已解缠绕区域的n个像素相位和生长块中的n个像素相位以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优整数补偿。

根据合并后的块，按照结合质量引导的区域增长方案和多项式曲面拟合方法对残余像素解缠绕的方法具体为：

步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个，那么从这些生长点中选取LD-PP值最小的像素为生长点。

步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19×19。

步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合并为新的已解缠绕区域，具体为：在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合窗内，假设有P个已解缠绕像素：

上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个P×(M+1)(N+1)大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x,y)的行向量为[1,x,y,xy,...,x^M,y^N,x^My^N]，M和N分别是沿着x和y方向的多项式指数。)。是含有多项式系数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵然后生长点的解缠绕相位可以估计为：

上式中，表示生长点相位估计值，是由生长点的多项式基构成。因此生长点的整数补偿可以计算为：

从执行残余像素解缠绕步骤后的解缠绕结果可以看出，存在缠绕的残余像素明显降低(请参阅图2)。

(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)。

对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、(4)、(5)、(6)和(7)从新解缠绕。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

(9)把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)解缠绕两次。

把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位图像上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这可以表示为：

上式中，表示添加滤波前与滤波后缠绕相位差的缠绕相位，代表非局部均值滤波后的缠绕相位，φ_NLM表示对滤波后的缠绕相位解缠绕后的结果。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，第一次我们根据的第二次我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。

从对滤波前缠绕相位图像的最终解缠绕结果和水脂分离情况可以看出，解缠绕结果很平滑，无论是在快速相位变化还是低信噪比区域都无明显的缠绕像素存在，从而在水脂分离图像中无明显水脂互换现象存在(请参阅图5)。

与现有技术相比，提出方法利用修改的复数非局部均值滤波器特点，缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。

本技术通过仿真实验与真实水脂分离磁共振成像实验对所提方法的性能进行了实验验证，实验结果表明提出方法能很好实现对缠绕相位图像解缠绕，水脂分离结果中无明显水脂互换现象存在。

实施例3。

通过本发明的方法对仿真和真实磁共振相位进行相位解缠绕实验，本实施例选取其中部分实验结果进行分析比较。

图3为具有信噪比从6等间隔变化到0.5的圆型仿真数据的相位解缠绕实验结果。未添加噪声的仿真数据真实相位为高斯相位曲面，相位变化区间分别为0到25和50弧度(rad)，即高度(Height)等于25和50弧度(rad)。添加均值为0，标准差为25随机噪声到的生成的复数仿真数据上以生成具有不同信噪比的复数仿真数据。其中第一行展示了幅值(magnitude)图像；第二行分别展示了高度为25rad的相关数据，分别是缠绕相位(Wrappedphase)图像，PRELUDE方法解缠绕结果，PRELUDE解缠绕结果与真实相位和添加噪声的差(Difference),提出方法的(Proposed method)解缠绕结果，提出方法的解缠绕结果与真实相位和添加噪声的差(Difference)；第三行分别展示了高度为50rad的相关数据，分别是缠绕相位图像，PRELUDE方法解缠绕结果，PRELUDE解缠绕结果与真实相位和添加噪声的差,提出方法解缠绕结果，提出方法解缠绕结果与真实相位和添加噪声的差。

比较两种方法的解缠绕结果，可以看出：采用现有经典算法PRELUDE的解缠绕结果错误像素较多，且有成块的错误出现，而采用本发明算法的解缠绕结果大大改善了这一情况。

图3为两种算法对具有相同信噪比但不同快速相位变化的仿真数据解缠绕结果。未添加噪声前真实相位变化区间分别为0到50、100、150和200rad，添加均值为0，标准差为10的随机噪声到复数仿真数据上。其中，第一列最上面图像为幅值图像，其余的为高度等于50、100、150和200rad的缠绕相位图像；第二列为PRELUDE方法解缠绕结果；第三列为PRELUDE解缠绕结果与真实相位和添加噪声的差；第四列为提出方法解缠绕结果；第五列为提出方法的解缠绕结果与真实相位和添加噪声的差。

可见，PRELUDE解缠绕结果(图4第二和第三列)随着缠绕复杂度的增加，错误像素逐渐增加，而当高度大于等于150rad后，在中心区域出现了明显的缠绕块。而在本发明算法的解缠绕结果中，虽然也有错误像素的出现，但是当高度等于150rad后，无明显的错误块出现(图4第四和第五列，与PRELUDE解缠绕结果相比，更低的相对误差也直观上说明了本发明的方法解缠绕结果更加准确。

图5、6和7为两种算法分别对人体不同部位(膝关节横断面、膝关节矢状面和踝关节矢状面)的相位解缠绕结果及应用到三点Dixon技术中进行水脂分离的结果。在图4、5和6中，第一行分别为缠绕相位(Wrapped phase)，反相位幅值(Out-phase magnitude)和同相位幅值(In-phase magnitude)；第二行分别是PRELUDE相位解缠绕结果，和使用PRELUDE相位解缠绕结果的水脂分离结果，包括水图(Water)，脂肪图(Fat)及脂肪分数图(FF)；第三行分别是提出方法相位解缠绕结果，和使用提出方法相位解缠绕结果的水脂分离结果，包括水图，脂肪图及脂肪分数图。

可见，PRELUDE在三种数据类型的解缠绕结果中(图6、7中第二行)都有缠绕残留，而残留的缠绕会导致水脂分离结果中出现了明显的互换现象(图6、7中第二行白色箭头处)。而在本发明算法的解缠绕结果中，而虽然偶尔也有缠绕残留(图7中第二行白色箭头处)，但与PRELUDE结果相比，互换面积明显减少。对掩模区域的脂肪分数计算无明显影响。因此，真实磁共振数据也说明了本发明的方法解缠绕结果更加准确。

通过大量仿真和真实数据实验证明，本发明的方法能够能利用缠绕在向量相位偏差和缠绕相位偏差的不同表现和修改的非局部均值滤波器实现掩模区内像素的有效分类，而利用局部多项式曲面拟合方法能有效的合并平滑的不连通块，最后再利用质量引导的区域增长和局部多项式曲面拟合方法能有效的实现对残余像素的真实相位估计，从而获得平滑解缠绕结果。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (9)

1.一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：包括如下步骤：

步骤(1)根据采集的复数磁共振数据计算其复向量并对采集的复数磁共振数据进行非局部均值滤波以获得滤波后的复向量

步骤(2)根据滤波后的复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤波后缠绕相位的差；

步骤(3)根据滤波后的缠绕相位图，计算每个像素的向量相位局部偏差和缠绕相位局部偏差；

步骤(4)根据向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据极点定义计算掩模内极点个数；

步骤(5)根据每个像素的缠绕相位局部偏差值，将掩模内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模内像素分类阈值；

步骤(6)根据像素分类阈值对滤波后的掩模内的像素分类：一类为相位平滑无缠绕的不连通块，一类为残余像素；

步骤(7)根据局部多项式曲面拟合方法解缠绕并合并所有不连通块，然后根据已解缠绕区域的相位，按照结合质量引导的区域增长方法和多项式曲面拟合方法对残余像素解缠绕；根据已解缠绕区域相位，按照结合质量引导的区域增长方法和多项式曲面拟合方法对残余像素解缠绕的方法具体为：

步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个，从这些生长点中选取LD-PP值最小的像素为生长点；

步骤b，选择残余像素解缠绕曲面拟合窗的大小；

步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点合并为新的已解缠绕区域，具体为：在拟合窗内，已解缠绕像素的真实相位值拟合成一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计，在拟合窗内，有P个已解缠绕像素：

上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个P×(M+1)(N+1)大小的矩阵，每一行是由每一个已解缠绕像素的多项式基，如在空间位置(x,y)的行向量为[1,x,y,xy,...,x^M,y^N,x^My^N]，M和N分别是沿着x和y方向的多项式指数；是含有多项式系数的列向量，E是含有残余误差的列向量；使用最小二乘方法求系数矩阵然后生长点的解缠绕相位估计为：

上式中，表示生长点相位估计值，是由生长点的多项式基构成，因此生长点的整数补偿计算为：

步骤(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，用前一循环计算出解缠绕相位替换步骤(3)中根据滤波后的缠绕相位，然后依次重复步骤(3)——(7)；

步骤(9)把步骤(2)中得到的滤波前缠绕相位与滤波后缠绕相位的差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)迭代解缠绕两次。

2.根据权利要求1所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：步骤(1)中利用非局部均值滤波器以相同参数同时对输入的磁共振数据实部和虚部滤波，获得滤波后的复向量上式中表示滤波后缠绕相位，其滤波参数分别是：t＝5；f＝3；k＝1.2,噪声标准差为0.3，i×i＝-1。

3.根据权利要求2所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(2)中计算出滤波后的缠绕相位，并计算滤波前缠绕相位与滤波后缠绕相位的相位差具体是：

根据滤波后的复向量对其做缠绕操作使得每个像素相位值都在区间(-π,π]内，从而获得滤波后的缠绕相位图像，同理获得滤波前的缠绕相位图像

根据前面计算的滤波前后缠绕相位图像，计算由滤波导致的缠绕相位值的变化即滤波前缠绕相位与滤波后缠绕相位的相位差。

4.根据权利要求3所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(3)中计算向量相位局部偏差和缠绕相位局部偏差具体是：

根据复向量计算出的向量相位局部偏差能把掩模内具有真实相位快速变换特点的像素更直观表现出来，具有较大向量相位局部偏差值的像素代表其所在的区域真实相位变化较快，因此，向量相位局部偏差可由下式获得：

上式中，N(x₀,y₀)表示以空间位置(x₀,y₀)为中心的八邻域像素空间坐标集合，为与空间位置(x₀,y₀)相邻的位于(x,y)处像素的复向量，为位于(x₀,y₀)处的复向量，为的复共轭；

由缠绕相位计算出的缠绕相位局部偏差不仅受到真实相位快速变化的影响，同时也受到缠绕影响，某像素的缠绕相位局部偏差值较大是由真实相位变化较快或缠绕引起的，缠绕相位局部偏差可由下式获得：

上式中，LD-WP表示缠绕相位局部偏差，为以空间位置(x₀,y₀)为中心的位于(x,y)处的缠绕相位，为位于(x₀,y₀)的缠绕相位。

5.根据权利要求4所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(4)中根据计算向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数和极点个数具体是：

M_w＝|(LD-WP)-(LD-PP)|>0 [3]

上式中，M_w表示缠绕附近像素个数，LD-WP和LD-PP分别是缠绕相位偏差和向量相位局部偏差；

N_P＝|Residues|>0 [4]

上式中，N_p代表极点个数，Residues表示根据极点定义计算出的所有极点。

6.根据权利要求5所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(5)中根据获得的缠绕相位局部偏差值，求得掩模区内像素分类阈值具体是：

根据每个像素的缠绕相位局部偏差值，将掩模内像素降序排列，然后根据缠绕附近像素个数和极点个数获得掩模内像素分类阈值；

对掩模内缠绕像素按照其缠绕相位局部偏差值进行降序排列，选择靠近序号(M_w+A*N_P)处的某个像素缠绕相位局部偏差值作为像素分类阈值Thre,如下式：

Thre_se＝round(M_w+A*N_P) [5]

上式中，A的取值为1，round(Z)表示一个计算距离Z最近整数的函数，Thre_se表示相位分类阈值的位置。

7.根据权利要求6所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(6)根据获得的像素分类阈值Thre对掩模区内像素分类：相位平滑无缠绕的不连通块和残余像素；

Blocks＝(LD-WP)>Thre [6]

上式中，Blocks表示无缠绕的不连通块，面积小于10个像素的不连通块被重新分类为残余像素，同时，如果不是面积最大的某不连通块距离其他块的边沿欧式距离都大于20个像素，那么这个块内像素也将被归类残余像素。

8.根据权利要求7所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(7)中利用局部多项式曲面拟合方法进行块间相位解缠绕并合并的方法具体为：

步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始块内相位已解缠绕；

步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素具有相同的相位补偿2kπ，k是整数，在确定生长块之前，计算已解缠绕区域与所有未解缠绕块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是0.9和0.1，距离最小的块被选为生长块；

步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块，并作为新的已解缠绕区域，具体为：已解缠绕区域内相位拟合成一个平滑的多项式曲面，从靠近生长块的已解缠绕区域中选取n个像素参加曲面拟合,从靠近已解缠绕区域的生长块中选取n个像素参加曲面拟合，同时假设生长块的整数补偿k的选取范围为一个较小的整数区间；

步骤d，根据步骤c选取的已解缠绕区域中的n个像素相位和生长块中的n个像素相位以及位于设定区间内的整数k进行曲面拟合，使多项式曲面拟合误差最小的k被选为最优整数补偿。

9.根据权利要求8所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：所述步骤(9)中把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，表示为：

上式中，表示了添加滤波前后缠绕相位差的缠绕相位，代表滤波后缠绕相位，φ_NLM表示对滤波后缠绕相位解缠绕后的结果；在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，第一次根据的第二次是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。