Through-focus phase retrieval and its connection to the spatial correlation for propagating fields

Through-focus phase retrieval methods aim to retrieve the phase of an optical field from its intensity distribution measured at different planes in the focal region. By using the concept of spatial correlation for propagating fields, for both the complex amplitude and the intensity of a field, we can infer which planes are suitable to retrieve the phase and which are not. Our analysis also reveals why all techniques based on measuring the intensity at two Fourier-conjugated planes usually lead to a good reconstruction of the phase. The findings presented in this work are important for aberration characterization of optical systems, adaptive optics and wavefront metrology. © 2013 Optical Society of America OCIS codes: (350.5500) Propagation; (100.5070) Phase retrieval; (010.7350) Wavefront sensing. References and links 1. R. Loudon, The Quantum Theory of Light (Oxford University Press, 2000). 2. R. W. Gerchberg, and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972). 3. J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758-2769 (1982). 4. H. M. Faulkner, and J. M. Rodenburg, “Movable aperture lensless transmission microscopy: a novel phase retrieval algorithm,” Phys. Rev. Lett. 93, 023903 (2004). 5. H. M. Faulkner, and J. M. Rodenburg, “Error tolerance of an iterative phase retrieval algorithm for moveable illumination microscopy,” Ultramicroscopy 103 153–164 (2005). 6. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400 342–344 (1999). 7. J. A. Rodrigo, T. Alieva, G. Cristbal, and M. L. Calvo, “Wavefield imaging via iterative retrieval based on phase modulation diversity,” Opt. Express 19 18621–18635 (2011). 8. M. R. Teague, “Image formation in terms of the transport equation,” J. Opt. Soc. Am. A 2, 2019–2026 (1985). 9. G. R. Brady, M. Guizar-Sicairos and J. R. Fienup, “Optical wavefront measurement using phase retrieval with transverse translation diversity,” Opt. Express 17 624−639 (2009). 10. R. A. Gonsalves, “Phase retrieval and diversity in adaptive optics,” Opt. Eng. 21, 829–832 (1982). 11. J. Miao, D. Sayre, and H. N. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). 12. B. Hanser, M. Gustafsson, D. Agard, and J. Sedat, “Phase retrieval for high-numerical-aperture optical systems,” Opt. Lett. 28, 801–803 (2003). #178863 $15.00 USD Received 1 Nov 2012; revised 12 Dec 2012; accepted 13 Dec 2012; published 27 Feb 2013 (C) 2013 OSA 11 March 2013 / Vol. 21, No. 5 / OPTICS EXPRESS 5550 13. J. Durnin, J. J. Miceli, and J. H. Eberly,“Diffraction-free beams,” Phys. Rev. Lett. 58, 1499–1501 (1987). 14. M. Born, and E. Wolf, Principles of Optics (Cambridge Univ. Press, 2007), 883–891.


Introduction
Measuring the phase of an optical field is a challenging issue.In fact, even the fastest detectors available nowadays have integration times that are several orders of magnitude larger than the temporal period of light oscillations [1].It follows that only the intensity of the field is directly accessible while any phase information is lost.This is at the basis of the well known phase problem in optics that is shared by many other disciplines as well, such as x-ray, electron or neutron diffraction.In optics, retrieving the phase of a field plays also an important role in microscopy or in the characterization of aberrations of an optical system.Although there are methods that allow to turn the phase information into an intensity one, such as interferometry, there exist specific applications that require solutions easier from an implementation point of view.For this reason, much attention has been paid to find alternative approaches to retrieve the phase of a wavefield, among which the most successful ones are probably the so-called iterative methods.In these methods, a set of intensity patterns is measured and used to retrieve the phase of the field everywhere in space by exploiting the constraints imposed by the intensity on the phase of a field.This stems from the fact that amplitude and phase of an optical field are not independent quantities but they are actually related via Maxwell's equations.Since the pivotal work by Gerchberg and Saxton of 1972 [2], several variations of the basic phase retrieval algorithm have been introduced by different authors [3][4][5][6][7][8][9][10].Among all the interesting works, we would like to recall the important contribution given by Miao et al. [11] to the understanding of the physical mechanisms behind a phase retrieval problem.These authors have suggested that when a detector endowed with N pixels is used to measure the intensity of a field at some plane, each pixel represents one equation of a (generally nonlinear) system of equations where 2N unknowns are present, the amplitude and phase of the field.This leads to the conclusion that the phase problem is undetermined by a factor of two and that at least two different intensity patterns must be used in a reconstruction algorithm.In most works appeared so far, these two planes are chosen to be Fourier-conjugated planes such as the near and far-field distributions of a wavefield.However, such minimum number of intensity measurements can be affected by two additional ingredients.The first one is the unavoidable presence of noise, which makes necessary to gather additional experimental data then the minimum required in order to successfully reconstruct the phase of a field.The second element is the availability of some a priori information, which can compensate the presence of noise in such a way that sometimes even a single intensity pattern can be sufficient to retrieve the phase of a field.An example of a typical a priori information is represented by the knowledge that the field (or its angular spectrum) is bounded to a limited support.In that case, one can set to zero the value of the retrieved complex field outside the boundaries of the said support.In this way, a consistent number of unknowns are removed from the problem.A special class of phase retrieval problems is represented by the through-focus ones.By this term it is intended a family of problems where the phase of a focused field needs to be determined by measuring its intensity in the focal region [12].For some reasons, that we are going to discuss in more detail later on in the paper, these problems present some additional difficulties compared to the classical phase retrieval problems investigated in the past and they require some extra care in their treatment.Here we present a theoretical model, based on the concept of spatial correlation for a propagated field, that allows us to properly predict the planes where the intensity has to be measured in order to successfully retrieve the phase of a focused field.Interestingly, our theoretical findings make also clear why choosing a couple of Fourier-conjugated planes is usually a good choice to set a phase retrieval method.Finally we apply our theoretical analysis to an experiment where we retrieve the phase of an optical field focused by an actual microscope objective of nominal numerical aperture NA = 0.4.The paper is organized as follows.In section 2 the theory is introduced and put into practise in section 3 by showing the experimental results of a phase retrieval problem for a through-focus configuration.Finally, in section 4 we summarize the main results of the work.

Spatial correlation for propagated fields
As anticipated in the introduction, in this section we need to develop some general facts on spatial correlations for propagated fields.The results of our analysis will allow us to define a criterion to choose which intensity measurement to use for phase retrieval in a through-focus configuration.In the next sections, we will apply the results obtained in the present section to the specific case of focused fields.Let us begin by considering a field U(x, y, z) which, for the time being, we will assume to be a scalar field, propagating in a given region of space.At the moment we do not need to specify whether U(x, y, z) is produced by a focusing system or not but we will simply assume that the field takes a given value on a reference plane that from now on we will label as z = 0 plane.Assuming that the field intensity is measured at the plane z 1 , we would like to give a quantitative answer to the following question: is there any additional information that can be extracted from the field by measuring its intensity at another plane z 2 ?In order to answer this question we need to quantify the similarities between two field distributions U(x, y, z 1 ) and U(x, y, z 2 ) originating from the same field distribution on the input plane z = 0, respectively.As common practice in other contexts, such as coherence theory, such similarity can be quantified by defining the degree of correlation between the two field distributions, as the normalized internal product which can be proven to be bounded, i.e., 0 ≤ |c(0, z)| ≤ 1.In Eq. (1) we have chosen, without loss of generality, z 1 = 0 and z 2 = z.When |c(0, z)| = 1 we say that the two fields are completely correlated.This happens when the two distributions U(x, y, z 1 ) and U(x, y, z 2 ) only differ by a trivial scaling factor.For instance, it is easy to check that in case the field belongs to the class of nondiffracting beams (such as a Bessel beam, for instance [13]), |c(0, z)| = 1 for any z.This of course just reflects the fact that the intensity distribution of a non-diffracting beam is invariant on propagation so that no additional information is gained by measuring the field intensity at different planes.
It is interesting to notice that c(0, z) can be expressed in terms of the angular spectrum of the field at z = 0 only as (Plancherel's theorem) where the angular spectrum at z = 0 is defined as and In Eq. (3), p and q represent the spatial frequencies in Fourier space.The integration domain D, that appears in the integrals in Eq. ( 2), is a subset of the homogeneous domain only.This means that we are neglecting any evanescent wave in the spectrum.This assumption is surely applicable when one deals with fields that propagate over distances longer than the wavelength λ .In case of a focused field, the domain D will coincide with the domain {(p, q) |p 2 + q 2 ≤ (NA/λ ) 2 }, with NA the numerical aperture of the optical system.The form of c(0, z) in terms of the angular spectrum of the field U(x, y, z) allows us to make some interesting observations.First of all, we point out that the correlation c(0, z) is a function of z.This indicates that, to what concerns the spatial correlation, different planes might be not all equivalent to each other.Another interesting observation to make is that the correlation c(0, z), which is a complex quantity, depends however only on the absolute value of the angular spectrum on z = 0.This observation plays an important role when the field U(x, y, z) is a focused field.In fact, in that case, the angular spectrum of the field in focus is directly related to the field distribution at the exit pupil of the focusing lens or objective.In most cases, the aberrations of the system affect only the phase of the field in the pupil which, however, does not change the form of c(0, z).Hence, under these circumstances the behavior of c(0, z) can be predicted without knowing the said phase aberrations.This is an appealing property since aberrations are usually unknown and they are typically what one would like to obtain by solving the phase problem.
As to the dependence on z, it is easy to see that This means that the correlation depends only on the distance between the two planes, and not on the absolute position of each single plane.Again, this has interesting consequences in case one deals with a focused field, where the field distributions at symmetrical positions before and after the focal plane, although eventually completely different, correlate exactly in the same way with the field in focus.In order to discuss the utility of c(0, z), let us consider an explicit 1. Degree of spatial correlation c(0, z) between the field in focus and the field in any other plane in the focal region, as function of the dimensionless variable z NA 2 /λ , with NA = 0.3.This curve corresponds to the case of a field with uniform angular spectrum A (0) (p, q) within NA. example that will also help us to explain how the theory can be applied to more general cases.Let us assume we have a field whose angular spectrum is uniform within the numerical aperture of a given optical system, namely A (0) (p, q) = 1, ∀(p, q) ∈ D. This case can well describe the situation of an expanded field with a diameter much larger than the diameter of a focusing objective.Interestingly, one can give the closed-form expression for the correlation c(0, z) in where in the integral in the numerator of Eq. ( 6) a simple integration by parts has been applied (the integral in the denominator is trivial and gives πNA 2 /λ 2 ).In order to solve the integrals, we have also made use of the circular coordinates (ρ, ϕ), with p = ρ cos ϕ and q = ρ sin ϕ.A simple second-order Taylor expansion shows that c(0, 0) = 1, as it should be.In Fig. 1 we plot the absolute value of the degree of spatial correlation c(0, z), as function of the dimensionless variable zNA 2 /λ , for NA = 0.3.One should notice that, at specific positions along the optical axis, the field at z = 0 and that at z tend to become orthogonal or uncorrelated, namely |c(0, z)| 0. Under these circumstances there is not so much common information shared by two fields U(x, y, 0) and U(x, y, z).For a phase retrieval method, this is a key observation since, from what we said in the introduction, one needs to feed into a reconstruction algorithm data sets that are as much independent as possible.A couple of planes showing very low values of spatial correlation c(0, z) are good candidates for this.It is also worth noting that for large values of z the envelope of the degree of spatial correlation goes to zero as 1/z.This can be actually proven by applying the method of stationary phase to Eq. ( 2) [14] which gives, for large z: Considered that for z → ∞ the field distribution tends to the far-field of the field at z = 0, Eq. ( 7) tells us that a Fourier-conjugated couple of planes, namely the near and the far-field, is probably always a good choice for a phase retrieval method, since the degree of correlation is always low in that case.This observation might be the physical explanation of the success of all phase retrieval methods appeared so far that are based on measuring the field intensity at the near and at the far field plane.So far we have considered the degree of spatial correlation for the complex fields but it has been pointed out how only the intensity of a field is actually accessible.For this reason, also the degree of spatial correlation between the corresponding intensity distributions needs to be evaluated.As common use, we define such degree of correlation, denoted by c int (0, z), as Unlike c(0, z), the expression of c int (0, z) cannot be simplified by passing to the Fourier domain.
For the reader's convenience we plot in Fig. 2 the absolute value of c int (0, z) as function of the dimensionless variable zNA 2 /λ , again for the field having A (0) (p, q) = 1, as in Fig. 1.We see that this correlation shows similar features as c(0, z), but the minima are less pronounced.Also, unlike c(0, z), c int (0, z) generally depends on the amplitude and the phase of the field  2. Degree of spatial correlation c int (0, z) between the intensity in focus and that in any other plane in the focal region, as function of the dimensionless variable zNA 2 /λ , with NA = 0.3.This curve corresponds to the case of a field with uniform angular spectrum A (0) (p, q) within NA.
angular spectrum A (0) (p, q).This means that its shape changes if aberrations are present.More specifically, any symmetry is generally broken in presence of aberrations.What we have learnt from this example is that given a field distribution in the exit pupil of a focusing system of interest, it is possible to find planes in the focal region that are poorly correlated with each other.These planes are the best candidates to be used for phase retrieval.Armed now with these theoretical findings, in the next section we will put the theory into practice by retrieving the phase of a focused field generated by a real, aberrated microscope objective of given nominal numerical aperture NA.

Through-focus phase retrieval: experimental implementation and reconstruction
In order to demonstrate the phase retrieval technique introduced in section 2, we built an experimental setup, where through focus intensity planes could be recorded.In this section, we show the experimental details followed by results and reconstruction.

Experimental setup
The schematic of the experimental setup is shown in Fig. 3.The beam from an intensity sta- placed on two translation stages.A plane in the focal region of the lens under test is imaged by the imaging system consisting of the imaging objective, the tube lens and the CCD camera.The aberrations of the test lens influence the intensity distribution in the planes across the focal region.The test lens is then moved along the optical axis by a mechanical translation stage and several approximately equidistant intensity planes in the vicinity of the focal plane are recorded in the camera.The movement along the optical axis is recorded with a displacement interferometer for the precise location of the intensity planes.The test lens in the present setup is a microscope objective with a numerical aperture of NA = 0.4.The imaging objective is of higher numerical aperture (NA = 0.9) in order to capture the complete angular spectrum of the test objective.Intensity planes in the vicinity of the focal plane are recorded in the camera with the same acquisition time.The pixel width of the 8-bit CCD camera is 3.75μm with 1600 × 1200 pixels in x and y directions respectively.The magnification of the measurement system has been calibrated by moving the test lens in the direction perpendicular to the optical axis with a piezo stage.The piezo movement has been calibrated with a displacement interferometer before the measurements were performed.A linear fit to the recorded positions of the focal spots reveals the magnification of the imaging system.In the present configuration of the experimental setup, the magnification is 50.4.The incident wavefront has been characterized through a Shack Hartmann sensor for high degree of collimation with rms λ /100 on a beam size of 8 mm, which corresponds to the diameter of the test objective.The aberrations of the measurement system have been calibrated for precise measurements.More specifically, the aberrations of the imaging objective have been characterized with an interferometer while the aberrations of the tube lens could not be measured with the interferometer due to parasitical interference fringes of the internal surfaces of the tube lens.However, the aberrations of the tube lens and imaging objective resulted to have no significant impact on the final reconstruction on the phase.Dark frames have been recorded before the through-focus measurements by blocking the illumination beam before the lens.These dark frames account for stray light in the setup and have been subtracted from the intensity measurements in the focal region.Each intensity plane is the average of 50 intensity maps at the same location.These dark frames account for stray light in the setup and have been subtracted from the intensity measurements in the focal region.Each intensity plane is the average of 50 intensity maps at the same location.The range of measurements is few microns (≈ 10μm) in both directions with respect to the focal plane.

Phase reconstruction
By means of the experimental setup described before, we have recorded several sets of throughfocus intensity distributions.Here we show the reconstruction in case of the focused field produced by a microscope objective of nominal numerical aperture NA = 0.4.By using the intensity distributions measured at different positions along the optical axis, we have first computed the (measured) degree of spatial correlation for the intensity c int (0, z).In Fig. 4 (panel a) we show c int (0, z) as obtained by using the measured through focus intensity distributions.As it can be easily seen, the curve is not symmetric, which is a manifestation of the presence of aberrations.In the same figure, panel b, we also plot the ideal degree of spatial correlation for both the intensity and for the field, in case of an ideal aberration-free system, uniformly illuminated.One can see that the position of the planes with minimal correlation does not change considerably from the aberrated to the aberration-free case.From the correlation for the intensities, we select the two measured distributions shown in Fig. 4, panel c and d, respectively.The first one corresponds to the field in focus, while the other is located 7.83 μm far from the best focus.The measured correlation between the intensities at those two planes resulted to be c int,meas (0, 7.83 μm) = 0.09, while the ideal one (no aberrations) for the fields and the intensities give |c(0, 7.83 μm)| = 0.03 and c int (0, 7.83 μm) = 0.12, respectively.The two chosen intensity measurements were then used in a standard iterative phase retrieval algorithm, that cycled continuously between the planes z = 0 and z = 7.83 μm.In each cycle of the algorithm, the intensity measured at each plane was used as constraint for the complex field that was propagated, back and forth, between the two planes.After few iterations (typically less than 50 − 60 already gave good results), the phase of the fields at both planes was obtained.In Fig. 5 we show such retrieved phases.In every phase retrieval method it is important to check whether the retrieved phase is correct or not.We have checked the correctness of our results in the following way.We have taken the complex field in focus (z = 0), consisting of the measured amplitude and the retrieved phase (shown in Fig. 5, first column), and we have propagated it along z.If the phase is the right one, then the intensity distributions obtained in this way should coincide with the intensities that have been measured.In Fig. 6 some of the measured and simulated field intensities are plotted where, the top row contains the measurements at different distances from the focal plane and the bottom row contains the simulations.As it is evident from the plots, the agreement between simulations and measurements is excellent, which confirms that the phase has been successfully retrieved.It is worth noting that even when the signal to noise ratio is low (such at z = 15.17 μm distance from the best focus) the calculated intensity shows a strong similarity with the measured one.In order to be more quantitative, we have computed the root mean square error (RMSE) between each simulated and measured intensity pattern shown in (C) 2013 OSA Fig. 6.We have found that at z = 4.04 μm the RMS deviation is 0.68%, at z = 7.83 μm is 0.38%, at z = 11μm is 0.47% and finally at z = 15.17 μm the deviation is 0.39%.

Role of SNR and a priori information
In the reconstructions shown so far, some a priori information has been used.More specifically, the nominal numerical aperture of the microscope objective NA = 0.4 has been included in the iterative algorithm to retrieve the phase of the field.Considered that a microscope objective in air can, at most, have NA = 1, this has allowed us to throw away a lot of data (actually more than half of the field in the exit pupil) that contained only noise but no useful signal.This makes the reconstruction much more robust and less sensitive to eventually high correlations between the two intensity patterns used in the phase retrieval algorithm.This means that, in this specific case, even choosing the field in focus (at the plane z = 0) and another field (at position z) with a high value of the correlation c int (0, z) could actually lead to an excellent reconstruction for the phase.In fact, while the increase in correlation is compensated by the said a priori information, there is the additional advantage that close to the focal plane also the SNR is generally high, which has beneficial effects on the phase reconstruction.Paradoxically, for extremely low (and known) values of NA the knowledge of only one measured intensity pattern would suffice to solve the phase problem.This would make unnecessary any further study based on the concept of spatial correlation introduced in the present work.On the other hand, for either high NA systems or for systems whose NA is unknown, making use of the concept of spatial correlation, as shown here, can be decisive.In order to clarify this aspect, let us perform the phase reconstruction by using the following sets of data: • case 1: Intensity of the field in focus (z 1 = 0) and intensity of the field on z 2 = 4.04μm, with c int (0, z 2 ) = 0.82.
• case 2: Intensity of the field in focus (z 1 = 0) and intensity of the field on z 3 = 7.83μm, with c int (0, z 3 ) = 0.09 (first minimum of the correlation curve).
• case 3: Intensity of the field in focus (z 1 = and intensity of the field on z 4 = 11μm, with c int (0, z 4 ) = 0.29 (first relative maximum of the correlation curve).
• case 4: Intensity of the field in focus (z 1 = 0) and intensity of the field on z 5 = 15.17μm, with c int (0, z 5 ) = 0.03 (second minimum of the correlation curve).
It is difficult to compare the reconstructions in all these cases, since the data are characterized by differences in their SNR values.However, we can still draw some interesting conclusion.After running a phase retrieval algorithm in all these cases, assuming to know the NA, we have computed the root mean square error (RMSE) between the simulated intensity and the measured one at the position z = −4.62 μm.We have obtained (NA = 0.4 known) • case 1: RMSE = 0.0577.
From these values we see that although in case 1 the correlation between the intensity patterns used in the reconstruction is higher than all the other three cases, the obtained RMSE is the smallest.This is due to the fact that the field intensity at z 2 = 4.04 μm is characterized by the highest SNR among the cases considered above.In order to be more quantitative, we have computed the SNR for the measured intensity distributions (cases 1-4).SNR has been defined as the ratio between the average intensity and the relative standard deviation of a set of 50 samples of the same intensity distribution, at a given plane.For the reader's commodity, in Table 1 we report the values for the intensity correlations along with the SNR, RMSE for the planes listed above.The role of the a priori information can be pointed out by performing, once again, the reconstruction but assuming that the NA of the system is not known.For instance, for the case 1 and 2, under this assumption (NA unknown), we obtain • case 1: RMSE = 0.2016.
As expected, both reconstructions becomes worse in this case, because of the lack of a priori information (NA is supposed to be unknown now).However, it is also evident that the situation is completely reversed now.In fact, since there is no a priori information available that can compensate for eventual high correlation values, the rms error of case 1 becomes much larger than case 2, while with the said a priori information available they showed very similar RMSE values (with case 1 actually slightly better, as we have seen already).All this confirms our theoretical predictions.We expect that for high-NA systems choosing the proper plane for a phase reconstruction is a central element of the whole phase problem, considered that the even knowledge of NA does not help in that case.

Conclusions
We have described through-focus phase retrieval methods and their connection to the concept of spatial correlation for propagated fields.By making use of the concept of degree of spatially correlation for complex fields and intensities, we could predict which intensity distributions can be fruitfully used to start a phase retrieval procedure in a through focus configuration.An experimental setup has been built, and the phase of an optical focused field, generated by a microscope objective with NA = 0.4, has been recovered with excellent accuracy.We firmly believe that the tools we have introduced here to approach a phase retrieval problem can be safely extended to many more situations, not only to through-focus configurations.

#(
Fig.2.Degree of spatial correlation c int (0, z) between the intensity in focus and that in any other plane in the focal region, as function of the dimensionless variable zNA 2 /λ , with NA = 0.3.This curve corresponds to the case of a field with uniform angular spectrum A (0) (p, q) within NA.

Fig. 3 .
Fig. 3. Schematic of the experimental setup for through-focus phase retrieval.

#Fig. 4 .
Fig. 4. Panel a: measured degree of spatial correlation c int (0, z) between the intensity in focus and that in any other plane in the focal region, as function of the dimensionless variable zNA 2 /λ .The curve has been generated by using actual through-focus intensity measurements of a field focused by a microscope objective with nominal NA = 0.4.In panel b we plot the ideal degree of spatial correlation for the intensity (blue line) and field (red line) for an aberration-free system of same NA illuminated by a uniform field distribution.Panel c and d contain the 2D plots of the two intensity distributions chosen for the phase reconstruction.

#Fig. 5 .Fig. 6 .
Fig. 5. Measured intensities (first row) and retrieved phases (second row) of the field U(x, y, 0) in focus and the field at position z = 7.83 μm that has a minimum degree of spatial correlation with U(x, y, 0).

Table 1 .
Positions z, intensity correlation (c int (0, z)), signal to noise ratio (SNR) and root mean square error (RMSE) in the reconstruction of the field at z = −4.62 μm for different intensity measurements with known NA