Vol. 24, No.

25 | 12 Dec 2016 | OPTICS EXPRESS 28713

Phase calibration unwrapping algorithm for

phase data corrupted by strong decorrelation
speckle noise
HAITING XIA,1,2,3 SILVIO MONTRESOR,3 RONGXIN GUO,1,2 JUNCHANG LI,5
FENG YAN,1,2 HEMING CHENG,1 AND PASCAL PICART3,4,*
1
Faculty of Civil Engineering and Mechanics, Kunming University of Science and Technology, Kunming
650500, China
2
Key Laboratory of Yunnan Higher Education Institutes for Mechanical Behavior and Microstructure
Design of Advanced Materials, Kunming University of Science and Technology, Kunming 650500,
China
3
Université du Maine, CNRS UMR 6613, LAUM, Avenue Olivier Messiaen, 72085 Le Mans Cedex 9,
France
4
École Nationale Supérieure d’Ingénieurs du Mans, rue Aristote, 72085 Le Mans Cedex 9, France
5
Faculty of Science, Kunming University of Science and Technology, Kunming 650500, China
*
pascal.picart@univ-lemans.fr

Abstract: Robust phase unwrapping in the presence of high noise remains an open issue.
Especially, when both noise and fringe densities are high, pre-filtering may lead to phase
dislocations and smoothing that complicate even more unwrapping. In this paper an approach
to deal with high noise and to unwrap successfully phase data is proposed. Taking into
account influence of noise in wrapped data, a calibration method of the 1st order spatial phase
derivative is proposed and an iterative approach is presented. We demonstrate that the
proposed method is able to process holographic phase data corrupted by non-Gaussian
speckle decorrelation noise. The algorithm is validated by realistic numerical simulations in
which the fringe density and noise standard deviation is progressively increased. Comparison
with other established algorithms shows that the proposed algorithm exhibits better accuracy
and shorter computation time, whereas others may fail to unwrap. The proposed algorithm is
applied to phase data from digital holographic metrology and the unwrapped results
demonstrate its practical effectiveness. The realistic simulations and experiments demonstrate
that the proposed unwrapping algorithm is robust and fast in the presence of strong speckle
decorrelation noise.
© 2016 Optical Society of America
OCIS codes: (100.5088) Phase unwrapping; (120.5050) Phase measurement; (110.4280) Noise in imaging systems;
(100.3175) Interferometric imaging; (100.3020) Image reconstruction-restoration; (100.2980) Image enhancement;
(090.1995) Digital holography.

References and links

1. D. C. Ghiglia and M. D. Pritt, Two-Dimensional Phase Unwrapping: Theory, Algorithms, and Software (Wiley,
1998).
2. H. Huang, “A fast multi-baseline and multi-frequency band phase-unwrapping algorithm,” Measurement 49,
401–406 (2014).
3. R. Schlögel, C. Doubre, J. P. Malet, and F. Masson, “Landslide deformation monitoring with ALOS/PALSAR
imagery: A D-InSAR geomorphological interpretation method,” Geomorphology 231, 314–330 (2015).
4. P. Andris and I. Frollo, “Simple and accurate unwrapping phase of MR data,” Measurement 42(5), 737–741
(2009).
5. P. Daga, T. Pendse, M. Modat, M. White, L. Mancini, G. P. Winston, A. W. McEvoy, J. Thornton, T. Yousry, I.
Drobnjak, J. S. Duncan, and S. Ourselin, “Susceptibility artefact correction using dynamic graph cuts:
application to neurosurgery,” Med. Image Anal. 18(7), 1132–1142 (2014).
6. B. N. Li, X. Shan, K. Xiang, N. An, J. Xu, W. Huang, and E. Kobayashi, “Evaluation of robust wave image
processing methods for magnetic resonance elastography,” Comput. Biol. Med. 54, 100–108 (2014).
7. F. Da and H. Huang, “A fast, accurate phase unwrapping method for wavelet-transform profilometry,” Opt.
Commun. 285(4), 421–432 (2012).

#278463 http://dx.doi.org/10.1364/OE.24.028713
Journal © 2016 Received 10 Oct 2016; revised 15 Nov 2016; accepted 15 Nov 2016; published 5 Dec 2016
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28714

8. L. Song, X. Dong, J. Xi, Y. Yu, and C. Yang, “A new phase unwrapping algorithm based on Three Wavelength
Phase Shift Profilometry method,” Opt. Laser Tech. 45, 319–329 (2013).
9. F. Chen and X. Su, “Phase-unwrapping algorithm for the measurement of 3D object,” Optik (Stuttg.) 123(24),
2272–2275 (2012).
10. B. Tayebi, M. R. Jafarfard, F. Sharif, Y. S. Song, D. Har, and D. Y. Kim, “Large step-phase measurement by a
reduced-phase triple-illumination interferometer,” Opt. Express 23(9), 11264–11271 (2015).
11. Z. Cheng, D. Liu, Y. Yang, T. Ling, X. Chen, L. Zhang, J. Bai, Y. Shen, L. Miao, and W. Huang, “Practical
phase unwrapping of interferometric fringes based on unscented Kalman filter technique,” Opt. Express 23(25),
32337–32349 (2015).
12. P. Girshovitz and N. T. Shaked, “Fast phase processing in off-axis holography using multiplexing with complex
encoding and live-cell fluctuation map calculation in real-time,” Opt. Express 23(7), 8773–8787 (2015).
13. S. Liu, W. Xiao, F. Pan, F. Wang, and L. Cong, “Complex-amplitude-based phase unwrapping method for
digital holographic microscopy,” Opt. Lasers Eng. 50(3), 322–327 (2012).
14. Y. T. Zhang, M. J. Huang, H. R. Liang, and F. Y. Lao, “Branch cutting algorithm for unwrapping photoelastic
phase map with isotropic point,” Opt. Lasers Eng. 50(5), 619–631 (2012).
15. J. C. de Souza, M. E. Oliveira, and P. A. dos Santos, “Branch-cut algorithm for optical phase unwrapping,” Opt.
Lett. 40(15), 3456–3459 (2015).
16. M. Zhao, L. Huang, Q. Zhang, X. Su, A. Asundi, and Q. Kemao, “Quality-guided phase unwrapping technique:
comparison of quality maps and guiding strategies,” Appl. Opt. 50(33), 6214–6224 (2011).
17. X. Su and W. Chen, “Reliability-guided phase unwrapping algorithm: a review,” Opt. Lasers Eng. 42(3), 245–
261 (2004).
18. S. Yuqing, “Robust phase unwrapping by spinning iteration,” Opt. Express 15(13), 8059–8064 (2007).
19. V. V. Volkov and Y. Zhu, “Deterministic phase unwrapping in the presence of noise,” Opt. Lett. 28(22), 2156–
2158 (2003).
20. J. J. Chyou, S. J. Chen, and Y. K. Chen, “Two-dimensional phase unwrapping with a multichannel least-mean-
square algorithm,” Appl. Opt. 43(30), 5655–5661 (2004).
21. Y. Guo, X. Chen, and T. Zhang, “Robust phase unwrapping algorithm based on least squares,” Opt. Lasers Eng.
63, 25–29 (2014).
22. M. Rivera, F. J. Hernandez-Lopez, and A. Gonzalez, “Phase unwrapping by accumulation of residual maps,”
Opt. Lasers Eng. 64, 51–58 (2015).
23. W. He, L. Xia, and F. Liu, “Sparse-representation-based direct minimum L (p) -norm algorithm for MRI phase
unwrapping,” Comput. Math. Methods Med. 2014, 134058 (2014).
24. F. Palacios, E. Gonçalves, J. Ricardo, and J. L. Valin, “Adaptive filter to improve the performance of phase-
unwrapping in digital holography,” Opt. Commun. 238(4), 245–251 (2004).
25. Q. Kemao, “Applications of windowed Fourier fringe analysis in optical measurement: a review,” Opt. Lasers
Eng. 66, 67–73 (2015).
26. J. C. Estrada, M. Servin, and J. A. Quiroga, “Noise robust linear dynamic system for phase unwrapping and
smoothing,” Opt. Express 19(6), 5126–5133 (2011).
27. M. A. Navarro, J. C. Estrada, M. Servin, J. A. Quiroga, and J. Vargas, “Fast two-dimensional simultaneous
phase unwrapping and low-pass filtering,” Opt. Express 20(3), 2556–2561 (2012).
28. J. F. Weng and Y. L. Lo, “Robust detection scheme on noise and phase jump for phase maps of objects with
height discontinuities--theory and experiment,” Opt. Express 19(4), 3086–3105 (2011).
29. J. Weng and Y. Lo, “Novel rotation algorithm for phase unwrapping applications,” Opt. Express 20(15), 16838–
16860 (2012).
30. H. Y. 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).
31. H. T. Xia, R. X. Guo, Z. B. Fan, H. M. Cheng, and B. C. Yang, “Non-invasive mechanical measurement for
transparent objects by digital holographic interferometry based on iterative least-squares phase unwrapping,”
Exp. Mech. 52(4), 439–445 (2012).
32. Th. Kreis, Holographic Interferometry – Pinciples and Methods (Akademie Verlag 1996).
33. J. W. Goodman, Speckle Phenomena in Optics (Roberts and Company Publishers, 2006).
34. T. C. Poon, Digital Holography and Three-Dimensional Display: Principles and Applications (Springer-Verlag,
2010).
35. P. Picart, New Techniques in Digital Holography (ISTE-Wiley 2015).
36. J. Poittevin, P. Picart, F. Gautier, and C. Pézerat, “Quality assessment of combined quantization-shot-noise-
induced decorrelation noise in high-speed digital holographic metrology,” Opt. Express 23(24), 30917–30932
(2015).
37. S. Montresor and P. Picart, “Quantitative appraisal for noise reduction in digital holographic phase imaging,”
Opt. Express 24(13), 14322–14343 (2016).
38. M. Karray, P. Slangen, and P. Picart, “Comparison between digital Fresnel holography and digital image-plane
holography: the role of the imaging aperture,” Exp. Mech. 52(9), 1275–1286 (2012).
39. Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. 43(13), 2695–2702 (2004).
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28715

1. Introduction
Phase unwrapping is a process of retrieving actual non-wrapped phase from the wrapped
phase obtained through numerical operations based on arctangent function [1]. Phase
unwrapping is an unavoidable process that can be encountered in many applications for which
the phase is of primary interest: synthetic aperture radar imaging (SAR) [2,3], magnetic
resonance imaging (MRI) [4–6], surface shape measurement [7–9], interferometry [10,11],
and digital holography [12,13]. In practical applications, phase unwrapping may be quite
difficult because of the presence of residues, noise or other phase singularities. In the last
decades, many phase unwrapping approaches were developed. As a general rule, these
methods can be classified as path-following [1,14–18] or minimum-norm methods [1,19–23].
In addition, branch-cut, quality map or masks are usually utilized to remove residues, noise or
holes with high reliability. Among them, quality-guided, Flynn’s minimum discontinuity and
minimum Lp-norm (L0) algorithm have exhibited the best performance in the presence of
noise [1]. In order to improve the performance of phase unwrapping algorithms in the
presence of high noise, spatial filtering can be applied so as to reduce noise and to clean the
wrapped phase map prior to unwrapping [24,25]. However, filtering algorithms may smear
the phase jumps or get rid of some useful information. So, phase unwrapping algorithms
which directly operate on the phase map corrupted by noise and holes have to be developed.
For example, a first-order feedback system in phase unwrapping was developed for directly
operating on noisy phase maps [26,27]. In addition, such an approach allows filtering some
noise since it is based on an infinite impulse response low-pass filter. However, this approach
is unstable in the presence of high noise. Weng and Lo recently proposed a noise and phase
jump detection scheme for noisy phase maps containing height discontinuities and a rotation
algorithm which operates on all the pixels in the wrapped phase map without prior filtering
[28,29]. In [30], Huang et al proposed a path-independent phase unwrapping method which
de-noises the phase gradient by total-variation minimization to separate the noise from the
genuine phase gradient map. The method is most suitable for unwrapping phase maps without
abrupt phase changes and its computation cost is expensive for the case of total-variation.
Recently, Xia et al presented a phase unwrapping algorithm based on least-squares and
iterations of unwrapped phase errors [31]. With this algorithm, phase unwrapping can be
performed in the presence of noise and holes and accurate unwrapped results were obtained.
In subsequent applications, however, it was found that this method fails when the local
standard deviation of noise is higher than 0.8rad.
In digital holographic metrology and speckle-based approaches, phase data is generally
obtained from an arctangent formula. Quantitative phase measurement is carried out by a
subtraction between the phase from the object under interest and a “reference” phase [32–37].
Unfortunately, a phenomena of speckle decorrelation occurs when the object is time-varying
or submitted to surface changes [36–38]. This decorrelation adds a high spatial frequency
noise to the useful signal [37]. As pointed out above, smoothing the phase map is an option,
but when the fringes are too close together, there exists a risk of generating phase
dislocations. So, the phase unwrapping algorithm which directly operates on the phase map
highly corrupted by noise has to be developed. In this paper, we present a robust algorithm
based on least-squares, iteration and calibration of the 1st order spatial phase derivative. The
proposed algorithm is validated by using simulated phase data which are highly corrupted by
realistic speckle decorrelation noise. The paper demonstrates that this method can retrieve the
phase with high accuracy in the presence of noise without any spatial filtering prior to the
phase unwrapping. In addition, comparison with other well-established algorithms exhibits its
accuracy and rapidity. This paper is organized as follows: in section 2, we describe the
theoretical basics for the calibrated phase unwrapping algorithm; section 3 discusses on the
realistic simulation of corrupted phase with high decorrelation noise. In section 4, we discuss
on the evaluation of the proposed approach with comparison with other established
unwrapping algorithms. Errors are quantitatively evaluated thanks to the database constituted
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28716

with the realistic simulation. Section 5 discusses on application of the proposed unwrapping
approach to experimental data including high decorrelation noise. Conclusions and
perspectives to the study are drawn in section 6.
2. Principle of proposed phase unwrapping algorithm
2.1 Calibration of phase derivatives
In this paper, φij represents the true phase and ψij (∈[−π, + π], that is modulo 2π) represents
the wrapped phase at the grid point (i,j) of a phase map. The wrapped phase is extracted from
a computation process based on the arctangent operator. This process is currently referred as
the “wrapping operator”. In this paper, the wrapping operator is symbolized as [1]:

W (φij ) = ψ ij ( i = 0,1..., M − 1; j = 0,1...N − 1) , (1)

where −π≤ψij≤π, M, N are respectively the number of grid points with respect to the i index
and the j index. The 1st order spatial wrapped phase derivatives are defined as:

(
Δ ijx = W ψ ( i +1) j −ψ ij ) ( i = 0,1..., M − 2; j = 0,1...N − 1)
Δ ijx = 0 otherwise
(2)
(
Δ ijy = W ψ i ( j +1) −ψ ij ) ( i = 0,1..., M − 1; j = 0,1...N − 2 )
Δ ijy = 0 otherwise,

where Δijx and Δijy are respectively the difference with respect to the i and the j indexes. The
phase derivative is an important amount for phase unwrapping. For the noise-free wrapped
phase, the wrapped phase derivatives of wrapped phase are the same as those of original
unwrapped phase. So the original phase can be unwrapped by the phase derivatives of
wrapped phase. However, the presence of noise will generate errors between noisy and noise-
free phase derivatives and makes unwrapping difficult, even failed. In order to illustrate the
influence of noise on phase unwrapping, Fig. 1 shows two different profiles of wrapped phase
derivatives extracted from phase maps from section 3. Figure 1(a) shows the noisy profile and
the noise-free profile. Figure 1(b) shows another noisy profile and its corresponding noise-
free profile. The red line corresponds to the standard deviation of the noisy phase derivatives.

Fig. 1. Examples of profiles of phase derivatives of noisy wrapped phase and noise-free
wrapped phase, (a) first example, blue points: noisy profile, dark line: noise-free profile, (b)
second example, blue points: noisy profile, dark line: noise-free profile; the red line correspond
to the standard deviation.

From these figures can be seen that there exist errors between noisy and noise-free phase
derivatives and phase derivatives are, as a general rule, mostly distributed near the noise-free
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28717

derivatives. These errors will make unwrapping very difficult or failing and have to be
removed. In [26,27,30], filtering is applied on the noisy phase derivatives to remove these
errors. However, in the presence of high noise these errors are difficult to be removed by
filtering. Here, we propose a calibration approach to calibrate the phase derivatives exhibiting
large errors:

Δ ijx = sgn ( Δ ijx ) Gx if Δ ijx ≥ Tx

Δ ijx = Δ ijx otherwise
(3)
Δ = sgn ( Δ ijy ) Gy
y
ij if Δ ijy ≥ Ty
Δ ijy = Δ ijy otherwise,

where sgn(…) is the signum function, Tx and Ty are the thresholds, Gx and Gy are the
calibrated phase derivatives. These parameters are defined as (E[…] means statistical
average):

(
Tx = E ( Δ ij )  − E  Δ ij 
 x 2 )
2
x

 , (4)
T = E ( Δ y )2  − E  Δ y 
( )
2

 y  ij   ij 

and,
 1 i = M −1 j = N −1 x
G
 x =   Δij
MN i = 0 j = 0

 . (5)
i = M −1 j = N −1
G = 1
 y MN   Δ ij
y

 i =0 j =0

This calibration approach means that the phase derivatives whose values are larger than the
standard deviation of phase derivatives are replaced by the average value of phase derivatives.
2.2 Phase unwrapping
The least-squares phase unwrapping algorithm seeks the unwrapped solution by minimizing
the differences between the discrete partial derivatives of the wrapped phase data and those of
the unwrapped solution. The least-squares solution is φij that minimizes quantity s such that
[1]:
M − 2 N −1 2 M −1 N − 2 2
s=   φ(
i =0 j =0
i +1) j
− φij − Δ ijx +   φi ( j +1) − φij − Δ ijy .
i =0 j =0
(6)

Equation (5) can be reduced to the discrete Poisson equation with Neumann boundary
conditions [1]:

(φ( i +1) j ) ( )
− 2φij + φ(i −1) j + φi ( j +1) − 2φij + φi ( j −1) = ρij , (7)

where

( ) (
ρij = Δijx − Δ(xi −1) j + Δijy − Δiy( j −1) . ) (8)

The discrete Poisson equation can be solved by many methods such as fast Fourier transform
(FFT), discrete cosine transform (DCT) or the multi-grid method. In this algorithm, the DCT
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28718

method is selected to solve the least-squares phase unwrapping problem. Performing the 2D
DCT to Eq. (7) yields the exact solution of φij in the DCT reciprocal space according to:
ρˆ ij
φˆij = , (9)
2 cos (π i / M ) + 2 cos (π j / N ) − 4

where φˆij and ρˆij are respectively the DCT of φij and ρij. Then computing the inverse 2D
DCT (named IDCT) to φˆij yields the least-squares unwrapped phase.

2.3 Unwrapping with iterations

Theoretically, the solution obtained from Eq. (9) is the exact one. However, there exist errors
between the unwrapped phase and the true phase due to noise and to the smoothing
performance of the least-squares method. In the proposed approach, iterations of unwrapped
phase errors are utilized to seek more accurate results. Assume that k is the iteration number.
At the k-th iteration, the least-squares unwrapped phase φij,k is calibrated to obtain the
calibrated unwrapped phase φij,k according to:

 φij , k −ψ ij 
ϕij , k = ψ ij + 2π round  , (10)
 2π 
where round(…) is the rounding operator. This equation can be rewritten as
ϕij , k −ψ ij  φij , k −ψ ij 
= round  . (11)
2π  2π 
Define the k-th phase error as:
Δψ ij , k = ϕij , k − φij , k . (12)

On one hand, from Eqs. (11) and (12) can be seen that the difference between the calibrated
unwrapped phase φij,k and the wrapped phase ψij is equal to an integer times 2π. If φij,k is
continuous and non-interrupt, the calibrated unwrapped phase φij,k is the true phase. On the
other hand, the least-squares unwrapped phase φij,k and the calibrated unwrapped phase φij,k
have the same wrap counts. So, the phase error Δψij,k is in the range [−π, + π]. Because the
least-squares unwrapped phase φij,k is continuous, if φij,k is continuous, then the error Δψij,k
must be also continuous, and, otherwise Δψij,k is wrapped into [−π, + π]. Therefore if Δψij,k is
wrapped, φij,k is also interrupted and not equal to the true phase. In this case, Δψij,k can be
unwrapped by least-squares algorithm and acquired results are added to the previous least-
squares unwrapped phase φij,k to obtain unwrapped results closer to the true phase. Therefore,
the average value of |Δψij,k−Δψij,k-1| is used to judge if the iteration process requires to
continue or to be stopped. If this value is smaller than a value fixed at Err, then iterations are
stopped and the final unwrapped phase is set to φij,k. Otherwise, phase unwrapping is
performed anew with Δψij,k to obtain the k + 1th unwrapped phase. The unwrapped phase
being the closest to the true phase can be obtained by:
φij , k +1 ⇐ φij , k + φij , k +1 . (13)

Substitute Eq. (13) to Eq. (10) to obtain the k + 1th calibrated unwrapped phase. Perform the
next step iteration until <|Δψij,k−Δψij,k-1|>≤Err (<…> stands for average value) or stop when
reaching the maximum number of iterations at value Nit. The block diagram of the proposed
phase unwrapping algorithm is summarized in Fig. 2. In the following, this approach is
referred as CPULSI (Calibrated Phase Unwrapping based on Least-Squares and Iterations).
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28719

Fig. 2. Flow chart of the unwrapping algorithm; DCT means Discrete Cosine Transform, IDCT
means Inverse Discrete Cosine Transform, <…> means average value.

3. Realistic simulation of noisy wrapped phase data

3.1 Principle
In order to evaluate the robustness of the proposed unwrapping approach, a realistic
numerical simulation was carried out and is based on that proposed in [37]. The goal of such
simulation is to produce phase maps corrupted by speckle decorrelation noise with the
adequate probability density function. In experiments, the amount of speckle phase
decorrelation is naturally controlled by the fringe density, which is mainly related to the
object surface modifications. The larger the surface deformation, the higher the number of
fringes, and hence the higher the phase decorrelation between the two states of the object
phase. The arrangement to produce speckle phase decorrelation was discussed in [37], but we
remind in this section the basics fundamentals. Suppose any object with rough surface
(natural object) illuminated in the visible range of light. Figure 3 shows the scheme, whereas
equations supporting the simulation are provided in [37]. Note that the value of Ru in Fig. 3 is
adjusted to control the speckle grain size in the image plane and to simulate realistic images.
The procedure is as follows: first calculate the convolution equation with the random surface
to produce a random “reference” speckle field in the image plane, second, calculate the
convolution equation with the random surface to which is added the surface deformation to
produce a random “deformed” speckle field in the image plane, last calculate the phase
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28720

difference between the two previous speckle phases from the two calculated optical fields.
Phase decorrelation occurs because of the spatial filtering produced by the diaphragm in the
back focal plane of the first lens in Fig. 3. Note that the standard deviation of phase noise
obtained in the phase difference image is related to the amplitude of the deformation, and thus
to the fringe density. Noise standard deviation and fringe density are closely linked as very
good friends. Increasing the fringe density leads to an increase of the noise standard
deviation. In the simulation, the noise standard deviation cannot be easily controlled. The
only things that can be quite controlled are the amplitude of surface changes and the size of
the speckle grains. The raw result of the simulation is that the speckle phase decorrelation
noise follows the same statistics as given in Eq. (4) in [37], that is, a non-Gaussian noise.

Fig. 3. Arrangement to produce speckle phase decorrelation in phase change due to surface
deformation.

3.2 Database for evaluation

Two sets including 10 phase maps each with increasing noise and fringe density were
computed. The data were obtained with an image size at 1024 × 1024 pixels, and a speckle
grain having a size equal to 4 pixels, which corresponds to what is usually encountered in
practical situation in speckle/holographic phase metrology [38]. Figure 4 exhibits the two sets
of modulo 2π noisy phase data. For each set, one observes a progressive increase of the fringe
density and phase noise (from left to right). The two first rows show the ten wrapped phases
of Set 1 and the two last rows that of Set 2. We define σt the total standard deviation of the
noise in the phase map, and σlmax the maximum value of the local noise standard deviation in a
small patch (31 × 31 pixels) in the phase map. Note that these two values cannot be equal
because speckle decorrelation noise is non stationary.
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28721

Fig. 4. Outputs from the numerical simulation; two sets of phase data with progressive increase
of fringe density and phase noise. Two first rows, Set 1, from left to right increase of fringe
density and phase noise, two last rows, Set 2, from left to right increase of fringe density and
phase noise.

The fringe density ρf is defined as the number of fringes per pixel which is calculated in a
small patch with K × L pixels by
i+K /2 j+L/2
1
  (Δ ) + (Δ )
2 2
ρ f ,ij = x y
. (14)
2π KL
pq pq
p =i − K / 2 q = j − L / 2

The maximum fringe density is denoted by ρfmax. In this paper, K = L = 100. The values of the
different parameters qualifying the fringe pattern are given in Table 1 for the two sets of
phase data. Parameters σt and σlmax are calculated from simulated noise and ρf (average fringe
density of the whole phase map) and ρfmax are calculated from simulated noise-free original
phase. Table 1 shows that σlmax and ρfmax in Set 2 are larger than for Set 1. This means that the
fringes of the phase maps in Set 2 are closer and noisier in few particular local areas.
Figure 5 shows the relationships between qualifying parameters. Figure 5(a) shows σlmax
versus σt and Fig. 5(b) shows ρfmax versus σt respectively for each phase map of the two sets.
As can be seen, the relation between the fringe density and the standard deviation is not linear
and depends on the fine structure of the fringe pattern. By using these two sets of phase maps,
the unwrapping algorithm is tested with a certain kind of “pattern density and structure
diversity”.
The next section discusses on the evaluation of the proposed unwrapping algorithm and
comparison with other established approaches.
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28722

Table 1. Standard deviation of noise and fringe density of the phase maps for the two
sets.
Set 1 Set 2
σt (rad) σlmax σt (rad) σlmax ρfmax
(rad) ρfmax ρfmax/ ρf (rad) (pixel−1 ρfmax/ ρf
N° (pixel−1) )
1 0.25 0.993 0.0016 2.34 0.25 1.036 0.0093 4.97
2 0.30 0.949 0.0053 2.34 0.30 0.824 0.0146 4.97
3 0.35 0.954 0.0084 2.34 0.35 0.947 0.0199 4.97

4 0.40 1.106 0.0105 2.34 0.40 1.009 0.0279 4.97

5 0.45 0.930 0.0158 2.34 0.45 1.220 0.0372 4.97

6 0.50 1.018 0.0200 2.34 0.50 1.248 0.0478 4.97

7 0.55 1.016 0.0253 2.34 0.55 1.328 0.0611 4.97

8 0.60 1.111 0.0321 2.34 0.60 1.521 0.0731 4.97

9 0.65 1.189 0.0395 2.34 0.65 1.761 0.0904 4.97
10 0.70 1.267 0.0474 2.34 0.70 1.805 0.1090 4.97

Fig. 5. Relationships between qualifying parameters, (a) σlmax versus σt for each phase map of
the two sets, (b) ρfmax versus σt for each phase map of the two sets.

4. Evaluation of the proposed algorithm

4.1 Methodology
In order to evaluate the performance of the proposed algorithm, five other different phase
unwrapping algorithms were selected: Goldstein branch cut algorithm (namely “Gold”) [1],
quality guided path following algorithm (“Quality”) [1], Flynn’s minimum-weight-
discontinuity algorithm (“Flynn”) [1], minimum Lp norm algorithm (“Lp”) [1] and phase
unwrapping based on least-squares and iteration (“PULSI”) [31]. These algorithms and
CPULSI were used to unwrap the phase maps of the two sets of data. The parameters in
CPULSI were set to Nit = 300, and Err = 0.001.
4.2 Unwrapped phase maps
In order to illustrate the results obtained from the six algorithms, the unwrapped phase maps
for fringe pattern n°10 for Set 1 and fringe pattern n°8 for Set 2 are shown in Figs. 6-9. Figure
6 shows the original phase and wrapped phase of Set 1 for fringe pattern n°10 (refer to Fig.
4).
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28723

Fig. 6. Phase maps n°10 of Set 1, (a) original unwrapped phase, (b) noisy wrapped phase
(colorbars given in radian units).

Figure 7 shows the unwrapped phase maps obtained from the six different algorithms.
Figure 7 shows that, with the naked eyes, results obtained from Flynn, Lp and CPULSI are
better than those from other algorithms.

Fig. 7. Unwrapped phase obtained for phase map n°10 of Set 1, (a) Gold, (b) Quality, (c)
Flynn, (d) Lp, (e) PULSI, (f) CPULSI (colorbars given in radian units).

Figure 8 shows the original phase and wrapped phase of Set 2 for fringe pattern n°8 (refer
to Fig. 4). Figure 9 shows the corresponding unwrapped results. It was found that the six
algorithms were not able to correctly unwrap all the fringe patterns from Set 2. Especially,
pattern n°7 is the last one that can be correctly unwrapped by Flynn, whereas, pattern n°8 is
the last one for CPULSI. Unwrapping can be considered to have failed for the other
algorithms below pattern n°7 for Set 2. From Fig. 9, one can appreciate that the unwrapped
results obtained from CPULSI are better than with other algorithms. These results
demonstrate that the proposed algorithm is robust and accurate in the presence of high speckle
decorrelation noise.
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28724

Fig. 8. Phase maps n°8 of Set 2, (a) original unwrapped phase, (b) noisy wrapped phase
(colorbars given in radian units).

4.3 Quantitative evaluation of phase errors

There are several ways to investigate the phase errors given by each algorithm. Especially,
from the numerical simulation of phase maps, one knows the original wrapped phase, the
original unwrapped phase, and the noise included in the wrapped noisy map. From the initial
noise and unwrapped original phase, one can get the theoretical unwrapped noisy phase. This
data can be compared to the raw unwrapped phase map from the different algorithms. If the
algorithm does not introduce error, then the theoretical unwrapped noisy phase should be
equal to the unwrapped noisy phase. This phase error was estimated and is provided in Fig. 10
and Fig. 11 for respectively Set 1, pattern n°10, and Set 2, pattern n°8. One can observe that
the phase errors for CPULSI are less than for other algorithms, and that they look as sparsely
points. Such point-like error can be easily reduced with low pass filtering applied to
unwrapped phase.

Fig. 9. Unwrapped phase obtained for phase map n°8 of Set 2, (a) Gold, (b) Quality, (c) Flynn,
(d) Lp, (e) PULSI, (f) CPULSI (colorbars given in radian units).
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28725

Fig. 10. Error maps of unwrapped phase for fringe pattern n°10 of Set 1, (a) Gold, (b) Quality,
(c) Flynn, (d) Lp, (e) PULSI, (f) CPULSI (colorbars given in radian units).

Fig. 11. Error maps of unwrapped phase for fringe pattern n°8 of Set 2, (a) Gold, (b) Quality,
(c) Flynn, (d) Lp, (e) PULSI, (f) CPULSI (colorbars given in radian units).

In order to evaluate the noise/error introduced by each algorithm, the standard deviation of
the phase difference between these two outputs is evaluated. If equal to 0, this means that the
unwrapped phase is quite similar to the theoretical noisy phase, and that the algorithm
unwraps the phase without introducing any phase error or noise filtering. Tables 2 and 3
summarize the standard deviations of phase differences between each unwrapped phase and
original noisy one, for each selected algorithm and each data set. Note that the standard
deviations for CPULSI are the least, and that they are correlated with the error maps of Figs.
10(f) and 11(f). For Set 2, when the initial noise standard deviation increases two much
(σt>0.6rad and σlmax>1.1rad for Set 2), there is no algorithm able to successfully unwrap the
phase maps. This can be explained by the high local density of fringes (or high maximum
local standard deviation of noise) in the low-right corner of the image (see Fig. 3).
However, Tables 2 and 3 clearly show that for both sets, the CPULSI approach provides
the best results since the estimated standard deviations are the smallest ones, and that it can
unwrap up to 8 phase maps in Set 2, whereas the Flynn approach can only deal with 7 phase
maps.
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28726

Table 2. Standard deviations of phase difference between unwrapped phase and original
noisy phase (rad) for Set n°1
Set 1
N° Gold Quality Flynn Lp PULSI CPULSI
1 0.173 0.175 0.165 0.130 0.164 0.080
2 0.192 0.166 0.169 0.136 0.183 0.097
3 0.236 0.194 0.202 0.161 0.428 0.127
4 0.328 0.240 0.253 0.186 2.944 0.150
5 0.406 0.260 0.292 0.219 3.987 0.178
6 0.542 0.317 0.350 0.259 6.026 0.218
7 0.804 0.382 0.428 0.299 8.676 0.256
8 1.343 0.471 0.498 0.340 12.800 0.315
9 2.108 0.656 0.567 0.393 17.417 0.351
10 5.320 0.785 0.660 0.472 23.785 0.419
Table 3. Standard deviations of phase difference between unwrapped phase and original
noisy phase (rad) for Set n°2
Set 2
N° Gold Quality Flynn Lp PULSI CPULSI
1 0.138 0.116 0.126 0.099 0.198 0.072
2 0.177 0.147 0.168 0.127 0.849 0.099
3 0.233 0.194 0.206 0.155 1.488 0.123
4 0.326 0.250 0.266 0.192 3.076 0.158
5 0.720 0.399 0.333 0.234 4.476 0.195
6 1.724 0.797 0.414 0.308 6.096 0.241
7 3.238 0.819 0.484 0.871 7.913 0.292
8 5.591 1.212 1.346 1.751 10.354 0.368
9 9.242 2.935 2.870 3.461 13.169 1.774
10 15.773 5.064 7.468 6.555 16.860 3.652
Figure 12 shows the plots of the estimated standard deviations of the phase errors versus
the initial input standard deviation of the decorrelation noise. The vertical scale is limited to
3rad. For Set 1 and Set 2, Fig. 12(a) and Fig. 12(b) confirm that CPULSI gives the best
results.

Fig. 12. Standard deviations of the errors between the unwrapped phase and the noisy original
phase: (a) for Set1, (b) for Set2.

Another way for qualifying the results is to compute the ratios between the peak-to-valley
phase error and the peak-to-valley of the original unwrapped noise-free phase. This helps to
understand how unwrapping produces hot spots in the unwrapped phase map. In case of
correct or successful unwrapping, this ratio would decreases to 0, meaning that hot spots are
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28727

quite weaker than the full phase dynamics. Figures 13(a) and 13(b) show such ratio obtained
for respectively Set 1 and Set 2, versus the initial noise standard deviation. One can observe
that the ratios given by CPULSI are the least ones.

Fig. 13. Ratio of peak-valley value of unwrapped phase errors divided by that of original
phase: (a) for Set 1, (b) for Set 2; PVphase: peak-to-valley of the original unwrapped phase.

Even though the standard deviations of the two sets of phase maps change in the same
way, the unwrapped errors exhibit different patterns (see Fig. 10 and Fig. 11). Figure 14
shows the standard deviation of unwrapping error versus the fringe density and the maximum
local standard deviation for the most successful algorithm Flynn, Lp and CPULSI. Figures
14(a) and 14(b) clearly show that the unwrapped error is closely related to the maximum
fringe density, or to maximum local standard deviation of noise.

Fig. 14. Evolution of the standard deviation of errors, (a) versus the maximum fringe density,
and (b) versus the maximum local standard deviation of noise, for the two sets of fringe
patterns and algorithms Flynn, Lp and CPULSI.

4.4 Computation time

The computation cost of the proposed approach was evaluated and compared with the other
selected algorithms. The computations were performed with MATLAB on a laptop equipped
with Intel Core i7-5500U CPU @2.4GHz and 8GB RAM. Table 4 provides the average
computation times of the different algorithms. From Table 3 can be seen that the computation
time of CPULSI is less than that of other algorithms, except for PULSI and Gold. However,
PULSI and Gold fail in presence of high noise (refer to Table 3). This means that the
proposed algorithm can be considered as a fast unwrapping algorithm in the presence of high
speckle noise. Note that the computation times in Table 4 are quite high. However, such
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28728

results do not prevent from optimized computation, such as with GPU for example, which
would be quite faster.
Table 4. Average computation times of the different algorithms.
Gold Quality Flynn Lp PULSI CPULSI
Computatio
122.8 564.5 783.7 485.0 5.9 145.9
n time (s)

5. Application to experimental holographic phase data

In order to test the proposed unwrapping approach on real data, two experimental phase maps
obtained from digital holography were used as inputs to the algorithm. Figure 15(a) shows the
phase map related to the stress field around a crack in a transparent plate tested by digital
transmission holographic interferometry [31]. The phase map includes 1000 × 750 pixels.
Figure 15(b) shows the phase map measured with digital Fresnel holography when submitting
a rough object to a mechanical loading [38]. This phase map includes 1360 × 1024 pixels. As
can be seen, there exists high speckle decorrelation noise and data-missing regions in both
phase maps. The noise was estimated for both phase maps. To get a robust estimation, the
windowed Fourier transform algorithm [39] was applied on the raw phase maps and, from the
results, noise was estimated by subtracting the raw phase. The fringe density was also
estimated. For Fig. 15(a), one gets σt = 0.515rad, σlmax = 1.429rad, ρfmax = 0.0435pixel−1, and
for Fig. 15(b), σt = 0.859rad, σlmax = 1.461rad, ρfmax = 0.0395pixel−1.

Fig. 15. Experimental wrapped phase maps, (a) phase map related to stress field around a crack
in a transparent plate tested by digital transmission holographic interferometry, (b) phase map
obtained to mechanical deformation of a rough sample measured with digital Fresnel
holography.

The phase maps were unwrapped by the different algorithms of Section 3. The parameters
for the CPULSI algorithm were set to Nit = 300 and Err = 0.001. Considering the different
algorithms, the unwrapped results for the two phase maps are shown in Figs. 16 and 17
respectively.
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28729

Fig. 16. Unwrapped phases of Fig. 15(a) obtained from (a) Gold, (b) Quality, (c) Flynn, (d) Lp,
(e) PULSI, (f) CPULSI.

Figures 16 and 17 show that the proposed approach yields very good unwrapped phases,
the other approaches failing, and this result is in good agreement with what was demonstrated
with the numerical simulation.

Fig. 17. Unwrapped phases of Fig. 15(b) obtained from (a) Gold, (b) Quality, (c) Flynn, (d) Lp,
(e) PULSI, (f) CPULSI.

6. Conclusion and perspectives

This paper proposes an approach to unwrap successfully phase data corrupted by high speckle
decorrelation noise. A calibration method of the 1st spatial phase derivative coupled with an
iterative approach is presented. The algorithm is validated by realistic numerical simulations
in which the fringe density and noise standard deviation is progressively increased.
Especially, two sets of phase data with a certain kind of pattern density and structure diversity
were computed to produce non-Gaussian speckle decorrelation noise. These two sets were
processed by the proposed approach. In order to demonstrate the powerfulness of the
algorithm, comparison with five other established algorithms was carried out. This
comparison shows that the proposed approach exhibit better accuracy and faster computation
speed, whereas others may fail to unwrap. From a practical point of view, the proposed
algorithm was applied to phase data from digital holographic metrology. The experimental
results demonstrate its robustness for difficult phase maps.
Future work will consider the case of phase dislocations which may exist in phase maps.
Such phase dislocations can be due to lacks in spatial resolution or too rapid temporal phase
 Vol. 24, No. 25 | 12 Dec 2016 | OPTICS EXPRESS 28730

changes. At the moment, there do not exist any unwrapping approaches able to deal with such
kind of phase singularity.
7. Funding
National Natural Science Foundation of China (NSFC) (11362007, 11462009); French
National Agency for Research (ANR) (ANR-14-ASTR-0005-01).