Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Phase calibration unwrapping algorithm for phase data corrupted by strong decorrelation speckle noise

Open Access Open Access

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

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 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...,M1;j=0,1...N1),
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...,M2;j=0,1...N1)Δijx=0otherwiseΔijy=W(ψi(j+1)ψij)(i=0,1...,M1;j=0,1...N2)Δijy=0otherwise,
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.

 figure: Fig. 1

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.

Download Full Size | PDF

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 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=ΔijxotherwiseΔijy=sgn(Δijy)|Gy|if|Δijy|TyΔijy=Δijyotherwise,
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[(Δijx)2](E[Δijx])2Ty=E[(Δijy)2](E[Δijy])2,
and,
{Gx=1MNi=0i=M1j=0j=N1ΔijxGy=1MNi=0i=M1j=0j=N1Δijy.
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]:

s=i=0M2j=0N1|ϕ(i+1)jϕijΔijx|2+i=0M1j=0N2|ϕi(j+1)ϕijΔijy|2.
Equation (5) can be reduced to the discrete Poisson equation with Neumann boundary conditions [1]:
(ϕ(i+1)j2ϕij+ϕ(i1)j)+(ϕi(j+1)2ϕij+ϕi(j1))=ρij,
where
ρij=(ΔijxΔ(i1)jx)+(ΔijyΔi(j1)y).
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 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=ρ^ij2cos(πi/M)+2cos(π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+2πround(ϕij,kψij2π),
where round(…) is the rounding operator. This equation can be rewritten as
φij,kψij2π=round(ϕij,kψij2π).
Define the k-th phase error as:
Δψij,k=φij,kϕij,k.
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.
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).

 figure: Fig. 2

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

Download Full Size | PDF

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 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.

 figure: Fig. 3

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

Download Full Size | PDF

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.

 figure: Fig. 4

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.

Download Full Size | PDF

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

ρf,ij=12πKLp=iK/2i+K/2q=jL/2j+L/2(Δpqx)2+(Δpqy)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.

Tables Icon

Table 1. Standard deviation of noise and fringe density of the phase maps for the two sets.

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”.

 figure: Fig. 5

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.

Download Full Size | PDF

The next section discusses on the evaluation of the proposed unwrapping algorithm and comparison with other established approaches.

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).

 figure: Fig. 6

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

Download Full Size | PDF

 figure: Fig. 7

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).

Download Full Size | PDF

 figure: Fig. 8

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

Download Full Size | PDF

 figure: Fig. 9

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).

Download Full Size | PDF

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.

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.

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.

 figure: Fig. 10

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).

Download Full Size | PDF

 figure: Fig. 11

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).

Download Full Size | PDF

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).

Tables Icon

Table 2. Standard deviations of phase difference between unwrapped phase and original noisy phase (rad) for Set n°1

Tables Icon

Table 3. Standard deviations of phase difference between unwrapped phase and original noisy phase (rad) for Set n°2

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.

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.

 figure: Fig. 12

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

Download Full Size | PDF

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 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.

 figure: Fig. 13

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.

Download Full Size | PDF

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.

 figure: Fig. 14

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.

Download Full Size | PDF

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 results do not prevent from optimized computation, such as with GPU for example, which would be quite faster.

Tables Icon

Table 4. Average computation times of the different algorithms.

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.

 figure: Fig. 15

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.

Download Full Size | PDF

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.

 figure: Fig. 16

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

Download Full Size | PDF

 figure: Fig. 17

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

Download Full Size | PDF

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.

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 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).

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). [CrossRef]  

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). [CrossRef]  

4. P. Andris and I. Frollo, “Simple and accurate unwrapping phase of MR data,” Measurement 42(5), 737–741 (2009). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

7. F. Da and H. Huang, “A fast, accurate phase unwrapping method for wavelet-transform profilometry,” Opt. Commun. 285(4), 421–432 (2012). [CrossRef]  

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). [CrossRef]  

9. F. Chen and X. Su, “Phase-unwrapping algorithm for the measurement of 3D object,” Optik (Stuttg.) 123(24), 2272–2275 (2012). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

17. X. Su and W. Chen, “Reliability-guided phase unwrapping algorithm: a review,” Opt. Lasers Eng. 42(3), 245–261 (2004). [CrossRef]  

18. S. Yuqing, “Robust phase unwrapping by spinning iteration,” Opt. Express 15(13), 8059–8064 (2007). [CrossRef]   [PubMed]  

19. V. V. Volkov and Y. Zhu, “Deterministic phase unwrapping in the presence of noise,” Opt. Lett. 28(22), 2156–2158 (2003). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

21. Y. Guo, X. Chen, and T. Zhang, “Robust phase unwrapping algorithm based on least squares,” Opt. Lasers Eng. 63, 25–29 (2014). [CrossRef]  

22. M. Rivera, F. J. Hernandez-Lopez, and A. Gonzalez, “Phase unwrapping by accumulation of residual maps,” Opt. Lasers Eng. 64, 51–58 (2015). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

25. Q. Kemao, “Applications of windowed Fourier fringe analysis in optical measurement: a review,” Opt. Lasers Eng. 66, 67–73 (2015). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]   [PubMed]  

29. J. Weng and Y. Lo, “Novel rotation algorithm for phase unwrapping applications,” Opt. Express 20(15), 16838–16860 (2012). [CrossRef]  

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). [CrossRef]   [PubMed]  

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). [CrossRef]  

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). [CrossRef]   [PubMed]  

37. S. Montresor and P. Picart, “Quantitative appraisal for noise reduction in digital holographic phase imaging,” Opt. Express 24(13), 14322–14343 (2016). [CrossRef]   [PubMed]  

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). [CrossRef]  

39. Q. Kemao, “Windowed Fourier transform for fringe pattern analysis,” Appl. Opt. 43(13), 2695–2702 (2004). [CrossRef]   [PubMed]  

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (17)

Fig. 1
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.
Fig. 2
Fig. 2 Flow chart of the unwrapping algorithm; DCT means Discrete Cosine Transform, IDCT means Inverse Discrete Cosine Transform, <…> means average value.
Fig. 3
Fig. 3 Arrangement to produce speckle phase decorrelation in phase change due to surface deformation.
Fig. 4
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.
Fig. 5
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.
Fig. 6
Fig. 6 Phase maps n°10 of Set 1, (a) original unwrapped phase, (b) noisy wrapped phase (colorbars given in radian units).
Fig. 7
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).
Fig. 8
Fig. 8 Phase maps n°8 of Set 2, (a) original unwrapped phase, (b) noisy wrapped phase (colorbars given in radian units).
Fig. 9
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).
Fig. 10
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
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).
Fig. 12
Fig. 12 Standard deviations of the errors between the unwrapped phase and the noisy original phase: (a) for Set1, (b) for Set2.
Fig. 13
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.
Fig. 14
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.
Fig. 15
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.
Fig. 16
Fig. 16 Unwrapped phases of Fig. 15(a) obtained from (a) Gold, (b) Quality, (c) Flynn, (d) Lp, (e) PULSI, (f) CPULSI.
Fig. 17
Fig. 17 Unwrapped phases of Fig. 15(b) obtained from (a) Gold, (b) Quality, (c) Flynn, (d) Lp, (e) PULSI, (f) CPULSI.

Tables (4)

Tables Icon

Table 1 Standard deviation of noise and fringe density of the phase maps for the two sets.

Tables Icon

Table 2 Standard deviations of phase difference between unwrapped phase and original noisy phase (rad) for Set n°1

Tables Icon

Table 3 Standard deviations of phase difference between unwrapped phase and original noisy phase (rad) for Set n°2

Tables Icon

Table 4 Average computation times of the different algorithms.

Equations (14)

Equations on this page are rendered with MathJax. Learn more.

W( ϕ ij )= ψ ij ( i=0,1...,M1;j=0,1...N1 ),
Δ ij x =W( ψ ( i+1 )j ψ ij )( i=0,1...,M2;j=0,1...N1 ) Δ ij x =0otherwise Δ ij y =W( ψ i( j+1 ) ψ ij )( i=0,1...,M1;j=0,1...N2 ) Δ ij y =0otherwise,
Δ ij x =sgn( Δ ij x )| G x |if| Δ ij x | T x Δ ij x = Δ ij x otherwise Δ ij y =sgn( Δ ij y )| G y |if| Δ ij y | T y Δ ij y = Δ ij y otherwise,
{ T x = E[ ( Δ ij x ) 2 ] ( E[ Δ ij x ] ) 2 T y = E[ ( Δ ij y ) 2 ] ( E[ Δ ij y ] ) 2 ,
{ G x = 1 MN i=0 i=M1 j=0 j=N1 Δ ij x G y = 1 MN i=0 i=M1 j=0 j=N1 Δ ij y .
s= i=0 M2 j=0 N1 | ϕ ( i+1 )j ϕ ij Δ ij x | 2 + i=0 M1 j=0 N2 | ϕ i( j+1 ) ϕ ij Δ ij y | 2 .
( ϕ ( i+1 )j 2 ϕ ij + ϕ ( i1 )j )+( ϕ i( j+1 ) 2 ϕ ij + ϕ i( j1 ) )= ρ ij,
ρ ij =( Δ ij x Δ ( i1 )j x )+( Δ ij y Δ i( j1 ) y ).
ϕ ^ ij = ρ ^ ij 2cos( πi/M )+2cos( πj/N )4 ,
φ ij,k = ψ ij +2πround( ϕ ij,k ψ ij 2π ),
φ ij,k ψ ij 2π =round( ϕ ij,k ψ ij 2π ).
Δ ψ ij,k = φ ij,k ϕ ij,k .
ϕ ij,k+1 ϕ ij,k + ϕ ij,k+1 .
ρ f,ij = 1 2πKL p=iK/2 i+K/2 q=jL/2 j+L/2 ( Δ pq x ) 2 + ( Δ pq y ) 2 .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.