Open Access
Add to BibSonomy
Abstract
Phase unwrapping is a problem to reconstruct true phase values from modulo 2π phase values measured using various phase imaging techniques. This procedure is essentially formulated as a discrete optimization problem. However, most energy minimization methods using continuous optimization techniques have ignored the discrete nature and solved it as a continuous minimization problem directly, leading to losing exactness of the algorithms. We propose a new minimum norm method that can yield the optimal solution of the discrete problem by minimizing a continuous energy function. In contrast to the graph-cuts method, which is state of the art in this field, the proposed method requires much less memory space and a very simple implementation. Therefore, it can be simply extended to 3D or 4D phase unwrapping problems.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
A Phase imaging has been widely used in various phase imaging modalities such as interferometric synthetic aperture radar (InSAR), phase-contrast X-ray imaging, and the water and fat separation problem in MRI etc. [1–6]. In these applications, a phase value at each image pixel is calculated from some trigonometric function, and we can only measure its principal value wrapped into the range $[{ - \pi ,\pi } )$. Phase unwrapping is a problem to recover the true phase values from the wrapped ones.
In the past, many unwrapping methods have been developed. A simple method is to add (or minus) some integer multiple of $2\pi $ to the measured phase at each pixel such that gradients (differences between two neighboring pixels) of the unwrapped phase is less than $\pi $. This method has been called the path-following method in the literature. If the gradients of true phases are actually less than $\pi $, it is known that the path-following method yields an exact solution [7]. Improved variations of the path-following method have been developed such as branch cuts method and quality maps method [8]. In general, these path-following methods work very fast. However, it is known that they are weak to noise in the measured data leading to providing incorrect results [9].
The other unwrapping method is minimum norm (energy minimization) method, which is based on minimizing the gradients of unwrapped phases at each pixel (differences between two neighboring pixels) in the 4 directions, i.e. top, bottom, left, and right, summed up for all the pixels. The gradients are evaluated according to ${L^p}({p \ge 0} )$-norm [10–13]. Unlike the path-following methods, the minimum norm method is a path-independent algorithm and is robust to the noise.
Alternative phase unwrapping methods based on respective different principles have been also developed. For example, these include the methods which perform denoising and unwrapping simultaneously [14–20], a method using multiple interferograms [21], methods using deep learning [22,23], and a large-scale phase unwrapping method [24].
Note that, in the phase unwrapping, the difference between a wrapped phase and an unwrapped one at each pixel is necessarily some integer multiple of $2\pi $. Therefore, in the minimum norm method, it causes a discrete minimization problem. However, most methods using continuous optimization techniques ignore this discrete nature, from which they cannot guarantee exactness of the solution. Such a difficulty can be overcome by solving the discrete minimization problem directly. As a typical method, the graph-cuts method [25] using max-flow/min-cut technique was developed to obtain a more accurate result. More specifically, the graph-cuts method can exactly minimize the norm according to ${L^1}$-norm under the above-mentioned discrete constraint. Since the minimum norm method with ${L^p}(p < 1)$-norm becomes an intractable non-convex optimization, the graph-cuts method can be considered state of art in this field. However, the graph-cuts method requires a large memory space to construct a graph with image pixels as nodes and cliques as edges into the memory. Therefore, it cannot be easily extended to higher-dimensional problems in 3-D or 4-D imaging [26–28].
In this paper, we propose a new phase unwrapping method. Similar to the most minimum norm methods, the proposed method also solves a continuous minimization (with ${L^1}$-norm) but can obtain an exact solution of the unwrapping problem under the discrete constraint. In contrast to the graph-cuts method, the proposed method requires less memory space and its implementation is very simple. The proposed method is based on minimizing a continuous convex function by using the dual ascent method followed by binarizing the solution to $\{{0,2\pi } \}$. Therefore, it can be easily extended to the higher-dimensional unwrapping problems.
This paper is organized as follows. In section 2, we briefly review the minimum norm method and explain the proposed phase unwrapping method. In section 3, we show results of simulation studies as well as an experiment with real data of electron holography. The conclusion is given in section 4.
2. Theory
Throughout this paper, we denote image pixel by $({i,j} )$. In many applications of phase imaging, a phase at each pixel is obtained from trigonometric functions, which yields a wrapped phase as
Fig. 1. (a) Original phase image ${\phi _{i,j}}$, (b) Wrapped phase image ${\psi _{i,j}}$, (c), (d) plots of profiles corresponding to (a) and (b) along central horizontal lines, respectively.
2.1 Energy function
In the minimum norm method, we wish to find the unwrapped phases ${\phi _{i,j}}$ from the measurements ${\psi _{i,j}}$ with the requirement that the unwrapped phases minimize the cost function which include the differences (between two neighboring pixels) in the 4 directions, i.e. top, bottom, left, and right, summed up for all the pixels, where the difference is evaluated according to the ${L^1}$-norm.
In [25], it was pointed out that this problem can be solved by finding the binary image ${l_{i,j}}$ in Eq. (1), where a value ${l_{i,j}}$ at each pixel can only take 0 or 1. In the following, we assume that ${l_{i,j}} \in \{{0,1} \}$. The intuitive meaning of ${l_{i,j}}$ is the indicator value whether the wrapping occurred or not at each pixel, where ${l_{i,j}} = 0$ in the absence of wrapping and ${l_{i,j}} = 1$ in the presence of wrapping. Denote the differences in the phase values in horizontal and vertical directions by
respectively. Then, the phase unwrapping problem can be formulated as the following binary minimization problem [25].The cost function in Eq. (4) means that the phase unwrapping formulated by the minimum norm method is a discrete optimization problem. A few novel unwrapping methods such as SNAPHU [29], Flynn's minimum discontinuity algorithm [30] and PUMA [21] can treat the discrete nature by utilizing algorithms of discrete optimization such as network flow or graph-cuts. Among them, PUMA is attractive because it can obtain an exact solution of the energy minimization problem using the graph-cuts. However, these methods based on discrete optimization tend to require a complicated implementation. Furthermore, the methods using the graph-cuts such as PUMA require a large memory space due to constructing the graph into the memory. If we are able to exactly solve the unwrapping problem, i.e. minimizing Eq. (4) under the binary constraint, using techniques of continuous optimization, the issue of memory space and complexity of implementation will be overcome by eliminating the necessity of graph data structure. However, most methods using continuous optimization techniques ignored this binary constraint, i.e. ${l_{i,j}}$ must be 0 or 1, because it is hard to exactly treat this constraint [11,12,13]. Without this constraint, however, it cannot be guaranteed that the resulting algorithm yields an exact solution. Based on the above rationale, this paper aims at developing a new unwrapping method based on continuous optimization which guarantees to exactly minimize Eq. (4) under the binary constraint. We also extend it to the case where the wrap count at each pixel takes one of one of $\{{0,1, \cdots ,L} \}$ with $L > 1$, which corresponds to a more general unwrapping problem. Since the implementation of the proposed method is simple with less memory space, we expect that it helps in overcoming the issues of memory usage and complexity of implementation in the other exact methods such as the graph-cuts method.
As an alternative issue, some unwrapping methods [14–20] perform unwrapping and denoising simultaneously by artificially dropping the binary constraint. Such methods are beyond the scope of this paper. This paper only deals with solving the unwrapping problem under the discrete constraint by using a simpler method based on continuous optimization.
2.2 Continuous optimization
As mentioned above, Eq. (4) is a discrete optimization problem and is difficult to solve compared to standard convex optimization problems involving only continuous variables. In this paper, we use a novel technique to convert the discrete minimization into a continuous minimization without the binary constraint in an exact way. This technique has appeared in the optimization literature several times by different researchers [31,32]. For the case of our problem, the most important result is summarized as the following theorem.
Theorem 2.1 Let $\{{x_{i,j}^{\textrm{min}}} \}$ is the (unique) solution of continuous optimization problem
Remarks on the practical implication of this theorem are as follows. First, in Eq. (7), the variables to optimize ${x_{i,j}}$ are continuous (not binary) so that Eq. (7) can be solved easily by using the standard convex optimization methods with the bound constraint on the variables ${x_{i,j}} \in [{0,1} ]$. More specifically, the structure of Eq. (7) possesses a large similarity to the Total Variation (TV) minimization so that most of TV minimization methods can be used to solve it. Second, Eq. (8) is a simple thresholding of the solution obtained by solving Eq. (7) with almost no computational cost. Therefore, Theorem 2.1 essentially implies that the discrete optimization of Eq. (4) can be solved by the continuous optimization of Eq. (7).
We provide a detailed mathematical explanation of Theorem 2.1 in Appendix A. As pointed out in [25], $E_{i,j}^h,\; E_{i,j}^v$ are submodular functions. By Proposition 3.1 ([31]), we see that the energy function in (7) is a convex function meaning that Eq. (7) is a convex optimization. Theorem 2.1 guarantees that the thresholding operation of Eq. (8) yields the optimal solution of Eq. (4). In fact, in some of the previous (continuous) ${L^p}({p \ge 0} )$-norm methods, we can also use similar operation as in Eq. (8) to approximately obtain a discrete solution for the unwrapping. In general, however, the operation like Eq. (8) may result in losing optimality of the solution. With the above formulation, it is mathematically exact.
2.3 Coordinate descent method
Note that the cost function of Eq. (7) is non-differentiable so that we use the dual method to solve it. Even if the primal problem is non-differentiable, the dual problem becomes differentiable, which is easier to solve. The typical example is the famous total variation minimization [33].
In Appendix C, we show the derivation of dual problem corresponding to Eq. [7]. We only need to solve the following simpler minimization problem.
In this work, we use Coordinate Descent (CD) method to solve Eq. (9) for its simplicity and fast convergence. The idea of CD method is very simple. To minimize a cost function of Eq. (9) having many variables, CD method successively minimizes the function with respect to only a single variable ${p_{i,j}}$ (or ${q_{i,j}}$) while other variables are fixed at current values, where the selection of variables to update is performed according to the raster-scan order $i = 1,2, \cdots ,M;j = 1,2, \cdots ,N$ (the horizontal variables ${p_{i,j}}$ are updated first and the vertical variables ${q_{i,j}}$ are updated next). More specifically, the final algorithm is summarized in the following form with k as the main iteration number.
[Initialization] Set $k = 0\; $, set the variables as $(p_{i,j}^{(0 )},q_{i,j}^{(0 )}) = ({0,0} )$.
[Main Iteration]
For $k = 1,2, \cdots $
For $i = 2,3, \cdots ,M - 1;\; j = 2,3, \cdots ,N - 1$ perform the following.
Update ${p_{i,j}}$ as
For $i = 2,3, \cdots ,M - 1;\; j = 2,3, \cdots ,N - 1$ perform the following.
Update ${q_{i,j}}$ as
2.4 Method of unwrapping with more than one wrap counts
Until now, we described how to solve the unwrapping problem when the wrap count value ${l_{i,j}}$ is 0 or 1 for every pixel. In many applications, however, severe wrappings with more than one wrap counts can occur. In this case, the wrap count variable ${l_{i,j}}$ can take one of $\{{0,1, \cdots ,L} \}$ with $L > 1$. The proposed method can be easily extended to the multiple wrap count case by using the iteration process described below [25]. The idea is simple and is explained as follows. Starting with the initial wrap count values ${l_{i,j}}\, = 0$ for every pixel, we increment the ${l_{i,j}}$ values at each pixel one by one, by successively solving the binary unwrapping problem using the algorithm in section 2.3. The final algorithm is summarized in the following form.
[Initialization] Set $n = 0$, set the variables as $l_{i,j}^{(0 )} = 0$.
[Main Iteration]
For $n = 0,1, \cdots ,L - 1$, perform the following two steps.
The final solution of the multiple count case is given by $l_{i,j}^{(L )}$, which is the output of the above algorithm. In the fields of discrete optimization, algorithms which repeatedly solve a simple problem in such a way that the cost value always decreases is called greedy algorithm. The above algorithm can be seen as a particular application of the greedy algorithm. When the general ${L^p}$-norm with $p \ne 1$ is used in the cost function, the above algorithm does not guarantee that the obtained solution $l_{i,j}^{(L )}$ coincides with the true optimum solution in the range ${l_{i,j}} \in \{{0,1, \cdots ,L} \}$ with $L > 1$. However, for the case of ${L^1}$-norm as treated in this paper, it is known that the above algorithm yields the optimal solution [25].
One important problem in implementing the above algorithm is how to determine the maximum wrap count L needed to terminate the iteration. A few approaches can be considered for its solution. Among them, we have adopted the following simple approach. We implement the above algorithm by assuming that L is a sufficiently large number. We calculate the energy value after each outer iteration, i.e. loop of n, is finished. If the energy value is not decreased, we terminate the iteration and finish the unwrapping procedure, because the non-decrease of energy value implies that pixels having the wrap count n+1 do not exist. The value of loop count n+1 when the iteration is terminated yields the maximum wrap count L. In our implementation studies, this approach worked well in most cases, although there is no theoretical guarantee. In Fig. 2, we show what will be happen when the wrapping was performed with a wrong value of L.
Fig. 2. (a) Original phase image (128 × 128 pixels), (b)Wrapped image, (c) unwrapped image with maximum wrap count 5, (d) unwrapped image with maximum wrap count 7, (e) unwrapped image with correct maximum wrap count 8.
3. Experiments
3.1 Simulation studies
In this section, we show results of some simulation studies to demonstrate performances of the proposed method. Figure 3(a) shows an original phase image of 256 × 256 (pixels) constructed by a two-dimensional Gaussian function $({ - 1 \le x,y \le 1} )$ with maximum value 15, mean $\mu = 0$, and standard deviation $\sigma = 0.1$. Figure 3(b) shows the corresponding wrapped image. For this example, we assumed that the wrapping is very severe with a maximum value of wrap count ${l_{i,j}} = 2$ and with no noise in the simulation. We used the method with more than one wrap count described in section 2.4, which applies the basic method of section 2.3 multiple times. Figure 3(c) shows the unwrapped image by the proposed method, where the exact unwrapping was achieved with this simple example.
Fig. 3. Example of unwrapping using a numerically constructed phase image with Gaussian intensity change consisting of 256 × 256 (pixels). (a) Original image, (b) wrapped image, (c) unwrapped image by the proposed method.
We also performed simulation studies by using real images consisting of 300 × 300 (pixels) shown in Fig. 4(a) and (e) to demonstrate validity of the proposed method. The used original images consisting of 321 × 481 (pixels) are available at Berkeley segmentation dataset at the following URL [34].
Fig. 4. (a) Mountains, and (e) Ski as the original phase images. Unwrapped images by (b), (f) Goldstein's method, (c), (g) Graph-cuts method, and (d), (h) the proposed method, respectively.
https://www2.eecs.berkeley.edu/Research/Projects/CS/vision/bsds.
We converted gray levels of these test images into the range $[{0,\;2.3\pi } )$ and then wrapped into the range $[{ - \pi ,\pi } )$ by applying the wrapping operation in Eq. (1). We compared the proposed method with Goldstein's method (belonging to the path-following method) and the graph-cuts method, which are other standard approaches of unwrapping.
Figure 4. show results of the unwrapping for Mountains, and Ski images. The results show that the two images were unwrapped with sufficient accuracy by both the graph-cuts method and the proposed method. It can be observed that the images unwrapped by both the graph-cuts (Figs. 4(c), (g)) and the proposed method (Figs. 4(d), (h)) are almost exact and are identical to each other. However, the unwrapped image by Goldstein’s method (Figs. 4(b), (f)) suffered from considerable errors. This result matches to the theory in the following sense. The graph-cuts method and the proposed method are exact methods using the same cost function, i.e. ${L^1}$-norm with the discrete constraint, so that their results should be almost identical except for the effects of numerical errors in the implementations.
We also performed simulation studies with noisy image. In general, for unwrapping images with noise, it is hard to distinguish good images from bad images. Therefore, we compared the corresponding label (wrap count) images $\{{l_{i,j}} \}$ to evaluate the results. We compared the proposed method with the four other algorithms, which are Goldstein's method, weighted least-squares phase unwrapping with preconditioned conjugate gradient technique (PCG), partial differential equation phase unwrapping with ${L^2}$-norm (L2-PDE), and the graph-cuts method. Codes corresponding to the compared first three algorithms are available in [8]. We used an image of penguin in this simulation, which is also available at the Berkeley segmentation dataset [34]. The gray level of image was rescaled into the range $[{0,\;2.3\pi } ]$, and Gaussian noise with mean $\mu = 0$, and standard deviation $\sigma = 4.5$ was added to the image and then wrapped into the range $[{ - \pi ,\pi } )$. Figure 5(a) and (b) show the noisy input image and the corresponding true label image, respectively. We can consider Fig. 5(b) as the correct solution in this simulation. Figures 5(c), (d), (e), (f), and (g) show label images obtained by Goldstein's method, PCG, L2-PDE, graph-cuts method, and the proposed method, respectively. The corresponding Root Mean Squared Error (RMSE) value of phase error for each unwrapped image is shown in the bottom-right corner. These results show that the path-following method and the minimum norm methods with the ${L^2}$-norm cost, i.e. PCG and L2-PDE, are affected by the noise much more compared to graph-cuts method and the proposed method. In this simulation, we note that PDE method with the ${L^0}$-norm may give a better result than with the ${L^2}$-norm [11]. However, as mentioned above, it is a non-convex optimization, which is beyond a scope of this paper. In addition, as in the noise-free simulation case, graph-cuts method and the proposed method produced almost same label images, which means that the proposed method can produce an exact minimum solution of the ${L^1}$-norm cost function similarly to graph-cuts method.
Fig. 5. From left to right: (a) Simulated noisy Penguin image, (b) true label image, (c) label image by Goldstein's method, (d) label image by PCG, (e) label image by L2-PDE, (f) label image by graph-cuts method, (g) label image by the proposed method. The corresponding RMSE value of phase error for each unwrapped image is shown in the bottom-right corner.
3.2 Unwrapping for more challenging cases
We also performed simulation studies for more challenging cases. It is well-known that accurate unwrapping becomes difficult when 1) phase images consist of dense and complex fringes or 2) phase images contain parts with discontinuous phase values. Furthermore, we also performed experimental studies using real data distributed in [8] aiming at making it easier to compare the results shown in other literature. These results are summarized in this section.
(1) Case of phase images with dense and complex fringes
For this experiment, we constructed two wrapped images shown in Fig. 6(b) and Fig. 6(d) according to the following steps. We prepared a phase image with 128 × 128 (pixels) by 2-D Gaussian function $({ - 1 \le x,y \le 1} )$ with maximum value 50, mean $\mu = 0$, and standard deviation $\sigma = 0.12$ (Fig. 6(a)). We applied the wrapping operator to Fig. 6(a) to obtain the corresponding wrapped image called Example 1 here (Fig. 6(b)). Furthermore, we calculated the difference with respect to the horizontal direction followed by the difference with respect to the vertical direction to Fig. 6(a) and multiplied each pixel value by 1200 to obtain an alternative phase image (Fig. 6(c)). The new phase image included some discontinuous part in which the corresponding difference from adjacent pixel exceeds $\pi $. We applied the wrapping operator to Fig. 6(c) to obtain the corresponding wrapped image called Example 2 here (Fig. 6(d)). We applied the proposed method to Example 1 and Example 2. The results are shown in Figs. 6(e) and 6(f), and exact unwrapping was achieved for both Example 1 and Example 2.
Fig. 6. Simulation results of unwrapping for interferograms with dense and complex fringes. (a) Original phase image, (b) wrapped image corresponding to (a) (Example 1), (c) alternative original phase image constructed by taking the difference along the horizontal direction followed by taking the difference along the vertical direction of (a) (each pixel value was multiplied by 1200), (d) wrapped image corresponding to (c) (Example 2), (e) unwrapped image by the proposed method for Example 1, (f) unwrapped image by the proposed method for Example 2.
(2) Case of phase images containing parts with discontinuous phase values
We also performed simulation studies for a phase image containing parts with discontinuous phase values. For this experiment, we constructed a wrapped image shown in Fig. 7(b) according to the following steps. First, we prepared a phase image consisting of 128 × 128 (pixels) shown in Fig. 7(a). It has a uniform value $4\pi $ in the left part, and the right part is a 2-D Gaussian function same as in Fig. 6(a). Note that this image has a discontinuous part, i.e. boundary, between the left and the right. Then, we added Gaussian noise with standard deviation $\sigma = 3.5$ followed by applying the wrapping into the range $[{ - \pi ,\pi } )$ to obtain the corresponding wrapped image shown in Fig. 7(b).
Fig. 7. (a) Original noisy phase image constructed by replacing the left part of Fig. 6(a) by a uniform value, (b) wrapped image corresponding to (a).
Figures 8(a), (c), (e), (g), and (i) show unwrapped images by Goldstein's method, PCG, L2-PDE, the graph-cuts method, and the proposed method, respectively. These results show that the path-following method and the minimum norm methods with the ${L^2}$-norm cost, i.e. PCG and L2-PDE, are affected by the discontinuous part and noise much more compared to the graph-cuts method and the proposed method. We also note that the RMSE value in Fig. 8(a) is significantly larger than the other methods, mainly due to the existence of large errors in the left uniform part. We expect that the result of Fig. 8(a) can be improved by applying the boundary condition after the unwrapping. In fact, if evaluate the RMSE value only for the right Gaussian part of Fig. 7(a), the RMSE value of Goldstein's method, PCG, L2-PDE, graph-cuts method, and the proposed method changed to 2.60, 2.84, 4,33, 0.67, and 0.33, respectively. However, in any case, the graph-cuts method and the proposed method produced better results than the other methods.
(3) Experimental studies using real data distributed in [8]
Fig. 8. Unwrapped images for the wrapped image of Fig. 7(b) and the corresponding histograms of phase error. (a) Unwrapped image by Goldstein's method, (b) corresponding histogram of phase error, (c) unwrapped image by PCG, (d) corresponding histogram of phase error, (e) unwrapped image by L2-PDE, (f) corresponding histogram of phase error, (g) unwrapped image by the graph-cuts method, (h) corresponding histogram of phase error, (i) unwrapped image by the proposed method, (j) corresponding histogram of phase error. RMSE value of phase error for each unwrapped image is also shown in the bottom-right corner.
We also applied the proposed method to a real image distributed in [8]. Figure 9(a) shows a real wrapped phase map obtained by Interferometric Synthetic Aperture Radar (IFSAR). Figures 9(b), (c), (d), (e), and (f) show the corresponding unwrapped images by Goldstein's method, PCG, L2-PDE, the graph-cuts method, and the proposed method, respectively. As in many real data cases, we do not know the correct image for this case so that it is hard to say which image is better than the others. However, we can still perform some qualitative visual inspection for the accuracy. For example, in Fig. 9(d) unwrapped by L2-PDE, there exist additional horizontal curved lines. From wrapped images, we expect that such horizontal lines should not be emerged. In addition, in Fig. 9(b) unwrapped by Goldstein's method, additional dark pattern, which does not exist in the other methods, emerged in many parts. PCG method gave a rather smooth image (although not very accurate), mainly because it solves a continuous optimization problem instead of exact discrete optimization problem under discrete constraint. The graph-cuts method and the proposed method again provided very similar results.
Fig. 9. (a) Real wrapped interferometric SAR image distributed in [8], (b) unwrapped image by Goldstein's method, (c) unwrapped image by PCG, (d) unwrapped image by L2-PDE, (e) unwrapped image by the graph-cuts method, and (f) unwrapped image by the proposed method.
3.3 Application to real image of electron holography
We applied the proposed method to two interferograms measured by a 300 kV transmission electron microscope. The two wrapped images are shown in Fig. 10(a1) and 10(a2). The magnification ratio in this measurement was approximately 270,000 times, and the image consisted of 800 × 900 (pixels). Furthermore, the image of Fig. 10(a2) was measured in shorter time than the image of Fig. 10(a1), so that Fig. 10(a2) contains more noise compared to Fig. 10(a1). We show results of unwrapping with various methods in Figs. 10(b1)-(f1) and Figs. 10(b2)-(f2). From the results, it can be observed that Goldstein's method is severely affected by the noise compared to the other methods. Furthermore, compared to the two ${L^2}$-norm methods (PCG and L2-PDE), the two ${L^1}$-norm methods (graph-cuts method and the proposed method) gave better unwrapped images which preserve fine structures. We believe that such a difference originates from the use of different cost function (${L^2}$ vs. ${L^1}$) together with the exact treatment of discrete nature of the unwrapping. Finally, we strengthen that graph-cuts method and the proposed method again gave almost same results, and the proposed method worked reasonably with the real data.
Fig. 10. Unwrapped results for real data of a Latex sphere measured by electron holography. Top row shows the case with lower noise, and bottom row shows the case with higher noise. From left to right: (a1), (a2) show wrapped images, (b1), (b2): unwrapped images by Goldstein's method, (c1), (c2): unwrapped images by PCG, (d1), (d2): unwrapped images by L2-PDE, (e1), (e2): unwrapped images by graph-cuts method, (f1), (f2): unwrapped images by the proposed method. The gray scale to show unwrapped images was set to 0.0-60.0.
3.4 Memory usage and convergence speed
In [9], it was pointed out that network flow method possesses a drawback of intensive memory usage and its large computational complexity. The graph-cuts method also possesses similar problems. Here, we make a discussion on the comparison of memory usage between the graph-cuts method and the proposed method. In implementing the graph-cuts method, a graph with each image pixel as node and each pair of neighboring pixels as edge needs to be constructed into the memory. We note that the construction of this graph cannot be avoided. On the other hand, the proposed method can be implemented using the array data structure representing the image. This difference of used data structure in implementation can produce a large difference in the memory usage, i.e. the proposed method requires less memory space.
Below, we show an example of memory usage in our implementation. Of course, the necessary memory space in each method can change dependent on details of implementation, programs, and use of special implementation techniques. So, we note that the comparison below is just an example in our case. In our implementation of graph-cuts method, we only converted MATLAB program (
http://www.lx.it.pt/∼bioucas/code.htm) into c++ program. Table 1 shows a comparison of necessary memory space between the graph-cuts method and the proposed method for three typical image sizes. It can be observed that the proposed method remarkably decreased the memory usage of graph-cuts method. Although this is just an example in our case, we believe that the memory usage of graph-cuts method cannot be smaller than that of the proposed method as long as the implementation of graph-cut method requires to construct the graph into the memory.
Table 1. Memory usage (Mb) for images with different sizes
In phase imaging applications, in higher-dimensions such as 3-D MRI and 4-D imaging, or in InSAR using high resolution satellite imagery, it is obvious that the memory usage will increase rapidly [35,36]. Therefore, to decrease the memory usage is an important issue. In addition, in higher dimensional problems, the graph used in graph-cuts method will become rather complex. However, the proposed method is expected to be still simple even in the higher dimensional problems thanks to its simple structure together with its continuous optimization structure.
We also investigated convergence speed of the proposed method using Ski image shown in Fig. 4(e). We converted gray levels of ski image into the range $[{0,\;4.6\pi } )$ and then wrapped into the range $[{ - \pi ,\pi } )$. Many minimum norm unwrapping methods involve an iteration procedure. The convergence speed is determined by the product of computation cost per iteration and iteration number. Figures 11(a1), 11(a2), and 11(a3) show unwrapped results when the number of iterations was only 1 by PCG, L2-PDE and the proposed method, respectively. Figures 11(b1), 11(b2), and 11(b3) show error images between the original Ski image and the unwrapped images. Table 2 show RMSE values for each unwrapped image with different number of iterations. The results show that all the investigated minimum norm methods converge quickly. However, the accuracy is rather different as follows. The error image by PCG shown in Fig. 11(b1) is not taking discrete values (see also the corresponding histogram of error image in Fig. 11(c1)), which is caused by solving the unwrapping problem without the discrete constrain unlike the proposed method. Due to this drawback, the RMSE value of phase error in PCG shown in Table 2 becomes larger. On the other hand, the RMSE value of phase error by the proposed method converged very quickly and stably than that of L2-PDE. Finally, we remark that the fast convergence by the proposed method is thanks to 1) the use of discrete constraint in an exact way and 2) the fast convergence property of CD method known in various continuous optimization problems.
Fig. 11. Unwrapped results (gray scale 0-16) together with the corresponding phase error images (gray scale -7-13) and histograms of phase error images at 1 iteration of PCG, L2PDE, and the proposed method. (a1), (b1), (c1): PCG, (a2), (b2), (c2): L2-PDE (a3), (b3), (c3): proposed method.
Table 2. RMSE values by PCG, L2-PDE, and the proposed method at 0, 1, 5, and 10 iterations.
4. Conclusion
In this paper, we proposed a new phase unwrapping method which belongs to the minimum norm method. We gave a complete mathematical derivation to prove that the proposed method yields a global optimal solution of the unwrapping problem with the ${L^1}$-norm cost function under the discrete constraint. We note that such an exact algorithm with the simple continuous optimization is still rare in the unwrapping problem. The simulation results showed that the proposed method and graph-cuts method provide almost identical results within a small difference arising from numerical errors. The reason is that they are exact solution methods minimizing the same cost function under the discrete constraint. However, the proposed method possesses a significant advantage in memory usage compared to graph-cuts method. This advantage would be very useful in 3-D or 4-D problems or large size problems in future. Furthermore, the implementation of proposed method is very simple compared to graph-cut method.
Appendix A:
This part is largely indebted to [31]. We first define a discrete Total Variation (TV) as a convex nonnegative function $J:{\mathrm{\mathbb{R}}^{M{\kern 1pt} N}} \to [{0, + \infty } ]$ satisfying a discrete co-area formula
where ${\chi ^{u \ge z}}$ is a vector such that $\chi _i^{u \ge z} = 1({i = 1,2, \cdots } )$ if ${u_i} \ge z$, and $\chi _i^{u \ge z} = 0$ if ${u_i} < z$.Proposition 3.3 ([31]): Let $J$ is a discrete TV, $X = \{{{x_{i,j}}} \}\in {\mathrm{\mathbb{R}}^{M{\kern 1pt} N}}$ is the (unique) solution of
Then, for all $z > 0$, the characteristic functions of the super level sets ${U_z} = \{{X \ge z} \}$ and $U_z^{\prime} = \{{X > z} \}$ (which are different only if $z \in \{{{x_{i,j}},i = 1,2, \cdots ,M,\,j = 1,2, \cdots ,N} \}$) are the largest and smallest minimizers of
Next, we use this proposition to solve the unwrapping problem. Define a discrete function as
Lemma A1: Let $\widetilde {J\; }({{x_{i,j}}} )$ is a discrete TV of $\widetilde {J\; }({{l_{i,j}}} )$, and $X \in {\mathrm{\mathbb{R}}^{M\,N}}$ is the (unique) solution of
Then, for any constant $0 < z < 1/({4M{\kern 1pt} N} )$, the characteristic functions of both super level sets ${U_z} = \{{X \ge z} \}$ and $U_z^{\prime} = \{{X > z} \}$ are the minimizers of
Proof: Assume that there exist a constant $0 < z < 1/({4M{\kern 1pt} N} )$ and level sets $U_z^{\prime}$ (or ${U_z}$) are not the minimum point of $\tilde{J}({{l_{i,j}}} )$. Denote the minimum point of $\tilde{J}({{l_{i,j}}} )$ by $l_{i,j}^{\min }$. Since the variable in (20) changes discretely (it can take only 1or 0), we have
On the other hand, by the definition of X and Proposition 3.3, we have
Therefore, for any constant $0 < z < 1/({4M{\kern 1pt} N} )$, we have
Equation (23) is a contradiction to Eq. (21).
QED.
Theorem 2.1 is an immediate consequence of this Lemma.
Appendix B:
In the following, we derive the specific form of $E_{i,j}^h({E_{i,j}^v} )({{x_{i,j}}} )$. For simplicity, we assume that $E_{i,j}^h({0,0} )= 0,E_{i,j}^h({1,1} )= 0$, which can be realized by replacing $E_{i,j}^h$ by $E_{i,j}^h - |{a_{i,j}^h} |$. Note that, after this operation, the minimum point of Eq. (4) is not changed. It is obvious that $a_{i,j}^h$ can only take one of three values, i.e. 0 or 1 or -1. Therefore, pixels comprising the image can be divided into three classes.
If $a_{i,j}^h = 0$, by Eq. (5), we see that $E_{i,j}^h({1,0} )= E_{i,j}^h({0,1} )= 1$.
Then, for ${x_{i,j}} < {x_{i - 1,j}}$
Similarly, for ${x_{i,j}} \ge {x_{i - 1,j}}$, we have $E_{i,j}^h = {\Delta ^h}{x_{i,j}}$. Therefore,
We can also calculate the same quantity as Eq. (25) for $a_{i,j}^h = 1$ as $E_{i,j}^h ={-} {\Delta ^h}{x_{i,j}}$, and for $a_{i,j}^h ={-} 1$ as $E_{i,j}^h = {\Delta ^h}{x_{i,j}}$, respectively.For the vertical terms, we can also obtain the specific form of $E^{v}_{i,j}$ according to the similar procedure. Finally, Eq. (7) can be rewritten as follows.
Appendix C:
Note that the energy function of Eq. (26) is non-differentiable so that standard differentiable convex optimization methods cannot be used directly. In the following, using the method of Lagrange multipliers (Lagrangian duality), we transform Eq. (26) into a dual problem, which is a minimization problem of differentiable cost function. For the sake of simplicity, we consider only the horizontal terms in Eq. (26) to show the procedure to derive the dual cost function, but we can do the same for the vertical terms. First, we rewrite Eq. (26) (with only horizontal terms) to an equivalent form as
By Karush-Kuhn-Tucker conditions, the optimality condition is expressed as
Substituting Eq. (31) and Eq. (32) into Eq. (30), we see that the dual problem takes the following simple form.
Combining $\lambda _i^1,\lambda _i^2 \ge 0$ with Eq. (32) yields the constraints for ${p_i}$ as follows.
This derivation process can be easily extended to the general case where both the horizontal terms and the vertical terms exist, from which we obtain the dual problem given by Eq. (9).
Funding
JST-ERATO project of Japan (Grant Number JPMJER1403); JST-CREST project of Japan (Grant Number JPMJCR1765).
Acknowledgments
Real IFSAR image used in the experimental studies of Section 3.2 were taken from the book (Ref. [8]).
We appreciate Prof. Yasukazu Murakami (Kyushu University) and Prof. Takehiro Tamaoka (Kyushu University) for providing real data of electron holography used in this work. We appreciate Dr. Katsuya Fujii (University of Tsukuba) and Mr. Keita Wada for useful discussions and helping with the experiments.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. D. Massonnet and K. L. Feigl, “Radar interferometry and its application to changes in the earth’s surface,” Rev. Geophys. 36(4), 441–500 (1998). [CrossRef]
2. P. A. Rosen, S. Hensley, I. R. Joughin, F. K. Li, S. N. Madsen, E. Rodriguez, and R. Goldstein, “Synthetic aperture radar interferometry,” Proc. IEEE. 88(3), 333–382 (2000). [CrossRef]
3. F. Cigna and D. Tapete, “Present-day land subsidence rates, surface faulting hazard and risk in Mexico City with 112161–2020 Sentinel-1 IW InSAR,” Remote Sens. Environ. 253, 112161 (2021). [CrossRef]
4. F. Pfeiffer, T. Weitkamp, B. Oliver, and D. Christian, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nat. Phys. 2(4), 258–261 (2006). [CrossRef]
5. T. Weitkamp, A. Diaz, C. David, D. Christian, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “Quantitative X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). [CrossRef]
6. P. Jezzard and R. S. Balaban, “Correction for geometric distortion in echo planar images from B0 field variations,” Magn. Reson. Med. 34(1), 65–73 (1995). [CrossRef]
7. K. Itoh, “Analysis of the phase unwrapping problem,” Appl. Opt. 21(14), 2470 (1982). [CrossRef]
8. D. Ghiglia and M. Pritt, “Two-Dimensional Phase Unwrapping. Theory, Algorithms, and Software,” (Wiley, New York), (1998).
9. L. Ying, “Phase unwrapping,” In: M. Akay, ed. Wiley encyclopedia of biomedical engineering. 1st Ed. (John Wiley & Sons; Hoboken, NJ, 2006).
10. D. Fried, “Least-squares fitting a wave-front distortion estimate to an array of phase-difference measurements,” J. Opt. Soc. Am. 67(3), 370–375 (1977). [CrossRef]
11. D. Ghiglia and L. A. Romero, “Minimum Lp-norm two-dimensional phase unwrapping,” J. Opt. Soc. Am. A 13(10), 1999–2013 (1996). [CrossRef]
12. HYH. Huang, L. Tian, Z. Zhang, Y. Liu, Z. Chen, and G. Barbastathis, “Path-independent phase unwrapping using phase gradient and total-variation (TV) denoising,” Opt. Express 20(13), 14075–14089 (2012). [CrossRef]
13. U. S. Kamilov, I. N. Papadopoulos, M. H. Shoreh, D. Psaltis, and M. Unser, “Isotropic inverse-problem approach for two-dimensional phase unwrapping,” J. Opt. Soc. Am. A 32(6), 1092–1100 (2015). [CrossRef]
14. O. Loffeld, H. Nies, S. Knedlik, and W. Yu, “Phase unwrapping for SAR interferometry: A data fusion approach by Kalman filtering,” IEEE Trans. Geosci. Remote Sens. 46(1), 47–58 (2008). [CrossRef]
15. X. Xie and Y. Li, “Enhanced phase unwrapping algorithm based on unscented Kalman filter, enhanced phase gradient estimator, and path-following strategy,” Appl. Opt. 53(18), 4049–4060 (2014). [CrossRef]
16. J. J. Martinez-Espla, T. Martinez-Marin, and J. M. Lopez-Sanchez, “A particle filter approach for InSAR phase filtering and unwrapping,” IEEE Trans. Geosci. Remote Sens. 47(4), 1197–1211 (2009). [CrossRef]
17. X. Xie, “Iterated unscented Kalman filter for phase unwrapping of interferometric fringes,” Opt. Express 24(17), 18872–18897 (2016). [CrossRef]
18. C. Magnard, M. Frioud, D. Small, T. Brehm, and E. Meier, “Analysis of a maximum likelihood phase estimation method for airborne multibaseline SAR interferometry,” IEEE J. Sel. Top. Appl. Earth Observations Remote Sensing 9(3), 1072–1085 (2016). [CrossRef]
19. H. Hongxing, J. M. Bioucas-Dias, and V. Katkovnik, “Interferometric phase image estimation via sparse coding in the complex domain,” IEEE Trans. Geosci. Remote Sens. 53(5), 2587–2602 (2015). [CrossRef]
20. J. Pineda, J. Bacca, J. Meza, L. A. Romero, H. Arguello, and A. G. Marrugo, “SPUD: simultaneous phase unwrapping and denoising algorithm for phase imaging,” Appl. Opt. 59(13), D81–D88 (2020). [CrossRef]
21. D. C. Ghiglia and D. E. Wahl, “Interferometric synthetic aperture radar terrain elevation mapping from multiple observations,” In Proceedings of IEEE 6th Digital Signal Processing Workshop (pp. 33–36). IEEE. (1994).
22. G. E. Spoorthi, R. K. S. S. Gorthi, and S. Gorthi, “PhaseNet 2.0: Phase unwrapping of noisy data based on deep learning approach,” IEEE Trans. on Image Process. 29, 4862–4872 (2020). [CrossRef]
23. F. Yang, T. A. Pham, N. Brandenberg, M. P. Lütolf, J. Ma, and M. Unser, “Robust phase unwrapping via deep image prior for quantitative phase imaging,” IEEE Trans. on Image Process. 30, 7025–7037 (2021). [CrossRef]
24. C. W. Chen and H. A. Zebker, “Phase unwrapping for large SAR interferograms: Statistical segmentation and generalized network models,” IEEE Trans. Geosci. Remote Sensing 40(8), 1709–1719 (2002). [CrossRef]
25. J. M. Bioucas-Dias and G. Valadão, “Phase unwrapping via graph cuts,” IEEE Trans. Image Process. 16(3), 698–709 (2007). [CrossRef]
26. M. Loecher, E Schrauben, KM Johnson, and O. Wieben, “Phase unwrapping in 4D MR Flow with a 4D single-step Laolacian Algorithm,” J. Magn. Reson. Imaging 43(4), 833–842 (2016). [CrossRef]
27. NDF. Campbell, G. Vogiatzis, C. Hernández, and R. Cipolla, “Automatic 3D object segmentation in multiple views using volumetric graph-cuts,” Image and Vision Computing. 28(1), 14–25 (2010). [CrossRef]
28. J. Dong, F. Chen, D. Zhou, T. Liu, Z. Yu, and Y. Wang, “Phase unwrapping with graph cuts optimization and dual decomposition acceleration for 3D high-resolution MRI data,” Magn Reson Med. 77(3), 1353–1358 (2017). [CrossRef]
29. C. W. Chen and H. A. Zebker, “Two-dimensional phase unwrapping with use of statistical models for cost functions in nonlinear optimization,” J. Opt. Soc. Am. A 18(2), 338–351 (2001). [CrossRef]
30. T. J. Flynn, “Two-dimensional phase unwrapping with minimum weighted discontinuity,” J. Opt. Soc. Am. A 14(10), 2692–2701 (1997). [CrossRef]
31. A. Chambolle and J. Darbon, “On Total Variation Minimization and Surface Evolution Using Parametric Maximum Flows,” Int J Comput Vis. 84(3), 288–307 (2009). [CrossRef]
32. L. Lovász, “Submodular functions and convexity,” In Mathematical programming: The state of the art, A. Bachem, M. Grötschel, and B. Korte, eds., (Springer-Verlag, Berlin), 235–257 (1983).
33. A. Chambolle, “An algorithm for total variation minimization and applications,” Journal of Mathematical imaging and vision. 20(1), 89–97 (2004). [CrossRef]
34. D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” In Proc. 8th Int. Conf. Comput. Vis. 2, 416–423 (2001).
35. H. Yu, Y. Lan, J. Xu, D. An, and H. Lee, “Large-scale L°-norm and L1-norm 2-D phase unwrapping,” IEEE Trans. Geosci. Remote Sens. 55(8), 4712–4728 (2017). [CrossRef]
36. H. Yu, Y. Lan, Z. Yuan, J. Xu, and H. Lee, “A Review on Phase Unwrapping in InSAR Signal Processing,” IEEE Geosci. and Remote Sens. Mag. 7(1), 40–58 (2019). [CrossRef]
