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

General iterative algorithm for phase-extraction from fringe patterns with random phase-shifts, intensity harmonics and non-uniform phase-shift distribution

Open Access Open Access

Abstract

Advanced iterative algorithm (AIA) is a flexible and effective phase-shifting algorithm (PSA) which can extract phase from fringe patterns with random unknown phase-shifts, making it attractive in the scenarios where phase-shifts are unknown or not accurate. However, accuracy of AIA degrades when intensity harmonics and/or phase-shift non-uniformity are presented. To solve this problem, multiple PSAs have been proposed, but they restrict their fringe model in one way or another, and thus sacrifice the immunity to certain error source(s). In this paper, a general iterative algorithm (GIA) which adopts a most general fringe model is proposed. In GIA, the many unknowns in the fringe pattern model are divided into three groups including: (i) the fringe amplitudes, (ii) the phase and (iii) the phase-shifts related parameters, and alternatively optimized through univariate search technique group by group to improve accuracy and convergence. The Levenberg-Marquart method is used for the optimization of each group of unknowns due to its excellent accuracy and robustness. GIA is shown to have better accuracies than all of its relevant competitors through both a large number of simulations as well as real experiments with a Fizeau interferometer.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Phase-shifting interferometry (PSI) as a simple and accurate technique has been widely used in optical metrology [14]. In PSI, a series of phase-shifted fringe patterns are generated and captured, from which a phase distribution relating to a measured quantity such as shape, deformation or refractive index [57] is extracted. In practice, phase extraction is affected by a variety of error sources [3] including (i) phase-shift error, (ii) random additive intensity noise, (iii) intensity harmonics and (iv) non-uniform phase-shift distribution. Among them, the phase-shift error due to mis-calibration or im-perfection of the phase-shifter is commonly encountered and has been intensively studied [815]. Meanwhile, the random noise induced by components of the optical system especially the detector is unavoidable [16]. Besides, the intensity harmonics due to the non-linearity of the detector or multi-beam reflection will affect the accuracy of high precision measurement [17,18]. Although intensity harmonics induced by detector non-linearity could be minimized with high-end hardware, it still affects interferometer with normal detector and/or multiple reflections, such as in Fizeau interferometer. Moreover, the phase-shift non-uniformity should not be ignored in the practical measurement with vibrations caused by rigid and non-rigid movement between the reference and specimen or air turbulence [19,20].

To accurately extract the phase with the presence of these error sources, many phase-shifting algorithms (PSA) have been developed. Most PSAs are with known phase-shifts [17]. The influence of the phase-shift error, random noise, intensity harmonics and non-uniform phase-shifts to these PSAs could be evaluated by analytical expansion method [8], characteristic polynomials [9], Fourier description method [10] or frequency transfer function [11].

Another type of PSAs, which is less widely used, adopts unknown phase-shifts. Typical examples include advanced iterative algorithm (AIA) based on an alternative linear least square fitting (LSF) framework [12], principal component analysis (PCA) algorithms [13] which are smart and popular with good accuracy and fast computation, and algorithms based on Euclidean matrix norm [14] and Gram–Schmidt orthonormalization [15]. These PSAs are naturally immune to the phase-shift errors. Among them, AIA is known for simplicity of the algorithm, flexibility of using arbitrary phase-shifts and frame numbers (≥3) and great immunity to random noise [21]. Therefore, a few AIA-based or related PSAs have subsequently been developed to suppress the influence of intensity harmonics and/or phase-shift non-uniformity, which can be further classified into three groups:

  • (i) To deal with intensity harmonics, Xu et al. proposed to directly incorporate them into the original alternative linear LSF framework of AIA [22], and Hoang et al. chose to combine LSF with exhaustive search [23]. We will show later that these two algorithms are sensitive to random noise and the phase-shift non-uniformity. In addition, Hoang’s algorithm is slow.
  • (ii) To deal with non-uniform phase-shift distribution, most works treat it as a tilt, i.e., a spatially linear phase-shift model. Then, Chen et al. replaced the linear LSF in AIA by non-linear LSF [24]; Duan et al. proposed a phase-tilt iteration (PTI) algorithm to extract phase-shift tilts [25]; Xu et al. and Liu et al. divided fringe patterns into small regions, and utilized the phase-shifts extracted from the small regions by AIA to determine the phase-shift tilts through linear fitting [26] and plane fitting [27], respectively. These algorithms did not consider intensity harmonics or spatially non-linear phase-shifts.
  • (iii) To deal with both intensity harmonics and non-uniform phase-shifts, Deck proposed a model-based phase-shifting interferometry (MPSI) [28] which estimates amplitude of harmonics through LSF and extracts phase-shift tilts through first order linearization in the deviation. MPSI does not consider the intensity harmonics caused by detector’s non-linearity and spatially non-linear phase-shifts. It also requires the phase-shifts to be pre-measured to initialize the algorithm and prefers phase-shifts who are evenly distributed among 2π. Thus, when the phase-shifts are arbitrary or not known in advance, as we are discussing in this paper, the performance of MPSI will decline.

Based on the above review and analysis, we are motivated to propose a PSA which is general enough to deal with all the mentioned error sources and unknown phase-shifts, simple enough to be easily realized, yet accurate and robust enough to be used in optical interferometry. To achieve this goal, we develop a general iterative algorithm (GIA) which adopts a most general fringe model to incorporate all the mentioned errors, and naturally, leads to many unknown parameters. Subsequently, GIA divides all the unknowns into three groups: (i) the background intensity and fringe amplitudes, (ii) the phase, and (iii) the phase-shift related parameters. These three groups of unknowns are alternatively optimized through univariate search technique [29] for better accuracy and convergence [30]. The Levenberg-Marquart method [31] is used for optimization of each group due to its excellent accuracy and robustness. In a nutshell, GIA only applies the Levenberg-Marquart method multiple times to reach the final convergence, making the algorithm simple in structure and easy for implementation. The good convergence and accuracy performance of GIA will be demonstrated by comparing with its most relevant competitors through a large number of simulations and real experiments from a Fizeau interferometer.

The rest of the paper is arranged as follows. In Sec. 2, a general fringe model is introduced first, and the most relevant works are briefly explained. The details of our proposed GIA are described in Sec. 3, followed by the simulation and experimental verification given in Sec. 4 and Sec. 5, respectively. This paper is concluded in Sec. 6.

2. General fringe pattern model and the related iterative algorithms

In this section, a general fringe pattern model which incorporates various error sources is provided, based on which, our research problem is stated. The most relevant algorithms are then reviewed. While keeping the introduction concise, sufficient details are given to differentiate and compare among the existing algorithms and our proposed GIA.

2.1 Fringe pattern model

Random noise is unavoidable in real measurement. Thus, a phase-shifted fringe pattern is first written as a combination of a noiseless part and a noise part as follows [32]

$$I({i;{j_{x,y}}} )= {I^t}({i;{j_{x,y}}} )+ n({i;{j_{x,y}}} ),$$
where I and It are the measured fringe intensity and the noiseless fringe intensity, respectively; n is the random noise which is modeled as additive white Gaussian noise (AWGN) in this paper; i=0,1, …, F-1 is the frame index with F as the total number of frames and jx,y=0, 1, …, N-1 is the pixel index with N as the total number of pixels. The index jx,y relates to the spatial location of the pixel (x, y) through jx,y=xNy + y where x=0, 1, …, Nx-1 and y=0, 1, …, Ny-1. Meanwhile, in later sections when the spatial location of the pixel is not critical, the subscripts x and y will be omitted.

Next, the noiseless phase-shifted fringe pattern with intensity harmonics is modeled as

$${I^t}({i;j} )= A({i;j} )+ \sum\nolimits_{k = 1}^P {{B_k}({i;j} )\cos [{k\varphi (j )+ k\delta ({i;j} )} ]} ,$$
where A(i; j) is the background intensity; k=1, 2, …, P is the order of harmonics with P as the highest order; Bk(i; j) is the fringe amplitude of the k-th order harmonic, and when P=1, we also write B1 as B; φ(j) is a simplified notation of φ(·; j) which indicates the phase information of the j-th pixel and is independent of the frame i. Such simplified notations will be frequently used for other parameters to maximize the similarity to the existing literature, when no confusion is raised; δ(i; j) is the phase-shift of the i-th frame at the j-th pixel.

Last, the phase-shift with a possible non-uniform distribution is written as a Maclaurin polynomial due to its simplicity and good accuracy in function approximations [33],

$$\delta ({i;j} )= \sum\nolimits_{u = 0}^Q {\sum\nolimits_{v = 0}^u {\frac{{{\alpha _{uv}}(i ){x^u}{y^{u - v}}}}{{{{({{N_x} - 1} )}^u}{{({{N_y} - 1} )}^{u - v}}}}} } ,$$
where αuv(i) are the coefficients of the Maclaurin polynomial; and u and u-v are the polynomial orders of x and y, respectively, with Q as the highest order. Since phase-shifts are completely unknown, in the first frame, we consider δ(0; j) = 0 in all the algorithms. It is worth noticing that, in the other frames, the α00(i) contains two parts, one is the true phase-shift introduced by phase-shifter, and the other is the zero-order term in the Maclaurin polynomial caused by undesired phase-shift error. However, as we extract the phase with unknown phase-shifts, the differentiation of these two parts is unnecessary.

The above fringe model is general and thus flexible to use. When the intensity harmonics are absent, we can simply set P=1; when the phase-shift non-uniformity is absent, we can set Q=0; the original AIA is developed considering both P=1 and Q=0. Generally, iteration-based algorithms minimize the following least squares error between the noiseless fringe intensity and measured intensity,

$$E = \sum\nolimits_{i,j} {{{[{{I^t}({i;j} )- I({i;j} )} ]}^2}} ,$$
where the summation is across both i and j, i.e. for all the pixels and frames. In some circumstances, the summation is pixel-wise and only across i,
$${E_P}(j )= \sum\nolimits_i {{{[{{I^t}({i;j} )- I({i;j} )} ]}^2}} ,$$
or frame-wise and only across j,
$${E_F}(i )= \sum\nolimits_j {{{[{{I^t}({i;j} )- I({i;j} )} ]}^2}} .$$

We now state our problem as to extract the phase (i) with high accuracy, (ii) without knowing phase-shift values, and (iii) with the presence of noise, intensity harmonics and non-unform phase-shifts in fringe pattern as modeled by Eqs. (13).

2.2 Advanced iterative algorithm (AIA)

AIA is introduced first for easy understanding of the other algorithms[12]. By setting P=1 and Q=0 in Eq. (2), AIA deals with the following fringe model,

$${I^t}({i;j} )= A({i;j} )+ B({i;j} )\cos [{\varphi (j )+ \delta (i )} ].$$

AIA iteratively estimates phase φ and phase-shifts δ, until a convergence condition is satisfied.

The first step estimates the phase in a pixel-wise manner. By treating the phase-shift δ(i) as known and assuming both the background intensity and fringe amplitude are constant across frames, Eq. (7) can be rewritten as

$${I^t}({i;j} )= a(j )+ b(j )\cos [{\delta (i )} ]+ c(j )\sin [{\delta (i )} ],$$
where a(j)=A(j); b(j)=B(j)cos[φ(j)] and c(j)=−B(j)sin[φ(j)]. For each fixed j, different i indices give multiple equations which are linear with respect to the unknowns a(j), b(j) and c(j). This linear equation system can be easily solved by LSF with an inverse of a $3 \times 3$ matrix. This technique will also be used in other algorithms and is called linear LSF with matrix inverse (LSF_MI) for convenience. Afterwards, the phase φ(j) can be calculated from b(j) and c(j).

The second step estimates phase-shifts in a frame-wise manner. By treating the phases of all pixels as known and assuming both background intensity and fringe amplitude are constant across pixels in each frame, Eq. (7) is rewritten as

$${I^t}({i;j} )= a^{\prime}(i )+ b^{\prime}(i )\cos [{\varphi (j )} ]+ c^{\prime}(i )\sin [{\varphi (j )} ],$$
where a′(i)=A(i); b′(i)=B(i)cos[δ(i)] and c′(i)= −B(i)sin([δ(i)]. Again, a′(i), b′(i) and c′(i) can be estimated by LSF_MI. Then, the phase-shifts δ(i) can be calculated from the b′(i) and c′(i).

However, neither the phase nor the phase-shifts are known at the beginning. AIA initializes phase-shift δ(i) as random values casted within [0, 2π] or as pre-defined values based on prior knowledge, and then iteratively execute the phase estimation step and the phase-shift estimation step. The stopping criteria is either the absolute phase-shift changes of all the frames are sufficiently small, or the maximum iteration number is reached, i.e.,

$$|{[{{\delta^m}(i )- {\delta^m}(0 )} ]- [{{\delta^{m - 1}}(i )- {\delta^{m - 1}}(0 )} ]} |< \varepsilon \textrm{, for all the frames,}$$
or
$$m \ge M\textrm{,}$$
where |·|takes the absolute value, ε is a pre-set error tolerance, the superscript m indicates the number of iterations, and M is the pre-set maximum iteration number. More details of AIA can be found in Ref. [12]. Setting ε=1×10−4 rad is usually suitable for high precision phase extraction [12,21,34] and can often be achieved with a few iterations. However, we set M=200 to secure possible slow convergence.

AIA performs well with noisy fringe patterns [21]. However, AIA does not consider intensity harmonics and non-uniform phase-shift distribution, resulting in worse performance if such errors exist. Thus, to improve the accuracy of the extracted phase, different PSAs has been proposed.

Note that we have proposed an enhanced AIA (eAIA) which improves the immunity of AIA to random noise through phase-shift control, frame number control and noise filtering. These factors are applicable to other methods as well and thus have been considered in our simulation conditions in Secs. 4.2-4.4. They are further discussed regarding convergence and accuracy for GIA in Sec. 4.5.

2.3 Algorithm for intensity harmonics rejection

Based on AIA, Xu et al. and Hoang et al. have proposed their intensity harmonic immune algorithms [22,23]. These two algorithms share the same fringe pattern model with P>1 and Q=0 in Eq. (2), i.e.,

$${I^t}({i;j} )= A({i;j} )+ \sum\nolimits_{k = 1}^P {{B_k}({i;j} )\cos [{k\varphi (j )+ k\delta (i )} ]} .$$

2.3.1 Xu’s algorithm [22]

Xu’s algorithm directly incorporates the intensity harmonics into the original alternative linear LSF framework of AIA. In the first step, by assuming that phase-shifts are known and the background intensity and fringe amplitude of each order of harmonics have no frame to frame variations, Eq. (12) can be rewritten as

$${I^t}({i;j} )= \sum\nolimits_{q = 0}^{2P} {{X_q}(j ){S_q}} (i ),$$
where X0(j)=A(j), X2k-1(j)=Bk(j)cos[(j)], X2k(j)=-Bk(j)sin[(j)], S0(i) = 1, S2k-1(i)=cos[(i)], and S2k(i)=sin[(i)]. The unknowns Xq(j) can be estimated by LSF_MI where the matrix has a size of (2P+1)×(2P+1) and its entry at s-th row and t-th column is $\mathop \sum \nolimits_{i = 0}^{F - 1} {S_s}(i ){S_t}(i )$. Subsequently, the phase can be calculated.

In the second step, by assuming that phase is known and the background intensity and fringe amplitude of each order of harmonics have no pixel to pixel variations, Eq. (12) can be rewritten as

$${I^t}({i;j} )= \sum\nolimits_{q = 0}^{2P} {{{X^{\prime}}_q}(i ){{S^{\prime}}_q}} (j ),$$
where X′0 (i)=A(i), X′2k-1(i)=Bk(i)cos[(i)], X′2k(i)=-Bk(i)sin[(i)], S′0 (j) = 1, S′2k-1(j)=cos[(j)], and S′2k(j)=sin[(j)]. Similarly, the unknowns X′q(i) can be estimated by LSF_MI where the matrix has a size of (2P+1)×(2P+1) and its entry at s-th row and t-th column is $\mathop \sum \nolimits_{j = 0}^{N - 1} S_s^{\prime}(j )S_t^{\prime}(j )$, and subsequently, the phase-shifts can be calculated.

Xu’s algorithm shares the same iterative execution, initialization and stopping criteria as AIA. Typically, it is sufficient to set P=3 for Xu’s algorithm, although a larger value can also be set [22].

By comparing Eqs. (13) and (14) with Eqs. (8) and (9), it can be seen clearly that Xu’s algorithm is a direct generalization of AIA, i.e., setting P=1 in Xu’s algorithm returns to the original AIA. However, this generalization naturally requires the inverse of two (2P+1)×(2P+1) least squares matrices which are greater than the corresponding 3×3 matrices used in AIA. Except when the phase-shifts are evenly distributed among 2π, the large size matrices will lead to larger condition numbers, resulting in worse immunity to the random noise [21,30]. Moreover, in a rare condition that the phase-shifts are evenly distributed as δ(i)=ilπ/k0 with k0$ {\mathbb{Z}}$+, k0P and l$ {\mathbb{Z}}$+, the entire (2k0)-th row becomes zero, making least squares matrix not inversible and hence failing the algorithm.

2.3.2 Hoang’s algorithm [23]

Hoang’s algorithm involves exhaustive search in the framework of AIA to extract the phase and phase-shifts. In the first step, assuming that the background intensity and fringe amplitude of each order of harmonics have no pixel to pixel variations, the fringe pattern in Eq. (12) can be represented as

$${I^t}({i;j} )= A(i )+ \sum\nolimits_{k = 1}^P {{B_k}(i )\cos [{k\varphi (j )+ k\delta (i )} ]} .$$

With the known phase and phase-shifts, the background intensity and all the fringe amplitudes can be estimated frame-wisely through LSF_MI with a (P+1)×(P+1) least squares matrix whose entry at s-th row and t-th column is $\mathop \sum \nolimits_{j = 0}^{N - 1} \{{\cos [{s\varphi (j )+ s\delta (i )} ]\cos [{t\varphi (j )+ t\delta (i )} ]} \}$. Next, the phase-shift in each frame is updated by minimizing the least squares error shown in Eq. (6) through exhaustive search (LSF_ES), using the estimated background intensity and fringe amplitudes, and the existing phase. It is worth noticing that the LSF here is non-linear, and thus cannot be solved by LSF_MI. As suggested by the authors [23], the phase-shift of each frame is estimated by exhaustive search within a range of π rads around the old phase-shift value and with a step of 0.0001 rad.

In the second step, assuming that the background intensity and fringe amplitude of each order of harmonics have no frame to frame variations, the fringe pattern in Eq. (12) can be represented as

$${I^t}({i;j} )= A(j )+ \sum\nolimits_{k = 1}^P {{B_k}(j )\cos [{k\varphi (j )+ k\delta (i )} ]} .$$

With the existing phase and updated phase-shifts, the background intensity and all the fringe amplitudes are estimated pixel-wisely through LSF_MI with a (P+1)×(P+1) least squares matrix whose entry at s-th row and t-th column as $\mathop \sum \nolimits_{i = 0}^{F - 1} \{{\cos [{s\varphi (j )+ s\delta (i )} ]\cos [{t\varphi (j )+ t\delta (i )} ]} \}$. Then, similar to the phase-shift optimization in the first step, the phase is updated by LSF_ES with exhaustive search within a range of $\pi$ rads with a step of 0.0001 rad.

Hoang’s algorithm starts with the phase and phase-shifts extracted by AIA and shares the same stopping criteria as AIA which are shown in Eqs. (10) and (11). We set P=3 for sufficient harmonics rejection.

The size of the least squares matrices in Hoang’s algorithm is (P+1)×(P+1) which is larger than 3×3 in AIA but smaller than (2P+1)×(2P+1) in Xu’s algorithm. Its random noise immunity also lies in between them [21,30]. A more significant problem of this algorithm is the slow computation speed due to the exhaustive search. In each iteration, π/0.0001≈31400 times of calculation are needed for the exhaustive search of phase of every pixel and phase-shift of every frame, which causes a huge computational cost. For example, to extract phase from 7 frames of fringe patterns with 256×256 pixels for 10 iterations, the algorithm requires 31400×(256×256 + 7)×10≈2×1010 times of calculations, which could cost hours or even longer to extract phase.

2.4 Algorithm for phase-shift tilts rejection

Based on AIA, multiple PSAs immune to phase-shift tilts have been proposed [2427]. In this sub-section, we introduce the phase-tilt iteration (PTI) algorithm [25] which has been demonstrated to be most effective for phase-shift tilts [24,25]. For PTI, the fringe pattern model used is a simplified case of Eq. (2) with P=1, Q=1, and an assumption that the background intensity and fringe amplitude have no frame to frame variation, i.e.,

$$\begin{aligned} {I^t}({i;j} )&= A(j )+ B(j )\cos [{\varphi (j )+ \delta ({i;j} )} ]\\ &= A(j )+ B(j )\cos [{\varphi (j )+ {\delta_0}(i )+ {\alpha_1}(i )x + {\beta_1}(i )y} ], \end{aligned}$$
where δ0, α1 and β1 are the parameters to represent the phase-shift tilts.

The first step of PTI is the same as that of AIA, except that the assumed known phase-shifts δ(i; j) is spatially tilted instead of a constant value. Through the same LSF_MI, the phase, background intensity and fringe amplitude can be estimated.

In the second step for phase-shift estimation, PTI first looks at a specific column (y = Y) and further re-writes the fringe signal as

$$\begin{aligned} {I^t}({i;{j_{x,Y}}} )&= A({{j_{x,Y}}} )+ B({{j_{x,Y}}} )\cos [{\varphi ({{j_{x,Y}}} )+ {\delta_{0Y}}(i )+ {\alpha_1}(i )x} ]\\ &= A({{j_{x,Y}}} )+ B({{j_{x,Y}}} )\cos [{{\Phi _r}({i;{j_{x,Y}}} )+ {p_r}(i )} ], \end{aligned}$$
where jx,Y is the index of a pixel (x,Y) in the Y-th column; Φr(i; jx,Y)=φ(jx,Y)+α1(i)(x-r) and pr(i)=δ0Y(i)+α1(i)r with r is an integer number. In this equation, Φr(i; jx,Y) is considered as known where the value of α1(i) uses the one from the last iteration, while the one in pr(i) is treated as unknown. Now, Eq. (18) can be linearized as follows,
$${I^t}({i;{j_{x,Y}}} )= A({{j_{x,Y}}} )+ {b_r}(i )B({{j_{x,Y}}} )\cos [{{\Phi _r}({i;{j_{x,Y}}} )} ]+ {c_r}(i )B({{j_{x,Y}}} )\sin [{{\Phi _r}({i;{j_{x,Y}}} )} ],$$
where Φr(i; jx,Y), A(jx,Y) and B(jx,Y) are known while br(i)=cos[pr(i)] and cr(i)=−sin[pr(i)] are the two unknowns that can be estimated through LSF_MI. Subsequently, pr(i) is calculated. It is noteworthy that the integer r can be changed so that a series of pr(i) is obtained. In PTI, r is scanned from 0 to Nx/4. Since pr(i)=δ0Y(i)+α1(i)r, both the δ0Y(i) and α1(i) can be calculated. The same process is repeated at a specific row (x = X) to obtain the δ0X(i) and β1(i). By “averaging” δ0X(i) and δ0Y(i), δ0(i) is determined [25].

The δ0(i), α1(i) and β1(i) are initialize by linear fitting of phase-shifts extracted by AIA from sub-segments of the fringe patterns and stops either when the maximum iteration number is reached as shown in Eq. (11) or when the change in δ0(i), α1(i) and β1(i) of all the frames are smaller than tolerance ε, ε1 and ε2, respectively. We set ε=1×10−3 rad and ε1= ε2=1×10−5 rad [25].

PTI can effectively solves the phase-shift tilts problem. However, when spatially non-linear phase-shift distribution or intensity harmonics are presented in the fringe patterns, the accuracy of the phase will be compromised.

2.5 Algorithm for both harmonic and phase-shift tilts rejection

Model based phase-shifting interferometry (MPSI) considers and incorporates all the mentioned error sources and has been implemented in the Zygo’s product, QPSI [28]. In general, the fringe pattern model used in MPSI can be considered as the one in Eq (2) with simplification of P>1 and Q=1, together with an assumption that the background intensity and fringe amplitude have no frame to frame variation. In addition, the intensity harmonics is modeled from multiple reflection of a Fizeau interferometer, resulting in their final model as follows,

$${I^t}({i;j} )= A(j )+ \sum\nolimits_{k = 1}^P {B(j ){{({ - g} )}^{k - 1}}\cos \{{k[{\varphi (j )+ {\delta_0}(i )+ {\alpha_1}(i )x + {\beta_1}(i )y} ]} \}} ,$$
where g is the parameter representing the intensity harmonics induced by multiple reflection.

MPSI classifies all the parameters to be estimated into two groups, those frame-independent, including A(j), B(j), g and φ(j), and those pixel-independent, including δ0(i), α1(i) and β1(i). In the first step to estimate the pixel-independent parameters, MPSI linearizes Eq. (20) by consider the small increments of these parameters and then estimate these increments by linear regression, with which, the pixel-independent parameters are updated.

The second step is in principle the same as the first step of Xu’s algorithm. With the calculated phase-shift δ(i; j) as the known input, all the frame-independent parameters including the phase are calculated.

We highlight that MPSI initializes φ(j) by measuring the phase using an existing PSA [35] and initializes A(j) and B(j) from the peaks and valleys of all the frames. Such initialization needs extra efforts but avoids a random guess. The above iteration stops when the change in normalized least squares error shown in Eq. (4) is smaller than the pre-defined tolerance ξ, or when maximum iteration number is reached as shown in Eq. (11). We set of ξ=0.01% and P=3 [28].

MPSI addresses both intensity harmonics and phase-shift tilts in the same time. As already mention in [28], for intensity harmonics, only multi-beam reflection is considered in MPSI, while other possible causes are ignored; for the non-uniform phase-shift distribution, only tilt is considered, while spatially nonlinear phase-shifts are not modeled. In addition, similar to Xu’s algorithm, the immunity of MPSI to random noise is poor due to the inverse of (2P+1)×(2P+1) least squares matrix in the second step [21,30]. Thus, in MPSI, phase-shifts evenly distributed among 2π are suggested. This treatment is reasonable in practical use, but it loses the capability of dealing with unknown and arbitrary phase-shifts, which is the context of this paper.

3. Proposed GIA

In this section, the principle and details of GIA which is intended for phase extraction from fringe pattern with unknown phase-shifts and all the mentioned possible error sources are described.

3.1 GIA strategy

To extract phase accurately with the presence of various error sources, generality and simplicity are pursued, where generality is meant for a capability to process fringe patterns with different error sources, while simplicity is meant for a solution without many hand-crafted manipulations.

3.1.1 Adopting a general fringe pattern model

As reviewed and discussed in Sec. 2, AIA, Xu’s algorithm, Hoang’s algorithm, PTI and MPSI restrict their fringe pattern model in one way or another, and thus sacrifice the immunity to certain error source(s). GIA aims at phase extraction from the general fringe model described by Eqs. (1)-(3) for good immunity to all the errors mentioned in this paper. Now, given that we have F frames and each frame has N pixels, adding up to F×N measured intensity values. However, there are F×N×(P+1) parameters related to background intensity and fringe amplitudes, N parameters related to phase distribution, (F-1)×(Q+1)×(Q+2)/2 parameters related to phase-shifts (under assumption that phase-shifts of first frame δ(0; j) = 0), and in total F×N×(P+1)+N+(F-1)×(Q+1)×(Q+2)/2 unknowns, making the problem theoretically underdetermined.

As did in all the algorithms reviewed in Sec. 2, we continue to assume that there is no frame to frame variations in the background intensity and fringe amplitudes. This assumption is reasonable in practice. Thus, the fringe pattern in Eq. (2) can be rewritten as

$${I^t}({i;j} )= A(j )+ \sum\nolimits_{k = 1}^P {{B_k}(j )\cos [{k\varphi (j )+ k\delta ({i;j} )} ]} .$$

As shown, the number of unknowns is reduced to N×(P+2)+(F-1)×(Q+1)×(Q+2)/2. These unknowns can be optimized if F>(P+2) since (F-1)×(Q+1)×(Q+2)/2 is typically much smaller than N. Conservatively, we can set F≥(2P+1) similar to other harmonic immune algorithms [22,23]. This fringe pattern model used in GIA is the most general one comparing with the other reviewed algorithms and needs no pre-measurement of the phase-shifts [12,2228].

3.1.2 Grouping unknowns for effective optimization

The adoption of the general fringe model should not sacrifice the accuracy and robustness of phase extraction; thus, GIA should be equipped with a powerful optimization solution. A natural choice is to simultaneously optimize all the unknowns by minimizing the objective function described in Eq. (4). However, it is difficult to converge and sensitive to noise due to the long unknown vector for optimization [30], which is consistent to our attempts in this phase extraction problem.

Another option, which has been used in all the algorithms mentioned in Sec. 2, is to divide the unknowns into two groups and use the univariate search method to iteratively optimize them. These two groups are the frame-independent parameters (background intensity, fringe amplitude and phase) and the pixel-independent parameters (background intensity, fringe amplitude and phase-shift related parameters). However, we find that (i) in AIA, Xu’s algorithm and Hoang’s algorithm, the background intensity and fringe amplitudes appear in both groups, resulting in a strong assumption that the background and fringe amplitudes are both frame- and pixel- independent; (ii) in Xu’s algorithm and MPSI where intensity harmonics are considered, their pixel-wise estimation step estimates the fringe amplitude and phase as combined variables through LSF_MI, where the inverse of a large matrix leads to poor immunity to random noise [21,30].

In GIA, we classify the unknowns into three groups based on the different natures of the parameters and their roles in Eq. (21): (i) Group B including the background intensity A(j) and fringe amplitude Bk(j) which are linear in Eq. (21) and have small spatial variations; (ii) Group P including the phase φ(j) which is non-linear in the fringe model and pixel-wisely independent; and (iii) Group S including phase-shift related parameters αuv(i) which is non-linear in the fringe model and frame-wisely independent. In other words, we keep the group of pixel-independent parameters as our Group S; meanwhile, we further divide the frame-independent parameters into two groups, group B and group P. Such division gives us two advantages: (i) as shown in Sec. 2.3 and 2.5, when group B and group P are considered together, the unknown vector length is 2P+1; when they are treated separately, the unknown vector lengths of groups B and P are P+1 and 1, respectively. The shorter lengths are more resilient to noise; (ii) since the parameters in group B have slow spatial variation, we can include neighboring pixels in the fringe amplitude optimization to avoid overfitting when using large value of P and to improve signal to noise ratio and hence the immunity to noise. The details are given in Sec. 3.2.2. Although the parameters could be further divided into more groups, such practice will not further improve the accuracy of GIA since the parameters in each group are independent to others in the same group.

Once all the parameters have been grouped, each group is estimated assuming that the other two groups of unknowns are fixed. Due to its known effectiveness and robustness, the Levenberg-Marquart method (briefly introduced in Appendix A) is selected for the estimation of all the three groups by optimizing their respective least squares errors, which is then called LSF_LM for convenience. In other words, LSF_LM is used for the optimization of all the three groups. We only mention that group B can also be optimized by LSF_MI, but LSF_LM performs better because the damper factor in the Levenberg-Marquart method can effectively avoid the matrix singularity problem.

3.2 GIA algorithm

This section provides the details of the algorithm for easy implementation.

3.2.1 Flowchart

Due to the simplicity of the GIA’s strategy, the flowchart of GIA is also simple, as shown in Fig. 1. With the initialization at the beginning, GIA sequentially optimizes groups B, P and S, and iterates the process until converged.

 figure: Fig. 1.

Fig. 1. Flowchart of GIA.

Download Full Size | PDF

3.2.2 Optimizing group B

The background intensity and fringe amplitudes are estimated by minimizing the objective function E shown in Eq. (4) involving the target pixel jx,y and its neighboring pixels within a (2w+1)×(2w+1) window through LSF_LM. The background intensity and fringe amplitudes are assumed to be constant within the window,

$$A({{j_{x + s,y + t}}} )= A({{j_{x,y}}} ),$$
$${B_k}({{j_{x + s,y + t}}} )= {B_k}({{j_{x,y}}} ),$$
where s = −w, …, 0, …, w and t = −w, …, 0, …, w. These two assumptions have been enforced to all the pixels in the entire frame in AIA, Xu’s and Hoang’s algorithms [12,22,23], which is relaxed to w=5∼10 in GIA. When w=0, the pixel-dependence is removed.

According to Appendix A, to use the Levenberg-Marquart method, we need to set up the vector of unknowns x with a size of 1×(P+1), the residual vector r(x) with a size of [(2w+1)×(2w+1)×F]×1 and the Jacobian matrix J(x) with a size of [(2w+1)×(2w+1)×F]×(P+1) for each target pixel jx,y:

$${x_k} = \left\{ \begin{array}{l} A({{j_{x,y}}} )\textrm{ if }k = 0,\\ {B_k}({{j_{x,y}}} )\textrm{ if }k = 1,\textrm{ }\ldots ,P, \end{array} \right.\textrm{ }$$
$${r_{({2w + 1} )\times ({2w + 1} )\times i + l}}({\mathbf x} )= {I^t}({i;{j_{x + s,y + t}}} )- I({i;{j_{x + s,y + t}}} ),$$
$${J_{({2w + 1} )\times ({2w + 1} )\times i + l,k}}({\mathbf x} )= \frac{{\partial {r_{({2w + 1} )\times ({2w + 1} )\times i + l}}({\mathbf x} )}}{{\partial {x_k}}} = \cos [{k\varphi ({{j_{x + s,y + t}}} )+ k\delta ({i;{j_{x + s,y + t}}} )} ],$$
where the subscript l=0, 1, …, (2w+1)×(2w+1)-1 is the index of pixels in the window.

3.2.3 Optimizing group P

This step optimizes the phase pixel-wisely by minimizing EP(j) using the LSF_LM. For a target pixel j, we have the following specifications for x (1×1), r(x) (F×1) and J(x) (F×1) of the Levenberg-Marquart method:

$${\mathbf x} = \varphi (j ),$$
$${r_i}({\mathbf x} )= {I^t}({i;j} )- I({i;j} ),$$
$${J_{i,1}}({\mathbf x} )= \frac{{\partial {r_i}({\mathbf x} )}}{{\partial {\mathbf x}}} = \sum\nolimits_{k = 1}^P { - k{B_k}(j )\sin [{k\varphi (j )+ k\delta ({i;j} )} ]} ,$$
with i=0, 1, …, F-1.

3.2.4 Optimizing group S

This step optimizes the phase-shift related parameters by minimizing EF(i) using the LSF_LM. The optimization is frame-wise and is executed for the second onward frames, as the phase-shifts of first frame is assumed to be 0. For a target frame i, we have the following specifications for x (1×[(Q+1)×(Q+2)/2]), r(x) (N×1) and J(x) (N×[(Q+1)×(Q+2)/2]) of the Levenberg-Marquart method:

$${x_{u({u + 1} )/2 + v}} = {\alpha _{uv}}(i ),$$
$${r_j}({\mathbf x} )= {I^t}({i;j} )- I({i;j} ),$$
$${J_{j,}}_{u({u + 1} )/2 + v}({\mathbf x} )= \frac{{\partial {r_j}({\mathbf x} )}}{{\partial {x_{u({u + 1} )/2 + v}}}} = \sum\nolimits_{k = 1}^P {\frac{{ - k{x^u}{y^{u - v}}{B_k}(j )\sin [{k\varphi (j )+ k\delta ({i;j} )} ]}}{{{{({{N_x} - 1} )}^u}{{({{N_y} - 1} )}^{u - v}}}}} ,$$
with j=0, 1, …, N-1, u=0, 1, …, Q and v=0, 1, …, u.

3.2.5 Convergence check

The three groups of parameters are iteratively optimized. GIA initializes phase, phase-shifts, background intensity and fringe amplitude by AIA [12] and set all higher order harmonics fringe amplitude Bk (k ≥ 2) and phase-shifts non-uniformity αuv (u≥1, v≥1) to zero. GIA stops until the maximum iteration number, M, is reached as shown in Eq. (11), or the absolute change of the least square error shown in Eq. (4) is sufficiently small, i.e.,

$${{|{{E^{m + 1}} - {E^m}} |} / {{E^m}}} < \zeta .$$

In later simulation and experiments, we set: (i) ξ=0.01%, (ii) η=1×10−4, (iii) M = MLM=200 (η and MLM are parameter of convergence check of the Levenberg-Marquart method as shown in Appendix A), (iv) w=8, (v) P=3 and (vi) Q=4 which will result in 15 coefficients of the Maclaurin polynomial αuv for each frame; such practice will be sufficient for phase-shifts non-uniformity removal and will be good enough for convergence and accuracy evaluation.

4. Simulation evaluation

In this section, the performance of GIA is examined by comparing with other algorithms through simulation. The simulation platform and algorithm specifications are introduced first, followed by fringe pattern generation. Then, two specific examples with random phase-shifts and regular phase-shifts are processed to exemplify the performances. Finally, more examples are processed to examine GIA’s overall convergence and accuracy. Although GIA is characterized by its generality and simplicity, its performance regarding convergence, accuracy and robustness is more crucial and carefully examined.

4.1 Simulation platform and algorithm specifications

The simulation is carried out with MATLAB 2020a on a workstation equipped with an Intel Xeon E5-2699 v3 processors (18 cores, 2.3Ghz main frequency) and 256 GB RAM. For the implementation of the existing algorithms used for comparison, we follow their descripts in the respective papers as strict as possible. Various parameters used in all the algorithms are summarized in Table 1, with the following additional parameters for GIA: (i) η=1×10−4 (ii) MLM=200 (iii) w=8 and (iv) Q=4.

Tables Icon

Table 1. Parameter settings for different algorithms.

4.2 Fringe pattern generation

We call it a complete example where the image size, frame number, phase are chosen and then fixed, but all other parameters vary in order to simulation different situations. The image size is Nx=Ny=128∼1024, i.e., the image size in one direction varies between 128 and 1024; the frame number is F=3∼19, i.e., the frame number is randomly generated between 3 and 19; and the phase distribution is created from the following three types,

$$\varphi (j )= {\omega _x}x + {\omega _y}y,$$
$$\varphi (j )= 2{\omega _p}\pi \times {{[{{{({x - {N_x}/2} )}^2} + {{({y - {N_y}/2} )}^2}} ]} / {[{{{({{N_x}/2} )}^2} + {{({{N_y}/2} )}^2}} ]}} + {\omega _x}x + {\omega _y}y,$$
$$\varphi (j )= {\omega _p} \times {{{\textrm peaks}\_s}}({{N_x},{N_y}} )+ {\omega _x}x + {\omega _y}y,$$
where peaks_s uses the MATLAB function “peaks”, linearly scaled it into the range of [0, 2π], with ωp=1∼10, ωx=0∼2 and ωy=0∼2. These simulated known phase distributions are called ground-truth phases, and will be used to evaluate the quality of the extracted phases by different algorithms.

In our simulation, one complete example includes ten sets of phase-shifted fringe patterns. Set I includes ideal phase-shifted fringe patterns as follows,

$${I_o}({i;j} )= A(j )+ B(j )\cos [{\varphi (j )+ {\delta_0}(i )} ],$$
with A(j) = 0.5 and B(j) = 0.4. Set I can be derived from Eq. (2) without intensity harmonics problem (i.e., P=1) or phase-shifts non-uniformity problem (i.e., Q=0). To increase the variety of the fringe patterns, multiple choices are given for the phase-shifts without tilt or non-uniformity, i.e., δ0, including randomly generated from [0, 2π], or specially set as δ0(i)=ilπ/k0 with k0$ {\mathbb{Z}}$+ and l$ {\mathbb{Z}}$+.

Based on Set I, we then generate eight more sets with different errors added as shown in Table 2, where a tick means this type of error is presented in this set. The added AWGN has a mean of zero and a standard deviation of σI ranging from 0.02∼0.4; the intensity harmonics are simulated using a gamma model, i.e., Ioγ is used to represent the non-linear fringe intensity [36] where γ=0.3∼3; the parameters of phase-shifts are set as random values with a mean of zero and a standard deviation of σα with α00=δ0, where σα=0.05∼0.5. In addition, Set X modifies Set IX by adding 50% spatial variations to the background intensity and fringe amplitude as

$$A(j )= 0.5 - 0.25\left[ {\frac{{{{({x - {N_x}/2} )}^2} + {{({y - {N_y}/2} )}^2}}}{{{{({{N_x}/2} )}^2} + {{({{N_y}/2} )}^2}}}} \right],$$
and B(j) = 0.8×A(j).

Tables Icon

Table 2. Errors simulated in fringe patterns in Group (i)-(x).

For the sake of easier and clearer discussions, these 10 sets of simulated fringe patterns are further classified into 4 categories as shown in first column of Table 2. Categories A-D highlight the presence of neither intensity harmonic nor phase-shift non-uniformity, only intensity harmonic, only phase-shift non-uniformity, and both of them, respectively.

Among 500 complete examples we tested, we specify two of them to show the detailed comparison results: (i) Nx=Ny=1024 (ii) F=8; (iii) phase distribution is based on Eq. (36) with ωp=10, ωx=0 and ωy=0.25; (iv) for the first example, the phase-shifts without tilt or non-uniformity are set as,

$${\delta _0} = \left[ {\begin{array}{llllllll} 0&{0.440} &{0.795}&{1.246}&{1.601}&{2.665}&{2.695}&{3.449} \end{array}} \right],$$
while for the second example, δ0(i)=iπ/4; (v) σI=0.04; (vi) γ=1.5; and (vii) σα=0.5. The first example is selected because the phase-shifts are “adversary”: they gather within half period (from 0 to π) and the 6-th and 7-th values are similar. On the contrary, the second example highlights regular phase-shifts evenly distributed among [0, 2π] as suggested by MPSI [28], which is most “friendly”. Even though the phase-shifts are simulated, they are assumed to be unknown for each algorithm, except for MPSI. In MPSI, the ground-truth phase, and phase-shifts without tilt and non-linearity, i.e. δ0, are used for initialization; if we use the AIA results to initialize MPSI, sometimes bad performances will be obtained.

Figure 2 shows some parameters of the first example: the background intensity expressed by Eq. (35) in (a), phase-shifts distribution of second frames in (b), the ground-truth phase expressed by Eq. (36) in (c), and first pattern from Set X in (d).

 figure: Fig. 2.

Fig. 2. Simulation Parameters (a) Background intensity; (b) Phase-shift of second frame; (c) Ground-truth phase; (d) First fringe pattern from Set X.

Download Full Size | PDF

4.3 Assessment through an example with random δ0

We use all the mentioned algorithms, AIA, Xu’s algorithm, Hoang’s algorithm, PTI, MPSI and the proposed GIA, to process the first example. The extracted phase by each algorithm is compared with the ground-truth phase, and their difference is recognized as the phase error. The root mean squares value of the obtained phase error (RMSE) is then computed. The RMSEs of all the algorithms are shown in Fig. 3 and Table 3. For clearer display, the 4 categories are set with different colors, i.e., pink for Category A, green for B, blue for C, and yellow for D. The observations from the simulation results are discussed below. In Category A, since Set I is ideal, all the algorithms achieve very low RMSEs (less than 1×10−4 rad), indicating their proper implementation and excellent performance. However, with the noise presenting in Set II, the RMSEs of Xu’s algorithm and MPSI become unacceptably high due to the large condition number (more than 20,000) of the matrix used in LSF_MI, which agrees with the analysis in Secs. 2.3 and 2.5. The RMSEs of the rest of the algorithms also increase but are less than 0.1 rad. GIA shows the least RMSE among the others.

 figure: Fig. 3.

Fig. 3. RMSEs of extracted phase from example with random δ0.

Download Full Size | PDF

Tables Icon

Table 3. RMSEs of extracted phase from example with random δ0 (Unit: rad).

In Category B, intensity harmonics are highlighted. For Set III, the harmonics-aware algorithms (Xu’s algorithm, Hoang’s algorithm, MPSI and GIA) perform better than the others (AIA and PTI). Set IV shows the noise influence again. GIA produces the least RMSE.

In Category C, phase-shift non-uniformity is highlighted. Set V shows that tilt-aware algorithms (PTI, MPSI and GIA) perform better than the others when phase-shift tilts are introduced, as expected. In Set VI where phase-shift nonlinearity presents, PTI and MPSI do not consider the problem and naturally perform worse, while GIA considers the problem and obtains excellent immunity. Besides, MPSI is worse than PTI due to the arbitrary phase-shifts. Set VII shows the noise influence.

In Category D, both intensity harmonics and phase-shift non-uniformity are highlighted. GIA performs better than all the other algorithms, which demonstrates the good immunity of GIA to all the mentioned error sources, even with the non-uniform background intensity and fringe amplitudes in Set X.

4.4 Assessment through an example with regular δ0

As already highlighted, the phase-shifts in the first example are “adversary”, although it could happen. We notice that MPSI does not work well for this example. Indeed, in the MPSI work [28], regular phase-shifts are suggested. We thus run the second complete example with regular δ0=iπ/4. The RMSEs are shown in Fig. 4 and Table 4.

 figure: Fig. 4.

Fig. 4. RMSEs of extracted phase from example with regular δ0.

Download Full Size | PDF

Tables Icon

Table 4. RMSEs of extracted phase from example with regular δ0 (Unit: rad).

The main observations and remarks are as follows: (i) all the algorithms perform better than the previous example, showing that regular phase-shifts are preferred if possible; (ii) GIA is able to outperform the others for both random and regular phase-shifts, indicating that pursuing generality does not affects its performance; (iii) comparing with the first example, MPSI benefits most from the selecting of regular phase-shifts, making it the second best algorithm. Although PTI does not consider the intensity harmonics, the problem is seen less significant than other error sources in this example, making its performance similar to that of MPSI. It is worth mentioning that in Sec. 4.3 and 4.4 we focus on the accuracy of the extracted phase. In the meantime, background intensity and fringe amplitude extracted by GIA are both better than those by the other algorithms.

4.5 Convergence and accuracy of GIA

In the above two examples, one with “adversary” random phase-shifts and the other with “friendly” uniformly distributed phase-shifts, GIA outperforms the other algorithms. Nevertheless, since GIA is an iterative algorithm with three Levenberg–Marquardt optimizations in each iteration, it is necessary to test its convergence on more examples. Thus, we carried out 500 complete examples (i.e., 5000 fringe sets) with different settings as described in Sec. 4.2. In all our simulations, GIA consistently converges with better accuracy than other algorithms (Hoang’s algorithm is excluded in these simulations due to its slow computation) if the density of the fringe patterns is more than one fringe over the entire frame. This requirement is actually the prerequisite of AIA’s convergence [21]. If good phase and phase-shift estimation can be obtained by other means to initialize GIA, then this requirement can be avoided. In Sec. 2, we mentioned that some special phase-shifts will cause singular matrices in Xu’s algorithm and MPSI, which has been tested not a restriction for GIA.

Although GIA has been tested to have the best accuracy among different algorithms, there are still a few factors affecting its accuracy:

  • (i) GIA is found to converge robustly against noise, but it does not have noise suppression capability. Figure 5 (a) shows RMSEs of the phase extracted by GIA from Set II with 8 frames of fringe pattern, regular phase-shifts δ0(i)=iπ/4 and different noise levels σI. The relationship between RMSE and the noise level is approximately linear.
  • (ii) The frame number will affect the accuracy of GIA. Set II is simulated with regular phase-shifts δ0(i) = 2iπ/F, fixed noise level σI=0.04 and different frame number F. The RMSEs of the phase extracted by GIA against different frame numbers is shown in Fig. 5 (b). Clearly, more fringe patterns should be acquired if lower RMSE is desired.
  • (iii) The phase-shifts may affect the accuracy of GIA. Set II is simulated with 8 frames of fringe patterns, fixed noise level σI=0.04, and different regular phase-shifts δ0(i)=ilπ/4. The RMSEs are shown in Fig. 5 (c). Although, GIA could converge with small phase-shift value between each two consecutive frames, the accuracy of phase will degrade. Hence, phase-shift control before an experiment starts is preferred and suggested, if possible, as indicated in Sec. 4.4.
  • (iv) The initialization may affect the accuracy of GIA. Set V is simulated with 8 frames of fringe pattern, regular phase-shifts δ0(i)=iπ/4, fixed noise level σI=0.04, and different values of σα. MPSI and GIA are used to extract the phase with initialization of AIA extracted and ground-truth phase and phase-shifts. The RMSEs are plotted in Fig. 5(d), which shows that, on one hand, initialization is essential for successful phase-extraction in both MPSI and GIA; on the other hand, both algorithms work well for σα≤0.6, which corresponds to a phase-tilt angle of more than 2π/3 between two consecutive frames. Comparing with MPSI, GIA is less sensitive to the initialization. In the meantime, we also noticed that a more accurate initialization will also lead to better performance for each step of GIA.
  • (v) Similar to all iterative algorithms, the accuracy of GIA is related to the convergence setting, ξ. A simulation is carried out with fringe pattern in Set X of the second example in Sec. 4.4 with a varying ξ, as shown in Fig. 5(e). As seen from the curve, ξ=0.01% is a good choice for MPSI and GIA, as we did in our simulations and later experiments.
  • (vi) We now discuss the setting of P and Q in GIA. Consider Categories A first, where there is nether intensity harmonics nor phase-shift nonlinearity. If we have such pre-knowledge, it is most reasonable to set P=1 and Q=0. However, if we do not have such knowledge, it is safer to set higher values for P and Q. The setting with pre-knowledge is expected to work better, but our simulation results show that two settings achieve the same accuracy. This is because GIA is able to converge well to the same pre-set convergence requirement from two different settings. Meanwhile, the convergence of the lower P and Q values is about 2.87 times faster. The same phenomenon happens to Categories B and C as well.

 figure: Fig. 5.

Fig. 5. Accuracy of GIA in different conditions (a) RMSEs with respect to different σI; (b) the RMSEs with respect to different F; (c) the RMSEs with respect to different l; (d) RMSEs with respect to different σα; (e) RMSEs with respect to ξ, the change in E.

Download Full Size | PDF

4.6 Computation time of GIA

The computation time of GIA is compared with other algorithms through simulation with Set X in Sec. 4.3. For fair comparison, we force the iteration number of every algorithm to be 10. The computation time of each algorithm, as an average of 10 executions, is given in Table 5, where Hoang’s algorithm needs much longer computation time while the others can extract the phase within 15 minutes with the MATLAB implementation. GIA requires about 10 minutes, and its acceleration by C++ and GPU is considered as our future work.

Tables Icon

Table 5. Computation time of each algorithm (Unit: second).

5. Experimental verification of GIA

To furtherly verify the effectiveness of GIA, we conducted serval optical testing experiments in which intensity harmonics and phase-shift non-uniformity are purposely introduced. The performance of GIA is evaluated and compared with all the other algorithms.

5.1 Experimental set up and ground-truth phase

The test is carried out with a Fizeau interferometer as sketched in Fig. 6. We used a wedged optical window as the reference and a piece of mirror as the test object. The phase-shifts were introduced by tuning the wavelength of the interferometer’s laser [37]. In the experiments, the frame size is Nx=Ny=1024.

 figure: Fig. 6.

Fig. 6. Sketch of the Fizeau interferometer.

Download Full Size | PDF

From this system, we first collected 17 frames of fringe patterns with phase-shifts of iπ/4 (i=0,1, …, F-1) and without experimentally introduced intensity harmonics and non-uniform phase-shift distribution. The first frame is shown in Fig. 7 (a), the phase extracted by AIA from 17 frames of fringe patterns is shown in Fig. 7 (b) and is used as the ground-truth phase of the test object.

 figure: Fig. 7.

Fig. 7. Fringe patterns and ground-truth phase. (a) One of acquired fringe patterns. (b) The ground-truth phase.

Download Full Size | PDF

5.2 Experiment results

Next, we capture two complete testing examples corresponding to the two simulation examples examined in Secs. 4.3 and 4.4, i.e., the phase-shifts are introduced randomly and controlled to be regular, respectively. However, in real experiments, random noise is unavoidable, making Sets I, III, V, VI and VIII impossible. Furthermore, we do not alter the distribution of background intensity and fringe amplitude, and thus Set X is excluded. Only four sets of experimental data are captured in each example: Sets II, IV, VII and IX, highlighting one error source (random noise), two error sources (intensity harmonic and random noise), two error sources (phase-shifts non-uniformity and random noise), and three error sources (intensity harmonic, phase-shifts non-uniformity and random noise), respectively.

In the experiments, the intensity harmonics are introduced by controlling the nonlinear response of the camera with γ=1.5 [38], and the phase-shift non-uniformity is introduced by a mechanical actuator to randomly excite the mounting of the interferometer’s Ref. [19]. The frame number is set as F=8. The phases extracted from the four sets of fringe patterns are compared with the ground-truth phase. The initialization of MPSI uses the ground-truth phase and the phase-shifts extracted by AIA. The phase, phase-shifts of second frame, and tilt removed phase-shifts of the same frame extracted by the GIA from Set IX with random phase-shifts are shown in Fig. 8. From Fig. 8(c), a high order phase-shifts distribution due to the non-rigid movement between the wedged window and the mirror are observed.

 figure: Fig. 8.

Fig. 8. Extracted phase and phase-shifts by GIA (unit: rads) (a) Phase; (b) Phase-shift; (c) Tilt removed phase-shift.

Download Full Size | PDF

Since all the parameters have been extracted, we can use them to reconstruct the fringe patterns and then compare them with the captured fringe pattern to check the royalty of the reconstruction. A smaller difference indicates better reconstruction and hence more accurate extraction of the parameters. The differences using all the algorithms are shown in Fig. 9, where the least difference is generated by GIA.

 figure: Fig. 9.

Fig. 9. Difference between the measured intensity and fringe pattern generated by (a) GIA; (b) AIA; (c) Xu’s algorithm; (d) Hoang’s algorithm; (e) PTI; (f) MPSI.

Download Full Size | PDF

The RMSEs of the extracted phases of the two complete examples are shown in Fig. 10 and Table 6. The following are highlighted: (i) for the four sets in both examples, the observations from the results are consistent with the examples in simulations as shown in Sec. 4.3 and 4.4; (ii) comparing with random phase-shifts, regular phase-shifts produce lower RMSEs; (iii) Comparing with other algorithms, GIA have better accuracy in both examples.

 figure: Fig. 10.

Fig. 10. RMSEs of extracted phase in experiment. (a) with random phase-shifts (b) with regular phase-shifts

Download Full Size | PDF

Tables Icon

Table 6. RMSEs of extracted phase in experiment (Unit: rad).

6. Conclusion

A general, simple and effective GIA is proposed for phase extraction from fringe patterns with unknown phase-shifts, random noise, intensity harmonics and phase-shifts non-uniformity. In GIA, a general fringe pattern model with many unknowns is adopted. For better accuracy and convergence, these unknowns are divided into three groups, i.e., the background intensity and fringe amplitudes, the phase, and the phase-shift related parameters, and are optimized through univariate search technique group by group through minimizing the least squares error between the noiseless fringe intensity and measured fringe intensity. The Levenberg-Marquart method is selected to optimize all the groups due to its excellent accuracy and robustness. Even with such a simple structure, GIA has been demonstrated to outperform five most relevant algorithms in both simulations and experiments.

Appendix A: The Levenberg-Marquart method for non-linear optimization

We briefly introduce the Levenberg-Marquart optimization method which is used in the proposed GIA. The Levenberg-Marquart method is very effective, widely used, and can be considered as a standard optimization module nowadays [31,32]. We thus only outline the important principle, leaving those details which do not affect the understanding of this paper to Ref. [31].

Given the following least squares objective function to be optimized regarding an unknown vector x as

$$f({\mathbf x} )= \frac{1}{2}{||{{\mathbf r}({\mathbf x} )} ||^2} = \frac{1}{2}\sum\nolimits_{i = 0}^{m - 1} {{{[{{r_i}({\mathbf x} )} ]}^2}} ,$$
where r(x) is a vector whose entries are the residuals, ri(x). The following approximation is made for r(x):
$${\mathbf r}({{\mathbf x + h}} )\approx {\mathbf r}({\mathbf x} )+ {\mathbf J}({\mathbf x} ){\mathbf h},$$
where J(x) is the Jacobian matrix, with each entry as
$${J_{i,j}}({\mathbf x} )= \frac{{\partial {r_i}}}{{\partial {x_j}}}({\mathbf x} ).$$

Substituting Eq. (41) to Eq. (40) yields

$$f({{\mathbf x + \textbf h}} )= f({\mathbf x} )+ {{\mathbf h}^T}{{\mathbf J}^T}({\mathbf x} ){\mathbf r}({\mathbf x} )+ \frac{1}{2}{{\mathbf h}^T}{{\mathbf J}^T}({\mathbf x} ){\mathbf J}({\mathbf x} ){\mathbf h}.$$

By letting the gradient of f with respect to h be zero, the Gauss–Newton method finds h by solving

$${{\mathbf J}^T}({\mathbf x} ){\mathbf J}({\mathbf x} ){\mathbf h} ={-} {{\mathbf J}^T}({\mathbf x} ){\mathbf r}({\mathbf x} ).$$

The Levenberg-Marquart method add a damper μ to Eq. (44) as

$$[{{{\mathbf J}^T}({\mathbf x} ){\mathbf J}({\mathbf x} )+ \mu {\mathbf I}} ]{{\mathbf h}_{\textrm{LM}}} ={-} {{\mathbf J}^T}({\mathbf x} ){\mathbf r}({\mathbf x} ).$$

A large value of μ gives a short step in the steepest descent direction which is similar to a method called gradient descent, while a small value of μ gives a result similar to the Gauss–Newton method. Thus, the choice of μ is essential, which should be updated with the progress of optimization [31]. The Levenberg-Marquart method stops when

$${||{{{\mathbf J}^T}({\mathbf x} ){\mathbf r}({\mathbf x} )} ||_\infty } < \eta, $$
or
$${m_{LM}} > {M_{LM}}, $$
where η are a pre-defined tolerance, mLM is the index of iterations and MLM is the maximum number of iterations.

Funding

Economic Development Board - Singapore (S17-1579-IPP-II).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. K. Creath, “Phase-measurement interferometry techniques,” in Prog. Opt., vol. 26E. Wolf, ed. (Elsevier, 1988), pp. 349–393. [CrossRef]  

2. Y. Kim, K. Hibino, R. Hanayama, N. Sugita, and M. Mitsuishi, “Multiple-surface interferometry of highly reflective wafer by wavelength tuning,” Opt. Express 22(18), 21145–21156 (2014). [CrossRef]  

3. D. Malacara, Optical shop testing (John Wiley & Sons, 2007), pp. 547–655, Vol. 59.

4. J. E. Greivenkamp, “Generalized data reduction for heterodyne interferometry,” Opt. Eng. 23(4), 234350 (1984). [CrossRef]  

5. J. Degrieck, W. Van Paepegem, and P. Boone, “Application of digital phase-shift shadow Moiré to micro deformation measurements of curved surfaces,” Optics and Lasers in Engineering 36(1), 29–40 (2001). [CrossRef]  

6. Y.-C. Chu, W.-Y. Chang, K.-H. Chen, J.-H. Chen, B.-C. Tsai, and K. Y. Hsu, “Full-field refractive index measurement with simultaneous phase-shift interferometry,” Optik 125(13), 3307–3310 (2014). [CrossRef]  

7. J. H. Bruning, D. R. Herriott, J. Gallagher, D. Rosenfeld, A. White, and D. Brangaccio, “Digital wavefront measuring interferometer for testing optical surfaces and lenses,” Appl. Opt. 13(11), 2693–2703 (1974). [CrossRef]  

8. K. Hibino, B. F. Oreb, D. I. Farrant, and K. G. Larkin, “Phase-shifting algorithms for nonlinear and spatially nonuniform phase shifts,” JOSA A 14(4), 918–930 (1997). [CrossRef]  

9. Y. Surrel, “Design of algorithms for phase measurements by the use of phase stepping,” Appl. Opt. 35(1), 51–60 (1996). [CrossRef]  

10. K. Freischlad and C. L. Koliopoulos, “Fourier description of digital phase-measuring interferometry,” JOSA A 7(4), 542–551 (1990). [CrossRef]  

11. M. Servin, J. C. Estrada, and J. A. Quiroga, “The general theory of phase shifting algorithms,” Opt. Express 17(24), 21867–21881 (2009). [CrossRef]  

12. Z. Wang and B. Han, “Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms,” Opt. Lett. 29(14), 1671–1673 (2004). [CrossRef]  

13. J. Vargas, J. A. Quiroga, and T. Belenguer, “Phase-shifting interferometry based on principal component analysis,” Opt. Lett. 36(8), 1326–1328 (2011). [CrossRef]  

14. J. Deng, H. Wang, D. Zhang, L. Zhong, J. Fan, and X. Lu, “Phase shift extraction algorithm based on Euclidean matrix norm,” Opt. Lett. 38(9), 1506–1508 (2013). [CrossRef]  

15. J. Vargas, J. A. Quiroga, C. Sorzano, J. Estrada, and J. Carazo, “Two-step demodulation based on the Gram–Schmidt orthonormalization method,” Opt. Lett. 37(3), 443–445 (2012). [CrossRef]  

16. Y. Surrel, “Additive noise effect in digital phase detection,” Appl. Opt. 36(1), 271–276 (1997). [CrossRef]  

17. R. Schödel, A. Nicolaus, and G. Bönsch, “Phase-stepping interferometry: methods for reducing errors caused by camera nonlinearities,” Appl. Opt. 41(1), 55–63 (2002). [CrossRef]  

18. K. Larkin and B. Oreb, “Design and assessment of symmetrical phase-shifting algorithms,” J. Opt. Soc. Am. A 9(10), 1740–1748 (1992). [CrossRef]  

19. P. J. De Groot, “Vibration in phase-shifting interferometry,” J. Opt. Soc. Am. A 12(2), 354–365 (1995). [CrossRef]  

20. N. Bobroff, “Residual errors in laser interferometry from air turbulence and nonlinearity,” Appl. Opt. 26(13), 2676–2682 (1987). [CrossRef]  

21. Y. Chen and Q. Kemao, “Advanced iterative algorithm for phase extraction: performance evaluation and enhancement,” Opt. Express 27(26), 37634–37651 (2019). [CrossRef]  

22. J. Xu, Q. Xu, and L. Chai, “An iterative algorithm for interferograms with random phase shifts and high-order harmonics,” J. Opt. A: Pure Appl. Opt. 10(9), 095004 (2008). [CrossRef]  

23. T. Hoang, Z. Wang, M. Vo, J. Ma, L. Luu, and B. Pan, “Phase extraction from optical interferograms in presence of intensity nonlinearity and arbitrary phase shifts,” Appl. Phys. Lett. 99(3), 031104 (2011). [CrossRef]  

24. Y.-C. Chen, P.-C. Lin, C.-M. Lee, and C.-W. Liang, “Iterative phase-shifting algorithm immune to random phase shifts and tilts,” Appl. Opt. 52(14), 3381–3386 (2013). [CrossRef]  

25. M. Duan, Y. Zong, R. Zhu, and J. Li, “Phase-tilt iteration: Accurate and robust phase extraction from random tilt-shift interferograms,” Optics and Lasers in Engineering 142, 106595 (2021). [CrossRef]  

26. J. Xu, Q. Xu, and L. Chai, “Iterative algorithm for phase extraction from interferograms with random and spatially nonuniform phase shifts,” Appl. Opt. 47(3), 480–485 (2008). [CrossRef]  

27. Q. Liu, Y. Wang, F. Ji, and J. He, “A three-step least-squares iterative method for tilt phase-shift interferometry,” Opt. Express 21(24), 29505–29515 (2013). [CrossRef]  

28. L. L. Deck, “Model-based phase shifting interferometry,” Appl. Opt. 53(21), 4628–4636 (2014). [CrossRef]  

29. W. Swann, “Direct search methods,” Numerical methods for unconstrained optimization 124, 13–28 (1972). [CrossRef]  

30. A. Pyzara, B. Bylina, and J. Bylina, “The influence of a matrix condition number on iterative methods’ convergence,” in 2011 Federated Conference on Computer Science and Information Systems (FedCSIS), (IEEE, 2011), 459–464.

31. K. Madsen and H. B. Nielsen, “Introduction to optimization and data fitting,” (2008), http://www2.imm.dtu.dk/pubdb/edoc/imm5938.pdf, retrieved June 2021, pp. 115-125.

32. Q. Kemao, Windowed fringe pattern analysis (Society of Photo-Optical Instrumentation Engineers, 2013), pp. 203–231.

33. R. Courant and F. John, Introduction to calculus and analysis I (Springer Science & Business Media, 2012), pp. 440–480.

34. J. L. Flores, A. Muñoz, S. Ordoñes, G. García-Torales, and J. A. Ferrari, “Color-fringe pattern profilometry using an efficient iterative algorithm,” Opt. Commun. 391, 88–93 (2017). [CrossRef]  

35. L. L. Deck, “Suppressing phase errors from vibration in phase-shifting interferometry,” Appl. Opt. 48(20), 3948–3960 (2009). [CrossRef]  

36. K. Liu, Y. Wang, D. L. Lau, Q. Hao, and L. G. Hassebrook, “Gamma model and its analysis for phase measuring profilometry,” J. Opt. Soc. Am. A 27(3), 553–562 (2010). [CrossRef]  

37. L. L. Deck and J. A. Soobitsky, “Phase-shifting via wavelength tuning in very large aperture interferometers,” in Optical manufacturing and testing III, 3782(International Society for Optics and Photonics, 1999), 432–442. [CrossRef]  

38. S. Zhang and S.-T. Yau, “Generic nonsinusoidal phase error correction for three-dimensional shape measurement using a digital video projector,” Appl. Opt. 46(1), 36–43 (2007). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

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 (10)

Fig. 1.
Fig. 1. Flowchart of GIA.
Fig. 2.
Fig. 2. Simulation Parameters (a) Background intensity; (b) Phase-shift of second frame; (c) Ground-truth phase; (d) First fringe pattern from Set X.
Fig. 3.
Fig. 3. RMSEs of extracted phase from example with random δ0.
Fig. 4.
Fig. 4. RMSEs of extracted phase from example with regular δ0.
Fig. 5.
Fig. 5. Accuracy of GIA in different conditions (a) RMSEs with respect to different σI; (b) the RMSEs with respect to different F; (c) the RMSEs with respect to different l; (d) RMSEs with respect to different σα; (e) RMSEs with respect to ξ, the change in E.
Fig. 6.
Fig. 6. Sketch of the Fizeau interferometer.
Fig. 7.
Fig. 7. Fringe patterns and ground-truth phase. (a) One of acquired fringe patterns. (b) The ground-truth phase.
Fig. 8.
Fig. 8. Extracted phase and phase-shifts by GIA (unit: rads) (a) Phase; (b) Phase-shift; (c) Tilt removed phase-shift.
Fig. 9.
Fig. 9. Difference between the measured intensity and fringe pattern generated by (a) GIA; (b) AIA; (c) Xu’s algorithm; (d) Hoang’s algorithm; (e) PTI; (f) MPSI.
Fig. 10.
Fig. 10. RMSEs of extracted phase in experiment. (a) with random phase-shifts (b) with regular phase-shifts

Tables (6)

Tables Icon

Table 1. Parameter settings for different algorithms.

Tables Icon

Table 2. Errors simulated in fringe patterns in Group (i)-(x).

Tables Icon

Table 3. RMSEs of extracted phase from example with random δ0 (Unit: rad).

Tables Icon

Table 4. RMSEs of extracted phase from example with regular δ0 (Unit: rad).

Tables Icon

Table 5. Computation time of each algorithm (Unit: second).

Tables Icon

Table 6. RMSEs of extracted phase in experiment (Unit: rad).

Equations (47)

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

I ( i ; j x , y ) = I t ( i ; j x , y ) + n ( i ; j x , y ) ,
I t ( i ; j ) = A ( i ; j ) + k = 1 P B k ( i ; j ) cos [ k φ ( j ) + k δ ( i ; j ) ] ,
δ ( i ; j ) = u = 0 Q v = 0 u α u v ( i ) x u y u v ( N x 1 ) u ( N y 1 ) u v ,
E = i , j [ I t ( i ; j ) I ( i ; j ) ] 2 ,
E P ( j ) = i [ I t ( i ; j ) I ( i ; j ) ] 2 ,
E F ( i ) = j [ I t ( i ; j ) I ( i ; j ) ] 2 .
I t ( i ; j ) = A ( i ; j ) + B ( i ; j ) cos [ φ ( j ) + δ ( i ) ] .
I t ( i ; j ) = a ( j ) + b ( j ) cos [ δ ( i ) ] + c ( j ) sin [ δ ( i ) ] ,
I t ( i ; j ) = a ( i ) + b ( i ) cos [ φ ( j ) ] + c ( i ) sin [ φ ( j ) ] ,
| [ δ m ( i ) δ m ( 0 ) ] [ δ m 1 ( i ) δ m 1 ( 0 ) ] | < ε , for all the frames,
m M ,
I t ( i ; j ) = A ( i ; j ) + k = 1 P B k ( i ; j ) cos [ k φ ( j ) + k δ ( i ) ] .
I t ( i ; j ) = q = 0 2 P X q ( j ) S q ( i ) ,
I t ( i ; j ) = q = 0 2 P X q ( i ) S q ( j ) ,
I t ( i ; j ) = A ( i ) + k = 1 P B k ( i ) cos [ k φ ( j ) + k δ ( i ) ] .
I t ( i ; j ) = A ( j ) + k = 1 P B k ( j ) cos [ k φ ( j ) + k δ ( i ) ] .
I t ( i ; j ) = A ( j ) + B ( j ) cos [ φ ( j ) + δ ( i ; j ) ] = A ( j ) + B ( j ) cos [ φ ( j ) + δ 0 ( i ) + α 1 ( i ) x + β 1 ( i ) y ] ,
I t ( i ; j x , Y ) = A ( j x , Y ) + B ( j x , Y ) cos [ φ ( j x , Y ) + δ 0 Y ( i ) + α 1 ( i ) x ] = A ( j x , Y ) + B ( j x , Y ) cos [ Φ r ( i ; j x , Y ) + p r ( i ) ] ,
I t ( i ; j x , Y ) = A ( j x , Y ) + b r ( i ) B ( j x , Y ) cos [ Φ r ( i ; j x , Y ) ] + c r ( i ) B ( j x , Y ) sin [ Φ r ( i ; j x , Y ) ] ,
I t ( i ; j ) = A ( j ) + k = 1 P B ( j ) ( g ) k 1 cos { k [ φ ( j ) + δ 0 ( i ) + α 1 ( i ) x + β 1 ( i ) y ] } ,
I t ( i ; j ) = A ( j ) + k = 1 P B k ( j ) cos [ k φ ( j ) + k δ ( i ; j ) ] .
A ( j x + s , y + t ) = A ( j x , y ) ,
B k ( j x + s , y + t ) = B k ( j x , y ) ,
x k = { A ( j x , y )  if  k = 0 , B k ( j x , y )  if  k = 1 ,   , P ,  
r ( 2 w + 1 ) × ( 2 w + 1 ) × i + l ( x ) = I t ( i ; j x + s , y + t ) I ( i ; j x + s , y + t ) ,
J ( 2 w + 1 ) × ( 2 w + 1 ) × i + l , k ( x ) = r ( 2 w + 1 ) × ( 2 w + 1 ) × i + l ( x ) x k = cos [ k φ ( j x + s , y + t ) + k δ ( i ; j x + s , y + t ) ] ,
x = φ ( j ) ,
r i ( x ) = I t ( i ; j ) I ( i ; j ) ,
J i , 1 ( x ) = r i ( x ) x = k = 1 P k B k ( j ) sin [ k φ ( j ) + k δ ( i ; j ) ] ,
x u ( u + 1 ) / 2 + v = α u v ( i ) ,
r j ( x ) = I t ( i ; j ) I ( i ; j ) ,
J j , u ( u + 1 ) / 2 + v ( x ) = r j ( x ) x u ( u + 1 ) / 2 + v = k = 1 P k x u y u v B k ( j ) sin [ k φ ( j ) + k δ ( i ; j ) ] ( N x 1 ) u ( N y 1 ) u v ,
| E m + 1 E m | / E m < ζ .
φ ( j ) = ω x x + ω y y ,
φ ( j ) = 2 ω p π × [ ( x N x / 2 ) 2 + ( y N y / 2 ) 2 ] / [ ( N x / 2 ) 2 + ( N y / 2 ) 2 ] + ω x x + ω y y ,
φ ( j ) = ω p × p e a k s _ s ( N x , N y ) + ω x x + ω y y ,
I o ( i ; j ) = A ( j ) + B ( j ) cos [ φ ( j ) + δ 0 ( i ) ] ,
A ( j ) = 0.5 0.25 [ ( x N x / 2 ) 2 + ( y N y / 2 ) 2 ( N x / 2 ) 2 + ( N y / 2 ) 2 ] ,
δ 0 = [ 0 0.440 0.795 1.246 1.601 2.665 2.695 3.449 ] ,
f ( x ) = 1 2 | | r ( x ) | | 2 = 1 2 i = 0 m 1 [ r i ( x ) ] 2 ,
r ( x + h ) r ( x ) + J ( x ) h ,
J i , j ( x ) = r i x j ( x ) .
f ( x + h ) = f ( x ) + h T J T ( x ) r ( x ) + 1 2 h T J T ( x ) J ( x ) h .
J T ( x ) J ( x ) h = J T ( x ) r ( x ) .
[ J T ( x ) J ( x ) + μ I ] h LM = J T ( x ) r ( x ) .
| | J T ( x ) r ( x ) | | < η ,
m L M > M L M ,
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.