Phase tomography from x-ray coherent diffractive imaging projections.

Coherent diffractive imaging provides accurate phase projections that can be tomographically combined to yield detailed quantitative 3D reconstructions with a resolution that is not limited by imaging optics. We present robust algorithms for post-processing and alignment of these tomographic phase projections. A simple method to remove undesired constant and linear phase terms on the reconstructions is given. Also, we provide an algorithm for automatic alignment of projections that has good performance even for samples with no fiducial markers. Currently applied to phase projections, this alignment algorithm has proven to be robust and should also be useful for lens-based tomography techniques that pursue nanoscale 3D imaging. Lastly, we provide a method for tomographic reconstruction that works on phase projections that are known modulo 2π, such that the phase unwrapping step is avoided. We demonstrate the performance of these algorithms by 3D imaging of bacteria population in legume root-nodule cells.


Introduction
Three-dimensional structure determination at nanoscale resolution is highly desirable in both biology and materials science applications. X-ray nanotomography [1][2][3][4][5] can provide this in a non-destructive way, allowing to preserve the sample for subsequent tests and reducing the possibility of introducing artifacts as may occur using slice-and-image methods. The high pen- etration depth of hard x-rays can be used to image sample volumes of several tens of microns of side, which are sought to obtain images that are statistically representative of bulk materials and to allow imaging of large biological specimens in their native environment, e.g., tissues and cellular networks. The smaller wavelength of hard x-rays, compared to soft x-rays, also increases the maximum sample thickness at which the projection approximation is valid for a fixed resolution.
Despite the significantly reduced absorption of hard x-ray energies, the information about the sample is available as a phase shift of the incident wave. This phase shift provides an imaging contrast mechanism that is not inherently linked to absorption or radiation damage, and makes phase tomography an attractive technique in this wavelength regime. Phase-contrast tomography with x-ray Zernike lenses [3] can be used for this purpose but is restricted to weakphase objects to remain quantitative and free of artifacts, and grating interferometer methods provide quantitative phase tomography [6] but remain difficult to scale down below a micron resolution.
A feasible route to retrieve the phase quantitatively is coherent diffractive imaging (CDI). CDI requires no imaging lenses or optics [7,8] but rather uses intensity measurements of x-ray diffraction patterns and iterative algorithms to reconstruct the transmissivity of the sample. The resolution of CDI is limited by the highest angle at which data is measured with a moderately good signal-to-noise ratio (SNR).
The transmissivity reconstructions are complex-valued, and the phase of the reconstruction corresponds to the integrated, or projected, phase advancement of the incident wavefront along the direction of propagation. This can be retrieved with high sensitivity [9] and is not limited to a weak-phase approximation. If the phase projections are measured at different sample orientations, they can be tomographically combined, as schematically shown in Fig. 1, to yield a three-dimensional map of the sample refractive index [5,10]. This allows the quantitative measurement of minute differences in electron density [5,11] and offers a path for quantitative tomography at the nanoscale.
Prior to the tomographic step a few processing steps need to be carried out. Depending on the CDI method used, the reconstructions may suffer from additive constant and linear phase terms that should be corrected. Also, the phase of the projections is computed from complex-valued reconstructions and is thus only known modulo 2π. This wrapped phase cannot be directly used in a filtered back projection (FBP) algorithm. In previous work this has been dealt with by a phase unwrapping step.
Additionally, the phase projections may be misaligned due to hardware limitations. This can be a problem for any nanoscale tomography technique where the axis of rotation and sample position needs to remain stable throughout the duration of the angular scan. Keeping the sample position to within tens of nanometers while rotating can be technologically challenging. This is especially true when full measurement times exceed a few hours, which is sometimes necessary for CDI to accumulate x-ray signal at high diffraction angles [4,5].
In this paper we describe effective and computationally efficient methods to deal with the above mentioned artifacts and limitations of phase projections from CDI. We provide a robust and accurate algorithm for alignment of tomographic phase projections when no fiducial markers are available. Additionally, we provide an algorithm for tomographic reconstruction that is insensitive to phase wrapping, thus avoiding the often problematic phase-unwrapping step.
Methods do exist that can avoid many or all of these projection post-processing steps by directly assembling the diffraction measurements in 3D and performing iterative reconstructions in 3D arrays [12,13]. However, they require isolated, small objects which makes sample preparation more difficult and limits their applicability to larger samples.

Experiment description
We demonstrate the performance of our tomographic methods by reconstructing the refractive index distribution of cells from the root-nodules of the mung bean plant Vigna radiata. The nodule cells host the rhizobia bacteria Bradyrhizobium japonicum, with which the plant sustains a symbiotic relationship. The legume plant provides the bacteria with a necessary anaerobic environment and nutrients, and receives in exchange nitrogen compounds that are synthesized by the rhizobia. This process provides the plant with natural fertilizer, and its study is motivated by the positive impact on the legume plant, which is an important source of proteins.
We have chosen this sample, originally measured for a correlative study, as it presented elements that made the processing of phase projections particularly challenging. Firstly, the small contrast of features in the projections made the alignment difficult with other alignment methods such as those based on fiducial markers or cross-correlation of projections. Additionally, to avoid radiation damage the sample had to be maintained at cryogenic temperature which induced thermal drifts.
To permit a subsequent correlative study, the plant tissue was prepared following a specially adapted protocol suitable for light and electron microscopy (EM) [14]. After high-pressure freezing, a freeze-substitution and subsequent resin-embedding (Lowicryl HM20) was performed in acetone, including fluorescent stains as well as uranyl acetate and osmium tetroxide. The latter two are fixatives that preserve the cellular ultra-structure and simultaneously work as heavy metal stains to provide contrast for EM. A 35 μm diameter cylindrical section was milled from the block of embedded tissue using a focused ion beam (FIB) and attached to a holder suitable for the tomography experiment. Sample preparation was carried out at Electron Microscopy ETH Zurich (EMEZ), Switzerland.
The measurements were performed at the cSAXS beamline at the Swiss Light Source, using the experimental setup for ptychographic tomography [5], at a photon energy of 6.2 keV. A coherent region of the incoming x-ray beam was selected using a 3 μm diameter pinhole, located 5 mm upstream of the sample. The complex-valued projections were obtained using ptychographic CDI [15], also referred to as scanning x-ray diffraction microscopy [16]. For each projection angle, the sample was translated to different positions in a direction transverse to the direction of x-ray propagation, covering a field of view of 50μm×30μm, and x-ray diffraction patterns were measured for each position at 7.2 meters from the sample using a PILATUS photon-counting detector with 172 μm pixels [17].
By ensuring a significant overlap of the x-ray beam for neighboring sample positions we provide diversity of measurements to the phase retrieval algorithm [18][19][20], which removes ambiguities and improves convergence and accuracy in the reconstruction. A total of 530 diffrac-tion patterns were measured for each projection angle, and we repeated this procedure at 360 sample orientations spaced by 0.5 degrees with a total measurement time of 14 hours. The reconstruction of the complex-valued projections and incident illumination was carried out using 128 × 128 pixels of the diffraction patterns with the difference-map algorithm [16,21,22], resulting in a reconstruction pixel size of 65 nm.

Processing and alignment of projections
For our experimental parameters the transmissivity of the sample is well approximated by a complex-valued projection, where (x, y) are the transverse Cartesian coordinates and z the direction of x-ray propagation, λ is the wavelength of the illumination and n(r) complex-valued index of refraction. The sample is rotated by an angle θ with the axis of rotation parallel to y, i.e. r = (x cos θ − z sin θ , y, z cos θ + x sin θ ).
The phase retrieval reconstructed projections, o(x, y, θ ), can additionally have a constant and linear phase terms, namely where a θ , b θ , and c θ are real-valued coefficients. The constant phase term, a θ , arises from an inherent ambiguity in phase retrieval reconstructions. Because we measure intensity patterns, any information about an overall phase shift provided by the sample is lost. Or in other words, the absolute phase of the sample transmissivity has no effect on the measurement and is thus undetermined upon reconstruction. The constant phase term a θ can then take any value. The linear phase terms, b θ and c θ , can be present in CDI reconstructions if the diffraction patterns are not centered in the computational window. Because a translation in Fourier domain corresponds to a linear phase in object domain, an indetermination of the center of the diffraction patterns leads to a consistent linear phase in the object domain reconstructions. The problem is worse when simultaneously reconstructing the illumination and the sample in ptychographic CDI. For the latter there is an inherent degree of freedom, as shown in Appendix A, that allows the object transmissivity to take arbitrary values of b θ and c θ .
The processing pipeline for the phase projections is depicted in Fig. 2. After the constant and linear phase terms are removed, a well-behaved region on the phase projections is selected and unwrapped for the alignment procedure. Once the translation errors are obtained and the projections aligned, we perform the tomographic step using an algorithm that works on the wrapped phase projections. In this way we can faithfully recover the tomogram even for regions that are difficult to unwrap.

Removal of constant and linear phase terms from the projections
For the removal of the constant and linear phase terms we rely on having measured a region outside the sample, thus knowing that for this region the phase should be zero. For tomography to be quantitative it is necessary to measure the whole sample along the x-direction, which guarantees that we should have such an empty area to work with. In sake of clarity and since we remove these terms individually from each projection, we will momentarily drop the dependance on θ . We define a window function w(x, y), that is unity in an empty region and zero elsewhere, and minimize with respect to a, u and v; where M and N are the dimensions of the array. For small phase errors σ 2 reduces to so that minimizing Eq. (3) is close to minimizing the root mean square (RMS) phase in w(x, y), or alternatively, to removing a least-squares fitted plane from the phase. Note that compared to a traditional least-squares fit, Eq. (3) has the advantage that it is insensitive to phase wrapping discontinuities. Following a derivation similar to that in [23] we analytically minimize σ 2 with respect to a and obtain min a,u,v where The problem is then reduced to finding the peak of this windowed Fourier transform within a small fraction of a pixel. Subpixel precision is necessary because a one-pixel error in the location of the peak of R(u, v) leads to a linear phase with 2π radians peak to valley. In practice we zero pad the windowed function to twice the number of pixels and compute a Fourier transform using the FFT algorithm, such that we locate the peak within half a pixel precision. We then compute the DFT in a small region around this initial estimate with an upsampling of 100. For this last step we use a selective upsampling matrix-multiply DFT [24,25], which allows us to upsample in a small neighborhood of the peak without the need to zeropad and without sacrificing accuracy in the computation. The computation time to perform this subpixel refinement is a small fraction of the time required to obtain the initial estimate [25] and the constant phase offset a is directly given by the phase at the aforementioned peak, After removal of these terms for all reconstructions, we can compute the phase of the transmissivity projections which are related to the real part of the 3D refractive index, δ (r), by a Radon transform.

Alignment of phase projections
Our attempts of alignment by cross-correlation registration of neighboring projections [26] gave unsatisfactory results for this sample. Although a subpixel cross-correlation algorithm [25] adequately aligned neighboring projections, we observed that a positioning error of several pixels accumulated over the 360 projections. The lack of markers or salient features in the projections did not allow fiducial marker tracking [27], and due to the large drifts induced by the cryogenic cooling the alignment for this sample was particularly difficult.
The alignment algorithm to be described in this section applies to unwrapped phase. The tomographic step, as detailed in Section 4, can work with wrapped phase such that currently phase unwrapping is only necessary for the alignment of projections. This gives us the freedom to choose for alignment "well-behaved" regions that are easiest to unwrap. The difficulties in phase unwrapping arise from the presence of residues [28], which are points in the phase that produce an inconsistency in the unwrapped phase and makes the result dependent on the selected unwrapping path. Selecting an appropriate unwrapping path is then the task of the algorithm in [28] and other unwrapping algorithms, such as the quality-guided map [29]. Reference [28] provides an easy and efficient method to locate the residues. A region with no residues is the easiest to unwrap as any path will lead to the same result, hence giving consistent results with any algorithm used. If there are residues present, the phase unwrapping can be made consistent (independent of unwrapping path) by defining branch cuts between pairs of "plus" and "minus" residues, and making sure to not cross these branch cuts in the unwrapping path [28]. Phase unwrapping becomes difficult when there is a big number of residues and their pairing is not so clear. When selecting the "well-behaved" region for alignment we look for areas in the reconstructions where no residues are present or where there are few residues appearing in close together pairs, such that finding the shortest branch cuts is easily automated.
The alignment of projections in x is performed using a center of mass computation within the selected well-behaved region, indicated with a red line in Fig. 3(a) for our plant tissue sample. For cases where the whole horizontal span of the sample is measured, the center of mass in a region of interest (ROI) of each projection, corresponds to the projection of the center of mass of the sample [30] and provides the means for estimating the position of the projections Δx θ . Here mass is to be understood as the integral of the projected phase. Although the center of mass argument also holds for the y-direction, in practice this does not guide the algorithm to a correct vertical alignment since the center of mass computations for the misaligned projections are strongly biased by effectively computing it in different regions. The latter is a consequence of not measuring the whole sample vertical span and for these cases we sought an alternative more robust alignment technique.
For the alignment of projections in the y-direction, we have implemented an iterative approach to maximize the correlation of vertical variations in the mass of the sample. The vertical mass distribution can be obtained in each projection by means of integration of the phase in the x-direction, namely M θ (y) is a function that, in ideal situations, should be the same for all projection angles, and provides the means for vertical alignment. We compute the function M θ (y) in the selected area from each projection and remove Legendre polynomial terms along the y direction using a leastsquares fit. We have obtained best results, in terms of stability and capture range, by removing a constant and linear polynomial term from M θ (y), thus arriving at the functions ψ θ (y), shown in Fig. 3(b). Removing a quadratic or higher order of polynomials allows to align functions ψ θ (y) that have more pronounced differences, due to imaging artifacts for example, but at the expense of a reduced capture range and increased sensitivity to noise. By removing polynomial terms, the functions ψ θ (y) emphasize small sample fluctuations that can be more robustly aligned, and the functions directly provide a visual indication of the sample misalignment. It is therefore convenient to select regions with larger vertical mass fluctuations whenever possible. Furthermore, by removing at least a constant and linear polynomial terms, the vertical align-ment becomes by definition independent of possible errors in the removal of the constant and linear phase terms, as described in Section 3.1. The estimate of the position for the (n + 1)-th iteration, Δy (n+1) θ , is found by minimizing the squared error with respect to the average ψ θ (y) where the average is taken over the projection angle θ . Note that the average is computed using the positions of the previous iteration. Hence, we can minimize E 2 individually for each projection with a line search. We perform this line search with a single pixel precision at first to minimize the need for interpolation.
In the iterative procedure we alternate between x and y alignment until a convergence criterion on the positions is met. On the first iterations ψ(y, θ ) is a smooth function with a broad capture range. With subsequent iterations ψ(y, θ ) is refined, showing more distinct features that allow improved alignment. Typically 10 to 15 iterations are required for the process to converge.
The iterative subpixel refinement step follows. For each projection we compute the contribution to E 2 at three subpixel displacement points, fit a parabola and set the best alignment as the minimum point. A handful of iterations are typically sufficient for convergence of this last step. Figure 3(c) shows the sample position offsets (Δx, Δy) and Fig. 3(d) shows the functions ψ(y, θ ) after alignment. For cases of drifts that are larger than the height of the selected region the algorithm can stagnate, in the case shown in Fig. 3 an initial vertical drift of −4.5 μm over the first 50 projections was obtained by inspection of Fig. 3(b) and given as an initial estimate.
The alignment of a single 1D function per projection makes this process fast, and since we align with respect to an average over all the available projections we find the process robust even with the lack of strong markers.
To further validate the results of our alignment technique we performed a similar experiment on a test sample, consisting of two spheres [11], which could be measured at room temperature. The sample holder was positioned on top of an optical reference mirror and the relative vertical position of the sample and pinhole was tracked with an optical interferometer. We performed two independent tomographic measurements, each with 180 projections spaced by 1 degree intervals. Figures 4(a) and 4(b) show the comparison of the alignment algorithm and interferometric measurement with an RMS difference of 55 nm and 40 nm, respectively. This error is within the expected accuracy of the interferometric setup and can be partially attributed to surface errors of the interferometer reference mirror, which has a manufacturer specified flatness of λ /10 ≈ 63 nm.

Computed tomography from derivative of wrapped phase
Probably the most widely used algorithm for tomographic reconstruction from projections is the filtered back projection (FBP). For this method a Fourier domain filter is applied on the projections which are then back-projected to arrive to the reconstruction [31]. For phase projections the output of the FBP is proportional to the real part of the refractive index, namely where Δ is the length of the side of the voxel, F x {·} is a one-dimensional Fourier transform from x to u coordinates, F −1 u {·} is the inverse FT from u to x coordinates, and |u| is the Fourier domain filter applied to the projections. However if the phase is wrapped, use of Eq. (11) produces strong artifacts. For beam deflection tomography [32] and x-ray differential phase contrast tomography [33], a modified Fourier filter has been introduced that allows direct tomographic reconstruction from the derivative of the phase, Using Eq. (12) the tomogram can be computed directly from measurements of the phase derivative, removing the need for an intermediate integration step that may introduce noise or artifacts. The derivation of Eq. (12) from Eq. (11) follows directly from the FT theorem For the case of wrapped phase, the phase derivative can be conveniently estimated independently of wrapping by where h is the step for the numerical differentiation, and we introduce to avoid interpolating artifacts that would otherwise occur if the wrapped phase was directly used to estimate the derivative. By using Eq. (13), the derivative is computed with the phase mapped to the unit circle, where the apparent wrapping discontinuities disappear, i.e. the function g(x, y, θ ) is well-behaved and insensitive to phase wrapping. Notice that this approach is still restricted to changes of the phase to be less than π radians from one pixel to the next, which is a requirement for adequate sampling of the phase profile in a complex-valued function. We used Eqs. (12) and (13) to perform tomography, with an isotropic voxel of 65 nm side, directly from the wrapped phase projections of the mung bean tissue sample. We included a Parzen filter to mitigate high-frequency noise [34] as it is a standard procedure in FBP reconstructions. Three-dimensional renderings of the mung-bean tissue sample are shown in Figs. 5(a) and 5(c). Tomographic slices are shown in Figs. 5(b) and 5(d). These slices are taken along the planes indicated in Fig. 5(a) and share color-map and scale.
The region of the sample that was imaged partially contains four cells that hosted the bacteria (B), and the reconstruction reveals the 3D distribution of the bacteria within the cells. A layer of a platinum compound (P) is shown on one side of the sample in white, as it has higher electron density (δ ≈ 25 × 10 −6 ). This layer was deposited using the FIB in order to attach the pillar to the micromanipulator and to the sample holder. Cell walls (W), round vesicles (V) and intercellular spaces can also be clearly identified. As the technique is quantitatively sensitive to changes in electron density, the contrast in this sample is enhanced due to the heavy-atom staining. This staining is an important contrast agent for electron microscopy (EM) and thus part of the sample preparation protocol. Since our reconstruction is proportional to the average electron density per voxel, it provides a quantitative map of the staining distribution among and within cell organelles and bacteria, which could be useful in studies of sample preparation protocols.

Summary and discussion
The reduced absorption at hard x-ray energies, compared to soft x-rays, is desirable to measure thick samples. For this wavelength regime the phase shift imparted on the x-ray field by the sample provides a contrast mechanism that is not inherently linked to radiation damage or absorption. Coherent diffractive imaging is a practical way of measuring this sample-induced phase shift with high sensitivity and with a resolution not limited by imaging optics.
We have presented methods for post-processing of CDI phase projections to remove constant and linear phase terms associated with complex-valued reconstructions, align projections, and perform tomography from wrapped phase. Although we have applied these methods to xray ptychographic tomography, they should prove useful in other tomography configurations, including for example x-ray holographic tomography [35] and curved-wavefront CDI [2,10].
The projection alignment algorithm, in particular, can also be of use for conventional absorption or phase-contrast tomography at optical, soft and hard x-ray wavelengths. For the case of x-ray nanotomography, monitoring and maintaining the sample position with nanometer precision while allowing a rotational degree of freedom is difficult, often requiring post-processing alignment. The vertical alignment procedure presented here has proven to be very robust to image artifacts, and has shown to work even for cases when Eqs. (9) are not completely satisfied, such as when only part of the sample is measured in some projections, or for images obtained with an x-ray Zernike phase-contrast microscope.

Appendix A. Linear phase terms in ptychographic CDI reconstructions
In conventional ptychographic CDI an extended object, o(x, y), is coherently illuminated by a spatially confined illumination function, p(x, y). The field leaving the sample, for an object transverse position (x n , y n ) is given by f n (x, y) = o(x − x n , y − y n )p(x, y). (A.1) The measured diffraction patterns are where F n (u, v) is the Fourier transform of f n (x, y). The object translations are such that the illumination overlaps the same region of the object for different measurements. This added diversity more effectively constrains the phase retrieval problem, removing ambiguities and increasing robustness and success rate [15,16,[18][19][20]. In many practical situations the illumination function p(x, y) is not known a priori and the phase retrieval problem is then to find an object and illumination that reproduces the measured intensity patterns. When reconstructing both the object and illumination an ambiguity on the linear phase term arises. An alternative solution that satisfies the object-domain constraints can be constructed by Notice that the fields f n (x, y) differ from f n (x, y) by constant phase term and thus have the same diffraction pattern intensity. Because both |F | 2 = |F| 2 and the overlap constraint is satisfied, the functions o (x, y) and p (x, y) are solutions for any value of b and c. This is an inherent ambiguity of the problem that does not hinder image quality, and is independent of the reconstruction algorithm. Hence, unless some a priori knowledge of the sample or the illumination is used, the object reconstruction can have an arbitrary linear phase provided that the illumination is reconstructed with a linear phase of equal magnitude but opposite sign.