- Email: [email protected]
Optimal phase-shifting algorithm for interferograms with arbitrary steps and phase noise
Optics and Lasers in Engineering 114 (2019) 129–135
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Optimal phase-shifting algorithm for interferograms with arbitrary steps and phase noise Gastón A. Ayubi∗, Ignacio Duarte, José A. Ferrari Instituto de Física, Facultad de Ingeniería, UdelaR, J. Herrera y Reissig 565, Montevideo, Uruguay
a r t i c l e
i n f o
Keywords: Phase retrieval Fringe analysis Phase measurements Interferometry
a b s t r a c t Phase-shifting is a technique for phase retrieval that requires a series of intensity measurements with certain phase steps. The purpose of the present work is to study the phase extraction from phase-shifted interferograms with arbitrary steps and random phase-noise, and to deduce an algorithm that minimizes the phase-retrieval error induced by weak phase-noise. Simulations and validation experiments are presented.
[ ( )] 𝐼𝑚 (𝑥, 𝑦) = 𝐼0 (𝑥, 𝑦) 1 + 𝐾(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝛿𝑚 + 𝑛(𝑚𝑝) (𝑥, 𝑦) + 𝑛(𝑚𝑎) (𝑥, 𝑦),
1. Introduction
(3) Phase-shifting interferometry (PSI) is a well-established technique for retrieving the phase of a test object (see, e.g., [1–3] and the references therein). PSI requires a series of intensity measurements, Im (x, y) (with 𝑚 = 1, 2, … , 𝑀), with certain phase steps 𝛿 m for retrieving the phase 𝜙(x, y) of the sample. The interferograms are given by [ ( )] 𝐼𝑚 (𝑥, 𝑦) = 𝐼0 (𝑥, 𝑦) 1 + 𝐾(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝛿𝑚 , (1) where (x, y) are Cartesian coordinates, I0 (x, y) is the mean intensity, and K(x, y) is the contrast of the fringes. In general, the algorithms for the phase recovering published in the literature can be written as a quotient of linear combinations of the Im , i.e., 𝑀 ∑
tan(𝜙) =
𝑚=1 𝑀 ∑ 𝑚=1
𝑏𝑚 𝐼𝑚 ,
(2)
where 𝑛(𝑚𝑝) (𝑥, 𝑦) is the phase noise and 𝑛(𝑚𝑎) (𝑥, 𝑦) is the additive noise. An algorithm that minimizes the propagation of additive noise when the interferogram have arbitrary phase-steps was proposed in [7,8]. The propagation of phase noise or random errors in the phase steps was previously studied for algorithms with constant phase-steps [15,16]. The purpose of the present work is to determinate the optimal phaseshifting algorithm that minimizes the phase-retrieval error under weak random phase noise in the case of arbitrary phase-steps. The plan of the paper is as follows. In Section 2 we will present the generalized procedure for generating M-frame arbitrarily spaced phaseshifting algorithms. In Section 3 the performance of new algorithms for phase extraction from interferograms with random phase noise is analyzed and finally simulations and validation experiments are presented in Sections 4 and 5 respectively.
𝑎𝑚 𝐼𝑚
with given (real) coefficients am and bm . Expression (2) presupposes that the phase steps are known. In practice, however, this requirement is often difficult to exactly meet because the effective reference phases are determined not only by the phase shifter, but also by any other phenomena that changes the relative optical path difference. Thus, the actual values of 𝛿 m are only approximately known a priori [4–8], or they are to be determined a posteriori by different methods from the experimentally obtained interferograms (see, e.g., [9–14]). Noise in interferograms can be mathematically modeled as a mixture of additive noise and phase noise [2], so an interferogram is given by
2. Phase-shifting algorithm with M arbitrary steps Let us consider that the coefficients am and bm appearing in Eq. (2) are the real and imaginary parts of a complex number cm , i.e., 𝑐𝑚 = 𝑎𝑚 + 𝑖𝑏𝑚 . Then, Eq. (2) can be rewritten as (𝑀 ) ∑ 𝜙 = arg 𝑐𝑚 𝐼𝑚 , (4) 𝑚=1
where arg denotes the argument function. By using Eq. (1), one obtains 𝑀 ∑ 𝑚=1
∗
𝑐𝑚 𝐼𝑚 = 𝐼0
Corresponding author. E-mail address: ayubi@fing.edu.uy (G.A. Ayubi).
https://doi.org/10.1016/j.optlaseng.2018.10.017 Received 11 July 2018; Received in revised form 17 September 2018; Accepted 27 October 2018 0143-8166/© 2018 Elsevier Ltd. All rights reserved.
𝑀 ∑ 𝑚=1
𝑐𝑚 +
𝑀 𝑀 ∑ 𝐼0 𝐾 𝑖𝜙 ∑ 𝐼 𝐾 𝑒 𝑐𝑚 𝑒𝑖𝛿𝑚 + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖𝛿𝑚 . 2 2 𝑚=1 𝑚=1
(5)
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Hence, it is clear that, given a set of values 𝛿 m , to secure the validity of Eq. (4) the conditions to be verified by the coefficients cm are [7] 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑐𝑚 = 0,
(6)
𝑐𝑚 𝑒−𝑖𝛿𝑚 = 0,
(7)
and 𝑀 ∑ 𝑚=1
𝑐𝑚 𝑒𝑖𝛿𝑚 = 𝛼,
(8)
with 𝛼 being an arbitrary positive real number. [It can be easily demonstrated that if Condition (8) is not satisfied, the retrieved phase will differ from the true phase by a global constant.] Equivalently, in terms of am and bm , Conditions (6)–(8) can be rewritten as 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑎𝑚 = 0,
(9)
𝑏𝑚 = 0,
(10)
𝑏𝑚 cos(𝛿𝑚 ) = 0,
Fig. 1. Phasor representation of Eq. (20). 𝜎 𝜙 /2 is the maximum deviation of the retrieved phase (𝜙′) with respect to the actual phase (𝜙).
(11)
𝑎𝑚 sin(𝛿𝑚 ) = 0,
This last expression can be simplified by considering that the phase noise is small, i.e., nm ≪ 1. Then, at first order in nm one has 𝑒𝑖𝑛𝑚 ≈ 1 + 𝑖𝑛𝑚 and 𝑒−𝑖𝑛𝑚 ≈ 1 − 𝑖𝑛𝑚 , so we get
(12)
𝑀 ∑
𝑎𝑚 cos(𝛿𝑚 ) = 𝛼∕2,
(13)
𝑏𝑚 sin(𝛿𝑚 ) = −𝛼∕2.
(14)
𝑚=1
(18)
By defining 𝑆 ′ = and 𝑁 − = −𝑖
Let consider the case of interferograms with random phase noise, i.e., [ ( )] 𝐼𝑚 = 𝐼0 1 + 𝐾 cos 𝜙 + 𝛿𝑚 + 𝑛(𝑚𝑝) , (15) where 𝑛(𝑚𝑝) (𝑥, 𝑦) is the phase noise. From here on we will omit the supraindex (p) and will write nm for the phase noise. Then, the phase retrieved using Eq. (2) – which we will denote as 𝜙′(x, y) – will not be exactly the true phase 𝜙(x, y). By using Eq. (15) we get 𝑀 ∑ 𝑚=1
𝑐𝑚 +
𝑀 𝑀 𝐼0 𝐾 ∑ 𝐼 𝐾 ∑ 𝑐𝑚 𝑒𝑖(𝜙+𝛿𝑚 +𝑛𝑚 ) + 0 𝑐 𝑒−𝑖(𝜙+𝛿𝑚 +𝑛𝑚 ) 2 𝑚=1 2 𝑚=1 𝑚
(16) or, equivalently (using Eq. (6)), 𝑀 ∑ 𝑚=1
𝑐𝑚 𝐼𝑚 =
𝑀 𝑀 ∑ 𝐼0 𝐾 𝑖𝜙 ∑ 𝐼 𝐾 𝑒 𝑐𝑚 𝑒𝑖(𝛿𝑚 +𝑛𝑚 ) + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖(𝛿𝑚 +𝑛𝑚 ) 2 2 𝑚=1 𝑚=1
𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑐𝑚 𝐼𝑚 , 𝑆 =
𝐼0 𝐾 𝑖𝜙 𝑒 𝛼, 2
𝑁+ = 𝑖
𝐼0 𝐾 𝑖𝜙 𝑒 2
𝑀 ∑ 𝑚=1
(19) 𝑛𝑚 𝑐𝑚 𝑒𝑖𝛿𝑚
𝑛𝑚 𝑐𝑚 𝑒−𝑖𝛿𝑚 one gets (20)
From above definitions it is clear that the actual phase is given by 𝜙 = arg(𝑆), while the phase obtained from the acquired interferograms is given by 𝜙′ = arg(𝑆 ′ ). Also it is not difficult to demonstrate that 𝐼 𝐾𝛼 |𝑆 | = 0 (21) 2 Fig. 1 shows a phasor representation of Eq. (20) on the complex plane; the circle represents the points where it is more probable to find S′. By calculating the radius of the circle it is possible to obtain the maximum deviation of the calculated phase with respect to its actual value. We will estimate the circle radius as the square root of the expected 2 value of ||𝑁 + + 𝑁 − || . Let us denote as 𝜎 𝜙 /2 the maximum deviation of the retrieved phase (𝜙′) with respect to the actual phase (𝜙). From Fig. 1 it is clear that ( 𝜎 ) √⟨|𝑁 + +𝑁 − |2 ⟩ 𝜙 sin 2 = , and thus, if 𝜎 𝜙 is small one gets |𝑆 | ⟨ + ⟩ |𝑁 + 𝑁 − |2 𝜎𝜙2 ≈ 4 . (22) |𝑆 |2 ⟨ ⟩ We can limit 𝜎 using the triangle inequality ( |𝑁 + + 𝑁 − |2 ⩽ ⟨ + 2⟩ ⟨ − 2⟩ 𝜙 |𝑁 | + |𝑁 | ), ⟨ + 2⟩ ⟨ − 2⟩ |𝑁 | + |𝑁 | 2 𝜎𝜙 ⩽ 4 . (23) |𝑆 |2
3.1. Estimation of the phase error induced by phase noise
𝑐𝑚 𝐼𝑚 = 𝐼0
𝐼0 𝐾 −𝑖𝜙 𝑒 2
𝑆′ = 𝑆 + 𝑁+ + 𝑁−
3. Phase retrieval from interferograms with random phase noise
𝑚=1
𝑀 𝑀 ∑ ( ) 𝐼 𝐾 ( ) 𝐼0 𝐾 𝑖𝜙 ∑ 𝑒 𝑐𝑚 𝑒𝑖𝛿𝑚 1 + 𝑖𝑛𝑚 + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖𝛿𝑚 1 − 𝑖𝑛𝑚 . 2 2 𝑚=1 𝑚=1
Using conditions (7) and (8), [ ] 𝑀 𝑀 𝑀 ∑ ∑ ∑ 𝐼0 𝐾 𝑖𝜙 𝑖𝜙 𝑖𝛿𝑚 −𝑖𝜙 −𝑖𝛿𝑚 𝑐𝑚 𝐼𝑚 ≈ 𝑒 𝛼+𝑒 𝑖 𝑛𝑚 𝑐𝑚 𝑒 − 𝑒 𝑖 𝑛𝑚 𝑐𝑚 𝑒 . 2 𝑚=1 𝑚=1 𝑚=1
Thus, when we are working with M phase steps, we have to set 2M real coefficients (am , bm ) that verify the six Eqs. (9)–(14), so that we actually have 2𝑀 − 6 free coefficients. For example, in three-frame PSI (i.e., for 𝑀 = 3) we do not have free parameters to choose [besides a trivial factor multiplying all coefficients (am , bm ), which is canceled in Eq. (2)], so that for a given set of values 𝛿 m there exists a unique algorithm for phase retrieval. But for M > 3, for each set of values 𝛿 m , the phase-retrieval algorithm is not unique, and thus, from the free coefficients it is possible to choose an optimal set that minimizes the phase errors caused by random phase noise, as shown in the next section.
𝑀 ∑
𝑐𝑚 𝐼𝑚 ≈
(17) 130
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
⟨ ⟩ )2 | ∑ |2 𝑀 | | 𝑖𝛿 Since 𝑛 𝑐 𝑒 | , we can apply the triangle = | |𝑚=1 𝑚 𝑚 | | | inequality once again obtaining ( ) 𝑀 ⟨ + 2⟩ 𝐼0 𝐾 2 ∑ ⟨ 2 ⟩| |2 (24) |𝑁 | ⩽ 𝑛𝑚 |𝑐𝑚 | . 2 𝑚=1 ⟨
|𝑁 + |2
(
⟩
and, from (28) and (32),
𝐼0 𝐾 2
𝜎𝑛2 𝑏𝑚 = 𝜆2 + 𝜆3 cos(𝛿𝑚 ) + 𝜆6 sin(𝛿𝑚 ),
for 𝑚 = 1, … , 𝑀. By substituting Eqs. (33) and (34) into Eq. (30) we get a 6 × 6 linear system to calculate the Lagrange multipliers.
Assuming that the expected value of the phase noise is zero, 𝜎𝑛2 = 𝑚 ⟩ , where 𝜎𝑛𝑚 is the standard deviation of nm . Recalling that 𝑐𝑚 = 𝑎𝑚 + 𝑖𝑏𝑚 we get, ⟨
⟨
⎛ 𝑆1 ⎜ 0 ⎜ ⎜ 0 ⎜ 𝑆2 ⎜ 𝑆 ⎜ 3 ⎝ 0
𝑛2𝑚
|𝑁 |
+ 2
⟩
( ⩽
𝐼0 𝐾 2
)2
𝑀 ∑ 𝑚=1
𝜎𝑛2 𝑚
(
𝑎2𝑚
+
𝑏2𝑚
) .
(25)
⟨ ⟩ In a similar way, for |𝑁 − |2 we have ⟨
⟩ |𝑁 − |2 ⩽
(
𝐼0 𝐾 2
)2
𝑀 ∑ 𝑚=1
𝜎𝑛2
(
𝑚
) 𝑎2𝑚 + 𝑏2𝑚 .
𝜎𝜙2 ⩽ 8
𝑚=1
𝜎𝑛2 𝑚
(
𝑎2𝑚
+
𝑏2𝑚
⎧ ⎪𝑆 ⎪ 1 ⎪ ⎪𝑆 ⎪ 2 ⎪ ⎪𝑆 ⎪ 3 ⎨ ⎪𝑆 ⎪ 4 ⎪ ⎪𝑆 ⎪ 5 ⎪ ⎪𝑆 ⎪ 6 ⎩
(26)
)
𝛼2
.
(27)
3.2. Minimization of the phase-retrieval error For a given set of phase steps 𝛿 m and 𝜎𝑛𝑚 , in order to minimize the noise-induced phase error, 𝜎 𝜙 , we have to choose properly the parameters am and bm to minimize the expression 𝑓 (𝑎1 , … , 𝑎𝑀 , 𝑏1 , … , 𝑏𝑀 ) = (1∕2)
𝑀 ∑ 𝑚=1
( ) 𝜎𝑛2 𝑎2𝑚 + 𝑏2𝑚 . 𝑚
⎧𝜆1 ⎪ ⎪𝜆2 ⎪𝜆3 ⎨𝜆 ⎪ 4 ⎪𝜆5 ⎪𝜆 ⎩ 6
𝐶= (29)
𝜕𝑔𝑘 𝜕𝑓 = 𝜆 , 𝜕𝑎𝑚 𝑘=1 𝑘 𝜕𝑎𝑚
(30)
⎧ ⎪𝜆 ⎪ 1 ⎪ ⎪𝜆 ⎪ 2 ⎪ ⎪𝜆 ⎪ 3 ⎨ ⎪𝜆 ⎪ 4 ⎪ ⎪𝜆 ⎪ 5 ⎪ ⎪𝜆 ⎪ 6 ⎩
(31)
and 6 ∑
𝜕𝑔𝑘 𝜕𝑓 = 𝜆 , 𝜕𝑏𝑚 𝑘=1 𝑘 𝜕𝑏𝑚
(32)
for 𝑚 = 1, … , 𝑀. So, from (28) and (31) we have 𝜎𝑛2 𝑎𝑚 = 𝜆1 + 𝜆4 sin(𝛿𝑚 ) + 𝜆5 cos(𝛿𝑚 ), 𝑚
𝑆3 0 0 𝑆5 𝑆4 0
0 𝑆2 𝑆5 0 0 𝑆6
⎞⎛𝜆1 ⎞ ⎛ 0 ⎞ ⎟⎜𝜆2 ⎟ ⎜ 0 ⎟ ⎟⎜ ⎟ ⎜ ⎟ ⎟⎜𝜆3 ⎟ + ⎜ 0 ⎟ = 0, ⎟⎜𝜆4 ⎟ ⎜ 0 ⎟ ⎟⎜𝜆 ⎟ ⎜−𝛼∕2⎟ ⎟⎜ 5 ⎟ ⎜ ⎟ ⎠⎝𝜆6 ⎠ ⎝ 𝛼∕2 ⎠
(35)
= = = = = =
𝑀 ∑
1 𝜎2 𝑚=1 𝑛𝑚 𝑀 ∑ sin(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos2 (𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos(𝛿𝑚 ) sin(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 2 ∑ sin (𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1
(36)
( ) = 𝐶 𝑆3 𝑆6 − 𝑆2 𝑆5 , ( ) = 𝐶 𝑆3 𝑆5 − 𝑆2 𝑆4 , ( ) = 𝐶 𝑆2 𝑆3 − 𝑆1 𝑆5 , ( ) = 𝐶 𝑆1 𝑆5 − 𝑆2 𝑆3 , ( 2 ) = 𝐶 𝑆2 − 𝑆1 𝑆6 , ( ) = 𝐶 𝑆1 𝑆4 − 𝑆32 ,
(37)
𝛼∕2 𝑆4 𝑆22 − 2𝑆2 𝑆3 𝑆5 + 𝑆6 𝑆32 + 𝑆1 𝑆52 − 𝑆1 𝑆4 𝑆6
.
(38)
Thus, given a set of values 𝛿 m and 𝜎𝑛𝑚 , from Eqs. (36) we can calculate 𝑆1−6 , and then, from Eqs. (37) one obtains 𝜆1−6 . Then, by substituting these values into Eqs. (33) and (34) we calculate the optimal parameter set am and bm that minimizes the phaseretrieval error. In the case that 𝜎𝑛2 are the same for all the interferograms, the ob𝑚 tained algorithm is the same one that minimizes the propagation of the additive random noise [7,8]. In this case, substituting (36) into (37) we obtain
with 𝑘 = 1, … , 6. Applying the Lagrange multipliers method we have to solve the set of equations 6 ∑
𝑆2 0 0 𝑆6 𝑆5 0
where
The constraints are 𝑔𝑘 (𝑎1 , … , 𝑎𝑀 , 𝑏1 , … , 𝑏𝑀 ) = 0,
0 𝑆3 𝑆4 0 0 𝑆5
The solutions of the system are
(28)
The target function given by Eq. (28) must be minimized subject to Conditions (9)–(14) which can be done for instance by using Lagrange multipliers and defining ⎧ 𝑀 ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑎 1 1 𝑀 1 𝑀 𝑚 ⎪ 𝑚=1 ⎪ 𝑀 ∑ ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = 𝑏𝑚 𝑀 1 𝑀 ⎪ 2 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑏 cos(𝛿 ) 𝑀 1 𝑀 𝑚 𝑚 ⎪ 3 1 𝑚=1 ⎨ 𝑀 ∑ ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = 𝑎𝑚 sin(𝛿𝑚 ) 𝑀 1 𝑀 ⎪ 4 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑎 cos(𝛿 ) − 𝛼∕2 𝑀 1 𝑀 𝑚 𝑚 ⎪ 5 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑏 sin(𝛿 ) + 𝛼∕2. 𝑀 1 𝑀 𝑚 𝑚 ⎪ 6 1 𝑚=1 ⎩
0 𝑆1 𝑆3 0 0 𝑆2
where
Now substituting Eqs. (21), (25) and (26) into Eq. (23) we obtain 𝑀 ∑
(34)
𝑚
(33) 131
=𝐶 =𝐶 =𝐶 =𝐶 =𝐶 =𝐶
𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1
] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) sin(𝛿𝑘 ), ] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) cos(𝛿𝑙 ), ] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) sin(𝛿𝑘 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) cos(𝛿𝑙 ).
(39)
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 2. (a)Mean intensity (I0 ) function; (b) modulation (I0 K) function; (c) phase function utilized to illustrate the performance of the algorithm.
Fig. 3. Reconstructed phases for the parameters shown in Table 1.
Table 1 Sets of values am and bm for generating the phase-retrieval algorithms whose performance is illustrated in Fig. 3 (a)–(c).
And substituting (39) in (33) and (34) differing only by a global constant one gets 𝑎𝑚 =
𝑀 ∑ 𝑀 ∑ [ ][ ] sin(𝛿𝑘 ) − sin(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
Numerator coefficients
𝑘=1 𝑙=1
(40)
𝑏𝑚 =
𝑀 ∑ 𝑀 ∑ 𝑘=1 𝑙=1
Fig.3 (a) b1 b2 b3 b4
[ ][ ] cos(𝛿𝑘 ) − cos(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 ) (41)
So from (2) 𝑀 𝑀 𝑀 ∑ ∑ ∑
tan(𝜙) =
[
][
𝐼𝑚 cos(𝛿𝑘 ) − cos(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
a1 a2 a3 a4
]
𝑘=1 𝑙=1 𝑚=1 𝑀 𝑀 𝑀
∑ ∑ ∑
[ ][ ] 𝐼𝑚 sin(𝛿𝑘 ) − sin(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
Fig.3 (b)
Fig.3 (c)
1 1 −0.1317 0.0752 0.1031 0.0602 −0.9714 −1.1355 Denominator coefficients
1 −1.1869 0.3217 −0.1347
Fig.3 (a)
Fig.3 (b)
Fig.3 (c)
6.7521 6.2800 −1.1537 −11.8784
4.5106 4.1025 −0.8070 −7.8060
0.7464 −0.4051 −0.5584 0.2171
𝑘=1 𝑙=1 𝑚=1
(42)
steps 𝛿1 = 0, 𝛿2 = 𝜋∕3, 𝛿3 = 8𝜋∕7, and 𝛿4 = 13𝜋∕6, each with an arbitrary amount of Gaussian phase noise. Figs. 3 shows the retrieved phase for three different sets of am and bm parameters shown in Table 1. All these parameters verify conditions (9)–(14) but only the third set (c) was deduced from Eqs. (40) and (41). As expected Fig. 3(c) shows the least error. Fig. 4 shows the difference between the test and reconstructed phases for the three different sets of parameters (am , bm ). It is evident that the parameters that show the least error are the deduced by Eqs. (40) and (41). The phase error was calculated using 𝜙𝐸𝑟𝑟𝑜𝑟 = 𝑊 (𝜙𝐼𝑑𝑒𝑎𝑙 − 𝜙𝐶𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒𝑑 ), where 𝜙Ideal is the ideal phase used to create the simulated interferograms, 𝜙Calculated is the calculated phase, and 𝑊 (𝑥) = arctan 2[sin(𝑥), cos(𝑥)]. Fig. 5 shows the RMSE of the retrieved phase as a function of the values of a1 and b1 , with a minimum at the calculated optimal values.
This last expression is quite general: It allows to retrieve the phase (𝜙) of the test object from M interferograms (Im ) acquired by shifting the phase in arbitrary steps 𝛿 m , and simultaneously the algorithm minimizes the phase-retrieval error induced by the phase-noise 𝑛(𝑚𝑝) . Additionally, when the variances 𝜎𝑛2 are the same for all the interferograms, the ob𝑚 tained algorithm is the same that minimizes the propagation of the additive random noise 𝑛(𝑚𝑎) (see Eq. (3), and [7,8]). Note that Eq. (42) seeks to represent a complete expression of the proposed phase shifting algorithms, but in order to reduce the computational load am and bm can be calculated first from Eqs. (40) and (41) and then substituted in (2). 4. Simulations In order to illustrate the proposed error-minimization procedure with unevenly spaced phase steps, let us consider an arbitrary set of phase steps, for example 𝛿1 = 0, 𝛿2 = 𝜋∕3, 𝛿3 = 8𝜋∕7, and 𝛿4 = 13𝜋∕6. To illustrate the performance of the generated algorithm, we utilized as test mean intensity (I0 ) the function shown in Fig. 2(a), as modulation (I0 K) the function shown in Fig. 2(b) and as test phase the function shown in Fig. 2(c). We then synthesized four interferograms with phase
5. Experimental results We performed a series of validation experiments using images that corresponds to an interferometer recently proposed in [17]. As test object we used a letter etched in a ito-layer. 132
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 4. Difference between original and reconstructed phase. (a) and (b) are the difference for the retrieved phase with arbitrary parameters and (c) is for the optimal parameters deduced by Eqs. (40) and (41).
Fig. 5. RMSE of the reconstructed phase as function of a1 and b1 . Dashed lines correspond to the coordinates of the calculated optimal values. Fig. 6. Acquired interferograms with phase-steps: (a) 0 rad; 0.525 rad; 1.225 rad; 2.8 rad; 4.025 rad; 5.6 rad. Table 2 Table of coefficient that minimizes the phase-retrieved error in case of phase-steps equal to 𝛿1 = 0𝑟𝑎𝑑, 𝛿2 = 0.525𝑟𝑎𝑑 , 𝛿3 = 1.225𝑟𝑎𝑑 , 𝛿4 = 2.8𝑟𝑎𝑑 , 𝛿5 = 4.025𝑟𝑎𝑑, 𝛿6 = 5.6𝑟𝑎𝑑. Numerator coefficients b1 b2 b3 b4 b5 b6
7.5617 10.3913 8.6019 −9.9652 −16.2360 −0.3537
Denominator coefficients a1 a2 a3 a4 a5 a6
−7.6854 2.8749 15.6064 15.2117 −8.1912 −17.8165
We acquired six patterns with phase-steps 0 rad, 0.525 rad, 1.225 rad, 2.8 rad, 4.025 rad, 5.6 rad. This interferograms are shown in Fig. 6. These images were acquired with a monochromatic 8 bit digital camera (Model DC310, Thorlabs) with 1024 × 768 pixels. Although the interferograms were acquired with arbitrary phasesteps and are corrupted by noise (see Fig. 6), the phase retrieved using the optimal parameters is reasonably free of noise and consistent with the original test object (see Fig. 7). We applied the proposed algorithm to the acquired interferograms. Fig. 7 shows the obtained phase which is consistent with an ito layer of 150 nm and refractive index of 1.7.
Fig. 7. Retrieved phase of the test object.
6. Conclusions In the present work we have considered the phase retrieval from M interferograms acquired by shifting the phase by arbitrary steps 𝛿 m , with random phase-noise nm . We have demonstrated that for each set of M values 𝛿 m , we actually have 2𝑀 − 6 free parameters from which we can 133
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 8. (a)Mean intensity (I0 ) function; (b) modulation (I0 K) function; (c) phase function utilized to illustrate the performance of the algorithm.
Fig. 9. Synthesized interferograms with 𝜎𝑛 = 0 rad.
choose an optimal set (am , bm ) in order to minimize phase-retrieval errors induced by the phase-noise. By applying the Lagrange multipliers method, we have explicitly calculated (see Eqs. (40) and (41)) the set of parameters (am , bm ) that minimizes the phase-retrieval error due to the phase-noise. As an illustration of the proposed error minimization algorithm, we compared simulations obtained with the optimal parameter set of am and bm (selected by the proposed algorithm) and other arbitrary non-optimal sets. Finally we tested the proposed algorithm with data from a validation experiment with satisfactory results. Acknowledgement The authors thank the Programa de Desarrollo de Ciencias Básicas (PEDECIBA, Uruguay), the Comisión Sectorial de Investigacin Científica (CSIC, UdelaR, Uruguay), the Comisión Académica de Posgrado (CAP, UdelaR, Uruguay) and the Agencia Nacional de Investigación e Innovación (ANII, Uruguay) for their financial support.
Fig. 10. RMSE of the reconstructed phase as function of the standard deviation of the phase noise (𝜎 n ).
noise ranges from 𝜎𝑛 = 0 rad to 𝜎𝑛 = 0.25 rad. Fig. 9(a)–(f) shows the synthesized interferograms with 𝜎𝑛 = 0 rad. From the utilized phase-steps, we obtained the coefficients of the optimal algorithm using the proposed method. Then we applied the obtained algorithm to the interferograms and reconstructed the phase for each set of interferograms. Fig. 10 shows the root mean square error (RMSE) as a function of the standard deviation of the phase noise (𝜎 n ), where the linearity in the function dependence can clearly be seen. Finally Fig. 11(a)–(c) shows the retrieved phase for three different sets of interferograms with different level of noise. The circle marks in Fig. 10 corresponds with the results in Fig. 11.
Appendix A In the following appendix we present additional simulations in order to show the dependence between phase noise levels and error of the reconstructed phase. With this purpose in mind, we utilized as test mean intensity (I0 ) the function shown in Fig. 8(a), as modulation (I0 K) the function shown in Fig. 8(b) and as test phase the function shown in Fig. 8(c) plus a linear term. We synthesized sets of six interferograms with different level of Gaussian phase noise, with phase-steps 0 rad, 0.1571 rad, 0.5411 rad, 2.7402 rad, 4.9393 rad, and 6.0912 rad. The standard deviation 𝜎 n of the phase 134
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 11. Reconstructed phases for different noise level.
References
[9] Okada K, Sato A, Tsujiuchi J. Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry. Opt Commun 1991;84(3):118–24. [10] Farrell CT, Player MA. Phase step measurement and variable step algorithms in phase-shifting interferometry. Meas Sci Technol 1992;3:953–8. [11] Wang Z, Han B. Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms. Opt Lett 2004;29:1671–3. [12] Gao P, Yao B, Lindlein N, Mantel K, Harder I, Geist E. Phase-shift extraction for generalized phase-shifting interferometry. Opt Lett 2009;34(22):3553–5. [13] Guo H, Zhang Z. Phase shift estimation from variances of fringe pattern differences. Appl Opt 2013;52(26):6572–8. [14] Ding H, Guo H. Estimating phase shifts from three fringe patterns by use of cross spectrum. Appl Opt 2017;56(4). [15] Rathjen C. Statistical properties of phase-shift algorithms. J Opt Soc Am A 1995;12(9):1997–2008. [16] Hibino K. Susceptibility of systematic error-compensating algorithms to random noise in phase-shifting interferometry. Appl Opt 1997;36(10):2084–93. [17] Ferrari JA, Flores JL, Ayubi GA, Frins E. Robust interferometer with external phase-shift control. Opt Express 2017;25(24):29965–70.
[1] Malacara D. Optical shop testing. Wiley; 2007. [2] Servin M, Quiroga JA, Padilla JM. Fringe pattern analysis for optical metrology. Wiley-VCH; 2014. [3] Servin M, Estrada JC, Quiroga JA. The general theory of phase shifting algorithms. Opt Express 2009;17:21867–81. [4] Lai G, Yatagai T. Generalized phase-shifting interferometry. J Opt Soc Am A 1991;8:822–7. [5] Téllez-Quinones A, Malacara-Doblado D. Inhomogeneous phase shifting: an algorithm for nonconstant phase displacements. Appl Opt 2010;49:6224–31. [6] Ayubi GA, Perciante CD, Flores JL, Martino JMD, Ferrari JA. Generation of phase-shifting algorithms with n arbitrarily spaced phase-steps. Appl Opt 2014;53:7168–76. [7] Ayubi GA, Perciante CD, Martino JMD, Flores JL, Ferrari JA. Generalized phase-shifting algorithms: error analysis and minimization of noise propagation. Appl Opt 2016;55(6):1461–9. [8] Ayubi GA, Perciante CD, Martino JMD, Flores JL, Ferrari JA. Generalized phase-shifting algorithms: error analysis and minimization of noise propagation: erratum. Appl Opt 2016;55(28):7763.
135
Contents lists available at ScienceDirect
Optics and Lasers in Engineering journal homepage: www.elsevier.com/locate/optlaseng
Optimal phase-shifting algorithm for interferograms with arbitrary steps and phase noise Gastón A. Ayubi∗, Ignacio Duarte, José A. Ferrari Instituto de Física, Facultad de Ingeniería, UdelaR, J. Herrera y Reissig 565, Montevideo, Uruguay
a r t i c l e
i n f o
Keywords: Phase retrieval Fringe analysis Phase measurements Interferometry
a b s t r a c t Phase-shifting is a technique for phase retrieval that requires a series of intensity measurements with certain phase steps. The purpose of the present work is to study the phase extraction from phase-shifted interferograms with arbitrary steps and random phase-noise, and to deduce an algorithm that minimizes the phase-retrieval error induced by weak phase-noise. Simulations and validation experiments are presented.
[ ( )] 𝐼𝑚 (𝑥, 𝑦) = 𝐼0 (𝑥, 𝑦) 1 + 𝐾(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝛿𝑚 + 𝑛(𝑚𝑝) (𝑥, 𝑦) + 𝑛(𝑚𝑎) (𝑥, 𝑦),
1. Introduction
(3) Phase-shifting interferometry (PSI) is a well-established technique for retrieving the phase of a test object (see, e.g., [1–3] and the references therein). PSI requires a series of intensity measurements, Im (x, y) (with 𝑚 = 1, 2, … , 𝑀), with certain phase steps 𝛿 m for retrieving the phase 𝜙(x, y) of the sample. The interferograms are given by [ ( )] 𝐼𝑚 (𝑥, 𝑦) = 𝐼0 (𝑥, 𝑦) 1 + 𝐾(𝑥, 𝑦) cos 𝜙(𝑥, 𝑦) + 𝛿𝑚 , (1) where (x, y) are Cartesian coordinates, I0 (x, y) is the mean intensity, and K(x, y) is the contrast of the fringes. In general, the algorithms for the phase recovering published in the literature can be written as a quotient of linear combinations of the Im , i.e., 𝑀 ∑
tan(𝜙) =
𝑚=1 𝑀 ∑ 𝑚=1
𝑏𝑚 𝐼𝑚 ,
(2)
where 𝑛(𝑚𝑝) (𝑥, 𝑦) is the phase noise and 𝑛(𝑚𝑎) (𝑥, 𝑦) is the additive noise. An algorithm that minimizes the propagation of additive noise when the interferogram have arbitrary phase-steps was proposed in [7,8]. The propagation of phase noise or random errors in the phase steps was previously studied for algorithms with constant phase-steps [15,16]. The purpose of the present work is to determinate the optimal phaseshifting algorithm that minimizes the phase-retrieval error under weak random phase noise in the case of arbitrary phase-steps. The plan of the paper is as follows. In Section 2 we will present the generalized procedure for generating M-frame arbitrarily spaced phaseshifting algorithms. In Section 3 the performance of new algorithms for phase extraction from interferograms with random phase noise is analyzed and finally simulations and validation experiments are presented in Sections 4 and 5 respectively.
𝑎𝑚 𝐼𝑚
with given (real) coefficients am and bm . Expression (2) presupposes that the phase steps are known. In practice, however, this requirement is often difficult to exactly meet because the effective reference phases are determined not only by the phase shifter, but also by any other phenomena that changes the relative optical path difference. Thus, the actual values of 𝛿 m are only approximately known a priori [4–8], or they are to be determined a posteriori by different methods from the experimentally obtained interferograms (see, e.g., [9–14]). Noise in interferograms can be mathematically modeled as a mixture of additive noise and phase noise [2], so an interferogram is given by
2. Phase-shifting algorithm with M arbitrary steps Let us consider that the coefficients am and bm appearing in Eq. (2) are the real and imaginary parts of a complex number cm , i.e., 𝑐𝑚 = 𝑎𝑚 + 𝑖𝑏𝑚 . Then, Eq. (2) can be rewritten as (𝑀 ) ∑ 𝜙 = arg 𝑐𝑚 𝐼𝑚 , (4) 𝑚=1
where arg denotes the argument function. By using Eq. (1), one obtains 𝑀 ∑ 𝑚=1
∗
𝑐𝑚 𝐼𝑚 = 𝐼0
Corresponding author. E-mail address: ayubi@fing.edu.uy (G.A. Ayubi).
https://doi.org/10.1016/j.optlaseng.2018.10.017 Received 11 July 2018; Received in revised form 17 September 2018; Accepted 27 October 2018 0143-8166/© 2018 Elsevier Ltd. All rights reserved.
𝑀 ∑ 𝑚=1
𝑐𝑚 +
𝑀 𝑀 ∑ 𝐼0 𝐾 𝑖𝜙 ∑ 𝐼 𝐾 𝑒 𝑐𝑚 𝑒𝑖𝛿𝑚 + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖𝛿𝑚 . 2 2 𝑚=1 𝑚=1
(5)
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Hence, it is clear that, given a set of values 𝛿 m , to secure the validity of Eq. (4) the conditions to be verified by the coefficients cm are [7] 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑐𝑚 = 0,
(6)
𝑐𝑚 𝑒−𝑖𝛿𝑚 = 0,
(7)
and 𝑀 ∑ 𝑚=1
𝑐𝑚 𝑒𝑖𝛿𝑚 = 𝛼,
(8)
with 𝛼 being an arbitrary positive real number. [It can be easily demonstrated that if Condition (8) is not satisfied, the retrieved phase will differ from the true phase by a global constant.] Equivalently, in terms of am and bm , Conditions (6)–(8) can be rewritten as 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑎𝑚 = 0,
(9)
𝑏𝑚 = 0,
(10)
𝑏𝑚 cos(𝛿𝑚 ) = 0,
Fig. 1. Phasor representation of Eq. (20). 𝜎 𝜙 /2 is the maximum deviation of the retrieved phase (𝜙′) with respect to the actual phase (𝜙).
(11)
𝑎𝑚 sin(𝛿𝑚 ) = 0,
This last expression can be simplified by considering that the phase noise is small, i.e., nm ≪ 1. Then, at first order in nm one has 𝑒𝑖𝑛𝑚 ≈ 1 + 𝑖𝑛𝑚 and 𝑒−𝑖𝑛𝑚 ≈ 1 − 𝑖𝑛𝑚 , so we get
(12)
𝑀 ∑
𝑎𝑚 cos(𝛿𝑚 ) = 𝛼∕2,
(13)
𝑏𝑚 sin(𝛿𝑚 ) = −𝛼∕2.
(14)
𝑚=1
(18)
By defining 𝑆 ′ = and 𝑁 − = −𝑖
Let consider the case of interferograms with random phase noise, i.e., [ ( )] 𝐼𝑚 = 𝐼0 1 + 𝐾 cos 𝜙 + 𝛿𝑚 + 𝑛(𝑚𝑝) , (15) where 𝑛(𝑚𝑝) (𝑥, 𝑦) is the phase noise. From here on we will omit the supraindex (p) and will write nm for the phase noise. Then, the phase retrieved using Eq. (2) – which we will denote as 𝜙′(x, y) – will not be exactly the true phase 𝜙(x, y). By using Eq. (15) we get 𝑀 ∑ 𝑚=1
𝑐𝑚 +
𝑀 𝑀 𝐼0 𝐾 ∑ 𝐼 𝐾 ∑ 𝑐𝑚 𝑒𝑖(𝜙+𝛿𝑚 +𝑛𝑚 ) + 0 𝑐 𝑒−𝑖(𝜙+𝛿𝑚 +𝑛𝑚 ) 2 𝑚=1 2 𝑚=1 𝑚
(16) or, equivalently (using Eq. (6)), 𝑀 ∑ 𝑚=1
𝑐𝑚 𝐼𝑚 =
𝑀 𝑀 ∑ 𝐼0 𝐾 𝑖𝜙 ∑ 𝐼 𝐾 𝑒 𝑐𝑚 𝑒𝑖(𝛿𝑚 +𝑛𝑚 ) + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖(𝛿𝑚 +𝑛𝑚 ) 2 2 𝑚=1 𝑚=1
𝑀 ∑ 𝑚=1 𝑀 ∑ 𝑚=1
𝑐𝑚 𝐼𝑚 , 𝑆 =
𝐼0 𝐾 𝑖𝜙 𝑒 𝛼, 2
𝑁+ = 𝑖
𝐼0 𝐾 𝑖𝜙 𝑒 2
𝑀 ∑ 𝑚=1
(19) 𝑛𝑚 𝑐𝑚 𝑒𝑖𝛿𝑚
𝑛𝑚 𝑐𝑚 𝑒−𝑖𝛿𝑚 one gets (20)
From above definitions it is clear that the actual phase is given by 𝜙 = arg(𝑆), while the phase obtained from the acquired interferograms is given by 𝜙′ = arg(𝑆 ′ ). Also it is not difficult to demonstrate that 𝐼 𝐾𝛼 |𝑆 | = 0 (21) 2 Fig. 1 shows a phasor representation of Eq. (20) on the complex plane; the circle represents the points where it is more probable to find S′. By calculating the radius of the circle it is possible to obtain the maximum deviation of the calculated phase with respect to its actual value. We will estimate the circle radius as the square root of the expected 2 value of ||𝑁 + + 𝑁 − || . Let us denote as 𝜎 𝜙 /2 the maximum deviation of the retrieved phase (𝜙′) with respect to the actual phase (𝜙). From Fig. 1 it is clear that ( 𝜎 ) √⟨|𝑁 + +𝑁 − |2 ⟩ 𝜙 sin 2 = , and thus, if 𝜎 𝜙 is small one gets |𝑆 | ⟨ + ⟩ |𝑁 + 𝑁 − |2 𝜎𝜙2 ≈ 4 . (22) |𝑆 |2 ⟨ ⟩ We can limit 𝜎 using the triangle inequality ( |𝑁 + + 𝑁 − |2 ⩽ ⟨ + 2⟩ ⟨ − 2⟩ 𝜙 |𝑁 | + |𝑁 | ), ⟨ + 2⟩ ⟨ − 2⟩ |𝑁 | + |𝑁 | 2 𝜎𝜙 ⩽ 4 . (23) |𝑆 |2
3.1. Estimation of the phase error induced by phase noise
𝑐𝑚 𝐼𝑚 = 𝐼0
𝐼0 𝐾 −𝑖𝜙 𝑒 2
𝑆′ = 𝑆 + 𝑁+ + 𝑁−
3. Phase retrieval from interferograms with random phase noise
𝑚=1
𝑀 𝑀 ∑ ( ) 𝐼 𝐾 ( ) 𝐼0 𝐾 𝑖𝜙 ∑ 𝑒 𝑐𝑚 𝑒𝑖𝛿𝑚 1 + 𝑖𝑛𝑚 + 0 𝑒−𝑖𝜙 𝑐𝑚 𝑒−𝑖𝛿𝑚 1 − 𝑖𝑛𝑚 . 2 2 𝑚=1 𝑚=1
Using conditions (7) and (8), [ ] 𝑀 𝑀 𝑀 ∑ ∑ ∑ 𝐼0 𝐾 𝑖𝜙 𝑖𝜙 𝑖𝛿𝑚 −𝑖𝜙 −𝑖𝛿𝑚 𝑐𝑚 𝐼𝑚 ≈ 𝑒 𝛼+𝑒 𝑖 𝑛𝑚 𝑐𝑚 𝑒 − 𝑒 𝑖 𝑛𝑚 𝑐𝑚 𝑒 . 2 𝑚=1 𝑚=1 𝑚=1
Thus, when we are working with M phase steps, we have to set 2M real coefficients (am , bm ) that verify the six Eqs. (9)–(14), so that we actually have 2𝑀 − 6 free coefficients. For example, in three-frame PSI (i.e., for 𝑀 = 3) we do not have free parameters to choose [besides a trivial factor multiplying all coefficients (am , bm ), which is canceled in Eq. (2)], so that for a given set of values 𝛿 m there exists a unique algorithm for phase retrieval. But for M > 3, for each set of values 𝛿 m , the phase-retrieval algorithm is not unique, and thus, from the free coefficients it is possible to choose an optimal set that minimizes the phase errors caused by random phase noise, as shown in the next section.
𝑀 ∑
𝑐𝑚 𝐼𝑚 ≈
(17) 130
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
⟨ ⟩ )2 | ∑ |2 𝑀 | | 𝑖𝛿 Since 𝑛 𝑐 𝑒 | , we can apply the triangle = | |𝑚=1 𝑚 𝑚 | | | inequality once again obtaining ( ) 𝑀 ⟨ + 2⟩ 𝐼0 𝐾 2 ∑ ⟨ 2 ⟩| |2 (24) |𝑁 | ⩽ 𝑛𝑚 |𝑐𝑚 | . 2 𝑚=1 ⟨
|𝑁 + |2
(
⟩
and, from (28) and (32),
𝐼0 𝐾 2
𝜎𝑛2 𝑏𝑚 = 𝜆2 + 𝜆3 cos(𝛿𝑚 ) + 𝜆6 sin(𝛿𝑚 ),
for 𝑚 = 1, … , 𝑀. By substituting Eqs. (33) and (34) into Eq. (30) we get a 6 × 6 linear system to calculate the Lagrange multipliers.
Assuming that the expected value of the phase noise is zero, 𝜎𝑛2 = 𝑚 ⟩ , where 𝜎𝑛𝑚 is the standard deviation of nm . Recalling that 𝑐𝑚 = 𝑎𝑚 + 𝑖𝑏𝑚 we get, ⟨
⟨
⎛ 𝑆1 ⎜ 0 ⎜ ⎜ 0 ⎜ 𝑆2 ⎜ 𝑆 ⎜ 3 ⎝ 0
𝑛2𝑚
|𝑁 |
+ 2
⟩
( ⩽
𝐼0 𝐾 2
)2
𝑀 ∑ 𝑚=1
𝜎𝑛2 𝑚
(
𝑎2𝑚
+
𝑏2𝑚
) .
(25)
⟨ ⟩ In a similar way, for |𝑁 − |2 we have ⟨
⟩ |𝑁 − |2 ⩽
(
𝐼0 𝐾 2
)2
𝑀 ∑ 𝑚=1
𝜎𝑛2
(
𝑚
) 𝑎2𝑚 + 𝑏2𝑚 .
𝜎𝜙2 ⩽ 8
𝑚=1
𝜎𝑛2 𝑚
(
𝑎2𝑚
+
𝑏2𝑚
⎧ ⎪𝑆 ⎪ 1 ⎪ ⎪𝑆 ⎪ 2 ⎪ ⎪𝑆 ⎪ 3 ⎨ ⎪𝑆 ⎪ 4 ⎪ ⎪𝑆 ⎪ 5 ⎪ ⎪𝑆 ⎪ 6 ⎩
(26)
)
𝛼2
.
(27)
3.2. Minimization of the phase-retrieval error For a given set of phase steps 𝛿 m and 𝜎𝑛𝑚 , in order to minimize the noise-induced phase error, 𝜎 𝜙 , we have to choose properly the parameters am and bm to minimize the expression 𝑓 (𝑎1 , … , 𝑎𝑀 , 𝑏1 , … , 𝑏𝑀 ) = (1∕2)
𝑀 ∑ 𝑚=1
( ) 𝜎𝑛2 𝑎2𝑚 + 𝑏2𝑚 . 𝑚
⎧𝜆1 ⎪ ⎪𝜆2 ⎪𝜆3 ⎨𝜆 ⎪ 4 ⎪𝜆5 ⎪𝜆 ⎩ 6
𝐶= (29)
𝜕𝑔𝑘 𝜕𝑓 = 𝜆 , 𝜕𝑎𝑚 𝑘=1 𝑘 𝜕𝑎𝑚
(30)
⎧ ⎪𝜆 ⎪ 1 ⎪ ⎪𝜆 ⎪ 2 ⎪ ⎪𝜆 ⎪ 3 ⎨ ⎪𝜆 ⎪ 4 ⎪ ⎪𝜆 ⎪ 5 ⎪ ⎪𝜆 ⎪ 6 ⎩
(31)
and 6 ∑
𝜕𝑔𝑘 𝜕𝑓 = 𝜆 , 𝜕𝑏𝑚 𝑘=1 𝑘 𝜕𝑏𝑚
(32)
for 𝑚 = 1, … , 𝑀. So, from (28) and (31) we have 𝜎𝑛2 𝑎𝑚 = 𝜆1 + 𝜆4 sin(𝛿𝑚 ) + 𝜆5 cos(𝛿𝑚 ), 𝑚
𝑆3 0 0 𝑆5 𝑆4 0
0 𝑆2 𝑆5 0 0 𝑆6
⎞⎛𝜆1 ⎞ ⎛ 0 ⎞ ⎟⎜𝜆2 ⎟ ⎜ 0 ⎟ ⎟⎜ ⎟ ⎜ ⎟ ⎟⎜𝜆3 ⎟ + ⎜ 0 ⎟ = 0, ⎟⎜𝜆4 ⎟ ⎜ 0 ⎟ ⎟⎜𝜆 ⎟ ⎜−𝛼∕2⎟ ⎟⎜ 5 ⎟ ⎜ ⎟ ⎠⎝𝜆6 ⎠ ⎝ 𝛼∕2 ⎠
(35)
= = = = = =
𝑀 ∑
1 𝜎2 𝑚=1 𝑛𝑚 𝑀 ∑ sin(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos2 (𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 ∑ cos(𝛿𝑚 ) sin(𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1 𝑀 2 ∑ sin (𝛿𝑚 ) 𝜎𝑛2𝑚 𝑚=1
(36)
( ) = 𝐶 𝑆3 𝑆6 − 𝑆2 𝑆5 , ( ) = 𝐶 𝑆3 𝑆5 − 𝑆2 𝑆4 , ( ) = 𝐶 𝑆2 𝑆3 − 𝑆1 𝑆5 , ( ) = 𝐶 𝑆1 𝑆5 − 𝑆2 𝑆3 , ( 2 ) = 𝐶 𝑆2 − 𝑆1 𝑆6 , ( ) = 𝐶 𝑆1 𝑆4 − 𝑆32 ,
(37)
𝛼∕2 𝑆4 𝑆22 − 2𝑆2 𝑆3 𝑆5 + 𝑆6 𝑆32 + 𝑆1 𝑆52 − 𝑆1 𝑆4 𝑆6
.
(38)
Thus, given a set of values 𝛿 m and 𝜎𝑛𝑚 , from Eqs. (36) we can calculate 𝑆1−6 , and then, from Eqs. (37) one obtains 𝜆1−6 . Then, by substituting these values into Eqs. (33) and (34) we calculate the optimal parameter set am and bm that minimizes the phaseretrieval error. In the case that 𝜎𝑛2 are the same for all the interferograms, the ob𝑚 tained algorithm is the same one that minimizes the propagation of the additive random noise [7,8]. In this case, substituting (36) into (37) we obtain
with 𝑘 = 1, … , 6. Applying the Lagrange multipliers method we have to solve the set of equations 6 ∑
𝑆2 0 0 𝑆6 𝑆5 0
where
The constraints are 𝑔𝑘 (𝑎1 , … , 𝑎𝑀 , 𝑏1 , … , 𝑏𝑀 ) = 0,
0 𝑆3 𝑆4 0 0 𝑆5
The solutions of the system are
(28)
The target function given by Eq. (28) must be minimized subject to Conditions (9)–(14) which can be done for instance by using Lagrange multipliers and defining ⎧ 𝑀 ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑎 1 1 𝑀 1 𝑀 𝑚 ⎪ 𝑚=1 ⎪ 𝑀 ∑ ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = 𝑏𝑚 𝑀 1 𝑀 ⎪ 2 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑏 cos(𝛿 ) 𝑀 1 𝑀 𝑚 𝑚 ⎪ 3 1 𝑚=1 ⎨ 𝑀 ∑ ⎪𝑔 ( 𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = 𝑎𝑚 sin(𝛿𝑚 ) 𝑀 1 𝑀 ⎪ 4 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑎 cos(𝛿 ) − 𝛼∕2 𝑀 1 𝑀 𝑚 𝑚 ⎪ 5 1 𝑚=1 ⎪ 𝑀 ⎪𝑔 (𝑎 , … , 𝑎 , 𝑏 , … , 𝑏 ) = ∑ 𝑏 sin(𝛿 ) + 𝛼∕2. 𝑀 1 𝑀 𝑚 𝑚 ⎪ 6 1 𝑚=1 ⎩
0 𝑆1 𝑆3 0 0 𝑆2
where
Now substituting Eqs. (21), (25) and (26) into Eq. (23) we obtain 𝑀 ∑
(34)
𝑚
(33) 131
=𝐶 =𝐶 =𝐶 =𝐶 =𝐶 =𝐶
𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1 𝑀 ∑ 𝑀 [ ∑ 𝑘=1 𝑙=1
] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) sin(𝛿𝑘 ), ] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) cos(𝛿𝑙 ), ] sin(𝛿𝑙 ) − sin(𝛿𝑘 ) sin(𝛿𝑘 ), ] cos(𝛿𝑙 ) − cos(𝛿𝑘 ) cos(𝛿𝑙 ).
(39)
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 2. (a)Mean intensity (I0 ) function; (b) modulation (I0 K) function; (c) phase function utilized to illustrate the performance of the algorithm.
Fig. 3. Reconstructed phases for the parameters shown in Table 1.
Table 1 Sets of values am and bm for generating the phase-retrieval algorithms whose performance is illustrated in Fig. 3 (a)–(c).
And substituting (39) in (33) and (34) differing only by a global constant one gets 𝑎𝑚 =
𝑀 ∑ 𝑀 ∑ [ ][ ] sin(𝛿𝑘 ) − sin(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
Numerator coefficients
𝑘=1 𝑙=1
(40)
𝑏𝑚 =
𝑀 ∑ 𝑀 ∑ 𝑘=1 𝑙=1
Fig.3 (a) b1 b2 b3 b4
[ ][ ] cos(𝛿𝑘 ) − cos(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 ) (41)
So from (2) 𝑀 𝑀 𝑀 ∑ ∑ ∑
tan(𝜙) =
[
][
𝐼𝑚 cos(𝛿𝑘 ) − cos(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
a1 a2 a3 a4
]
𝑘=1 𝑙=1 𝑚=1 𝑀 𝑀 𝑀
∑ ∑ ∑
[ ][ ] 𝐼𝑚 sin(𝛿𝑘 ) − sin(𝛿𝑙 ) cos(𝛿𝑘 ) sin(𝛿𝑙 ) + cos(𝛿𝑚 ) sin(𝛿𝑘 ) + cos(𝛿𝑙 ) sin(𝛿𝑚 )
Fig.3 (b)
Fig.3 (c)
1 1 −0.1317 0.0752 0.1031 0.0602 −0.9714 −1.1355 Denominator coefficients
1 −1.1869 0.3217 −0.1347
Fig.3 (a)
Fig.3 (b)
Fig.3 (c)
6.7521 6.2800 −1.1537 −11.8784
4.5106 4.1025 −0.8070 −7.8060
0.7464 −0.4051 −0.5584 0.2171
𝑘=1 𝑙=1 𝑚=1
(42)
steps 𝛿1 = 0, 𝛿2 = 𝜋∕3, 𝛿3 = 8𝜋∕7, and 𝛿4 = 13𝜋∕6, each with an arbitrary amount of Gaussian phase noise. Figs. 3 shows the retrieved phase for three different sets of am and bm parameters shown in Table 1. All these parameters verify conditions (9)–(14) but only the third set (c) was deduced from Eqs. (40) and (41). As expected Fig. 3(c) shows the least error. Fig. 4 shows the difference between the test and reconstructed phases for the three different sets of parameters (am , bm ). It is evident that the parameters that show the least error are the deduced by Eqs. (40) and (41). The phase error was calculated using 𝜙𝐸𝑟𝑟𝑜𝑟 = 𝑊 (𝜙𝐼𝑑𝑒𝑎𝑙 − 𝜙𝐶𝑎𝑙𝑐𝑢𝑙𝑎𝑡𝑒𝑑 ), where 𝜙Ideal is the ideal phase used to create the simulated interferograms, 𝜙Calculated is the calculated phase, and 𝑊 (𝑥) = arctan 2[sin(𝑥), cos(𝑥)]. Fig. 5 shows the RMSE of the retrieved phase as a function of the values of a1 and b1 , with a minimum at the calculated optimal values.
This last expression is quite general: It allows to retrieve the phase (𝜙) of the test object from M interferograms (Im ) acquired by shifting the phase in arbitrary steps 𝛿 m , and simultaneously the algorithm minimizes the phase-retrieval error induced by the phase-noise 𝑛(𝑚𝑝) . Additionally, when the variances 𝜎𝑛2 are the same for all the interferograms, the ob𝑚 tained algorithm is the same that minimizes the propagation of the additive random noise 𝑛(𝑚𝑎) (see Eq. (3), and [7,8]). Note that Eq. (42) seeks to represent a complete expression of the proposed phase shifting algorithms, but in order to reduce the computational load am and bm can be calculated first from Eqs. (40) and (41) and then substituted in (2). 4. Simulations In order to illustrate the proposed error-minimization procedure with unevenly spaced phase steps, let us consider an arbitrary set of phase steps, for example 𝛿1 = 0, 𝛿2 = 𝜋∕3, 𝛿3 = 8𝜋∕7, and 𝛿4 = 13𝜋∕6. To illustrate the performance of the generated algorithm, we utilized as test mean intensity (I0 ) the function shown in Fig. 2(a), as modulation (I0 K) the function shown in Fig. 2(b) and as test phase the function shown in Fig. 2(c). We then synthesized four interferograms with phase
5. Experimental results We performed a series of validation experiments using images that corresponds to an interferometer recently proposed in [17]. As test object we used a letter etched in a ito-layer. 132
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 4. Difference between original and reconstructed phase. (a) and (b) are the difference for the retrieved phase with arbitrary parameters and (c) is for the optimal parameters deduced by Eqs. (40) and (41).
Fig. 5. RMSE of the reconstructed phase as function of a1 and b1 . Dashed lines correspond to the coordinates of the calculated optimal values. Fig. 6. Acquired interferograms with phase-steps: (a) 0 rad; 0.525 rad; 1.225 rad; 2.8 rad; 4.025 rad; 5.6 rad. Table 2 Table of coefficient that minimizes the phase-retrieved error in case of phase-steps equal to 𝛿1 = 0𝑟𝑎𝑑, 𝛿2 = 0.525𝑟𝑎𝑑 , 𝛿3 = 1.225𝑟𝑎𝑑 , 𝛿4 = 2.8𝑟𝑎𝑑 , 𝛿5 = 4.025𝑟𝑎𝑑, 𝛿6 = 5.6𝑟𝑎𝑑. Numerator coefficients b1 b2 b3 b4 b5 b6
7.5617 10.3913 8.6019 −9.9652 −16.2360 −0.3537
Denominator coefficients a1 a2 a3 a4 a5 a6
−7.6854 2.8749 15.6064 15.2117 −8.1912 −17.8165
We acquired six patterns with phase-steps 0 rad, 0.525 rad, 1.225 rad, 2.8 rad, 4.025 rad, 5.6 rad. This interferograms are shown in Fig. 6. These images were acquired with a monochromatic 8 bit digital camera (Model DC310, Thorlabs) with 1024 × 768 pixels. Although the interferograms were acquired with arbitrary phasesteps and are corrupted by noise (see Fig. 6), the phase retrieved using the optimal parameters is reasonably free of noise and consistent with the original test object (see Fig. 7). We applied the proposed algorithm to the acquired interferograms. Fig. 7 shows the obtained phase which is consistent with an ito layer of 150 nm and refractive index of 1.7.
Fig. 7. Retrieved phase of the test object.
6. Conclusions In the present work we have considered the phase retrieval from M interferograms acquired by shifting the phase by arbitrary steps 𝛿 m , with random phase-noise nm . We have demonstrated that for each set of M values 𝛿 m , we actually have 2𝑀 − 6 free parameters from which we can 133
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 8. (a)Mean intensity (I0 ) function; (b) modulation (I0 K) function; (c) phase function utilized to illustrate the performance of the algorithm.
Fig. 9. Synthesized interferograms with 𝜎𝑛 = 0 rad.
choose an optimal set (am , bm ) in order to minimize phase-retrieval errors induced by the phase-noise. By applying the Lagrange multipliers method, we have explicitly calculated (see Eqs. (40) and (41)) the set of parameters (am , bm ) that minimizes the phase-retrieval error due to the phase-noise. As an illustration of the proposed error minimization algorithm, we compared simulations obtained with the optimal parameter set of am and bm (selected by the proposed algorithm) and other arbitrary non-optimal sets. Finally we tested the proposed algorithm with data from a validation experiment with satisfactory results. Acknowledgement The authors thank the Programa de Desarrollo de Ciencias Básicas (PEDECIBA, Uruguay), the Comisión Sectorial de Investigacin Científica (CSIC, UdelaR, Uruguay), the Comisión Académica de Posgrado (CAP, UdelaR, Uruguay) and the Agencia Nacional de Investigación e Innovación (ANII, Uruguay) for their financial support.
Fig. 10. RMSE of the reconstructed phase as function of the standard deviation of the phase noise (𝜎 n ).
noise ranges from 𝜎𝑛 = 0 rad to 𝜎𝑛 = 0.25 rad. Fig. 9(a)–(f) shows the synthesized interferograms with 𝜎𝑛 = 0 rad. From the utilized phase-steps, we obtained the coefficients of the optimal algorithm using the proposed method. Then we applied the obtained algorithm to the interferograms and reconstructed the phase for each set of interferograms. Fig. 10 shows the root mean square error (RMSE) as a function of the standard deviation of the phase noise (𝜎 n ), where the linearity in the function dependence can clearly be seen. Finally Fig. 11(a)–(c) shows the retrieved phase for three different sets of interferograms with different level of noise. The circle marks in Fig. 10 corresponds with the results in Fig. 11.
Appendix A In the following appendix we present additional simulations in order to show the dependence between phase noise levels and error of the reconstructed phase. With this purpose in mind, we utilized as test mean intensity (I0 ) the function shown in Fig. 8(a), as modulation (I0 K) the function shown in Fig. 8(b) and as test phase the function shown in Fig. 8(c) plus a linear term. We synthesized sets of six interferograms with different level of Gaussian phase noise, with phase-steps 0 rad, 0.1571 rad, 0.5411 rad, 2.7402 rad, 4.9393 rad, and 6.0912 rad. The standard deviation 𝜎 n of the phase 134
G.A. Ayubi, I. Duarte and J.A. Ferrari
Optics and Lasers in Engineering 114 (2019) 129–135
Fig. 11. Reconstructed phases for different noise level.
References
[9] Okada K, Sato A, Tsujiuchi J. Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry. Opt Commun 1991;84(3):118–24. [10] Farrell CT, Player MA. Phase step measurement and variable step algorithms in phase-shifting interferometry. Meas Sci Technol 1992;3:953–8. [11] Wang Z, Han B. Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms. Opt Lett 2004;29:1671–3. [12] Gao P, Yao B, Lindlein N, Mantel K, Harder I, Geist E. Phase-shift extraction for generalized phase-shifting interferometry. Opt Lett 2009;34(22):3553–5. [13] Guo H, Zhang Z. Phase shift estimation from variances of fringe pattern differences. Appl Opt 2013;52(26):6572–8. [14] Ding H, Guo H. Estimating phase shifts from three fringe patterns by use of cross spectrum. Appl Opt 2017;56(4). [15] Rathjen C. Statistical properties of phase-shift algorithms. J Opt Soc Am A 1995;12(9):1997–2008. [16] Hibino K. Susceptibility of systematic error-compensating algorithms to random noise in phase-shifting interferometry. Appl Opt 1997;36(10):2084–93. [17] Ferrari JA, Flores JL, Ayubi GA, Frins E. Robust interferometer with external phase-shift control. Opt Express 2017;25(24):29965–70.
[1] Malacara D. Optical shop testing. Wiley; 2007. [2] Servin M, Quiroga JA, Padilla JM. Fringe pattern analysis for optical metrology. Wiley-VCH; 2014. [3] Servin M, Estrada JC, Quiroga JA. The general theory of phase shifting algorithms. Opt Express 2009;17:21867–81. [4] Lai G, Yatagai T. Generalized phase-shifting interferometry. J Opt Soc Am A 1991;8:822–7. [5] Téllez-Quinones A, Malacara-Doblado D. Inhomogeneous phase shifting: an algorithm for nonconstant phase displacements. Appl Opt 2010;49:6224–31. [6] Ayubi GA, Perciante CD, Flores JL, Martino JMD, Ferrari JA. Generation of phase-shifting algorithms with n arbitrarily spaced phase-steps. Appl Opt 2014;53:7168–76. [7] Ayubi GA, Perciante CD, Martino JMD, Flores JL, Ferrari JA. Generalized phase-shifting algorithms: error analysis and minimization of noise propagation. Appl Opt 2016;55(6):1461–9. [8] Ayubi GA, Perciante CD, Martino JMD, Flores JL, Ferrari JA. Generalized phase-shifting algorithms: error analysis and minimization of noise propagation: erratum. Appl Opt 2016;55(28):7763.
135