Speckle-based X-ray microtomography via preconditioned Wirtinger flow

Sample preparation

Glass beads (≤106 µm, 2.5 g cm−3, G8893, Merck KGaA) are attached to a toothpick with an adhesive for the sample shown in Figs. 2 and 3. The samples shown in Fig. 4 (a cumin seed, a dried anchovy, a dried shrimp, and a piece of cork) are procured from a local store. The piece of cork is obtained from a wine cork.

Experimental setup

The experiments are performed at the 6 C beamline of PLS-II in Korea. A wiggler source is used (500 µm × 30 µm full width at half maximum), followed by a double multilayer monochromator (DMM) tuned to 10 keV (∆E/E = 2.1% measured at 16 keV). The sample is placed 36 m after the source. The diffuser is four layers of P3000 grit sandpaper (991 A, Starcke GmbH & Co. KG) placed L1 = 3 mm after the sample. The scintillator is 50-µm thick LuAG:Ce placed L2 = 20 mm downstream of the diffuser. Note that the shortest possible L1 and L2 are chosen to minimize sample-induced decoherence (Supplementary Text). Specifically, L1 is chosen to ensure that it does not physically interfere with the sample rotation stage, and L2 is determined to ensure that the spatial bandwidth of the speckle pattern is well sampled without aliasing. Note that the speckle grains become coarser as L2 increases due to the partial coherence (Fig. S2). The image on the scintillator is observed with an optical microscope equipped with an objective lens (NA = 0.25, LMPLFLN 10X, Olympus Corp.) and a sCMOS camera (6.5 µm, 2560 × 2160, pco.edge 5.5, Excelitas PCO GmbH). The corresponding spatial sampling period and area are 650 nm and 1.664 mm × 1.404 mm, respectively.

Data acquisition

For each sample, data are acquired in the following order: (i) dark frame, (ii) flat field for the beam image, (iii) reference speckle from the diffuser, (iv) sample projection images. In step (iv), the samples are rotated continuously from 0° to 180°, while the projection images are taken on the fly at each 0.225° rotation angle, resulting in a total of 801 images. The camera exposure time is kept at 100 ms throughout the measurements. For speckle-tracking methods, which typically require 10 to 20 independent measurements27,29, the acquisition process is repeated 12 times at different lateral positions of the diffuser. In all subsequent reconstruction steps, all speckle images are considered to have been pre-processed with dark frame subtraction and flat-field correction. The entire acquisition process, including all stage movements and pauses, is completed in approximately 4 min (see Video S6).

Intensity optical transfer function (IOTF)

The coherence length at the detection plane is calibrated from the autocorrelation function of the reference speckle pattern based on the Siegert relation

while \({g}^{(1)}(\Delta x)\) and \({g}^{\left(2\right)}\left(\Delta x\right)\) are first- and second-order normalized autocorrelation functions of a speckle pattern61. Note the \({\left|{g}^{(1)}(\Delta x)\right|}^{2}\) term is closely related to the spatial coherence of light62. Due to the inherent asymmetry of the synchrotron source, horizontal and vertical coherence lengths (xcoh and ycoh) are estimated separately. For example, for xcoh, we compute the one-dimensional normalized autocorrelation function along the x-axis, average it along the y-axis, and estimate the coherence length by fitting it to the Gaussian model \({g}^{(2)}(\Delta x)=1+\exp \left(-\pi \Delta {x}^{2}/{x}_{{\rm{coh}}}^{2}\right)\)62. Corresponding intensity point spread function (IPSF) becomes

where r is the index of vectorized image in (x, y) space. By the Fourier transform of Eq. (4), we can calculate the intensity optical transfer function (IOTF),

where u and v are the spatial frequencies of x and y, respectively, and k is the index of vectorized image in (u, v).

Diffuser transmission function

We derived the transmission function of the diffuser (tr) from the reference speckle. First, we removed the spatial incoherence effect using Wiener deconvolution with the calibrated IOTF (Eq. 5). This step estimates the ‘coherent’ reference speckle intensity at the camera plane. Then, the tr can be obtained by the numerical backpropagation (L2) of the deconvoluted reference speckle. One problem here is that its phase is unknown.

We originally attempted to retrieve the sample and diffuser fields simultaneously, as done in ptychographic iterations47,48. However, we found that the obtained diffuser phase is usually unreliable (e.g., slowly varying, spatially inhomogeneous) and changes drastically depending on the initial guess (Fig. S12). Despite this instability, the sample field was retrieved in a stable manner with indiscernible tomogram results. It is noteworthy that this result agrees with near-field ptychography results32,63,64. This may suggest that there is only one sample field that satisfies both the sample and reference speckle simultaneously, regardless of the phase of the reference speckle. In other words, the sample field can be retrieved with any of the possible diffuser phase. Based on this, we simply fix the phase of the reference speckle to zero.

This is a significant advantage over conventional speckle-based imaging methods, which usually require a well-defined complex diffuser transmission function34,50. We expect this advantage to hold only for weakly scattering samples with a near-field setup that does not induce significant overlap between speckle grains. Further investigation is needed to determine these conditions.

Preconditioning filter

The phase retrieval problem in speckle-based X-ray microtomography is well conditioned on the phase gradient rather than the phase. As such, the gradient with respect to the phase gradient would have a better convergence property. Thus, at each iteration, we want to perform the gradient descent step to the phase gradient such as

where L is the loss function defined in Eq. (2), \({\phi }_{r}={\rm{Im}}({\psi }_{r})\) is the sample phase, and η is a step size. Applying additional partial derivative to both sides of Eq. (6) and adding them together, we get

where \({\nabla }_{\perp }^{2}={\partial }^{2}/\partial {x}^{2}+{\partial }^{2}/\partial {y}^{2}\) is a two-dimensional Laplacian. The second term on the right-hand side can be simplified by the chain rule

which leads

where the factor of 2 is omitted here as it is coupled with the step size. Applying inverse Laplacian to both sides of Eq. (9), we obtain

Compared to the conventional gradient-based phase iteration, \({\phi }_{r}\leftarrow {\phi }_{r}-\eta \partial {\mathcal{L}}/\partial {\phi }_{r}\) there is an additional inverse Laplacian term before the gradient term, which is the preconditioner. The preconditioner is introduced by applying the inverse quadratic preconditioning filer \({P}_{k}^{-1}:=1/({u}^{2}+{v}^{2})\) to the update term,

where u and v are the spatial frequencies of x and y, and \(\mathcal{F}\left\{ \cdot \right\}\) is the Fourier transform. For the zero frequency (u = 0 and v = 0) \({P}_{k}^{-1}\) is set to zero, and is determined separately later by zeroing the background phase value. Since \({\phi }_{r}={\rm{Im}}({\psi }_{r})\), the preconditioning filter is applied only to the imaginary part of \({\psi }_{r}\). The detailed calculation can be found in Algorithm S2.

Regularization window

In this work, the number of measured spatial modes is defined by the IOTFk, which effectively limits the measurable spatial bandwidth. Since Eq. (5) gives the effective bandwidths of \(\sqrt{2}{x}_{{\rm{coh}}}^{-1}\) and \(\sqrt{2}{y}_{{\rm{coh}}}^{-1}\) 62, we can estimate the number of measured spatial modes from the space-bandwidth product (SBP) of the measured images50,

where FOVx and FOVy are the horizontal and vertical extents of acquired images (i.e., field of view). Note that the factor 1/2 is introduced here to correctly estimate the SBP of speckle field from the SBP of intensity speckle image. Assuming the Gaussian angular spectrum, the intensity speckle would have greater bandwidth by the factor of \(\sqrt{2}\), which results in the SBP increase by the factor of 2. The number of reconstructed spatial modes (N) can be defined from given oversampling ratio γ and Eq. (12),

Since the reconstructed fields have the same spatial extent as the measured images, Eq. (13) directly leads to additional constraints on their Fourier space by the bandwidth of \({x}_{{\rm{coh}}}^{-1}/\sqrt{\gamma }\) and \({y}_{{\rm{coh}}}^{-1}/\sqrt{\gamma }\). Assuming a Gaussian profile, we can define the sample transfer function (STFk),

and subsequently the regularization window,

The used oversampling ratio (γ = 1) may seem inconsistent with previous works that empirically report γ ≥ 4 for robust field reconstruction33,43. The difference is that here we use a Gaussian STF, which provides a soft boundary, while previous works used a well-defined domain, which is equivalent to binary STFs. We find that the soft boundary relaxes the dependence of the reconstruction fidelity on γ (Fig. S4), and we prefer to use a smaller γ to obtain the best possible resolution.

Due to the Gaussian shape of the STFk, it does not have a distinct limit on the resolution; rather, it depends heavily on the fineness of the sample structure and the noise level. For instance, if the sample bandwidth is retrievable up to STFk = 0.1 boundary, the expected spatial resolution can be calculated \(\sqrt{\pi }{x}_{{\rm{coh}}}/\left(2\log 10\right)\approx 0.385{x}_{{\rm{coh}}}\) from Eq. (14). Using the measured coherence lengths in Fig. S2, the expected horizontal and vertical resolution becomes 1.34 µm and 1.66 µm, respectively. The root-mean-square resolution is 1.51 µm, which agrees with the acquired resolution in FSC analysis (Fig. 3f).

Reconstruction procedure

The PWF iteration is stopped when the normalized correlation between the retrieved field of the current iteration and the previous iteration is greater than 10−0.00001. Typically, 700 to 1300 iterations are required to meet the criterion, which takes 5 s to 6 s on a personal computer equipped with a graphics processing unit (GPU; GeForce RTX 4090, NVIDIA Corp.) using MATLAB software (The MathWorks, Inc.). Thus, for each sample, it took 70 min to 80 min to process all 801 projection angles. Note that all the computations in Algorithm S1 are element-wise operations and Fourier transforms, which were effectively accelerated by the GPU.

In both PWF and LCS, the root-mean-square error (RMSE) image is calculated from the reconstructed sample fields. By applying reconstructed sample fields to the forward model, we calculate the estimated speckle patterns, \(f({x}_{r})\) and subtract it from the measured speckle patterns, \({y}_{r}\). In PWF, there is only one speckle pattern to compare (K = 1), so the RMSE image simply is \(\left|{y}_{r}-f({x}_{r})\right|\). In LCS, we have twelve speckle patterns to compare (K = 12), so the RMSE image becomes \(\sqrt{\sum_{m=1}^{K}{\left({y}_{r}^{m}-{f}_{{\rm{LCS}}}\left({x}_{r}^{m}\right)\right)}^{2}}\), where fLCS is the forward model of LCS. The RMSE images are normalized with the mean intensity of the reference speckle and present in % (Fig. 2d, h).

The tomogram results are reconstructed using the standard filtered back projection (FBP) algorithm with Ramachandran-Lakshminarayanan (“Ram–Lak”) filter. The iradon function of MATLAB is independently applied at each vertical position of the sample. A constant \(\lambda /p/(2\pi )\) is then multiplied to convert the tomogram results to \(\delta ({\bf{r}})\) and \(\beta ({\bf{r}})\), where p is a pixel size of the reconstructed sample tomogram. We often experience low-frequency background curvature due to the intrinsic lower sensitivity on slow phase variation of phase-gradient sensing techniques54. We mitigate such artifacts by fitting the background area with a two-dimensional quadratic function after tomographic reconstruction.

Fourier shell correlation (FSC)

For FSC, we utilize two tomograms of the toothpick sample (Fig. 3) that are reconstructed from the measurements at different lateral positions of the diffuser. Defining either the real or imaginary parts of two tomograms as \({g}_{1}({\bf{r}})\) and \({g}_{2}({\bf{r}})\), the FSC can be calculated from the normalized cross-correlation of their Fourier transforms over the shell,

where \({\tilde{g}}_{1,2}({\bf{k}})\) is the Fourier transform of \({g}_{1,2}({\bf{r}})\), and κ is the radius of the shell in the 3D Fourier space65. For PWF, a total of 12 tomograms are independently reconstructed from the K = 12 measurements, and two are randomly selected for FSC calculation. Five FSC are calculated from different pairs of tomograms and averaged. For LCS, the K = 12 measurements are randomly divided into two sets of K = 6, and the tomograms reconstructed from the two sets are used for the FSC. Similarly, five FSC are calculated from different random splits and averaged.

Regarding the resolution criterion, we followed the derivation in Ref. 65, which suggests the 0.143 (or 1/7) criterion for the “full-set averaged” result. For example, if there are K repetitive measurements, the FSC is calculated between the two half-set averaged results of K/2 non-overlapping measurements. To extrapolate this FSC result to the resolution of the full-set averaged result, an improvement in signal-to-noise ratio (SNR) by \(\sqrt{2}\) was assumed. Unfortunately, this result cannot be applied directly to our case for two reasons: first, PWF requires only a single measurement (K = 1), which cannot be divided into two sets; and more importantly, the SNR may not be simply proportional to K because noise propagation in PWF is more complex due to the nonlinear reconstruction process. Therefore, we redo the derivation without the full set extrapolation.

Let us divide the measured tomogram into the sample signal and noise terms, \({g}_{1,2}({\bf{r}})={s}_{r}({\bf{r}})+{n}_{1,2}({\bf{r}})\). Note that \({s}_{r}({\bf{r}})\) is invariant across measurements, while \({n}_{1,2}({\bf{r}})\) fluctuates randomly. Here we can consider the result to be reliable if the normalized correlation between the data and the ideal sample signal is greater than 1/2,

where the \({{||}{\bf{k}}{||}}_{2}=\kappa\) condition of the average operator, and the argument of \(\tilde{g}({\bf{k}})\), \(\tilde{s}({\bf{k}})\), and \(\tilde{n}({\bf{k}})\) are dropped for brevity. Note that the \(\left\langle {\tilde{n}}^{* }\tilde{s}\right\rangle\) and \(\left\langle {\tilde{s}}^{* }\tilde{n}\right\rangle\) terms have disappeared since the noise should be random and uniform over space, which leads \(\left\langle {\tilde{s}}^{* }\tilde{n}\right\rangle =\left\langle {\tilde{s}}^{* }\right\rangle \left\langle \tilde{n}\right\rangle\) and \(\left\langle \tilde{n}\right\rangle =0\) for the non-zero spatial frequencies (\(\kappa > 0\)), respectively. Similarly, we can rewrite Eq. (16) as

Substituting Eq. (18) into Eq. (17), we get \({\rm{FSC}}(\kappa)\geqq 1/4\), which suggests 1/4 as the resolution criterion. Again, the only difference from the derivation in Ref. 65 is the full set extrapolation that introduces an additional 1/2 factor for \(\left\langle {\tilde{n}}^{* }\tilde{n}\right\rangle\) in Eq. (17).