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

一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法
【专利摘要】一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括：(1)对输入磁共振数据进行非局部均值滤波，计算滤波前后复向量；(2)计算滤波前后缠绕相位和滤波前后缠绕相位差；(3)计算每个像素的向量相位局部偏差和缠绕相位局部偏差；(4)计算缠绕附近像素个数和极点个数；(5)将掩模内像素按其缠绕相位局部偏差值降序排列，计算像素分类阈值；(6)对掩膜内像素分类：平滑无缠绕不连通块和残余像素；(7)解缠绕并合并所有不连通块，根据已解缠绕区域相位对残余像素解缠绕；(8)对解缠绕结果从新解缠绕，依次重复(3)——(7)；(9)把与差添加到解缠绕结果上，重复(8)迭代解缠绕两次。
【专利说明】
一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法
技术领域
[0001] 本发明属于磁共振成像技术领域，具体涉及一种在图像域内基于缠绕识别和局部 多项式曲面拟合的二维磁共振相位解缠绕方法。
【背景技术】
[0002] 磁共振成像(Magnetic resonance imaging，MRI)由于其无电离福射、能够获得丰 富的组织对比度信息以及具有非入侵式检测等优点，已广泛应用于临床医学影像检查。在 磁共振成像中，在每个像素上获得的信号数据都是复数，同时包含幅值和相位两部分信息。 在大多数临床应用中，仅幅值图像被用于疾病的诊断与研究，而忽视了相位信息。也存在一 些磁共振成像技术，如化学偏移成像、磁敏感度成像、磁共振相位对比血流成像、血流测量、 人脑相位成像等，利用采集到相位数据去获得对疾病诊断有价值的信息。这些基于相位的 磁共振成像技术通常需要获得准确的相位信息。然而，从复数数据中获得的相位被定义在 (-31,31]这一区间内。如果真实相位值不在该区间内，那么它将被缠绕到这个区间内，引起相 位混叠，这种结果性模糊通常被称为相位缠绕。理论上，真实相位与缠绕相位的差为2k3I，k 是整数。可以表示为如下形式：
[0003]
[11]
[0004] 上式中，Φ (x，y)和炉(U)分别表不位于空间位置（x，y)处的真实相位和缠绕相 位。从缠绕相位中准确可靠地恢复出真实像位，对于和相位相关的磁共振成像是非常重要 的。
[0005] 相位解缠绕通常假设相邻像素的真实相位差小于π，如果该条件在所有像素上都 得到满足，获得准确的解缠绕结果是比较容易的。然而，当待处理的相位图像中含有强烈的 噪声、相位快速变化、或不连通区域时，相邻像素的真实相位差有可能大于I这种情况下， 从缠绕相位图像中准确地恢复出真实相位是非常困难的，尤其对二维及以上的数据来说。 在过去很多年里，提出了许多相位解缠绕方法。这些方法主要分成两类：（1)路径跟踪法， (2)最小范数法。
[0006] 路径追踪法是利用相邻像素的相位梯度进行路径积分来实现相位解缠绕的。如果 存在极点，则解缠绕结果依赖于积分路径。大多数路径追踪法尝试通过优化积分路径来处 理这种由于极点导致的不连续。例如，Goldstein's枝切算法先识别出缠绕相位图像中的极 点并使用枝切连接它们，然后沿着不穿过枝切线的任何路径进行积分以获得准确的解缠绕 相位。基于区域增长的相位解缠绕方法是另外一类著名的路径追踪方法。区域增长解缠绕 方法选取相位相对均匀区域的某个像素作为起始的种子点，然后根据已解缠绕相邻像素的 相位信息来预测将要被解缠绕像素的相位值。该方法相位解缠绕的结果依赖于区域增长 的路径，因此需要鲁棒的区域增长准则来确保相位解缠绕沿着一个可靠的路径。Xu和 Cumming提出通过检测相位预测的连贯性来引导解缠绕路径沿着最可靠方向。质量图被提 出用来指导区域增长路径:按照从质量由好到坏的顺序依次进行逐像素的相位解缠绕。局 部一阶相位偏差和二阶相位偏差曾作为质量来指导区域增长的路径和种子点选取。最近， Witoszynskyj等利用幅值信息来指导解缠绕路径和相位局部相关信息来选取多个种子点。 然而上述路径追踪方法中，整个路径积分或区域增长中的任何错误都可能累计至后续像素 的解缠绕过程，导致最终解缠绕结果出现严重的缠绕残留。
[0007] 最小范数法通过最小化未知的真实相位局部偏差与已知的缠绕相位局部偏差的 差值来实现相位解缠绕。最简单易用的最小范数相位解缠绕方法可能是最小二乘法，最小 二乘法最小化缠绕图像梯度与估计图像梯度的差平方和。更常用的最小范数法是假设真实 相位是平滑的，并可以用经验的数学模型表示，如多项式模型、truncated泰勒级、MarkoV模 型等，从而把把相位解缠绕问题转化为模型的参数估计问题。而以模型为基础最小范数方 法的优点是相位解缠绕结果一般不受局部区域噪声或相位不连贯的影响，甚至是对某些被 较大无信号（signal voids)区域分割开的区域，该方法仍能获得准确的解缠绕结果。然而 该方法主要局限是经验的模型很难准确的描述快速或复杂的相位变化。
[0008] 把缠绕相位图像分成若干个块，先进行块内解缠绕，再利用以块为基本单元的区 域增长策略并结合最小二乘方法进行块间解缠绕，直至合并所有的块，即"split and merge"，被表明是一种有效准确鲁棒的相位解缠绕方案。其中Jenkinson提出的PRELUDE方 法首先把(_m]这一区间划分成偶数个等间隔子区间，然后利用这些子区间把相位图像分 成若干块，而不需块内解缠绕，因此相位解缠绕只需要进行块与块之间的解缠绕并合并。当 信噪比较低时，相邻像素的原始相位差可能会大于2jt，PRELUDE方法获得的块将存在块内缠 绕，导致该方法不能成功地估计出真实相位。此外，该方法利用块间相邻的像素进行块合 并，如果这些相邻的像素位于真实相位存在快速变化的区域，如组织间隙等，基于这些低质 量的像素进行块合并可能会降低相位解缠绕的准确性与可靠性。
[0009] 因此，针对现有技术不足，提出方法利用修改的复数非局部均值滤波器特点，根据 缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证 块内相位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解 缠绕后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块 与逐像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。本技术通过仿真 实验与真实磁共振水脂分离实验对所提方法的性能进行了可靠性验证。

【发明内容】

[0010] 本发明的目的在于针对现有技术不足，提供一种基于缠绕识别和局部曲面拟合磁 共振相位解缠绕方法，该方法能够较好的在图像域中对缠绕图像实现解缠绕，获得平滑解 缠绕结果，应用在三点Dixon技术中能获得较好的水脂分离结果。
[0011] 本发明的上述目的通过如下技术方案实现：
[0012] -种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于：依次包 括如下步骤：
[0013] 步骤（1)根据采集的复数磁共振数据计算其复向量&，并对采集的复数磁共振数 据进行非局部均值滤波以获得滤波后的复向量
[0014] 步骤(2)根据滤波后的复向量《计算出滤波后的缠绕相位图，并计算滤波前缠 绕相位与滤波后缠绕相位的差；
[0015] 步骤(3)根据滤波后的缠绕相位图，计算每个像素的向量相位局部偏差和缠绕相 位局部偏差；
[0016] 步骤(4)根据向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并 根据极点定义计算感兴趣区A内掩模内（用户兴趣的区域，通过人为设定)极点个数；
[0017] 步骤(5)根据每个像素的缠绕相位局部偏差值，将掩模感兴趣区内像素降序排列， 然后根据缠绕附近像素个数和极点个数获得掩模感兴趣区内像素分类阈值；
[0018] 步骤(6)根据计算的像素分类阈值对滤波后的掩模感兴趣区内的像素分类:一类 为相位平滑无缠绕的不连通块和，一类为残余像素；
[0019] 步骤(7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解 缠绕区域的相位值，按照结合质量引导的局域增长方案法和多项式曲面拟合方法对残余像 素解缠绕；
[0020] 步骤(8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤 (3)、（4)、（5)、（6)和(7);
[0021] 步骤(9)把步骤(2)中得到的滤波前缠绕相位与滤波后缠绕相位的差滤波前后缠 绕相位差添加到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠 绕相位图像仍为缠绕相位图像，然后依次重复步骤(8)迭代解缠绕两次。解缠绕两次。
[0022] 上述步骤（1)中利用非局部均值滤波器以相同参数同时对输入的磁共振数据实部 和虚部滤波，获得滤波后的复P
良示滤波后缠绕相位， 其滤波参数分别是:t = 5;f = 3;k=1.2,噪声标准差为0.3，i X i = -1。
[0023]上述步骤(2)中计算出滤波后的缠绕相位，并计算滤波前缠绕相位与滤波后缠绕 相位的相位差具体是：
[0024]根据滤波后的复向量，对其做缠绕操作使得每个像素相位值都在区间（-13I]内， 从而获得滤波后的缠绕相位图像，同理获得滤波前的缠绕相位图像P ;
[0025] 根据前面计算的滤波前后缠绕相位图像，计算由滤波导致的缠绕相位值的变化；
[0026] 上述步骤(3)中计算向量相位局部偏差和缠绕相位局部偏差具体是：
[0027] 根据复向量计算出的向量相位局部偏差(LD-PP)能把掩模内具有真实相位快速变 换特点的像素更直观表现出来，即在LD-PP图中，具有较大值得像素代表其所在的区域真实 相位变化较快，因此，LD-PP可由下式获得：
[0028]
[1]
[0029]上式中，N(XQ，yQ)表示以空间位置（XQ，yQ)为中心的八邻域像素空间坐标集合， 4〇^')为与空间位置(^7〇)相邻的位于&，7)处像素的复向量，以.1^'())为位于(如7())处 的复向量。6 (.W)为巧(X刃的复共辄。
[0030] 由缠绕相位计算出的缠绕相位局部偏差(LD-WP)不仅受到真实相位快速变化的影 响，同时也受到缠绕影响，即在LD-WP图中，某像素的值较大可能由真实相位变化较快或有 缠绕引起的。LD-WP图可由下式获得：
[0031]
[2]
[0032] 上式中，p(x,.V)为以空间位置(xo，yQ)为中心的位于(x，y)处的缠绕相位，(?，.V0) 为位于(XQ，yo)的缠绕相位。
[0033] 上述步骤(4)中根据计算向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像 素个数和极点个数具体是：
[0034] Mw= I (LD-WP)-(LD-PP) | >0 [3]
[0035] 上式中，Mw表示缠绕附近像素个数，LD-WP和LD-PP分别是缠绕相位偏差和向量相 位局部偏差。
[0036] Np= I Residues I >0 [4]
[0037]上式中，Np代表极点个数，Residues表示根据极点定义计算出的所有极点。
[0038] 上述步骤(5)中根据获得的缠绕相位局部偏差值，求得掩模区内像素分类阈值具 体是：
[0039] 对掩模内缠绕像素按照其缠绕相位局部偏差值进行降序排列，选择靠近序号(Mw+ A*NP)处的某个像素缠绕相位局部偏差值作为像素分类阈值Thre，如下式：
[0040] Threse = round (Mw+A*Np) [5]
[0041 ]上式中，A的取值为I，r〇Und(Z)表示一个计算距离Z最近整数的函数，Threse表示相 位分类阈值的位置。
[0042]上述步骤(6)根据获得的像素分类阈值Thre对掩模区内像素分类:相位平滑无缠 绕的不连通块和残余像素；
[0043] Blocks = (LD-ffP) >Thre [6]
[0044] 上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连 通块被重新分类为残余像素，同时，如果不是面积最大的某不连通块距离其他块的边沿欧 式距离都大于20个像素，那么这个块内像素也将被归类残余像素。
[0045] 上述步骤(7)中利用局部多项式曲面拟合方法进行块间相位解缠绕并合并的方法 具体为：
[0046] 步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始 块内相位已解缠绕。
[0047] 步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素 具有相同的相位补偿2k3i，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕 块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是〇. 9和0.1。距离最小 的块被选为生长块。
[0048]步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块，并作为新的已 解缠绕区域，具体为:我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近 生长块的解缠绕区域中选取η个像素参加曲面拟合，从靠近解缠绕区域的生长块中选取η个 像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的 整数区间，如[_3,3]。
[0049] 步骤d，根据步骤c选取的已解缠绕区域的η个像素相位和生长块中的η个像素相位 以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优 整数补偿。
[0050] 上述步骤(7)中根据已解缠绕区域相位，按照结合质量引导的局域增长方案和多 项式曲面拟合方法对残余像素解缠绕的方法具体为：
[0051] 步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素 中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个， 那么从这些生长点中选取LD-PP值最小的像素为生长点。
[0052] 步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有 不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19X19。 [0053]步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合 并为新的已解缠绕区域，具体为:在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成 一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合 窗内，假设有P个已解缠绕像素：
[0054]
[7]
[0055]上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个PX (M+1)(N+1) 大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x，y)的行向量为 [1，^7，^，...，/， 7'/71肩和~分别是沿着1和7方向的多项式指数。）。6是含有多项式系 数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵匕然后生长点 的解缠绕相位可以估计为：
[0056]
[8]
[0057] 上式中，η )表示生长点相位估计
是由生长点的多项式基构成。因此生长点的整数补偿&可以计算为：
[0058]
[9]〇
[0059] 上述步骤(8)中对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重 复步骤(3)、（4)、（5)、（6)和(7)从新解缠绕。在重复步骤3中，计算向量相位偏差和缠绕相位 偏差时，我们是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。
[0060] 上述步骤(9)中把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解 缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这 可以表5?.
[0061: [1:0]
[0062] 上式中，％'表示了添加滤波前后缠绕相位差的缠绕相位，Am代表滤波后缠绕相 位，Φ Νυ?表示对滤波后缠绕相位解缠绕后的结果。在重复步骤(3)中，计算向量相位偏差和 缠绕相位偏差时，第一次我们根据的?%，第二次我们是根据前一循环计算出解缠绕相位，而 不再是滤波后的缠绕相位。
[0063] 本发明的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括如 下步骤：（1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行 非局部均值滤波以获得滤波后的复向量；（2)根据滤波后复向量计算出滤波后的缠绕相位 图，并计算滤波前缠绕相位与滤波后缠绕相位的差；（3)根据滤波后缠绕相位图，计算向量 相位局部偏差和缠绕相位局部偏差；（4)根据计算的向量相位局部偏差和缠绕相位局部偏 差估计缠绕附近像素个数，并根据极点定义计算掩模区内极点个数；（5)根据计算的缠绕相 位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕附近像素个数和极点个数获得 掩模区内像素分类阈值；（6)根据计算的像素分类阈值对滤波后的掩模区内像素分类:相位 平滑无缠绕的不连通块和残余像素；（7)根据局部多项式曲面拟合方法解缠绕合并所有不 连通块，然后根据已解缠绕区域的相位值，按照结合质量引导的局域增长方案和多项式曲 面拟合方法对残余像素解缠绕；（8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然 后依次重复步骤(3)、（4)、（5)、（6)和(7); (9)把步骤(2)中得到的滤波前后缠绕相位差添加 到步骤(8)得到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍 为缠绕相位图像，然后依次重复步骤(8)解缠绕两次。因此，针对现有技术不足，提出方法利 用修改的复数非局部均值滤波器特点，缠绕在LD-PP与LD-WP中的不同表现和极点补偿的方 式自动获得像素分类阈值，有效地保证块内相位信息的平滑，从而保证块内相位解缠绕的 可靠性。首先将图像中平滑的块进行解缠绕后，再用逐像素区域增长方法对靠近变化剧烈 区域的残余像素解缠绕，可以综合分块与逐像素区域增长两种方法的优点，有助于实现更 稳定准确的相位估计。本技术通过仿真实验与真实磁共振水脂分离实验对所提方法的性能 进行了可靠性验证。
【附图说明】
[0064]利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限 制。
[0065] 图1为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中计算得 到的幅值图像，滤波前后缠绕相位图像，滤波前后缠绕相位差，LD-PP，LD-WP，LD-WP与LD-PP 之差和像素分类结果；
[0066] 图2为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法的流程 图；
[0067] 图3为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用具 有不同信噪比和不同快速相位变化的仿真数据集1来比较验证提出方法；
[0068] 图4为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用具 有相同信噪比但不同快速相位变化的仿真数据集2来判断提出方法对具有相同信噪比但不 同快速相位变换的数据解缠绕结果示意图；
[0069] 图5为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真 实人体膝关节横断面磁共振数据判断提出方法应用在三点Dixon水脂分离技术中的效果示 意图；
[0070] 图6为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真 实人体膝关节矢状位磁共振数据判断提出方法应用在三点Dixon水脂分离技术的效果示意 图；
[0071] 图7为本发明一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法中使用真 实人体踝关节矢状位磁共振数据判断提出方法应用在三点Dixon水脂分离技术中的效果示 意图；
【具体实施方式】
[0072]下面结合具体的实施例对本发明进行详细描述。
[0073] 实施例1
[0074] -种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括如下步骤：
[0075] (1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行 非局部均值滤波以获得滤波后的复向量；
[0076] 滤波前后复向量采用如下方式计算获得：
[0077] 滤波前复向量$可由获取的复数磁共振数据除以其幅值获得。滤波后复向量/丨〃 由滤波后的复数磁共振数据除以相应幅值获得。其中非局部均值滤波参数分别是:t = 3;f = l;k=1.0;噪声标准差为0.3,以相同滤波参数同时对磁共振数据的实部和虚部进行滤 波。
[0078] (2)根据滤波后复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤 波后缠绕相位的差；
[0079] 计算出滤波后的缠绕相位Ziw，滤波前缠绕相位口与滤波后缠绕相位的相位 差具体是：
[0080] 根据计算出的滤波前后复向量戽和，分别对其做缠绕操作使得每个像素相位 值都在区间(-^]内，从而获得滤波前后的缠绕相位图像P和;
[0081] 根据前面计算的滤波前后缠绕相位图像供和^，计算由滤波导致的缠绕相位变
[0082] (3)根据滤波后缠绕相位图，计算向量相位的局部偏差和缠绕相位的局部偏差； [0083]计算向量相位的局部偏差(LD-PP)和缠绕相位的局部偏差(LD-WP)具体是：
[0084]根据复向量计算出的向量相位局部偏差LD-PP能把掩模内具有真实相位快速变换 特点的像素更直观表现出来，即在LD-PP图中，具有较大值的像素代表其所在的区域真实相
位变4丨7*六丄+| T F\ r>r>^1=TT 々曰
[0085： [1]
[0086] 上式中，N( XO，yo)表示以空间位置（XO，yo)为中心的八邻域像素空间坐标集合， Α(λ%.ν)为与空间位置(Xhyo)相邻的位于( X，y)处像素的复向量，/；(XQ，.VQ)为位于(XQ，y 0) 处的复向量。为P/X，.V)的复共辄。
[0087] 由缠绕相位计算出的缠绕相位局部偏差LD-WP不仅受到真实相位快速变化的影 响，同时也受到缠绕影响，即在LD-WP图中，如果某像素的LD-WP值较大，这可能是由真实相 位变化较快或缠绕引起的。LD-WP图可由下式获得：
[0088]
[2]
[0089] 上式中，利χ,.ν)为与空间位置(XQ，yo)相邻的位于(x，y)处的缠绕相位，供(4，>'。）为 位于(XQ，yo)的缠绕相位。
[0090] (4)根据计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数， 并根据极点定义计算掩模区内极点个数；
[0091] 根据计算向量相位局部偏差LD-PP和缠绕相位局部偏差LD-WP估计缠绕附近像素 个数和极点个数具体是：
[0092] Mw= I (LD-WP)-(LD-PP) | >0 [3]
[0093] 上式中，MW表示缠绕附近像素个数。
[0094] Np= I Residues I >0 [4]
[0095]上式中，Np代表极点个数，Residues表示根据极点定义计算出的所有极点。
[0096] (5)根据计算的缠绕相位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕 附近像素个数和极点个数获得掩模区内像素分类阈值；
[0097] 根据缠绕相位局部偏差值的大小，求得掩模区内像素分类阈值具体是：
[0098]对掩模内缠绕像素按照其LD-WP值的大小进行降序排列，选择靠近(Mw+A*Np)处的 某个像素的LD-WP值作为像素分类阈值Thre，如下式：
[0099] Thre Se = round (Mw+A*Np) [5]
[0100] 上式中，A的取值为I，r〇und(Z)表示一个计算距离Z最近整数的函数，Threse表示相 位分类阈值的位置。
[0101] (6)根据计算的像素分类阈值对滤波后的掩模区内像素分类:相位平滑无缠绕的 不连通块和残余像素；
[0102] 根据像素分类阈值Thre对掩模区内像素分类:相位平滑无缠绕的不连通块和残余 像素；
[0103] Blocks = (LD-ffP) <Thre [6]
[0104] 上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连 通块被重新分类为残余像素。如果某个不是面积最大的不连通块距离其他块的最小边沿欧 式距离都大于20个像素，那么这个块被重新归类为残余像素。
[0105] (7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕 区域的相位值，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠 绕；
[0106] 其中利用局部多项式曲面拟合方法进行块间相位解缠绕的方法具体为：
[0107] 步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始 块内相位已解缠绕。
[0108] 步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素 具有相同的相位补偿2k3i，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕 块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是〇. 9和0.1。距离最小 的块被选为生长块。
[0109] 步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块为新的已解缠 绕区域，具体为:我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近生长 块的已解缠绕区域中选取η个像素参加曲面拟合，从靠近已解缠绕区域的生长块中选取η个 像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的 整数区间，如[_3,3]。
[0110]步骤d，根据步骤c选取的已解缠绕区域的η个像素相位和生长块中的η个像素相位 以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优 整数补偿。
[0111]根据合并后的块，按照结合质量引导的区域增长方案和多项式曲面拟合方法对残 余像素解缠绕的方法具体为：
[0112] 步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素 中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个， 那么从这些生长点中选取LD-PP值最小的像素为生长点。
[0113] 步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有 不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19X19。
[0114] 步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合
并为新的已解缠绕区域，具体为:在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成 一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合 窗内，假I力"PAP 像素：
[0115] [7]
[0116] 上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个PX (Μ+1)(Ν+1) 大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x，y)的行向量为 [l，x，y，xy，. . .，XM，yN，xMyN]，M和N分别是沿着X和y方向的多项式指数。）。《$是含有多项式系 数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵I然后生长点 的解缠绕相位可以估计为：
[0117]
[8]
[0118] 上式中，^(?? ,u)表示生长点相位估计值，
是由生长点的多烦式某构成。闵此牛长点的整数补偿气%,Λ)可以计算为：
[0119]
[9]
[0120] (8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、 (4)、（5)、（6)和(7);
[0121] 对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、（4)、 (5 )、（6)和(7)从新解缠绕。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，我们 是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。
[0122] (9)把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位 上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重 复步骤(8)解缠绕两次。
[0123] 把步骤(2)中得到的滤波前与滤波后缠绕相位差添加到步骤(8)得到的解缠绕相 位图像上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这可 以表示为：
[0124：
[10]
[0125] 上式中，仍v表示添加滤波前与滤波后缠绕相位差的缠绕相位，外代表非局部均 值滤波后的缠绕相位，Φ?Μ表示对滤波后的缠绕相位解缠绕后的结果。在重复步骤3中，计 算向量相位偏差和缠绕相位偏差时，第一次我们根据的，第二次我们是根据前一循环计 算出解缠绕相位，而不再是滤波后的缠绕相位。
[0126] 与现有技术相比，提出方法利用修改的复数非局部均值滤波器特点，根据缠绕在 LD-PP与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相 位信息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕 后，再用逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐 像素区域增长两种方法的优点，有助于实现更稳定准确的相位估计。
[0127] 本技术通过仿真实验与真实磁共振水脂分离实验对所提方法的性能进行了实验 验证，实验结果表明提出的方法能很好实现对缠绕相位图像解缠绕，水脂分离结果中无明 显水脂互换现象存在。
[0128] 实施例2
[0129] -种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，依次包括以下步骤： [0130] (1)根据采集的复数磁共振数据计算其复向量，并对采集的复数磁共振数据进行 非局部均值滤波以获得滤波后的复向量。在图1中，第一行分别显示了一个膝关节横断面磁 共振数据的幅值(Magnitude)图像，缠绕相位(Wrapped phase)图像，非局部均值滤波后获 得缠绕相位(Filtered wrapped phase)和滤波前后缠绕相位差。第二行显示了由滤波前磁 共振数据获得的缠绕相位局部偏差（LD-WP)、向量相位局部偏差（LD-PP)、缠绕附近像素 (Wraps locations)和像素分类后的不连通块(Blocks)和残余像素(Residual pixels);第 三行显示了由滤波后磁共振数据获得的缠绕相位局部偏差、向量相位局部偏差、缠绕附近 像素和像素分类后的不连通块和残余像素，
[0131] 在对掩模区像素分类中，由滤波前磁共振数据获得的缠绕附近像素和极点个数分 别是3374和268,在缠绕相位局部偏差获得的分块阈值为1.997,像素分类结果如图1所示。 由滤波后磁共振数据获得的缠绕附近像素和极点个数分别是2185和0,在缠绕相位局部偏 差图获得的分块阈值为〇. 5679，分块结果如图1所示。
[0132] 因此，利用非局部均值滤波能有效降低极点对像素分类的影响，获得较好的像素 分类结果。
[0133] (2)根据滤波后复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相位与滤 波后缠绕相位的差。利用步骤（1)中计算的滤波前后复向量，再利用Matlab函数库中angle 0函数，可以分别获得滤波前后缠绕相位及滤波前后缠绕相位差(请参阅图1)。
[0134] (3)根据滤波缠绕相位图，计算向量相位局部偏差和缠绕相位局部偏差。
[0135] 分别利用步骤(1)和(2)中计算的滤波后复向量和缠绕相位，可以分别获得向量相 位局部偏差(LD-PP)和缠绕相位局部偏差(LD-WP)(请参阅图1)。其中计算向量相位局部偏 差和缠绕相位局部偏差具体是：
[0136] 根据复向量计算出的向量相位局部偏差LD-PP能把掩模内具有真实相位快速变换 特点的像素更直观表现出来，即在LD-PP图中，具有较大值的像素代表其所在的区域真实相 位变化较快，LD-PP可由下式获得：
[0137]
[1]
[0138] 上式中，N(XQ，yQ)表示以空间位置（XQ，yQ)为中心的八邻域像素空间坐标集合， A(U')为与空间位置(x Q，yQ)相邻的位于(x，y)处像素的复向量，AU5,凡)为位于(x〇，y〇)处 的复向量。为A(U')的复共辄。如图1所示，LD-PP的较大值发生在真实相位快速变 换区域。
[0139] 由缠绕相位计算出的缠绕相位局部偏差LD-WP不仅受到真实相位快速变化的影 响，同时也受到缠绕影响，即在LD-WP图中，如果某像素的LD-WP值较大，这可能是由真实相 位变化较快或缠绕引起的。LD-WP图可由下式获得：
[0140]
[2]
[0141] 上式中，列.V)为与空间位置(XQ，yο)相邻的位于(X，y)处的缠绕相位，舛, .V。)为 位于(XQ，yo)的缠绕相位。如图1所示，LD-WP的较大值出现在缠绕发生区域和/或相位快速变 换区域。
[0142] 因此，在本实施方案中，根据缠绕在计算LD-WP和LD-PP中的不同表现，可以有效的 估计缠绕附近像素的位置。
[0143] (4)根据计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数， 并根据极点定义计算掩模区内极点个数。
[0144] 根据步骤(3)计算的向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个 数和极点个数具体是：
[0145] Mw= I (LD-WP)-(LD-PP) | >0 [3]
[0146] 上式中，]\^表示缠绕附近像素个数。在本例中]\^为2185(请参阅图1)。
[0147] Np= I Residues I >0 [4]
[0148]上式中，Np代表极点个数，Res idues表示根据极点定义计算出的所有极点。在本例 中^为0(请参阅图1)。
[0149] (5)根据计算的缠绕相位局部偏差值，使得掩模区内像素降序排列，然后根据缠绕 附近像素个数和极点个数获得掩模区内像素分类阈值。
[0150] 根据步骤(3)计算的缠绕相位局部偏差值的大小，求得掩模区内像素分类阈值具 体是：
[0151] 对掩模内缠绕像素按照其缠绕相位局部偏差值的大小进行降序排列，选择最接近 序号(Mw+A*Np)处的某个像素的LD-WP值作为像素分类阈值Thre，如下式：
[0152] Threse=Mw+Np [5]
[0153] 上式中，A的取值为l，r〇und(Z)表示一个计算距离Z最近整数的函数，Threse表示相 位分类阈值的位置。在本例中Thres e和Thre分别等于2185和0.5679(请参阅图1)。
[0154] (6)根据计算的像素分类阈值对滤波后的掩模区内像素分类:相位平滑无缠绕的 不连通块和残余像素。
[0155] 根据步骤(5)计算的像素分类阈值Thre对掩模区内像素分类，具体如下式：
[0156] Blocks = (LD-ffP) <Thre [6]
[0157] 上式中，Blocks表示无缠绕的不连通块。为了便于计算，面积小于10个像素的不连 通块被重新分类为残余像素。如果某个不是面积最大的不连通块距离其他块的最小边沿欧 式距离都大于20个像素，那么这个块被重新归类为残余像素。
[0158] (7)根据局部多项式曲面拟合方法解缠绕合并所有不连通块，然后根据已解缠绕 区域的相位值，按照结合质量引导的局域增长方案和多项式曲面拟合方法对残余像素解缠 绕。
[0159] 其中利用局部多项式曲面拟合方法进行块间相位解缠绕的方法具体为：
[0160]步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始 块内相位已解缠绕。
[0161] 步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素 具有相同的相位补偿2k3i，k是整数。在确定生长块之前，计算已解缠绕区域与所有未解缠绕 块之间的边沿欧式距离与中心欧式距离之和，这两类距离权重分别是〇. 9和0.1。距离最小 的块被选为生长块。
[0162] 步骤c，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块为新的已解缠 绕区域，具体为:我们假设已解缠绕区域内相位可以拟合成一平滑多项式曲面。从靠近生长 块的已接缠绕区域中选取η个像素参加曲面拟合，从靠近已接缠绕区域的生长块中选取η个 像素参加曲面拟合，如100个像素。同时假设生长块的整数补偿k的选取范围为一个较小的 整数区间，如[_3,3]。
[0163] 步骤d，根据步骤c选取的已解缠绕区域的η个像素相位和生长块中的η个像素相位 以及位于设定区间内的整数k进行曲面拟合，使得多项式曲面拟合误差最小的k被选为最优 整数补偿。
[0164] 根据合并后的块，按照结合质量引导的区域增长方案和多项式曲面拟合方法对残 余像素解缠绕的方法具体为：
[0165] 步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素 中最靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个， 那么从这些生长点中选取LD-PP值最小的像素为生长点。
[0166] 步骤b，选择残余像素解缠绕曲面拟合窗的大小。不同缠绕复杂度对曲面拟合窗有 不同的要求，如在所有真实数据实验中，统一设定残余像素解缠绕的拟合窗为19X19。
[0167] 步骤c，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点为合 并为新的已解缠绕区域，具体为:在拟合窗内，假设已解缠绕像素的真实相位值可以拟合成 一平滑多项式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计。在拟合 窗内，假设有P个已解缠绕像素：
[0168] Φ = Jf 十茗 [7]
[0169] 上式中，Φ是由这P个已解缠绕像素相位值构成的列向量，X是一个PX (M+1)(N+1) 大小的矩阵，每一行是由每一个已解缠绕像素的多项式基(如在空间位置(x，y)的行向量为 [l，x，y，xy，. . .，XM，yN，xMyN]，M和N分别是沿着X和y方向的多项式指数。）。6是含有多项式系 数的列向量，E是含有残余误差的列向量。使用最小二乘方法可以求系数矩阵θ，然后生长点 的解缠绕相位可以估计为：
[0170]
[8]
[0171] 上式中，表示生长点相位估计值
是由生长点的多项式基构成。因此生长点的整数补偿可以计算为：
[0172]
[9]
[0173] 从执行残余像素解缠绕步骤后的解缠绕结果可以看出，存在缠绕的残余像素明显 降低(请参阅图2)。
[0174] (8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、 (4)、（5)、（6)和(7)。
[0175] 对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤(3)、（4)、 (5 )、（6)和(7)从新解缠绕。在重复步骤(3)中，计算向量相位偏差和缠绕相位偏差时，我们 是根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。
[0176] (9)把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位 上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，然后依次重 复步骤(8)解缠绕两次。
[0177] 把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到的解缠绕相位图像 上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图像，这可以表示 为：
[0178] [10]
[0179] 上式中，表示添加滤波前与滤波后缠绕相位差的缠绕相位，代表非局部均 值滤波后的缠绕相位，表示对滤波后的缠绕相位解缠绕后的结果。在重复步骤(3)中， 计算向量相位偏差和缠绕相位偏差时，第一次我们根据的，第二次我们是根据前一循环 计算出解缠绕相位，而不再是滤波后的缠绕相位。
[0180] 从对滤波前缠绕相位图像的最终解缠绕结果和水脂分离情况可以看出，解缠绕结 果很平滑，无论是在快速相位变化还是低信噪比区域都无明显的缠绕像素存在，从而在水 脂分离图像中无明显水脂互换现象存在(请参阅图5)。
[0181 ]与现有技术相比，提出方法利用修改的复数非局部均值滤波器特点，缠绕在LD-PP 与LD-WP中的不同表现和极点补偿的方式自动获得像素分类阈值，有效地保证块内相位信 息的平滑，从而保证块内相位解缠绕的可靠性。首先将图像中平滑的块进行解缠绕后，再用 逐像素区域增长方法对靠近变化剧烈区域的残余像素解缠绕，可以综合分块与逐像素区域 增长两种方法的优点，有助于实现更稳定准确的相位估计。
[0182] 本技术通过仿真实验与真实水脂分离磁共振成像实验对所提方法的性能进行了 实验验证，实验结果表明提出方法能很好实现对缠绕相位图像解缠绕，水脂分离结果中无 明显水脂互换现象存在。
[0183] 实施例3。
[0184] 通过本发明的方法对仿真和真实磁共振相位进行相位解缠绕实验，本实施例选取 其中部分实验结果进行分析比较。
[0185] 图3为具有信噪比从6等间隔变化到0.5的圆型仿真数据的相位解缠绕实验结果。 未添加噪声的仿真数据真实相位为高斯相位曲面，相位变化区间分别为O到25和50弧度 (rad)，即高度(Height)等于25和50弧度(rad)。添加均值为0,标准差为25随机噪声到的生 成的复数仿真数据上以生成具有不同信噪比的复数仿真数据。其中第一行展示了幅值 (magnitude)图像;第二行分别展示了高度为25rad的相关数据，分别是缠绕相位(Wrapped phase)图像，PRELUDE方法解缠绕结果，PRELUDE解缠绕结果与真实相位和添加噪声的差 (Difference)，提出方法的(Proposed method)解缠绕结果，提出方法的解缠绕结果与真实 相位和添加噪声的差(Difference);第三行分别展示了高度为50rad的相关数据，分别是缠 绕相位图像，PRELUDE方法解缠绕结果，PRELUDE解缠绕结果与真实相位和添加噪声的差，提 出方法解缠绕结果，提出方法解缠绕结果与真实相位和添加噪声的差。
[0186] 比较两种方法的解缠绕结果，可以看出：采用现有经典算法PRELUDE的解缠绕结果 错误像素较多，且有成块的错误出现，而采用本发明算法的解缠绕结果大大改善了这一情 况。
[0187] 图3为两种算法对具有相同信噪比但不同快速相位变化的仿真数据解缠绕结果。 未添加噪声前真实相位变化区间分别为0到50、100、150和200rad，添加均值为0,标准差为 10的随机噪声到复数仿真数据上。其中，第一列最上面图像为幅值图像，其余的为高度等于 50、100、150和200rad的缠绕相位图像；第二列为PRELUDE方法解缠绕结果；第三列为 PRELUDE解缠绕结果与真实相位和添加噪声的差;第四列为提出方法解缠绕结果;第五列为 提出方法的解缠绕结果与真实相位和添加噪声的差。
[0188]可见，PRELUDE解缠绕结果（图4第二和第三列）随着缠绕复杂度的增加，错误像素 逐渐增加，而当高度大于等于150rad后，在中心区域出现了明显的缠绕块。而在本发明算法 的解缠绕结果中，虽然也有错误像素的出现，但是当高度等于150rad后，无明显的错误块出 现（图4第四和第五列，与PRELUDE解缠绕结果相比，更低的相对误差也直观上说明了本发明 的方法解缠绕结果更加准确。
[0189]图5、6和7为两种算法分别对人体不同部位(膝关节横断面、膝关节矢状面和踝关 节矢状面)的相位解缠绕结果及应用到三点Dixon技术中进行水脂分离的结果。在图4、5和6 中，第一行分别为缠绕相位(Wrapped phase)，反相位幅值(Out-phase magnitude)和同相 位幅值（In-phase magnitude);第二行分别是PRELUDE相位解缠绕结果，和使用PRELUDE相 位解缠绕结果的水脂分离结果，包括水图(Water)，脂肪图(Fat)及脂肪分数图（FF);第三行 分别是提出方法相位解缠绕结果，和使用提出方法相位解缠绕结果的水脂分离结果，包括 水图，脂肪图及脂肪分数图。
[0190] 可见，PRELUDE在三种数据类型的解缠绕结果中（图6、7中第二行)都有缠绕残留， 而残留的缠绕会导致水脂分离结果中出现了明显的互换现象（图6、7中第二行白色箭头 处）。而在本发明算法的解缠绕结果中，而虽然偶尔也有缠绕残留（图7中第二行白色箭头 处），但与PRELUDE结果相比，互换面积明显减少。对掩模区域的脂肪分数计算无明显影响。 因此，真实磁共振数据也说明了本发明的方法解缠绕结果更加准确。
[0191]通过大量仿真和真实数据实验证明，本发明的方法能够能利用缠绕在向量相位偏 差和缠绕相位偏差的不同表现和修改的非局部均值滤波器实现掩模区内像素的有效分类， 而利用局部多项式曲面拟合方法能有效的合并平滑的不连通块，最后再利用质量引导的区 域增长和局部多项式曲面拟合方法能有效的实现对残余像素的真实相位估计，从而获得平 滑解缠绕结果。
[0192]最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护 范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理 解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和 范围。
【主权项】
1. 一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法，其特征在于:包括如下 步骤： 步骤(1)根据采集的复数磁共振数据计算其复向量馬，并对采集的复数磁共振数据进行 非局部均值滤波W获得滤波后的复向量? 步骤(2)根据滤波后的复向量计算出滤波后的缠绕相位图，并计算滤波前缠绕相 位与滤波后缠绕相位的差； 步骤(3)根据滤波后的缠绕相位图，计算每个像素的向量相位局部偏差和缠绕相位局 部偏差； 步骤(4)根据向量相位局部偏差和缠绕相位局部偏差估计缠绕附近像素个数，并根据 极点定义计算掩模内极点个数； 步骤(5)根据每个像素的缠绕相位局部偏差值，将掩模内像素降序排列，然后根据缠绕 附近像素个数和极点个数获得掩模内像素分类阔值； 步骤(6)根据像素分类阔值对滤波后的掩模内的像素分类:一类为相位平滑无缠绕的 不连通块，一类为残余像素； 步骤(7)根据局部多项式曲面拟合方法解缠绕并合并所有不连通块，然后根据已解缠 绕区域的相位，按照结合质量引导的局域增长方法和多项式曲面拟合方法对残余像素解缠 绕； 步骤（8)把解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次重复步骤 (3)--(7)； 步骤(9)把步骤(2)中得到的滤波前缠绕相位与滤波后缠绕相位的差添加到步骤(8)得 到的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图 像，然后依次重复步骤(8)迭代解缠绕两次。2. 根据权利要求1所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:步骤（1)中利用非局部均值滤波器W相同参数同时对输入的磁共振数据实部 和虚部滤波，获得滤波后的复向量=cxp(i/χv^"")，上式中口》w表示滤波后缠绕相位， 其滤波参数分别是:t = 5;f = 3;k=l .2,噪声标准差为0.3，i X i = -l。3. 根据权利要求2所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(2)中计算出滤波后的缠绕相位，并计算滤波前缠绕相位与滤波后缠 绕相位的相位差具体是： 根据滤波后的复向量对其做缠绕操作使得每个像素相位值都在区间（-η，π]内， 从而获得滤波后的缠绕相位图像，同理获得滤波前的缠绕相位图像0; 根据前面计算的滤波前后缠绕相位图像，计算由滤波导致的缠绕相位值的变化即滤波 前缠绕相位与滤波后缠绕相位的相位差。4. 根据权利要求3所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(3)中计算向量相位局部偏差和缠绕相位局部偏差具体是： 根据复向量计算出的向量相位局部偏差能把掩模内具有真实相位快速变换特点的像 素更直观表现出来，具有较大向量相位局部偏差值的像素代表其所在的区域真实相位变化 较快，因此，向量相位局部偏差可由下式获得：[1] 上式中，N(x〇，y〇)表示W空间位置(x〇，y〇)为中屯、的八邻域像素空间坐标集合，/；片，>'） 为与空间位置(x〇，y〇)相邻的位于(X，y)处像素的复向量，ッ?，.v。；!为位于(χo，yo)处的复向 量。C化.V')为沁y)的复共辆。； 由缠绕相位计算出的缠绕相位局部偏差不仅受到真实相位快速变化的影响，同时也受 到缠绕影响，某像素的缠绕相位局部偏差值较大可能是由真实相位变化较快或缠绕引起 的，缠绕相位局部偏差可由下式获得：口] 上式中，LD-WP表示缠绕相位局部偏差，口 (Xv)为W空间位置(x〇，y〇)为中屯、的位于(X， y)处的缠绕相位，抑而，乃）为位于(X0，yo)的缠绕相位。5. 根据权利要求4所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(4)中根据计算向量相位局部偏差和缠绕相位局部偏差估计缠绕附 近像素个数和极点个数具体是： Mw=I(LD-WP)-(LD-PP)|>0 [3] 上式中，Mw表示缠绕附近像素个数，LD-WP和LD-PP分别是缠绕相位偏差和向量相位局部 偏差； Np= I Residues I >0 [4] 上式中，Np代表极点个数，Residues表示根据极点定义计算出的所有极点。6. 根据权利要求5所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(5)中根据获得的缠绕相位局部偏差值，求得掩模区内像素分类阔值 具体是： 根据每个像素的缠绕相位局部偏差值，将掩模内像素降序排列，然后根据缠绕附近像 素个数和极点个数获得掩模内像素分类阔值； 对掩模内缠绕像素按照其缠绕相位局部偏差值进行降序排列，选择靠近序号(Mw+A*Np) 处的某个像素缠绕相位局部偏差值作为像素分类阔值化re，如下式： Τ'虹 ese = :round(Mw+A*Np) [5] 上式中，A的取值为l，round(Z)表示一个计算距离Z最近整数的函数，Threse表示相位分 类阔值的位置。7. 根据权利要求6所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(6)根据获得的像素分类阔值化re对掩模区内像素分类:相位平滑无 缠绕的不连通块和残余像素； Blocks=(LD-WP)yrhre [6] 上式中，Blocks表示无缠绕的不连通块，面积小于10个像素的不连通块被重新分类为 残余像素，同时，如果不是面积最大的某不连通块距离其他块的边沿欧式距离都大于20个 像素，那么运个块内像素也将被归类残余像素。8.根据权利要求7所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(7)中利用局部多项式曲面拟合方法进行块间相位解缠绕并合并的 方法具体为： 步骤a，选取起始块，面积最大的块，即像素个数最多的块被选为起始块，认为起始块内 相位已解缠绕； 步骤b，选择生长块，下一个将要被解缠绕块被命名为生长块，生长块内所有像素具有 相同的相位补偿化n，k是整数，在确定生长块之前，计算已解缠绕区域与所有未解缠绕块之 间的边沿欧式距离与中屯、欧式距离之和，运两类距离权重分别是0.9和0.1，距离最小的块 被选为生长块； 步骤C，根据局部多项式曲面拟合方法合并已解缠绕区域和生长块，并作为新的已解缠 绕区域，具体为:已解缠绕区域内相位可W拟合成一个平滑的多项式曲面，从靠近生长块的 已解缠绕区域中选取η个像素参加曲面拟合，从靠近已解缠绕区域的生长块中选取η个像素 参加曲面拟合，同时假设生长块的整数补偿k的选取范围为一个较小的整数区间； 步骤d，根据步骤C选取的已解缠绕区域中的η个像素相位和生长块中的η个像素相位W 及位于设定区间内的整数k进行曲面拟合，使多项式曲面拟合误差最小的k被选为最优整数 补偿； 所述步骤(7)中根据已解缠绕区域相位，按照结合质量引导的区域增长方案和多项式 曲面拟合方法对残余像素解缠绕的方法具体为： 步骤a，选择生长点，生长点选取根据两个标准，第一个是选取待解缠绕残余像素中最 靠近已解缠绕区域的残余像素为生长点；第二个是如果标准一选取的生长点有多个，从运 些生长点中选取LD-PP值最小的像素为生长点； 步骤b，选择残余像素解缠绕曲面拟合窗的大小； 步骤C，利用局部多项式曲面拟合方法把已解缠绕区域和解缠绕后的生长点合并为新 的已解缠绕区域，具体为:在拟合窗内，已解缠绕像素的真实相位值可W拟合成一平滑多项 式曲面，待解缠绕像素相位值可由拟合的局部多项式曲面去接近估计，在拟合窗内，有P个 已解缠绕像素：[7] 上式中，Φ是由运P个已解缠绕像素相位值构成的列向量，X是一个PX (M+1KN+1)大小 的矩阵，每一行是由每一个已解缠绕像素的多项式基，如在空间位置(x，y)的行向量为[1， x，y，xy，. . .，χM，/，χMy^^]，M和N分别是沿着χ和y方向的多项式指数;c^是含有多项式系数的 列向量，E是含有残余误差的列向量;使用最小二乘方法可W求系数矩阵含，然后生长点的解 缠绕相位估计为：[8] 上式中，<i。，如表示生长点相位估计值，义(V.V。)=化為，拘，為说，…，而M，沁V，而 生长点的多项式基构成，因此生长点的整数补偿^^<^^.，.",可^计算为：脚。9. 根据权利要求8所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(9)中把步骤(2)中得到的滤波前后缠绕相位差添加到步骤(8)得到 的解缠绕相位上，并认为添加了滤波前后缠绕相位差的解缠绕相位图像仍为缠绕相位图 像，运可W表示为：閲 上式中，巧V表示了添加滤波前后缠绕相位差的缠绕相位，稱&W代表滤波后缠绕相位， Φνιμ表示对滤波后缠绕相位解缠绕后的结果;在重复步骤(3)中，计算向量相位偏差和缠绕 相位偏差时，第一次我们根据的嘶，第二次我们是根据前一循环计算出解缠绕相位，而不再 是滤波后的缠绕相位。10. 根据权利要求8所述的一种基于缠绕识别和局部曲面拟合磁共振相位解缠绕方法， 其特征在于:所述步骤(8)中对解缠绕后相位图像认为是仍存在缠绕的相位图像，然后依次 重复步骤(3)、（4)、（5)、（6)和(7)从新解缠绕，在重复步骤(3)中，计算向量相位偏差和缠绕 相位偏差时，根据前一循环计算出解缠绕相位，而不再是滤波后的缠绕相位。
【文档编号】G06T7/00GK105844626SQ201610159078
【公开日】2016年8月10日
【申请日】2016年3月18日
【发明人】程军营, 刘晓云, 陈武凡, 王常青
【申请人】电子科技大学