Print

发布时间: 2018-06-16
摘要点击次数:
全文下载次数:
DOI: 10.11834/jig.170515
2018 | Volume 23 | Number 6




    医学图像处理    




  <<上一篇 




  下一篇>> 





相位分区和局部多项式曲面拟合的相位解缠绕
expand article info 程军营1,2, 王常青1,2, 冯衍秋1, 陈武凡1,2
1. 电子科技大学自动化工程学院, 成都 611731;
2. 南方医科大学生物医学工程学院, 广州 510515

摘要

目的 准确的解缠绕相位是两点或三点Dixon技术等在磁共振临床应用的前提和关键,然而当相位图像中存在严重噪声、快速相位变换或不连通区域时,当前许多已经提出的相位解缠绕算法将会失败。为此本文提出一种基于相位分区和局部多项式曲面拟合的相位解缠绕新方法,该新方法在相位图像存在严重噪声、快速相位变换或不连通区域的情况下仍可以稳定可靠的工作。方法 首先将获得的相位图像分成连通块,块内相位都在给定的相位区间内,把像素个数小于给定阈值的块归类为残余像素。然后利用局域增长的局部多项式曲面拟合方法依次进行块与块之间相位解缠绕和残余像素相位解缠绕。最后使用仿真与真实磁共振Dixon数据来评价提出方法,并与PRELUDE(phase region expanding labeler for unwrapping discrete estimates)方法进行了比较。结果 在不同信噪比、快速相位变换或存在不连通区域的仿真实验中,即使当数据中存在信噪比为0.5、相邻相位变换大于π弧度或不连通区域时,提出方法的平均错误率不大于0.51%。对于100层真实的磁共振膝关节和踝关节水与脂肪分离图像,提出方法生成结果中发生明显解缠绕错误及水与脂肪互换的比率为6.00%,而PRELUDE却为42.00%。结论 本文提出了一种磁共振相位解缠绕算法,利用相位分区方法,可靠的实现相位图像分块;利用局部多项式曲面拟合方法,准确的实现相位解缠绕。提出方法能够更加鲁棒的实现相位解缠绕,这将有益于相位相关的磁共振临床的应用,如两点和三点Dixon水脂分离技术、磁敏感加权成像和人脑相位成像等。

关键词

磁共振成像; 相位解缠绕; 相位分区; 局部多项式曲面拟合; 水与脂肪分离

New phase-unwrapping method based on phase partition and local polynomial surface fitting
expand article info Cheng Junying1,2, Wang Changqing1,2, Feng Yanqiu1, Chen Wufan1,2
1. School of Automation Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China;
2. School of Biomedical Engineering, Southern Medical University, Guangzhou 510515, China
Supported by: National Natural Science Foundation of China (61471188)

Abstract

Objective An accurate phase map is crucial for magnetic resonance imaging (MRI)-based clinical applications, such as two-point and three-point Dixon techniques, susceptibility-weighted imaging, and human-brain phase imaging. However, the phase calculated from the complex MR signal is commonly wrapped back into the range of (-π, π] for the arctangent operation. Phase unwrapping is required to recover the underlying actual phase. A certain number of phase unwrapping methods have been proposed under the assumption that the underlying actual phase difference between adjacent pixels should not be larger than π. If the wrapped phase map contains severe noise, rapid phase change, or even disconnected regions in the interesting region, then the underlying actual phase difference between adjacent pixels may be larger than π; thus, this assumption becomes invalid. A new phase-unwrapping method is proposed; this method is based on phase partition and local polynomial surface fitting, which works robustly under the situation of severe noise, rapid phase changes, and disconnected regions. Method The proposed method first exploits the phase partition method to split the acquired phase map into the connected blocks. The phase inside each block remains within a given interval, and the wrapped phase difference between adjacent pixels is less than π. Therefore, no phase wrap exists inside each block. To reduce the affection of noise, the blocks with pixel number of less than a given threshold (such as 50) are clustered into residual pixels. Then, the proposed method sequentially performs inter-block phase unwrapping and residual-pixel phase unwrapping by a region-growing local polynomial surface fitting approach. Three simulated datasets were separately generated with signal-to-noise ratios (SNRs) varying from 6.5 to 0.5, phase height changing from 100 rad to 200 rad, and disconnected regions to test the performance of the proposed method on the data under the situation of severe noise, different phase change levels and disconnected regions. The three-point Dixon knee and ankle data of five healthy adult human volunteers were acquired on a 0.35 T permanent magnet MR scanner to test the performance of the proposed method on the in vivo MR data (XGY-OPER, Ningbo Xingaoyi, Ningbo, China). In addition, the two-point Dixon knee data of one healthy volunteer were acquired on a 3.0 T MR scanner (Philips, Achieva, Netherlands). The simulation and in vivo two-point and three-point water-fat Dixon data were used to evaluate the proposed method and compared with phase-region expanding labeler for unwrapping discrete estimates (PRELUDE). The unwrapped error ratio was calculated as the number of inaccurate unwrapped pixels divided by the total number of pixels to quantitatively evaluate the performance of the proposed method. Each simulation was repeated 20 times, and the corresponding means and standard deviations (SDs) of unwrapped error ratio (%) were calculated. The unwrapped results were implemented via Dixon technique to produce water and fat images to evaluate the performance of the proposed method on in vivo data. If the phase-unwrapping method does not acquire an accurate phase map, then the water and fat images will contain swaps. The water-fat separation results of every slice were evaluated by a blinded board-certified radiologist according to the following four-point scale: 1) many swaps (slices), 2) few swaps (slices), 3) total swap slices, and 4) error ratio. Result 1) In the simulation experiment of different signal-to-noise ratio (SNR) levels, the phase map unwrapped by PRELUDE contains evident wrapping residues in low SNR regions. However, the unwrapped error is substantively reduced in the results generated by the proposed method. The proposed method consistently produced a substantially lower mean and SD of unwrapped error ratio than those of PRELUDE when SNR is below 2.5. 2) In the simulation experiment of different phase change levels, the unwrapped results produced by PRELUDE contain evident wraps when the phase height increased to 200 rad. On the contrary, the proposed method consistently produced accurate unwrapped results for different phase heights. 3) In the simulation experiment of existing discontented regions, PRELUDE generated results with serious wrapping residues, and the mean and SD of the unwrapped error ratio was 12.79±0.67 (%). The proposed method produced an accurate unwrapped phase with the mean and SD of unwrapped error ratio of 0.10±0.05 (%). 4) In the in vivo 0.35 T datasets of five volunteers and the 3.0 T dataset of one volunteer (a total of 100 slices, consisting of 75 sagittal knees and 25 sagittal ankles) water-fat separation experiments, the results generated by PRELUDE had nine times many swaps and 33 times few swaps out of 100 outputs. Meanwhile, the proposed method only produced six times few swaps. The total error ratio of PRELUDE was 42% and that of the proposed method was only 6%. Conclusion A new phase-unwrapping method, which first splits the acquired phase map into connected blocks by exploiting the phase partition method, is presented. The blocks with the pixel number of less than a given threshold are clustered into the residual pixels. Then, the proposed method sequentially performs inter-block phase unwrapping and residual-pixel phase unwrapping by using a region-growing local polynomial surface fitting approach. The simulation results demonstrate that the proposed method can achieve robust and accurate phase unwrapping even under serious noise, rapid phase changes, and disconnected regions. The application on two-point and three-point water-fat Dixon MRI data shows that the proposed method can successfully separate water and fat with few swaps. Therefore, the proposed method is beneficial to phase-related MRI applications, such as two-point and three-point Dixon techniques, susceptibility-weighted imaging, and human-brain phase imaging.

Key words

magnetic resonance imaging; phase unwrapping; phase partition; local polynomial surface fitting; water-fat separation

0 引言

从光学到雷达干涉测量技术,再到磁共振成像等学科,获得的相位都需要解缠绕后才能进一步处理。其中,对于许多医学磁共振临床应用来说,要求对其获得的相位图像解缠绕,例如:场图成像[1],化学偏移成像[2-3]和流量测量[4]等。当一个测量的相位信号仅仅在(-π, π]范围内取值以至于 ${\varphi _{{\rm{measured}}}} = {\rm{mod}}\; {\phi _{{\rm{original}}}}$ 时,缠绕就会发生。同时缠绕发生区域边缘将会引入2π跳转[5]。这时,为了重建一个空间上连续的相位图像,需要进行相位解缠绕操作来移除所获得相位图像中的2π跳转。相位解缠绕算法通常假设潜在真实相位在空间上是平滑或连续的,即相邻像素相位差不大于π[6]。如果该假设得到满足,则很容易获得真实相位图。但是当缠绕相位图像中含有严重噪声、快速相位变换或不连通区域时,相邻像素真实相位差有可能大于π。这时,准确的恢复真实相位变得具有挑战性,得到的相位解缠绕结果有可能是错误的。

在过去的三十多年里,曾经提出了很多相位解缠绕算法来改善解缠绕结果的准确性。这些算法可以分类[7]为:全局算法、局部算法和路径追踪方法。全局算法是指每一个像素的相位解缠绕依赖于整个图像相位特征,把相位解缠绕问题转化为一个全局函数的最小化求解[8]。所有已知的全局算法虽然是鲁棒的但计算代价高昂。局部算法则是指每个像素的相位解缠绕依赖于周围像素的相位特征,把整个相位图像分成一系列子部分,然后对每个子部分相位解缠绕,最后合并解缠绕后的所有子部分[9-10]。局部算法具有既保持了一定鲁棒性又降低了计算强度的优点。路径追踪算法是根据Itoh理论[11],利用相邻像素相位梯度来实现在整个相位图中的线性积分解缠绕。如果相位图中存在极点[12],则解缠绕结果依赖于积分路径。大多数路径追踪方法是通过优化积分路径来解决由极点引起的相位不一致。例如,在Goldstein枝切相位解缠绕算法[12]中,缠绕相位图像中的极点先被识别,并利用枝切连接它们,然后沿着一条不穿过任意枝切线的路径通过线性积分来计算解缠绕相位。基于局域增长方法的相位解缠绕可以被认为是另一类路径追踪方法[5]。在局域增长方法中,相位相对均匀区域中某个像素被选为起始点,然后利用已解缠绕区域的相位信息来确定增长点的位置和预测它的正确相位值[9-10]。质量引导的局域增长方法的表现依赖于种子点和生长点的选取及解缠绕方式。目前常见的质量图主要包括局部一阶相位导数方差、局部二阶相位偏导数方差和幅值信息等。在路径追踪相位解缠绕方法中,在积分或局域增长路径的任何位置发生的错误都有可能传递到后面的待解缠绕像素并累积。

基于分块与合并策略的相位解缠绕算法是一类高效的方法[13]。该类方法把整个缠绕相位图分成块并执行块内解缠绕,然后执行块与块之间的相位解缠绕来合并所有已经解缠绕的块。该类方法中最典型的是PRELUDE (phase region expanding labeler for unwrapping discrete estimates)[14]。在该方法中,首先,设置相位分区(例如{(-π, -2π/3] rad, (-2π/3, -π/3] rad,…, (2π/3, π/] rad }),将相位空间数据分成多个块,并认为每个块内无相位缠绕。每个块作为一个整体,计算交界的两两块的代价函数以决定哪两个块优先被合并融合。然后将最佳匹配的两个块优先解缠绕合并作为一个新的块,并重新计算交界的两两块的代价函数,重复这一步,直至只剩下一个块。PRELUDE假设感兴趣区域中相邻像素的真实相位变化不大于π。当存在相邻像素真实相位变化大于π的情况时,PRELUDE会失败。而现实采集的磁共振成像数据中,通常因含有噪声、快速相位变换或不连通区域而破坏了该假设,从而导致PRELUDE解缠绕结果含有缠绕残余。

相位解缠绕是基于相邻像素真实相位变换不大于π的假设,因此使用相位分区方法对相位图像进行分块,块内相邻像素的相位变化小于给定相位间隔。如设置相位分区为{(-π, -2π/3] rad, (-2π/3, -π/3] rad,…, (2π/3, π/] rad },则任一块内相邻像素的缠绕相位变化必小于π/3。因此块内不存在相位缠绕的像素。而多项式拟合方法对噪声不太敏感,在局域增长解缠绕过程中,利用临近的已解缠绕像素的相位信息对生长点解缠绕更能准确的估计出其真实相位值。另外,局部多项式曲面模型更能精确的描述真实相位的复杂状况和快速变化。

为此,本文提出了一种基于相位分区和局部多项式曲面拟合的相位解缠绕新方法。具体来说,利用相位分区方法把获得的相位图像分成多个连通块,块内像素相位都在给定的相位区间内,块内无相位缠绕的像素存在。为了降低噪声的影响,把像素个数小于给定阈值的块归类为残余像素。然后利用局域增长的局部多项式曲面拟合方法依次进行块与块之间相位解缠绕和残余像素相位解缠绕。本文使用仿真和真实水与脂肪Dixon磁共振成像实验来评价提出方法,并与PRELUDE进行了对比。

1 方法

1.1 相位解缠绕问题

在磁共振成像中,获得的复数信号 $S$ 可以写成 $S = M \times {{\rm{e}}^{i\phi }}$ $M$ 代表幅值信息, $\phi $ 为真实相位。因此缠绕相位 $\varphi $ 计算公式为

$ \varphi = \angle S $ (1)

式中,∠(·)是从一个复数中获得角度在(-π, π] rad的函数,因此缠绕相位与真实相位之间的关系可以写成

$ \phi = \varphi + 2k\pi $ (2)

$k$ 是整数。相位解缠绕就是为相位图中的每一个像素都找到一个正确的相位补偿 $k$ 来恢复真实相位。

1.2 局部多项式曲面拟合

在许多相位解缠绕方法中,一个关键步骤就是利用一个平滑函数模型来逼近真实相位[15-17]。本文利用局部多项式曲面模型来拟合局部窗内的真实相位,即

$ {\phi _{(x,y)}} = \sum\limits_{m = 0}^M {\sum\limits_{n = 0}^N {{C_{m,n}}{x^m}{y^n} + {e_{(x,y)}}} } $ (3)

式中, ${{C_{m, n}}}$ 表示多项式拟合系数, $M$ $N$ 为多项式拟合函数沿着 $x$ $y$ 方向的阶, ${{e_{(x, y)}}}$ 表示多项式模型的逼近错误。

因此,利用多项式曲面模型拟合真实相位,可以有效地把相位解缠绕问题转化为一个参数估计问题。利用已解缠绕像素来拟合等式(3)的函数模型从而确定拟合系数 ${{C_{m, n}}}$ 。而在局部多项式拟合中,假定在一个以生长像素 $({x_0}, {y_0})$ 为中心的局部窗内有 $P$ 个已解缠绕像素,因此等式(3)以矩阵形式可以表示为

$ \mathit{\pmb{\Phi }} = \mathit{\boldsymbol{X\hat c}} + \mathit{\boldsymbol{E}} $ (4)

式中, $\mathit{\pmb{\Phi }}$ 为列向量,包含了 $P$ 个已解缠绕像素的真实相位, $\mathit{\boldsymbol{X}}$ $P \times \left({M + 1} \right)\left({N + 1} \right)$ 的矩阵,每一行对应的是一个已解缠绕像素的多项式基, ${\mathit{\boldsymbol{\hat c}}}$ 为列向量,包含多项式拟合系数, $\mathit{\boldsymbol{E}}$ 为包含拟合错误的列向量。通过等式(3)拟合 $P$ 个已解缠绕像素的真实相位值,可以计算得到拟合系数 ${\mathit{\boldsymbol{\hat c}}}$ 。这时,生长像素的真实相位可以计算为

$ {\phi _{({x_0},{y_0})}} = R\left( {\frac{{{X_{({x_0},{y_0})}}\hat c - {\varphi _{({x_0},{y_0})}}}}{{2\pi }}} \right) \times 2\pi + {\varphi _{({x_0},{y_0})}} $ (5)

式中, ${{X_{({x_0}, {y_0})}}}$ 是生长像素的多项式基, $R$ (·)是四舍五入取整的函数。

1.3 提出方法

整个算法流程图如图 1所示。有许多不同的方法可以用来分块。采用PRELUDE中的分块方法,该方法把整个相位图分成多个连通块,而每一个块内任一像素相位值都在给定相位间隔内。也就是说,把缠绕相位区间(-π, π]偶等分(如六等分),然后对一个相位间隔,创建一个掩模。如果图像中某个像素的缠绕相位值在该相位间隔内,则掩模内相同位置的像素值设置为1。确定掩模内块个数,并对每一个块做标记。把相素个数小于给定阈值(如50)的块划分为残余像素,然后对块重新标记。对其他相位间隔做重复操作。相位解缠绕过程是基于相邻像素的真实相位差不大于π的假设。因此相位分区方法确保了单个块内相邻像素的相位差不大于相应的相位间隔尺寸(如π/3),若相邻像素真实相位变换不大于π,相位分区方法获得的块内无缠绕发生。另外,相位间隔划分的越小,任一块内也就越不可能包含缠绕。

图 1 提出算法流程图
Fig. 1 Flowchart of the proposed algorithm

块与块之间相位解缠绕采用了局域增长的策略,并把每一个块当成一个基本单元,块内所有像素的相位都有相同的整数补偿。感兴趣区域内像素最多的块被选为起始块,并假设块内相位已解缠绕。将要被解缠绕的块被定义为生长块。距离已解缠绕区域最近的块被选为生长块。既要考虑已解缠绕区域与未解缠绕块的质心距离,又要考虑二者最近像素之间的欧氏距离。

生长块的解缠绕等价于为生长块内所有像素找到一个最优的整数补偿 $k$ 。选取位于生长块内并距离已解缠绕区域最近的 ${Q_{\rm{g}}}$ 个像素和位于已解缠绕区域内并距离生长块最近的 ${Q_{\rm{u}}}$ 个像素参加曲面拟合。假定生长块的最优整数补偿在区间 $\left[ { - K, K} \right]$ (如 $K$ 等于5)内,使得多项式曲面拟合错误最小的整数被选取为生长块的最优整数补偿。

当所有的块已完成解缠绕并合并,使用基于局域增长的多项式曲面拟合方法逐像素对感兴趣区域内的残余像素解缠绕。具体解缠绕过程如下:

1) 选取最靠近已解缠绕区域的残余像素;

2) 对步骤1)获得残余像素,根据其幅值大小按降序排列;

3) 根据步骤2)确定的先后顺序,根据式(4)(5),使用以生长点为中心的局部窗内已解缠绕像素的相位来对生长点相位解缠绕;

4) 重复步骤1)2)3)直至所有的残余像素完成解缠绕。

2 实验

2.1 实验数据

利用标准差为50个像素的高斯函数生成一个尺寸为256 × 256的相位图。仿真数据的幅值图像为一个12等分的扇形圆盘,幅值从120以幅度为10的节奏递减到10[17-18]。添加均值为0,标准差为20的随机噪声到生成复数数据的实部和虚部上。因此,拟合数据的信噪比从6以幅度0.5的节奏递减到0.5。为了测试提出方法在不同情况下的表现,分别生成了3个仿真数据集。每一个数据集有不同的特性,具体如下:

1) 为了测试本文方法在不同信噪比下的表现,生成了相位高度为40 rad的仿真数据,如图 1所示。

2) 为了测试本文方法在不同相位变换水平下的表现,生成了相位高度分别为100 rad、150 rad和200 rad的仿真数据,如图 2所示。

图 2 提出方法与PRELUDE在不同信噪比下的解缠绕结果
Fig. 2 Phase-unwrapping results of PRELUDE and the proposed method under different SNRs ((a) reference phase; (b) wrapped phase; (c) magnitude; (d) and (e) the unwrapped results by PRELUDE and the proposed method; (f) and (j) the error maps by PRELUDE and the proposed method)

3) 为了测试提出方法在存在不连通区域情况下的表现,生成了相位高度为60 rad的仿真数据,3个预先设定的区域设置为0,如图 3所示。

图 3 使用不同相位变换水平的数据来比较本文方法与PRELUDE(均值±标准差/%)
Fig. 3 Comparison of the proposed method with PRELUDE using the dataset with different phase variations levels ((a) phase height equal to 100 rad; (b) phase height equal to 150 rad; (c) phase height equal to 200 rad) (means ± standard deviations/%)

为了评价仿真实验结果,拟合相位与噪声引起相位变换的和作为参考相位。从解缠绕相位中减去参考相位,如果某个像素相位差的绝对值大于π/10 rad,则认为该像素解缠绕失败。感兴趣区域内,错误率定义为解缠绕错误像素个数占总像素个数的比率,并以此来定量评价提出方法。每一类仿真实验重复了20次,并计算相应错误率的均值和方差。

2.2 真实实验数据

使用一台0.35 T永磁磁共振机器(XGY-OPER,宁波鑫高益公司,中国)获得了5个健康成年人的膝关节和踝关节数据。采用三点Dixon序列收集可用于水与脂肪分离的数据。获取数据参数为:脉冲重复时间/回波时间=580/28 ms,矩阵大小为256×256,视野窗=240×240 mm2,每组膝关节数据和踝关节数据分别包含了10层和5层。

另外使用一台3.0 T磁共振机器(Achieva,飞利浦,荷兰)获得了一个成年人的膝关节数据。获取双回波梯度成像数据的参数为:脉冲重复时间/回波时间1/回波时间2=152/2.30/3.43 ms,矩阵大小为192×192, 视野窗=180×180 mm2,包含了25层。所有的志愿者在实验之前告知了相关机构审查委员会的政策。

为了评价本文方法在真实数据上的表现,应用真实数据的相位解缠绕结果到两点或三点Dixon技术中生成水和脂肪图像。如果解缠绕结果出错,则生成的水和脂肪图像会发生互换,即水图中出现了脂肪,反之亦然。评价了每一层的水与脂肪分离结果,并按以下标准分类[17, 19]:1)很多互换;2)较少互换;3)总的互换层数;4)错误率。错误率定义为发生水与脂肪互换层的个数与总的层个数的百分比。

实验平台为HP EliteDesk 800(Intel i7 3.6 GHz RAM: 4 GB),程序采用Matlab2014实现。在局部多项式曲面拟合中,涉及多项式拟合函数阶 $M$ $N$ 和局部拟合窗的大小。多项式阶 $M$ $N$ 同时设置为2,因为2阶多项式函数能够准确地描述局部窗内的真实相位,而增加阶数,虽然稍微提高计算精度但增加了计算复杂度。在块与块之间合并解缠绕时,选取距离已解缠绕区域最近的位于生长块内50个像素和距离生长块最近的位于已解缠绕区域内50个像素参加拟合。选取50个像素参与拟合是由经验得到的,并等于最小块的尺寸。基于观察仿真实验,残余像素解缠绕拟合窗的大小设置为15,并应用到接下来的所有实验中去。

3 结果与讨论

3.1 仿真实验

图 2显示了提出方法与PRELUDE在不同信噪比下的解缠绕结果。在低信噪比区域,PRELUDE生成的解缠绕结果包含明显缠绕,而提出方法的解缠绕错误则明显降低。表 1提供了两种方法重复20次的错误率均值和标准差。当信噪比大于2.5时,两种方法都获得了准确的解缠绕结果,两种方法结果的错误率均值和标准差都低于0.01%。当信噪比低于2.5时,提出方法生成结果的错误率均值和标准差都比PRELUDE的相应结果要低。

表 1 重复20次后两种方法的错误率均值和标准差
Table 1 Means and standard deviations of errors ration using the two methods over 20 repetitions under different SNRs

下载CSV
/%
信噪比 PRELUDE 提出方法
1 6 0.00±0.00 0.00±0.00
2 5.5 0.00±0.00 0.00±0.00
3 5 0.00±0.00 0.00±0.00
4 4.5 0.00±0.00 0.00±0.00
5 4 0.00±0.00 0.00±0.00
6 3.5 0.00±0.00 0.00±0.00
7 3 0.00±0.00 0.00±0.00
8 2.5 0.01±0.01 0.00±0.00
9 2 0.06±0.03 0.01±0.01
10 1.5 0.35±0.12 0.09±0.04
11 1 1.81±0.30 0.19±0.10
12 0.5 6.43±1.32 0.51±0.24

图 3比较了本文方法与PRELUDE在不同相位变换水平下的表现。无论相位高度是100 rad还是200 rad,PRELUDE都未能获得准确的解缠绕结果,其解缠绕错误率均值随着相位高度的增加而快速增加。而本文方法均获得了满意的解缠绕结果,解缠绕错误率均值虽然也随着相位高度的增加而增加,但都低于0.36%。

图 4显示了本文方法与PRELUDE在出现不连通区域时的表现。PRELUDE生成的结果都包含严重的解缠错误,其解缠绕错误率的均值和标准差为12.79±0.67%。而本文方法能够生成准确的解缠绕结果,解缠绕错误率的均值和标准差为0.10±0.05%。

图 4 使用存在不连通区域的数据来比较本文方法与PRELUDE(均值±标准差/%)
Fig. 4 Comparison of the proposed method with PRELUDE in the presence of disconnected regions ((a) reference phase; (b) wrapped phase; (c) magnitude; (d) and (e) the unwrapped results by PRELUDE and the proposed method; (f) and (j) the error maps by PRELUDE and the proposed method) (means ± standard deviations/%)

仿真实验结果表明,PRELUDE对噪声、快速相位变换和不连通区域敏感,而本文方法则能获得满意的解缠绕结果,这是因为本文方法利用周围已解缠绕的相位信息来对低信噪比区域解缠绕,有效地降低了增长路径中误差的传递和累积。同时,本文方法利用局部多项式曲面拟合方法代替PRELUDE中的代价函数法来合并块和对残余像素相位解缠绕,由于多项式曲面模型对噪声和不连通区域不敏感,因此能够获得准确的解缠绕结果,甚至当相邻像素相位差大于π时。

3.2 真实磁共振成像实验

图 5图 6比较了本文方法与PRELUDE在0.35 T(Tesla)真实水与脂肪Dixon磁共振数据上的表现。图 5展示了一个典型的膝关节层面相位解缠绕和水与脂肪分离结果。图 6展示了一个典型的踝关节层面相位解缠绕和水与脂肪分离结果。如箭头所指,PRELUDE生成结果包含严重的缠绕残余。该现象导致了水与脂肪互换发生。而本文方法生成的结果中并无明显的解缠绕错误和水与脂肪互换现象发生。

图 5 一个0.35 T膝关节层面的相位解缠绕和水与脂肪分离结果
Fig. 5 Phase-unwrapping and water-fat separation results of a slice in the 0.35 T knee dataset ((a) wrapped phase; (b) out-phase magnitude; (c) in-phase magnitude; (d) the unwrapped results by PRELUDE and the proposed method; (e) water maps by PRELUDE and the proposed method; (f) fat maps by PRELUDE and the proposed method)
图 6 一个0.35 T踝关节层面的相位解缠绕和水与脂肪分离结果
Fig. 6 Phase-unwrapping and water-fat separation results of a slice in the 0.35 T ankle dataset ((a) wrapped phase; (b) out-phase magnitude; (c) in-phase magnitude; (d) the unwrapped results by PRELUDE and the proposed method; (e) water maps by PRELUDE and the proposed method; (f) fat maps by PRELUDE and the proposed method)

图 7比较了本文方法与PRELUDE在3.0 T真实水与脂肪Dixon磁共振数据上的表现。如白色箭头所指,PRELUDE生成结果存在缠绕残余,并导致了水与脂肪互换发生。而本文方法生成结果中并无明显的解缠绕错误和水与脂肪互换。

图 7 一个3.0 T膝关节层面的相位解缠绕和水与脂肪分离结果
Fig. 7 Phase-unwrapping and water-fat separation results of a slice in the 3.0 T knee dataset ((a) wrapped phase; (b) out-phase magnitude; (c) in-phase magnitude; (d) the unwrapped results by PRELUDE and the proposed method; (e) water maps by PRELUDE and the proposed method; (f) fat maps by PRELUDE and the proposed method)

表 2提供了所有0.35 T和3.0 T水与脂肪分离的统计结果(共100层:75层膝关节和25层踝关节)。总的来说,在100个结果中,PRELUDE生成结果中包含了9次很多互换和33次较少互换。而本文方法无许多互换发生,只发生了6次较少互换。PRELUDE总错误率为42.00%,而本文方法总错误率仅为6.00%。

表 2 6位志愿者水与脂肪分离数据集中互换发生层的统计结果
Table 2 Statistical results of slices with swaps in the water-fat separation datasets of 6 volunteers

下载CSV
数据集 方法 很多互换 较少互换 错误层数 错误率/%
膝关节
(0.35 T)
PRELUDE 2 7 9 18.00
本文 0 2 2 4.00
踝关节
(0.35 T)
PRELUDE 1 11 12 48.00
本文 0 3 3 12.00
膝关节
(3.0 T)
PRELUDE 6 15 21 84.00
本文 0 1 1 4.00
总数 PRELUDE 9 33 42 42.00
本文 0 6 6 6.00

磁共振真实数据表明,本文方法明显优于原始的PRELUDE方法。PRELUDE失败的原因可能是由于在采集的膝关节和踝关节数据中,受多种组织和空气的影响下,相位发生了急剧的变化,甚至有相邻像素真实相位变换大于π的情况发生,而使用代价函数合并的前提假设却是相邻相位变换不大于π。违背前提条件将使得解缠绕结果发生错误,并使最终的水与脂肪分离结果出错。而本文方法使用阈值法把较小的块归类于残余像素,这是因为信噪比越低区域,相邻像素相位波动较大,相位分区划分的块也通常越小。另外,在完成块与块之间合并解缠绕之后再利用局域增长的局部多项式拟合方法对残余像素相位解缠绕,降低了低信噪比区域像素点提前出现解缠绕路径中可能性,降低了解缠绕错误的传递和累积。

4 结论

针对存在严重噪声、快速相位或不连通块的磁共振相位图像,本文提出了一种基于相位分区和局部多项式曲面拟合的相位解缠绕新方法。本文方法首先利用相位分区方法把相位图像分成连通块,块内相邻像素相位变换不大于给定的相位间隔,因此块内不含相位缠绕。由于低信噪比区域相位变化较紊乱,分的块较小,为了降低噪声的影响,利用阈值法把连通块分为块和残余像素两类。最后再利用局部曲面拟合方法依次进行块与块之间和残余像素解缠绕。这是由于局部多项式曲面模型对噪声不太敏感,在局域增长过程中利用周围相位信息来估计待解缠绕像素真实相位值,提高了解缠绕结果的鲁棒性。另外当相邻真实相位变换大于π,利用多项式曲面模型能够拟合出真实的相位变化。

本文方法可有效地克服PRELUDE算法中对噪声、快速相位变换或存在不连通区域敏感的缺点,提高了解缠绕结果的准确性。当获取相位图像存在严重噪声、快速相位变换或不连通区域时,仿真实验结果表明本文方法仍能获得准确的解缠绕相位。水与脂肪Dixon磁共振成像数据实验表明本文方法可以获得准确的解缠绕相位,并应用到Dixon技术中去实现水与脂肪的可靠分离。因此,本文方法有益于那些需要相位解缠绕的医学磁共振成像应用,如水脂分离,磁敏感权重成像和人脑相位成像等。

与PRELUDE相比,本文方法中有3个参数需要经验选取,如块与块之间解缠绕过程中曲面拟合点个数、残余像素解缠绕半径和生长点的选取等,这将给本文方法的适用性带来了困难。如何可靠鲁棒的解决参数选取问题以及把提出算法扩展到3维将是下一步工作的研究方向。

参考文献

  • [1] Munger P, Crelier G R, Peters T M, et al. An inverse problem approach to the correction of distortion in EPI images[J]. IEEE Transactions on Medical Imaging, 2000, 19(7): 681–689. [DOI:10.1109/42.875186]
  • [2] Glover G H, Schneider E. Three-point dixon technique for true water/fat decomposition with ${B_0}$ inhomogeneity correction[J]. Magnetic Resonance in Medicine, 1991, 18(2): 371–383. [DOI:10.1002/mrm.1910180211]
  • [3] Coombs B D, Szumowski J, Coshow W. Two-point Dixon technique for water-fat signal decomposition with ${B_0}$ inhomogeneity correction[J]. Magnetic Resonance in Medicine, 1997, 38(6): 884–889. [DOI:10.1002/mrm.1910380606]
  • [4] Nayler G L, Firmin D N, Longmore D B. Blood flow imaging by cine magnetic resonance[J]. Journal of Computer Assisted Tomography, 1986, 10(5): 715–722. [DOI:10.1097/00004728-198609000-00001]
  • [5] Ghiglia D C, Pritt M D. Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software[M]. New York: Wiley, 1998.
  • [6] Ma J F. Dixon techniques for water and fat imaging[J]. Journal of Magnetic Resonance Imaging, 2008, 28(3): 543–558. [DOI:10.1002/jmri.21492]
  • [7] Zhu Y J, Liu L R, Luan Z, et al. A reliable phase unwrapping algorithm based on the local fitting plane and quality map[J]. Journal of Optics A:Pure and Applied Optics, 2006, 8(6): 518–523. [DOI:10.1088/1464-4258/8/6/005]
  • [8] Ghiglia D C, Romero L A. Minimum ${L_p}$ -norm two-dimensional phase unwrapping[J]. Journal of the Optical Society of America A, 1996, 13(10): 1999–2013. [DOI:10.1364/JOSAA.13.001999]
  • [9] Xu W, Cumming I. A region-growing algorithm for InSAR phase unwrapping[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(1): 124–134. [DOI:10.1109/36.739143]
  • [10] Witoszynskyj S, Rauscher A, Reichenbach J R, et al. Phase unwrapping of MR images using ΦUN-a fast and robust region growing algorithm[J]. Medical Image Analysis, 2009, 13(2): 257–268. [DOI:10.1016/j.media.2008.10.004]
  • [11] Itoh K. Analysis of the phase unwrapping algorithm[J]. Applied Optics, 1982, 21(14): 2470. [DOI:10.1364/AO.21.002470]
  • [12] Goldstein R M, Zebker H A, Werner C L. Satellite radar interferometry:two-dimensional phase unwrapping[J]. Radio Science, 1988, 23(4): 713–720. [DOI:10.1029/RS023i004p00713]
  • [13] Strand J, Taxt T, Jain A K. Two-dimensional phase unwrapping using a block least-squares method[J]. IEEE Transactions on Image Processing, 1999, 8(3): 375–386. [DOI:10.1109/83.748892]
  • [14] Jenkinson M. Fast, automated, $N$ -dimensional phase-unwrapping algorithm[J]. Magnetic Resonance in Medicine, 2003, 49(1): 193–197. [DOI:10.1002/mrm.10354]
  • [15] Liang Z P. A model-based method for phase unwrapping[J]. IEEE Transactions on Medical Imaging, 1996, 15(6): 893–897. [DOI:10.1109/42.544507]
  • [16] Liu W T, Tang X, Ma Y J, et al. 3D phase unwrapping using global expected phase as a reference:application to MRI global shimming[J]. Magnetic Resonance in Medicine, 2013, 70(1): 160–168. [DOI:10.1002/mrm.24448]
  • [17] Cheng J Y, Mei Y J, Liu B S, et al. A novel phase-unwrapping method based on pixel clustering and local surface fitting with application to Dixon water-fat MRI[J]. Magnetic Resonance in Medicine, 2018, 79(1): 515–528. [DOI:10.1002/mrm.26647]
  • [18] Maier F, Fuentes D, Weinberg J S, et al. Robust phase unwrapping for MR temperature imaging using a magnitude-sorted list, multi-clustering algorithm[J]. Magnetic Resonance in Medicine, 2015, 73(4): 1662–1668. [DOI:10.1002/mrm.25279]
  • [19] Cui C, Wu X D, Newell J D, et al. Fat water decomposition using globally optimal surface estimation (GOOSE) algorithm[J]. Magnetic Resonance in Medicine, 2015, 73(3): 1289–1299. [DOI:10.1002/mrm.25193]