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 aswhere 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
(1) 
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 illposed 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, subsampling 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 pathfollowing 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 wrappedphase 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 wrappedphase 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 2D phase unwrapping is defined as ^{9}
(2) 
where and is the domain of integration. To obtain the solution of the problem expressed in Eq. (2
), the firstorder optimality condition or EulerLagrange equation has to be derived, resulting in the following partial differential equation (PDE)
(3) 
with boundary conditions
where
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
respectively.
Hence the numerical approximation of the Eq. (3) is given by
(4) 
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
(5)  
with boundary conditions
where for simplicity and without loss of generality we consider The terms y are defined as
(6)  
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
(7)  
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.
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}:

where
is the identity matrix; with this matrix, PCG becomes CG
^{10}. 
where is the diagonal of matrix this technique is known as Jacobi preconditioning.

the incomplete LU factorization with no fillin.

the incomplete Cholesky factorization with no fillin.

the successive overrelaxation factorization.
The main reason for selecting these preconditioning techniques was their no fillin 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) 64bit 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.
matrix  

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 
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}:
(8) 
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.
References
 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 phaseunwrapping 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.: Twodimensional 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.: TwoDimensional Phase Unwrapping: Theory, Algorithms, and Software. WileyInterscience, New York (1998)
 9 Ghiglia, D.C., Romero, L.A.: Minimum L pnorm twodimensional 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/s1066501698497
 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/02665611/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/{~}quakepapers/painlessconjugategradient.pdf
 19 Tupin, F., Inglada, J., Nicolas, J.M.: Remote Sensing Imagery. WileyISTE (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)