Sparse-Representation-Based Direct Minimum L p-Norm Algorithm for MRI Phase Unwrapping

A sparse-representation-based direct minimum L p-norm algorithm is proposed for a two-dimensional MRI phase unwrapping. First, the algorithm converts the weighted-L p-norm-minimization-based phase unwrapping problem into a linear system problem whose system (coefficient) matrix is a large, symmetric one. Then, the coefficient-matrix is represented in the sparse structure. Finally, standard direct solvers are employed to solve this linear system. Several wrapped phase datasets, including simulated and MR data, were used to evaluate this algorithm's performance. The results demonstrated that the proposed algorithm for unwrapping MRI phase data is reliable and robust.


Introduction
From the MRI complex data, the phase information can be extracted with a restricted interval (− , ]. That is, the phase value is wrapped. We call it a wrapped phase or the principal value . This relationship between the wrapped phase and its corresponding true phase can be described by = { } = + 2 , where is an integer and {⋅} is the wrapping operator, forcing the value of its argument inside the curly braces into the range (− , ] by adding or subtracting an integral multiple of 2 radians from its argument. However, what is needed is the true unknown phase , because this relates to certain properties of interest, such as the velocity of the moving spins, the main 0 field inhomogeneity, and the magnetic susceptibility variations. Given the wrapped phase, phase unwrapping is applied to restore the true phase, obtaining the unwrapped estimate . This technique is an important tool in many MRI applications, for example, threepoint Dixon water and fat separation [1], MR venography [2,3], motion tracking in a tagged cardiac MRI [4], and field mapping in EPI [5]. Should the wrapped phases have no inconsistencies, the process of the phase unwrapping will be merely integrating the phase gradients over a path that covers the whole domain of interest. This process is quite simple and path independent. No matter which path is followed, the results will be the same regardless of the constant offset. However, in practice, there are always inconsistencies owing to the presence of the noise, undersampling, and/or object discontinuities. Consequently, phase unwrapping becomes intractable and path dependent. The inconsistency is usually called "residue" and is detected by summing the wrapped phase gradients around each 2 × 2 closed loop in the two-dimensional (2D) array, as shown in Figure 1. Consider where Δ , and Δ , are defined to be the gradients of the wrapped phases: In the literature, numerous two-dimensional (2D) phase unwrapping approaches have been developed. They can be classified into four categories: path-following [6][7][8][9][10][11], minimum-norm [12][13][14][15][16][17][18], Bayesian/regularization [19,20], and 2 Computational and Mathematical Methods in Medicine i+1,j+1 i+1,j parametric modeling [21] methods. Our method belongs to the second category; thus, we briefly introduce the minimumnorm methods.
The minimum-norm methods estimate the unwrapped phases by minimizing the -norm of the differences between the gradients of the wrapped and the unwrapped phases. With = 2, this results in least-square algorithms. These employ the FFT/DCT transforms [18,22] or/and iterative techniques [22] to reach an approximation for the least-square solution. The exact least-square solution is obtained by applying network programming techniques in [15]. Nevertheless, the least-square minimization tends to smooth the discontinuities, unless these discontinuities are given in advance as binary weights. The minimum 1 -norm algorithms [13,23] have a better ability than the 2 -norm ones to preserve the discontinuities. As approaches zero, the minimumnorm method tends to obtain a more reliable result [24,25]. Thus, the 0 -norm minimization is widely accepted as the most desirable in practice. However, the 0 -norm minimization is a nondeterministic polynomial-time hard (NP-hard) problem [14] with only approximate, not exact, solutions being developed in [12,14]. The conventional minimum -norm method [12] converts the -norm minimization problem into a generalized matrix equation with a flexible option of . It yields more accurate solutions than other methods, except for the cases including some residues closer to the periphery than to the other residues with opposite polarity. Additionally, because it is implemented in a dualiterative structure (an iteration structure embedded in an outer one), it can be computationally intense.
In this work, a sparse-representation-based direct minimum -norm (SDM ) algorithm is proposed. It introduces the user-defined weights into the generalized matrix equation of the conventional minimum -norm method to improve the performance of phase unwrapping, because, in this way, the discontinuities in the unwrapped phase surface can be confined to the low-quality or zero-weight regions [22]. On the other hand, the sparse representations of the matrices in the modified matrix equation and the direct solvers are exploited in this algorithm to significantly reduce the computational time of the conventional minimum -norm method for MRI phase unwrapping.

Mathematics Foundation.
In general, for an × 2D wrapped phase array, the weighted minimum -norm phase unwrapping problem is expressed by [14,17] where , and , are the user-defined weights for corresponding differences and indicate where the phase values are more reliable than others. The user-defined weights are given by 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 = 3. , is relatively higher in the areas where the phase changes smoothly. Through the analogous derivation process in [12,26], the minimum -norm solution of (3) eventually entails the solution of the following linear system of equations: , = , ⋅ , +1 − , − Δ , −2 , = 1, 2, . . . , ; = 1, 2, . . . , − 1, That is to say, when the weighted minimum -norm phase unwrapping solution is desired, it can be simpler to find the solution , of (6) instead.
For the setting of the values of data-dependent weights , and , empirically, the weights are always normalized in the interval [0, 1] in practice. Hence, (7) and (8) are modified to incorporate this constraint. To avoid | ⋅ | −2 with 0 ≤ ≤ 2 being extremely great or even infinite, the normalization is defined by Here, setting to be 0.01 radians is a compromise between accuracy and efficiency [12,22].

Sparse
Representation of the Objective. The estimated phase distribution is shown in Figure 1. To simplify the algorithm, concatenating the columns of the phase matrix yields a vector of length MN: The superscript refers to a matrix transpose. The wrapped phase differences along the two directions, Δ and Δ , are also matrices. Their sizes are ( − 1) × and × ( − 1), respectively. As above, we concatenate the columns of each matrix and then merge these two vectors vertically into a single array: The length of is ( − 1) + ( − 1).
Analogous to the matrix equation used in the weighted least-squares phase unwrapping algorithm [18], the system of equations defined by (6) can be represented in matrix form as follows: where Q denotes an × matrix given by A is a matrix consisting of an ( − 1) × upper partition and an ( − 1) × lower partition: where D 1 is an ( − 1) × matrix given by and I is an × identity matrix.
To make the expansion of (13) equal to (6), the matrix W W is constructed as where diag{⋅} puts its arguments in order on the main diagonal. It is easy to see that W W is an If these foregoing matrices are stored and operated using a standard matrix structure, this will consume a significant amount of memory and become unfeasible in standard computers. For instance, given a 256 × 256 phase image, the sizes of A, W W, and Q are 130560×65536, 130560×130560, and 65536 × 65536, respectively. These matrices are obviously too large (exceeding 2 30 ) to be stored and manipulated in a program. Fortunately, all these matrices are sparse and can be represented with a sparse structure. We store only nonzero entries of the matrix together with their indexes. That is, a 3-tuple, ( , , ), is applied to uniquely identify a nonzero entry of the sparse matrix [27,28], where and are the row 4 Computational and Mathematical Methods in Medicine and column indexes, respectively, and denotes the value of the nonzero entry located in ( , ). Moreover, to make the subsequent computation more efficient, the nonzero entries are ordered first by columns and then by rows. Finally, the size of the original sparse matrix should also be stored.
It should be noted that, henceforth, A, W W, and Q will still refer to the original sparse matrices but will be stored and operated, respectively, in the corresponding sparse structures. (14) implies that Q is a large sparse, real, symmetric matrix. Thus, (13) becomes a linear system of equations involving the large sparse symmetric coefficient matrix Q. To solve this system of equations, considerable effort has been devoted over the past three decades [29]. Among these algorithms, the direct solvers that rely on the explicit factorization of the coefficient matrix are widely used owing to their generality and robustness. Especially for a symmetric coefficient matrix, the direct solvers prefer to adopt the factorization [30], a variant of Gaussian elimination that can also be considered as an alternative form of the Cholesky factorization without extracting the square roots [31]. Under this factorization, the matrix Q is decomposed into a lower triangular matrix L, a diagonal matrix D with 1 × 1 and 2 × 2 blocks, and the conjugate transpose of L. Consider

Direct Solver. Equation
Then, given the right-hand side , the estimated phases in (13) can be obtained through a forward elimination followed by the backward substitution [28]. The time complexity is O((MN) 1.5 ) and the memory complexity is O (MNlog(MN)). For more details of the direct solvers, please see [32].

Implementation of the Algorithm.
During the last two decades a number of software packages that implement direct solvers have been developed [33]. If the coefficient matrix is positive definite, CHOLMOD [34] is adopted, because of its rapid computation and relatively small amounts of memory demand. If not, MA57 [35] is strongly recommended.
To accelerate the convergence, the termination condition is set to no residues for the distinctions [12]. First, the distinction is defined by Then, we treat these distinctions similarly to the phase data. Based on the theory that any wrapped phase image can be uniquely unwrapped (with an arbitrary constant offset) if it does not have residues [36], we check whether the , have residues. Given no residues, we unwrap , by a floodfill algorithm [37,38] with a centre start point and no branch cuts.
Step 1. Choose the start pixel as the point in the location of = round( /2) and = round( /2), where round(⋅) rounds its argument to the nearest integer. Its phase value is stored as an unwrapped phase value in the solution matrix.
The four neighbouring pixels are next unwrapped and their unwrapped phase values are placed in the solution matrix. These four pixels are inserted in the unwrapped list.
Step 2. Pick a pixel from the unwrapped list and then eliminate it from the unwrapped list. Unwrap the phase values of its four neighbouring pixels, avoiding pixels that have been unwrapped. Insert these pixels in the unwrapped list and put their unwrapped phase values in the solution matrix.
Step 3. Repeat Step 2 until the unwrapped list becomes empty.
In summary, the detailed procedure of the SDM phase unwrapping algorithm is described as below.
Step 1. Set the iteration time = 0. Set the value of . Initialize the estimated phases.
Step 4. Compute the distinction , by (20). Check , for residues. If there are no residues, continue. If > max , end. Otherwise, set = + 1 and go to Step 2.
Step 5. Unwrap , by the flood-fill algorithm, and form the final estimated phases by (21).

Evaluation.
The following weighted 0 measure is used to evaluate the quality of an unwrapped solution: This counts the ratio of pixels where the gradients of the unwrapped solution mismatch the wrapped phase gradients, which is what we have also minimized with = 0. Therefore, the lower weighted 0 measure indicates the better performance.

Results and Discussion
Several 2D datasets are used to evaluate the performance of the proposed SDM algorithm. The actual number of the iterations of the SDM algorithm varies, depending on the loops after which the distinctions of each example become residue free. The limited maximum number of iterations max is set to be 20. The performance of the proposed algorithm is compared with conventional minimum -norm algorithm [12] and other two widely used 2D algorithms, PUMA [16] and PHUN [10]. Note that we choose = 0 for the first two methods, because 0 -norm minimization is more well behaved and most desired in practice (explained in Section 1). We merely display the best result we can obtain for the PHUN method in each example, because the behaviour of this algorithm is controlled by many parameters and the optimum values of the parameters vary for each dataset. All these methods are implemented in MATLAB (MathWorks, Natick, MA) on a PC (Intel 2 Quad CPU 2.39 GHz).

Simulated Data.
We begin with a 128 × 128 phase dataset with one shear line located horizontally along the 13th row below the upper border. The wrapped phase image is shown in Figure 2(a). Figure 2(b), where the black border is added, shows its residues marked as black dots.
The proposed method, SDM , made no unwrapping errors for this case, producing a solution (see Figure 2(c)) exactly the same as the true phase. However, while it obtained the smallest weighted 0 measure for this example, it converged more slowly than the PUMA and PHUN methods (see Figure 6). The conventional minimum -norm and PUMA methods both yielded extra short vertical shear lines marked with red arrows in the upper parts of their results (see Figures 2(d) and 2(e)). Figure 2(f) shows the solution achieved by the PHUN method with five unexpected crooked shear lines and a black spot of outliers whose values are very low (all marked with red arrows). Figure 3(a) shows a 257 × 257 wrapped phase shear image with Gaussian noise (the signal-to-noise ratio is 0.8379 dB). Its residue distribution is shown in Figure 3(b).
The result of the SDM algorithm in Figure 3(c) is composed of a two-planar surface tearing along the median shear line as expected, except for some protrusions on this line where the phases are severely corrupted by both the noise and the object discontinuities. The conventional minimum -norm method offers a slightly worse result, shown in Figure 3 conventional minimum -norm methods are almost the same (see Figure 6(a)). However, it is noted in this case that the former converged much faster than the latter. The unwrapped phase image of the PUMA method in Figure 3(e) has some small-scale, anomalous "layer" artifacts around the median shear line. The PHUN method generated an incorrect unwrapped phase image with the left half full of "layer" artifacts that propagated from the median line to the surrounding regions or even the image border. In summary, these two simulation examples, both having object discontinuities lying in the shear line, demonstrate the discontinuity-preserving ability of the SDM algorithm. Additionally, the SDM algorithm did not produce extraundesired discontinuities that appeared as shear lines or "layer" artifacts. . Figures 4(a) and 4(b) show the magnitudes and phases of a 44 × 44 displacement encoded MR heart dataset [39], respectively.

MR Data
The results of the SDM , the conventional minimum -norm, and the PUMA algorithms are shown in Figures   4(d), 4(e), and 4(f). There is no significant visible difference between these three images. In these unwrapped phase images, the rough shape of the heart (compared to the magnitude picture in Figure 4(a)) is dimly visible. We then examine the corresponding discontinuity maps (added black borders) that are the distributions of the discontinuities where pixels (marked in black) differ from a neighbouring pixel by more than radians. These maps in Figures 4(h), 4(i), and 4(j) have little differences. In addition, the weighted 0 measures, depicted in Figure 6(a), of these three methods are similar. The extremely fast algorithm, PHUN, offered a solution involving many black spots (see Figure 4(g)) that correspond roughly to the locations of the residues (see Figure 4(c)). Also, it created a large number of undesired discontinuities (see Figure 5(j)).
An MR head example [12] is shown in Figure 5. The SDM , the conventional minimum -norm, and the PUMA algorithms all produced plausible unwrapped phase images, shown in Figures 5(c), 5(d), and 5(e). Moreover, the discontinuity maps (see Figures 5(g), 5(h), and 5(i)) and weighted 0 measures (see Figure 6(a)) of these three Computational and Mathematical Methods in Medicine From the horizontal comparison of the execute time, it can be concluded that the execute time of the SDM algorithm depended partly on the size of the phase dataset. For example, the SDM algorithm converged in 0.5 seconds for the 44×44 heart dataset, while, for the 257×257 simulated dataset, it took 31 seconds.
The final example is the transverse section of a 5-slice MR knee dataset, as shown in Figure 7. The size of each slice is 256 × 256. The measurements were made with a 0.5T MRI scanner (Ningbo Xingaoyi Co., LTD., China).
Excellent results with no phase errors, shown in the third row of Figure 7, are derived from the unwrapping by the SDM algorithm. The conventional minimum -norm algorithm worked well, except for generating two undesirable shear lines in the fourth slice. The PUMA algorithm yielded undesired results in the third and fourth slices with some shear lines but successfully unwrapped the other slices. The PHUN method made some unwrapping errors in all slices, shear lines in the first, third, and fourth slices, and white patches of outliers, whose values are very high, in the second and fifth slices. All the truncations and outliers are marked with yellow and green arrows, respectively, in Figure 7.
As above, the weighted 0 measures and execute time of all the methods for this multislice dataset are compared in Figures 8(a) and 8(b), respectively. The SDM method returned the smallest weighted 0 measures for all slices. Additionally, it yielded these results in 19.62 ± 0.99 seconds, slower than the PUMA and PHUN methods.

Conclusions
In this work, we developed a sparse-representation-based direct minimum -norm (SDM ) algorithm for the 2D phase unwrapping. The user-defined weights are introduced in the objective function to improve the discontinuity-preserving ability of the SDM algorithm. Furthermore, the sparse structures are used to represent the matrices involved in the objective function to accelerate the computation and decrease the memory space. Finally, the SDM algorithms computed effectively and efficiently by employing direct solvers.
The proposed algorithm does produce excellent, reliable results with a very small weighted 0 measure; it even allows phase images with large discontinuities through the whole phases to be unwrapped correctly. Moreover, benefiting from using the sparse representation and well-developed direct solvers, the SDM method converges much faster than the conventional minimum -norm method. However, it is not fast enough. Further research can be devoted to reducing the execution time.