We develop a deterministic algorithm for coherent diffractive imaging (CDI) that employs a modified Fourier transform of a Fraunhofer diffraction pattern to quantitatively reconstruct the complex scalar wavefield at the exit surface of a sample of interest. The sample is placed in a uniformly-illuminated rectangular hole with dimensions at least two times larger than the sample. For this particular scenario, and in the far-field diffraction case, our non-iterative reconstruction algorithm is rapid, exact and gives a unique analytical solution to the inverse problem. The efficacy and stability of the algorithm, which may achieve resolutions in the nanoscale range, is demonstrated using simulated X-ray data.
© 2007 Optical Society of America
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 , the CDI problem continues to be attacked using iterative algorithms that generalize the pioneering work of Gerchberg and Saxton . 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 ; (ii) the experimental realization of CDI with visible light  using Fienup’s iterative reconstruction algorithms; (iii) the demonstration of Miao et al. , building on seminal work by Sayre  and Miao et al. , 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 ; (v) three-dimensional CDI ; (vi) CDI using intense, pulsed beams .
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 . 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 , unknown but finite wave-field support , atomicity , positivity of real-space phase , 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.
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 I0. The sample is assumed to lie within a square support of length and width no greater than 2Δ0, this support lying entirely within a rectangular-shaped hole with centre at (x0, y0) and sizes of 2Δx≥4Δ0 and 2Δy≥4Δ0 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, ψ0, 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 ID will be (see Appendix):
with a tilde indicating Fourier transformation with respect to x and y, and (qs,qy) being Fourier-space coordinates dual to (x, y). Here z is the distance between the object and the detector, the pair (xd, yd) 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 ψ0 from ID, 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 . 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, (x0, y0), and the dimensions, (2Δx,2Δy), of the rectangular hole.
It should be noted that the reconstructions of the function ψ0(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 : a transverse spatial shift ψ0(x+x1, y+y1), a constant additional phase ψ0(x, y)eiθand a complex conjugation ψ*0(-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 ψ0(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., ] in the intensity function ID(qx=kxd/z,qy=kyd/z) (i.e. the simulated Fraunhofer diffraction pattern) using the following procedure. First, we calculate the noise-free intensity ID(qx=kxd/z,qx=kyd/z) [see Eq. (2)] for each detector pixel position, which is determined by the appropriate coordinates (xd, yd) in the detector plane. Second, we normalise this 2D map by multiplying by a constant Nmax/max[ID], where Nmax is the number of photons collected in the pixel registering the highest intensity max[ID]. This produces a 2D map, N(xd,yd), of the registered number of photons without the effects of noise. We then replace each of these numbers, N(xd,yd), by a random deviate, N̂(xd,yd), produced by a Poisson generator, for which 〈N̂(xd,yd)〉. Finally, we renormalise the new 2D map, N〈(xd,yd), to the new “noisy” intensity map ÎD(qx=kxd/z,qy=kyd/z), by multiplying by the factor max[ID]/Nmax. 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 Nmax. Typically, this point will be in the vicinity of the centre of the pattern.
We use two levels for Nmax in our simulations, namely 109 and 1010 photons. This corresponds to SNR=10log10(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 , namely a normalised root-mean-square (RMS) error criterion, defined as: , and a normalized absolute difference, , where Aidealjk and Arecjk are ideal and reconstructed two-dimensional images, respectively, and 〈Aideal〉 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 Nmax=109); compare to 0.063 and 0.016 (for Nmax=1010). The reconstructed phase images, φ(x, y), have d and r criteria of 0.189 and 0.050 (for Nmax=109); compare to 0.044 and 0.012 (for Nmax=1010), respectively. This is only marginally worse than for the amplitude images, a(x, y), for the lower SNR, and better for the higher SNR.
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  where Abbe’s theory  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.,  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 212–214 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 ].
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 qx = 0 or qy=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 (qx=0) and horizontal (qy=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  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” , 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.
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.
Here we give a detailed derivation of the key results presented in the main text of this paper. First, to obtain the far-field intensity, ID, [see Eq. (2)], we calculate the propagated wavefield using Fresnel diffraction theory as :
Note that, in line 2 of Eq. (A4) we used a tabulated integral [see e.g., p. 14 in ]. The remaining terms, i.e. B 1,2,3,4 and C 1,2,3,4 are straight-forward results of the Fourier shift theorem [see e.g., p. 8 in ]. For instance,
The authors acknowledge funding from the Australian Research Council.
References and links
2. 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). [CrossRef]
3. 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). [CrossRef]
4. S. Marchesini, H. He, H. N. Chapman, S. P. Hau-Riege, A. Noy, M. R. Howells, U. Weierstall, and J. C. H. Spence, “X-ray image reconstruction from a diffraction pattern alone,” Phys. Rev. B 68, 140101(R) (2003). [CrossRef]
7. I. K. Robinson, I. A. Vartanyants, G. J. Williams, M. A. Pfeifer, and J. A. Pitney, “Reconstruction of the shapes of gold nanocrystals using coherent x-ray diffraction,” Phys. Rev. Lett. 87, 195505 (2001). [CrossRef] [PubMed]
9. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Physics 2, 258–261 (2006). [CrossRef]
10. W. PauliH. Geiger and K. Scheel, “Die allgemeinen Prinzipien der Wellenmechanik” in Handbuch der Physik, ed. (Springer, Berlin, 1933) 24(1), 83–272.
11. R. W. Gerchberg and W.O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik , 35, 237–246 (1972).
13. D. Sayre, “Some implications of a theorem due to Shannon,” Acta Cryst. 5, 843–843 (1952). [CrossRef]
14. H. N. Chapman, A. Barty, M. J. Bogan, S. Boutet, M. Frank, S. P. Hau-Riege, S. Marchesini, B. W. Woods, S. Bajt, W. H. Benner, R. A. London, E. Plönjes, M. Kuhlmann, R. Treusch, S. Düesterer, T. Tschentscher, J. R. Schneider, E. Spiller, T. Möller, C. Bostedt, M. Hoener, D. A. Shapiro, K. O. Hodgson, D. van der Spoel, F. Burmeister, M. Bergh, C. Caleman, G. Huldt, M. M. Seibert, F. R. N.C. Maia, R. W. Lee, A. Szöke, N. Timneanu, and J. Hajdu, “Femtosecond diffractive imaging with a soft-X-ray free-electron laser,” Nature Physics 2, 839–843 (2006). [CrossRef]
15. V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Am. A 20, 40–55 (2003). [CrossRef]
16. P. A. M. Dirac, “Quantised singularities in the electromagnetic field,” Proc. R. Soc. A 133, 60–72 (1931). [CrossRef]
17. M. V. BerryR. Balian “Singularities in waves and rays,” et al. (eds), Les Houches Lecture Series, session XXXV, Physics of defects (North Holland, Amsterdam, 1981) pp. 453–543
18. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, New York, 2006). [CrossRef]
19. G. A. Korn and T. M. Korn, Mathematical handbook for scientists and engineers, 2nd ed (McGraw-Hill Book Company, 1968).
20. R. H. T. Bates, “Fourier phase problems are uniquely solvable in more than one dimension. I: Underlying theory,” Optik 61(3), 247–262 (1982).
21. S. G. Podorov, G. Hölzer, E. Förster, and N. N. Faleev, “Fourier analysis of X-ray rocking curves from superlattices,” Phys. Stat. Sol. (b) 213, 317–324 (1999). [CrossRef]
22. S. G. Podorov, G. Hölzer, E. Förster, and N. N. Faleev, “Semidynamical solution of the inverse problem of X-ray Bragg diffraction on multilayered crystals,” Phys. Stat. Sol. (a) 169, 9–16 (1998). [CrossRef]
23. J. W. Goodman, Statistical optics (John Wiley & Sons, Inc., 2000).
24. G. T. Herman, Image reconstruction from projections (Academic Press, 1980).
25. G. J. Williams, H. M. Quiney, B. B. Dhal, C. Q. Tran, K. A. Nugent, A. G. Peele, D. Paterson, and M. D. de Jonge, “Fresnel coherent diffractive imaging,” Phys. Rev. Lett. 97, 025506 (2006). [CrossRef] [PubMed]
26. M. Born and E. Wolf, Principles of Optics, 7th edition (Cambridge University Press, Cambridge, 1999).
27. Q. Shen, I. Bazarov, and P. Thibault, “Diffractive imaging of nonperiodic materials with future coherent sources,”. J. Synchr. Rad. 11, 432–438 (2004). [CrossRef]
29. J. T. Winthrop and C. R. Worthington, “X-ray microscopy by successive Fourier transformation,” Phys. Lett. 15, 124–126 (1965). [CrossRef]
32. G. W. Stroke, “Lensless Fourier-transform method for optical holography,” Appl. Phys. Lett. 6(10), 201–203 (1965). [CrossRef]
33. G. W. Stroke and D. G. Falconer, “Attainment of high resolution in wavefront-reconstruction imaging,” Phys. Lett. 13, 306–309 (1964). [CrossRef]
34. G. W. Stroke, An introduction to coherent optics and holography, (Academic Press, New York, London, 1966).
35. J. W. Goodman, Introduction to Fourier optics, 3rd edition (Roberts & Company, Englewood, Colorado, 2005).