2.1 Fast phase unwrapping

The common observation in phase unwrapping methods is that the phase has to be corrected when an abrupt change occurs in the analyzed data (e.g., wrapping over a 2π range). One of the strategies used to correct the phase wraps is to calculate the phase derivatives and then compute the integral of the result. This leads to the estimation of the unwrapped phase, as will be explained in detail in a subsequent section. Depending on the data dimensionality, one, two, or three partial spatial derivatives have to be calculated. In the approach presented by Schofield et al. [43], a two-dimensional (2D) phase derivative calculation was performed using 2D Laplace operators. The key observation in this popular algorithm was that the use of Fourier transformations simplifies and speeds up the computations. Therefore, the name fast phase unwrapping (FPU) was used. Herein, we outline the FPU algorithm for data with an arbitrary number of spatial dimensions before we introduce the optimizations in the next section. The dataflow in the algorithm is depicted in Fig. 1.

Graphical representation of the data flow in the fast phase unwrapping (FPU) algorithm. The wrapped phase distribution, φ (r), is the input data. It is used to calculate the phase estimate ψ_est (r). The phase correcting function Q(r) is computed as the difference between the phase estimate and the input phase scaled by 2π and rounded to the nearest integer numbers, i.e., it is an integer map of 2π phase wraps. The phase wrap map is multiplied by 2π and added to the input phase to yield the output unwrapped phase distribution ψ (r).

If the spatial phase distribution φ (r) in the acquired data has phase wraps, a correcting function Q(r) needs to be identified to obtain the unwrapped phase distribution ψ (r):

where r denotes the spatial coordinates, and Q(r) is a function used to correct the wrapped phase using an integer number of 2π radians. The estimate of Q(r) can be calculated as,

Q′(r) is a function of real numbers (can be fractional, positive, and negative). Accordingly, ψ_est (r) is the unwrapped phase estimate which should to be identified using only wrapped phase information. This can be done by integrating the derivatives of the phase, and in practice, by applying the Laplace operator ∇² and its inverse ∇⁻²:

Schofield proposed to solve this equation for ∇²ψ_est (r) by defining a function P(r) such that

where j is the imaginary unit. This function has a key property in that it has the same values for wrapped and unwrapped phase, i.e., it is insensitive to phase wraps. The Laplacian of the P(r) function can be expressed as,

which can be used to find ∇²ψ_est (r):

In the above formula, the Laplacian of the unwrapped phase distribution is expressed with the use of function P(r) which can be obtained from the known (measured) wrapped phase distribution. In the original approach, Schofield et al. [43] used Eq. 4 and Euler’s formula to derive the Laplacian of the unwrapped phase,

Application of the inverse Laplace operator to both sides of this equation leads to an expression for the unwrapped phase estimate,

The Laplace operators of an arbitrary function g(r) can be calculated using Fourier transformations (FT) [54],

where K is Fourier space conjugate of the vector r, and F{} and F−1{}are forward and inverse Fourier transformation operators. With this observation, Eq. 8 can be rewritten as,

The above version of Schofield’s algorithm (Eq. (10)) requires eight Fourier transformations, and constitutes a computational challenge if the phase distribution is multidimensional. A straightforward improvement in the numerical efficiency arises from the simple rearrangement of Eq. (8), as proposed by Jeught et al. [37]:

which leads to a simplified version of Eq. (10) with six (instead of eight) Fourier transformations,

If the periodic boundary conditions required for the existence of the Fourier transforms are not fulfilled, phase unwrapping errors arise, thus preventing direct use of the estimate ψ_est (r) as the final unwrapped phase distribution. To overcome this issue, Schofield proposed the use of the phase correcting function defined as a difference between unwrapped and wrapped phase distributions, scaled by 2π, and rounded to the nearest integer values,

Q(r) can be thought of as an integer map of 2π phase wraps. The phase correcting function was subsequently used in Eq. (1) to calculate the unwrapped phase.