On the performance of preconditioned methods to solve L^p-norm phase unwrapping

by   Ricardo Legarda-Saenz, et al.

In this paper, we analyze and evaluate suitable preconditioning techniques to improve the performance of the L^p-norm phase unwrapping method. We consider five preconditioning techniques commonly found in the literature, and analyze their performance with different sizes of wrapped-phase maps. Keywords.- Phase unwrapping, L^p-norm based method, Preconditioning techniques.


page 11

page 14


Estimating Absolute-Phase Maps Using ESPIRiT and Virtual Conjugate Coils

Purpose: To develop an ESPIRiT-based method to estimate coil sensitiviti...

A Multiphase Image Segmentation Based on Fuzzy Membership Functions and L1-norm Fidelity

In this paper, we propose a variational multiphase image segmentation mo...

Optimal convex lifted sparse phase retrieval and PCA with an atomic matrix norm regularizer

We present novel analysis and algorithms for solving sparse phase retrie...

2D-Phase Unwrapping via Balanced Spanning Forests

Phase unwrapping is the process of recovering a continuous phase signal ...

High-fidelity quantitative differential phase contrast deconvolution using dark-field sparse prior

Differential phase contrast (DPC) imaging plays an important role in the...

A new family of nonconforming elements with H(curl)-continuity for the three-dimensional quad-curl problem

We propose and analyze a new family of nonconforming finite elements for...

Optimization of Geometric Constellation Shaping for Wiener Phase Noise Channels with Varying Channel Parameters

We present a novel method to investigate the effects of varying channel ...

1 Introduction

There has been a growing interest in the development of new techniques for the processing of coherent signals. This kind of signals is generated by measurement techniques like synthetic aperture radar (SAR), magnetic resonance imaging (MRI), and interferometry among others. The objective of this processing is to estimate the phase term

, also known as phase map, from a signal which general model could be expressed as

where is the position and is the signal amplitude. The estimation of the phase term becomes very relevant, given that this term can be related to different physical quantities such as geographical topography in the case of SAR, or the optical path difference in the case of optical interferometry 6, 19.

However, as it can be seen from the signal model, the estimation of the phase term is not straightforward. Instead, the estimated phase map from the signal is defined as


where is a function that bounds the values to The term is known as wrapped phase and is a nonlinear function of This term is not useful for measurements because just offers the principal values of the phase term so it is necessary estimate from this wrapped phase . This process is called phase unwrapping 8.

Phase unwrapping is an ill-posed problem 1, 20. The unwrapping process consists of integrating the gradient field of the wrapped phase map 8. Even in ideal conditions, phase estimation is not trivial due to the non linearity of caused by cyclic discontinuities. In real conditions, the unwrapping process becomes very difficult, where noise signal, sub-sampling or differences larger than (real or not) generate ambiguities hard to process not allowing the accurate recovering of the phase map

In the literature there are two main strategies to solve the unwrapping problem 8, 7

: the first one, known as path-following or local methods, consists of integrating the differences of the wrapped phase over a path that covers the entire phase map. The second one, considers the problem globally and the solution is presented in the form of an integral over the wrapped region. In the global strategy, there are two approaches: the first consists of using Green functions together with a numerical solution based on the Fourier transform 

4, 5, 14, 13.

The second approach expresses the solution as a -norm minimization problem, resulting on the solution of weighted differential equations. The numerical solution of these differential equations leads to a nonlinear system which has to be solved iteratively with great computational cost 9, 12, 11. Typically, conjugate gradient (CG) or multigrid method are the methods of choice for this kind of problems. In the case of the conjugate gradient, matrix is expected to be well conditioned, otherwise the convergence will be slow. However, very often one comes across with ill conditioned matrices when working with these nonlinear systems. The solution is to precondition the matrix ; this is, instead of solving the original system we solve the preconditioned system The matrix , called a preconditioner for the matrix , is chosen to improve the condition number of the matrix In most cases, this preconditioning matrix is problem dependent.

The goal of this paper is to analyze and evaluate suitable preconditioning techniques to improve the performance of the -norm phase unwrapping method 9, 8. We consider five preconditioning techniques commonly found in the literature, and analyze their performance with different sizes of wrapped-phase maps. The organization of this paper is as follows: first, we describe the -norm phase unwrapping method and their numerical solution. Then, the performance of the selected preconditioning techniques is evaluated by numerical experiments with different sizes of a synthetic wrapped-phase map. Finally, we discuss our results and present some concluding remarks.

2 Methodology

The -norm based method proposed by D. C. Ghiglia and M.D. Pritt for 2-D phase unwrapping is defined as 9


where and is the domain of integration. To obtain the solution of the problem expressed in Eq. (2

), the first-order optimality condition or Euler-Lagrange equation has to be derived, resulting in the following partial differential equation (PDE)


with boundary conditions


denotes the unit outer normal vector to the boundary.

2.1 Numerical solution

Let to denote the value of a function at point defined on where the sampling points are with and are the number of points in the discrete grid of points. We use to represent any of the variables and defined in the previous equations. Derivatives are approximated using standard forward and backward finite difference schemes

The gradient and the divergence are approximated as


Hence the numerical approximation of the Eq. (3) is given by


with boundary conditions

where the terms y are defined as

where is the wrapping operator 8. Applying the previous discrete approximations, we have that the numerical solution of Eq. (3) is given by


with boundary conditions

where for simplicity and without loss of generality we consider The terms y are defined as


to force the values of the terms to lie in the range (0,1). This helps the stability and convergence of the numerical solution 2, 9. Usually, is selected for phase unwrapping process 9, 8.

Discretization of Eq. (3) leads to a nonlinear PDE because the terms and are functions of the input data and the solution. This is solved using the following iterative procedure 2, 17, 9: first, given an initial value the terms and are computed; then, the terms and are held fixed and Eq. (5) is solved using preconditioned conjugate gradient 18, 10. With the current solution the terms and are updated and a new solution is computed. This process is repeated until convergence.

Now, we present the algorithm to solve the discretization of Eq. (3). First we arrange Eq. (5) in matrix form as


Notice that is a sparse matrix which depends on the terms and so it needs to be constructed at each iteration. The preconditioned conjugate gradient (PCG) used in our work is the implementation proposed in Ref. 18, and the explicit structure of the algorithm is given in Algorithm 1.

Data: the wrapped phase and
Result: the unwrapped phase
1 while  do
2       compute and , Eq. (6) solve Eq. (7) using PCG: begin
3             construct and estimate preconditioning matrix from while   do
4                   if  is divisible by  then
6                  else
8                   end if
10             end while
12       end
15 end while
Algorithm 1 -norm phase unwrapping algorithm.

The key point of the performance of the algorithm is the proper selection of the preconditioning matrix , step 1 of Algorithm 1. can be defined in many different ways but it should meet the following requirements: a) the preconditioned system should be easy to solve, and b) the preconditioning matrix should be computationally cheap to construct and apply 16.

Here we test suitable preconditioning techniques to estimate matrix in Algorithm 1 and analyze the performance of the -norm phase unwrapping method. For this purpose, we selected the following commonly found preconditioning techniques in literature 21, 16, 3:

  1. where

    is the identity matrix; with this matrix, PCG becomes CG 


  2. where is the diagonal of matrix this technique is known as Jacobi preconditioning.

  3. the incomplete LU factorization with no fill-in.

  4. the incomplete Cholesky factorization with no fill-in.

  5. the successive over-relaxation factorization.

The main reason for selecting these preconditioning techniques was their no fill-in property. That is, the preconditioning matrices preserve the sparse structure of matrix and require just the same memory space used by matrix

3 Numerical experiments

To illustrate the performance of the selected preconditioning techniques, we carried out the numerical experiments using a Intel® Core i7 @ 2.40 GHz laptop with Debian GNU/Linux 10 (buster) 64-bit and 16 GB of memory. For our experiments, we programmed all the functions using C language, GNU g++ 8.3 compiler and Intel® MKL 2019 library. It is important to highlight that all the functions were programmed from scratch and were programmed from the basic algorithms, without using modifications that improve their performance.

We use the wrapped phase map shown in Fig. 1 as the data term of the Algorithm 1. This wrapped phase map has a resolution of 640x480 pixels, and this resolution was used as our reference scale.

We generated several scaled wrapped phase map from the one shown in Fig. 1, and used them as data in Algorithm 1. The scaled image sizes and the resultant sizes of sparse matrix in our experiments are shown in Table 1.

scale image size variables density (%)
0.25 160 x 120 95440 0.0250
0.50 320 x 240 382880 0.0064
0.75 480 x 360 862320 0.0028
1.00 640 x 480 1533760 0.0016
1.25 800 x 600 2397200 0.0010
1.50 960 x 720 3452640 0.0007
1.75 1120 x 840 4700080 0.0005
2.00 1280 x 960 6139520 0.0004
Table 1: Image sizes used in the experiments.

For each scaled wrapped phase map, we tested the selected five preconditioning matrices. The stopping criteria used in Algorithm 1 were Fig. 2 shows the iterations needed to obtain the solution for the different scaled wrapped phase maps, and Fig. 3 shows the computational time consumed for each preconditioning technique. It is worth the value to remark that the time taken to construct preconditioning matrix : except for the incomplete Cholesky factorization, all the preconditioning techniques employed between 1% and 2% of the total processing time. However, the incomplete Cholesky factorization took more than 90%. This is why in Figs 2 and 3, we only show the first four experiments with this technique.

Finally, we use a normalized error to compare the unwrapping estimation; this error is defined as 15:


where and are the signals to be compared. The normalized error values vary between zero (for perfect agreement) and one (for perfect disagreement). For all the cases, we found that the normalized error was around . An example of the resultant unwrapped phase map is shown in Fig. 4.

4 Discussion of results and conclusions

From the results obtained, we have the following remarks. First, any of the techniques used had no impact on the obtained results, since for all the experiments the normalized error Q was approximately the same, . The differences are evident when we analyze the performance in terms of iterations and computational time. If we only consider the number of iterations, the incomplete LU factorization and incomplete Cholesky factorization clearly show their advantage over the other three methods. However, this perception changes when we also consider the computational time used. Clearly, the incomplete Cholesky factorization consumed a lot of time which makes it unviable to be considered as a preconditioning technique. In general, numerical results show the incomplete LU factorization as the best choice to be used as preconditioning technique in Algorithm 1.

As work in the future, we are going to analyze the use of dedicated libraries for estimating the preconditioning matrix, where we hope to obtain better results.


  • 1 Bertero, M., Boccacci, P.: Introduction to Inverse Problems in Imaging. Institute of Physics Publishing, Bristol (1998)
  • 2 Bloomfield, P., Steiger, W.: Least Absolute Deviations: Theory, Applications and Algorithms. Birkhäuser Boston (1983)
  • 3 Chen, K.: Matrix Preconditioning Techniques and Applications. Cambridge University Press, Cambridge (2005)
  • 4 Fornaro, G., Franceschetti, G., Lanari, R.: Interferometric SAR phase unwrapping using Green’s formulation. IEEE Transactions on Geoscience and Remote Sensing 34(3), 720–727 (1996). https://doi.org/10.1109/36.499751
  • 5 Fornaro, G., Franceschetti, G., Lanari, R., Sansosti, E.: Robust phase-unwrapping techniques: a comparison. Journal of the Optical Society of America A 13(12), 2355–2366 (1996). http://dx.doi.org/10.1364/JOSAA.13.002355
  • 6 Gasvik., K.J.: Optical metrology. John Wiley & Sons, Chichester, 3rd edn. (2002)
  • 7 Gens, R.: Two-dimensional phase unwrapping for radar interferometry: Developments and new challenges. International Journal of Remote Sensing 24(4), 703–710 (2003). https://doi.org/10.1080/0143116021000016725
  • 8 Ghiglia, D.C., Pritt, M.D.: Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software. Wiley-Interscience, New York (1998)
  • 9 Ghiglia, D.C., Romero, L.A.: Minimum L p-norm two-dimensional phase unwrapping. Journal of the Optical Society of America A 13(10), 1999–2013 (1996). http://dx.doi.org/10.1364/JOSAA.13.001999
  • 10 Golub, G.H., Van Loan, C.F.: Matrix Computations. Johns Hopkins University Press, Baltimore, third edn. (1996)
  • 11 Guo, Y., Chen, X., Zhang, T.: Robust phase unwrapping algorithm based on least squares. Optics and Lasers in Engineering 63, 25–29 (2014). http://dx.doi.org/10.1016/j.optlaseng.2014.06.007
  • 12 Hooper, A., Zebker, H.A.: Phase unwrapping in three dimensions with application to InSAR time series. Journal of the Optical Society of America A 24(9), 2737–2747 (2007). http://dx.doi.org/10.1364/JOSAA.24.002737
  • 13 Lyuboshenko, I.: Unwrapping circular interferograms. Applied Optics 39(26), 4817–4825 (2000). http://dx.doi.org/10.1364/AO.39.004817
  • 14

    Lyuboshenko, I., MaiTre, H.: Phase unwrapping for interferometric synthetic aperture radar by use of Helmholtz equation eigenfunctions and the first Green’s identity. Journal of the Optical Society of America A

    16(2),  378 (1999). https://doi.org/10.1364/JOSAA.16.000378
  • 15 Perlin, M., Bustamante, M.D.: A robust quantitative comparison criterion of two signals based on the Sobolev norm of their difference. Journal of Engineering Mathematics 101(1), 115–124 (dec 2016). http://dx.doi.org/10.1007/s10665-016-9849-7
  • 16 Saad, Y.: Iterative methods for sparse linear systems. SIAM, second edn. (2003)
  • 17 Scales, J.A., Gersztenkorn, A.: Robust methods in inverse theory. Inverse Problems 4(4), 1071–1091 (oct 1988). https://doi.org/10.1088/0266-5611/4/4/010
  • 18 Shewchuk, J.R.: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. Tech. rep., Carnegie Mellon University, Pittsburgh (1994), https://www.cs.cmu.edu/{~}quake-papers/painless-conjugate-gradient.pdf
  • 19 Tupin, F., Inglada, J., Nicolas, J.M.: Remote Sensing Imagery. Wiley-ISTE (2014)
  • 20 Vogel, C.R.: Computational Methods for Inverse Problems. SIAM (2002)
  • 21 van der Vorst, H.A.: Iterative Krylov Methods for Large Linear Systems. Cambridge University Press, Cambridge (2003)


Figure 1: Wrapped phase map used in the experiments.
Figure 2: Iterations employed for the different test.
Figure 3: Computational time used for the different test.
Figure 4: Resultant unwrapped phase map using Algorithm 1. The phase maps was wrapped for purposes of illustration.