Abstract
A new method based on discrete particle swarm optimization (dPSO) algorithm is proposed to solve the branch-cut phase unwrapping problem of MRI data. In this method, the optimal order of matching the positive residues with the negative residues is first identified by the dPSO algorithm, then the branch cuts are placed to join each pair of the opposite polarity residues, and in the last step phases are unwrapped by flood-fill algorithm. The performance of the proposed algorithm was tested on both simulated phase image and MRI wrapped phase data sets. The results demonstrated that, compared with conventionally used branch-cut phase unwrapping algorithms, the dPSO algorithm is rather robust and effective.
1. Introduction
In magnetic resonance imaging (MRI), the complex signal contains both the magnitude and phase parts. Usually the magnitude of the MRI signal has been mainly considered. However, the phase of MRI signal offers very important information on the velocity of the moving spins, and can also be used to deduce useful information about the main field inhomogeneity and the magnetic susceptibility variations [1]. In MRI, the phase information is usually obtained from a complex MRI dataset through some mathematical operations, and the value always lies in the principal interval of , consequently producing a wrapped phase . This relationship can be described by , where is an integer and defines a wrapping operator that forces all values of its argument into the range by adding or subtracting an integral multiple of radians from its argument. Phase unwrapping is the process of estimating the true phase from the wrapped phase . As an important tool, it can not only be used for the three-point Dixon water and fat separation, but also be applied to increase the dynamic range of phase contrast MR velocity measurements [2]. If the true phase gradients (i.e., the differences of ) between contiguous pixels are less than π radians in magnitude in the entire space, the true phase can be unwrapped in a straightforward manner by just integrating the wrapped phase gradients [3]. However, the presence of the noise, undersampling, and/or object discontinuities often makes this condition unavailable. Therefore, the problem of phase unwrapping becomes complex in practice and difficult to solve, although significant amount of research effort has been devoted to date. In the literature, there are quite a few existing phase unwrapping algorithms [4], which can be grouped into two categories: path-following and minimum-norm methods [5].
The branch-cut phase unwrapping method is one kind of the path-following methods. Unlike the minimum-norm methods, the branch-cut phase unwrapping technique offers correct phase unwrapping with no solution approximations [4]. In the operation, it locates residues, joins the residues by branch cuts, and then unwraps all the pixels avoiding those branch cuts. In the algorithm, residues are defined as local inconsistencies, which mark the starting and end of discontinuities. Corresponding to formula (1), that is, the value is 1 or −1 in a 2 × 2 closed-loop of the wrapped phase gradients [4], as shown in Figure 1:
These wrapped phase gradients are computed counterclockwise:
In formula (1), when , we label it a residue with a positive polarity. And when , the residue is labeled with a negative polarity. Otherwise it indicates there is no residue when . In this algorithm, the total length of branch cuts must been minimized, resulting a decrease of the amount of good pixels used as the branch cuts. This provides more unwrapping paths in dense residue areas, leading to a smoother result [5].
To achieve a minimum total length of branch cuts, various techniques have been developed for the branch-cut phase unwrapping method. These techniques include, for example, the Goldstein’s algorithm [5, 6], the nearest-neighbor algorithm [7], the minimum-cost matching (MCM) algorithm [8], and the hybrid genetic algorithm (HGA) [4].
The nearest-neighbor algorithm is very efficient, but it utilizes local heuristics, thus causing some long branch cuts embedded in the phase image. Therefore, the distribution of the branch cuts does not achieve the optimum, and the resultant phase image is lack of the smoothness. The MCM is a graph theory-based algorithm, which uses the Hungarian algorithm to find the minimum total length of branch cuts. Although powerful, it is computationally expensive. The HGA employs a combination of global and local search methods, and its solution is usually good. However, the complexity of this algorithm tends to be a problem with the increase of residues.
Compared with these three methods, which place the branch cuts to connect pairs of residues of opposite polarity (called dipoles), the Goldstein’s method joins the residues in clumps instead of pairs [5]. The Goldstein’s algorithm is very efficient, but often forms some isolated patches.
In this paper, we propose a new discrete particle-swarm-optimization- (dPSO-) based branch-cut algorithm for phase unwrapping. The new dPSO algorithm is used to find the best way in which the negative polarity residues match with the positive ones, so that the overall length of the branch cuts is minimized. The performance of the new dPSO algorithm is compared with the Goldstein’s and MCM algorithms.
2. Phase Unwrapping Using the Proposed Algorithm
2.1. Overview of the Basic PSO Algorithm
As an artificial intelligent algorithm, the particle swarm optimization (PSO) [9, 10], is easy for implementation (only few parameters to be adjusted) and converges fast.
In the PSO algorithm, the swarm consists of several particles, and each contains elements. Then each particle is viewed as a point in an -dimensional space. The th particle of swam is represented as . All the particles share their information and move to find the global optima. represents that the local best position that the th particle has reached, and is the global best position. The velocity of the th particle is . Each particle of the swarm updates its velocity and position using the following formulas: where denotes the iteration number, and are learning factors (nonnegative constants), controlling (or regulating) the influence of and , the function and generate a random number ([0~1]), and is the inertia weight factor.
A problem-specific fitness function (symbolized by ) is employed to measure the performance of each particle. Thereby, for a minimization problem, and can be found in current iteration as follows:
To date the PSO technique has been well developed for the continuous problem, but not in discrete domain.
2.2. The dPSO Algorithm for Phase Unwrapping
2.2.1. Particle and Swarm Initialization
Any problem adopting the PSO algorithm has to be interpreted into PSO particle form. For phase unwrapping, every particle should be composed of some elements corresponding to the indexes of residues. If we calculate all the residues in the wrapped phase image at one stroke, sometimes the size of the particle may be too large and the swarm size required may be enlarged accordingly. This may result in a poor solution and/or extending the convergence time. To avoid this, the whole image is divided into sub regions and therefore the residues are set into different small groups.
The process can be described as follows.(1)The image is partitioned on the basis of its phase derivative variance map. The phase derivative variance is defined as follows [5]: where for each sum the indexes cover over the window centered at the pixel . The terms and are the wrapped phase gradients in the windows, and and are the averages of these wrapped phase gradients. In this paper . (2)Based on an appropriate threshold, the phase derivative variance map can be converted into a binary one, where the low phase derivative variance values turns out to be 0 and the high ones becomes to be 1. In this way the image is divided into separate areas. To obtain a proper threshold, a classical approach called Otsu’s method [11] has been adopted in this study. (3)It is observed that the majority of residues cluster in the patches of value 1. Thus the residues are grouped by taking the above steps.
In each region the indexes of the residues are inserted in two arrays regardless of the order. One is a positive polarity residue array, and the other is a negative polarity residue array. Provided there are positive polarity residues and negative polarity residues in one region, respectively, the positive residue array and the negative one are accordingly denoted as and . The former will be fixed throughout all generations and serves as a reference. And the later acts as a particle of the initial swarm. The rest particles of the swarm are generated via arranging the order of the elements in in a random manner.
2.2.2. Fitness Estimate
From the abovementioned, it can be easily seen that the dPSO algorithm is used to find out the best matched order of the elements in the particle with the reference array.
In dPSO, the quality of the current solution is judged by the fitness function. Since the branch-cut phase unwrapping must minimize the total length of branch cuts, the corresponding fitness function is obviously for calculating the total length of branch cuts in the wrapped phase image. Here we employ Euclidian distance to assess the total cut length: where and denote the residue’s -coordinate and -coordinate.
2.2.3. Velocity Update
Owing to the attribute of dPSO algorithm for branch-cut phase unwrapping, the iterative velocity in formula (3) should be a set of permutation operators rather than a usual vector. Various permutation operators have been introduced for discrete particle swarm optimization. Here we choose the adjustment operator [12], not the swap operator [13, 14], as the permutation operator. This is because compared with the swap operator, the adjustment operator avoids returning to the previous position [12].
In addition, during all the iterations, w is set to be linearly varied ([0.9~0.4]) [10]. This setting considers the global searching capacity and convergence rate of the optimization: where is the maximal iteration times.
To make formula (3) suitable for the dPSO operation, some concepts are given as follows.
Definition 1. The adjustment operator is defined as deleting the element in the th position and popping it in the th position in the array.
For example, acting on the array gets a result .
Definition 2. One or more adjustment operators make up an adjustment sequence (AS). That is .
Acting an AS on an array means that every adjustment operator of the sequence acts on the array in turn. Therefore, in the AS, it is critical to properly arrange the order of adjustment operators.
Definition 3. The plus sign “+” between ASs has its new meaning. It is defined as forming a new longer AS by putting the latter behind the former.
For example, , , . This operation does not satisfy the commutative law.
Definition 4. The minus sign “−” between two arrays means constructing an AS which can act on the array after minus sign to obtain the array before.
Supposing there are two arrays, according to the array before minus sign from left to right, this AS can be obtained by adjusting the order of the array after. For example, in and . , so the first adjustment operator of AS is , and the first new array . Then , so the second adjustment operator is , and the second new array is . W(4) = , so the third adjustment operator is , and we finally obtain the same array as . Thus the calculation ends and .
Definition 5. The multiplication sign “” between two real numbers is simply the operation of multiplication.
Definition 6. The sign “” between a real number and an AS means reserving a certain amount of the adjustment operators in the AS in turn when the real number is in the range (0, 1).
Given that the number of adjustment operators in the AS denotes as , the number of retaining adjustment operators is , which rounds the product of and to the nearest integer.
2.2.4. Position Update
Similarly, the formula (4) has to be redeclared to fit the requirement of the dPSO.
Definition 7. The plus sign “+” between an array and an AS means acting the adjustment operators of the AS on the array in order. For example, given and .
It is obvious that the operations defined in Definitions 4 and 7 are inverse to each other.
2.2.5. The Procedure Description of dPSO Algorithm for Phase Unwrapping
The process of using the dPSO algorithm for phase unwrapping is summarized as follows.(1)According to Section 2.2.1, the residues in the image are divided into several clusters. Set .(2)Do the following steps in the th group of residues.(3)Set values to the parameters of dPSO, which include the learning factors (), the maximal iteration times (), and the number of particles in swarm, and also the termination criteria.(4)Initialize the swarm according to Section 2.2.1. Each particle has its random velocity, that is, AS. Set .(5)Evaluate the fitness of every particle according to formula (8) and find the current by formula (5), (6), respectively. Calculate the current inertia weight factor according to formula (9).(6)Set . Use formula (3) to get the new velocity . Then calculate the new position according to formula (4).(7)Repeat (5)-(6) until or meeting the termination criteria.(8) is the best indexes order of the negative polarity residues matched with the positive ones in this group.(9)Set . Repeat (2)–(8) until equals the number of the residues groups plus 1.
2.3. Branch Cuts and Unwrapping
Once the best match in every group has been found by dPSO, each pair of two matching opposite polarity residues are connected by the branch cuts. It is worth mentioning that, owing to that the number of opposite polarity residues is not always equal to each other, there are usually one or more residues left in each group. Then the nearest-neighbor algorithm [7] is employed to place branch cuts to balance these remaining residues. So far all the residues have been balanced by branch cuts.
Finally, the phase data can be unwrapped by flood-fill algorithm [15, 16], without crossing the branch cuts as follows:(1)Choose a start pixel, whose phase value is stored as an unwrapped phase value in the solution matrix. The four neighboring pixels are unwrapped next and their unwrapped phase values are placed in the solution matrix. These four pixels are inserted in the unwrapped list. (2)Pick (and then eliminate) a pixel from the unwrapped list. Unwrap the phase values of its four neighboring pixels, avoiding pixels that have been unwrapped. Insert these pixels in the unwrapped list and put their unwrapped phase values in the solution matrix.(3)Repeat (2) until the unwrapped list becomes empty.
In fact, it does not always mean that all the pixels have been unwrapped when the unwrapped list becomes empty. Because sometimes there are some pixels in the image encircled by the branch cuts, they cannot be unwrapped if not crossing the branch cuts.
2.4. Weighted Measure
Weighted measure, the most general/practical error measure to consider [5], is used to evaluate the quality of an unwrapped solution: where and are the number of rows and columns, and are user-defined weights, and the norm measures a count of the number of pixels at which the gradients of the unwrapped solution mismatch the wrapped phase gradients. In this paper the weights adopted are derived from the quality map mentioned above, not just omitted (i.e., equal to 1).
3. Results and Discussion
We have tested the performance of the proposed algorithm on both simulated and MRI phase data on a PC (Intel 2 Quad CPU 2.39 GHz, MATLAB). We set both learning factors () used in dPSO to be 2. The results were compared with the well-known Goldstein’s and the MCM algorithms.
3.1. Simulation Results
The proposed algorithm was implemented on a simulated wrapped phase image with salt and pepper noise (the signal-to-noise ratio is 7.58 dB), which has 2460 residues. Figure 2(a) shows the simulated wrapped phase image. And its residue distribution is shown in Figure 2(b), where the positive and negative residues are marked as white and black pixels, respectively.
(a)
(b)
The resultant unwrapped phase image in Figure 3(a) was achieved by dPSO using a swarm of 300 particles and . Figures 3(b) and 3(c) depict the corresponding unwrapped phase images obtained by Goldstein’s and MCM algorithms, respectively. We rewrapped these resultant solutions and subtracted the original wrapped phase data from them. The corresponding difference maps, which are plotted in 3D visualization, are shown in Figures 3(d)–3(f). As mentioned in Section 1, the difference between the wrapped phase and the true phase is an integer multiples of . As a result, these curve surfaces of difference maps give a visual representation of the deviations between the true phases and the results, which can illustrate the accuracy of the unwrapped solutions directly. The average value of the difference map is called average difference.
(a)
(b)
(c)
(d)
(e)
(f)
The performance of dPSO with respect to the other two algorithms can clearly be seen in Table 1. In terms of weighted measure, average difference, and total cuts length, the dPSO algorithm is better than the Goldstein’s, but not as good as the MCM. On the other hand, the execution time of dPSO is much less than that of the MCM, but not comparable to that of the Goldstein’s.
3.2. Results of MRI Data
The proposed algorithm was also executed on a displacement encoded MRI heart phase data set [17] with 396 residues. The wrapped phase image and its corresponding residue distribution are shown in Figure 4.
(a)
(b)
Figures 5(a)–5(c) depict the branch-cut distribution achieved by the three algorithms, where the black pixels mark the branch cuts. The dPSO result in Figure 5(d) was obtained by using a swarm of 300 particles. And the results of the other algorithms are shown in Figures 5(e) and 5(f), respectively. In Figure 5(e) several patches are isolated, two large ones in the upper part, two small ones in the middle, and a very large one in the lower right part. Compared with Figures 4(b) and 5(b), it is easy to observe that these patches are completely isolated by branch cuts, which would lead to an incorrect unwrapping. In addition, the isolated areas tend to arise in the regions with dense residues, because in such regions the branch cuts are often close to each other and have more possibility to encircle some pixels. However, in Figures 5(d) and 5(f) the isolated patches are much smaller and less. The dipole branch-cut methods appear to be less likely to isolate regions in the phase image by branch cuts, since the branch cuts balance the residues in pairs not in clumps. The difference maps of the three approaches are then generated like in Section 3.1, shown in Figures 5(g)–5(i). Intuitively, both the dPSO and MCM produce more desirable results than the Goldstein’s.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
As shown in Table 2, the Goldstein’s algorithm is extremely fast. Neither dPSO nor MCM is comparable to it. But the proposed approach has the smallest weighted measure and average difference. In addition, its total cuts length is the shortest.
Another example is the MRI head phase data. The wrapped phase image is shown in Figure 6(a). Figure 6(b) depicts its residues distribution. The unwrapped phase images achieved by the dPSO, Goldstein’s, and MCM algorithms are displayed in Figures 7(a)–7(c), respectively. Obviously, in Figure 7(b) the surrounding area was not unwrapped at all. According to the previous analysis, the Goldstein’s method isolates this area by branch cuts. In comparison to Figure 6(a), the surrounding part roughly is the back ground with dense noise. However, due to the property of dipole branch-cut phase unwrapping method the branch cuts could hardly enclose the surrounding area, which certainly caused unwrapping phases in this area. The same inferences can be also made according to the difference maps of three methods shown in Figures 7(d)–7(f).
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
The weighted measure, average difference and total cuts length are calculated over the whole image whether the pixel is inside the region of interest (ROI) or background. Thereby in these three respects, as shown in Table 3, the dPSO and MCM algorithms are much better than the Goldstein’s. Though the dPSO method does not get a better solution than that of the MCM, there is little difference between them. That is, dPSO is comparable to MCM. Furthermore the former converges nearly 68% faster than the latter.
Viewing dPSO’s performance on these three examples, we can find that the dPSO algorithm takes more time to achieve an optimal solution when plenty of residues uniformly scatter throughout large areas. This is because the pixels in each of these areas often have similar quality and then the residues in each area can hardly be separated into more than one group, which results in the increase of the particle size for every group.
4. Conclusions
We have presented a new branch-cut phase unwrapping method based on dPSO algorithm in this paper. Both simulated and real wrapped phase data were used to test the performance of the proposed algorithm. The results of dPSO were compared with the Goldstein’s and the MCM algorithms. It was found that the dPSO method is better than Goldstein’s algorithm in terms of weighted measure, average difference and total branch cuts length. Moreover, the dPSO is much faster than the MCM algorithm in getting a global optimum solution while it is comparable to the latter in terms of weighted measure, average difference and total branch cuts length. Generally speaking, it has been demonstrated to be robust, effective for the phase unwrapping application.
In addition, it is capable of dealing with large branch-cut problem with thousands of residues. The complexity of the dPSO algorithm increases when the number of residues in a group increases, as the length of the particle extends which requires a larger swarm size. Future research will make this algorithm to be more efficiently operated for the phase unwrapping study.
Acknowledgment
This work was supported by the 973 National Basic Research and Development Program of China (2010CB732502).