Fast Eikonal Phase Retrieval for High-Throughput Beamlines

TL;DR

This paper tackles the challenge of robust, high-throughput phase retrieval in propagation-based X-ray phase-contrast imaging by extending the Eikonal Phase Retrieval (EPR) framework. It introduces a fast, second-order near-field model that retains the full $O(L^2)$ content of a WKB expansion (WKB0 plus the leading WKB1 correction) while employing FFT-diagonal preconditioning to accelerate convergence; to remain robust beyond sub-pixel shifts, it adds a non-local, multi-pixel shift solver that transports intensity via a mass-conserving WKB0 mapping and uses a discrete adjoint for back-projection. The method also accommodates polychromatic data through a compact spectral discretisation, enabling energy-dependent transport and inversion without sacrificing GPU/FFT efficiency. Compared with the original EPR, the new approach achieves roughly a 2.8-order-of-magnitude speedup per volume for similar hardware, enabling routine nonlinear phase retrieval in high-throughput HiP-CT workflows, with demonstrated robustness in challenging regimes (e.g., strong gradients, sub-pixel and multi-pixel shifts). Overall, the work provides a practical, scalable framework that integrates near-field, wave-optics corrections into fast, spectrally-resolved phase retrieval, significantly expanding the operational envelope of PPC-μCT for large, heterogeneous samples.

Abstract

We introduce a fast Eikonal Phase Retrieval (EPR) formulation that accelerates eikonal phase retrieval by more than two orders of magnitude while retaining controlled accuracy. The method is derived from a second-order asymptotic expansion in the propagation distance $L$ and complemented by the leading Wentzel--Kramers--Brillouin (WKB) wave-optics correction, yielding an efficient iterative correction scheme preconditioned by FFT-diagonal, energy-dependent inverse operators (Paganin-type filters). To ensure robustness across practical experimental regimes, we combine two complementary solvers: (i) a local $O(L^2)$ closure that is accurate when eikonal shifts remain sub-pixel, and (ii) a non-local formulation for multi-pixel shifts, in which intensity is propagated through an explicit eikonal ray mapping using a mass-conserving bilinear redisribution on the detector grid, and detector residuals are transferred back to the object grid by the corresponding adjoint (transpose), implemented as bilinear interpolation, before applying an approximate FFT-diagonal preconditioner to accelerate convergence. The same framework supports polychromatic data through a compact spectral discretisation, allowing energy-dependent transport and inversion while keeping the iteration GPU/FFT efficient. Overall, this unified approach enables accurate and computationally efficient phase retrieval across propagation conditions relevant to high-throughput PPC-$μ$CT experiments.

Figures (6)

• Figure 1: Full reconstructed slice of the sheep-head HiP-CT dataset (same specimen and acquisition context as in the original EPR paper). The yellow rectangle marks the region of interest (ROI) used for the detailed comparisons in Figs. \ref{['figgridone']}--\ref{['figgridtwo']}. Among the multiple zoomed regions discussed in the original EPR work, the ROI chosen here corresponds to the brain region (Inset A in Fig. 3 of the original EPR paper), where soft-tissue contrast coexists with strong nearby gradients induced by adjacent bone structures.
• Figure 2: ROI comparison (brain region) for the four forward models. Top row: monochromatic models; bottom row: polychromatic models. Left column: first-order model $L^{1}$; right column: second-order model $L^{2}$. (Here $L^{1}$ mono is equivalent to the Paganin single-distance forward model in the homogeneous-object setting pag2002.)
• Figure 3: Non-local: ROI comparison (brain region) for the four forward models. Top row: monochromatic models; bottom row: polychromatic models. Left column: first-order model $L^{1}$; right column: second-order non-local model $L^{2}$. (Here $L^{1}$ mono is equivalent to the Paganin single-distance forward model in the homogeneous-object setting pag2002.)
• Figure 4: Local vs non-local phase retrieval on the same ROI. Left: local ($L^{2}$) solver (oversampling factor $2$). Centre: non-local solver with oversampling factor $2$, showing residual high-frequency noise. Right: non-local solver with oversampling factor $4$, where this noise is strongly reduced.
• Figure 5: Convergence in the polychromatic case. A large improvement is already achieved when switching from $L^{1}$ poly to $L^{2}$ poly with a single iteration. Subsequent iterations produce only minor (often barely visible) refinements, indicating rapid convergence of the $L^{2}$ polychromatic inversion in this ROI.