30
OPTICS LETTERS / Vol. 39, No. 1 / January 1, 2014
Optimum plane selection criteria for single-beam phase retrieval techniques based on the contrast transfer function Konstantinos Falaggis,* Tomasz Kozacki, and Malgorzata Kujawinska Institute of Micromechanics and Photonics, Warsaw University of Technology, 8Sw. A. Boboli St., 02-525 Warsaw, Poland *Corresponding author: [email protected] Received September 25, 2013; revised November 13, 2013; accepted November 14, 2013; posted November 18, 2013 (Doc. ID 198340); published December 17, 2013 In this work, an optimum plane selection methodology is reported that can be applied to a wide range of single-beam phase retrieval techniques, based on the contrast transfer function. It is shown that the optimum measurement distances obtained by this method form a geometric series that maximizes the range of spatial frequencies to be recovered using a minimum number of planes. This allows a noise-robust phase reconstruction that does not rely on regularization techniques, i.e., an extensive search for a regularization parameter is avoided. Measurement systems that employ this optimization criteria give an instant deterministic noise-robust phase reconstruction with higher accuracy, and enable the phase retrieval of the entire object spectrum, including lower frequency components. © 2013 Optical Society of America OCIS codes: (100.5070) Phase retrieval; (110.7440) X-ray imaging; (120.0120) Instrumentation, measurement, and metrology; (120.3940) Metrology; (120.5050) Phase measurement; (150.3045) Industrial optical metrology. http://dx.doi.org/10.1364/OL.39.000030
Single-beam phase retrieval techniques (PRTs) have gained increased interest [1–5] because they allow wavefront reconstructions in single beam devices, such as optical microscopes. Iterative PRTs (IPRTs), e.g., the Yan–Gu algorithm [6], allow phase reconstructions from two intensities, but make assumptions on the properties of the object. Variants of single-beam multiple-intensity reconstruction (SBMIR) techniques [5,7] do not require such assumptions, but the convergence of the solver requires a high diversity [5] of measured intensities. For the Fresnel regime, there are two widely known deterministic PRTs available. First, techniques based on the transport of intensity equation (TIE) [1,2] (valid in the limit of small object-to-detector distances) and, second, those based on the contrast transfer function (CTF) [3,4] (valid for objects with slowly varying phase and weak absorption). Reference [3] reported a mixed TIE–CTF approach [3] that combines the advantages of both solvers. TIE-based techniques relate the phase of a wave to the axial intensity derivative, where the latter is estimated using a series of intensity measurements. However, errors in the estimation of the axial derivative directly propagate into the calculated phase. A suitable plane selection for noise-robust measurements depends on prior knowledge of the object, i.e., in the form of the second [1] or third [2] axial intensity derivative. On the contrary, CTF-based solvers and mixed approaches [3] do not rely on such prior knowledge, because the limitation lies in the quest overcoming the roots and minima of the CTF when applying the inversion. Commonly, the inversion either involves a limitation of the spatial frequencies via low-pass filtering, or the use of regularization techniques (RT) that limit the ability to measure low spatial frequency components. Nevertheless, both options affect the accuracy of the reconstruction, as they do not prevent the presence of these minima. In this Letter, an analysis of the mixed TIE–CTF approach is presented, with a focus on the distribution of recovered spatial frequency content. It is shown that 0146-9592/14/010030-04$15.00/0
low values, or zeros, of the CTF create difficulty in the CTF inversion. However, this can be overcome for a wide range of spatial frequencies when using an optimum plane selection. The approach presented here is based on a geometric series [8], and gives a robust reconstruction, using a minimum number of planes. In this way, noise is reduced and the need of RT is overcome, such that the entire object spectrum, including lower frequencies, can be recovered. Let us consider the case of a coherent wave illuminating an object with a 3D complex refractive index distribution, where the wave–object interaction can be described by the transmittance function Tx; y expBx; y iϕx; y, where x and y are the spatial coordinates in the transverse plane and Bx; y and ϕx; y are the absorption and phase shift induced by the object, respectively. Both Bx; y and ϕx; y are given as: Z Bx; y 2π∕λ
βx; y; zdz; Z φx; y −2π∕λ δx; y; zdz;
where βx; y; z and δx; y; z are related to the refractive index distribution nx; y; z as nx; y; z 1 − δx; y; z iβx; y; z. The mixed approach [3] recovers the phase using following semi-iterative procedure: FfI 0
ϕv1 xg
PN
n1
aDn ; f I~ Dn f − bv Dn ; f ; 2W f ε (1)
where I~ D f is the Fourier-transform of the captured intensity, N is the number of planes, Dn is the defocus distance at the nth plane, v is the iteration step, and φv x is the phase at the iteration step v with φ0 x 0 © 2014 Optical Society of America
January 1, 2014 / Vol. 39, No. 1 / OPTICS LETTERS
W f
N X
31
wn Dn ; f ; wn Dn ; f aD; f 2 ;
n1
aD; f sinπλDf 2 ; bv D; f I~ D f ϕ0 cosπλDf 2
λD Ffφv x∇I 0 g: 2π
(2)
The problem of the phase retrieval in Eq. (1) is the division by the factor 2W f . Spatial frequencies having a small value of W f give rise to strong noise amplification. The simplest approach to avoid this problem is to use filtering techniques, i.e., low-pass filters, or, alternatively, applying a regularization parameter (RP) ε > 0, so that 1∕2W ε does not exceed a certain limit. The function W f is periodic because it comprises a total of N periodic functions wn . In the spatial discrete case, the spatial frequencies are given as f x mx Δf and f y my Δf where Δf M⋅p−1 , mx and my are integers, p is the pixel pitch, and M is the number of pixels of the measured intensities in the x− or y− direction. Hence, wn can be expressed as wn sinπλDn mx Δf 2 my Δf 2 2 sinπq∕T Dn 2 , where T Dn is the period of wn with T Dn 1∕λDn Δf 2 , and q is a positive integer with q m2x m2y . Equation (2) shows that the first minimum of W f occurs when all individual terms wn are close to zero, i.e., all individual products q⋅T Dn coincide and are close to an integer. A similar problem is encountered in multiwavelength interferometry (MWI) for phase unwrapping problems [8], where, e.g., a series of synthetic wavelengths with various periods is created. The quest is to choose wavelengths that maximize the measurement range and ensure a given reliability [8]. It is possible to apply the same methodology to the terms wn , so it is ensured that the range of spatial frequencies of W f is above a certain threshold. One strategy that can be adopted is based on the general optimum multifrequency (GOMF) approach [8], which gives an optimal solution for Dn . GOMF strategies are often referred to as hierarchical approaches, because they form a hierarchy of beats with a constant geometric progression factor. Consider the example in Fig. 1, where the T Dn are chosen such that: T Dn gr N−n ;
(3)
where r is a geometric progression factor calculated as r DN ∕D1 1∕N−1 , and g is determined by the largest defocus distance DN as: g 1∕λDN Δf 2 :
(4)
The parameters r and g allow working with normalized quantities. The geometric series in Eq. (4) creates a set of periodic patterns of wn , which ensures that W mx; my ≥ sin π∕g2 . The principle of this methodology can be understood as follows: the periodic pattern with the smallest period wN is ≥ sin π∕g2 for all spatial frequencies, except for values of q near integer multiples of g, i.e., q NINTk⋅g, k 0; 1; 2; 3…. The particular choice of T Dn in Eq. (4) ensures that the local minima of wN coincide with values of wN−1 that are ≥ sin π∕g2 . This can be applied successfully throughout
Fig. 1. (a) Value of W f and wn for the hierarchical approach with T Dn 10n , and (b) the comparison of W f for N 4 for the case of the hierarchical approach (100, 794, 6300, and 50 mm, r 166.67, and g 530), and the case of equidistant planes (12.5, 25, 37.5, and 50 mm).
the whole frequency interval, but does not account for values of q at q0 NINTk⋅g⋅r. Similarly, additional terms of wn increase the range of spatial frequencies where ≥ sin π∕g2 , i.e., the addition of a third, fourth, or fifth measurement plane can account for the local minima of wN−2 , wN−3 , and wN−4 , respectively. This can be done continuously until one reaches the measurement plane closest to the object. The range of spatial frequencies that can be covered by this approach is equal to q < NINTgr N−1 − r N−2 , where the maximum values of q are given by q≤M 2 ∕2. Figure 1 shows an example of an “hierarchical structured” W f composed of three wi with T D1 1000, T D2 100, and T D3 10. The term w3 ensures that W f is larger than the threshold within the interval q 1; 9. The sum of the terms w3 and w2 ensures that W f is above the threshold for the interval q 1; 99, but zero for values of q that are a multiples of T D2 , i.e., q 100, 200, 300, etc. These minima are covered by the term w1 , which ensures that W f is above the threshold within the interval q 1; 899. Local minima of W f are found near small integer multiples of T D1 and T D2 and are in the order of the threshold. These considerations have a number of properties: (1) The range of spatial frequencies covered by the hierarchical approach grows exponentially with N. (2) The minimum value of W is equal to sin π∕g2 , where the actual value of g depends on the largest defocus distance and the field of view (FoV). (3) A smaller value of r provides robustness against noise, but requires more measurement planes. (4) The values of r and g can be reduced by increasing the defocus distance or reducing the FoV. (5) The value of the RP depends on g, as it is related to the minima of W f . If sin π∕g2 is within the computational accuracy, then the use of RT is not required.
32
OPTICS LETTERS / Vol. 39, No. 1 / January 1, 2014
(6) The amplitudes of low frequency signal components are damped by the RT (similar to a high-pass filter); hence, a small value of ε is desired. (7) The remaining minimum of W f at q → 0 reflects the difficulty of this type of solver to recover very low spatial frequencies. This can be understood as an illconditioned nature of the phase retrieval problem. However, it is encountered as well in other methods. TIE-based approaches [1,2] have strong noise amplification at f → 0 when applying the inverse Laplacian f 2x f 2y −1 , and the IPRTs struggle to remove low frequency artifacts, due to stagnation [5]. The geometric series in Eq. (4) provides an optimal solution for Dn for N measurement planes, where W f ≥ sin π∕g2 . If there is no preference for the measurement plane D1 (plane closest to the object), then one can choose r g and a GOMF series [8] is formed, solely determined by the largest defocus distance and covering the largest range of spatial frequencies with a minimum number of planes. The example in Fig. 1(b) shows the case of N 4 for the hierarchical approach with Dn at 100, 794, 6300, and 50 mm, (r 166.67, g 530). The choice of r < g provides additional noise robustness. In contrast, the case of equidistant measurement planes (at 12.5, 25, 37.5, and 50 mm) has a significantly larger number of low values of W f . Hence, a large number of additional planes is needed to provide additional support for spatial frequencies with low W f . Alternatively, it is possible to apply RT, using knowledge of the object. A comparison of the performance of this hierarchical plane selection strategy is shown in Fig. 2 for a detector with M 512, p 8 μm, and a signal-to-noise ratio (SNR) of 10 dB. The values of Dn are equal to the case of Fig. 1(b). Figure 2(a) shows the simulation results for a U. S. Air Force (USAF) 1951 resolution test chart that is used as a pure phase object with a phase difference of
π∕2. A subregion of this USAF chart is shown in Fig. 2(b), as well as the reconstruction results for the hierarchical approach and the case of equidistant planes in Figs. 2(c) and 2(d), respectively. Low-frequency artifacts, due to noise amplification near f → 0, are suppressed, using a high-pass filter with a corner frequency of 15Δf . For better comparison and visualization, the same high-pass filter is applied to the reference phase in Fig. 2. The results indicate that the hierarchical approach with no regularization gives a qualitatively good result in the presence of high noise (SNR 10 dB). The case of equidistant measurement planes, with a RP of ε 0.1, provides a visually good result; however, the search for an appropriate RP ε required an additional computational effort. A similar behavior is also obtained for the speckle pattern with an object feature size of 32 μm, a maximum phase difference of π∕2, and a SNR of 15 dB, shown in Fig. 2(e). Similar to the case of the USAF chart, the plane selection criteria, based on the hierarchical approach in Fig. 2(g), provides an instantly good result with ε 0; whereas, in the case of equidistant planes in Fig. 2(h), an extensive search for an appropriate RP is required to achieve a visually good result. The computation time for experimental data on a computer with an Intel Core i7-2600 Processor with 3.4 GHz and an ASRock p67 extreme4 gen3 motherboard with 8 GB RAM using MATLAB was ∼60 s for M 2048 and N 11. The search for an optimum ε increases the computation time, i.e., a number of 100 trials for ε results in a computation time of 6000 s. The chosen value of ε maximizes the IQI [9], which measures losses due to correlation, luminance distortion, and contrast distortion. However, this search for ε requires the ideal reference phase to be known, which cannot be applied in experiments. Figure 3 shows a quantitative comparison of those strategies by means of the image quality index (IQI)
Fig. 2. (a)–(d) Phase reconstructions [rad] of a USAF chart (SNR 10 dB) and, (e)–(h) a speckle patterns with an object feature size of 32 μm, M 512, p 8 μm (SNR 15 dB). Figures (b)–(d) and (f)–(h) correspond to subregions of the ideal phase object in (a) and (e), respectively. Figures (b) and (f) correspond to the ideal phase object; (c) and (g) correspond to the hierarchical approach; and (d) and (h) correspond to the case of equidistant planes. The values of Dn are equal to the case of Fig. 1(b).
January 1, 2014 / Vol. 39, No. 1 / OPTICS LETTERS
Fig. 3. (a) IQI [9] and (b) RMS for the example in Fig. 2(e) for various levels of SNR.
and root mean square (RMS) for the object in Fig. 2(e). The results indicate that the hierarchical approach improves accuracy. In this work, a new criterion has been introduced for optimum plane selection in single-beam PRTs, based on the TIE–CTF approach. The presented plane selection criterion generates a geometric series of T Dn for wi , which has the ability to cover all spatial frequencies with a minimum number of planes r g. In addition, the choice of r < g provides further noise robustness, at the expense of additional planes. The simulations carried out for the case of N 4 indicate that the set of Dn of the hierarchical approach provides improved accuracy compared with the case of equidistant planes, despite the fact that, for the hierarchical approach, no RT has been applied. Hence, the knowledge of an optimum RP [3] is not required, which is advantageous in the measurement of unknown objects. Low frequency artifacts are a result of noise amplification due to the minimum of W f at q → 0, and are related to the largest defocus distance that influences the parameter g. Figures 4(a)–4(c) show examples where the choice of DN can effectively suppress low-frequency artifacts in the reconstruction, giving this method advantages compared with TIE approaches [1,2], which are not immune to such errors for the case of no regularization. Figures 4(d)–4(f) show the same case, but highlight one major limitation of RT that inevitably introduces low frequency phase errors. The use of ε in Eq. (1) damps the low frequency signal components because 1∕2W ε cannot exceed a certain limit; this is similar to the case of a nonlinear high-pass filter. Hence, the measurements of slowly varying object features cannot be accomplished when using RT, as they require a very small value of ε or a hierarchical strategy with ε 0.
33
Fig. 4. (a) Reconstructed phase of a random object using the GOMF approach (feature size 128 μm, phase difference π∕2, M 512, p 8 μm), for SNR 25 dB, N 11, and DN 5 mm, and (b) DN 50 mm. The cross sections (c) show the ideal phase object (green curve) and the cases of DN 5 mm (blue) and DN 50 mm (red), respectively. Figures (d)–(f) correspond to the cases of (a)–(c), but using ε 0.1. The phases in (a)–(f) were not high-pass filtered and include all low-frequency artifacts.
In summary, this Letter presents an optimization scheme for the mixed TIE–CTF approach [3] that is based on a geometric series [8]. This plane selection method improves the measurement accuracy with a minimum number of planes, increases the immunity to noise, does not require RT, and recovers the entire object spectrum, including low frequency components. The research leading to the described results was realized within program TEAM/2011-7/7 of the Foundation for Polish Science, cofinanced by the European Funds of Regional Development. We acknowledge the support of the statutory funds of Warsaw University of Technology. References 1. M. Soto and E. Acosta, Appl. Opt. 46, 7978 (2007). 2. D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, J. Microsc. 214, 51 (2004). 3. J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, Opt. Lett. 32, 1617 (2007). 4. M. Langer, P. Cloetens, and F. Peyrin, IEEE Trans. Image Process. 19, 2428 (2010). 5. K. Falaggis, T. Kozacki, and M. Kujawinska, Opt. Lett. 38, 1660 (2013). 6. G. Z. Yang, B. Z. Dong, B. Y. Gu, J. Y. Zhuang, and O. K. Ersoy, Appl. Opt. 33, 209 (1994). 7. G. Pedrini, W. Osten, and Y. Zhang, Opt. Lett. 30, 833 (2005). 8. C. E. Towers, D. P. Towers, and J. D. C. Jones, Opt. Lett. 29, 1348 (2004). 9. A. C. Bovik, IEEE Signal Process. Lett. 9, 81 (2002).
OPTICS LETTERS / Vol. 39, No. 1 / January 1, 2014
Optimum plane selection criteria for single-beam phase retrieval techniques based on the contrast transfer function Konstantinos Falaggis,* Tomasz Kozacki, and Malgorzata Kujawinska Institute of Micromechanics and Photonics, Warsaw University of Technology, 8Sw. A. Boboli St., 02-525 Warsaw, Poland *Corresponding author: [email protected] Received September 25, 2013; revised November 13, 2013; accepted November 14, 2013; posted November 18, 2013 (Doc. ID 198340); published December 17, 2013 In this work, an optimum plane selection methodology is reported that can be applied to a wide range of single-beam phase retrieval techniques, based on the contrast transfer function. It is shown that the optimum measurement distances obtained by this method form a geometric series that maximizes the range of spatial frequencies to be recovered using a minimum number of planes. This allows a noise-robust phase reconstruction that does not rely on regularization techniques, i.e., an extensive search for a regularization parameter is avoided. Measurement systems that employ this optimization criteria give an instant deterministic noise-robust phase reconstruction with higher accuracy, and enable the phase retrieval of the entire object spectrum, including lower frequency components. © 2013 Optical Society of America OCIS codes: (100.5070) Phase retrieval; (110.7440) X-ray imaging; (120.0120) Instrumentation, measurement, and metrology; (120.3940) Metrology; (120.5050) Phase measurement; (150.3045) Industrial optical metrology. http://dx.doi.org/10.1364/OL.39.000030
Single-beam phase retrieval techniques (PRTs) have gained increased interest [1–5] because they allow wavefront reconstructions in single beam devices, such as optical microscopes. Iterative PRTs (IPRTs), e.g., the Yan–Gu algorithm [6], allow phase reconstructions from two intensities, but make assumptions on the properties of the object. Variants of single-beam multiple-intensity reconstruction (SBMIR) techniques [5,7] do not require such assumptions, but the convergence of the solver requires a high diversity [5] of measured intensities. For the Fresnel regime, there are two widely known deterministic PRTs available. First, techniques based on the transport of intensity equation (TIE) [1,2] (valid in the limit of small object-to-detector distances) and, second, those based on the contrast transfer function (CTF) [3,4] (valid for objects with slowly varying phase and weak absorption). Reference [3] reported a mixed TIE–CTF approach [3] that combines the advantages of both solvers. TIE-based techniques relate the phase of a wave to the axial intensity derivative, where the latter is estimated using a series of intensity measurements. However, errors in the estimation of the axial derivative directly propagate into the calculated phase. A suitable plane selection for noise-robust measurements depends on prior knowledge of the object, i.e., in the form of the second [1] or third [2] axial intensity derivative. On the contrary, CTF-based solvers and mixed approaches [3] do not rely on such prior knowledge, because the limitation lies in the quest overcoming the roots and minima of the CTF when applying the inversion. Commonly, the inversion either involves a limitation of the spatial frequencies via low-pass filtering, or the use of regularization techniques (RT) that limit the ability to measure low spatial frequency components. Nevertheless, both options affect the accuracy of the reconstruction, as they do not prevent the presence of these minima. In this Letter, an analysis of the mixed TIE–CTF approach is presented, with a focus on the distribution of recovered spatial frequency content. It is shown that 0146-9592/14/010030-04$15.00/0
low values, or zeros, of the CTF create difficulty in the CTF inversion. However, this can be overcome for a wide range of spatial frequencies when using an optimum plane selection. The approach presented here is based on a geometric series [8], and gives a robust reconstruction, using a minimum number of planes. In this way, noise is reduced and the need of RT is overcome, such that the entire object spectrum, including lower frequencies, can be recovered. Let us consider the case of a coherent wave illuminating an object with a 3D complex refractive index distribution, where the wave–object interaction can be described by the transmittance function Tx; y expBx; y iϕx; y, where x and y are the spatial coordinates in the transverse plane and Bx; y and ϕx; y are the absorption and phase shift induced by the object, respectively. Both Bx; y and ϕx; y are given as: Z Bx; y 2π∕λ
βx; y; zdz; Z φx; y −2π∕λ δx; y; zdz;
where βx; y; z and δx; y; z are related to the refractive index distribution nx; y; z as nx; y; z 1 − δx; y; z iβx; y; z. The mixed approach [3] recovers the phase using following semi-iterative procedure: FfI 0
ϕv1 xg
PN
n1
aDn ; f I~ Dn f − bv Dn ; f ; 2W f ε (1)
where I~ D f is the Fourier-transform of the captured intensity, N is the number of planes, Dn is the defocus distance at the nth plane, v is the iteration step, and φv x is the phase at the iteration step v with φ0 x 0 © 2014 Optical Society of America
January 1, 2014 / Vol. 39, No. 1 / OPTICS LETTERS
W f
N X
31
wn Dn ; f ; wn Dn ; f aD; f 2 ;
n1
aD; f sinπλDf 2 ; bv D; f I~ D f ϕ0 cosπλDf 2
λD Ffφv x∇I 0 g: 2π
(2)
The problem of the phase retrieval in Eq. (1) is the division by the factor 2W f . Spatial frequencies having a small value of W f give rise to strong noise amplification. The simplest approach to avoid this problem is to use filtering techniques, i.e., low-pass filters, or, alternatively, applying a regularization parameter (RP) ε > 0, so that 1∕2W ε does not exceed a certain limit. The function W f is periodic because it comprises a total of N periodic functions wn . In the spatial discrete case, the spatial frequencies are given as f x mx Δf and f y my Δf where Δf M⋅p−1 , mx and my are integers, p is the pixel pitch, and M is the number of pixels of the measured intensities in the x− or y− direction. Hence, wn can be expressed as wn sinπλDn mx Δf 2 my Δf 2 2 sinπq∕T Dn 2 , where T Dn is the period of wn with T Dn 1∕λDn Δf 2 , and q is a positive integer with q m2x m2y . Equation (2) shows that the first minimum of W f occurs when all individual terms wn are close to zero, i.e., all individual products q⋅T Dn coincide and are close to an integer. A similar problem is encountered in multiwavelength interferometry (MWI) for phase unwrapping problems [8], where, e.g., a series of synthetic wavelengths with various periods is created. The quest is to choose wavelengths that maximize the measurement range and ensure a given reliability [8]. It is possible to apply the same methodology to the terms wn , so it is ensured that the range of spatial frequencies of W f is above a certain threshold. One strategy that can be adopted is based on the general optimum multifrequency (GOMF) approach [8], which gives an optimal solution for Dn . GOMF strategies are often referred to as hierarchical approaches, because they form a hierarchy of beats with a constant geometric progression factor. Consider the example in Fig. 1, where the T Dn are chosen such that: T Dn gr N−n ;
(3)
where r is a geometric progression factor calculated as r DN ∕D1 1∕N−1 , and g is determined by the largest defocus distance DN as: g 1∕λDN Δf 2 :
(4)
The parameters r and g allow working with normalized quantities. The geometric series in Eq. (4) creates a set of periodic patterns of wn , which ensures that W mx; my ≥ sin π∕g2 . The principle of this methodology can be understood as follows: the periodic pattern with the smallest period wN is ≥ sin π∕g2 for all spatial frequencies, except for values of q near integer multiples of g, i.e., q NINTk⋅g, k 0; 1; 2; 3…. The particular choice of T Dn in Eq. (4) ensures that the local minima of wN coincide with values of wN−1 that are ≥ sin π∕g2 . This can be applied successfully throughout
Fig. 1. (a) Value of W f and wn for the hierarchical approach with T Dn 10n , and (b) the comparison of W f for N 4 for the case of the hierarchical approach (100, 794, 6300, and 50 mm, r 166.67, and g 530), and the case of equidistant planes (12.5, 25, 37.5, and 50 mm).
the whole frequency interval, but does not account for values of q at q0 NINTk⋅g⋅r. Similarly, additional terms of wn increase the range of spatial frequencies where ≥ sin π∕g2 , i.e., the addition of a third, fourth, or fifth measurement plane can account for the local minima of wN−2 , wN−3 , and wN−4 , respectively. This can be done continuously until one reaches the measurement plane closest to the object. The range of spatial frequencies that can be covered by this approach is equal to q < NINTgr N−1 − r N−2 , where the maximum values of q are given by q≤M 2 ∕2. Figure 1 shows an example of an “hierarchical structured” W f composed of three wi with T D1 1000, T D2 100, and T D3 10. The term w3 ensures that W f is larger than the threshold within the interval q 1; 9. The sum of the terms w3 and w2 ensures that W f is above the threshold for the interval q 1; 99, but zero for values of q that are a multiples of T D2 , i.e., q 100, 200, 300, etc. These minima are covered by the term w1 , which ensures that W f is above the threshold within the interval q 1; 899. Local minima of W f are found near small integer multiples of T D1 and T D2 and are in the order of the threshold. These considerations have a number of properties: (1) The range of spatial frequencies covered by the hierarchical approach grows exponentially with N. (2) The minimum value of W is equal to sin π∕g2 , where the actual value of g depends on the largest defocus distance and the field of view (FoV). (3) A smaller value of r provides robustness against noise, but requires more measurement planes. (4) The values of r and g can be reduced by increasing the defocus distance or reducing the FoV. (5) The value of the RP depends on g, as it is related to the minima of W f . If sin π∕g2 is within the computational accuracy, then the use of RT is not required.
32
OPTICS LETTERS / Vol. 39, No. 1 / January 1, 2014
(6) The amplitudes of low frequency signal components are damped by the RT (similar to a high-pass filter); hence, a small value of ε is desired. (7) The remaining minimum of W f at q → 0 reflects the difficulty of this type of solver to recover very low spatial frequencies. This can be understood as an illconditioned nature of the phase retrieval problem. However, it is encountered as well in other methods. TIE-based approaches [1,2] have strong noise amplification at f → 0 when applying the inverse Laplacian f 2x f 2y −1 , and the IPRTs struggle to remove low frequency artifacts, due to stagnation [5]. The geometric series in Eq. (4) provides an optimal solution for Dn for N measurement planes, where W f ≥ sin π∕g2 . If there is no preference for the measurement plane D1 (plane closest to the object), then one can choose r g and a GOMF series [8] is formed, solely determined by the largest defocus distance and covering the largest range of spatial frequencies with a minimum number of planes. The example in Fig. 1(b) shows the case of N 4 for the hierarchical approach with Dn at 100, 794, 6300, and 50 mm, (r 166.67, g 530). The choice of r < g provides additional noise robustness. In contrast, the case of equidistant measurement planes (at 12.5, 25, 37.5, and 50 mm) has a significantly larger number of low values of W f . Hence, a large number of additional planes is needed to provide additional support for spatial frequencies with low W f . Alternatively, it is possible to apply RT, using knowledge of the object. A comparison of the performance of this hierarchical plane selection strategy is shown in Fig. 2 for a detector with M 512, p 8 μm, and a signal-to-noise ratio (SNR) of 10 dB. The values of Dn are equal to the case of Fig. 1(b). Figure 2(a) shows the simulation results for a U. S. Air Force (USAF) 1951 resolution test chart that is used as a pure phase object with a phase difference of
π∕2. A subregion of this USAF chart is shown in Fig. 2(b), as well as the reconstruction results for the hierarchical approach and the case of equidistant planes in Figs. 2(c) and 2(d), respectively. Low-frequency artifacts, due to noise amplification near f → 0, are suppressed, using a high-pass filter with a corner frequency of 15Δf . For better comparison and visualization, the same high-pass filter is applied to the reference phase in Fig. 2. The results indicate that the hierarchical approach with no regularization gives a qualitatively good result in the presence of high noise (SNR 10 dB). The case of equidistant measurement planes, with a RP of ε 0.1, provides a visually good result; however, the search for an appropriate RP ε required an additional computational effort. A similar behavior is also obtained for the speckle pattern with an object feature size of 32 μm, a maximum phase difference of π∕2, and a SNR of 15 dB, shown in Fig. 2(e). Similar to the case of the USAF chart, the plane selection criteria, based on the hierarchical approach in Fig. 2(g), provides an instantly good result with ε 0; whereas, in the case of equidistant planes in Fig. 2(h), an extensive search for an appropriate RP is required to achieve a visually good result. The computation time for experimental data on a computer with an Intel Core i7-2600 Processor with 3.4 GHz and an ASRock p67 extreme4 gen3 motherboard with 8 GB RAM using MATLAB was ∼60 s for M 2048 and N 11. The search for an optimum ε increases the computation time, i.e., a number of 100 trials for ε results in a computation time of 6000 s. The chosen value of ε maximizes the IQI [9], which measures losses due to correlation, luminance distortion, and contrast distortion. However, this search for ε requires the ideal reference phase to be known, which cannot be applied in experiments. Figure 3 shows a quantitative comparison of those strategies by means of the image quality index (IQI)
Fig. 2. (a)–(d) Phase reconstructions [rad] of a USAF chart (SNR 10 dB) and, (e)–(h) a speckle patterns with an object feature size of 32 μm, M 512, p 8 μm (SNR 15 dB). Figures (b)–(d) and (f)–(h) correspond to subregions of the ideal phase object in (a) and (e), respectively. Figures (b) and (f) correspond to the ideal phase object; (c) and (g) correspond to the hierarchical approach; and (d) and (h) correspond to the case of equidistant planes. The values of Dn are equal to the case of Fig. 1(b).
January 1, 2014 / Vol. 39, No. 1 / OPTICS LETTERS
Fig. 3. (a) IQI [9] and (b) RMS for the example in Fig. 2(e) for various levels of SNR.
and root mean square (RMS) for the object in Fig. 2(e). The results indicate that the hierarchical approach improves accuracy. In this work, a new criterion has been introduced for optimum plane selection in single-beam PRTs, based on the TIE–CTF approach. The presented plane selection criterion generates a geometric series of T Dn for wi , which has the ability to cover all spatial frequencies with a minimum number of planes r g. In addition, the choice of r < g provides further noise robustness, at the expense of additional planes. The simulations carried out for the case of N 4 indicate that the set of Dn of the hierarchical approach provides improved accuracy compared with the case of equidistant planes, despite the fact that, for the hierarchical approach, no RT has been applied. Hence, the knowledge of an optimum RP [3] is not required, which is advantageous in the measurement of unknown objects. Low frequency artifacts are a result of noise amplification due to the minimum of W f at q → 0, and are related to the largest defocus distance that influences the parameter g. Figures 4(a)–4(c) show examples where the choice of DN can effectively suppress low-frequency artifacts in the reconstruction, giving this method advantages compared with TIE approaches [1,2], which are not immune to such errors for the case of no regularization. Figures 4(d)–4(f) show the same case, but highlight one major limitation of RT that inevitably introduces low frequency phase errors. The use of ε in Eq. (1) damps the low frequency signal components because 1∕2W ε cannot exceed a certain limit; this is similar to the case of a nonlinear high-pass filter. Hence, the measurements of slowly varying object features cannot be accomplished when using RT, as they require a very small value of ε or a hierarchical strategy with ε 0.
33
Fig. 4. (a) Reconstructed phase of a random object using the GOMF approach (feature size 128 μm, phase difference π∕2, M 512, p 8 μm), for SNR 25 dB, N 11, and DN 5 mm, and (b) DN 50 mm. The cross sections (c) show the ideal phase object (green curve) and the cases of DN 5 mm (blue) and DN 50 mm (red), respectively. Figures (d)–(f) correspond to the cases of (a)–(c), but using ε 0.1. The phases in (a)–(f) were not high-pass filtered and include all low-frequency artifacts.
In summary, this Letter presents an optimization scheme for the mixed TIE–CTF approach [3] that is based on a geometric series [8]. This plane selection method improves the measurement accuracy with a minimum number of planes, increases the immunity to noise, does not require RT, and recovers the entire object spectrum, including low frequency components. The research leading to the described results was realized within program TEAM/2011-7/7 of the Foundation for Polish Science, cofinanced by the European Funds of Regional Development. We acknowledge the support of the statutory funds of Warsaw University of Technology. References 1. M. Soto and E. Acosta, Appl. Opt. 46, 7978 (2007). 2. D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, J. Microsc. 214, 51 (2004). 3. J. P. Guigay, M. Langer, R. Boistel, and P. Cloetens, Opt. Lett. 32, 1617 (2007). 4. M. Langer, P. Cloetens, and F. Peyrin, IEEE Trans. Image Process. 19, 2428 (2010). 5. K. Falaggis, T. Kozacki, and M. Kujawinska, Opt. Lett. 38, 1660 (2013). 6. G. Z. Yang, B. Z. Dong, B. Y. Gu, J. Y. Zhuang, and O. K. Ersoy, Appl. Opt. 33, 209 (1994). 7. G. Pedrini, W. Osten, and Y. Zhang, Opt. Lett. 30, 833 (2005). 8. C. E. Towers, D. P. Towers, and J. D. C. Jones, Opt. Lett. 29, 1348 (2004). 9. A. C. Bovik, IEEE Signal Process. Lett. 9, 81 (2002).
