- Email: [email protected]
Advanced spatial spectrum fitting algorithm for significantly improving the noise resistance ability of self-calibration phase shifting interferometry
Optics and Lasers in Engineering 112 (2019) 170–181
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Advanced spatial spectrum fitting algorithm for significantly improving the noise resistance ability of self-calibration phase shifting interferometry Shengqi Cao, Yi Wang, Xiaoxu Lu, Liyun Zhong∗ Guangdong Provincial Key Laboratory of Nanophotonic Functional Materials and Devices, South China Normal University, Guangzhou 510006, China
a r t i c l e Keywords: Fringe analysis Interferometry Metrology Phase retrieval Phase measurement
i n f o
a b s t r a c t We propose a novel phase demodulation algorithm for phase shifting interferometry (PSI) named as advanced spatial spectrum fitting (ASSF) algorithm. By utilizing dual-operation of singular value decomposition (SVD) and hybrid temporally and spatially estimating the background intensity, the retrieved phase distribution will become more robust to the noise distribution compared with the traditional PSI. At the same time, with no particular temporal phase shifts distribution requirement, it realize stable self-calibration phase retrieval with few interferograms containing less than one fringe. Therefore, even when the interferograms are seriously corrupted by the random noise or the fringe number is small, the high quality phase retrieval is possible. Both simulation and experimental results will be used to demonstrate the outstanding performance of the proposed ASSF algorithm.
1. Introduction After the tested object being illuminated by a light beam, the spatial phase distribution encoded optical wavefront will change corresponding to the micro-nano structure of the object [1]. Various important quantitative applications of it have been presented, including profile testing [2,3], fluid monitoring [4], refractive index tomography [5,6] and biological imaging [7,8]. Although the direct detection of phase distribution is still impossible due to lack of devices with ultra-high temporal resolution, by encoding the phase information into optical intensity, a time-average quantity, we can achieve it indirectly. To date, some spatial phase encoding and demodulation methods have been reported [9– 14]. Among them the optical interferometry is regarded as an important candidate, in which the phase information is directly encoded into the interference fringe patterns in a linear way, and then the phase demodulation procedure becomes mathematically convenient and deterministic. In current optical interferometry, the off-axis configurations [15,16] are usually utilized to perform dynamic analysis and the phase-shifting configurations [17] are used to make the full utilization of space-bandwidth product of the camera. In phase shifting interferometry (PSI), we need to ensure the accurate estimation of phase shifts for accurate phase retrieval, otherwise, the retrieved phase distribution will suffer from a serious rippled error along the fringe orientations [18] . Initially, people have simplified the computational model by assuming that the phase shifts are known [19]. or the phase-shifter produce equal phase steps which are uniformly distributed in the range of [0, 2𝜋] [20]. Subsequently, it is found that this assump-
∗
tion is not valid in most cases. Under this model, although by increasing the number of captured interferograms we can release the rippled-like error [21,22], more time will be needed, and the vibrational robustness of the interferometric system will be correspondingly reduced. To solve this problem, the self-calibration phase-shifting algorithms are proposed and developed, in which the blind phase shifts are available [23–33] . Some of them are designed to work with 3 or less phase-shifting interferograms to reduce the necessary image capturing time, then some reasonable hypotheses of these interferograms are necessary to balance the equations and unknown parameters in the mathematical model. Generally speaking, in the condition that the interference-induced fringes number and the number of phase-shifting interferograms are both few, or the illumination field is nonuniform, it is very difficult to achieve accurate self-calibration phase retrieval because these hypotheses are no longer valid. This is a significant problem in PSI because the number reduction of interference fringes is usually appeared when we use white light to suppress the coherent artifact [8], and the reduction of the fringes number induced by the aberration is indeed helpful for more sufficient utilization of the space-bandwidth of the camera [34]. Another challenge in PSI is that the phase shifts between interferograms cannot be too small, otherwise the retrieved phase distribution will be greatly sensitive to the noise, which will lead to a restriction during practical measurement. Recently, a method is reported to realize self-calibration phase retrieval from only 3-frame phase-shifting interferograms containing less than one fringe [30]. In this method, the mid-frequency spatial spectrum of the fringe patterns is utilized to determine the phase shifts, so the large amplitude of mid-frequency spatial spectrum of the interfero-
Corresponding author. E-mail addresses: [email protected] (Y. Wang), [email protected] (L. Zhong).
https://doi.org/10.1016/j.optlaseng.2018.09.007 Received 1 June 2018; Received in revised form 8 August 2018; Accepted 10 September 2018 Available online 27 September 2018 0143-8166/© 2018 Elsevier Ltd. All rights reserved.
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
number of phase-sifting interferograms is more than three (N > 3). For 𝑁 ∑ a group of coefficients𝑘1 , 𝑘2 ......𝑘𝑁 , we can define 𝑘𝑛 𝐼𝑛 as a combina-
grams is preferred. However, in the ultra-sparse fringe patterns (USFP) contain less than one fringe, the amplitude of the mid-frequency spatial spectrum is weaker than the fringe patterns which contain many fringes. This lead to the accuracy decrease when we intend to use the USFPs to make phase retrieval. What’s more, in this algorithm the number of interferograms is fixed at 3, and we will clarify in the next section, to guarantee the accuracy, the phase shifts between arbitrary 2 interferograms cannot close to𝜋. Apparently, these restrictions may bring inconvenience in its practical applications. We have also noticed that although this method is not sensitive to the difference of fringe patterns, when the phase shifts are small, the noise will still be amplified. It has been reported that pre-removal of background intensity by spatial filtering instead of temporal subtraction is useful [35] for improving the signal-to-noise ratio (SNR) of phase retrieval in the case that phase shifts are small, but it is also found that the high-pass filtering function will lead to a serious signal distortion [31] due to the overlap of interference and background intensity. We will demonstrate later that the direct high-pass spatial filtering of these interferograms can both reduce the noise and signal amplitude, which indicates that the improvement of noise-resistance ability in Ref [35] is sensitive to the filtering window. In this manuscript, we proposed a phase-shifting algorithm named as advanced spatial spectrum fitting (ASSF) algorithm that can realize self-calibration phase retrieval from few phase shifting interferograms containing less than one fringe with no strict temporal phase shifts requirement and further boost the noise resistance ability in PSI. In the first step, we summarize the PSI algorithm by introducing a combination modes (CM) theory. According to this theory, we optimize the SNR of retrieved signals by SVD. Next, by using an advanced background intensity removal procedure, we achieve the interference signals with further improved SNR. Finally, by replying the SVD and performing parameters estimation using the spatial spectrum characters of the interferograms, we realize phase retrieval with strong self-calibration ability and ultrahigh SNR. Next, we will discuss the proposed algorithm in detail.
𝑛=1
tion mode (CM) for the N-frame interferograms. Apparently, a group of OCC contains three CMs, and in the right hand of Eqs. 2b and 2c, the first and second term represents the interference signal and the noise, respectively. Apparently, the SNR of CMs containing in different OCCs may be different as well. Our first goal is to search for the CMs that have maximum SNR. If the phase shifts are known, the CMs with the maximum SNR can be easily achieved by temporal least square fitting algorithm. However, the phase shifts cannot be determined precisely before calculation in practical. In this situation, we need a method to realize the blind optimization of SNR. Recently, researchers tend to extract the sub-vectors of the correlation matrix formed by the intensities of interferograms by using singular value decomposition (SVD) operation to achieve high SNR phase retrieval [26,32,36-38]. As we will discuss later, out of the perspective of the combination modes theory, these methods cannot always help us to achieve the optimum SNR. This is because that without optimization, the temporal average subtraction for background intensity A(x, y) removal operation embedded in these methods are blind. The optimum SNR can be achieved by this procedure only when the phase shifts are distributed in [0, 2𝜋] in a quasi-linear way. Only in this case, the temporal average of N-frame interferograms is accidentally corresponding to the first CM in the optimum OCC, i.e., 𝑎1 𝐼1 + .....𝑎𝑁 𝐼𝑁 . Before beginning the description of the proposed method, we briefly review the basic properties of the SVD operation. The SVD operation transforms a real and symmetric correlation matrix A constructed from N input data vectors to a diagonal matrix, and then output N eigenvectors and their eigenvalues. If N input data vectors are I1 , I2 ......IN , then the elements of matrix A can be described as ∑ 𝐴𝑚𝑛 = 𝐼𝑚 (𝑥, 𝑦) ⋅ 𝐼𝑛 (𝑥, 𝑦) 𝑚, 𝑛 = 1......𝑁. (3) 𝑥,𝑦
2. Methods
Two important properties of SVD operation used in this work are listed as follow
2.1. Combination modes of PSI and its optimization
a) The matrix A has N eigenvectors with length of N if the number of interferograms is N, i.e., 𝐶 𝐶𝑚 = (𝑘1𝑚 , 𝑘2𝑚 ......𝑘𝑁𝑚 ), 𝑚 = 1......𝑁. For these eigenvectors we have
First, we introduce a general theory to describe PSI algorithms named as combination mode theory. In PSI, the intensity of the nth phase-shifting interferogram can be described as 𝐼𝑛 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝜃𝑛 ] + 𝑜𝑛 (𝑥, 𝑦),
𝑛 = 1......𝑁.
(𝑘1𝑚 )2 + ...... + (𝑘𝑁𝑚 )2 = 1
(1)
where x and y denotes the transverse and longitude coordinates, respectively; A(x, y) and B(x, y) represents the background intensity and modulation amplitude, respectively; 𝜙(x, y)is the encoded phase distribution; 𝜃 n represents the phase shifts and 𝑛 = 1......𝑁are the index numbers of phase-shifting interferograms; on (x, y)denotes the noise distribution in the nth interferogram. Noticed that the goal of PSI algorithms are in fact to determine a set of orthogonal combination coefficients (OCC) for each interferograms to temporally combine them into B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y)as 𝑁 ∑ 𝑛=1 𝑁 ∑ 𝑛=1 𝑁 ∑ 𝑛=1
𝑎𝑛 𝐼𝑛 = 𝐴(𝑥, 𝑦) +
𝑁 ∑ 𝑛=1
𝑎𝑛 𝑜𝑛 .
𝑏𝑛 𝐼𝑛 = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) +
𝑐𝑛 𝐼𝑛 = 𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦) +
𝑛 𝑘𝑚 𝑘𝑛 + 𝑘𝑚 𝑘𝑛 + ...... + 𝑘𝑚 𝑁 𝑘𝑁 = 0 1 1 2 2
𝑏𝑛 𝑜𝑛 .
(2b)
𝑐𝑛 𝑜𝑛 .
(2c)
𝑛=1 𝑁 ∑ 𝑛=1
𝑚, 𝑛 = 1......𝑁.
(4a)
(4b)
(𝑘1𝑚 𝐼1 + 𝑘2𝑚 𝐼2 + ...... + 𝑘𝑁𝑚 𝐼𝑁 )(𝑘1𝑛 𝐼1 + 𝑘2𝑛 𝐼2 + ...... + 𝑘𝑁𝑛 𝐼𝑁 ) = 0 𝑚, 𝑛 = 1......𝑁.
(4c)
b) If there is an arbitrary vector with length of N, i.e., 𝐶 𝐶 ′ = (𝑐1 ′, 𝑐2 ′......𝑐𝑁 ′), and satisfy the condition of 𝑐 ′21 + 𝑐 ′22 + ...... + 𝑐′2𝑁 = 1, the function 𝑇 = (𝑐 ′1 𝐼1 + 𝑐 ′2 𝐼2 + ...... + 𝑐 ′𝑁 𝐼𝑁 )2 will reach its first maximum 𝜆1 in the case that 𝐶 𝐶 ′ = 𝐶 𝐶1 , then the second maximum 𝜆2 in the case that 𝐶 𝐶 ′ = 𝐶 𝐶2 and so on (𝜆1 > 𝜆2 > ...... > 𝜆N ).
(2a) 𝑁 ∑
𝑚 = 1......𝑁.
Each interferogram will contain three linear and independent components, i.e. A(x, y) B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y). If the spatial noise distributions 𝑜𝑛 (𝑥, 𝑦), 𝑛 = 1......𝑁of every interferograms have the same variance 𝜎 2 and are random and independent with each other, then according to the statistical theory, the variance of noise fluctuations in Eqs .(2b) and (2c) are (𝑏1 2 + 𝑏2 2 + ......𝑏𝑁 2 )𝜎 2 and (𝑐1 2 + 𝑐2 2 + ......𝑐𝑁 2 )𝜎 2 respectively. If we define the SNR of two interference signals i.e., B(x, y)cos 𝜙(x, y) and B(x, y)sin 𝜙(x, y), obtained by Eq. (2) as ∑ 2 𝐵 (𝑥, 𝑦)cos2 [𝜙(𝑥, 𝑦) + 𝜃𝑛 ]
Here, we definea1 ......aN , b1 ......bN and c1 ......cN as a group of OCC. If the OCC is known, we can achieve 𝜙(x, y)by performing the arctangent operation to the division ofB(x, y)sin 𝜙(x, y)andB(x, y)cos 𝜙(x, y). The OCC are directly related with the phase shifts, but it is not unique if the
𝑆𝑁 𝑅1 = 171
𝑥,𝑦
(𝑏1 2 + 𝑏2 2 + ......𝑏𝑁 2 )𝜎 2
.
(5a)
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
∑ 𝑆𝑁 𝑅2 =
𝑥,𝑦
As discussed in Section 2.1, to achieve accurate interference signals, it is needed to perform the background removal in PSI. Traditionally, this procedure can be realized by the temporal subtraction operation, but both the effective phase shifts and modulation amplitude may be reduced after this, and thus the noise level is amplified. The background removal can also be realized using spatial convolution, and this idea has directly given birth to the 2-step PSI. Recently, it is found that in the case that phase shifts are small, the noise amplification induced by the temporal subtraction background removal can be released by using spatial convolution instead [35]. To clarify its principle briefly, we let N = 3 and 𝜃1 = 0, 𝜃2 = 𝜋6 , 𝜃3 = 𝜋3 in Eq. (1), then we have
𝐵 2 (𝑥, 𝑦)sin2 [𝜙(𝑥, 𝑦) + 𝜃𝑛 ]
(𝑐1 2 + 𝑐2 2 + ......𝑐𝑁 2 )𝜎 2
.
(5a)
Then apparently, according to the properties of SVD above mentioned, 3 CMs with the maximum SNR, in which A(x, y), B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y)are contained in them should correspond to the first 3 eigenvectors or their orthogonal transformations. In the other CMs, we can easily understand that they contain the redundant components that can be described as 𝑁 ∑ 𝑛=1
𝑐𝑛 𝐼𝑛 = 0 +
𝑁 ∑ 𝑛=1
𝑐𝑛 𝑜𝑛 .
(6)
𝐼1 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝑜1 (𝑥, 𝑦).
In which only the noise terms are left. The SVD operation is in fact to help us excluding these redundant CMs. After the 3 optimized CMs with the same noise level (guaran𝑁 ∑ teed by the property a) of SVD operation), i.e., 𝑎1 𝐼1 + ...... + 𝑎𝑁 𝐼𝑁 , 𝑁 ∑ 𝑛=1
𝑏1 𝐼1 + ...... + 𝑏𝑁 𝐼𝑁 ,
𝑁 ∑ 𝑛=1
𝑛=1
𝑎𝑚
𝑛=1
𝑎𝑛
−
𝑐𝑚
𝑁 ∑
𝑛=1
𝜋 ] + 𝑜3 (𝑥, 𝑦). 3
(8c)
𝜋 + ] + 𝑜1 (𝑥, 𝑦) − 𝑜3 (𝑥, 𝑦). 6
(7b)
𝑎𝑚
𝑁 ∑
𝑎𝑛
−
𝑏𝑚
𝑁 ∑
𝑛=1
𝜋 𝐵(𝑥, 𝑦) sin[𝜙(𝑥, 𝑦) 12
𝜋 ] + 𝑜1 (𝑥, 𝑦) − 𝑜2 (𝑥, 𝑦). 12
𝐼2′ (𝑥, 𝑦) = 𝐼3 (𝑥, 𝑦) − 𝐼1 (𝑥, 𝑦) = 2 sin
𝑛=1
𝑁 ∑
𝐼3 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
+
𝑁 ∑
𝑚=1
𝑚=1
(8b)
(7a)
√ √𝑁 √∑ according to the Eq. (2). Here, two denominators √ ( √ √ √𝑁 √∑ and √ ( √
𝜋 ] + 𝑜2 (𝑥, 𝑦). 6
𝐼1′ (𝑥, 𝑦) = 𝐼1 (𝑥, 𝑦) − 𝐼2 (𝑥, 𝑦) = 2 sin
√ √ 2 √ ⎛ ⎞ √ 𝑎𝑛 𝐼𝑛 𝑐𝑛 𝐼𝑛 √ 𝑁 ⎜ ⎟ √ 𝑐 𝑛=1 𝑛=1 √ ∑ ⎜ 𝑎𝑚 𝑃2 = ( − )∕√ − 𝑚 ⎟ . √ 𝑁 𝑁 𝑁 𝑁 ⎜ ⎟ ∑ ∑ ∑ ∑ √𝑚=1 ⎜ 𝑎𝑛 𝑐𝑛 𝑎 𝑐 ⎟ ⎝ 𝑛=1 𝑛 𝑛=1 𝑛 ⎠ 𝑛=1 𝑛=1 𝑁 ∑
𝐼2 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
In which the variance of o1 (x, y), o2 (x, y), o3 (x, y) is equal to 𝜎 2 . The phase shifts of Eq. (8) are apparently small. If we remove A(x, y) by subtracting I1 (x, y) from I2 (x, y) and I3 (x, y), we have
𝑐1 𝐼1 + ...... + 𝑐𝑁 𝐼𝑁 are obtained, we can sub-
tract the A(x, y) components from them by √ √ 2 𝑁 𝑁 ∑ ⎛∑ ⎞ √ ⎛ ⎞ √ 𝑎 𝐼 𝑏 𝐼 𝑁 𝑛 𝑛 𝑛 𝑛⎟ √ ⎜ ⎜ ⎟ √ ∑ 𝑎𝑚 𝑏 𝑛=1 𝑛=1 ⎟∕ √ ⎜ 𝑃1 = ⎜ − − 𝑚 ⎟ . √ 𝑁 𝑁 𝑁 𝑁 ⎜ ∑ ⎟ √ ⎜ ⎟ ∑ ∑ ∑ √𝑚=1 ⎜ ⎜ 𝑎 𝑏𝑛 ⎟ 𝑎 𝑏 ⎟ ⎝ 𝑛=1 𝑛 ⎠ ⎝ 𝑛=1 𝑛 𝑛=1 𝑛 ⎠ 𝑛=1
(8a)
(9a)
𝜋 𝐵(𝑥, 𝑦) sin[𝜙(𝑥, 𝑦) 6 (9b)
The variance of noise terms are increased from 𝜎 2 in Eq. (8) to 2𝜎 2 in Eq. (9) according to the statistical theory. the modulation amplitude 𝜋 is reduced from B(x, y) to 2 sin 12 𝐵(𝑥, 𝑦) and 2 sin 𝜋6 𝐵(𝑥, 𝑦); the effective 𝜋 𝜋 𝜋 phase shifts are decreased from 6 and 𝜋3 to 𝜋6 − 12 = 12 . However, if A(x, y) is instead removed by the spatial filtering and there is no spatial distortion of the interference terms by this filtering function, which is an ideal condition, then we have
2
)
𝑏𝑛
2
) are employed to normalize the noise variance
𝐼1′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝑜1 (𝑥, 𝑦).
𝑐𝑛
of Eqs. (7a) and (7b). Then we can use self-calibration algorithms [25, 26, 29] to obtain the two interference terms, i.e., B(x, y)sin 𝜙(x, y)and B(x, y)cos 𝜙(x, y) from P1 and P2 and the spatial phase distribution can be obtained by the antitangent of their division. In comparison, due to the disadvantage of background removal procedure, the sub-space optimization method used before are not in the best way because the blind background removal by temporal subtraction may include the redundant CMs, for example, the temporal average subtraction in Ref [32] and the subtraction of the first interferograms in Ref [26] . To clarify this, 30-frame phase shifting interferograms are captured and one of them is shown in Fig. 1a. Next, we use these interferograms and advanced iterative algorithm (AIA) method to demodulate the encoded phase and regard it as the ground truth., Then, we choose first 5 phase shifting interferograms to make phase retrieval. Fig. 1c and d respectively present the corresponding phase errors by using background removal procedure adopted by APCA method [26] in which the first interferograms is subtracted by another, and the proposed CMs optimization operation. The comparison between Fig. 1c – e have clarified the importance of CM optimization.
(10a)
𝐼2′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
𝜋 ] + 𝑜2 (𝑥, 𝑦). 6
(10b)
𝐼3′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
𝜋 ] + 𝑜3 (𝑥, 𝑦). 3
(10c)
Compared with Eq. (9), the solution of Eq. (10) is apparently with higher SNR because the larger modulation amplitude, smaller noise variance and increase of equations number. However, different with the ideal condition, there is a spatial phase reconstruction effect if the spatial filtering operation is directly applied to the interferograms. In following analysis, for simplicity, we will neglect the y dependence. In general, if a spatial filtering operation is utilized, corresponding to a convolution operation (H(x)⊗) in spatial domain or a multiply operation (H(fx ) · ) in frequency domain, the distortion of the interference term can be expressed as [31] ∑ 𝐼̃𝑛 (𝑥) = 𝐻(𝑥) ⊗ 𝐵(𝑥) cos[𝜙(𝑥) + 𝛿𝑛 ] = 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos[𝜙(𝑥 − 𝜏) + 𝛿𝑛 ] 𝜏 ∑ = 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos[𝜙(𝑥) + 𝑓 (𝑥, 𝜏) + 𝛿𝑛 ] 𝜏 ∑ = cos[𝜙(𝑥) + 𝛿𝑛 ] 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏) − sin[𝜙(𝑥) 𝜏 ∑ +𝛿𝑛 ] 𝐻(𝜏)𝐵(𝑥 − 𝜏) sin 𝑓 (𝑥, 𝜏). √ 𝜏 √ = 𝑇 (𝑥) cos[𝜙(𝑥) + 𝛿𝑛 + 𝜂(𝑥)] = 𝑇 (𝑥) cos[𝜙′ (𝑥) + 𝛿𝑛 ]
2.2. Advanced noise suppression strategy Out of the traditional perspective, after three optimized CMs are retrieved, we have already achieved the signals with the highest SNR. However, when the phase shifts are small, the error induced by the noise may still too large even the optimized CMs are used. In these cases, we need a more efficient noise suppression strategy.
(11) in which 𝑓 (𝑥, 𝜏) = 𝜙(𝑥 − 𝜏) − 𝜙(𝑥) 172
(12a)
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 1. (a) Interferogram of a single macrophage cell; (b) The retrieved phase with aberration phase removed using 30-frame phase shifting interferograms and AIA algorithm which is regarded as the ground truth; The errors obtained with (c) the reported APCA background removal procedure; (d) the proposed CMs optimization operation. (e) The 200th row of (c) and (d).
𝑇 (𝑥) = [
∑
𝐻(𝜏)𝐵(𝑥 − 𝜏) sin 𝑓 (𝑥, 𝜏)]2 + [
∑
𝜏
be utilized only in the harsh condition with very small phase shifts. When the phase shifts become large or the noise level is moderate, the error induced by this distortion will become larger relative to the error induced by noise and the reduction of modulation amplitude may lead to an overall decrease of SNR instead of increase. Next, we will demonstrate an advanced strategy to remove A(x, y) without distorting the phase distribution or amplifying the noise level. Due to there are 3 components in the interferograms, i.e. A(x, y), B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y), there are 3 optimum CMs corresponding to the first 3 eigenvectors of the correlation matrix as follow
𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏)]2
𝜏
(12b) ∑ 𝜂(𝑥) = arccos(
𝜏
𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏) √
).
𝑇 (𝑥)
(12c)
𝜙′ (𝑥) = 𝜙(𝑥) + 𝜂(𝑥).
(12d)
From Eq. (12d), we can observe that the phase distribution is changed from the original 𝜙(x)to 𝜙′(x)with an addictive distortion term𝜂(x), and this distortion will bring in a serious error especially in the case that the fringes density is relative small. It is worth noting that although the interference terms have more extended spatial spectrum relative to the background term, they are both centered at the zero frequency and thus overlapped. That is to say, when high-pass filtering operation is applied to filter out the background, it will also reduce the ∑√ ∑ amplitude of interference term, i.e. 𝑇 (𝑥) < 𝐵(𝑥), and contribute 𝑥
𝑇1 (𝑥, 𝑦) = 𝐶 𝑀1 = 𝑘11 𝐼1 + 𝑘12 𝐼2 + ......𝑘1𝑁 𝐼𝑁 𝑁 ∑ = 𝑡1 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿1 ] + 𝑘1𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13a)
𝑇2 (𝑥, 𝑦) = 𝐶 𝑀2 = 𝑘21 𝐼1 + 𝑘22 𝐼2 + ......𝑘2𝑁 𝐼𝑁 𝑁 ∑ = 𝑡2 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿2 ] + 𝑘2𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13b)
𝑇3 (𝑥, 𝑦) = 𝐶 𝑀3 = 𝑘31 𝐼1 + 𝑘32 𝐼2 + ......𝑘3𝑁 𝐼𝑁 𝑁 ∑ = 𝑡3 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿3 ] + 𝑘3𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13c)
𝑛=1
𝑛=1
𝑥
to the reduction of SNR consequently. Therefore, the above method can
𝑛=1
173
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 2. Spatial spectrums of three CMs, in which the background and interference terms are mixed. (f2 − , f 1 − ) and (f 1 + , f 2 + ) denote the spatial frequency band where mainly occupied by the interference term. (f 1 − , f 1 + ) denotes the frequency band where mainly occupied by background intensity A(x,y).
The basic idea is that, if we want to filtered out A(x, y) without the interference terms distortion, we need to avoid the overlap of them. There is a beneficial property can be helpful, which is, in PSI, the phase contrast is mainly induced by the spatial phase difference between the reference and object wave. That is to say, the background A(x, y) and modulation amplitude B(x, y) changes spatially slower than B(x, y)cos 𝜙(x, y) and B(x, y)sin 𝜙(x, y). As shown in Fig. 2, we can separate the spatial spectrum of the fringe patterns into three parts: (f2 − , f 1 − ) and (f 1 + , f 2 + ) are mainly occupied by the interference terms and (f 1 − , f 1 + ) is mainly occupied by A(x,y). Next, we will use this spectrum dividing properties of fringe patterns to separate A(x,y) and the interference terms. It is easy to understand that there exists k1 and k2 that satisfy 𝑘1 𝑇1 (𝑥, 𝑦) + 𝑘2 𝑇2 (𝑥, 𝑦) + 𝑇3 (𝑥, 𝑦) = 𝐿0 (𝑥, 𝑦)
and Eq. (14) as 𝐴(𝑥, 𝑦) =
Here, 𝐿̃ 0 (𝑥, 𝑦)denotes the low-pass function of L0 (x, y). By subtracting A(x, y) described by Eq. (16) from every interferograms described in Eq. (1), we have that 𝐼̃𝑛 (𝑥, 𝑦) = 𝐼𝑛 (𝑥, 𝑦) − 𝐴(𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝜃𝑛 ] + 𝑜𝑛 (𝑥, 𝑦).
(17)
Clearly, there are only two independent components in 𝐼̃𝑛 (𝑥, 𝑦)left now. Then, we can perform the procedure of SVD again, and the correlation matrix elements are changed from Eq. (3) to ∑ 𝐴̃ 𝑚𝑛 = 𝐼̃𝑚 (𝑥, 𝑦) ⋅ 𝐼̃𝑛 (𝑥, 𝑦) 𝑚, 𝑛 = 1......𝑁. (18)
(14)
𝑥,𝑦
In which the interference terms is compensated. Here, L0 (x, y) is a quasi-constant in spatial domain and be proportional to A(x, y). Now, we intend to determine the values of k1 and k2 . It has demonstrated [30] that by using a low-frequency suppression spatial filter, the accuracy of A(x, y) subtraction is improved relative to the direct spatial average subtraction. Therefore, we firstly apply the high-pass filtering operation to both sides of Eq. (14) and then L0 (x, y) can be well removed, and we have 𝑘1 𝑇̃1 (𝑥, 𝑦) + 𝑘2 𝑇̃2 (𝑥, 𝑦) + 𝑇̃3 (𝑥, 𝑦) = 0
𝐿̃ 0 (𝑥, 𝑦) 𝑘1 (𝑘11 + ...... + 𝑘1𝑁 ) + 𝑘2 (𝑘21 + ...... + 𝑘2𝑁 ) + (𝑘31 + ...... + 𝑘3𝑁 ) (16)
Thus, 2 instead 3 CMs with the highest SNR are achieved as 𝑇̃1 (𝑥, 𝑦) = 𝑘̃ 11 𝐼1 + 𝑘̃ 12 𝐼2 + ......𝑘̃ 1𝑁 𝐼𝑁 = 𝑡1 ′𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿1 ′]
(19a)
𝑇̃2 (𝑥, 𝑦) = 𝑘̃ 21 𝐼1 + 𝑘̃ 22 𝐼2 + ......𝑘̃ 2𝑁 𝐼𝑁 = 𝑡2 ′𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿2 ′]
(19b)
In which A(x, y) is removed and only the interference intensities are left. Till now, the effective phase shifts and modulation amplitude in Eq. (19) has been optimized based on the subtraction operation of A(x, y) through using the spatial property of interferogram in advance. The flow chart of the noise-suppressing strategy described in this section is given in Fig. 3 as follow Unfortunately, we cannot achieve the phase distribution due to 𝑡1 ′ , 𝑡2 ′ , 𝛿2 ′ − 𝛿1 ′ are still unknown. Therefore, it is needed to perform a self-calibration procedure for accurate spatial phase extraction. Next, we will demonstrate a rapid and accurate self-calibration algorithm using Eq. (19).
(15)
To solve the spatial least square problem (𝑘1 , 𝑘2 ) = ∑ 2 min{ [𝑘1 𝑇̃1 (𝑥, 𝑦) + 𝑘2 𝑇̃2 (𝑥, 𝑦) + 𝑇̃3 (𝑥, 𝑦)] }, we can determine k1 and k2 . 𝑥,𝑦
Then, we back substitute k1 and k2 to Eq. (14) to achieve L0 (x, y). The high-frequency spatial spectrum in the obtained L0 (x, y) in this step is corresponding to the noise fluctuations and the residual interference terms. Next, by utilizing a low-pass filtering operation to L0 (x, y) to obtain A(x, y) instead of direct high-pass filtering operation for each interferogram [35], we can avoid the distortion of interference terms. Another advantage of this procedure is that the filtering procedure is performed after the compensation of interference terms, so the final resolution will not be reduced, which is a significant disadvantage in traditional noise suppression method with low-pass filtering [39]. The relationship between A(x, y) and L0 (x, y) can be determined by Eq. (13)
2.3. Self-calibration phase retrieval from two optimum CMs We know that not only A(x, y) but also B(x, y) are quasi-constant in spatial domain. To sufficiently use both of these two properties, we will separate B(x, y) from the interference terms firstly. Neglecting the 174
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 3. The flow chart of the advanced noise suppressing strategy.
in which 𝐾̃ 1 , 𝐾̃ 2 and 𝐾̃ 3 denotes the high-pass filtering results of K1 , K2 and K3 , and M0 (x,y) is eliminated by the filtering operation. By solving the linear least-square problem, i.e. (𝑙1 , 𝑙2 ) = ∑ 2 min[ |𝑙1 𝐾̃ 1 (𝑥, 𝑦) + 𝑙2 𝐾̃ 2 (𝑥, 𝑦) + 𝐾̃ 3 (𝑥, 𝑦)| ], we can achieve l1 and l2 . Con-
piston phase shifts, we have 𝛿1 ′ = 0, 𝛿2 ′ = 𝛿, and then 𝐾1 (𝑥, 𝑦) = 𝑇̃1 ⋅ 𝑇̃2 = 𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿] = =
𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥,𝑦) {cos[2𝜙(𝑥, 𝑦) + 𝛿] + cos 𝛿} 2 𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥,𝑦) [cos 2𝜙(𝑥, 𝑦) cos 𝛿 − sin 2𝜙(𝑥, 𝑦) sin 𝛿 2
𝐾2 (𝑥, 𝑦) = 𝑇̃12 = 𝑡1 ′2 𝐵 2 (𝑥, 𝑦)cos2 𝜙(𝑥, 𝑦) =
(20a)
𝑥,𝑦
+ cos 𝛿].
sequently, the accuracy of the blind calibration is directly determined by the accuracy of the achieved l1 and l2 . Compared with the previous MSSM method [30], we can conclude the advantages of proposed algorithm as following: First, if the phase shifts is 𝜋, the MSSM algorithm is invalid, but in this method, this restriction can be canceled by adding a spatial quasi-constant assumption of B(x, y) compared with Ref [30]. Second, the interference terms B2 (x, y)cos 2𝜙(x, y)and B2 (x, y)sin 2𝜙(x, y) in Eq .(13) have stronger mid-frequency amplitude than the interference terms B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y) used in MSSM algorithm, this will bring stronger blind calibration ability for this algorithm. Thus, in the case of USFP, it will reveal more obvious advantage. Finally, the MSSM algorithm is specific for the case of 3 phase shifting interferograms are captured and this method can make full use of the N-frame (N ≥ 3) interferograms.
𝑡1 ′2 𝐵 2 (𝑥, 𝑦) [cos 2𝜙(𝑥, 𝑦) + 1]. 2 (20b)
𝐾3 (𝑥, 𝑦) = 𝑇̃22 = 𝑡2 2 𝐵 2 (𝑥, 𝑦)cos2 [𝜙(𝑥, 𝑦) + 𝛿] = =
𝑡1 2 𝐵 2 (𝑥,𝑦) {cos[2𝜙(𝑥, 𝑦) + 2𝛿] + 1} 2 𝑡1 ′ 2 𝐵 2 (𝑥,𝑦) [cos 2𝜙(𝑥, 𝑦) cos 2𝛿 − sin 2𝜙(𝑥, 𝑦) sin 2𝛿 2
(20c) + 1].
There are three components B2 (x, y), B2 (x, y)cos 2𝜙(x, y)and B2 (x, y)sin 2𝜙(x, y) in K1 , K2 and K3 . In general, B2 (x, y) is a quasi-constant in spatial domain, then there will exist l1 and l2 that satisfy 𝑙1 𝐾1 (𝑥, 𝑦) + 𝑙2 𝐾2 (𝑥, 𝑦) + 𝐾3 (𝑥, 𝑦) = 𝑀0 (𝑥, 𝑦). in which M0 (x, y) is proportional to (𝑇̃2 +
B2 (x,
𝑙1 𝑙 2 𝑇̃1 )2 + (𝑙2 − 1 )𝑇̃12 = 𝑃0 (𝑥, 𝑦). 2 4
(21) 3. Simulations and discussions
y). Then we have
Following, we will present the performance of proposed ASSF algorithm through using two simulations. The size of simulated fringe pattern is 500 × 500 pixels with pixel size of 0.01 mm. 5-frame of fringe patterns with phase shifts of 0 rad, 0.2 rad, 0.3 rad, 0.4 rad, 0.5 rad are generated for each simulation; In the first simulation, the phase model that induced several fringes is set as 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ). In the second simulation, the phase model induced few (less than one) fringe is set as 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500); the background and modulation intensity is A(x,y) = 50exp[− 0.2(x2 + y2 )] + 50 and B(x,y) = 30exp[− 0.2(x2 + y2 )] + 30, respectively. The SNR of all these fringe patterns are fixed at 50 dB. Fig. 4 shows the results of phase retrieval achieved with the proposed ASSF, MSSM [30], APCA [26] and PCAN [35] algorithms. Note that the noise fluctuation in fringe patterns will lead to a serious error when the phase shifts are small. Even when the fringes number is large enough, the APCA and MSSM algorithms only can achieve noisy phase distribution while the ASSF algorithm can achieve accurate one, in which both the rippled error induced by the miss-calibration of phase shifts and the error induced by the random noise can be well-suppressed. Although the PCAN algorithm can achieve smooth profile of phase dis-
(22)
P0 (x, y) is also a quasi-constant in spatial domain. If the piston phase shifts is neglected, the first and second terms in Eq. (22) are corresponding toB(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y), respectively. Apparently, if l1 and l2 are determined, we have 𝑇̃2 +
𝑙1 𝑇̃ = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦). 2 1
(23a)
√ 𝑙2 −
𝑙1 2 𝑇̃ = 𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦). 4 1
(23b)
And then 𝜙(x, y)can be determined by 𝜙(𝑥, 𝑦) = arctan[
𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦) ]. 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦)
(24)
To achieve l1 and l2 , we perform the spatial-high pass filtering operation to both two sides of Eq. (21), then 𝑙1 𝐾̃ 1 (𝑥, 𝑦) + 𝑙2 𝐾̃ 2 (𝑥, 𝑦) + 𝐾̃ 3 (𝑥, 𝑦) = 0.
(25) 175
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 4. (a1) (a2) the theoretical phase distribution of 𝜙1 (x, y) and one of the 5-frame simulated phase-shifting interference patterns, respectively; the spatial phase distributions of 𝜙1 (x, y) achieved with different algorithms: (a3) MSSM (a5)APCA (a7)ASSF (a9) PCAN; (a4) (a6) (a8) (a10) is the differences between 𝜙1 (x, y) and (a3) (a5) (a7) and (a9), respectively; (b1) (b2) the theoretical phase distribution of 𝜙2 (x, y)and one of 5-frame simulated phase-shifting interference patterns, respectively; the spatial phase distributions 𝜙2 (x, y)achieved by different algorithms: (b3) MSSM (b5) APCA (b7) ASSF (b9) PCAN; (b4) (b6) (b8) (b10) is the differences between 𝜙1 (x, y) and (b3) (b5) (b7) and (b9), respectively.
tribution which indicates less noisy result, its phase distortion becomes very serious, especially in the center area of interference pattern with sparse fringe, as shown in Fig. 4(a10). In the second simulation, the fringe number is decreased and the influence of noise fluctuation is still large due to the small phase shifts. We can see that only the proposed ASSF algorithm is available in this case. The small range of the phase distribution has violated the prior-assumptions of APCA algorithm and the sparsity of the fringe patterns makes the high-pass filtering of the interferograms be invalid in PCAN algorithm. The strong noise also has overwhelmed the achieved spatial phase distribution except with ASSF algorithm. These results prove that the proposed ASSF algorithm can
be utilized at the most difficult situations in PSI, in which the fringes number is small and the phase shifts are small at the same time. Note that in both the ASSF algorithm and the PCAN algorithm, the improvement of noise-resistance ability is dependent on the quality of background removal, which is directly related with the accuracy of phase retrieval. As shown in Fig. 5, the A(x, y) achieved by the proposed ASSF algorithm reveals obvious accuracy advantage compared with the PCAN algorithm. Like the ASSF algorithm [30], the proposed MSSM algorithm also can work well in the case that the fringes number is less than one. However, when the phase shifts between arbitrary two interferograms
176
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 5. (a) theoretical spatial distribution of A(x, y); A(x, y) achieved by the proposed ASSF algorithm in the case that (b) 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ); (c) 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500); (d) (e) the differences between (b) (c) and (a), respectively; A(x, y) achieved by the PCAN algorithm in the case that (f) 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ); (g) 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500), respectively; (h) (i) the differences between (f) (g) and (a), respectively.
is equal to𝜋, the MSSM algorithm becomes invalid if the total number of interferograms N = 3. This is because in this case, the spatial quasiconstant assumption of A(x, y) can no longer support the determination of phase shifts mathematically in the MSSM algorithm. However, in the ASSF algorithm, by introducing an another assumption, i.e. the slowly spatial variation of B(x, y), we can get rid of the above restriction. Fig. 6 presents the RMSEs of phase retrieval achieved with the ASSF and MSSM algorithms when the phase shifts are changed from 2.9 rad to 3.3 rad. According to Fig. 6, although the ASSF and MSSM algorithm can both realize self-calibration phase retrieval from such a sparse fringe pattern, we can clearly observe a sharp peak near 𝜋 in the blue curve, which corresponds to the invalidity of the MSSM algorithm when the phase shift is close to 𝜋. The green curve changes slow with different phase
Fig. 6. The RMSE of phase retrieval using MSSM and ASSF when the phase shifts near of 𝜋.
Fig. 7. The accuracy of phase retrieval using different algorithm described by the corresponding RMSEs with different (a) phase shifts; (b) SNR of fringe patterns; (c) fringe density; (d) variation rate of A(x, y) and B(x, y). 177
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 8. Experimental interferogram encoded with a circular spatial phase distribution; (b) the reference phase achieved by the temporal Fourier transform algorithm; the spatial phase distributions achieved by different algorithms (c) APCA; (d) MSSM; (e) ASSF; (f) (g) (h) the differences between (b) and (c) (d) (e), respectively;.
The phase model is given by 𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500)+0.2(𝑥2 + 𝑦2 ), the phase shifts between the adjacent interferograms are 0.2 rad and the SNR of the interferograms are 50 dB, respectively. When t changes from 0.01 to 0.3, the RMSEs are plotted in Fig. 7d. The results of the PCAN algorithm are not shown in these curves because it cannot output a reasonable phase distribution under some values of these parameters due to the serious distortion of the interference terms by the filtering bank. Clearly, from Fig. 7, we can see several significant advantages of proposed ASSF algorithm: First, we can see that if k < 0.15 rad, which represents small phase shifts, the accuracy of phase retrieval with the ASSF algorithm is almost unchanged while the results with the APCA and MSSM algorithms are greatly reduced (Fig. 7a). Second, it has no strict requirement of the fringe density (Fig. 7c). Third, although the accuracy with the ASSF algorithm is decreased with the SNR of fringe pattern (Fig. 7b), but the noise-suppressing strategy can work in each noise level, so the corresponding result is more reliable compared with APCA and MSSM algorithms. Specially, when t value is large, which indicates that A(x, y) and B(x, y) are nonuniform, the accuracy of ASSF algorithm decreases slower compared with the MSSM algorithm due to the former has stronger mid-frequency spatial spectrum (Fig. 7d).
shifts, indicates that the ASSF algorithm has no such a restriction. Next, we will present the influence of the phase shifts, noise level, fringe density and uniformity of A(x, y), B(x, y) on the accuracy of phase retrieval with the ASSF algorithm. First, we let A(x, y) = 80exp[− 0.05(x2 + y2 )] and B(x, y) = 80exp[− 0.05(x2 + y2 )], 𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500), the SNR is fixed at 50 dB, the phase shift between the adjacent interferograms are set as k rad and 5-frame interferograms are generated. The RMSEs are plotted when k changes from 0.1 to 0.4, as Fig. 7a shown. Next, we let the phase shift between the adjacent interferograms be 0.2 rad, let the SNR of the fringe patterns changes from 35 dB to 65 dB and the RMSEs are plotted in Fig. 7b. Furthermore, we change the phase model to𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500)+𝑚(𝑥2 + 𝑦2 ), in which an addictive quadratic phase component is added. In practical, this quadratic phase component is often induced when the optical system is asymmetric in the reference and objective arms, and consequently increase the density of the fringes. We change the value of m from 0 to 3 and keep the phase shifts between the adjacent interferograms be 0.2 rad and SNR of 50 dB. The RMSEs are plotted in Fig. 7c. At last, we know that in most cases the background and modulation intensity are both nonuniform, to investigate the accuracy robustness of the proposed ASSF algorithm under different spatial changing rate of the A(x, y) and B(x, y), we let 𝐴(𝑥, 𝑦) = 80 exp[−𝑡(𝑥2 + 𝑦2 )] + 80, 𝐵(𝑥, 𝑦) = 80 exp[−𝑡(𝑥2 + 𝑦2 )]. 178
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 9. Experimental interferogram; (b) the reference phase achieved by the temporal Fourier transform algorithm; the unwrapped distributions achieved by three different algorithms: (c) APCA; (d)MSSM; (e)ASSF; (f) (g) (h) represents the phase errors obtained by the 3 algorithms compared with the reference phase: (f) APCA; (g) MSSM; (h) ASSF; (i) the cross-section curves of the 300th row in (f) (g) (h), respectively.
3. Experiments
wavelength of 632.8 nm is utilized as the light source, and a CCD with pixels size of 10𝜇m × 10𝜇m is employed to capture phase-shifting interferograms. 300-frame phase-shifting interferograms with total phase shifts of about 24 rad are captured. The PZT we used is WDT-2D20 made by CETC company, its voltage-phase shifts relationship is not linear, so
Next, we will utilize the proposed ASSF algorithm to perform phase retrieval from the experimental phase-shifting interferograms. A MachZehnder type PSI system is built, in which a He-Ne stabilized laser with 179
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
References
a self-calibration phase retrieval is needed. First, these 300-frame interferograms are utilized to achieve the reference spatial phase distribution with temporal Fourier transform algorithm [21], in which both the self-calibration and noise-resistance ability can be guaranteed by a large number of interferograms with a quasi-uniform phase shifts. Two significant properties of PSI algorithm, i.e. the self-calibration accuracy and noise-resistance ability are respectively verified. Generally speaking, for almost all the self-calibration PSI algorithms, the self-calibration ability can be guaranteed by the large fringes. When the fringes density decreases, especially when the fringes number in interferogram is less than one, the self-calibration accuracy will decrease. As shown in Fig. 8, in the first-group of experiment, 3-frame phase-shifting interferograms, containing few (less than one) fringe are chosen to achieve self-calibration phase retrieval by the APCA, MSSM and ASSF algorithms, respectively. It is found that when the fringes number is less than one, the APCA algorithm becomes invalid while the MSSM and ASSF algorithms still can work. However, we can still see an apparent circular distortion in Fig. 8g corresponds to the MSSM algorithm when it is compared with Fig. 8h corresponds to the ASSF algorithm. This is because the spatial spectrum of B2 (x, y)cos 2𝜙(x, y), B2 (x, y)sin 2𝜙(x, y) is more expanded in Fourier domain relative to B(x, y)cos 𝜙(x, y), B(x, y)sin 𝜙(x, y), which makes the self-calibration ability of ASSF algorithm be stronger relative to the MSSM algorithm, as discussed in above theoretical analysis. Subsequently, we will verify the noise-resistance ability of the ASSF algorithm. When the phase shifts are small, the error induced by the noise amplitude will be amplified. Thus, we choose 3-frame phase-shifting interferograms, i.e. the first, third, and fifth frame from 300-frame interferograms, approximately corresponding to the phase shifts of 0 rad, 0.16 rad and 0.32 rad, respectively. To reduce the unnecessary fringe density, the configuration of reference arm and objective arm is set as a symmetric. An objective with NA of 0.5 and magnification of 20 × are placed in the reference arm and objective arm. A light guide panel with the circular cavities array is chosen as the measured sample. Fig. 9c–e respectively show the phase distribution achieved by APCA, MSSM and ASSF algorithm after TIE unwrapping [40]. Fig. 9f–h show the difference between the wrapped phase distribution achieved by these 3 different algorithms with the reference phase. By comparing Fig. 9f and g with Fig. 9h, the obvious noise-resistance advantage of ASSF algorithm is observed. In the unwrapped spatial phase distribution, the cloud-like error is well-suppressed in Fig. 9e but be apparent in Fig. 9c and d. From Fig. 9i, it is also found that using the ASSF algorithm, the spatial phase fluctuation induced by the noise has been greatly suppressed.
[1] Kim KMyung. Digital holographic microscopy: principles, techniques, and applications. Springer; 2011. [2] Malacara D, Servín M, Malacara Z. Interferogram analysis for optical testing. Second Edn. CRC; 2005. [3] Bruning JH, Herriott DR, Gallagher JE, Rosenfeld DP, White AD, Brangaccio DJ. Digital wavefront measuring interferometer for testing optical surfaces and lenses. Appl Opt 1974;13:2693–703. [4] Bishara W, Zhu H, Ozcan A. Holographic opto-fluidic microscopy. Opt Express 2010;18:27499–510. [5] Choi W, Fang-Yen C, Badizadegan K, Oh S, Lue N, Dasari RR, Feld MS. Tomographic phase microscopy. Nat Methods 2007;4:717. [6] Kim T, Zhou R, Mir M, Babacan SD, Carney PS, Goddard LL, Popescu G. White-light diffraction tomography of unlabelled live cells. Nat Photonics 2014;8:256. [7] Marquet P, Rappaz B, Magistretti PJ, Cuche E, Emery Y, Colomb T, Depeursinge C. Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy. Opt Lett 2005;30:468–70. [8] Wang Z, Millet L, Mir M, Ding H, Unarunotai S, Rogers J, Gillette MU, Popescu G. Spatial light interference microscopy (SLIM). Opt Express 2011;19:1016–26. [9] Paganin D, Nugent KA. Noninterferometric phase imaging with partially coherent light. Phys Rev Lett 1998;80:2586–9. [10] Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72:156–60. [11] Arnison MR, Larkin KG, Sheppard CJR, Smith NI, Cogswell CJ. Linear phase imaging using differential interference contrast microscopy. J Microsc 2004;214:7–12. [12] Zhang T, Yamaguchi I. Three-dimensional microscopy with phase-shifting digital holography. Opt Lett 1998;23:1221–3. [13] Tian L, Waller L. Quantitative differential phase contrast imaging in an LED array microscope. Opt Express 2015;23:11394–403. [14] Zheng G, Horstmeyer R, Yang C. Wide-field, high-resolution Fourier ptychographic microscopy. Nat Photonics 2013;7:739. [15] Debnath SK, Park Y. Real-time quantitative phase imaging with a spatial phase-shifting algorithm. Opt Lett 2011;36:4677–9. [16] Wang Y, Qiu X, Xiong J, Li B, Zhong L, Liu S, Zhou Y, Tian J, Lu X. General spatial phase-shifting interferometry by optimizing the signal retrieving function. Opt Express 2017;25:7170–80. [17] Stern A, Javidi B. Space-bandwidth conditions for efficient phase-shifting digital holographic microscopy. J Opt Soc Am A 2008;25:736–41. [18] Servín M, Quiroga JA, Padilla JM. Fringe pattern analysis for optical metrology: theory, algorithms, and applications. Wiley-VCH; 2014. [19] Morgan CJ. Least-squares estimation in phase-measurement interferometry. Opt Lett 1982;7:368–70. [20] Larkin KG, Oreb BF. Design and assessment of symmetrical phase-shifting algorithms. J Opt Soc Am A 1992;9:1740–8. [21] Hao Q, Zhu Q, Hu Y. Random phase-shifting interferometry without accurately controlling or calibrating the phase shifts. Opt Lett 2009;34:1288–90. [22] Vargas J, Sorzano COS. Quadrature component analysis for interferometry. Opt Laser Eng 2013;51:637–41. [23] Gao P, Yao B, Lindlein N, Mantel K, Harder I, Geist E. Phase-shift extraction for generalized phase-shifting interferometry. Opt Lett 2009;34:3553–5. [24] Deng J, Wang H, Zhang D, Zhong L, Fan J, Lu X. Phase shift extraction algorithm based on Euclidean matrix norm. Opt Lett 2013;38:1506–8. [25] Liu F, Wu Y, Wu F. Correction of phase extraction error in phase-shifting interferometry based on Lissajous figure and ellipse fitting technology. Opt Express 2015;23:10794–807. [26] Deng J, Wang K, Wu D, Lv X, Li C, Hao J, Qin J, Chen W. Advanced principal component analysis method for phase reconstruction. Opt Express 2015;23:12222–31. [27] Wang Z, Han B. Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms. Opt Lett 2004;29:1671–3. [28] Guo H. Blind self-calibrating algorithm for phase-shifting interferometry by use of cross-bispectrum. Opt Express 2011;19:7807–15. [29] Farrell CT, Player MA. Phase step measurement and variable step algorithms in phase-shifting interferometry. Meas Sci Technol 1992;3:953. [30] Wang Y, Lu X, Liu Y, Tian J, Zhong L. Self-calibration phase-shifting algorithm with interferograms containing very few fringes based on Fourier domain estimation. Opt Express 2017;25:29971–82. [31] Wang Y, Li B, Zhong L, Tian J, Lu X. Spatial dual-orthogonal (SDO) phase-shifting algorithm by pre-recomposing the interference fringe. Opt Express 2017;25:17446–56. [32] Vargas J, Quiroga JA, Belenguer T. Phase-shifting interferometry based on principal component analysis. Opt Lett 2011;36:1326–8. [33] Vargas J, Quiroga JA, Belenguer T. Analysis of the principal component algorithm in phase-shifting interferometry. Opt Lett 2011;36:2215–17. [34] Weijuan Q, Yingjie Y, Choo CO, Asundi A. Digital holographic microscopy with physical phase compensation. Opt Lett 2009;34:1276–8. [35] Deng J, Wu D, Wang K, Vargas J. Precise phase retrieval under harsh conditions by constructing new connected interferograms. Sci Rep 2016;6:24416. [36] Yatabe K, Ishikawa K, Oikawa Y. Improving principal component analysis based phase extraction method for phase-shifting interferometry by integrating spatial information. Opt Express 2016;24:22881–91. [37] Yatabe K, Ishikawa K, Oikawa Y. Hyper ellipse fitting in subspace method for phase-shifting interferometry: practical implementation with automatic pixel selection. Opt Express 2017;25:29401–16. [38] Yatabe K, Ishikawa K, Oikawa Y. Simple, flexible, and accurate phase re-
4. Conclusions In this study, we focus on two significant challenges in PSI algorithm, i.e. the low noise phase retrieval from few interferograms with small phase shifts and accurate self-calibration phase retrieval from USFP. By analyzing the general properties of PSI method by introducing the combination mode theory, for the first challenge, we propose an advanced noise-resistance method for PSI algorithm. To overcome the second challenge, for the first time, we propose an self-calibration phase retrieval algorithm which can realize phase retrieval from fringe patterns contain less than one fringe without strict requirement of phase shifts. Based on the above procedures, the proposed ASSF algorithm can significantly improve the quality of phase retrieval with PSI method. The Matlab codes can be downloaded at Ref [41].
Acknowledgment This work is supported by National Natural Science Foundation of China grants (61475048, 61275015 and 61177005) 180
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
trieval method for generalized phase-shifting interferometry. J Opt Soc Am A 2017;34:87–96. [39] Gonzalez RC, Woods RE. Digital image processing. 3rd Edition. Prentice-Hall; 2007.
[40] Pandey N, Ghosh A, Khare K. Two-dimensional phase unwrapping using the transport of intensity equation. Appl Opt 2016;55:2418–25. [41] https://pan.baidu.com/s/1htwqUuK
181
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Advanced spatial spectrum fitting algorithm for significantly improving the noise resistance ability of self-calibration phase shifting interferometry Shengqi Cao, Yi Wang, Xiaoxu Lu, Liyun Zhong∗ Guangdong Provincial Key Laboratory of Nanophotonic Functional Materials and Devices, South China Normal University, Guangzhou 510006, China
a r t i c l e Keywords: Fringe analysis Interferometry Metrology Phase retrieval Phase measurement
i n f o
a b s t r a c t We propose a novel phase demodulation algorithm for phase shifting interferometry (PSI) named as advanced spatial spectrum fitting (ASSF) algorithm. By utilizing dual-operation of singular value decomposition (SVD) and hybrid temporally and spatially estimating the background intensity, the retrieved phase distribution will become more robust to the noise distribution compared with the traditional PSI. At the same time, with no particular temporal phase shifts distribution requirement, it realize stable self-calibration phase retrieval with few interferograms containing less than one fringe. Therefore, even when the interferograms are seriously corrupted by the random noise or the fringe number is small, the high quality phase retrieval is possible. Both simulation and experimental results will be used to demonstrate the outstanding performance of the proposed ASSF algorithm.
1. Introduction After the tested object being illuminated by a light beam, the spatial phase distribution encoded optical wavefront will change corresponding to the micro-nano structure of the object [1]. Various important quantitative applications of it have been presented, including profile testing [2,3], fluid monitoring [4], refractive index tomography [5,6] and biological imaging [7,8]. Although the direct detection of phase distribution is still impossible due to lack of devices with ultra-high temporal resolution, by encoding the phase information into optical intensity, a time-average quantity, we can achieve it indirectly. To date, some spatial phase encoding and demodulation methods have been reported [9– 14]. Among them the optical interferometry is regarded as an important candidate, in which the phase information is directly encoded into the interference fringe patterns in a linear way, and then the phase demodulation procedure becomes mathematically convenient and deterministic. In current optical interferometry, the off-axis configurations [15,16] are usually utilized to perform dynamic analysis and the phase-shifting configurations [17] are used to make the full utilization of space-bandwidth product of the camera. In phase shifting interferometry (PSI), we need to ensure the accurate estimation of phase shifts for accurate phase retrieval, otherwise, the retrieved phase distribution will suffer from a serious rippled error along the fringe orientations [18] . Initially, people have simplified the computational model by assuming that the phase shifts are known [19]. or the phase-shifter produce equal phase steps which are uniformly distributed in the range of [0, 2𝜋] [20]. Subsequently, it is found that this assump-
∗
tion is not valid in most cases. Under this model, although by increasing the number of captured interferograms we can release the rippled-like error [21,22], more time will be needed, and the vibrational robustness of the interferometric system will be correspondingly reduced. To solve this problem, the self-calibration phase-shifting algorithms are proposed and developed, in which the blind phase shifts are available [23–33] . Some of them are designed to work with 3 or less phase-shifting interferograms to reduce the necessary image capturing time, then some reasonable hypotheses of these interferograms are necessary to balance the equations and unknown parameters in the mathematical model. Generally speaking, in the condition that the interference-induced fringes number and the number of phase-shifting interferograms are both few, or the illumination field is nonuniform, it is very difficult to achieve accurate self-calibration phase retrieval because these hypotheses are no longer valid. This is a significant problem in PSI because the number reduction of interference fringes is usually appeared when we use white light to suppress the coherent artifact [8], and the reduction of the fringes number induced by the aberration is indeed helpful for more sufficient utilization of the space-bandwidth of the camera [34]. Another challenge in PSI is that the phase shifts between interferograms cannot be too small, otherwise the retrieved phase distribution will be greatly sensitive to the noise, which will lead to a restriction during practical measurement. Recently, a method is reported to realize self-calibration phase retrieval from only 3-frame phase-shifting interferograms containing less than one fringe [30]. In this method, the mid-frequency spatial spectrum of the fringe patterns is utilized to determine the phase shifts, so the large amplitude of mid-frequency spatial spectrum of the interfero-
Corresponding author. E-mail addresses: [email protected] (Y. Wang), [email protected] (L. Zhong).
https://doi.org/10.1016/j.optlaseng.2018.09.007 Received 1 June 2018; Received in revised form 8 August 2018; Accepted 10 September 2018 Available online 27 September 2018 0143-8166/© 2018 Elsevier Ltd. All rights reserved.
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
number of phase-sifting interferograms is more than three (N > 3). For 𝑁 ∑ a group of coefficients𝑘1 , 𝑘2 ......𝑘𝑁 , we can define 𝑘𝑛 𝐼𝑛 as a combina-
grams is preferred. However, in the ultra-sparse fringe patterns (USFP) contain less than one fringe, the amplitude of the mid-frequency spatial spectrum is weaker than the fringe patterns which contain many fringes. This lead to the accuracy decrease when we intend to use the USFPs to make phase retrieval. What’s more, in this algorithm the number of interferograms is fixed at 3, and we will clarify in the next section, to guarantee the accuracy, the phase shifts between arbitrary 2 interferograms cannot close to𝜋. Apparently, these restrictions may bring inconvenience in its practical applications. We have also noticed that although this method is not sensitive to the difference of fringe patterns, when the phase shifts are small, the noise will still be amplified. It has been reported that pre-removal of background intensity by spatial filtering instead of temporal subtraction is useful [35] for improving the signal-to-noise ratio (SNR) of phase retrieval in the case that phase shifts are small, but it is also found that the high-pass filtering function will lead to a serious signal distortion [31] due to the overlap of interference and background intensity. We will demonstrate later that the direct high-pass spatial filtering of these interferograms can both reduce the noise and signal amplitude, which indicates that the improvement of noise-resistance ability in Ref [35] is sensitive to the filtering window. In this manuscript, we proposed a phase-shifting algorithm named as advanced spatial spectrum fitting (ASSF) algorithm that can realize self-calibration phase retrieval from few phase shifting interferograms containing less than one fringe with no strict temporal phase shifts requirement and further boost the noise resistance ability in PSI. In the first step, we summarize the PSI algorithm by introducing a combination modes (CM) theory. According to this theory, we optimize the SNR of retrieved signals by SVD. Next, by using an advanced background intensity removal procedure, we achieve the interference signals with further improved SNR. Finally, by replying the SVD and performing parameters estimation using the spatial spectrum characters of the interferograms, we realize phase retrieval with strong self-calibration ability and ultrahigh SNR. Next, we will discuss the proposed algorithm in detail.
𝑛=1
tion mode (CM) for the N-frame interferograms. Apparently, a group of OCC contains three CMs, and in the right hand of Eqs. 2b and 2c, the first and second term represents the interference signal and the noise, respectively. Apparently, the SNR of CMs containing in different OCCs may be different as well. Our first goal is to search for the CMs that have maximum SNR. If the phase shifts are known, the CMs with the maximum SNR can be easily achieved by temporal least square fitting algorithm. However, the phase shifts cannot be determined precisely before calculation in practical. In this situation, we need a method to realize the blind optimization of SNR. Recently, researchers tend to extract the sub-vectors of the correlation matrix formed by the intensities of interferograms by using singular value decomposition (SVD) operation to achieve high SNR phase retrieval [26,32,36-38]. As we will discuss later, out of the perspective of the combination modes theory, these methods cannot always help us to achieve the optimum SNR. This is because that without optimization, the temporal average subtraction for background intensity A(x, y) removal operation embedded in these methods are blind. The optimum SNR can be achieved by this procedure only when the phase shifts are distributed in [0, 2𝜋] in a quasi-linear way. Only in this case, the temporal average of N-frame interferograms is accidentally corresponding to the first CM in the optimum OCC, i.e., 𝑎1 𝐼1 + .....𝑎𝑁 𝐼𝑁 . Before beginning the description of the proposed method, we briefly review the basic properties of the SVD operation. The SVD operation transforms a real and symmetric correlation matrix A constructed from N input data vectors to a diagonal matrix, and then output N eigenvectors and their eigenvalues. If N input data vectors are I1 , I2 ......IN , then the elements of matrix A can be described as ∑ 𝐴𝑚𝑛 = 𝐼𝑚 (𝑥, 𝑦) ⋅ 𝐼𝑛 (𝑥, 𝑦) 𝑚, 𝑛 = 1......𝑁. (3) 𝑥,𝑦
2. Methods
Two important properties of SVD operation used in this work are listed as follow
2.1. Combination modes of PSI and its optimization
a) The matrix A has N eigenvectors with length of N if the number of interferograms is N, i.e., 𝐶 𝐶𝑚 = (𝑘1𝑚 , 𝑘2𝑚 ......𝑘𝑁𝑚 ), 𝑚 = 1......𝑁. For these eigenvectors we have
First, we introduce a general theory to describe PSI algorithms named as combination mode theory. In PSI, the intensity of the nth phase-shifting interferogram can be described as 𝐼𝑛 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝜃𝑛 ] + 𝑜𝑛 (𝑥, 𝑦),
𝑛 = 1......𝑁.
(𝑘1𝑚 )2 + ...... + (𝑘𝑁𝑚 )2 = 1
(1)
where x and y denotes the transverse and longitude coordinates, respectively; A(x, y) and B(x, y) represents the background intensity and modulation amplitude, respectively; 𝜙(x, y)is the encoded phase distribution; 𝜃 n represents the phase shifts and 𝑛 = 1......𝑁are the index numbers of phase-shifting interferograms; on (x, y)denotes the noise distribution in the nth interferogram. Noticed that the goal of PSI algorithms are in fact to determine a set of orthogonal combination coefficients (OCC) for each interferograms to temporally combine them into B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y)as 𝑁 ∑ 𝑛=1 𝑁 ∑ 𝑛=1 𝑁 ∑ 𝑛=1
𝑎𝑛 𝐼𝑛 = 𝐴(𝑥, 𝑦) +
𝑁 ∑ 𝑛=1
𝑎𝑛 𝑜𝑛 .
𝑏𝑛 𝐼𝑛 = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) +
𝑐𝑛 𝐼𝑛 = 𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦) +
𝑛 𝑘𝑚 𝑘𝑛 + 𝑘𝑚 𝑘𝑛 + ...... + 𝑘𝑚 𝑁 𝑘𝑁 = 0 1 1 2 2
𝑏𝑛 𝑜𝑛 .
(2b)
𝑐𝑛 𝑜𝑛 .
(2c)
𝑛=1 𝑁 ∑ 𝑛=1
𝑚, 𝑛 = 1......𝑁.
(4a)
(4b)
(𝑘1𝑚 𝐼1 + 𝑘2𝑚 𝐼2 + ...... + 𝑘𝑁𝑚 𝐼𝑁 )(𝑘1𝑛 𝐼1 + 𝑘2𝑛 𝐼2 + ...... + 𝑘𝑁𝑛 𝐼𝑁 ) = 0 𝑚, 𝑛 = 1......𝑁.
(4c)
b) If there is an arbitrary vector with length of N, i.e., 𝐶 𝐶 ′ = (𝑐1 ′, 𝑐2 ′......𝑐𝑁 ′), and satisfy the condition of 𝑐 ′21 + 𝑐 ′22 + ...... + 𝑐′2𝑁 = 1, the function 𝑇 = (𝑐 ′1 𝐼1 + 𝑐 ′2 𝐼2 + ...... + 𝑐 ′𝑁 𝐼𝑁 )2 will reach its first maximum 𝜆1 in the case that 𝐶 𝐶 ′ = 𝐶 𝐶1 , then the second maximum 𝜆2 in the case that 𝐶 𝐶 ′ = 𝐶 𝐶2 and so on (𝜆1 > 𝜆2 > ...... > 𝜆N ).
(2a) 𝑁 ∑
𝑚 = 1......𝑁.
Each interferogram will contain three linear and independent components, i.e. A(x, y) B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y). If the spatial noise distributions 𝑜𝑛 (𝑥, 𝑦), 𝑛 = 1......𝑁of every interferograms have the same variance 𝜎 2 and are random and independent with each other, then according to the statistical theory, the variance of noise fluctuations in Eqs .(2b) and (2c) are (𝑏1 2 + 𝑏2 2 + ......𝑏𝑁 2 )𝜎 2 and (𝑐1 2 + 𝑐2 2 + ......𝑐𝑁 2 )𝜎 2 respectively. If we define the SNR of two interference signals i.e., B(x, y)cos 𝜙(x, y) and B(x, y)sin 𝜙(x, y), obtained by Eq. (2) as ∑ 2 𝐵 (𝑥, 𝑦)cos2 [𝜙(𝑥, 𝑦) + 𝜃𝑛 ]
Here, we definea1 ......aN , b1 ......bN and c1 ......cN as a group of OCC. If the OCC is known, we can achieve 𝜙(x, y)by performing the arctangent operation to the division ofB(x, y)sin 𝜙(x, y)andB(x, y)cos 𝜙(x, y). The OCC are directly related with the phase shifts, but it is not unique if the
𝑆𝑁 𝑅1 = 171
𝑥,𝑦
(𝑏1 2 + 𝑏2 2 + ......𝑏𝑁 2 )𝜎 2
.
(5a)
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
∑ 𝑆𝑁 𝑅2 =
𝑥,𝑦
As discussed in Section 2.1, to achieve accurate interference signals, it is needed to perform the background removal in PSI. Traditionally, this procedure can be realized by the temporal subtraction operation, but both the effective phase shifts and modulation amplitude may be reduced after this, and thus the noise level is amplified. The background removal can also be realized using spatial convolution, and this idea has directly given birth to the 2-step PSI. Recently, it is found that in the case that phase shifts are small, the noise amplification induced by the temporal subtraction background removal can be released by using spatial convolution instead [35]. To clarify its principle briefly, we let N = 3 and 𝜃1 = 0, 𝜃2 = 𝜋6 , 𝜃3 = 𝜋3 in Eq. (1), then we have
𝐵 2 (𝑥, 𝑦)sin2 [𝜙(𝑥, 𝑦) + 𝜃𝑛 ]
(𝑐1 2 + 𝑐2 2 + ......𝑐𝑁 2 )𝜎 2
.
(5a)
Then apparently, according to the properties of SVD above mentioned, 3 CMs with the maximum SNR, in which A(x, y), B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y)are contained in them should correspond to the first 3 eigenvectors or their orthogonal transformations. In the other CMs, we can easily understand that they contain the redundant components that can be described as 𝑁 ∑ 𝑛=1
𝑐𝑛 𝐼𝑛 = 0 +
𝑁 ∑ 𝑛=1
𝑐𝑛 𝑜𝑛 .
(6)
𝐼1 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝑜1 (𝑥, 𝑦).
In which only the noise terms are left. The SVD operation is in fact to help us excluding these redundant CMs. After the 3 optimized CMs with the same noise level (guaran𝑁 ∑ teed by the property a) of SVD operation), i.e., 𝑎1 𝐼1 + ...... + 𝑎𝑁 𝐼𝑁 , 𝑁 ∑ 𝑛=1
𝑏1 𝐼1 + ...... + 𝑏𝑁 𝐼𝑁 ,
𝑁 ∑ 𝑛=1
𝑛=1
𝑎𝑚
𝑛=1
𝑎𝑛
−
𝑐𝑚
𝑁 ∑
𝑛=1
𝜋 ] + 𝑜3 (𝑥, 𝑦). 3
(8c)
𝜋 + ] + 𝑜1 (𝑥, 𝑦) − 𝑜3 (𝑥, 𝑦). 6
(7b)
𝑎𝑚
𝑁 ∑
𝑎𝑛
−
𝑏𝑚
𝑁 ∑
𝑛=1
𝜋 𝐵(𝑥, 𝑦) sin[𝜙(𝑥, 𝑦) 12
𝜋 ] + 𝑜1 (𝑥, 𝑦) − 𝑜2 (𝑥, 𝑦). 12
𝐼2′ (𝑥, 𝑦) = 𝐼3 (𝑥, 𝑦) − 𝐼1 (𝑥, 𝑦) = 2 sin
𝑛=1
𝑁 ∑
𝐼3 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
+
𝑁 ∑
𝑚=1
𝑚=1
(8b)
(7a)
√ √𝑁 √∑ according to the Eq. (2). Here, two denominators √ ( √ √ √𝑁 √∑ and √ ( √
𝜋 ] + 𝑜2 (𝑥, 𝑦). 6
𝐼1′ (𝑥, 𝑦) = 𝐼1 (𝑥, 𝑦) − 𝐼2 (𝑥, 𝑦) = 2 sin
√ √ 2 √ ⎛ ⎞ √ 𝑎𝑛 𝐼𝑛 𝑐𝑛 𝐼𝑛 √ 𝑁 ⎜ ⎟ √ 𝑐 𝑛=1 𝑛=1 √ ∑ ⎜ 𝑎𝑚 𝑃2 = ( − )∕√ − 𝑚 ⎟ . √ 𝑁 𝑁 𝑁 𝑁 ⎜ ⎟ ∑ ∑ ∑ ∑ √𝑚=1 ⎜ 𝑎𝑛 𝑐𝑛 𝑎 𝑐 ⎟ ⎝ 𝑛=1 𝑛 𝑛=1 𝑛 ⎠ 𝑛=1 𝑛=1 𝑁 ∑
𝐼2 (𝑥, 𝑦) = 𝐴(𝑥, 𝑦) + 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
In which the variance of o1 (x, y), o2 (x, y), o3 (x, y) is equal to 𝜎 2 . The phase shifts of Eq. (8) are apparently small. If we remove A(x, y) by subtracting I1 (x, y) from I2 (x, y) and I3 (x, y), we have
𝑐1 𝐼1 + ...... + 𝑐𝑁 𝐼𝑁 are obtained, we can sub-
tract the A(x, y) components from them by √ √ 2 𝑁 𝑁 ∑ ⎛∑ ⎞ √ ⎛ ⎞ √ 𝑎 𝐼 𝑏 𝐼 𝑁 𝑛 𝑛 𝑛 𝑛⎟ √ ⎜ ⎜ ⎟ √ ∑ 𝑎𝑚 𝑏 𝑛=1 𝑛=1 ⎟∕ √ ⎜ 𝑃1 = ⎜ − − 𝑚 ⎟ . √ 𝑁 𝑁 𝑁 𝑁 ⎜ ∑ ⎟ √ ⎜ ⎟ ∑ ∑ ∑ √𝑚=1 ⎜ ⎜ 𝑎 𝑏𝑛 ⎟ 𝑎 𝑏 ⎟ ⎝ 𝑛=1 𝑛 ⎠ ⎝ 𝑛=1 𝑛 𝑛=1 𝑛 ⎠ 𝑛=1
(8a)
(9a)
𝜋 𝐵(𝑥, 𝑦) sin[𝜙(𝑥, 𝑦) 6 (9b)
The variance of noise terms are increased from 𝜎 2 in Eq. (8) to 2𝜎 2 in Eq. (9) according to the statistical theory. the modulation amplitude 𝜋 is reduced from B(x, y) to 2 sin 12 𝐵(𝑥, 𝑦) and 2 sin 𝜋6 𝐵(𝑥, 𝑦); the effective 𝜋 𝜋 𝜋 phase shifts are decreased from 6 and 𝜋3 to 𝜋6 − 12 = 12 . However, if A(x, y) is instead removed by the spatial filtering and there is no spatial distortion of the interference terms by this filtering function, which is an ideal condition, then we have
2
)
𝑏𝑛
2
) are employed to normalize the noise variance
𝐼1′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝑜1 (𝑥, 𝑦).
𝑐𝑛
of Eqs. (7a) and (7b). Then we can use self-calibration algorithms [25, 26, 29] to obtain the two interference terms, i.e., B(x, y)sin 𝜙(x, y)and B(x, y)cos 𝜙(x, y) from P1 and P2 and the spatial phase distribution can be obtained by the antitangent of their division. In comparison, due to the disadvantage of background removal procedure, the sub-space optimization method used before are not in the best way because the blind background removal by temporal subtraction may include the redundant CMs, for example, the temporal average subtraction in Ref [32] and the subtraction of the first interferograms in Ref [26] . To clarify this, 30-frame phase shifting interferograms are captured and one of them is shown in Fig. 1a. Next, we use these interferograms and advanced iterative algorithm (AIA) method to demodulate the encoded phase and regard it as the ground truth., Then, we choose first 5 phase shifting interferograms to make phase retrieval. Fig. 1c and d respectively present the corresponding phase errors by using background removal procedure adopted by APCA method [26] in which the first interferograms is subtracted by another, and the proposed CMs optimization operation. The comparison between Fig. 1c – e have clarified the importance of CM optimization.
(10a)
𝐼2′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
𝜋 ] + 𝑜2 (𝑥, 𝑦). 6
(10b)
𝐼3′′ (𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) +
𝜋 ] + 𝑜3 (𝑥, 𝑦). 3
(10c)
Compared with Eq. (9), the solution of Eq. (10) is apparently with higher SNR because the larger modulation amplitude, smaller noise variance and increase of equations number. However, different with the ideal condition, there is a spatial phase reconstruction effect if the spatial filtering operation is directly applied to the interferograms. In following analysis, for simplicity, we will neglect the y dependence. In general, if a spatial filtering operation is utilized, corresponding to a convolution operation (H(x)⊗) in spatial domain or a multiply operation (H(fx ) · ) in frequency domain, the distortion of the interference term can be expressed as [31] ∑ 𝐼̃𝑛 (𝑥) = 𝐻(𝑥) ⊗ 𝐵(𝑥) cos[𝜙(𝑥) + 𝛿𝑛 ] = 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos[𝜙(𝑥 − 𝜏) + 𝛿𝑛 ] 𝜏 ∑ = 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos[𝜙(𝑥) + 𝑓 (𝑥, 𝜏) + 𝛿𝑛 ] 𝜏 ∑ = cos[𝜙(𝑥) + 𝛿𝑛 ] 𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏) − sin[𝜙(𝑥) 𝜏 ∑ +𝛿𝑛 ] 𝐻(𝜏)𝐵(𝑥 − 𝜏) sin 𝑓 (𝑥, 𝜏). √ 𝜏 √ = 𝑇 (𝑥) cos[𝜙(𝑥) + 𝛿𝑛 + 𝜂(𝑥)] = 𝑇 (𝑥) cos[𝜙′ (𝑥) + 𝛿𝑛 ]
2.2. Advanced noise suppression strategy Out of the traditional perspective, after three optimized CMs are retrieved, we have already achieved the signals with the highest SNR. However, when the phase shifts are small, the error induced by the noise may still too large even the optimized CMs are used. In these cases, we need a more efficient noise suppression strategy.
(11) in which 𝑓 (𝑥, 𝜏) = 𝜙(𝑥 − 𝜏) − 𝜙(𝑥) 172
(12a)
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 1. (a) Interferogram of a single macrophage cell; (b) The retrieved phase with aberration phase removed using 30-frame phase shifting interferograms and AIA algorithm which is regarded as the ground truth; The errors obtained with (c) the reported APCA background removal procedure; (d) the proposed CMs optimization operation. (e) The 200th row of (c) and (d).
𝑇 (𝑥) = [
∑
𝐻(𝜏)𝐵(𝑥 − 𝜏) sin 𝑓 (𝑥, 𝜏)]2 + [
∑
𝜏
be utilized only in the harsh condition with very small phase shifts. When the phase shifts become large or the noise level is moderate, the error induced by this distortion will become larger relative to the error induced by noise and the reduction of modulation amplitude may lead to an overall decrease of SNR instead of increase. Next, we will demonstrate an advanced strategy to remove A(x, y) without distorting the phase distribution or amplifying the noise level. Due to there are 3 components in the interferograms, i.e. A(x, y), B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y), there are 3 optimum CMs corresponding to the first 3 eigenvectors of the correlation matrix as follow
𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏)]2
𝜏
(12b) ∑ 𝜂(𝑥) = arccos(
𝜏
𝐻(𝜏)𝐵(𝑥 − 𝜏) cos 𝑓 (𝑥, 𝜏) √
).
𝑇 (𝑥)
(12c)
𝜙′ (𝑥) = 𝜙(𝑥) + 𝜂(𝑥).
(12d)
From Eq. (12d), we can observe that the phase distribution is changed from the original 𝜙(x)to 𝜙′(x)with an addictive distortion term𝜂(x), and this distortion will bring in a serious error especially in the case that the fringes density is relative small. It is worth noting that although the interference terms have more extended spatial spectrum relative to the background term, they are both centered at the zero frequency and thus overlapped. That is to say, when high-pass filtering operation is applied to filter out the background, it will also reduce the ∑√ ∑ amplitude of interference term, i.e. 𝑇 (𝑥) < 𝐵(𝑥), and contribute 𝑥
𝑇1 (𝑥, 𝑦) = 𝐶 𝑀1 = 𝑘11 𝐼1 + 𝑘12 𝐼2 + ......𝑘1𝑁 𝐼𝑁 𝑁 ∑ = 𝑡1 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿1 ] + 𝑘1𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13a)
𝑇2 (𝑥, 𝑦) = 𝐶 𝑀2 = 𝑘21 𝐼1 + 𝑘22 𝐼2 + ......𝑘2𝑁 𝐼𝑁 𝑁 ∑ = 𝑡2 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿2 ] + 𝑘2𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13b)
𝑇3 (𝑥, 𝑦) = 𝐶 𝑀3 = 𝑘31 𝐼1 + 𝑘32 𝐼2 + ......𝑘3𝑁 𝐼𝑁 𝑁 ∑ = 𝑡3 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿3 ] + 𝑘3𝑛 [𝐴(𝑥, 𝑦) + 𝑜𝑛 (𝑥, 𝑦)]
(13c)
𝑛=1
𝑛=1
𝑥
to the reduction of SNR consequently. Therefore, the above method can
𝑛=1
173
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 2. Spatial spectrums of three CMs, in which the background and interference terms are mixed. (f2 − , f 1 − ) and (f 1 + , f 2 + ) denote the spatial frequency band where mainly occupied by the interference term. (f 1 − , f 1 + ) denotes the frequency band where mainly occupied by background intensity A(x,y).
The basic idea is that, if we want to filtered out A(x, y) without the interference terms distortion, we need to avoid the overlap of them. There is a beneficial property can be helpful, which is, in PSI, the phase contrast is mainly induced by the spatial phase difference between the reference and object wave. That is to say, the background A(x, y) and modulation amplitude B(x, y) changes spatially slower than B(x, y)cos 𝜙(x, y) and B(x, y)sin 𝜙(x, y). As shown in Fig. 2, we can separate the spatial spectrum of the fringe patterns into three parts: (f2 − , f 1 − ) and (f 1 + , f 2 + ) are mainly occupied by the interference terms and (f 1 − , f 1 + ) is mainly occupied by A(x,y). Next, we will use this spectrum dividing properties of fringe patterns to separate A(x,y) and the interference terms. It is easy to understand that there exists k1 and k2 that satisfy 𝑘1 𝑇1 (𝑥, 𝑦) + 𝑘2 𝑇2 (𝑥, 𝑦) + 𝑇3 (𝑥, 𝑦) = 𝐿0 (𝑥, 𝑦)
and Eq. (14) as 𝐴(𝑥, 𝑦) =
Here, 𝐿̃ 0 (𝑥, 𝑦)denotes the low-pass function of L0 (x, y). By subtracting A(x, y) described by Eq. (16) from every interferograms described in Eq. (1), we have that 𝐼̃𝑛 (𝑥, 𝑦) = 𝐼𝑛 (𝑥, 𝑦) − 𝐴(𝑥, 𝑦) = 𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝜃𝑛 ] + 𝑜𝑛 (𝑥, 𝑦).
(17)
Clearly, there are only two independent components in 𝐼̃𝑛 (𝑥, 𝑦)left now. Then, we can perform the procedure of SVD again, and the correlation matrix elements are changed from Eq. (3) to ∑ 𝐴̃ 𝑚𝑛 = 𝐼̃𝑚 (𝑥, 𝑦) ⋅ 𝐼̃𝑛 (𝑥, 𝑦) 𝑚, 𝑛 = 1......𝑁. (18)
(14)
𝑥,𝑦
In which the interference terms is compensated. Here, L0 (x, y) is a quasi-constant in spatial domain and be proportional to A(x, y). Now, we intend to determine the values of k1 and k2 . It has demonstrated [30] that by using a low-frequency suppression spatial filter, the accuracy of A(x, y) subtraction is improved relative to the direct spatial average subtraction. Therefore, we firstly apply the high-pass filtering operation to both sides of Eq. (14) and then L0 (x, y) can be well removed, and we have 𝑘1 𝑇̃1 (𝑥, 𝑦) + 𝑘2 𝑇̃2 (𝑥, 𝑦) + 𝑇̃3 (𝑥, 𝑦) = 0
𝐿̃ 0 (𝑥, 𝑦) 𝑘1 (𝑘11 + ...... + 𝑘1𝑁 ) + 𝑘2 (𝑘21 + ...... + 𝑘2𝑁 ) + (𝑘31 + ...... + 𝑘3𝑁 ) (16)
Thus, 2 instead 3 CMs with the highest SNR are achieved as 𝑇̃1 (𝑥, 𝑦) = 𝑘̃ 11 𝐼1 + 𝑘̃ 12 𝐼2 + ......𝑘̃ 1𝑁 𝐼𝑁 = 𝑡1 ′𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿1 ′]
(19a)
𝑇̃2 (𝑥, 𝑦) = 𝑘̃ 21 𝐼1 + 𝑘̃ 22 𝐼2 + ......𝑘̃ 2𝑁 𝐼𝑁 = 𝑡2 ′𝐵(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿2 ′]
(19b)
In which A(x, y) is removed and only the interference intensities are left. Till now, the effective phase shifts and modulation amplitude in Eq. (19) has been optimized based on the subtraction operation of A(x, y) through using the spatial property of interferogram in advance. The flow chart of the noise-suppressing strategy described in this section is given in Fig. 3 as follow Unfortunately, we cannot achieve the phase distribution due to 𝑡1 ′ , 𝑡2 ′ , 𝛿2 ′ − 𝛿1 ′ are still unknown. Therefore, it is needed to perform a self-calibration procedure for accurate spatial phase extraction. Next, we will demonstrate a rapid and accurate self-calibration algorithm using Eq. (19).
(15)
To solve the spatial least square problem (𝑘1 , 𝑘2 ) = ∑ 2 min{ [𝑘1 𝑇̃1 (𝑥, 𝑦) + 𝑘2 𝑇̃2 (𝑥, 𝑦) + 𝑇̃3 (𝑥, 𝑦)] }, we can determine k1 and k2 . 𝑥,𝑦
Then, we back substitute k1 and k2 to Eq. (14) to achieve L0 (x, y). The high-frequency spatial spectrum in the obtained L0 (x, y) in this step is corresponding to the noise fluctuations and the residual interference terms. Next, by utilizing a low-pass filtering operation to L0 (x, y) to obtain A(x, y) instead of direct high-pass filtering operation for each interferogram [35], we can avoid the distortion of interference terms. Another advantage of this procedure is that the filtering procedure is performed after the compensation of interference terms, so the final resolution will not be reduced, which is a significant disadvantage in traditional noise suppression method with low-pass filtering [39]. The relationship between A(x, y) and L0 (x, y) can be determined by Eq. (13)
2.3. Self-calibration phase retrieval from two optimum CMs We know that not only A(x, y) but also B(x, y) are quasi-constant in spatial domain. To sufficiently use both of these two properties, we will separate B(x, y) from the interference terms firstly. Neglecting the 174
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 3. The flow chart of the advanced noise suppressing strategy.
in which 𝐾̃ 1 , 𝐾̃ 2 and 𝐾̃ 3 denotes the high-pass filtering results of K1 , K2 and K3 , and M0 (x,y) is eliminated by the filtering operation. By solving the linear least-square problem, i.e. (𝑙1 , 𝑙2 ) = ∑ 2 min[ |𝑙1 𝐾̃ 1 (𝑥, 𝑦) + 𝑙2 𝐾̃ 2 (𝑥, 𝑦) + 𝐾̃ 3 (𝑥, 𝑦)| ], we can achieve l1 and l2 . Con-
piston phase shifts, we have 𝛿1 ′ = 0, 𝛿2 ′ = 𝛿, and then 𝐾1 (𝑥, 𝑦) = 𝑇̃1 ⋅ 𝑇̃2 = 𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) cos[𝜙(𝑥, 𝑦) + 𝛿] = =
𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥,𝑦) {cos[2𝜙(𝑥, 𝑦) + 𝛿] + cos 𝛿} 2 𝑡1 ′ 𝑡2 ′ 𝐵 2 (𝑥,𝑦) [cos 2𝜙(𝑥, 𝑦) cos 𝛿 − sin 2𝜙(𝑥, 𝑦) sin 𝛿 2
𝐾2 (𝑥, 𝑦) = 𝑇̃12 = 𝑡1 ′2 𝐵 2 (𝑥, 𝑦)cos2 𝜙(𝑥, 𝑦) =
(20a)
𝑥,𝑦
+ cos 𝛿].
sequently, the accuracy of the blind calibration is directly determined by the accuracy of the achieved l1 and l2 . Compared with the previous MSSM method [30], we can conclude the advantages of proposed algorithm as following: First, if the phase shifts is 𝜋, the MSSM algorithm is invalid, but in this method, this restriction can be canceled by adding a spatial quasi-constant assumption of B(x, y) compared with Ref [30]. Second, the interference terms B2 (x, y)cos 2𝜙(x, y)and B2 (x, y)sin 2𝜙(x, y) in Eq .(13) have stronger mid-frequency amplitude than the interference terms B(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y) used in MSSM algorithm, this will bring stronger blind calibration ability for this algorithm. Thus, in the case of USFP, it will reveal more obvious advantage. Finally, the MSSM algorithm is specific for the case of 3 phase shifting interferograms are captured and this method can make full use of the N-frame (N ≥ 3) interferograms.
𝑡1 ′2 𝐵 2 (𝑥, 𝑦) [cos 2𝜙(𝑥, 𝑦) + 1]. 2 (20b)
𝐾3 (𝑥, 𝑦) = 𝑇̃22 = 𝑡2 2 𝐵 2 (𝑥, 𝑦)cos2 [𝜙(𝑥, 𝑦) + 𝛿] = =
𝑡1 2 𝐵 2 (𝑥,𝑦) {cos[2𝜙(𝑥, 𝑦) + 2𝛿] + 1} 2 𝑡1 ′ 2 𝐵 2 (𝑥,𝑦) [cos 2𝜙(𝑥, 𝑦) cos 2𝛿 − sin 2𝜙(𝑥, 𝑦) sin 2𝛿 2
(20c) + 1].
There are three components B2 (x, y), B2 (x, y)cos 2𝜙(x, y)and B2 (x, y)sin 2𝜙(x, y) in K1 , K2 and K3 . In general, B2 (x, y) is a quasi-constant in spatial domain, then there will exist l1 and l2 that satisfy 𝑙1 𝐾1 (𝑥, 𝑦) + 𝑙2 𝐾2 (𝑥, 𝑦) + 𝐾3 (𝑥, 𝑦) = 𝑀0 (𝑥, 𝑦). in which M0 (x, y) is proportional to (𝑇̃2 +
B2 (x,
𝑙1 𝑙 2 𝑇̃1 )2 + (𝑙2 − 1 )𝑇̃12 = 𝑃0 (𝑥, 𝑦). 2 4
(21) 3. Simulations and discussions
y). Then we have
Following, we will present the performance of proposed ASSF algorithm through using two simulations. The size of simulated fringe pattern is 500 × 500 pixels with pixel size of 0.01 mm. 5-frame of fringe patterns with phase shifts of 0 rad, 0.2 rad, 0.3 rad, 0.4 rad, 0.5 rad are generated for each simulation; In the first simulation, the phase model that induced several fringes is set as 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ). In the second simulation, the phase model induced few (less than one) fringe is set as 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500); the background and modulation intensity is A(x,y) = 50exp[− 0.2(x2 + y2 )] + 50 and B(x,y) = 30exp[− 0.2(x2 + y2 )] + 30, respectively. The SNR of all these fringe patterns are fixed at 50 dB. Fig. 4 shows the results of phase retrieval achieved with the proposed ASSF, MSSM [30], APCA [26] and PCAN [35] algorithms. Note that the noise fluctuation in fringe patterns will lead to a serious error when the phase shifts are small. Even when the fringes number is large enough, the APCA and MSSM algorithms only can achieve noisy phase distribution while the ASSF algorithm can achieve accurate one, in which both the rippled error induced by the miss-calibration of phase shifts and the error induced by the random noise can be well-suppressed. Although the PCAN algorithm can achieve smooth profile of phase dis-
(22)
P0 (x, y) is also a quasi-constant in spatial domain. If the piston phase shifts is neglected, the first and second terms in Eq. (22) are corresponding toB(x, y)cos 𝜙(x, y)and B(x, y)sin 𝜙(x, y), respectively. Apparently, if l1 and l2 are determined, we have 𝑇̃2 +
𝑙1 𝑇̃ = 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦). 2 1
(23a)
√ 𝑙2 −
𝑙1 2 𝑇̃ = 𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦). 4 1
(23b)
And then 𝜙(x, y)can be determined by 𝜙(𝑥, 𝑦) = arctan[
𝐵(𝑥, 𝑦) sin 𝜙(𝑥, 𝑦) ]. 𝐵(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦)
(24)
To achieve l1 and l2 , we perform the spatial-high pass filtering operation to both two sides of Eq. (21), then 𝑙1 𝐾̃ 1 (𝑥, 𝑦) + 𝑙2 𝐾̃ 2 (𝑥, 𝑦) + 𝐾̃ 3 (𝑥, 𝑦) = 0.
(25) 175
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 4. (a1) (a2) the theoretical phase distribution of 𝜙1 (x, y) and one of the 5-frame simulated phase-shifting interference patterns, respectively; the spatial phase distributions of 𝜙1 (x, y) achieved with different algorithms: (a3) MSSM (a5)APCA (a7)ASSF (a9) PCAN; (a4) (a6) (a8) (a10) is the differences between 𝜙1 (x, y) and (a3) (a5) (a7) and (a9), respectively; (b1) (b2) the theoretical phase distribution of 𝜙2 (x, y)and one of 5-frame simulated phase-shifting interference patterns, respectively; the spatial phase distributions 𝜙2 (x, y)achieved by different algorithms: (b3) MSSM (b5) APCA (b7) ASSF (b9) PCAN; (b4) (b6) (b8) (b10) is the differences between 𝜙1 (x, y) and (b3) (b5) (b7) and (b9), respectively.
tribution which indicates less noisy result, its phase distortion becomes very serious, especially in the center area of interference pattern with sparse fringe, as shown in Fig. 4(a10). In the second simulation, the fringe number is decreased and the influence of noise fluctuation is still large due to the small phase shifts. We can see that only the proposed ASSF algorithm is available in this case. The small range of the phase distribution has violated the prior-assumptions of APCA algorithm and the sparsity of the fringe patterns makes the high-pass filtering of the interferograms be invalid in PCAN algorithm. The strong noise also has overwhelmed the achieved spatial phase distribution except with ASSF algorithm. These results prove that the proposed ASSF algorithm can
be utilized at the most difficult situations in PSI, in which the fringes number is small and the phase shifts are small at the same time. Note that in both the ASSF algorithm and the PCAN algorithm, the improvement of noise-resistance ability is dependent on the quality of background removal, which is directly related with the accuracy of phase retrieval. As shown in Fig. 5, the A(x, y) achieved by the proposed ASSF algorithm reveals obvious accuracy advantage compared with the PCAN algorithm. Like the ASSF algorithm [30], the proposed MSSM algorithm also can work well in the case that the fringes number is less than one. However, when the phase shifts between arbitrary two interferograms
176
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 5. (a) theoretical spatial distribution of A(x, y); A(x, y) achieved by the proposed ASSF algorithm in the case that (b) 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ); (c) 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500); (d) (e) the differences between (b) (c) and (a), respectively; A(x, y) achieved by the PCAN algorithm in the case that (f) 𝜙1 (𝑥, 𝑦) = 0.5𝜋(𝑥2 + 𝑦2 ); (g) 𝜙2 (𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500), respectively; (h) (i) the differences between (f) (g) and (a), respectively.
is equal to𝜋, the MSSM algorithm becomes invalid if the total number of interferograms N = 3. This is because in this case, the spatial quasiconstant assumption of A(x, y) can no longer support the determination of phase shifts mathematically in the MSSM algorithm. However, in the ASSF algorithm, by introducing an another assumption, i.e. the slowly spatial variation of B(x, y), we can get rid of the above restriction. Fig. 6 presents the RMSEs of phase retrieval achieved with the ASSF and MSSM algorithms when the phase shifts are changed from 2.9 rad to 3.3 rad. According to Fig. 6, although the ASSF and MSSM algorithm can both realize self-calibration phase retrieval from such a sparse fringe pattern, we can clearly observe a sharp peak near 𝜋 in the blue curve, which corresponds to the invalidity of the MSSM algorithm when the phase shift is close to 𝜋. The green curve changes slow with different phase
Fig. 6. The RMSE of phase retrieval using MSSM and ASSF when the phase shifts near of 𝜋.
Fig. 7. The accuracy of phase retrieval using different algorithm described by the corresponding RMSEs with different (a) phase shifts; (b) SNR of fringe patterns; (c) fringe density; (d) variation rate of A(x, y) and B(x, y). 177
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 8. Experimental interferogram encoded with a circular spatial phase distribution; (b) the reference phase achieved by the temporal Fourier transform algorithm; the spatial phase distributions achieved by different algorithms (c) APCA; (d) MSSM; (e) ASSF; (f) (g) (h) the differences between (b) and (c) (d) (e), respectively;.
The phase model is given by 𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500)+0.2(𝑥2 + 𝑦2 ), the phase shifts between the adjacent interferograms are 0.2 rad and the SNR of the interferograms are 50 dB, respectively. When t changes from 0.01 to 0.3, the RMSEs are plotted in Fig. 7d. The results of the PCAN algorithm are not shown in these curves because it cannot output a reasonable phase distribution under some values of these parameters due to the serious distortion of the interference terms by the filtering bank. Clearly, from Fig. 7, we can see several significant advantages of proposed ASSF algorithm: First, we can see that if k < 0.15 rad, which represents small phase shifts, the accuracy of phase retrieval with the ASSF algorithm is almost unchanged while the results with the APCA and MSSM algorithms are greatly reduced (Fig. 7a). Second, it has no strict requirement of the fringe density (Fig. 7c). Third, although the accuracy with the ASSF algorithm is decreased with the SNR of fringe pattern (Fig. 7b), but the noise-suppressing strategy can work in each noise level, so the corresponding result is more reliable compared with APCA and MSSM algorithms. Specially, when t value is large, which indicates that A(x, y) and B(x, y) are nonuniform, the accuracy of ASSF algorithm decreases slower compared with the MSSM algorithm due to the former has stronger mid-frequency spatial spectrum (Fig. 7d).
shifts, indicates that the ASSF algorithm has no such a restriction. Next, we will present the influence of the phase shifts, noise level, fringe density and uniformity of A(x, y), B(x, y) on the accuracy of phase retrieval with the ASSF algorithm. First, we let A(x, y) = 80exp[− 0.05(x2 + y2 )] and B(x, y) = 80exp[− 0.05(x2 + y2 )], 𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500), the SNR is fixed at 50 dB, the phase shift between the adjacent interferograms are set as k rad and 5-frame interferograms are generated. The RMSEs are plotted when k changes from 0.1 to 0.4, as Fig. 7a shown. Next, we let the phase shift between the adjacent interferograms be 0.2 rad, let the SNR of the fringe patterns changes from 35 dB to 65 dB and the RMSEs are plotted in Fig. 7b. Furthermore, we change the phase model to𝜙(𝑥, 𝑦) = 0.2 ⋅ 𝑝𝑒𝑎𝑘𝑠(500)+𝑚(𝑥2 + 𝑦2 ), in which an addictive quadratic phase component is added. In practical, this quadratic phase component is often induced when the optical system is asymmetric in the reference and objective arms, and consequently increase the density of the fringes. We change the value of m from 0 to 3 and keep the phase shifts between the adjacent interferograms be 0.2 rad and SNR of 50 dB. The RMSEs are plotted in Fig. 7c. At last, we know that in most cases the background and modulation intensity are both nonuniform, to investigate the accuracy robustness of the proposed ASSF algorithm under different spatial changing rate of the A(x, y) and B(x, y), we let 𝐴(𝑥, 𝑦) = 80 exp[−𝑡(𝑥2 + 𝑦2 )] + 80, 𝐵(𝑥, 𝑦) = 80 exp[−𝑡(𝑥2 + 𝑦2 )]. 178
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
Fig. 9. Experimental interferogram; (b) the reference phase achieved by the temporal Fourier transform algorithm; the unwrapped distributions achieved by three different algorithms: (c) APCA; (d)MSSM; (e)ASSF; (f) (g) (h) represents the phase errors obtained by the 3 algorithms compared with the reference phase: (f) APCA; (g) MSSM; (h) ASSF; (i) the cross-section curves of the 300th row in (f) (g) (h), respectively.
3. Experiments
wavelength of 632.8 nm is utilized as the light source, and a CCD with pixels size of 10𝜇m × 10𝜇m is employed to capture phase-shifting interferograms. 300-frame phase-shifting interferograms with total phase shifts of about 24 rad are captured. The PZT we used is WDT-2D20 made by CETC company, its voltage-phase shifts relationship is not linear, so
Next, we will utilize the proposed ASSF algorithm to perform phase retrieval from the experimental phase-shifting interferograms. A MachZehnder type PSI system is built, in which a He-Ne stabilized laser with 179
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
References
a self-calibration phase retrieval is needed. First, these 300-frame interferograms are utilized to achieve the reference spatial phase distribution with temporal Fourier transform algorithm [21], in which both the self-calibration and noise-resistance ability can be guaranteed by a large number of interferograms with a quasi-uniform phase shifts. Two significant properties of PSI algorithm, i.e. the self-calibration accuracy and noise-resistance ability are respectively verified. Generally speaking, for almost all the self-calibration PSI algorithms, the self-calibration ability can be guaranteed by the large fringes. When the fringes density decreases, especially when the fringes number in interferogram is less than one, the self-calibration accuracy will decrease. As shown in Fig. 8, in the first-group of experiment, 3-frame phase-shifting interferograms, containing few (less than one) fringe are chosen to achieve self-calibration phase retrieval by the APCA, MSSM and ASSF algorithms, respectively. It is found that when the fringes number is less than one, the APCA algorithm becomes invalid while the MSSM and ASSF algorithms still can work. However, we can still see an apparent circular distortion in Fig. 8g corresponds to the MSSM algorithm when it is compared with Fig. 8h corresponds to the ASSF algorithm. This is because the spatial spectrum of B2 (x, y)cos 2𝜙(x, y), B2 (x, y)sin 2𝜙(x, y) is more expanded in Fourier domain relative to B(x, y)cos 𝜙(x, y), B(x, y)sin 𝜙(x, y), which makes the self-calibration ability of ASSF algorithm be stronger relative to the MSSM algorithm, as discussed in above theoretical analysis. Subsequently, we will verify the noise-resistance ability of the ASSF algorithm. When the phase shifts are small, the error induced by the noise amplitude will be amplified. Thus, we choose 3-frame phase-shifting interferograms, i.e. the first, third, and fifth frame from 300-frame interferograms, approximately corresponding to the phase shifts of 0 rad, 0.16 rad and 0.32 rad, respectively. To reduce the unnecessary fringe density, the configuration of reference arm and objective arm is set as a symmetric. An objective with NA of 0.5 and magnification of 20 × are placed in the reference arm and objective arm. A light guide panel with the circular cavities array is chosen as the measured sample. Fig. 9c–e respectively show the phase distribution achieved by APCA, MSSM and ASSF algorithm after TIE unwrapping [40]. Fig. 9f–h show the difference between the wrapped phase distribution achieved by these 3 different algorithms with the reference phase. By comparing Fig. 9f and g with Fig. 9h, the obvious noise-resistance advantage of ASSF algorithm is observed. In the unwrapped spatial phase distribution, the cloud-like error is well-suppressed in Fig. 9e but be apparent in Fig. 9c and d. From Fig. 9i, it is also found that using the ASSF algorithm, the spatial phase fluctuation induced by the noise has been greatly suppressed.
[1] Kim KMyung. Digital holographic microscopy: principles, techniques, and applications. Springer; 2011. [2] Malacara D, Servín M, Malacara Z. Interferogram analysis for optical testing. Second Edn. CRC; 2005. [3] Bruning JH, Herriott DR, Gallagher JE, Rosenfeld DP, White AD, Brangaccio DJ. Digital wavefront measuring interferometer for testing optical surfaces and lenses. Appl Opt 1974;13:2693–703. [4] Bishara W, Zhu H, Ozcan A. Holographic opto-fluidic microscopy. Opt Express 2010;18:27499–510. [5] Choi W, Fang-Yen C, Badizadegan K, Oh S, Lue N, Dasari RR, Feld MS. Tomographic phase microscopy. Nat Methods 2007;4:717. [6] Kim T, Zhou R, Mir M, Babacan SD, Carney PS, Goddard LL, Popescu G. White-light diffraction tomography of unlabelled live cells. Nat Photonics 2014;8:256. [7] Marquet P, Rappaz B, Magistretti PJ, Cuche E, Emery Y, Colomb T, Depeursinge C. Digital holographic microscopy: a noninvasive contrast imaging technique allowing quantitative visualization of living cells with subwavelength axial accuracy. Opt Lett 2005;30:468–70. [8] Wang Z, Millet L, Mir M, Ding H, Unarunotai S, Rogers J, Gillette MU, Popescu G. Spatial light interference microscopy (SLIM). Opt Express 2011;19:1016–26. [9] Paganin D, Nugent KA. Noninterferometric phase imaging with partially coherent light. Phys Rev Lett 1998;80:2586–9. [10] Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am 1982;72:156–60. [11] Arnison MR, Larkin KG, Sheppard CJR, Smith NI, Cogswell CJ. Linear phase imaging using differential interference contrast microscopy. J Microsc 2004;214:7–12. [12] Zhang T, Yamaguchi I. Three-dimensional microscopy with phase-shifting digital holography. Opt Lett 1998;23:1221–3. [13] Tian L, Waller L. Quantitative differential phase contrast imaging in an LED array microscope. Opt Express 2015;23:11394–403. [14] Zheng G, Horstmeyer R, Yang C. Wide-field, high-resolution Fourier ptychographic microscopy. Nat Photonics 2013;7:739. [15] Debnath SK, Park Y. Real-time quantitative phase imaging with a spatial phase-shifting algorithm. Opt Lett 2011;36:4677–9. [16] Wang Y, Qiu X, Xiong J, Li B, Zhong L, Liu S, Zhou Y, Tian J, Lu X. General spatial phase-shifting interferometry by optimizing the signal retrieving function. Opt Express 2017;25:7170–80. [17] Stern A, Javidi B. Space-bandwidth conditions for efficient phase-shifting digital holographic microscopy. J Opt Soc Am A 2008;25:736–41. [18] Servín M, Quiroga JA, Padilla JM. Fringe pattern analysis for optical metrology: theory, algorithms, and applications. Wiley-VCH; 2014. [19] Morgan CJ. Least-squares estimation in phase-measurement interferometry. Opt Lett 1982;7:368–70. [20] Larkin KG, Oreb BF. Design and assessment of symmetrical phase-shifting algorithms. J Opt Soc Am A 1992;9:1740–8. [21] Hao Q, Zhu Q, Hu Y. Random phase-shifting interferometry without accurately controlling or calibrating the phase shifts. Opt Lett 2009;34:1288–90. [22] Vargas J, Sorzano COS. Quadrature component analysis for interferometry. Opt Laser Eng 2013;51:637–41. [23] Gao P, Yao B, Lindlein N, Mantel K, Harder I, Geist E. Phase-shift extraction for generalized phase-shifting interferometry. Opt Lett 2009;34:3553–5. [24] Deng J, Wang H, Zhang D, Zhong L, Fan J, Lu X. Phase shift extraction algorithm based on Euclidean matrix norm. Opt Lett 2013;38:1506–8. [25] Liu F, Wu Y, Wu F. Correction of phase extraction error in phase-shifting interferometry based on Lissajous figure and ellipse fitting technology. Opt Express 2015;23:10794–807. [26] Deng J, Wang K, Wu D, Lv X, Li C, Hao J, Qin J, Chen W. Advanced principal component analysis method for phase reconstruction. Opt Express 2015;23:12222–31. [27] Wang Z, Han B. Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms. Opt Lett 2004;29:1671–3. [28] Guo H. Blind self-calibrating algorithm for phase-shifting interferometry by use of cross-bispectrum. Opt Express 2011;19:7807–15. [29] Farrell CT, Player MA. Phase step measurement and variable step algorithms in phase-shifting interferometry. Meas Sci Technol 1992;3:953. [30] Wang Y, Lu X, Liu Y, Tian J, Zhong L. Self-calibration phase-shifting algorithm with interferograms containing very few fringes based on Fourier domain estimation. Opt Express 2017;25:29971–82. [31] Wang Y, Li B, Zhong L, Tian J, Lu X. Spatial dual-orthogonal (SDO) phase-shifting algorithm by pre-recomposing the interference fringe. Opt Express 2017;25:17446–56. [32] Vargas J, Quiroga JA, Belenguer T. Phase-shifting interferometry based on principal component analysis. Opt Lett 2011;36:1326–8. [33] Vargas J, Quiroga JA, Belenguer T. Analysis of the principal component algorithm in phase-shifting interferometry. Opt Lett 2011;36:2215–17. [34] Weijuan Q, Yingjie Y, Choo CO, Asundi A. Digital holographic microscopy with physical phase compensation. Opt Lett 2009;34:1276–8. [35] Deng J, Wu D, Wang K, Vargas J. Precise phase retrieval under harsh conditions by constructing new connected interferograms. Sci Rep 2016;6:24416. [36] Yatabe K, Ishikawa K, Oikawa Y. Improving principal component analysis based phase extraction method for phase-shifting interferometry by integrating spatial information. Opt Express 2016;24:22881–91. [37] Yatabe K, Ishikawa K, Oikawa Y. Hyper ellipse fitting in subspace method for phase-shifting interferometry: practical implementation with automatic pixel selection. Opt Express 2017;25:29401–16. [38] Yatabe K, Ishikawa K, Oikawa Y. Simple, flexible, and accurate phase re-
4. Conclusions In this study, we focus on two significant challenges in PSI algorithm, i.e. the low noise phase retrieval from few interferograms with small phase shifts and accurate self-calibration phase retrieval from USFP. By analyzing the general properties of PSI method by introducing the combination mode theory, for the first challenge, we propose an advanced noise-resistance method for PSI algorithm. To overcome the second challenge, for the first time, we propose an self-calibration phase retrieval algorithm which can realize phase retrieval from fringe patterns contain less than one fringe without strict requirement of phase shifts. Based on the above procedures, the proposed ASSF algorithm can significantly improve the quality of phase retrieval with PSI method. The Matlab codes can be downloaded at Ref [41].
Acknowledgment This work is supported by National Natural Science Foundation of China grants (61475048, 61275015 and 61177005) 180
S. Cao et al.
Optics and Lasers in Engineering 112 (2019) 170–181
trieval method for generalized phase-shifting interferometry. J Opt Soc Am A 2017;34:87–96. [39] Gonzalez RC, Woods RE. Digital image processing. 3rd Edition. Prentice-Hall; 2007.
[40] Pandey N, Ghosh A, Khare K. Two-dimensional phase unwrapping using the transport of intensity equation. Appl Opt 2016;55:2418–25. [41] https://pan.baidu.com/s/1htwqUuK
181