# Sparse-Representation-Based Direct MinimumL

Authors: Wei He; Ling Xia; Feng Liu
Journal: Computational and Mathematical Methods in Medicine
 (2014)

---

## Abstract

A sparse-representation-based direct minimumL
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.


---

## Body

## 1. 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
k
π, where k 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 B
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, three-point 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 each2
×
2 closed loop in the two-dimensional (2D) array, as shown in Figure 1. Consider
(1)
Δ
ψ
i
,
j
x
-
Δ
ψ
i
,
j
y
-
Δ
ψ
i
,
j
+
1
x
+
Δ
ψ
i
+
1
,
j
y
=
±
2
π
,
where Δ
ψ
i
,
j
x and Δ
ψ
i
,
j
y are defined to be the gradients of the wrapped phases:
(2)
Δ
ψ
i
,
j
x
=
ω
{
ψ
i
+
1
,
j
-
ψ
i
,
j
}
,
Δ
ψ
i
,
j
y
=
ω
{
ψ
i
,
j
+
1
-
ψ
i
,
j
}
.Figure 1
Residue calculation.In the literature, numerous two-dimensional (2D) phase unwrapping approaches have been developed. They can be classified into four categories: path-following [6–11], minimum-norm [12–18], Bayesian/regularization [19, 20], and parametric modeling [21] methods. Our method belongs to the second category; thus, we briefly introduce the minimum-norm methods.The minimum-norm methods estimate the unwrapped phases by minimizing theL
p-norm of the differences between the gradients of the wrapped and the unwrapped phases. With p
=
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 L
1-norm algorithms [13, 23] have a better ability than the L
2-norm ones to preserve the discontinuities. As p approaches zero, the minimum L
p-norm method tends to obtain a more reliable result [24, 25]. Thus, the L
0-norm minimization is widely accepted as the most desirable in practice. However, the L
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 L
p-norm method [12] converts the L
p-norm minimization problem into a generalized matrix equation with a flexible option of p. 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 dual-iterative structure (an iteration structure embedded in an outer one), it can be computationally intense.In this work, a sparse-representation-based direct minimumL
p-norm (SDM
L
P) algorithm is proposed. It introduces the user-defined weights into the generalized matrix equation of the conventional minimum L
p-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 L
p-norm method for MRI phase unwrapping.


## 2. Materials and Methods

### 2.1. Mathematics Foundation

In general, for anM
×
N 2D wrapped phase array, the weighted minimum L
p-norm phase unwrapping problem is expressed by [14, 17]
(3)
argmin
φ
∑
i
=
1
M
-
1
∑
j
=
1
N
w
i
,
j
x
·
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
p
+
∑
i
=
1
M
∑
j
=
1
N
-
1
w
i
,
j
y
·
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
p
,
where w
i
,
j
x and w
i
,
j
y 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
(4)
h
h
h
h
h
h
h
w
i
,
j
x
=
min
⁡
(
w
i
+
1
,
j
,
w
i
,
j
)
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
h
h
h
h
h
h
h
w
i
,
j
y
=
min
⁡
(
w
i
,
j
+
1
,
w
i
,
j
)
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
,
where
(5)
w
i
,
j
=
max
⁡
(
z
i
,
j
)
-
z
i
,
j
max
⁡
(
z
i
,
j
)
-
min
⁡
(
z
i
,
j
)
,
z
i
,
j
=
∑
(
Δ
ψ
i
,
j
x
-
Δ
ψ
m
,
n
x
¯
)
2
+
∑
(
Δ
ψ
i
,
j
y
-
Δ
ψ
m
,
n
y
¯
)
2
l
2
.
For each sum, the indexes (
i
,
j
) cover over the l
×
l window centered at the pixel (
m
,
n
). The terms Δ
ψ
i
,
j
x and Δ
ψ
i
,
j
y are the wrapped phase gradients in the l
×
l windows and Δ
ψ
m
,
n
x
¯ and Δ
ψ
m
,
n
y
¯ are the averages of these wrapped phase gradients. In this paper l
=
3. w
i
,
j is relatively higher in the areas where the phase changes smoothly.Through the analogous derivation process in [12, 26], the minimum L
p-norm solution of (3) eventually entails the solution of the following linear system of equations:
(6)
R
i
,
j
(
φ
i
+
1
,
j
-
φ
i
,
j
)
+
C
i
,
j
(
φ
i
,
j
+
1
-
φ
i
,
j
)
-
R
i
,
j
(
φ
i
,
j
-
φ
i
-
1
,
j
)
-
C
i
,
j
(
φ
i
,
j
-
φ
i
,
j
-
1
)
=
ρ
i
,
j
,
where
(7)
R
i
,
j
=
w
i
,
j
x
·
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
p
-
2
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
(8)
C
i
,
j
=
w
i
,
j
y
·
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
p
-
2
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
,
(9)
ρ
i
,
j
=
R
i
,
j
Δ
ψ
i
,
j
x
-
R
i
,
j
Δ
ψ
i
-
1
,
j
x
+
C
i
,
j
Δ
ψ
i
,
j
y
-
C
i
,
j
Δ
ψ
i
,
j
-
1
y
.That is to say, when the weighted minimumL
p-norm phase unwrapping solution is desired, it can be simpler to find the solution φ
i
,
j of (6) instead.For the setting of the values of data-dependent weightsR
i
,
j and C
i
,
j 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 |
·
|
p
-
2 with 0
≤
p
≤
2 being extremely great or even infinite, the normalization is defined by
(10)
R
i
,
j
=
w
i
,
j
x
·
α
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
2
-
p
+
α
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
C
i
,
j
=
w
i
,
j
y
·
α
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
2
-
p
+
α
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
.
Here, setting α to be 0.01 radians is a compromise between accuracy and efficiency [12, 22].


### 2.2. Sparse-Representation-Based Direct MinimumL
p-Norm (SDM
L
P) Algorithm

#### 2.2.1. Sparse Representation of the Objective

The estimated phase distribution is shown in Figure1. To simplify the algorithm, concatenating the columns of the phase matrix yields a vector of lengthMN:
(11)
ϕ
=
[
φ
1,1
,
φ
2,1
,
…
,
φ
M
,
1
,
φ
1,2
,
φ
2,2
,
…
,
φ
M
,
2
,
…
,
φ
1
,
N
,
…
,
φ
M
,
N
]
T
.
The superscript T refers to a matrix transpose.The wrapped phase differences along the two directions,Δ
ψ
x and Δ
ψ
y, are also matrices. Their sizes are (
M
-
1
)
×
N and M
×
(
N
-
1
), respectively. As above, we concatenate the columns of each matrix and then merge these two vectors vertically into a single array:
(12)
d
=
[
Δ
ψ
1,1
x
,
…
,
Δ
ψ
M
-
1,1
x
,
…
,
Δ
ψ
1
,
N
x
,
…
,
Δ
ψ
M
-
1
,
N
x
,
Δ
ψ
1,1
y
,
…
,
Δ
ψ
M
,
1
y
,
…
,
Δ
ψ
1
,
N
-
1
y
,
…
,
Δ
ψ
M
,
N
-
1
y
]
T
.
The length of d is (
M
-
1
)
N
+
M
(
N
-
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:
(13)
Q
ϕ
=
s
,
where Q denotes an M
N
×
M
N matrix given by
(14)
Q
=
A
T
W
T
W
A
,
(15)
s
=
A
T
W
T
W
d
.A is a matrix consisting of an (
M
-
1
)
N
×
M
N upper partition and an M
(
N
-
1
)
×
M
N lower partition:
(16)
A
=
[
D
1
0
0
⋯
0
0
0
D
1
0
⋯
0
0
0
0
D
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
0
D
1
I
-
I
0
⋯
0
0
0
I
-
I
⋯
0
0
0
0
I
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
I
-
I
]
,
where D
1 is an (
M
-
1
)
×
M matrix given by
(17)
D
1
=
[
1
-
1
0
⋯
0
0
0
1
-
1
⋯
0
0
0
0
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
1
-
1
]
and I is an M
×
M identity matrix.To make the expansion of (13) equal to (6), the matrix W
T
W is constructed as
(18)
W
T
W
=
diag
⁡
{
R
1,1
,
…
,
R
M
-
1,1
,
…
,
R
1
,
N
,
…
,
R
M
-
1
,
N
,
000000000000
C
1,1
,
…
,
C
M
,
1
,
…
C
1
,
N
-
1
,
…
,
C
M
,
N
-
1
}
,
where diag
⁡
{
·
} puts its arguments in order on the main diagonal. It is easy to see that W
T
W is an [
(
M
-
1
)
×
N
+
M
×
(
N
-
1
)
]
2 matrix.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 a256
×
256 phase image, the sizes of A, W
T
W
, and Q are 130560
×
65536, 130560
×
130560, and 65536
×
65536, respectively. These matrices are obviously too large (exceeding 230) 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, (
i
,
j
,
a
i
j
), is applied to uniquely identify a nonzero entry of the sparse matrix [27, 28], where i and j are the row and column indexes, respectively, and a
i
j denotes the value of the nonzero entry located in (
i
,
j
). 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
T
W, and Q will still refer to the original sparse matrices but will be stored and operated, respectively, in the corresponding sparse structures.


#### 2.2.2. Direct Solver

Equation (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 L
D
L
T 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
(19)
Q
=
L
D
L
T
.Then, given the right-hand sides, 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].


#### 2.2.3. 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
(20)
E
i
,
j
=
ω
{
ψ
i
,
j
-
φ
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.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 E
i
,
j have residues. Given no residues, we unwrap E
i
,
j by a flood-fill algorithm [37, 38] with a centre start point and no branch cuts.Step 1.
Choose the start pixel as the point in the location ofi = round(M
/
2) and j = round(N
/
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 Step2 until the unwrapped list becomes empty.UnwrappingE
i
,
j is signified by ω
-
1
{
E
i
,
j
}, which is added back to the estimated phase to obtain the final result. Consider
(21)
φ
i
,
j
final
=
φ
i
,
j
+
ω
-
1
{
E
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.In summary, the detailed procedure of theSDM
L
P phase unwrapping algorithm is described as below.Step 1.
Set the iteration timek
=
0. Set the value of p. Initialize the estimated phases.Step 2.
Compute the weightsR
i
,
j and C
i
,
j from (10).Step 3.
ComputeQ and s from (14) and (15); then solve (13) by the direct solver.Step 4.
Compute the distinctionE
i
,
j by (20). Check E
i
,
j for residues. If there are no residues, continue. If k
>
k
max
⁡, end. Otherwise, set k
=
k
+
1 and go to Step 2.Step 5.
UnwrapE
i
,
j by the flood-fill algorithm, and form the final estimated phases by (21).


#### 2.2.4. Evaluation

The following weightedL
0 measure is used to evaluate the quality of an unwrapped solution:
(22)
ρ
=
1
M
N
[
∑
i
=
1
M
-
1
∑
j
=
1
N
w
i
,
j
x
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
0
h
h
h
h
+
∑
i
=
1
M
∑
j
=
1
N
-
1
w
i
,
j
y
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
0
]
.
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 p
=
0. Therefore, the lower weighted L
0 measure indicates the better performance.


## 2.1. Mathematics Foundation

In general, for anM
×
N 2D wrapped phase array, the weighted minimum L
p-norm phase unwrapping problem is expressed by [14, 17]
(3)
argmin
φ
∑
i
=
1
M
-
1
∑
j
=
1
N
w
i
,
j
x
·
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
p
+
∑
i
=
1
M
∑
j
=
1
N
-
1
w
i
,
j
y
·
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
p
,
where w
i
,
j
x and w
i
,
j
y 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
(4)
h
h
h
h
h
h
h
w
i
,
j
x
=
min
⁡
(
w
i
+
1
,
j
,
w
i
,
j
)
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
h
h
h
h
h
h
h
w
i
,
j
y
=
min
⁡
(
w
i
,
j
+
1
,
w
i
,
j
)
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
,
where
(5)
w
i
,
j
=
max
⁡
(
z
i
,
j
)
-
z
i
,
j
max
⁡
(
z
i
,
j
)
-
min
⁡
(
z
i
,
j
)
,
z
i
,
j
=
∑
(
Δ
ψ
i
,
j
x
-
Δ
ψ
m
,
n
x
¯
)
2
+
∑
(
Δ
ψ
i
,
j
y
-
Δ
ψ
m
,
n
y
¯
)
2
l
2
.
For each sum, the indexes (
i
,
j
) cover over the l
×
l window centered at the pixel (
m
,
n
). The terms Δ
ψ
i
,
j
x and Δ
ψ
i
,
j
y are the wrapped phase gradients in the l
×
l windows and Δ
ψ
m
,
n
x
¯ and Δ
ψ
m
,
n
y
¯ are the averages of these wrapped phase gradients. In this paper l
=
3. w
i
,
j is relatively higher in the areas where the phase changes smoothly.Through the analogous derivation process in [12, 26], the minimum L
p-norm solution of (3) eventually entails the solution of the following linear system of equations:
(6)
R
i
,
j
(
φ
i
+
1
,
j
-
φ
i
,
j
)
+
C
i
,
j
(
φ
i
,
j
+
1
-
φ
i
,
j
)
-
R
i
,
j
(
φ
i
,
j
-
φ
i
-
1
,
j
)
-
C
i
,
j
(
φ
i
,
j
-
φ
i
,
j
-
1
)
=
ρ
i
,
j
,
where
(7)
R
i
,
j
=
w
i
,
j
x
·
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
p
-
2
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
(8)
C
i
,
j
=
w
i
,
j
y
·
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
p
-
2
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
,
(9)
ρ
i
,
j
=
R
i
,
j
Δ
ψ
i
,
j
x
-
R
i
,
j
Δ
ψ
i
-
1
,
j
x
+
C
i
,
j
Δ
ψ
i
,
j
y
-
C
i
,
j
Δ
ψ
i
,
j
-
1
y
.That is to say, when the weighted minimumL
p-norm phase unwrapping solution is desired, it can be simpler to find the solution φ
i
,
j of (6) instead.For the setting of the values of data-dependent weightsR
i
,
j and C
i
,
j 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 |
·
|
p
-
2 with 0
≤
p
≤
2 being extremely great or even infinite, the normalization is defined by
(10)
R
i
,
j
=
w
i
,
j
x
·
α
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
2
-
p
+
α
,
i
=
1,2
,
…
,
M
-
1
;
j
=
1,2
,
…
,
N
,
C
i
,
j
=
w
i
,
j
y
·
α
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
2
-
p
+
α
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
-
1
.
Here, setting α to be 0.01 radians is a compromise between accuracy and efficiency [12, 22].


## 2.2. Sparse-Representation-Based Direct MinimumL
p-Norm (SDM
L
P) Algorithm

### 2.2.1. Sparse Representation of the Objective

The estimated phase distribution is shown in Figure1. To simplify the algorithm, concatenating the columns of the phase matrix yields a vector of lengthMN:
(11)
ϕ
=
[
φ
1,1
,
φ
2,1
,
…
,
φ
M
,
1
,
φ
1,2
,
φ
2,2
,
…
,
φ
M
,
2
,
…
,
φ
1
,
N
,
…
,
φ
M
,
N
]
T
.
The superscript T refers to a matrix transpose.The wrapped phase differences along the two directions,Δ
ψ
x and Δ
ψ
y, are also matrices. Their sizes are (
M
-
1
)
×
N and M
×
(
N
-
1
), respectively. As above, we concatenate the columns of each matrix and then merge these two vectors vertically into a single array:
(12)
d
=
[
Δ
ψ
1,1
x
,
…
,
Δ
ψ
M
-
1,1
x
,
…
,
Δ
ψ
1
,
N
x
,
…
,
Δ
ψ
M
-
1
,
N
x
,
Δ
ψ
1,1
y
,
…
,
Δ
ψ
M
,
1
y
,
…
,
Δ
ψ
1
,
N
-
1
y
,
…
,
Δ
ψ
M
,
N
-
1
y
]
T
.
The length of d is (
M
-
1
)
N
+
M
(
N
-
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:
(13)
Q
ϕ
=
s
,
where Q denotes an M
N
×
M
N matrix given by
(14)
Q
=
A
T
W
T
W
A
,
(15)
s
=
A
T
W
T
W
d
.A is a matrix consisting of an (
M
-
1
)
N
×
M
N upper partition and an M
(
N
-
1
)
×
M
N lower partition:
(16)
A
=
[
D
1
0
0
⋯
0
0
0
D
1
0
⋯
0
0
0
0
D
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
0
D
1
I
-
I
0
⋯
0
0
0
I
-
I
⋯
0
0
0
0
I
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
I
-
I
]
,
where D
1 is an (
M
-
1
)
×
M matrix given by
(17)
D
1
=
[
1
-
1
0
⋯
0
0
0
1
-
1
⋯
0
0
0
0
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
1
-
1
]
and I is an M
×
M identity matrix.To make the expansion of (13) equal to (6), the matrix W
T
W is constructed as
(18)
W
T
W
=
diag
⁡
{
R
1,1
,
…
,
R
M
-
1,1
,
…
,
R
1
,
N
,
…
,
R
M
-
1
,
N
,
000000000000
C
1,1
,
…
,
C
M
,
1
,
…
C
1
,
N
-
1
,
…
,
C
M
,
N
-
1
}
,
where diag
⁡
{
·
} puts its arguments in order on the main diagonal. It is easy to see that W
T
W is an [
(
M
-
1
)
×
N
+
M
×
(
N
-
1
)
]
2 matrix.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 a256
×
256 phase image, the sizes of A, W
T
W
, and Q are 130560
×
65536, 130560
×
130560, and 65536
×
65536, respectively. These matrices are obviously too large (exceeding 230) 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, (
i
,
j
,
a
i
j
), is applied to uniquely identify a nonzero entry of the sparse matrix [27, 28], where i and j are the row and column indexes, respectively, and a
i
j denotes the value of the nonzero entry located in (
i
,
j
). 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
T
W, and Q will still refer to the original sparse matrices but will be stored and operated, respectively, in the corresponding sparse structures.


### 2.2.2. Direct Solver

Equation (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 L
D
L
T 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
(19)
Q
=
L
D
L
T
.Then, given the right-hand sides, 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].


### 2.2.3. 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
(20)
E
i
,
j
=
ω
{
ψ
i
,
j
-
φ
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.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 E
i
,
j have residues. Given no residues, we unwrap E
i
,
j by a flood-fill algorithm [37, 38] with a centre start point and no branch cuts.Step 1.
Choose the start pixel as the point in the location ofi = round(M
/
2) and j = round(N
/
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 Step2 until the unwrapped list becomes empty.UnwrappingE
i
,
j is signified by ω
-
1
{
E
i
,
j
}, which is added back to the estimated phase to obtain the final result. Consider
(21)
φ
i
,
j
final
=
φ
i
,
j
+
ω
-
1
{
E
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.In summary, the detailed procedure of theSDM
L
P phase unwrapping algorithm is described as below.Step 1.
Set the iteration timek
=
0. Set the value of p. Initialize the estimated phases.Step 2.
Compute the weightsR
i
,
j and C
i
,
j from (10).Step 3.
ComputeQ and s from (14) and (15); then solve (13) by the direct solver.Step 4.
Compute the distinctionE
i
,
j by (20). Check E
i
,
j for residues. If there are no residues, continue. If k
>
k
max
⁡, end. Otherwise, set k
=
k
+
1 and go to Step 2.Step 5.
UnwrapE
i
,
j by the flood-fill algorithm, and form the final estimated phases by (21).


### 2.2.4. Evaluation

The following weightedL
0 measure is used to evaluate the quality of an unwrapped solution:
(22)
ρ
=
1
M
N
[
∑
i
=
1
M
-
1
∑
j
=
1
N
w
i
,
j
x
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
0
h
h
h
h
+
∑
i
=
1
M
∑
j
=
1
N
-
1
w
i
,
j
y
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
0
]
.
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 p
=
0. Therefore, the lower weighted L
0 measure indicates the better performance.


## 2.2.1. Sparse Representation of the Objective

The estimated phase distribution is shown in Figure1. To simplify the algorithm, concatenating the columns of the phase matrix yields a vector of lengthMN:
(11)
ϕ
=
[
φ
1,1
,
φ
2,1
,
…
,
φ
M
,
1
,
φ
1,2
,
φ
2,2
,
…
,
φ
M
,
2
,
…
,
φ
1
,
N
,
…
,
φ
M
,
N
]
T
.
The superscript T refers to a matrix transpose.The wrapped phase differences along the two directions,Δ
ψ
x and Δ
ψ
y, are also matrices. Their sizes are (
M
-
1
)
×
N and M
×
(
N
-
1
), respectively. As above, we concatenate the columns of each matrix and then merge these two vectors vertically into a single array:
(12)
d
=
[
Δ
ψ
1,1
x
,
…
,
Δ
ψ
M
-
1,1
x
,
…
,
Δ
ψ
1
,
N
x
,
…
,
Δ
ψ
M
-
1
,
N
x
,
Δ
ψ
1,1
y
,
…
,
Δ
ψ
M
,
1
y
,
…
,
Δ
ψ
1
,
N
-
1
y
,
…
,
Δ
ψ
M
,
N
-
1
y
]
T
.
The length of d is (
M
-
1
)
N
+
M
(
N
-
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:
(13)
Q
ϕ
=
s
,
where Q denotes an M
N
×
M
N matrix given by
(14)
Q
=
A
T
W
T
W
A
,
(15)
s
=
A
T
W
T
W
d
.A is a matrix consisting of an (
M
-
1
)
N
×
M
N upper partition and an M
(
N
-
1
)
×
M
N lower partition:
(16)
A
=
[
D
1
0
0
⋯
0
0
0
D
1
0
⋯
0
0
0
0
D
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
0
D
1
I
-
I
0
⋯
0
0
0
I
-
I
⋯
0
0
0
0
I
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
I
-
I
]
,
where D
1 is an (
M
-
1
)
×
M matrix given by
(17)
D
1
=
[
1
-
1
0
⋯
0
0
0
1
-
1
⋯
0
0
0
0
1
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
0
0
0
⋯
1
-
1
]
and I is an M
×
M identity matrix.To make the expansion of (13) equal to (6), the matrix W
T
W is constructed as
(18)
W
T
W
=
diag
⁡
{
R
1,1
,
…
,
R
M
-
1,1
,
…
,
R
1
,
N
,
…
,
R
M
-
1
,
N
,
000000000000
C
1,1
,
…
,
C
M
,
1
,
…
C
1
,
N
-
1
,
…
,
C
M
,
N
-
1
}
,
where diag
⁡
{
·
} puts its arguments in order on the main diagonal. It is easy to see that W
T
W is an [
(
M
-
1
)
×
N
+
M
×
(
N
-
1
)
]
2 matrix.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 a256
×
256 phase image, the sizes of A, W
T
W
, and Q are 130560
×
65536, 130560
×
130560, and 65536
×
65536, respectively. These matrices are obviously too large (exceeding 230) 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, (
i
,
j
,
a
i
j
), is applied to uniquely identify a nonzero entry of the sparse matrix [27, 28], where i and j are the row and column indexes, respectively, and a
i
j denotes the value of the nonzero entry located in (
i
,
j
). 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
T
W, and Q will still refer to the original sparse matrices but will be stored and operated, respectively, in the corresponding sparse structures.


## 2.2.2. Direct Solver

Equation (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 L
D
L
T 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
(19)
Q
=
L
D
L
T
.Then, given the right-hand sides, 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].


## 2.2.3. 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
(20)
E
i
,
j
=
ω
{
ψ
i
,
j
-
φ
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.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 E
i
,
j have residues. Given no residues, we unwrap E
i
,
j by a flood-fill algorithm [37, 38] with a centre start point and no branch cuts.Step 1.
Choose the start pixel as the point in the location ofi = round(M
/
2) and j = round(N
/
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 Step2 until the unwrapped list becomes empty.UnwrappingE
i
,
j is signified by ω
-
1
{
E
i
,
j
}, which is added back to the estimated phase to obtain the final result. Consider
(21)
φ
i
,
j
final
=
φ
i
,
j
+
ω
-
1
{
E
i
,
j
}
,
i
=
1,2
,
…
,
M
;
j
=
1,2
,
…
,
N
.In summary, the detailed procedure of theSDM
L
P phase unwrapping algorithm is described as below.Step 1.
Set the iteration timek
=
0. Set the value of p. Initialize the estimated phases.Step 2.
Compute the weightsR
i
,
j and C
i
,
j from (10).Step 3.
ComputeQ and s from (14) and (15); then solve (13) by the direct solver.Step 4.
Compute the distinctionE
i
,
j by (20). Check E
i
,
j for residues. If there are no residues, continue. If k
>
k
max
⁡, end. Otherwise, set k
=
k
+
1 and go to Step 2.Step 5.
UnwrapE
i
,
j by the flood-fill algorithm, and form the final estimated phases by (21).


## 2.2.4. Evaluation

The following weightedL
0 measure is used to evaluate the quality of an unwrapped solution:
(22)
ρ
=
1
M
N
[
∑
i
=
1
M
-
1
∑
j
=
1
N
w
i
,
j
x
|
φ
i
+
1
,
j
-
φ
i
,
j
-
Δ
ψ
i
,
j
x
|
0
h
h
h
h
+
∑
i
=
1
M
∑
j
=
1
N
-
1
w
i
,
j
y
|
φ
i
,
j
+
1
-
φ
i
,
j
-
Δ
ψ
i
,
j
y
|
0
]
.
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 p
=
0. Therefore, the lower weighted L
0 measure indicates the better performance.


## 3. Results and Discussion

Several 2D datasets are used to evaluate the performance of the proposedSDM
L
P algorithm. The actual number of the iterations of the SDM
L
P algorithm varies, depending on the loops after which the distinctions of each example become residue free. The limited maximum number of iterations k
max
⁡ is set to be 20.The performance of the proposed algorithm is compared with conventional minimumL
p-norm algorithm [12] and other two widely used 2D algorithms, PUMA [16] and PHUN [10]. Note that we choose p
=
0 for the first two methods, because L
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).


### 3.1. Simulated Data

We begin with a128
×
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.A128
×
128 phase dataset with shear, (a) its wrapped phase, and (b) its residue distribution including 6 residues. Unwrapped phase image using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods. Unwrapping errors are marked with red arrows.
(a)
(b)
(c)
(d)
(e)
(f)The proposed method,SDM
L
P, 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 L
0 measure for this example, it converged more slowly than the PUMA and PHUN methods (see Figure 6). The conventional minimum L
p-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).Figure3(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).A257
×
257 noisy phase dataset with shear, (a) its wrapped phase, and (b) its residue distribution including 150 residues. Unwrapped phase image using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)The result of theSDM
L
P 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 L
p-norm method offers a slightly worse result, shown in Figure 3(d), where minor errors (marked with red arrows) occurred near the top and bottom of the median shear line. So the weighted L
0 measure values of the SDM
L
P and conventional minimum L
p-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 theSDM
L
P algorithm. Additionally, the SDM
L
P algorithm did not produce extra-undesired discontinuities that appeared as shear lines or “layer” artifacts.


### 3.2. MR Data

Figures4(a) and 4(b) show the magnitudes and phases of a 44
×
44 displacement encoded MR heart dataset [39], respectively.A44
×
44 displacement encoded MR heart dataset, (a) its magnitudes, (b) its wrapped phases, and (c) its residue distribution including 210 residues. Unwrapped phase image using the (d) SDM
L
P, (e) conventional minimum L
p-norm, (f) PUMA, and (g) PHUN methods. Discontinuities in the corresponding unwrapped phase images got by the (h) SDM
L
P, (i) conventional minimum L
p-norm, (j) PUMA, and (k) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)The results of theSDM
L
p, the conventional minimum L
p-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 L
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)).A256
×
256 MR head dataset, (a) its wrapped phases, and (b) its residue distribution including 1929 residues. Unwrapped phase images using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods. Discontinuities in the corresponding unwrapped phase images got by the (g) SDM
L
P, (h) conventional minimum L
p-norm, (i) PUMA, and (j) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)A comparison of the (a) weightedL
0 measures and (b) execute time of the SDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods in unwrapping the preceding four datasets.
(a)
(b)An MR head example [12] is shown in Figure 5. The SDM
L
p, the conventional minimum L
p-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 L
0 measures (see Figure 6(a)) of these three methods are almost the same. The PHUN method failed to unwrap this dataset correctly (see Figure 5(f)). Its result has a large number of undesired discontinuities (see Figure 5(j)).The weightedL
0 measures and execute time of all methods in the preceding four examples are compared in Figures 6(a) and 6(b), respectively. For a better visual effect, both of the ordinate axes are in a logarithmic scale. The SDM
L
p algorithm achieved the smallest weighted L
0 measures for all four examples above. It significantly reduced the computational time compared to the conventional minimum L
p-norm method. However, the algorithm did not converge fast enough compared with the PUMA and PHUN algorithms. From the horizontal comparison of the execute time, it can be concluded that the execute time of the SDM
L
p algorithm depended partly on the size of the phase dataset. For example, the SDM
L
P 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 Figure7. The size of each slice is 256
×
256. The measurements were made with a 0.5T MRI scanner (Ningbo Xingaoyi Co., LTD., China).Figure 7
The transverse sections of a 5-slice MR knee dataset. The magnitude and wrapped phase images are shown in the first and second rows, respectively. Unwrapping results of theSDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods are followed in a top-down order. Undesirable shear lines are marked with yellow arrows. The patches of outliers are marked with green arrows.Excellent results with no phase errors, shown in the third row of Figure7, are derived from the unwrapping by the SDM
L
P algorithm. The conventional minimum L
p-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 weightedL
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
L
P method returned the smallest weighted L
0 measures for all slices. Additionally, it yielded these results in 19.62
±
0.99 seconds, slower than the PUMA and PHUN methods.A comparison of the (a) weightedL
0 measures and (b) execute time of the SDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods in unwrapping the multislice dataset.
(a)
(b)


## 3.1. Simulated Data

We begin with a128
×
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.A128
×
128 phase dataset with shear, (a) its wrapped phase, and (b) its residue distribution including 6 residues. Unwrapped phase image using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods. Unwrapping errors are marked with red arrows.
(a)
(b)
(c)
(d)
(e)
(f)The proposed method,SDM
L
P, 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 L
0 measure for this example, it converged more slowly than the PUMA and PHUN methods (see Figure 6). The conventional minimum L
p-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).Figure3(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).A257
×
257 noisy phase dataset with shear, (a) its wrapped phase, and (b) its residue distribution including 150 residues. Unwrapped phase image using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)The result of theSDM
L
P 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 L
p-norm method offers a slightly worse result, shown in Figure 3(d), where minor errors (marked with red arrows) occurred near the top and bottom of the median shear line. So the weighted L
0 measure values of the SDM
L
P and conventional minimum L
p-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 theSDM
L
P algorithm. Additionally, the SDM
L
P algorithm did not produce extra-undesired discontinuities that appeared as shear lines or “layer” artifacts.


## 3.2. MR Data

Figures4(a) and 4(b) show the magnitudes and phases of a 44
×
44 displacement encoded MR heart dataset [39], respectively.A44
×
44 displacement encoded MR heart dataset, (a) its magnitudes, (b) its wrapped phases, and (c) its residue distribution including 210 residues. Unwrapped phase image using the (d) SDM
L
P, (e) conventional minimum L
p-norm, (f) PUMA, and (g) PHUN methods. Discontinuities in the corresponding unwrapped phase images got by the (h) SDM
L
P, (i) conventional minimum L
p-norm, (j) PUMA, and (k) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)The results of theSDM
L
p, the conventional minimum L
p-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 L
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)).A256
×
256 MR head dataset, (a) its wrapped phases, and (b) its residue distribution including 1929 residues. Unwrapped phase images using the (c) SDM
L
P, (d) conventional minimum L
p-norm, (e) PUMA, and (f) PHUN methods. Discontinuities in the corresponding unwrapped phase images got by the (g) SDM
L
P, (h) conventional minimum L
p-norm, (i) PUMA, and (j) PHUN methods.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)A comparison of the (a) weightedL
0 measures and (b) execute time of the SDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods in unwrapping the preceding four datasets.
(a)
(b)An MR head example [12] is shown in Figure 5. The SDM
L
p, the conventional minimum L
p-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 L
0 measures (see Figure 6(a)) of these three methods are almost the same. The PHUN method failed to unwrap this dataset correctly (see Figure 5(f)). Its result has a large number of undesired discontinuities (see Figure 5(j)).The weightedL
0 measures and execute time of all methods in the preceding four examples are compared in Figures 6(a) and 6(b), respectively. For a better visual effect, both of the ordinate axes are in a logarithmic scale. The SDM
L
p algorithm achieved the smallest weighted L
0 measures for all four examples above. It significantly reduced the computational time compared to the conventional minimum L
p-norm method. However, the algorithm did not converge fast enough compared with the PUMA and PHUN algorithms. From the horizontal comparison of the execute time, it can be concluded that the execute time of the SDM
L
p algorithm depended partly on the size of the phase dataset. For example, the SDM
L
P 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 Figure7. The size of each slice is 256
×
256. The measurements were made with a 0.5T MRI scanner (Ningbo Xingaoyi Co., LTD., China).Figure 7
The transverse sections of a 5-slice MR knee dataset. The magnitude and wrapped phase images are shown in the first and second rows, respectively. Unwrapping results of theSDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods are followed in a top-down order. Undesirable shear lines are marked with yellow arrows. The patches of outliers are marked with green arrows.Excellent results with no phase errors, shown in the third row of Figure7, are derived from the unwrapping by the SDM
L
P algorithm. The conventional minimum L
p-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 weightedL
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
L
P method returned the smallest weighted L
0 measures for all slices. Additionally, it yielded these results in 19.62
±
0.99 seconds, slower than the PUMA and PHUN methods.A comparison of the (a) weightedL
0 measures and (b) execute time of the SDM
L
P, conventional minimum L
p-norm, PUMA, and PHUN methods in unwrapping the multislice dataset.
(a)
(b)


## 4. Conclusions

In this work, we developed a sparse-representation-based direct minimumL
p-norm (SDM
L
P) algorithm for the 2D phase unwrapping.The user-defined weights are introduced in the objective function to improve the discontinuity-preserving ability of theSDM
L
P 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
L
P algorithms computed effectively and efficiently by employing direct solvers.The proposed algorithm does produce excellent, reliable results with a very small weightedL
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
L
P method converges much faster than the conventional minimum L
p-norm method. However, it is not fast enough. Further research can be devoted to reducing the execution time.



---

*Source: 134058-2014-03-26.xml*