Coherent diffractive imaging: a new statistically regularized amplitude constraint

Statistical information about measurement errors is incorporated in an algorithm that reconstructs the image of an object from x-ray diffraction data. The distribution function of measurement errors is included directly into reconstruction processes using a statistically based amplitude constraint. The algorithm is tested using simulated and experimental data and is shown to yield high-quality reconstructions in the presence of noise. This approach can be generalized to incorporate experimentally determined measurement error functions into image reconstruction algorithms.


Introduction
Coherent diffractive imaging (CDI) is a method in which an object is illuminated using a coherent beam. The diffraction intensities recorded from the object are used to extract its complex transmission function [1]. The method shows considerable promise as a means to achieve very high spatial resolution imaging without the need for a high-resolution lens. As such, it is receiving considerable attention in the synchrotron [2]- [5] and x-ray freeelectron laser [6,7] communities, as well as for imaging applications using the diffraction of electrons [8,9].
The fundamental idea is to measure the diffracted intensity and to use iterative algorithms to find an object distribution that is consistent both with the measured intensity distribution and with a priori information that is known about the object, such as its physical extent, or support. Typically, iterative algorithms adapted from electron microscopy are used to obtain the solution [10]. If a solution that is consistent with the measured intensity and the known support is found, then it may be assumed to a high degree of certainty that the solution is the correct one [11].
The successful application of any iterative algorithm depends on the availability of an accurate model for the data acquisition process. Recent work has investigated issues associated with deviations from perfect spatial coherence [12]- [15] and has established that even small deviations from the ideal conditions implicit in the formulation of CDI algorithms can have significantly adverse effects on the quality of the reconstruction. A more common deviation from the ideal case arises because of the impact of noise or low photon counts on the diffraction data. Simulation studies typically investigate the impact of noise on the quality of reconstructions [16] and explore the limit when the iterative scheme fails to converge to a solution that is consistent with the simulated data. The impact of the low count statistics on the reconstruction of the phase of the complex-valued electron density of nanocrystals using Bragg CDI has also been pointed out recently [17]. The authors of [17] indicated, in a context related to the one considered here, that measured data points corresponding to low intensities may lead to difficulties in applying the modulus constraint in phase retrieval algorithms. There has been comparatively little work, however, on the development of approaches that include a knowledge of the statistics of the data within the reconstruction process. Indeed, while many 3 areas of modern image analysis incorporate the statistical properties of the data as an integral component of the methods, this area of development has not, to our knowledge, been applied consistently to the reconstruction of images from coherent diffraction data. The present paper represents a step in this direction.
In this paper, we present a new reconstruction approach involving a modification of the intensity constraint in the CDI algorithm that incorporates information about the distribution of intensity measurement errors.

Constraints and coherent diffractive imaging
The essential feature of the coherent diffractive imaging reconstruction algorithm is to determine a coherent exit wave leaving the sample that is consistent with both the measured far-field intensity distribution and the a priori information that is known about the object. A number of possible constraints exist for the plane of the sample, the most fundamental of which is some knowledge of its physical extent and shape, including the use of known illumination [18].
The present paper is concerned with the constraint that is applied using the experimental data. A coherent wave is described by both its amplitude and its phase. The iterative methods used to recover the exit wave distribution are often considered to be phase recovery methods because the measured intensity distribution provides the data from which the wave amplitude is extracted.
The method used to recover coherent x-ray images enforces consistency with the experimental data by imposing the measured amplitude after each iteration of the algorithm. Expressed symbolically, the amplitude of the complex wave (q) at scattering vector q is replaced at each iteration using the operation where l = 1, 2, . . . is the iteration number and A Est (q) is the estimate of the amplitude obtained from the experimental data. The key issue here concerns optimal estimates of the amplitude from imperfect measurements of the intensity.
In the next section, we explore how to estimate the wavefield amplitude, A Est (q), from experimentally observed intensities that include measurement errors.

Estimating amplitude from intensity measurements
In practice, the experimentally measured intensity is a mean value deduced from many independent measurements. For example, in the experimental data discussed in section 5, the far-field intensity distribution is formed by summing a large number of relatively short exposures. In this case, the measured intensity is more properly written as I (q) and the amplitude constraint, equation (1), is where the estimate of the wavefield amplitude is 4 In the limit of a large number of intensity measurements, I (q) is the expectation value of the intensity, and is given by where the probability density function, P I (q), is uniformly normalized, bounded for all I (q) and vanishes for I (q) → ∞ and I (q) < 0.
In this paper, we question the prevailing paradigm contained in equation (3) that assumes that equation (3) represents the best estimate of the wave amplitude from the measured intensity.
Here we seek a better estimate of the amplitude of the wavefield from the set of I (q), which we assume to be a continuous random variable. In this paper, we conjecture that the expected value of the wavefield amplitude,Â(q), measured at scattering vector q and derived from this set of measurements is better given by [19] For random variables X and Y, the Cauchy-Schwartz inequality states that Hence, equating X with the wave amplitude and Y with the identity operator, and noting that the wavefield and intensities are both intrinsically non-negative quantities, it follows that indicating that, given equation (5), the square root of the average intensity systematically overestimates the scattered wave amplitude in the presence of noise. It should be readily apparent that the inequality defined by equation (7) will only be significant for data with a poor signal-to-noise ratio, which will be the case for the high angle parts of the diffraction pattern and regions around intensity zeroes. Earlier work on the effect of partial coherence has shown that the regions around intensity zeroes can be of particular importance in the reconstruction process [15], an issue we discuss again in section 5. In the next section we estimate the significance of the distribution of errors on the estimate of the amplitude.

The probability distribution function
In this section, we analytically calculate the magnitude of the error in the amplitude estimate using a Gaussian probability distribution function. We use this distribution as an illustrative example and as a tool to explore the generic properties of this approach. We recognize that this will not always be the best model for the experimental statistics. We present some numerical results for other distributions at the end of the section.
Assume that the intensity at a given q has a mean value I (q) and a variance σ I (q) with a Gaussian probability distribution function that is therefore described for I (q) > 0 by The expectation value for the square root of the intensity is, therefore, described bŷ where I (q) = I (q) − I (q) . Clearly, equations (9) and (10) are valid only for I (q) > 0, and so we assume that any contributions to the integrand of equation (10) are strongly suppressed for I (q) < 0. This is consistent with prevailing practice in experimental CDI where negative intensities that appear after background subtraction are set to zero. Accordingly, the limits of integration have been extended to ±∞. We may form the expansion By symmetry, the odd powers vanish. Retaining the first two even terms, we find which is easily evaluated, using equation (3), to becomê If we define the signal-to-noise ratio for the measurement at scattering vector q as then we obtain the result for Gaussian statistics that where Note that the requirement for negligible negative intensity measurements broadly requires that S I (q) > 1. Figure 1 plots the curve E I (q) for the Gaussian distribution, along with numerical evaluations of equation (4) for a number of other distributions. Table 1 presents analytical forms for these curves. Across various error distribution functions, P I (q), E I (q) adopts the form where α and µ are real and positive constants as displayed in table 1.  1.199 Data set consisted of single-shot measurements or low-intensity level Exponential √ π 2 Rare random events a Calculated using a continuous interpolation of the discrete distribution and fitted using the y(x) = 1 − α · x µ function.  The form of equation (15) enables the amplitude constraint, equation (1), to be better written as This modified form of the wavefield amplitude constraint allows the incorporation of experimentally acquired knowledge about the distribution of measurement errors in the experimental data within the reconstruction processes.

Examples
This section provides examples of reconstruction from simulated and experimental data using the conventional amplitude constraint, equation (2), and the modified, statistically regularized, amplitude constraint, equation (18).

Numerical example
The modified wavefield amplitude constraint was tested on simulated data obtained from the object illustrated in figure 2 where N is the total number of scattered photons. To analyze the influence of measurement errors on the quality of reconstruction, several data sets for different N , model (a) 1 × 10 8 (dynamic range 1 : 1 × 10 6 ), model (b) 5 × 10 6 (dynamic range 1 : 6 × 10 4 ), model (c) 1 × 10 6 (dynamic range 1 : where I k N [m, n] is the intensity distribution at the kth data frame. To analyze the accuracy of two different estimates of the wavefield amplitude, the relative amplitudes,   (21) represents strong evidence that the initial hypothesis encapsulated by equation (5) is correct.
The object was reconstructed using a support function, Su, selected as a 128 × 128dimensional rectangle. The support was fixed during the reconstruction processes. Reconstruction results are presented in figure 3, where column (I) corresponds to reconstruction using equation (2) and column (II) corresponds to reconstruction using equation (18). The Poisson-type correction, E I (q), was adopted for the reconstruction in this case.
Although a reasonable quality image can still be obtained using CDI with a small level of noise, models (a) and (b), the quality of the reconstruction degrades with increasing noise level much more rapidly than in the case where the modified constraint is used.  (2) for the amplitude constraint, (II) reconstruction using equation (18) for the amplitude constraint. The red rectangle on images indicates the support assumed in the reconstruction.

Experimental application
The statistically regularized constraint is used here for image reconstruction from experimental data. The diffraction experiment was performed at the 2-ID-B beamline at the Advanced Photon Source. The CDI experiment used the sample shown as an inset of figure 4 and consisted of a set of apertures milled into a thick gold substrate, which is opaque for 1.4 keV x-rays. The sample was placed 11 m from the 5 µm exit slit of the monochromator and 1.4 m from the detector. As the sample is primarily opaque, the contribution by the undiffracted beam is negligible; hence the diffraction data can be collected without a beam stop. The diffraction data were obtained using a charge-coupled device (CCD) camera containing a 2048 × 2048 array of 13.5 µm pixels. The data set consisted of 600 frames. Each data frame was corrected for dark field. Then the average intensity distribution, I [m, n] , was calculated using where M is the number of data frames and I k [m, n] is the intensity distribution at the kth data frame. Figure 4 shows the resulting diffraction pattern. Consequently, the variance for each pixel in the data was calculated using (23) Figure 5 shows the normalized intensity distribution function for three selected data points with different average intensities. The center 1024 × 1024-dimensional area of the discretely sampled diffraction pattern was selected for reconstruction. The support function was chosen to be a 342 × 342-dimensional rectangle, fixed during the reconstruction processes. Experimental data were analyzed using two estimates for wavefield amplitudes, equations (2) and (18). The Gaussian-type probability distribution function was used for the statistically regularized constraint, for whicĥ Reconstruction results are presented in figure 6. The conventional CDI approach, which assumes fully coherent illumination of the object and error-free diffraction data, is able to reconstruct an image of the object but displays non-uniform distribution of the magnitude of the reconstructed wavefield, and a significant number of spurious artifacts, figure 6(a). Figure 6(b) shows the reconstruction using the modified constraint. One can see that the resulting image displays a more uniform distribution of the magnitude of the reconstructed wavefield than is obtained using the conventional CDI algorithm. Note that an analysis of these data has also been presented elsewhere [15]. Interestingly, comparison with the results in this paper demonstrates that the correction factor applied here can serve to moderate the impact of partial coherence by suppressing the effect of noise around intensity zeroes [14], but at the cost of reduced resolution and a degree of distortions of the features in the reconstruction-note, for example, that the circular feature is not reconstructed with a great deal of fidelity.   The current algorithm's advantage, in addition to those obtained in [15], is that it is able to accommodate the effects of measurement noise, even if the illumination is fully coherent. The algorithm described in [15] is designed primarily to accommodate partial coherence in data that are largely free of measurement noise. A combination of both approaches can be envisaged in which the reconstruction of noisy data in the presence of partial coherence is improved. This will be the subject of further work.

Conclusion
In this paper, we consider the error-based CDI constraint for image reconstruction from noisy data. The key idea of this method is to use statistical information about the uncertainties in intensity measurements in the reconstruction process. The results presented in section 5 demonstrate that the modified CDI algorithm is able to reconstruct images more successfully from noisy data than the conventional CDI approach. We do not claim that the prescription of equation (5) is optimal; however, we find that we obtained superior results with both simulated and real experimental data using it.