An iterative method for robust in-line phase contrast imaging

We present an iterative near-field in-line phase contrast method that allows quantitative determination of the thickness of an object given the refractive index of the sample material. The iterative method allows for quantitative phase contrast imaging in regimes where the contrast transfer function (CTF) and transport of intensity equation (TIE) cannot be applied. Further, the nature of the iterative scheme offers more flexibility and potentially allows more high-resolution image reconstructions when compared to TIE method and less artefacts when compared to the CTF method. While, not addressed here, extension of our approach in future work to broadband illumination will also be straightforward as the wavelength dependence of the refractive index of an object can be readily incorporated into the iterative approach.


Introduction
Propagation based x-ray phase contrast imaging exploits the change in phase that occurs when x-rays propagate through a material [1,2]. It allows for improved imaging of objects that exhibit poor x-ray absorption contrast, or which are sensitive to x-ray damage.
Various analytical techniques-such as transport of intensity (TIE) [3] and contrast transfer function (CTF) [4,5] methods-have been developed for propagation-based x-ray phase contrast imaging. These methods rely on different approximations [1,2,5] and show different susceptibility to noise [6,7] and to other sources of imaging artefacts. Typically all these approaches are applied in the Fresnel regime, which can be referred to as the near field consistent in this context of phase contrast imaging.
In recent years, approaches described as coherent diffractive imaging (CDI) [8] have been successful in recovering the complex transmission function of a sample through iterative methods that rely on oversampling in the far field [9]. A common algorithm often used in CDI is Fienup's error reduction (ER) algorithm [10]. A range of iterative methods have allowed for the inclusion of complex constraints [11], variable supports through the iteration sequence [12] and models for the partial coherence of the illuminating wavefield [13]. The effects of significant wavefront curvature have also been incorporated in, for example, Fresnel CDI (FCDI) [14].
Following the success of far-field iterative methods, iterative near-field methods have also been explored where the Fresnel approximation now applies due to curvature in the propagated wavefield rather than purely from the incident wavefield as in the case of FCDI [15][16][17][18][19][20][21].
In this work, we demonstrate the iterative implementation of an analytical phase retrieval approach that was developed for use with compositionally uniform objects [22]. As with other near-field or Fresnel iterative approaches, the different approximations involved offer improved robustness to noise and artefacts compared to the analytical methods. The uniform material approximation also offers an advantage in that it readily allows the thickness of an object to be recovered, which is a quantity that is often straightforward to validate.

Algorithm
The standard ER algorithm [10] can be written in operator notation as where the ith iterate of the exit surface wavefield of the sample, y, is defined as the product of the incident wavefield, inc y , and the transmission function, T . The modulus constraint is defined as where P and P 1 represent, separately, the forward and reverse propagation operators of the paraxial wavefield.
The modulus constraint, m p , operating on the detector plane wavefield, ŷ , is defined as where I is the measured intensity. The support constraint is where r is a two-dimensional position vector in the plane of the wavefield perpendicular to the direction of propagation and S is the area of the support, which contains the sample.
In the numerical implementation of the algorithm the wavefield must be oversampled in the detector plane [9] and the algorithm may be commenced with a guess for 0 y .
Convergence of the algorithm is assessed by evaluating the quantity where the sum is performed over all N pixels in the numerically sampled quantities.
To obtain information about the sample, the final iterate of the wavefield can be divided by the known or recovered [23] incident wavefield to obtain the complex transmission function of the sample where the amplitude, A, and phase, f, can be written in terms of the complex refractive index, n 1 i d b = -+ . The real part of the refractive index, δ, defines the phase shift that occurs, while the magnitude of imaginary component, β, determines the absorption.
It can be seen that T can be expressed as a linear projection through the sample. Therefore if sufficient angular projections through the sample can be recovered the three-dimensional distribution of the complex refractive index can be recovered using standard tomographic methods [1,2].
In the case of a uniform sample of known composition the recovered transmission function can be used to calculate the thickness of the sample directly via This can be incorporated into the sample plane support constraint so that the current estimate of the sample thickness is used to calculate the next iterate of the transmission function and the wavefield so that information from both the phase and amplitude are combined [11]. Accordingly, in our implementation the sample plane support constraint may be written as In simulation tests of the algorithm the error metric for the known sample can also be monitored via  Figure 1 shows the results of a numerical simulation that demonstrates the effectiveness of the algorithm in obtaining a quantitative reconstruction of an object in the near field at 9 keV. Our test object, shown in figure 1(a), has a uniform composition with 2.33 10 6 b =´and 2.69 10 5 d =´-, and has thickness from 0 (black) to 1 μm (white) according to the linear grayscale in figure 1(a). A point source 0.76 m from the sample and a sample-to-detector distance of 0.24 m were used with a pixel size of 342 nm in the sample plane resulting in a detector pixel size of 450 nm in the 710×710 pixel image shown in figure 1(b). The Fresnel propagator in its convolution form [1,24] was used to propagate between the sample and detector planes. It can be seen in figures 1(c) and (d) that the algorithm converges smoothly to an accurate estimate for the thickness.
For the iterative method to work, the sampling requirements for the relevant propagator must be satisfied [24]. For a uniform sample those requirements are set out in table 1. In the case of the TIE and CTF analytic approaches, the influence of the sample transmission function-which, in the projection approximation [25], is readily related to the sample thickness-can also be incorporated into the validity conditions, as set out in table 1.
The validity conditions for the TIE and CTF are defined by the need to linearise the relevant expressions for the propagated wavefield [5]. The iterative method instead requires that the relevant combination of the exit wavefield from the sample and the propagator term is properly sampled. This is most easily written in the case of the single fast Fourier transform (FFT) method where the requirement is For the convolution and angular spectrum methods the expression is a little more complicated, but a good rule of thumb from our simulation work is that provided both the propagator term and phase of the sample transmission function are well sampled-that is D is the pixel size-then the iterative approach will provide an accurate answer. Accordingly, due to the different validity conditions, samples that do not satisfy the phase gradient conditions required for TIE and CTF may still be recovered using an iterative approach.
Using the same simulation as the previous example, figure 2 shows a demonstration of the robustness of the iterative method when compared to the CTF (with a Tikhonov regularisation [26] of 0.01 a = ) and TIE. In this geometry the phase gradient validity condition for the CTF is violated. It can be seen that only the iterative approach provides a high fidelity reconstruction of the input object.  (7) and (9); and (d) a line out (blue line in figure 1(a)) comparing the object (blue) and the reconstruction (green)-an image of the iterative reconstruction is shown in figure 2(a). x y r , = ( ), z eff is the effective propagation distance, λ is the wavelength and other variables are as defined previously.

Technique
Validity condition Reference Experiment A proof-of-concept experiment was undertaken at diamond light source's B16 Test Beamline [27]. A conventional in-line phase imaging geometry was used with an energy of 9 keV. A point source was created through the use of a Kirkpatrick-Baez mirror, where the horizontal and vertical full width half maximum (FWHM) of the focus was 0.837 μm and 0.607 μm respectively. The source to sample distance, z 1 , was 0.76 m and the sample to detector distance, z 2 , was 0.24 m. The area imaged contained an object that was constructed from Au deposited on a Si N 3 4 window that was larger than 700 μm×700 μm. The design of this object can be seen in figure 3(a) and was fabricated via direct current sputtering. The object consisted of two overlaying deposits, with the deposits having a thickness of 350 or 650 nm (blue and green shades) respectively. The regions where the deposits overlapped (red) resulted in a thickness of 1000 nm. The sample was contained within an area of 105 μm×105 μm of the Si N 3 4 window. A cooled 14-bit charge-coupled device camera was used in conjunction with a 20×objective lens to image the distribution of optical photons produced in a YAG:Ce 35 μm thick scintillator. The resulting detector pixel size was 450 nm. As a result, the effective pixel size at the sample plane was 342 nm. A total of 62 sample phasecontrast and sample-free white-field images were acquired. The exposure time for each of these images was 75 s.
Under these conditions the experimental geometry satisfies the validity conditions for the iterative and TIE methods, but in some regions violated the phase conditions for the CTF. From simulations the CTF violation was not expected to cause significant artefacts. The projection approximation is also satisfied.

Results
Variations in the illumination and non-uniformity in the scintillator were removed by dividing the sample phase contrast image with the white-field image taken with no sample present. The field of view used in the analysis is shown in figure 3(b). The 710×710 pixel area resulted in an oversampling factor equal to four. Despite the background normalisation step, the image contained noticeable artefacts (seen in figure 3(b)) likely due to beam instability over the acquisition time.  The β and δ values required for our algorithm were calculated from the relevant scattering form factors, obtained from Henke et al [28]. To obtain the correct sample density, the average intensities transmitted through the sample for various regions were calculated and used in conjunction with Beer's Law. While the phase contrast redistributes the intensity in the form of some fringing, an average over the whole area of the sample that includes all fringes should include all the transmitted power and give an accurate measure of the average intensity. We found the density of the deposited gold to be 13.8±2.0 g cm −3 compared to bulk density of 19.3 g cm −3 . This is consistent with previous work using similarly deposited samples [29,30].
Using the Fresnel propagator in its convolution form, 1000 iterations of our algorithm were applied to the phase contrast image in order to test the quality of the solution and its convergence properties. During this process a very loose square support was applied. Using a shrink-wrap support [12] gives results that are not significantly different for the sample and in the results shown we have instead used a loose support in order to better compare regions outside the sample with those recovered using the TIE and CTF algorithms. The error metric stabilised around the 500th iteration. Figure 3(c) shows the reconstruction of the sample thickness using our iterative algorithm.
We compared our thickness reconstruction to that of both the TIE and CTF as shown in figures 4 and 5. A Tikhonov regularisation [26] was applied to the CTF with 0.01 a = . The mean of the region outside the sample was used to define zero sample thickness for the CTF and TIE, as these methods include an arbitrary constant phase offset in the result.
Qualitatively, all three algorithms show similar non-uniformities in the sample reconstruction for regions that are of uniform thickness most of which are attributable to instabilities in the illumination and consequent non-uniform flat-fielding of the data. Assuming the manufactured regions are of constant height (see figure 6) standard deviation of the recovered thickness in each of the regions is shown in figure 4 is taken as an estimate of the error due to the reconstruction process. The resulting thickness and uncertainty for each method is shown in the table in figure 4.
The residual values inside the loose support but outside the sample provides an additional analysis of the estimated error. The standard deviation of these values, which was 46 nm, may be attributed to the level of uncertainty in the measured intensities at the detector due to sources other than the sample, such as scatter from components in the imaging system, beam variation and noise in the detector. It can be seen that this is a relatively large part of the overall uncertainty and that taking care to reduce these sources of error may allow more quantitatively accurate results to be obtained.  While the iterative and CTF methods offer better spatial resolution than the TIE approach, they also introduce a cross-hatching image artefact. This artefact is consistent with the loss of phase contrast information in the measured image due to effects such as finite source size, noise and detector resolution as demonstrated in figure 7. In our case the scintillator system employed means that the detector point spread function will be the limiting factor. In addition, the sampling requirements may be affected by issues such as detector resolution and, consequently, the choice of algorithm for obtaining the desired information will depend on the interplay of system parameters and the actual sample.

Conclusion
We have experimentally demonstrated an iterative near-field in-line phase contrast method that allows quantitative determination of the thickness of an object when monochromatic x-ray illumination is used. The image reconstruction and quantitative results were validated against existing analytic phase retrieval methods. The iterative approach can be used under conditions where the analytic methods are not valid and does not need to be adjusted for an arbitrary constant phase offset.
For this demonstration we employed hard x-ray illumination and isolated single material objects. However, the algorithm should be viable at other wavelengths and while it is also compatible with ptychographic or other scanning techniques it does offer a way for single-shot recovery when multiple exposure approaches are not practical.
Extension of our approach to broadband illumination will also be straightforward as the wavelength dependence of the refractive index of an object can be readily incorporated into the iterative approach.