1. Introduction

The non-crystallographic phase problem of coherent diffractive imaging, namely the reconstruction of a two-dimensional aperiodic complex scalar wavefield from the known squared modulus of its Fourier transform, has been subject to intensive research in recent years [1–9]. With historical roots that date back at least as far as the “Pauli problem” of determining a complex spatial wavefunction given knowledge of its probability density in both real-space and momentum-space [10], the CDI problem continues to be attacked using iterative algorithms that generalize the pioneering work of Gerchberg and Saxton [11]. Key developments include: (i) Fienup’s realization, that the Pauli–Gerchberg–Saxton assumption that the real-space modulus is known, could be replaced with the much weaker assumption that the support of the real-space wavefunction is known [12]; (ii) the experimental realization of CDI with visible light [1] using Fienup’s iterative reconstruction algorithms; (iii) the demonstration of Miao et al. [3], building on seminal work by Sayre [13] and Miao et al. [2], giving the first experimental realization of CDI for X-ray data given both a momentum-space intensity and a low-resolution real-space intensity; (iv) the experimental realisation of CDI from diffraction data alone [4]; (v) three-dimensional CDI [9]; (vi) CDI using intense, pulsed beams [14].

All of these approaches attack the problem of CDI from an iterative perspective, using what are in essence modified forms of the Gerchberg–Saxton–Fienup algorithm [12]. Such classes of algorithm iterate between real and momentum space, imposing known constraints in each of these spaces. Typical momentum-space constraints include a known intensity of the momentum-space wavefield (Fraunhofer diffraction pattern) [10, 11], with typical real-space constraints including known wave-field support [12], unknown but finite wave-field support [12], atomicity [15], positivity of real-space phase [2], and known coupling between real-space phase and real-space amplitude.

Such iterative approaches have reached a high state of accomplishment, with many reconstructions now having been reported in the literature [1–9, 14]. However, given that the CDI phase retrieval algorithms are iterative, the question naturally arises as to whether the problem of CDI is amenable to a closed-form analytic solution. Such a question of deterministic CDI has proven elusive, in part due to the fact that the phase of a CDI pattern is typically riddled with a complex distribution of branch points for screw-type phase singularities [16–18]. Notwithstanding this difficulty, it remains important to pursue deterministic solutions to the CDI problem, both for (i) the foundational and conceptual clarity that deterministic algorithms provide, and (ii) the possibility of enhanced speed of reconstruction, afforded by eschewing iteration and thereby realizing “single shot” phase–amplitude retrieval.

Here we report a particular analytical solution of the inverse problem of CDI, for a uniformly illuminated object placed within a rectangular hole with transverse dimensions at least two times larger than the object. The solution employs a modified Fourier transform [see Eq. (3)] of the far-field diffracted intensity. After application of our reconstruction algorithm to the Fraunhofer pattern, we may obtain eight independent phase/amplitude reconstructions of the original object, depending on the relative position of the object and the rectangular hole, as well as on the size of the hole.

 

Fig. 1. Setup for non-iterative Coherent Diffractive Imaging.

Download Full Size | PDF

2. Method

We assume the sample and rectangular hole, illustrated in Fig. 1, to be coherently illuminated with normally incident monochromatic scalar plane waves of uniform intensity I₀. The sample is assumed to lie within a square support of length and width no greater than 2Δ₀, this support lying entirely within a rectangular-shaped hole with centre at (x₀, y₀) and sizes of 2Δ_x≥4Δ₀ and 2Δ_y≥4Δ₀ in the horizontal and vertical directions, respectively. We emphasize that the constraining rectangular hole must be at least twice the size of the object, in both lateral directions, to avoid overlapping of the reconstructed object functions [see Eq. (4)].

One can write the spatial part of the wavefield just after the screen as

where, in the case of a unit-intensity incident plane monochromatic wave, ψ₀, defined within the hole, is the difference between the coherent complex scalar field at the exit surface of the object and unity, and ψ_S is a ‘support function’ that is by definition equal to unity within the entire area of the rectangular hole, and equal to zero outside the rectangular hole. Further, (x, y) denotes Cartesian coordinates in the object plane. Note that harmonic time dependence, of the coherent scalar field, is suppressed throughout.

Assume that a two-dimensional detector is placed at a sufficient distance downstream to yield a far-field diffraction pattern, with the surface of the detector being parallel to the screen. Neglecting the effects of noise, this far-field intensity I_D will be (see Appendix):

with a tilde indicating Fourier transformation with respect to x and y, and (q_s,q_y) being Fourier-space coordinates dual to (x, y). Here z is the distance between the object and the detector, the pair (x_d, y_d) denotes Cartesian coordinates in the detector plane, and the wavenumber k is related to the radiation wavelength λ via k=2π/λ.

To proceed further we introduce the auxiliary function U(x, y), which will be seen to provide a closed-form solution to the inverse problem of reconstructing the amplitude and phase of ψ₀ from I_D, via the definition:

Here, ÎD(qx,qy) is the intensity with noise. The procedure of introducing noise in the simulated intensity is described below. For real experimental data this will correspond to the registered intensity.

After some algebra (see Appendix) we can represent U(x, y) in the following form:

where ere sgn(x) is the signum function, and θ(x) is the step function [19]. Equations (3) and (4) are the main result of the article. The first term in Eq. (4) is the derivative of the object auto-correlation function, while the last term in Eq. (4) is the auto-correlation of the rectangular hole. Both of these functions are placed in the centre of the reconstructed image [see Fig.1]. The remaining terms represent eight laterally displaced reconstructions of the desired object function, half of which are complex conjugated, and at least six of which are separated from the object auto-correlation. The lateral displacements of the phase–amplitude reconstructions are determined by the position, (x₀, y₀), and the dimensions, (2Δ_x,2Δ_y), of the rectangular hole.

 

Fig. 2. Images corresponding to 45 dB SNR for the brightest pixel (N_max=10⁹): 2a (noisy Fraunhofer pattern), 2b - reconstructed function a(x, y) with d=0.255 and r=0.066, 2c - reconstructed function φ(x, y) with d=0.189 and r=0.050.

Download Full Size | PDF

 

Fig. 3. Images corresponding to 50 dB SNR for the brightest pixel (N_max=10¹⁰): 3a (noisy Fraunhofer pattern), 3b - reconstructed function a(x, y) with d=0.063 and r=0.016, 3c - reconstructed function φ(x, y) with d=0.044 and r=0.012.

Download Full Size | PDF

It should be noted that the reconstructions of the function ψ₀(x, y) shown in Eq. (4) are unique and unambiguous, apart from three “trivial” ambiguities irresolvable when only the modulus of the wavefunction can be registered [20]: a transverse spatial shift ψ₀(x+x₁, y+y₁), a constant additional phase ψ₀(x, y)e^iθand a complex conjugation ψ*₀(-x,-y). Note also that Podorov et al. have already used the technique of differentiation of the Fourier transform [see Eq. (3)], in a different context, to analyse the strain distribution in thin surface layers [21–22].

3. Numerical results

We now turn to numerical modelling of the reconstruction algorithm described above. First we simulate the Fraunhofer diffraction pattern from a test object with non-negative real amplitude a(x, y) and phase φ(x, y), so that ψ₀(x,y)=0(x,y)exp(iφ(x,y))-1. We use the “Lena” image as the phase function, with 0≤φ(x, y)≤3 radians, and a photo of Einstein as the amplitude function, with 0≤a(x, y)≤1 [see Figs 2, 3].

To demonstrate the robustness of our algorithm, we incorporate pseudo-random Poisson noise [see, e.g., [23]] in the intensity function I_D(q_x=kx_d/z,qy=ky_d/z) (i.e. the simulated Fraunhofer diffraction pattern) using the following procedure. First, we calculate the noise-free intensity I_D(qx=kx_d/z,qx=ky_d/z) [see Eq. (2)] for each detector pixel position, which is determined by the appropriate coordinates (x_d, y_d) in the detector plane. Second, we normalise this 2D map by multiplying by a constant N_max/max[I_D], where N_max is the number of photons collected in the pixel registering the highest intensity max[I_D]. This produces a 2D map, N(x_d,y_d), of the registered number of photons without the effects of noise. We then replace each of these numbers, N(x_d,y_d), by a random deviate, N̂(x_d,y_d), produced by a Poisson generator, for which 〈N̂(x_d,y_d)〉. Finally, we renormalise the new 2D map, N〈(x_d,y_d), to the new “noisy” intensity map Î_D(q_x=kx_d/z,qy=ky_d/z), by multiplying by the factor max[I_D]/N_max. Note that the signal-to-noise ratio (SNR) of the resulting image will vary with position in the detector plane, with the maximum SNR corresponding to the point in the diffraction pattern that has the maximum intensity N_max. Typically, this point will be in the vicinity of the centre of the pattern.

We use two levels for N_max in our simulations, namely 10⁹ and 10¹⁰ photons. This corresponds to SNR=10log₁₀(Signal Noise), for the brightest pixel, of 45 dB and 50 dB, respectively. The “noisy” simulated Fraunhofer diffraction patterns along with reconstructed images of the amplitude, a(x, y), and phase distribution, φ(x, y), for both levels of noise, are shown in Figs. 2(a), 2(b), 2(c) and Figs. 3(a), 3(b), 3(c), respectively. To estimate errors in the reconstructed images we employ two different metrics [24], namely a normalised root-mean-square (RMS) error criterion, defined as: $d=\sqrt{\sum {\left(A_{\mathrm{jk}}^{\mathrm{rec}}-A_{\mathrm{jk}}^{\mathrm{ideal}}\right)}^2⁄\sum {\left(A_{\mathrm{jk}}^{\mathrm{ideal}}-〈A^{\mathrm{ideal}}〉\right)}^2}$, and a normalized absolute difference, $r=\sqrt{\sum \mid A_{\mathrm{jk}}^{\mathrm{rec}}-A_{\mathrm{jk}}^{\mathrm{ideal}}\mid ⁄\sum \mid A_{\mathrm{jk}}^{\mathrm{ideal}}\mid }$, where A^ideal_jk and A^rec_jk are ideal and reconstructed two-dimensional images, respectively, and 〈A^ideal〉 is the mean of the original image, with all sums being taken over the integer pixel coordinates (j,k). The reconstruction results [Figs. 2(b), 2(c) and 3(b), 3(c)] show that a high-quality reconstruction requires a sufficiently intense registered signal, which, in the context of coherent X-ray optics, will typically require a filtered synchrotron source. Note that the reconstruction quality of the amplitude function, on, a(x, y), and the phase function, φ(x, y), is comparable to one another regardless of the level of noise. For instance, the d and r criteria for the a(x, y) images are 0.255 and 0.066 (for N_max=10⁹); compare to 0.063 and 0.016 (for N_max=10¹⁰). The reconstructed phase images, φ(x, y), have d and r criteria of 0.189 and 0.050 (for N_max=10⁹); compare to 0.044 and 0.012 (for N_max=10¹⁰), respectively. This is only marginally worse than for the amplitude images, a(x, y), for the lower SNR, and better for the higher SNR.

4. Discussion

It is important to elucidate the major effects that determine the resolving power of our method, i.e. the ability to resolve two neighbouring object points. There is a well-known theoretical criterion to the resolution limit Γ of an optical instrument, e.g., a microscope, determined by Abbe’s theory in the case of coherent illumination [see, e.g., reference [25] where Abbe’s theory [26] is used]: Γ=0.82λ/NA. Here, NA is the numerical aperture, i.e. half of the detector acceptance angle. Our simulations performed using the Fast Fourier Transform (FFT) algorithm, as well as the analytical form of the solution itself, show that there are several factors which can significantly deteriorate the quality of the reconstructed images. They are as follows: (i) the signal-to-noise ratio determined (apart from the object characteristics) by the intensity of the incident beam [see, e.g., [27] and references therein]. This effect is illustrated by Figs.2(b), 2(c) and 3(b), 3(c); (ii) the detector efficiency and dynamic range. For instance, a typical 2¹²–2¹⁴ dynamic range of CCD detectors necessitates a large number of subsequent images to be collected and averaged to get an acceptable SNR as well as using a beam stop to avoid saturation of the central part of the Fraunhofer pattern; (iii) edge smoothness and rectangularity of hole in the screen, which is important to obtain the second and third terms in Eq. (2) in the current form; (iv) the distance between the object and the detector, the detector size and number of pixels in the detector [see e.g., reference [25]].

Our theoretical simulations show that nano-scale resolution might be achievable. For instance, if one chooses an X-ray energy of 8 keV to collect the far-field pattern of an object contained within a 10µm×10µm square hole, using a detector with 50 µm pixel size placed at a distance of 7 meters from the object, the sampling step-size in the reconstructed real space object will be about 22 nm and the resolution Γ ≈ 36 nm.

It should be also noted that our method requires that the transverse coherence area covers the entire rectangular hole, not only the object itself. This implies the limitation that the area of reconstruction in the real space be at least two times larger than (and, therefore, the sampling spacing in the reciprocal space should be two times finer than) that required by current Gerchberg–Saxton-type algorithms (which typically consider tighter supports).

As one can see from Eq. (3) the part of the diffraction pattern with either q_x = 0 or q_y=0 can be easily excluded from the reconstruction process, because these terms are multiplied by 0. This means that a beam stop can cover the central vertical (q_x=0) and horizontal (q_y=0) lines of the detector pixels without affecting the reconstruction process. Alternatively, the analysis may proceed by “stitching together” two images: (i) the Fraunhofer pattern obtained with a beam stop, and (ii) a short-exposure image of the central portion of the Fraunhofer pattern.

Finally, let us compare our method with two well-established holography techniques, namely inline (on-axis) holography [26, 28] and Fourier holography [26, 29–34]. Regarding the comparison of our technique to inline holography, which corresponds to the limit case Δ_x,Δ_y→∞, we see that the introduction of a confining rectangular hole allows us to remove the twin-image problem [28] of the latter method. The physical origin of this twin-image removal is that the “on axis” reference wave of the Gabor theory (i.e. the unscattered plane wave) is replaced with the “off axis” Young–Maggi–Rubinowicz “boundary wave” [26], scattered from the four edges of the rectangular hole. Regarding the comparison of our technique with Fourier holography: in Fourier holography the reference wave is provided by the radiation transmitted through a small hole or scattered from a small point, with either the small hole or the point scatterer being laterally displaced from the sample, while in our approach the reference wave is provided by the “boundary wave” scattered by the edge of the rectangular hole. A further point of difference: in Fourier holography the resolution is limited by the size of the reference hole or reference scatterer, which, in the ideal case scenario, should be point-like, while in our method the resolution is limited by the sharpness of the edges of the confining rectangular hole. Note, in this context, that our method is also applicable to the case where the object wave does not overlap with the rectangular hole.

5. Conclusions

In summary, we proposed and successfully employed a deterministic reconstruction algorithm for coherent diffractive imaging, utilizing simulated data. This involves collecting a far-field diffraction pattern from a normally-illuminated object of finite support, that lies within the field of view of a rectangular hole of transverse dimensions at least twice that of the object. Our robust closed-form algorithm allows us to unambiguously reconstruct the scalar complex wavefield (in nanoscale range) just after the object using a single far-field diffraction pattern.