Reconstruction of wave front and object for inline holography from a set of detection planes

We illustrate the errors inherent in the conventional empty beam correction of full field X-ray propagation imaging, i.e. the division of intensities in the detection plane measured with an object in the beam by the intensity pattern measured without the object, i.e. the empty beam intensity pattern. The error of this conventional approximation is controlled by the ratio of the source size to the smallest feature in the object, as is shown by numerical simulation. In a second step, we investigate how to overcome the flawed empty beam division by simultaneous reconstruction of the probing wavefront (probe) and of the object, based on measurements in several detection planes (multi-projection approach). The algorithmic scheme is demonstrated numerically and experimentally, using the defocus wavefront of the hard X-ray nanoprobe setup at the European Synchrotron Radiation Facility (ESRF). © 2014 Optical Society of America OCIS codes: (340.7440) X-ray imaging; (100.5070) Phase retrieval. References and links 1. P. Cloetens, W. Ludwig, J. Baruchel, D. van Dyck, J. van Landuyt, J. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. Lett. 75, 2912–2914 (1999). 2. D. M. Paganin, Coherent x-ray optics (New York: Oxford University Press, 2006). 3. K. A. Nugent, “The measurement of phase through the propagation of intensity: An introduction,” Contemp. Phys. 52, 55–69 (2011). 4. H. M. Quiney, “Coherent diffractive imaging using short wavelength light sources,” J. Mod. Opt. 57, 1109–1149 (2010). 5. D. R. Luke, J. V. Burke, and R. G. Lyon, “Optical wavefront reconstruction: Theory and numerical methods,” SIAM Rev. 44, 169–224 (2002). 6. T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. 133, 339–346 (1997). 7. M. Krenkel, M. Bartels, and T. Salditt, “Transport of intensity phase reconstruction to solve the twin image problem in holographic x-ray imaging,” Opt. Express 21, 2220–2235 (2013). 8. E. Maire, J.-Y. Buffière, L. Salvo, J. J. Blandin, W. Ludwig, and J. M. Letang, “On the application of x-ray microtomography in the field of materials science,” Adv. Eng. Mater. 3, 539–546 (2001). 9. P. Tafforeau, R. Boistel, E. Boller, A. Bravin, M. Brunet, Y. Chaimanee, P. Cloetens, M. Feist, J. Hoszowska, and J.-J. Jaeger, “Applications of x-ray synchrotron microtomography for non-destructive 3D studies of paleontological specimens,” Appl. Phys. A Mater. 83, 195–202 (2006). #203649 $15.00 USD Received 10 Jan 2014; revised 25 Mar 2014; accepted 30 Mar 2014; published 6 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.011552 | OPTICS EXPRESS 11552 10. B. D. Metscher, “Biological applications of x-ray microtomography: Imaging microanatomy, molecular expression and organismal diversity,” Microsc. Anal. (Am Ed) 27, 13–16 (2013). 11. T. van de Kamp, P. Vagovič, T. Baumbach, and A. Riedel, “A biological screw in a beetle’s leg,” Science 333, 52–52 (2011). 12. J. Moosmann, A. Ershov, V. Altapova, T. Baumbach, M. S. Prasad, C. LaBonne, X. Xiao, J. Kashef, and R. Hofmann, “X-ray phase-contrast in vivo microtomography probes new aspects of xenopus gastrulation,” Nature 497, 374 – 377 (2013). 13. M. Bartels, V. H. Hernandez, M. Krenkel, T. Moser, and T. Salditt, “Phase contrast tomography of the mouse cochlea at microfocus x-ray sources,” Appl. Phys. Lett. 103, 083703 (2013). 14. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68, 2774–2782 (1997). 15. P. Cloetens, R. Mache, M. Schlenker, and S. Lerbs-Mache, “Quantitative phase tomography of Arabidopsis seeds reveals intercellular void network,” Proc. Natl. Acad. Sci. USA 103, 14626–14630 (2006). 16. M. Bartels, M. Priebe, R. Wilke, S. Kruger, K. Giewekemeyer, S. Kalbfleisch, C. Olendrowitz, M. Sprung, and T. Salditt, “Low-dose three-dimensional hard x-ray imaging of bacterial cells,” Opt. Nanoscopy, 1-10 (2012). 17. K. Giewekemeyer, S. Krüger, S. Kalbfleisch, M. Bartels, C. Beta, and T. Salditt, “X-ray propagation microscopy of biological cells using waveguides as a quasipoint source,” Phys. Rev. A 83, 023804 (2011). 18. C. Olendrowitz, M. Bartels, M. Krenkel, A. Beerlink, R. Mokso, M. Sprung, and T. Salditt, “Phase-contrast x-ray imaging and tomography of the nematode Caenorhabditis elegans,” Phys. Med. Biol. 57, 5309–5323 (2012). 19. G. Martinez-Criado, R. Tucoulou, P. Cloetens, P. Bleuet, S. Bohic, J. Cauzid, I. Kieffer, E. Kosior, S. Laboure, S. Petitgirard, A. Rack, J. A. Sans, J. Segura-Ruiz, H. Suhonen, J. Susini, and J. Villanova, “Status of the hard x-ray microprobe beamline ID22 of the European Synchrotron Radiation Facility,” J. Synchrotron Radiat. 19, 10–18 (2012). 20. S. Kalbfleisch, H. Neubauer, S. P. Krüger, M. Bartels, M. Osterhoff, D. D. Mai, K. Giewekemeyer, B. Hartmann, M. Sprung, and T. Salditt, “The Göttingen holography endstation of beamline P10 at PETRA III/DESY,” AIP Conf. Proc. 1365, 96–99 (2011). 21. R. Barrett, R. Baker, P. Cloetens, Y. Dabin, C. Morawe, H. Suhonen, R. Tucoulou, A. Vivo, and L. Zhang, “Dynamically-figured mirror system for high-energy nanofocusing at the ESRF,” Proc. SPIE pp. 813904–813912 (2011). 22. K. Giewekemeyer, “A study on new approaches in coherent x-ray microscopy of biological specimens,” Ph.D. thesis, Universität Göttingen (2011). 23. P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy 109, 338–343 (2009). 24. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109, 1256–1262 (2009). 25. M. Stockmar, P. Cloetens, I. Zanette, B. Enders, M. Dierolf, F. Pfeiffer, and P. Thibault, “Near-field ptychography: Phase retrieval for inline holography using a structured illumination,” Sci. Rep. 3, 1–6 (2013). 26. H. M. Quiney, A. G. Peele, Z. Cai, D. Paterson, and K. A. Nugent, “Diffractive imaging of highly focused x-ray fields,” Nature 2, 101–104 (2006). 27. F. Döring, A.-L. Robisch, C. Eberl, M. Osterhoff, A. Ruhlandt, T. Liese, F. Schlenkrich, S. Hoffmann, M. Bartels, T. Salditt, and H. U. Krebs, “Sub-5 nm hard x-ray point focusing by a combined Kirkpatrick-Baez mirror and multilayer zone plate,” Opt. Express 21, 19311–19323 (2013). 28. T. B. Edo, D. J. Batey, A. M. Maiden, C. Rau, U. Wagner, Z. D. Pešić, T. A. Waigh, and J. M. Rodenburg, “Sampling in x-ray ptychography,” Phys. Rev. A 87, 053850 (2013). 29. D. L. Misell, “An examination of an iterative method for the solution of the phase problem in optics and electron optics: I. Test calculations,” J. Phys. D: Appl. Phys. 6, 2200–2216 (1973). 30. B. M. Hanser, M. G. L. Gustafsson, D. A. Agard, and J. W. Sedat, “Phase-retrieved pupil functions in wide-field fluorescence microscopy,” J. Microsc. 216, 32–48 (2004). 31. L. J. Allen, W. McBride, and M. P. Oxley, “Exit wave reconstruction using soft x-rays,” Opt. Commun. 233, 77–82 (2004). 32. L. J. Allen, W. McBride, N. O’Leary, and M. P. Oxley, “Exit wave reconstruction at atomic resolution,” Ultramicroscopy 100, 91–104 (2004). 33. D. R. Luke, “Relaxed averaged alternating reflections for diffraction imaging,” Inverse Probl. 21, 37–50 (2005). 34. S. C. Kramer and J. Hagemann, “SciPAL: Expression templates and composition closure objects for high performance computational physics with CUDA and OpenMP,” ACM Trans. on Parallel Comput. (submitted). 35. J. N. Clark, C. T. Putkunz, M. A. Pfeifer, A. G. Peele, G. J. Williams, B. Chen, K. A. Nugent, C. Hall, W. Fullagar, S. Kim, and I. McNulty, “Use of a complex constraint in coherent diffractive imaging,” Opt. Express 18, 1981– 1993 (2010). 36. L. J. Allen, H. M. L. Faulkner, K. A. Nugent, M. P. Oxley, and D. Paganin, “Phase retrieval from images in the presence of first-order vortices,” Phys. Rev. E 63, 037602 (2001). 37. L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. 199, 65–75 (2001). #203649 $15.00 USD Received 10 Jan 2014; revised 25 Mar 2014; accepted 30 Mar 2014; published 6 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.011552 | OPTICS EXPRESS 11553 38. A.-L. Robisch and T. Salditt, “Phase retrieval for object and probe using a series of defocus near-field images,” Opt. Express 21, 23345–23357 (2013). 39. C. Homann, J. Hagemann, T. Salditt, and T. Hohage, “Remarks on the product approximation in the detector plane,” (in prepartion).


Introduction
X-ray propagation full field imaging and tomography have become powerful imaging techniques, based on successful development of phase retrieval algorithms [1][2][3][4][5]. In the detection plane, the contrast for both absorbing and phase shifting features is generated by the inherent near-field diffraction behind a sample which is illuminated (probed) by an extended x-ray wavefront P of sufficiently high spatial coherence. The exit wavefront and hence the complex-valued sample transmission function of the object O can be reconstructed from the intensity distribution in the detection plane at some (defocus) distance z behind the object plane, based, for example, on the transport-of-intensity equation (TIE) or various iterative algorithms [1,6,7]. A wide spectrum of applications ranging from material science [8], paleontology [9], histology [10], small arthropod biology [11] has emerged, taking advantage of the unprecedented 3D imaging capability, including fast tomography applications of dynamic processes [12]. While this development has been driven by the brilliance and coherence of synchrotron radiation, to some extent the method can be translated to compact laboratory-scale µ-CT instruments, preserving fully quantitative reconstruction [13]. Furthermore, in order to achieve higher spatial resolution, propagation imaging and tomography can be extended from the parallel beam case to quasi-spherical beams, allowing for adjustable geometric magnification [14][15][16]. With this approach, projection micrographs and tomograms of biological cells [16,17], tissues and small organisms [18] can be recorded with quantitative density contrast and pixel size below 50 nm.
The necessary diverging spherical wavefronts for high magnification are created by nanofocus x-ray optics, such as Kirkpatrick-Baez (KB) mirror systems [19,20]. For such optical systems, already small imperfections like mirror figure errors result in strong phase aberrations due to the small wavelength of x-rays [21]. The empty beam far-field patterns are therefore strongly affected and never correspond to the clean homogeneous flat field distribution ideally assumed. This aberrated empty beam, however, constitutes the probing field P for O. According to current wisdom, the problem is dealt with by simple division of the recorded image with object by that of the empty beam [17,22]. In fact, the raw images are hardly ever presented in publications. Hence, significant errors can be expected in the subsequent phase retrieval process. While the simple flat field correction scheme by division is actually a crude approximation, it still works in some cases, whereas going terribly wrong in others. Figure 1 illustrates the basic problem investigated here. The inhomogeneities in P cause distortions which cannot be divided out. Strong artifacts of the mirror system are still observed in the 'corrected' hologram of a polystyrene sphere. The goal of the present work is to investigate the errors of the common flat field correction. In which cases are they strong, and in which cases are they negligible? In a second step, algorithmic approaches beyond the simple empty beam correction are presented, and tested numerically as well as experimentally.
Simultaneous reconstruction of probing beam P and object O have been achieved by ptychographic algorithms [23,24], based on the redundancy generated by lateral scanning with partial overlap. The method was first developed for far-field coherent diffraction data collected from a compact probe in the object plane [23,24], but some recent generalizations to nearfield (Fresnel) diffraction data have been demonstrated [25]. However, to this end P had to be further scrambled by a wavefront diffuser in order to generate enough diversity. Further, for the case of an illumination by a Fresnel zone plate, special reconstruction schemes were devised [26], based on constraints in several planes including the pupil of the zone plate to define the probe's support. Unfortunately, both ptychography based on lateral scanning and the multiplane Fresnel-CDI scheme [26,27] turn out to be unsuited for the problem of the empty beam reconstruction in full field propagation imaging, since the necessary redundancy/diversity is generally not provided by lateral scanning. For a discussion of redundancy/diversity as well as of the sampling criteria relevant for ptychography, see also [28]. This is why we turn to an approach based on multiple detection planes in this work in order to reconstruct P. A multi-plane reconstruction scheme was first introduced in optics and electron optics in 1973 [29]. Since then this scheme and variations thereof have been used in microscopy to recover the complex valued point-spread-function (PSF) [30], and also in the field of x-ray imaging to reconstruct an exit wavefield P · O [31,32]. However, in [31,32] the multi-plane detection has not been used to recover P itself, simultaneously with the object O, which is provided here. We therefore extend near-field imaging in a similar manner as ptychography has previously done for far-field coherent imaging [23,24]. The manuscript is organized as follows: Section 2 further defines the problem of the wavefront distortions and the associated correction schemes. The notation and the basic concepts of the work are introduced. Section 3 presents the experimental setup and parameters. Section 4 presents a single distance approach for the reconstruction of P, which we denote as divide&update. Finally, we turn to the reconstruction of P based on multiple detection planes in Section 5, before the paper closes with conclusions and outlook in Section 6.

How wrong is the standard empty beam correction?
Before we turn to reconstructions of P let us address the validity of the standard empty beam correction, i.e. the approximation  where D z denotes the Fresnel propagation over the propagation distance z, with coordinates defined in Fig. 2. D z is given by where F denotes the Fourier transform, λ the wavelength and k x,y the spatial frequencies.
The correction according to Eq. (1) can be justified for point sources, or equivalently plane waves, see also appendix A. Only in idealized cases, the measured |D z (P · O)| 2 can be written as |D z (P)| 2 |D z (O)| 2 . In all other cases in particular for P resulting from focusing the x-rays to a finite spot size, |D z (P · O)| 2 corresponds to a convolution of P and O, and cannot be 'deconvoluted' by applying Eq. (1). Note that propagation of parallel and spherical wavefronts can be treated in the same way, after application of a simple variable transform introducing an effective defocus distance z → z eff = z 1 z 2 /(z 1 + z 2 ) as well as a geometric magnification of the pixel size M = (z 1 + z 2 )/z 1 (Fresnel scaling theorem) [2]. For numerical calculation, the divergent (cone) beam case can therefore be easily transformed to an equivalent parallel beam case, see Fig. 2(b). By rearranging Eq. (1) we can introduce a metric δ to measure the error of the flat field correction  Since δ yields a matrix of errors, we use for our purposesδ for the mean error over all pixels in the image in the simulations below:δ where N is the number of pixels in one direction. In order to study the dependence ofδ on the geometric parameters, we perform a numerical simulation. The relevant parameters are the source size a generating the probe P and the characteristic length of the object ∆x O . To be able to work in the (numerically convenient) equivalent geometry of parallel beams, a correlation length of the probe ∆x P is introduced, characterizing the smallest length scale on which P shows variations. In other words ∆x P quantifies the smoothness of the image representing P cf. Fig. 2(c). Next, in view of the physical interpretation, we can relate ∆x P to a source size a, since the correlation length is limited by diffraction In the simulation we prescribe different ∆x P by a Fourier filtering of P according to . Different filter kernels can be used, such as Gaussian, hat profile, etc. In fact the simulation results below are found to be very similar. Finally, we chose a rectangular filter with the cutoff parameter ξ as the control variable. ξ is related to the effective source size (in the divergent beam geometry) according to a(ξ ) = z 1 λ ξ 2π . Stated differently, the filter is used to blur the probe P(x, y), in order to make it appear more like a plane wave on local scales. With ∆x P as the second control parameter, the errorδ can be studied for any chosen example of P and objects O of varied feature size in pixel units, such as a sine grating with periodicity ∆x O . Note that the filter kernel is not introduced to take the effect of partial coherence into account, which are of course relevant but not considered here. We always assume a fixed phase relationship between different points on the wavefront, i.e. a spatially fully coherent but aberrated probe.
In the simulation, we then evaluateδ for the probe phantom shown in Fig. 2(c), denoted as the mandrill-Dürer-probe, where the amplitudes of P are given by the Dürer carving, and the phases by a mandrill monkey. The test object is chosen as a purely phase shifting object consisting of a sine grating with a grid constant (periodicity) ∆x O . Both probe and object are represented by a N × N matrices with N = 512. The error (of the conventional empty beam correction according to Eq. (4) is then evaluated for a constant Fresnel number F = ∆x 2 /(λ z eff ) = 0.0015, as a function of ∆x O in the range from 16 px to 128 px and as a function of the filter cutoff ξ , which was varied in the range from 2π/(N∆x) to π/∆x, with ∆x the pixel size in the plane of the object. Note that this formulation in terms of the natural units of pixel units and Fresnel number has been chosen in view of the generality. The value of F corresponds in particular to the following (experimental) parameters: λ = 0.07 nm, ∆x = 25 nm, z 1 = 6 mm, z 2 = 519 mm and z eff = 5.93 mm. Next, we note that the varied range of the filter cutoff parameter can be transformed to a range of a from 0.34 px to 349.92 px, according to Eq. (5), evaluated for the above (experimental) parameters.
The expectation is thatδ increases when either the source size a is increased, or the feature size ∆x O is decreased. Figure 3 shows the simulation results forδ for different ∆x O as function of ξ /ξ max , with ξ max = π/∆x. The alternative x-axis label on the top indicates the corresponding a. As expectedδ decreases if the bandwidth decreases, i.e. if P is filtered more strongly with a smaller cut-off frequency ξ . Furthermore, the errorδ is always larger for smaller grid parameters than for larger grid parameters. The smaller the structure, the more stringent is the requirement to keep the source size a small or the correlation length of the wavefront aberrations ∆x P large. From this we can infer that small structures in the object will not be transferred correctly after application of the empty beam division, while the contrast of larger objects may not suffer so much. Interestingly, the experimental example shown in Fig. 1 seems to be in contradiction to the simulation results. The hologram of a 10 µm sphere which is illuminated by a beam emanating from a KB-focus with 100 nm spot size should not be affected by artifacts due to empty beam division. The paradox is resolved by the fact that while the full width at half maximum (FWHM) of the focal spot is well below 100 nm, the KB focus is characterized by long tails which extend over several µm. It is these tails which render the effective source size too large for a 'clean' empty beam division. This conclusion has been validated experimentally by placing a pinhole with 1.4 µm diameter in the focal plane, which then results in a less disturbed, i.e. cleaner empty beam, see Fig. 4.

Experimental setup
The experiment was performed at the nano-focusing beamline ID22NI at ESRF (Grenoble) at a photon energy of 17 keV (0.0729 nm). The source size (FWHM of the focus of the KB-mirrorsystem [21]) was determined to be 67 nm(vertical) and 73 nm(horizontal), as measured by a fluorescence knife edge scan. The sample and detector stage can be moved within a defocus range of 0 − 300 mm and 300 − 1500 mm, respectively. A FReLoN HD 4m detector (ESRF) equipped with a LSO:Tb (20 µm) scintillator and an optical magnification scheme was used, with an effective pixel size of 756 nm. In the first configuration the detector was placed at z 1 +z 2 = 525 mm. The variation of z eff is achieved by placing the sample at different z 1 . Figure 4 shows the measured intensity patterns recorded in this setup. The left column corresponds to the empty beam, the center to the raw hologram of a grid and the right column to the empty beam corrected hologram, for the data recorded without (top row) and with (bottom row) a 'cleaning' pinhole of 1.4 µm diameter positioned in the focal plane. P and P · O were recorded quickly after each other, so that the flat field correction is not affected by temporal drifts. In the second configuration the detector is moved and the sample (gold grid) is kept at a constant z 1 (25 mm). This approach also yields different z eff but the relative variation in z eff is much smaller if only z 2 is changed. The detector was placed at z 2 = {500, 501, 510, 600} mm. For each distance 100 images with and without O were recorded over a total of 40 minutes. Figure 5 shows data recorded using this configuration. Here, the time span between the measurements of P and P · O was 20 minutes. This delay has visible influence on the quality of the flat field correction (c). Thus the primary prerequisite for the empty beam division is of course a stationary illumination field, or a correspondingly small sampling interval, in order to avoid temporal drift. While this trivial reason for possible artifacts is not our main concern, it does affect some of the experimental data. For this reason phase reconstructions run on data as Fig. 5(c) yield poor results as depicted in Fig. 5(d). The multiple detection plane reconstruction (Section 5) suffers also from this issue. For comparison, using data, collected over a considerably smaller time interval, as shown in Fig. 4(c) to run the phase reconstruction yields Fig. 7(c), as discussed further below. All relevant experimental parameters for the measurements shown are summarized in Tab. 1.

Reconstructions from a single distance
In order to overcome the errors involved with the conventional empty beam division, as illustrated above, we now turn to reconstruction strategies of P and O. To this end, we use the common class of projection algorithms for the reconstruction [5], such as Gerchberg-Saxton (GS), Hybrid-Input-Output (HIO), Relaxed Averaged Alternate Reflection (RAAR) [33] etc. The algorithms were implemented using the SciPAL library [34] and are based on two alternating projections in two different planes: the projection P M adapts the current iterate of the wave where I denotes the measured intensities. The second projection P S acts in the plane of O and can imply different constraints on O. For example positivity of the electron density variation, meaning O implies negative phase shifts φ with respect to air/vacuum. P S could thus keep the negative φ and set all positive φ values of the iterate to zero. It can be formulated as where x denotes the coordinates in the plane perpendicular to the propagation axis and A the (reconstructed) amplitude at a given x. Combinations of different constraints are also possible for example negative phase as well as a support constraint [17]. However, none of these constraints is sufficient for reconstruction of P from a single intensity measurement, since (i) P fills the whole field of view so that no compact support can be used, (ii) amplitudes and phases of P are in the general case not proportional (coupled) or zero in any plane, and (iii) positivity may apply to the electron density which determines O, but figure errors of the mirrors can induce both positive and negative phases. Thus more information in form of additional measurements has to be used. The first idea which comes to mind is to bring a well known O into the beam and take images with and without O. To this end, in the experiment a 200 nm gold grid with 3 µm period and 50% duty cycle (referred to as grid below) was used, see Figs. 4(a) and 4(b) for the measurements used in the reconstruction discussed further below. The basic principle of the algorithm for coupled reconstruction of O and P is sketched in Fig. 6(a). We will denote it as the divide&update (d&u) algorithm. The right side of the sketch shows P M (Ψ P·O n ) and P M (Ψ P n ) where n denotes the iteration index. The left side shows the update of P and O. In order to obtain a separation of the exit wave field P · O, we need an initial guess for O. A guess for O is easily obtained by running a phase reconstruction step on the empty beam corrected measurements. One iteration of d&u starts with applying P M (Ψ P·O n ) and P M (Ψ P n ) on P and O. P 0 and (P · O) 0 (the initializations) are chosen with uniform amplitude 1 and phase 0. After application of P M (Ψ P·O n ) we can use the guess of O to obtain an update for P (blue arrow). We likewise use P M (Ψ P n ) to update O (red arrow). The updated versions O n+1 and P n+1 are then used to compute the new exit wave (P · O) n+1 (green arrows). This concludes one iteration of d&u. In practice, a few enhancements were added to the algorithm, and several additional constraints were tested. First, in the spirit of hybrid input-output, the update P n+1 is calculated as an average of P n and the update, same holds for O n+1 . Second, O n+1 and P n+1 were always projected on negative phases, see Eq. (8). We can consider this negativity (or positivity) constraint to be based on a physical relationship, such as positivity of the electron density (implying negativity of the object phase with respect to the vacuum wave). Therefore, one may be tempted to conclude that it would not hold for scattering length distributions which can change sign (as for neutrons), or for a sample embedded in a matrix. In the latter case the density of the object could be larger or lower than the matrix. Since a global phase offset is irrelevant, translating or reflecting all pixel values to negative or positive sign thus is not restrictive in our case, but is found to stabilize the algorithms by 'aligning' the phases. Finally, we projected O to unit amplitude 1, implying the constraint for a purely phase shifting O. to the initial guess (d), which was computed from the standard empty beam corrected hologram by a GS-algorithm. The reconstruction in the d&u scheme is clearly better then the GS-solution which is spoiled by the wrongful empty beam division. However, the d&u result still shows residual cloudy artifacts, which show that the separation of O and P is not sufficient. If we use a purely phase shifting P in the simulation and imply the corresponding constraint, the reconstruction quality is excellent, see the reconstructions of P in (f), and of O in (g), respectively. The amplitudes of P are not shown but exhibit ±15% (max.) variations around 1. This variation is due to the fact, that the final reconstruction is generated after a last application of P M , we use this procedure for every reconstruction. Figure 7 shows the d&u reconstruction of the experimental data recorded for the gold grid. Close inspection of the reconstructed phases of P (a) show again regular structures which are artifacts stemming from the gold grid. The amplitudes of P (b) look reasonable, with no bigger imprints of the grid. Here, the reconstruction quality cannot be improved significantly by assuming a constraint of coupled real and imaginary component of the index of refraction (dipersion and absorption correction), such as proposed in [35]. The results also do not improve if we take a second z 1 position for the same object into account such that two positions in the same P are probed.

Reconstructions from multiple detector distances
The previous section has shown that images of the object and probe at a single distance do not fully constrain the reconstruction of P. We thus turn to an approach of using multiple measurements of the beam in parallel. By moving the detector coordinate z 2 , the effective propagation distance z eff is modified both for P·O and the empty beam P, which is measured at each detector distance. This experimental setup is sketched in Fig. 8. For the geometry of high magnification, z 2 z 1 , z eff changes only very little with the variation of z 2 , since for z 2 → ∞, z eff → z 1 . We use the multi-distance (multi-magnitude) projection algorithm proposed by Allen et al. [36] to reconstruct the wave fields. The algorithm is denoted as mmp (multi-magnitude projection) and is rather straight forward. It uses only the measurements in the different detection planes as constraints, with no constraint in the sample plane. Using 4 detection distances, we have Ψ n+1 = P M 4 (P M 3 (P M 2 (P M 1 (Ψ n )))) , where P M i denotes the projection on the i-th measurement. We stress that no constraints on O are used in the multi-plane algorithm. The algorithm has been used before in the literature to reconstruct the exit wave in x-ray imaging [37]. Instead of the sequential scheme implemented here, one could alternatively also update via an average over all detection planes, i.e. in a parallel manner [32]. mmp was demonstrated to be robust against noise, and can deal with first order phase jumps. But to our knowledge the algorithm has not been used for the reconstruction of P. We mention two small deviations from the original algorithm (i) the measurements are not equally spaced (ii) P M brings Ψ always back in the plane of O before projecting on the next sample multiple detector positions Fig. 8. Detector distance variations. In this approach z 2 is varied and z 1 is kept constant. measurement. The original implementation projects directly on the next measurement.

Simulation
The simulated reconstruction was carried out for the same four detector distances as in the experimental case discussed below, but in the equivalent (parallel beam) geometry with z eff = {23.810, 23.812, 23.832, 24} mm. Poisson noise was added to each pixel of the simulated holograms of P and P · O for an expected number of 10 9 detector photons. For P the mandrill-Dürer-probe was used as above. For the phase of O, the image of a stone sculpture was used with phases in the range of [-0.4 0] rad, see Fig. 2(e). Amplitudes of O were set to 1 (purely phase shifting object). Figure 9 shows the simulation results for 1000 mmp iterations. The reconstructed phases and amplitudes of P, as shown in (a) and (b), respectively, are in very good agreement with the phantom, see Fig. 2 where N 2 is the number of pixels, and is given in the titles of the sub-figures for each reconstruction. We notice that ∆ is about 3 times larger for the phase reconstructions than for the amplitude reconstructions. Also the reconstruction of P · O has a higher error value than the one of P. The reconstruction results and error values are already quite good, despite the small changes in z eff (corresponding to the experimental values). As expected, the simulations show that by choosing propagation distances with more variation in z eff (translational diversity) the errors can be further decreased by several orders of magnitude. The challenge is therefore clearly to cover a sufficient z eff range in the experimental implementation.

Experimental data
Motivated by the successful numerical results, the next step was to apply the mmp algorithm to experimental data. The data was recorded in the second configuration as described in Section 3, exemplary data is depicted in Fig. 5. Due to the cone beam geometry and different magnifications the recorded images must be rescaled and aligned. The images were cropped to the smallest field of view. We have reconstructed the two data sets with and without O independently. Figure 10 shows the results of the reconstruction of experimental data. The reconstructions of phase and amplitude of P are shown in (a) and (b), respectively. The result looks similar but more plausible than the reconstruction of P in Figs. 7(a) and 7(b). Comparing the reconstruction of P and P · O in (c) and (d) we can also identify the same structures in both, which indicates that the reconstructions are consistent. However, when extracting O from the reconstruction of P · O by dividing by P, certain artifacts remain, and the result looks worse than the result for O reconstructed from a single distance as before in the d&u approach, shown in Fig. 7(d). This is attributed to the fact that this data set is more affected by temporal drift of P, simply because the delay between measuring P and P · O was much too long, namely about 20 min. It is surprising that the reconstruction of P is still possible at this point. For an optimized scan, a time interval between subsequent images of less than a minute is realistic, in view of the necessary detector movement. This could sufficiently avoid drift of P and enable reconstructions of similar quality as those in the numerical simulation.

Summary and Conclusion
The first aim of this work was to illustrate the problems associated with the standard flat field correction, both numerically and with experimental data. To this end, we have used an error metric δ Eq. (4). We have investigated, under which conditions the division of intensity images of P and P · O in the detection plane is a sufficient approximation, in terms of object and source size. From the simulation we can conclude that the reconstruction of the probe P is an important step towards higher resolution in propagation imaging, since otherwise good data quality can only be expected for feature sizes much larger than the source size.
As a first attempt, we have presented 2 algorithms to reconstruct the empty beam P in x-ray propagation imaging. The first algorithm d&u, devised for a single distance data set, uses a nested scheme to update P and O. This scheme works well for purely phase modulated P and O, but it breaks down for the general case of an amplitude and phase modulated P. We have also inferred from this approach that it is not sufficient to use one image of P and P · O each, even if O is known. The second algorithm mmp is already known in literature but has not been used to reconstruct P. The algorithm requires multiple distance measurements of P. mmp yields good reconstructions in simulations even for noisy data with small variation in z eff . In order to probe P at different distances in the experiment, the sample-detector distance z 2 was changed, which is intrinsically different from changing the sample position z 1 . Analyzing this data set showed that in order to get a good reconstruction of P the time between measurements should be short. This is of course also well known for the standard empty beam correction. However, in order to cover a sufficient diversity (range in Fresnel numbers), the detector has to be moved over a rather wide range, which makes this approach more prone to temporal drift. Flight tubes are another experimental difficulty in implementing this approach, unless the detector is contained in the flight tube. For the experimental strategy, the number of detector positions and the shift δ z 2 must be optimized in view of sufficient diversity while keeping the total measurement time as small as possible in view of drift elimination. Interestingly, a recent generalization of ptychography to longitudinal scanning (OAP, on-axis ptychography) has also shown, that a defocus scan with constant detector position, i.e. moving only the sample, is not sufficient for simultaneous reconstruction of P and O [38]. However, numerical simulations have yielded good results already for two detector distances. This still awaits experimental verification, but it can be anticipated that more than one algorithm will be suitable to perform the simultaneous which differs from the heuristic expression given in Eq. (4). The numerator of δ measures the difference of D z eff (P · O) and the product D z eff (O)D z eff (P). This difference is then normalized by the modulus of the propagated P and O. It is shown in [39] that the supremum of Eq. (11) is bounded by √ 2 z eff λ 8 π ∆x P ∆x O Figure 11 shows the results obtained with the metric Eq. (11). The simulation was carried out, following the same procedure as described in Section 2 also using the same set of parameters. Both errors Fig. 3 and Fig. 11 show the same behavior, regarding the strength of filtering (ξ /ξ max varied) and the variation of object periodicity ∆x O . In fact it is possible to define several other metrics, showing the same general behavior as here, but without the mathematical foundation as Eq. (11).