Non-iterative method for phase retrieval and coherence characterization by focus variation using a fixed star-shaped mask

: A novel non-iterative phase retrieval method is proposed and demonstrated with a proof-of-principle experiment. The method uses a ﬁxed specially designed mask and through-focus intensity measurements. It is demonstrated that this method is robust to spatial partial coherence in the illumination, making it suitable for coherent diﬀractive imaging using spatially partially coherent light, as well as for coherence characterization.


Introduction
Coherent Diffractive Imaging (CDI) has emerged as a popular method for lensless imaging and phase imaging. In CDI, one illuminates a thin sample with coherent light, and one or multiple intensity patterns of the diffracted transmitted field are measured (for example in the far field). With this set of intensity measurements, the complex-valued transmission function of the sample is reconstructed. Several reconstruction schemes have been developed, and these can be divided in two categories: iterative reconstruction algorithms, and non-iterative methods.
The most basic iterative reconstruction algorithms use a single far-field intensity measurement and a support constraint that assumes we know the shape of the object a priori. Examples of such algorithms are the Error Reduction and Hybrid Input-Output algorithms [1,2]. Other algorithms use multiple intensity measurements. For example, in ptychography the object is illuminated with a spot of light referred to as 'the probe'. By shifting the probe to different overlapping positions (or equivalently, by moving the object), different far-field intensity patterns are obtained, and the object is reconstructed from this set of diffraction patterns [3]. Another method is to take set of through-focus intensity measurements and reconstruct the object from these patterns [4].
Non-iterative methods are for example phase-shifting holography (also sometimes referred to as quantitative Zernike phase contrast imaging [5]) and Fourier holography [6]. In these methods, the transmitted field is modulated in a single point (ideally) so that in the far-field the diffracted field interferes with a reference plane wave. The object is then reconstructed using methods that are equivalent to phase-shifting interferometry and off-axis holography. Other examples of non-iterative methods include single-parameter phase modulation [7], focus-variation that assumes a high zeroth diffraction order [8], solving the Transport of Intensity Equation [9], and using shifting Gaussian filters [10].
The reconstruction methods that were just mentioned all assume that the illumination is perfectly coherent. It may be relevant to consider how the algorithms perform or how they must be adapted when the illumination is spatially partially coherent. For the iterative solutions, the solution usually involves mode decomposition. It can be shown that the mutual coherence function can be decomposed into modes that propagate coherently, but are mutually incoherent [11]. In the completely coherent case only one mode is required to describe the transmitted field, and as the coherence decreases, more modes are required. If it is known a priori that the mutual coherence function is shift invariant, the diffracted intensity can be calculated more easily using a convolution [12]. For several non-iterative methods such as Fourier holography and phase-shifting holography, it has been demonstrated [13] that these methods also work when the illumination is partially coherent, and that in the absence of an object these methods can be used to reconstruct the coherence function of the illumination.
In this article, a method for a non-iterative object reconstruction from through-focus intensity measurements is presented. The method works by using a fixed specially designed mask. This method is robust to partial coherence, and can also be used to characterize the coherence structure of a field. If it can be assumed that the coherence function and/or the object transmission function is real-valued, only half of the focal field needs to be scanned, which means that it is also possible to reconstruct the field using just free-space propagation, thus eliminating the need for focusing optics. Since it is a non-iterative method, it is computationally inexpensive. The method does not require mode decomposition (as opposed to the iterative methods), which therefore avoids the arbitrary truncation of the number of modes that are used to describe the coherence function. Focus-variation is achieved by simply moving the detector, or by modulating the wave front with a spherical or quadratic phase pattern, which is a relatively simple modulation that can be introduced easily (using e.g. a liquid lens for optical wavelengths, or electromagnetic lenses for electrons), especially compared to other phase modulation methods [7].

Method
Suppose we have a mutual coherence function J(x 1 , x 2 ) of a quasi-monochromatic spatially partially coherent field of which we want to reconstruct information using a through-focus intensity scan. The intensity in the observation plane is given by where A is a defocus parameter and is given by (see Appendix A) where z f is the distance to the back focal plane of a lens with focal length f , and λ is the wavelength of the light. The reciprocal space coordinate u is related to the lateral focal field coordinates (x f , y f ) as Inverse Fourier transforming this intensity pattern with respect to both u and A gives If for the sake of notation we write x 2 = y, then we get Our goal is to use the delta function in this integral to filter out information about J(x 1 , x 2 ) directly. To do this, we choose ω A = |x| 2 + 2x · P (see Appendix B for details), where P is a fixed point, the choice of which will be explained later. We get Due to this delta function, the integrand contributes to the integral only when x = 0 or when x is perpendicular to P − y. Let us assume in the following that x 0. What we want to do now, is to make sure that the integrand contributes only when y = P, so that we can directly reconstruct the object modulated by the mutual coherence function of the illumination Let us for now consider the case of full coherence in which case we can write J(x + P, can be interpreted as the transmission function of a sample that we want to reconstruct. If we want a nonzero contribution to the integral only if y = P, we require the following O(x + y) = 0 or O(y) * = 0 when P − y is perpendicular to x and y P.
We can make sure that these requirements hold by imposing a certain mask on the object, and we can check these requirements visually as follows: • Sketch O(x) and check that O(P) * 0, i.e. P should lie within the domain of the object (which can be defined by a mask). This fulfills requirement 1.
• Sketch O(x + P) as a function of x and see for which x it holds that O(x + P) 0. These are the x we are interested in. This fulfills requirement 2.
• For all the directions of x that we are interested in, draw lines perpendicular to x through P, and lines perpendicular to x through x + P. These lines make up the collection of y and x + y respectively for which x · (P − y) = 0. Verify that there is no y except for y = P for which both O(x + y) and O(y) * are nonzero. This fulfills requirement 3.
In Fig. 1 this procedure is used to demonstrate that a star-shaped mask is suitable for non-iterative focus-variation reconstruction: because of the sharp protrusions of the star, the line that is perpendicular to any relevant x and goes through P will intersect O(x) in only a small region, even if the line has a finite thickness due to the finite sampling range of A (as will be explained in Section 2.1). Therefore, only a small region of y around P will contribute to the integral in Eq. (6). On a practical note, it is important to recall that in the reconstruction formula in Eq. (7), Fig. 1. Demonstration that a star-shaped mask is suitable for non-iterative focus-variation reconstruction. If P is chosen to be one of the star's points, then we see that for any x for which O(x + P) is non-zero, all y for which P − y is perpendicular to x intersect the object only in P.
F −1 u, A {·} denotes the three-dimensional inverse Fourier transform with respect to both u and A. In practice however, it would be unnecessarily cumbersome to compute a three-dimensional inverse Fourier transform and extract the desired two-dimensional surface from it. Therefore, we instead perform the reconstruction by summing two-dimensional inverse Fourier transforms with respect to u for all A J(x + P, P) ∝ where in this case, F −1 u {·} denotes the two-dimensional inverse Fourier transform with respect to u.

Finite scanning range of A
When performing the Fourier transform with respect to A as in Eq. (4), note that we can scan A only over a finite range. As a result, we will not have a sharp delta peak in ω A , but rather we have a function a finite width ∆ω A . We can mitigate the effects of a finite sampling range by introducing a sampling function H(A) If we take a rectangular sampling window H(A) = 1, taking the inverse Fourier transform with respect to A gives a factor proportional to where W is the sampling range of A. This function gets narrower when W or |x| gets larger, but it always has significant sidelobes. By choosing H(A) to be a Hamming window, these sidelobes can be suppressed, as is shown in Fig. 2. Another problem that may occur is that if one varies A by varying the propagation distance in the Fresnel approximation, rather than focusing the field, then only positive A can be scanned. In this case, A is given by and u is related to the real-space coordinates (x, y) as where λ is the wavelength and z is the propagation distance. If it is known that the mutual coherence function is real-valued (i.e. J(x 1 , x 2 ) = J(x 1 , x 2 ) * ), the intensity measurements obtained for positive A can be used to infer the intensity patterns for negative A.
This extrapolation method as well as others are discussed in [14].

Experiment
To test the method in a proof-of-principle experiment, we use the setup as shown in Fig. 3. In this setup, we use green laser light (λ = 532nm) which is expanded in the beam expander. The beam is made partially coherent by focusing it on a rotating ground-glass disk, after which the beam is approximately collimated again with a second lens. By varying the position of the first lens, we can vary the spot size on the rotating disk, thereby changing the degree of coherence. The resulting beam has a Gaussian correlation structure. The partially coherent beam is incident on a reflective liquid crystal phase-only Spatial Light Modulator (SLM, HOLOEYE-GAEA-VIS, 3840×2160 resolution, 3.74 µm pixel size), on which a pattern is assigned that serves as the object that is to be reconstructed. The radius of the object on the SLM is R = 2.62mm, see Fig. 4. Using a third lens with focal length f = 15cm, the modulated light is focused, and in the back focal plane of the lens, the intensity pattern is recorded with a CCD camera. In order to generate the through-focus data set, quadratic phase factors are added to the pattern assigned to the SLM. The defocus parameter AR 2 is varied from -7 to 7 in 100 steps, which according to Eq. (2) means that for the given parameters this is equal to physically scanning through the focal field from z f = −2.45cm to z f = 2.45cm. In case one knows a priori that the sample and the coherence function of the illumination are real-valued, the third lens is in principle not necessary, since the intensity patterns obtained with free-space propagation is already sufficient for the reconstruction. With this experimental setup, we take two data sets each consisting of 100 through-focus images, where for each data set a different degree of coherence has been used. The coherence widths of the Gaussian-correlated beam with constant amplitude at the SLM-plane are σ = 2.3mm (high coherence) and σ = 0.5mm (low coherence). These values have been obtained using coincidence measurements [15].
In Fig. 5 the non-iterative reconstructions are plotted, and compared to the simulated amplitudes which are obtained by multiplying the amplitude of O(x + P) with the Gaussian correlation function with respect to the reference point P:  It is thus shown that the coherence structure of the illuminating field limits the field of view of the reconstruction, and conversely, the degree of coherence can be inferred directly from the non-iterative reconstruction.

Extending the field of view
For each P, a coherence function J 0 (x) = J(0, x) 'illuminates' a different part of the object O(x). In the example of Eq. (14), we have J 0 (x) = e − |x | 2 /2σ 2 . So by shifting P around, we 'illuminate' different parts of the object, thereby extending the field of view. In Fig. 6 it is shown how the field of view can be extended by considering multiple reference points P. In order to synthesize an object reconstruction O(x) with an extended field of view from the set of reconstructions for different P, we use a factorization method that was originally used for ptychography [16]. Let us define These ψ j (x) are the reconstructions we obtain using the non-iterative method as shown in Fig. 6, except that we have shifted them by P to the center. Also, we have multiplied the reconstructions with a window function to eliminate the artifacts of the non-iterative reconstruction, as shown in the top two rows of Fig. 7. We denote the collection of all ψ j (x) as Ψ. Next, we define the (a) High coherence, σ = 2.3mm.
where is a small constant to prevent division by 0. To put it in words, the operator is known. Since in our case both O(x) and J 0 (x) are unknown, we have to apply an iterative procedure to reconstruct O(x) and J 0 (x). Denoting the k th estimate of O(x) and J 0 (x) as O (k) (x) and J (k) 0 (x) respectively, one iteration of the procedure is defined as follows Here µ is a step size that should be chosen between 0 and 1. In order to obtain the results in the bottom row of Fig. 7, we first applied 10 iterations with µ = 1, then 10 iterations with µ = 0.1. It must be noted that even though this is technically an iterative method, it is a computationally inexpensive one: in ptychography, this iterative method is applied in each iteration of the algorithm.

Characterization of a Laguerre-Gaussian correlated Schell-model beam
To illustrate how the non-iterative method can be used to characterize the coherence of a partially coherent beam, we generate a Laguerre-Gaussian correlated Schell-model beam (LGCSM) by inserting a spiral phase plate of topological charge n = 2 before the lens L1 in Fig. 3. The resulting coherence function is given by [15,17] where n = 2 is the topological charge of the spiral phase plate, and L 0 n denotes the Laguerre polynomial of mode order n and 0. The parameter we aim to find is the coherence width σ. To do this, we perform a non-iterative through-focus reconstruction, and take two cross-sections, taken from reconstructions obtained with different P, as shown in Fig. 8. These two cross-sections are then fitted to Eq. (18) using a least-squares method, which gives a value of σ = 0.80mm. With this parameter value, one can compute the coherence function as prescribed by Eq. (18).

Conclusion
In this article we have proposed a non-iterative phase retrieval method using a fixed specially designed mask and a through-focus data set of intensity patterns. We have demonstrated mathematically how the method works, and that it is robust to partial spatial coherence in the illumination. With a proof-of-principle experiment using visible laser light and a spatial light modulator to create an object and scan through focus, we tested the proposed algorithm for different degrees of coherence of the illumination. The non-iterative reconstruction directly reveals information about the coherence structure of the illuminating field. Focus-variation can be introduced relatively easily, especially when the object is known to be real-valued, in which case only half of the focal field needs to be scanned, which means one does not necessarily need to focus the transmitted field, but rather let it propagate freely. Therefore, this method can be a valuable contribution to the field of diffractive imaging and coherence characterization.

A. Fresnel propagation of the focal field
The focal field is given byÛ where (x f , y f ) denote the coordinates in the back focal plane, λ is the wavelength, f is the focal length of the lens, and (x, y) are the coordinates in the front focal plane. We propagate this field by z f using the Angular Spectrum method. The spectrum ofÛ(x f , y f ) is (20) In the Fresnel approximation the propagator is given by Thus, for the propagated field we get x f λ f , y f λ f .

(22)
If we define R to be the size of the object, so that (X, Y ) = x R , y R are normalized coordinates, we can writeÛ

B. Conditions for ω A
We have where v = µx − 2y Our goal is to use the delta function in this integral to filter out information about J(x) directly.

B.1. Condition for ω A given y = P
We want to choose ω A such that the delta function gives a contribution for y = P. We require for any constant λ. We also know by definition (Eq. (25)), and given that y = P Combining Eq. (27) and (28) gives us the matrix equation Solving this equation gives and from Eq. (25) we then get B.2. Condition for ω A given x + y = P We want to choose ω A such that the delta function gives a contribution for y + x = P. In this case we have to substitute in Eq. (28) P with P − x, which gives v = (µ + 2) We can quickly see from Eq. (33) that in this case and from Eq. (25) we then get ω A = −|x| 2 + 2P · x.