Advanced iterative algorithm for phase extraction: performance evaluation and enhancement

Yuchi Chen; Qian Kemao

doi:10.1364/OE.27.037634

1. Introduction

Phase-shifting interferometry has been widely used to precisely measure physical quantities such as shape, deformation and refractive index [1–3]. Phase-shifted fringe patterns can be represented as [4]

(1)$$I_{ij}^t = {A_{ij}} + {B_{ij}}\cos ({\varphi _j} + {\delta _i}),i = 1,2, \ldots ,M;j = 1,2, \ldots ,N,$$

where the subscript i represents the i^th fringe pattern with M as the total number of fringe patterns; the subscript j represents the j^th pixel with N as the total number of pixels; ${I_{ij}^{t}}$ is the theoretical fringe intensity; A_ij and B_ij represent the background intensity and fringe amplitude, respectively; φ_j is the phase information of the j^th pixel and is assumed to be independent of i; and δ_i is the phase-shift of the i^th frame with δ₁= 0, and is assumed to be independent of j. Among these variables, the phase information φ_j is directly related to the measured quantities [1–3] and thus need to be extracted from fringe patterns. In practice, as noise is unavoidable, the measured fringe patterns are normally presented as

(2)$${I_{ij}} = I_{ij}^t + {n_{ij}},$$

where n_ij is the noise of the j^th pixel at i^th fringe pattern.

Due to the extreme importance of phase extraction, many phase-shifting algorithms have been developed. Based on the prior knowledge about the phase-shifts, they can be classified into two categories [4]:

(i) Algorithms with known phase-shifts: the least square algorithm (LSA) [1], the N + 1 bucket algorithm [5] and the windowed digital Fourier transform algorithms [6] are such examples. Performance evaluation of these algorithms is relatively easy and has been intensively carried out [6–9];
(ii) Algorithms with unknown phase-shifts: the Carre’s algorithm [10], the Hariharan algorithm [11], the advanced iterative algorithm (AIA) [12], the algorithms based on principal component analysis [13], Euclidean matrix norm [14] and Gram–Schmidt orthonormalization [15] belong to this category. Naturally, the unknown phase-shifts bring uncertainty and challenge for performance evaluation.

Among all these algorithms, AIA is of particular interest to us. AIA iteratively estimates the phase and phase-shifts until the stopping criterion is satisfied. It offers the following merits: (i) the algorithmic simplicity and elegance, (ii) the flexibility of using arbitrary phase-shifts, and frame numbers (≥3), and (iii) the good performance demonstrated in the literature and our own testing. However, we also noticed that there is a lack of comprehensive algorithm analysis and performance evaluation with respect to (w.r.t) different phase distributions, phase-shifts, frame numbers, and noise, which leaves an obstacle for using AIA in industrial applications. Specifically, the following aspects of the AIA performance are unclear and should be answered: (i) Is AIA applicable to any set of phase-shifted fringe patterns and successfully converge? (ii) How is the performance affected by different phase-shifts, frame numbers and noise levels, and can the performance be further improved? (iii) Can the accuracy of AIA be predicted? These questions motivate us to investigate the AIA’s performance thoroughly, through which we can answer all these questions and propose various improvements accordingly. Finally, by integrating all our findings and solutions, an enhanced AIA (eAIA) is proposed.

The rest of the paper is arranged as follows. In Section 2, the AIA is briefly described and analyzed. Condition numbers are introduced to characterize the least squares matrices used in AIA. In Section 3, simulations of fringe patterns to evaluate AIA are set up. The requirement on fringe density for the effective use of AIA is also obtained and suggested. In Section 4, the performance of AIA w.r.t different variables is thoroughly evaluated, resulting in an integrated general error model. In Section 5, individual improvements w.r.t. different variables in AIA are proposed. In Section 6, practical error estimation methods are proposed to predict the performance of AIA for a given set of phase-shifted fringe patterns. In Section 7, the eAIA is finally proposed and verified in an experiment. Section 8 concludes the paper.

2. Advanced iterative algorithm (AIA) and its analysis

In this section, the AIA is described and analyzed. For clarity, the representations and symbols largely follow the original AIA in [12].

2.1 Algorithm description

AIA iteratively estimates phase and phase-shifts, until a convergence condition is satisfied. First, the phase is estimated by treating the phase-shifts as known; and both fringe background and amplitude are assumed to be constant across frames. Equation (1) is then rewritten as

(3)$$I_{ij}^t = {a_j} + {b_j}\cos {\delta _i} + {c_j}\sin {\delta _i},$$

where ${a_j} = {A_{ij}}$; ${b_j} = {B_{ij}}\cos {\varphi _j}$ and ${c_j} ={-} {B_{ij}}\sin {\varphi _j}$. The following loss function between the measured intensity I_ij and the theoretical model ${I_{ij}^{t}}$ should be minimized in the least squares sense:

(4)$${S_j} = \sum\nolimits_{i = 1}^M {(I_{ij}^t - } {I_{ij}}{)^2} = \sum\nolimits_{i = 1}^M {({a_j} + {b_j}\cos {\delta _i} + {c_j}\sin {\delta _i} - } {I_{ij}}{)^2}.$$

Letting $\frac{{\partial {S_j}}}{{\partial {a_j}}} = \frac{{\partial {S_j}}}{{\partial {b_j}}} = \frac{{\partial {S_j}}}{{\partial {c_j}}} = 0$ yields:

(5)$${\textbf{A}_\textrm{p}}{\textbf{X}_{\textrm{p},j}} = {\textbf{B}_{\textrm{p},j}},$$

where

(6)$${\textbf{A}_\textrm{p}} = \left[ {\begin{array}{ccc} M&{\sum\nolimits_{i = 1}^M {\cos {\delta_i}} }&{\sum\nolimits_{i = 1}^M {\sin {\delta_i}} }\\ {\sum\nolimits_{i = 1}^M {\cos {\delta_i}} }&{\sum\nolimits_{i = 1}^M {{{\cos }^2}{\delta_i}} }&{\sum\nolimits_{i = 1}^M {\cos {\delta_i}\sin {\delta_i}} }\\ {\sum\nolimits_{i = 1}^M {\sin {\delta_i}} }&{\sum\nolimits_{i = 1}^M {\cos {\delta_i}\sin {\delta_i}} }&{\sum\nolimits_{i = 1}^M {{{\sin }^2}{\delta_i}} } \end{array}} \right],$$

${\textbf{X}_{\textrm{p},j}} = {\left[ {\begin{array}{ccc} {{a_j}}&{{b_j}}&{{c_j}} \end{array}} \right]^T}$, and ${\textbf{B}_{\textrm{p},j}} = {\left[ {\begin{array}{ccc} {\sum\nolimits_{i = 1}^M {{I_{ij}}} }&{\sum\nolimits_{i = 1}^M {{I_{ij}}\cos {\delta_i}} }&{\sum\nolimits_{i = 1}^M {{I_{ij}}\sin {\delta_i}} } \end{array}} \right]^T}$. It is clear that A_p is a function of phase-shifts. If A_p is nonsingular, we have

(7)$${\textbf{X}_{\textrm{p},j}} = \textbf{A}_\textrm{p}^{ - 1}{\textbf{B}_{\textrm{p},j}},$$

where the subscript − 1 indicates the matrix inverse. Subsequently, the phase can be extracted as

(8)$${\varphi _j} = {\tan ^{ - 1}}\left( { - \frac{{{c_j}}}{{{b_j}}}} \right).$$

Similarly, the phase-shifts are estimated by treating the phases of all pixels as known and both fringe background and amplitude are assumed to be constant across pixels in each frame. Equation (1) is then rewritten as

(9)$$I_{ij}^t = a_i^{\prime} + b_i^{\prime}\cos {\varphi _j} + c_i^{\prime}\sin {\varphi _j},$$

where $a_i^{\prime} = {A_{ij}}$; $b_i^{\prime} = {B_{ij}}\cos {\delta _i}$ and $c_i^{\prime} ={-} {B_{ij}}\sin {\delta _i}$. The following loss function between the measured intensity I_ij and the theoretical model ${I_{ij}^{t}}$ should be minimized in the least squares sense:

(10)$$S_i^{\prime} = \sum\nolimits_{j = 1}^N {(I_{ij}^t - } {I_{ij}}{)^2} = \sum\nolimits_{j = 1}^N {(a_i^{\prime} + b_i^{\prime}\cos {\varphi _j} + c_i^{\prime}\sin {\varphi _j} - } {I_{ij}}{)^2}.$$

Letting $\frac{{\partial S_i^{\prime}}}{{\partial a_i^{\prime}}} = \frac{{\partial S_i^{\prime}}}{{\partial b_i^{\prime}}} = \frac{{\partial S_i^{\prime}}}{{\partial c_i^{\prime}}} = 0$ yields:

(11)$${\textbf{A}_{\textrm{ps}}}{\textbf{X}_{\textrm{ps},i}} = {\textbf{B}_{\textrm{ps},i}},$$

where

(12)$${\textbf{A}_{\textrm{ps}}} = \left[ {\begin{array}{ccc} N&{\sum\nolimits_{j = 1}^N {\cos {\varphi_j}} }&{\sum\nolimits_{j = 1}^N {\sin {\varphi_j}} }\\ {\sum\nolimits_{j = 1}^N {\cos {\varphi_j}} }&{\sum\nolimits_{j = 1}^N {{{\cos }^2}{\varphi_j}} }&{\sum\nolimits_{j = 1}^N {\cos {\varphi_j}\sin {\varphi_j}} }\\ {\sum\nolimits_{j = 1}^N {\sin {\varphi_j}} }&{\sum\nolimits_{j = 1}^N {\cos {\varphi_j}\sin {\varphi_j}} }&{\sum\nolimits_{j = 1}^N {{{\sin }^2}{\varphi_j}} } \end{array}} \right],$$

${\textbf{X}_{\textrm{ps},i}} = {\left[ {\begin{array}{ccc} {a_i^{\prime}}&{b_i^{\prime}}&{c_i^{\prime}} \end{array}} \right]^T}$ and ${\textbf{B}_{\textrm{ps},i}} = {\left[ {\begin{array}{ccc} {\sum\nolimits_{j = 1}^N {{I_{ij}}} }&{\sum\nolimits_{j = 1}^N {{I_{ij}}\cos {\varphi_j}} }&{\sum\nolimits_{j = 1}^N {{I_{ij}}} } \end{array}\sin {\varphi_j}} \right]^T}.$ The A_ps is a function of phase distribution. If A_ps is nonsingular, we have

(13)$${\textbf{X}_{\textrm{ps},i}} = \textbf{A}_{\textrm{ps}}^{ - 1}{\textbf{B}_{\textrm{ps},i}},$$

from which, the phase-shifts are calculated,

(14)$${\delta _i} = {\tan ^{ - 1}}\left( { - \frac{{c_i^{\prime}}}{{b_i^{\prime}}}} \right).$$

The above phase calculation and phase-shifts calculation are iteratively executed. AIA starts with random guessed phase-shifts and stops when the following condition is satisfied for all i:

(15)$$|(\delta _i^k - \delta _1^k) - (\delta _i^{k - 1} - \delta _1^{k - 1})|< {\varepsilon} ,$$

where $\varepsilon $ is a pre-set error tolerance and the superscript k is the number of iteration.

2.2 Condition numbers of least squares matrices

As can be seen from Eqs. (7) and (13), the performance of AIA directly relies on the two least squares matrices A_p and A_ps, thus phase distribution and phase-shifts. AIA works poorly when A_p and A_ps are singular or nearly singular. Since condition number is a good indicator of matrix singularity [16], we propose to adopt it to qualify these matrices. Given a linear equation system Hx = F, the condition number of matrix H is calculated as follows [17],

(16)$$\kappa (\textbf{H} )= \frac{{{{\sigma}_{\max }}(\textbf{H})}}{{{{\sigma}_{\min }}(\textbf{H})}},$$

where ${{\sigma}_{\min }}(\textbf{H} )$ and ${{\sigma}_{\max }}(\textbf{H} )$ are the minimal and maximal singular values of H, respectively.

Consider matrix A_ps first. For any vector x with real components, the dot product of A_psx and x is

(17)$$\left\langle {{\textbf{A}_{\textrm{ps}}}\textbf{x},\textbf{x}} \right\rangle = {\textbf{x}^T}{\textbf{A}_{\textrm{ps}}}\textbf{x} = \sum\nolimits_{j = 1}^N {{{({x_1} + \cos {\varphi _j}{x_2} + \sin {\varphi _j}{x_3})}^2}} \ge 0,$$

which implies that A_ps is positive semidefinite [17]. Since A_ps is a symmetrical positive semidefinite real matrix, based on the property of the Rayleigh quotient [18] and the fact that eigenvalues of A_ps equal to its singular values [17], for any vector x with nonzero real component, we have,

(18)$${{\sigma}_{\min }}({{\textbf{A}_{\textrm{ps}}}} )\le \frac{{{\textbf{x}^T}{\textbf{A}_{\textrm{ps}}}\textbf{x}}}{{{\textbf{x}^T}\textbf{x}}} \le {{\sigma}_{\max }}({{\textbf{A}_{\textrm{ps}}}} ).$$

Substituting different values of x, e.g. ${\left[ {\begin{array}{ccc} 1&0&0 \end{array}} \right]^T}$, ${\left[ {\begin{array}{ccc} 0&1&0 \end{array}} \right]^T}$ and ${\left[ {\begin{array}{ccc} 0&0&1 \end{array}} \right]^T}$, into Eq. (18) yields

(19)$${{\sigma}_{\max }}({{\textbf{A}_{\textrm{ps}}}} )\ge N,$$

(20)$${{\sigma}_{\min }}({{\textbf{A}_{\textrm{ps}}}} )\le \min \left( {\sum\nolimits_{j = 1}^N {{{\cos }^2}{\varphi_j}} ,\sum\nolimits_{j = 1}^N {{{\sin }^2}{\varphi_j}} } \right).$$

Thus following Eq. (16), we have

(21)$${\kappa}({{\textbf{A}_{\textrm{ps}}}} )\ge \frac{N}{{\min \left( {\sum\nolimits_{j = 1}^N {{{\cos }^2}{\varphi_j}} ,\sum\nolimits_{j = 1}^N {{{\sin }^2}{\varphi_j}} } \right)}}.$$

Based on Pythagorean trigonometric identity,

(22)$${\cos ^2}\varphi + {\sin ^2}\varphi = 1,$$

it is not difficult to find that

(23)$$\min \left( {\sum\nolimits_{j = 1}^N {{{\cos }^2}{\varphi_j}} ,\sum\nolimits_{j = 1}^N {{{\sin }^2}{\varphi_j}} } \right) \le \frac{N}{2}.$$

Hence the inequality Eq. (21) becomes

(24)$${\kappa}({{\textbf{A}_{\textrm{ps}}}} )\ge \textrm{2},$$

where ${\kappa}({{\textbf{A}_{\textrm{ps}}}} ) = 2$ is achievable when ${{\sigma}_{\min }}({{\textbf{A}_{\textrm{ps}}}} )= \sum\nolimits_{j = 1}^N {{{\cos }^2}{\varphi _j}} = \sum\nolimits_{j = 1}^N {{{\sin }^2}{\varphi _j}} = \frac{N}{2}$, i.e., the phase evenly distributed among 2π. In the meantime, A_ps could be singular in special cases, e.g. φ_j = 0, for j = 1, 2, …, N. Thus, we have 2≤κ(A_ps)<∞. Similarly, we also have 2≤κ(A_p)<∞, where ${\kappa}({{\textbf{A}_\textrm{p}}} ) = 2$ is achievable when δ_i = (i-1)×2π/M.

3. Fringe pattern simulation and fringe density requirement for AIA

In this section, the fringe pattern simulation strategies which are used throughout the paper are introduced first. Based on simulations, a fringe density requirement for using AIA is established.

3.1 Fringe pattern simulation

Our main tool for AIA evaluation in this paper is simulation, which provides us the ground true values for comparison and enables us to run tests thousands of times easily. Nevertheless, an experimental example will be provided at the end of the paper.

The fringe patterns are simulated according to Eq. (2) with the following settings:

(i) M = 3 ∼ 40;
(ii) N = N_x (horizontal) × N_y (vertical) = 256×256;
(iii) the background intensity A and fringe amplitude B are constant across both frames and pixels;
(iv) three typical phase distributions, linear for straight fringe patterns (φ¹), quadratic for circular fringe patterns (φ²), and a “peaks” function for complex fringe patterns (φ³), are used as represented below, $(25)$${\varphi ^1}(x,y) = \frac{{2{N_p}{\pi}}}{{{N_x}}}x,$$$ $(26)$${\varphi ^2}(x,y) = \frac{{2{N_p}{\pi}}}{{{{\left( {\frac{{{N_x}}}{2}} \right)}^2} + {{\left( {\frac{{{N_y}}}{2}} \right)}^2}}}{r^2},$$$ $(27)$${\varphi ^3}(x,y) = {N_p} \times \textrm{peaks}_{\textrm{s}}({N_x},{N_y}) + {\omega _c}x$$$ where N_p is number of fringes over ROI, ${r^2} = {\left( {x - \frac{{{N_x}}}{2}} \right)^2} + {\left( {y - \frac{{{N_y}}}{2}} \right)^2}$, “peaks_s” uses the MATLAB function “peaks”, linearly scaled it into the range of [0, 2π] and ω_c represents the carrier frequency;
(v) phase-shifts are either randomly drawn from [0, 2π] or specially set as δ_i = (i-1)×2π/M so that κ(A_ps) is minimal;
(vi) additive white noises have a zero mean and different standard deviations relative to the fringe amplitude B, i.e., σ/B = 0% ∼ 100% [19].

As examples, a noiseless straight fringe pattern with N_p = 1, a circular fringe pattern with N_p = 2 and additive noise with σ/B = 10% and a complex fringe pattern with N_p = 3, ω_c = 0 and additive noise with σ/B = 20% are shown in Figs. 1(a)–1(c), respectively.

Fig. 1. Typical fringe patterns used for simulation. (a) A noiseless straight fringe pattern with N_p = 1. (b) A circular fringe pattern with N_p= 2 and additive noise with σ/B = 10%. (c) A complex fringe pattern with N_p = 3 and ω_c = 0 and additive white noise with σ/B = 20%

Abstract

1. Introduction

2. Advanced iterative algorithm (AIA) and its analysis

2.1 Algorithm description

2.2 Condition numbers of least squares matrices

3. Fringe pattern simulation and fringe density requirement for AIA

3.1 Fringe pattern simulation

3.2 Fringe density requirement

4. Performance evaluation of AIA

4.1 Influence of phase-shifts

4.2 Influence of frame number

4.3 Influence of noise

4.4 The integrated general error model of AIA

5. Improvements of AIA

5.1 Improvement by phase-shifts control

5.2 Improvement by frame number control

5.3 Improvement by fringe denoising

6. Practical error estimation

6.1 Estimation of phase-shifts

6.2 Estimation of fringe amplitude

6.3 Estimation of noise standard deviation

6.4 Practical overall error estimation

7. Enhanced AIA (eAIA)

7.1 Structure of eAIA

7.2 Experimental verification of eAIA

8. Conclusions

Funding

References

Cited By

Figures (15)

Equations (41)

Optics Express