🔬 optics

Fast Eikonal Phase Retrieval for High-Throughput Beamlines

This paper introduces a fast, unified Eikonal Phase Retrieval framework that accelerates computation by over two orders of magnitude through a hybrid local/non-local solver and FFT-based preconditioning, enabling accurate, energy-dependent phase retrieval for high-throughput beamline experiments.

Original authors: Alessandro Mirone, Theresa Urban, Joseph Brunet, Claire L. Walsh, Peter D. Lee, Paul Tafforeau

Published 2026-02-02

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Alessandro Mirone, Theresa Urban, Joseph Brunet, Claire L. Walsh, Peter D. Lee, Paul Tafforeau

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to take a perfect photograph of a complex, 3D object (like a sheep's head or a piece of bamboo) using X-rays. In a standard setup, the X-rays pass through the object and hit a detector. However, to see the soft, delicate details that don't absorb much light, scientists move the detector further away. This creates a "shadow" effect where the light bends and interferes, creating a blurry, distorted image full of streaks and artifacts.

To fix this, scientists use a mathematical trick called "Phase Retrieval" to reverse-engineer the original object from the distorted shadow. For years, the best method was like using a simple, linear map: it worked well for gentle hills but failed miserably when the terrain had steep cliffs (strong gradients), causing the map to tear apart and produce garbage images.

This paper introduces a super-fast, upgraded version of that map. Here is how it works, broken down into everyday concepts:

1. The Problem: The "Too-Long" Shadow

When X-rays travel a long distance to the detector, they don't just travel in straight lines; they bend.

The Old Way (Linear): Imagine trying to describe a winding river by only looking at the water's speed at one specific point. If the river curves sharply, your description fails. This is what older methods did. They assumed the bending was small and gentle. When the object had sharp edges (like bone meeting air), the math broke, creating long, ugly streaks across the image.
The New Way (Eikonal Phase Retrieval): The authors realized they needed a better map that accounts for the river curving sharply. They developed a new formula that includes a "second-order" correction. Think of it as upgrading from a flat, 2D paper map to a 3D topographical model that understands steep slopes and sharp turns.

2. The Two Tools in the Toolbox

The authors realized that sometimes the "bending" of the X-rays is tiny (less than the size of a single pixel on the camera), and sometimes it's huge (jumping over several pixels). They built a system with two different tools to handle these situations:

Tool A: The Local Solver (For Gentle Bends)
- Analogy: Imagine you are walking on a smooth path. You can predict your next step by just looking at the ground right under your feet.
- How it works: If the X-rays only shift a tiny bit, the computer uses a quick, local calculation. It's like doing a quick mental math problem. This is incredibly fast.
- Speed: It's so fast that it can process thousands of images in minutes, whereas the old method took hours or even days.
Tool B: The Non-Local Solver (For Wild Bends)
- Analogy: Imagine you are walking on a path that suddenly jumps over a canyon. You can't just look at your feet; you have to look at where you are landing on the other side.
- How it works: When the X-rays bend so much that they jump over multiple pixels, the "local" math fails. This tool explicitly tracks where every single ray of light lands on the detector, like a courier delivering a package to a specific house rather than just dropping it in the neighborhood. It then works backward to figure out where the package came from.
- Robustness: This tool is slower but much stronger. It doesn't break when the image has extreme distortions.

3. The "Polychromatic" Upgrade

Real X-ray beams aren't just one color (energy); they are a mix, like white light containing all colors of the rainbow.

The Old Problem: Previous methods often pretended the beam was just one color to save time, which introduced subtle errors.
The New Solution: This new method treats the beam as a mix of different colors simultaneously. It calculates how each "color" of X-ray bends differently and then combines them. It's like having a choir sing in harmony rather than a single voice, resulting in a much clearer, more accurate picture of the object.

4. The Result: Speed Meets Accuracy

The most impressive claim of the paper is the speed.

Before: Using the original version of this advanced method on a large dataset (like a sheep's head) took about 870 minutes (nearly 15 hours) on a powerful computer.
Now: With this new "Fast Eikonal" method, the same job takes about 1.5 minutes on a single graphics card.
The Takeaway: They made the process roughly 580 times faster. This means scientists can now use this high-accuracy method in real-time, high-throughput workflows (like scanning hundreds of samples a day) without waiting days for the computer to finish the math.

Summary

The authors have built a "smart, fast, and flexible" calculator for X-ray imaging. It knows when to use a quick shortcut (for simple bends) and when to do the heavy lifting (for complex bends). It handles the full spectrum of X-ray colors and, most importantly, it does it so fast that it turns a previously impractical, slow process into a routine, everyday tool for seeing the invisible details of complex 3D objects.

Technical Summary: Fast Eikonal Phase Retrieval for High-Throughput Beamlines

Problem Statement
Propagation-based phase-contrast micro-tomography (PPC-µCT) at fourth-generation synchrotron sources enables high-sensitivity imaging of weakly refracting features by utilizing long sample-to-detector propagation distances. However, this regime introduces significant challenges for image reconstruction. Standard single-distance phase retrieval methods (e.g., Paganin filters) rely on linearized transport-of-intensity assumptions. In heterogeneous specimens containing both weak soft tissues and strong interfaces (e.g., bone–air or bone–fluid), steep phase gradients generate strong intensity modulations that violate these linear assumptions. This leads to characteristic nonlinear artifacts, such as long-range streaks, which degrade image quality and compromise downstream analysis.

While full Fresnel wave-optics formulations can model these effects, they require explicit pixel-wise representation of the complex wavefield. This imposes strict sampling constraints; when phase variations are rapid (as caused by sharp interfaces), standard Fresnel propagation becomes computationally prohibitive due to aliasing and massive memory requirements. Conversely, the previously introduced Eikonal Phase Retrieval (EPR) algorithm successfully modeled strong refraction using an eikonal approximation but suffered from high computational costs (tens of hours for large datasets), rendering it incompatible with the high-throughput acquisition rates of modern beamlines.

Methodology
The authors introduce a "Fast Eikonal Phase Retrieval" (Fast EPR) formulation that accelerates the original EPR algorithm by more than two orders of magnitude while retaining controlled accuracy. The methodology is built on a second-order asymptotic expansion of the Fresnel propagator in the propagation distance $L$ , complemented by the leading Wentzel–Kramers–Brillouin (WKB) wave-optics correction.

The approach employs a unified framework with two complementary solvers to handle different physical regimes:

Local Solver (Sub-pixel shifts): For regimes where the eikonal-induced transverse ray shift $\Delta(r)$ is less than one pixel ( $\eta < 1$ ), the method uses a local $O(L^2)$ closure. This involves a second-order differential expansion of the intensity transport, including the quadratic eikonal term (WKB0) and the leading diffraction-pressure term (WKB1). This local model is preconditioned by FFT-diagonal, energy-dependent inverse operators (Paganin-type filters) to accelerate convergence.
Non-local Solver (Multi-pixel shifts): When shifts exceed one pixel ( $\eta > 1$ ), the local Taylor expansion breaks down. The method switches to a non-local formulation where intensity is propagated through an explicit eikonal ray mapping. This is implemented as a mass-conserving bilinear redistribution on the detector grid. The detector residuals are transferred back to the object grid using the corresponding adjoint (transpose) operator, implemented as bilinear interpolation, before applying the FFT-diagonal preconditioner.

Polychromatic Support
The framework supports polychromatic data through a compact spectral discretization. The spectrum is modeled as an incoherent sum of $N_s$ spectral components. All per-line quantities (phase/absorption parameters and shifts) are treated energy-dependently, allowing for energy-dependent transport and inversion while maintaining FFT efficiency.

Key Contributions

Algorithmic Acceleration: The new formulation reduces the computational cost of phase retrieval by approximately 580 times (2.8 orders of magnitude) compared to the original EPR implementation when normalized to the same GPU resources.
Hybrid Solver Strategy: The introduction of a dual-solver approach (local $O(L^2)$ for sub-pixel shifts and non-local ray mapping for multi-pixel shifts) ensures robustness across a wide range of experimental conditions, from weakly refracting soft tissues to strong interfaces.
Second-Order Physics: By retaining the complete $O(L^2)$ content of the WKB expansion (including the diffraction-pressure term), the method suppresses nonlinear streak artifacts that persist in first-order linear models.
Polychromatic Efficiency: The method enables accurate reconstruction of polychromatic data without sacrificing the computational efficiency required for high-throughput workflows.

Results
The authors validated the method using a sheep-head HiP-CT dataset and a bamboo sample stress test:

Artifact Reduction: Comparisons on the sheep-head dataset demonstrated that the second-order ( $L^2$ ) model significantly reduces residual nonlinear propagation artifacts compared to first-order ( $L^1$ ) models, providing cleaner rendering of soft-tissue contrast near strong gradients.
Convergence: The algorithm converges rapidly; the majority of the improvement is achieved in the first iteration once the second-order physics are enabled.
Robustness: On the bamboo sample, where ray shifts exceeded the detector pixel size, the local solver failed (diverging after the first iteration), whereas the non-local solver remained stable and accurate.
Performance: For a dataset of 8,000 radiographs ( $256 \times 3104$ pixels), the Fast EPR method completed phase retrieval in approximately 1.1 to 1.9 minutes on a single GPU (for monochromatic data) and 1.47 to 1.9 minutes for polychromatic data (5 spectral lines). This contrasts with the original EPR benchmark, which required approximately 870 minutes for a comparable volume.

Significance
The paper claims that this unified approach enables accurate and computationally efficient phase retrieval across propagation conditions relevant to high-throughput PPC-µCT experiments. By reducing reconstruction times to be compatible with beamline acquisition rates (one sample every few hours), the method facilitates the routine use of nonlinear phase retrieval in high-throughput workflows. Furthermore, the availability of a fast second-order near-field model with explicit wave-optics correction allows researchers to operate closer to the near-field boundary and explore regimes with stronger phase gradients (e.g., lower energy or higher resolution) where standard linear approaches fail, without incurring the prohibitive costs of full Fresnel propagation. The authors note that the algorithmic complexity of the new method allows it to be hidden behind tomographic back-projection in concurrent hardware setups, further integrating it into practical reconstruction pipelines.