Sparse recovery of undersampled intensity patterns for coherent diffraction imaging at high X-ray energies

Coherent X-ray photons with energies higher than 50 keV offer new possibilities for imaging nanoscale lattice distortions in bulk crystalline materials using Bragg peak phase retrieval methods. However, the compression of reciprocal space at high energies typically results in poorly resolved fringes on an area detector, rendering the diffraction data unsuitable for the three-dimensional reconstruction of compact crystals. To address this problem, we propose a method by which to recover fine fringe detail in the scattered intensity. This recovery is achieved in two steps: multiple undersampled measurements are made by in-plane sub-pixel motion of the area detector, then this data set is passed to a sparsity-based numerical solver that recovers fringe detail suitable for standard Bragg coherent diffraction imaging (BCDI) reconstruction methods of compact single crystals. The key insight of this paper is that sparsity in a BCDI data set can be enforced by recognising that the signal in the detector, though poorly resolved, is band-limited. This requires fewer in-plane detector translations for complete signal recovery, while adhering to information theory limits. We use simulated BCDI data sets to demonstrate the approach, outline our sparse recovery strategy, and comment on future opportunities.

In a BCDI experiment, increasing the beam energy from a value currently suitable for CDI (about 9 keV) to one necessary for HEDM (≥50 keV) compresses reciprocal space by at least 80% in each dimension. Under these conditions, the solid angle subtended by each detector pixel covers a substantial portion of the diffraction pattern. As a result signal features (diffraction fringes) critical to inversion of the diffraction pattern are no longer resolved by a typical pixelated X-ray area detector. This is demonstrated in Fig. 1(b) which shows a coherent diffraction intensity pattern depicted with well-resolved and under-resolved fringes, mimicking measurements made at lower and higher X-ray energies respectively. The subsequent loss of diffraction feature visibility at higher energies renders conventional phase retrieval ineffective.
For the recovery of fine feature detail in the intensity pattern, additional independent measurements (upsamples) of the wave field intensity are therefore necessary. In this work, we consider the case where these additional measurements come from translating the area detector perpendicular to the exit beam in sub-pixel steps, as has been explored previously 9 . Each pixel measurement in such a data set imposes a constraint on a particular region of the scattered intensity (Fig. 2). In theory, the fine detail at a desired sub-pixel resolution could be recovered if one were to acquire an adequate number of these constraints, taken from sufficiently small detector offsets. The difficulties in fulfilling this requirement are described in detail in Section 1 of the Supplementary Material. In this paper we utilise compressed sensing 10,11 to demonstrate signal recovery with significantly fewer of such pixel constraints than there are fine pixels in the desired image. The required number of constraints is dictated by the information content when the desired image is expressed in Fourier space. The reduced number of constraints (or measurements, in compressed sensing parlance) is made possible by the fact that the BCDI intensity pattern of a compact single crystal is necessarily band-limited. In the compressed sensing approach employed here, one needs to make approximately as many coarse pixel measurements of the poorly-resolved signal as there are non-zero components in the Fourier representation of the well-resolved signal. The specific representation used in our method is the cosine basis. Recovery of this condensed representation of the signal is perfectly suited Three-dimensional reciprocal space is sampled in finite steps as shown. Two of the three steps are in the detector plane (one of which is perpendicular to the figure plane, denoted by '°') and the third dictates the migration of the plane. This set of sampling vectors is fixed by the crystallographic orientation of the scatterer with respect to the beam, and the manner of its rotation. (b) Acquisition of more data by detector translations perpendicular to k f in case of coarse detector resolution. Application of sparse recovery techniques greatly reduces the number of additional measurements required. to mathematical optimization algorithms that specialize in sparse arrays. Throughout this paper, we refer to the upsampled intensity patterns resulting from this sparse recovery scheme as the 'recovered' images or patterns. This paper is organised in the following manner: Beginning with simplifying assumptions, we describe our BCDI simulations at a range of energies from 9 keV through to 54 keV, the latter corresponding to a high-energy experiment conceivable at a high-brightness synchrotron. This includes establishing a suitable ground truth intensity pattern with well-resolved intensity fringes consistent with currently feasible measurements, against which all recovered intensity patterns are benchmarked. We describe the detector sub-pixel translation method and the additional constraints provided by the intrinsic band-limit of the intensity pattern, followed by the compressed sensing algorithm employed to recover each slice of the 3D data set independently. We demonstrate the recovery of qualitatively different intensity patterns, with simulations that represent different beam energies and varying sub-pixel detector translations. Finally, we show three-dimensional reconstructions of a simulated nanocrystal obtained by applying conventional BCDI phase retrieval to simulated high-energy data sets recovered with our approach. We describe how our method can emulate a variety of physical modifications to a BCDI experiment that may be difficult to realise experimentally. We close with comments on new experiments that this method can potentially enable at next-generation synchrotrons.

High-energy BCDI Simulations
BCDI measurements query the wave intensity in the Fraunhofer regime 12 in which the scattered wave front is the Fourier transform (FT) of the illuminated compact crystal. As the scatterer is rotated through the Bragg condition the intensity pattern is measured, slice by slice, on the area detector. For a crystal of size ~300 nm illuminated with hard X-rays of energy ~9 keV, the required range of this rotation is sufficiently small (about 0.7°) that the successive slices are approximated as parallel in reciprocal space 4 .
We simulate the compact crystal as the set of grid points inside a faceted volume, at the centre of a three-dimensional complex array. The interior points are complex-valued with magnitude 1 and a spatially varying phase that mimics a strain field. The exterior points have magnitude 0. The corresponding diffraction signal is obtained via the three-dimensional FT. In the 3D array of this FT, two dimensions represent pixel coordinates of the area detector while the third represents successive images acquired by rotating the crystal through the Bragg condition. In our simulations, the crystal resides inside a 22 × 24 × 22 box within the simulation array of size 128 × 128 × 70, denoting 70 angular steps with a 128 × 128 pixel detector. The ground truth intensity pattern is taken to be the squared modulus of the FT of this array. We nominally associate a beam energy of 9 keV to this ground truth. A data set of this kind could conceivably be collected at an existing BCDI facility.
Simulation of an overbinned diffraction signal at higher energies is now straightforward: for each 2D slice in the ground truth, blocks of pixels are summed to a single intensity value (Fig. 2). This binning operation mimics a high-energy BCDI experiment since photons that would have spread over a larger solid angle at lower x-ray energies now aggregate into fewer pixels at higher energies because of the compression of reciprocal space. Equivalently, one may imagine the overbinning to arise from a 9 keV measurement with proportionately larger pixels. The ratio of pixel sizes ('pixel binning factor' of PBF) is equal to the ratio of the beam energies. For example, the diffraction features at 9 keV contained in every 6 × 6 block of pixels is squeezed into a single pixel when the beam energy is 54 keV. In this case, PBF = 6.
In this manner binning effectively reduces the feature visibility of a BCDI intensity measurement ( Fig. 1(b)), which we redress through additional information collected from sub-pixel detector translations. A similar method for ptychography has been proposed by Batey et al. 13 . In general the finite size of the physical pixels makes it impossible to sufficiently constrain the intensity features of the ground truth that are lost to binning, regardless of how finely the detector is translated in its plane. This assertion is rigorously proved in Section 1 of the Supplementary Material. As mentioned earlier, our approach instead is to demonstrate sub-pixel-scale feature recovery with fewer measurements, which is enabled by the sparsity of the Fourier representation of the scattered intensity. This is the focus of Section 3 in which we discuss our binned data sets in the context of a compressed sensing measurement and motivate the use of sparse recovery techniques.
We point out that in a real-world BCDI experiment at energies >9 keV, a given detector would span a larger q-space aperture that would not be entirely covered at 9 keV. In our simulations we consider a fixed q-space aperture corresponding to the 128 × 128 detector space of the ground truth. We focus on the recovery of fine features in the original ground truth simulation from binned data sets that subtend this q range. We justify this by pointing out that in typical BCDI experiments at 9 keV the edge pixels of the detector capture relatively low intensities compared to the central pixels (by a few orders of magnitude). The scattered intensities outside this aperture contribute negligibly to the Fourier representation of the original scatterer.

Sparse recovery: Mathematical Details
We seek a method to reverse the binning process described in Section 2 and obtain the original ground-truth diffraction pattern that features well-resolved fringes. In our method, recovery of the fine detail from a limited set of binned pixel measurements hinges on representing each two-dimensional slice in the discrete cosine basis, which we utilise as a numerically convenient variant of the FT. The two-dimensional discrete cosine transform (DCT) is defined for an N × N-sized image A ij as the linear transformation 14 : The three-dimensional diffraction pattern from a finite crystal is inherently a band-limited function, since its FT is the auto-correlation of the scattering factor of a finite crystal (the Patterson function, Fig. 3(b)). From this it follows that each two-dimensional slice of the diffraction pattern is also band-limited (this is proved rigorously in Section 2 of the Supplementary Material). Since each slice contains essentially the same signal information as its two-dimensional DCT, one would require approximately as many pixel measurements as the number of non-zero components in this latter representation in order to fully determine it. This well-known principle from information theory underlies all compressed sensing/sparse recovery algorithms. This approach is of immense importance in scenarios where a scant number of physical measurements seemingly fail to mathematically constrain a system. After sparse recovery, the target intensity pattern we ultimately seek is obtained simply by inverting this transform.
To avoid information loss in our ground truth simulation, the Patterson function of the simulated crystal has to be fully captured within the three dimensional simulation array. In other words, the buffer width (populated with zero-valued pixels) around the scatterer should be at least the span of the scatterer in each dimension (this is equivalent to the Nyquist sampling criterion 17 ). The simulation sizes described in Section 2 go well above this minimum requirement to ensure that the Patterson function is not only fully captured, but also sparse in the array.
We define the sparse recovery problem in terms of (1) The unknown (target) diffraction slice I of the desired fine-pixel resolution, which is sparse in the DCT representation: I ≡ Bx (where x is sparse and the columns of matrix B are the inverse DCT basis vectors) and (2) a set of measurements made on this signal, represented by a linear operation A resulting in measured (binned) values I (measured) = AI = A(Bx). We illustrate these operations using the example intensity pattern shown in Fig. 2 that has a relatively small number of pixels. In this example, the coarse pixelation of a fine grid of size 28 × 28 = 784 pixels results in a total of 4 × 4 = 16 measurements, and therefore A is a 16 × 784 matrix. Additional measurements are obtained by translating the detector in the plane perpendicular to the exit beam as shown. Each detector offset provides a new set of pixel measurements (fewer than 16 per detector translation, since we ignore the pixels that fall outside the original q range of interest). Each pixel measurement corresponds to a row of the matrix A. As we have proved in Section 1 of the Supplementary Material, it is not possible to obtain 784 independent rows of the matrix A through detector shifts, and we address this problem with compressed sensing. All compressed sensing techniques solve the system of equations ABx = I (measured) by enforcing that x be sparse.
While there exist various algorithms for sparse recovery 18-20 , we adopt the LASSO regression method 21 common in machine learning applications: for some small α > 0 (set to 2 × 10 −4 throughout this paper). The |x| penalty imposition on the objective function in Equation (2) explicitly enforces sparsity on the unknown x ( 1  -optimisation). Briefly, the optimisation converges to that solution x which has the fewest non-zero components and simultaneously satisfies the set of constraints ABx = I (measured) . The target image is recovered by inverting the sparsifying transform (I = Bx optimal ). We point out that this recovery scheme does not strictly enforce the constraint of non-negativity on I. This results in recovered images with a few negative-valued pixels, which we simply threshold to zero ( Fig. 4(g,h)). We show in Section 4 that this thresholding still results in a working approximation of the continuous signal, from a phase retrieval point of view.
With regard to the number of independent coarse pixel measurements M, compressed sensing theory provides a more precise success threshold: for an image of size N × N pixels, M ≥ Klog(N 2 /K) where K is the number of significant (nonzero) components in the sparse representation 11 . In general, results of the recovery when M is well above this threshold are negligibly different from the original image itself. We note that K depends on the basis of the sparse representation (in our case, the discrete cosine basis). Given the simulated particle described in Section 2, a crude estimate of K for the DCT can be made in the following manner. Since the central slice of the particle is contained within a 22 × 22-pixel box at the centre of a two-dimensional array, any projection of the Patterson function parallel to this plane is contained in a 44 × 44-pixel box (from the projection-slice theorem 22 ). Since our simulated object has facets, this projection may be simply modeled as a 'diamond'-shaped area whose vertices lie on the midpoints of the sides of the 44 × 44-pixel square. The 'diamond' has an area of 44 2 /2 = 968 pixels. Therefore the estimated number of non-zero components in this projection (which is the Fourier transform of a single detector measurement) is approximately 2 × 968 = 1936, after accounting for comlpex numbers. Given the close relation between the Fourier and cosine transforms, it is reasonable to assume a similar number of components in the cosine basis as well. In other words, K ≃ 1936.
Our estimate of ≃1500 comes from simply counting the number of significant components in Fig. 3(c) after suitable thresholding. In a real experiment it is impossible to estimate K without making a measurement. Figure 4 shows a visual comparison between a benchmark ground truth diffraction slice (4(b)), and the corresponding fine-pixel diffraction slice recovered from a set of simulated intensity patterns at 54 keV obtained with sub-pixel detector translation (4(e)). The coarsely binned diffraction slice has size 20 × 20 pixels (4(d)), while the recovered diffraction slice has size 120 × 120 (fine) pixels, corresponding to a PBF of 6.

Recovery Results
For a better quantitative picture, we examine the fidelity of two different recovered diffraction patterns to their respective 9 keV ground truth benchmark slices, as a function of the beam energy and degree of upsampling. The benchmark images represent qualitatively different intensity distributions: the "on-Bragg" central slice (Fig. 5(a)) has one strong, highly localised peak while the intensity distribution in the "off-Bragg" terminal slice (Fig. 5(b)) is weaker and more spread out. The fidelity is quantified by the sparse recovery transfer function (SRTF), which we define for a single recovered image as recovered groundtruth where the indices (i, j) run over the pixels of the finely pixelated diffraction pattern. The SRTF is analogous to the phase retrieval transfer function common in phase retrieval literature 23 . For a perfect recovery, SRTF = 1 for all pixels. Our results are expressed in terms of the mean μ of the SRTF and the standard deviation spread μ ± σ around the mean, evaluated over the recovered (upsampled) images.  With higher beam energy, there is a greater variance in the SRTF. This is due to the appearance of artefacts from an insufficiently constrained optimisation, as can most prominently be seen in Fig. 6(a) and (b). In Fig. 5(c) and (d), on the other hand, the slight deviation of the SRTF from 1 at high PBF stems from the small penalty α|x| on the objective function in Equation (2). Finally in Fig. 5(f-h) we show the results of the standard BCDI phase retrieval applied to diffraction intensity patterns upsampled from simulations of 36, 45 and 54 keV data sets (PBF = 4, 5, 6) with sub-pixel translations. Each upsampled intensity pattern data set was recovered from the original detector position and 12 additional detector offsets along the pixel diagonals. The reconstructions are in agreement with the original structure 5(e). Section 1 of the Supplementary Material describes the rationale for offsets along both pixel diagonals, as opposed to purely horizontal or vertical offsets.

Discussion
We have outlined a framework by which fine features in high-energy coherent diffraction intensity patterns from individual compact crystals can be recovered using compressed sensing and sparse recovery. The success of these methods is often justified in literature by demonstrating recovery with fewer measurements than the size of the actual signal array. Thus far we have also described our methodology in terms of these coarse 'pixel measurements' when quantifying information content in a signal. In this section we translate this phrasing into quantities that are directly related to the design of CDI experiments. Specifically, we show how our sparse recovery technique can be exploited to emulate smaller pixel sizes or larger sample-detector distances. We thereby demonstrate the advantage of putting the burden of signal detection on numerical algorithms, as opposed to substantially modifying an existing experimental setup or constructing massive experimental enclosures to enable large detector distances at large Bragg angles.
Firstly, Table 1a shows the number of pixel measurement constraints in our simulations as a function of the X-ray energy multiplier and the number of detector positions. The maximum number of pixel measurements at each PBF is computed using Equation 2 of the Supplementary Material. The two regions of interest are: (1) the numbers in ()-brackets which fall below the theoretical limit of M = Klog(N 2 /K), and (2) the numbers in []-brackets which show that further detector offsets along the diagonals do not give additional binning constraints due to the grid periodicity. Here, we aimed to recover the ground truth within a 120 × 120 pixel region, such that N = 120. The number of unknown sub-pixels to recover is therefore N 2 = 14400. Also, K ≃ 1500 for the simulated crystal (i.e. the number of significant components in the DCT of the central diffraction slice, Fig. 4(f)). The numbers in region (1) result from an inadequate number of detector positions. In this regime, the acquired signal information is truly deficient and sparse recovery is not possible. Region (2) denotes the limitations of the chosen detector translation strategy (in this case, diagonal offsets only). For a given crystal size and detector pixel size, appropriate choices of sample-detector distance and detector offset strategy ensure the existence of the intermediate region between 1 and 2 where compressed sensing can successfully recover the diffracted intensity pattern.
The connection between M and the experimental parameters is made if we imagine that M binning constraints can also be obtained from a single detector image of × M M pixels that queries the same region of interest in reciprocal space. Stated differently, the set of offset detector images with a larger pixel size is conceptually equivalent to a single detector image with a smaller pixel size. From an information theoretical point of view, the two sets of constraints are equivalent descriptions of the scattered wave intensity. This immediately suggests a reinterpretation of Table 1a in terms of experimental parameters.
For example, at a beam energy of 54 keV, we have PBF = 6 relative to the same measurement at 9 keV. The original 120 × 120-pixel q-space region of interest in the ground truth simulation is now squeezed into a 20 × 20-pixel region in the centre of the detector due to reciprocal space compression. The effective coarse pixel size over the original region of interest is simply Δx = 1/20 = 0.05 units. But with 7 detector positions at this energy, Table 1a shows that we have access to a total of 2566 pixel measurements, effectively emulating a × 2566 2566 pixel grid over the same region of interest. The effective pixel size is now scaled by a factor of , or equivalently, we are able to resolve spatial frequencies better by a factor of  = . 2566 20 2 533, which we label f.
In the Fraunhofer approximation, the quantity f can also be interpreted as an effective increase in sample-detector distance, given a fixed pixel size. The spatial frequency step is related to the pixel size Δx and sample-detector distance z by Δk ∝ Δx/z. A scaled pixel size of Δx/f implies an effective sample-detector distance of fz. Table 1b shows Table 1a reinterpreted in terms of the distance multiplier f.
With a fixed beam energy, the limit of f is reached when all possible binning constraints are accounted for and further detector offsets do not give new constraints. We have shown in Section 1 of the Supplementary Material that in this case, the number of pixel measurements that fall fully within the q-space region of interest Figure 7 shows the maximum values of f under this condition, for a wider range of beam energies. We see that maximum upsampling nearly undoes the effect of pixel binning at high energies, through a proportionate increase of f, effectively emulating the sample-detector distance that would otherwise be needed to resolve intensity fringes. The difference is accounted for by the fact that pixels that are not fully contained within the original q-range are not used as measurement constraints.

Summary
We have described a signal recovery technique for single-crystal BCDI data sets acquired at X-ray energies typically suited to HEDM applications (>50 keV). Our methodology relies on a modification of the conventional BCDI setup in which additional data is acquired by translating the detector perpendicular to the diffracted beam. We have described a fundamental incompleteness in the data collected by this method of upsampling, which calls for the incorporation of some form of extra information about the signal. With this necessity we have motivated the use of techniques that recover not the signal directly, but sparse representations of it. We have based our methodology on a rigorous proof of the existence of such a representation for all coherent scattering from a fully-illuminated compact single crystal. We have quantified image recovery by this compressed sensing process We note that the upsampling of coarsely pixelated data in some transformed form is in general a valuable alternative to direct measurement for many high-energy diffraction experiments. The subsequent signal recovery is readily achieved by identifying an appropriate sparse basis for a given measurement.
Information theory-based analysis of detector-space upsampling shows us that smaller pixels and larger sample-detector distances can be emulated, with minimal changes to the experimental setup. This capability has the potential to significantly influence the design of space-constrained experiments at high energy coherent scattering beamlines. When incorporated into existing near-and far-field HEDM workflows, high-energy BCDI could be the nanoscale component of a generalised experiment for imaging and characterisation of polycrystalline samples at multiple length scales. Coupled with the development of new computational methodologies, such multiscale characterisation capabilities could go a long way in validating existing and new models of materials physics, as well as informing the creation and processing of engineering materials. Data availability. The ground truth simulation and the Python code for sparse recovery are available upon reasonable request. Please contact the corresponding author. For a fixed pixel size, the distance multiplier asymptotically approaches the PBF in the limit of large reciprocal space aperture (a wider detector).