Propagation-Based Phase Contrast Imaging

• Propagation-based phase contrast imaging is a technique that exploits free-space Fresnel diffraction to transform phase shifts into measurable intensity modulations.
• It employs deterministic and learning-based phase retrieval methods, such as the Paganin approach, to improve signal-to-noise ratio and enable quantitative imaging.
• The method spans X-ray and neutron modalities, relying on coherence and model assumptions to balance spatial resolution, noise, and material-specific contrast.

Propagation-based phase-contrast imaging is the in-line, free-space-propagation form of phase-contrast imaging in which a coherent wavefield traverses a specimen, acquires attenuation and phase shifts, and is then recorded after a finite propagation distance so that Fresnel diffraction converts otherwise invisible phase variations into measurable intensity modulations (Paganin et al., 2019). In X-ray implementations, the sample is described by the complex refractive index $n=1-\delta+i\beta$, with $\delta$ governing phase shift and $\beta$ absorption; in neutron implementations, equivalent formulations often use the bound coherent scattering length $b$ and total cross section $\sigma$ directly in the forward and retrieval models (Paganin et al., 2019, Østergaard et al., 2022). The subject spans radiography and tomography, deterministic and learned phase retrieval, spectral and material-selective inversion, neutron and X-ray variants, and a substantial mathematical literature on uniqueness, stability, detector truncation, and model validity.

1. Physical basis and image formation

Propagation-based phase contrast is generated by free-space propagation itself. Under the projection approximation, a specimen modifies the exiting wave by projected attenuation and projected phase. For X-rays, the projected phase shift and attenuation are written as

$\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$

$$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$

The propagated field is then

$$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$

with Fresnel propagation operator

$$\mathcal{D}_{\Delta} = e^{ik\Delta}\, \mathcal{F}^{-1} \exp\left[ \frac{-i\Delta(k_x^2+k_y^2)}{2k} \right] \mathcal{F}.$$

The corresponding transport-of-intensity equation is

$$-\nabla_{\perp}\cdot\left[I(x,y,z)\nabla_{\perp}\phi(x,y,z)\right] = k\frac{\partial I(x,y,z)}{\partial z},$$

which makes explicit that longitudinal intensity evolution is driven by transverse phase structure (Paganin et al., 2019).

This mechanism is especially useful for weakly attenuating materials, because propagation converts phase curvature and refraction into edge-enhanced intensity structure. In the small-defocus regime,

$\delta$0

so interfaces and projected density curvature become visible even when attenuation contrast is weak (Paganin et al., 2019).

Spatial coherence is a necessary experimental condition. In the X-ray tutorial literature, finite source size produces penumbral blurring

$\delta$1

so increasing object-to-detector distance strengthens propagation contrast but also increases source-induced blur. In neutron imaging, the same trade-off is expressed through beam divergence $\delta$2 and blur width $\delta$3; stronger collimation improves phase sensitivity but reduces flux (Paganin et al., 2019, Paganin et al., 2019).

Propagation-based imaging is therefore distinct from analyzer-based, grating, or edge-illumination methods in both hardware and transfer characteristics. Within the TIE viewpoint, propagation-based methods are associated mainly with the $\delta$4 term, whereas gradient-based methods are linked to $\delta$5 (Paganin et al., 2019).

2. Deterministic phase retrieval and homogeneous-object models

The most influential deterministic inversion is the single-distance homogeneous-object method of Paganin et al. For a single material of projected thickness $\delta$6, the tutorial gives

$\delta$7

The denominator is a Lorentzian low-pass filter, which explains why the method is fast, deterministic, and noise robust, but also resolution reducing (Paganin et al., 2019).

In multi-material settings, the homogeneous assumption becomes restrictive. A widely used two-material generalization is written as

$\delta$8

This reconstructs the chosen interface correctly, but only that interface. In samples containing weakly attenuating material adjacent to dense material, retrieval tuned for the weak material gives strong denoising but can over-blur the dense boundary; tuning to the dense boundary preserves the interface but leaves more remnant noise (Pollock et al., 2023).

A direct remedy is 3D Masked Phase Retrieval (3DMPR), a non-iterative volumetric workflow for samples such as brain tissue enclosed by skull. It segments the high-$\delta$9 component in a boundary-correct reconstruction, replaces those voxels in the raw phase-contrast CT volume by the theoretical attenuation coefficient of the low-$\beta$0 material, applies strong 3D homogeneous retrieval to the modified volume, and then interleaves the high-$\beta$1 regions back from the interface-specific reconstruction. In the rabbit kitten brain dataset, this yielded a 6.9-fold improvement in brain-tissue SNR compared to the standard phase-retrieval procedure without over-smoothing. In the aluminium-water phantom, it provided a 4.2-fold SNR boost while preserving the boundary resolution at $\beta$2, compared to $\beta$3 in conventional phase retrieval (Pollock et al., 2023).

A different deterministic extension addresses unknown multi-material samples by exploiting interface distortions produced by incorrect retrieval parameters. In that approach, the reconstructed interface is fitted by an error-function model plus a Gaussian-weighted correction term,

$\beta$4

and the fitted coefficient $\beta$5 is used to infer the correct interface ratio $\beta$6. Applied to breast tissue, the method recovered the refractive component $\beta$7 with 0.6–2.4% accuracy compared to theoretical values while remaining single-distance, with one exposure per tomographic angle (Alloo et al., 2021).

3. Tomographic formulations and computational inversion

Propagation-based CT inherits the projection-space forward model but introduces additional choices about where retrieval is applied. In conventional pipelines, one retrieves each projection and then reconstructs. A mathematically equivalent alternative is post-reconstruction 3D retrieval. For homogeneous TIE-Hom PB-CT, the reconstructed refractive decrement can be written as

$\beta$8

where $\beta$9 is filtered backprojection and $b$0. This allows phase retrieval to be applied directly in the reconstructed 3D volume, and even only to localized regions of interest. In multi-material PB-CT, that avoids repeated 2D retrieval and repeated full CT reconstruction for each material choice of $b$1 (Thompson et al., 2018).

The associated 3D TIE-Hom filter has Fourier transfer function

$b$2

and real-space kernel

$b$3

This makes explicit that post-reconstruction phase retrieval is a 3D low-pass inversion whose support must be accommodated by the chosen ROI (Thompson et al., 2018).

A second line of work replaces continuous Fourier calculus by a discrete derivation matched to pixelated detectors. The generalized Paganin method substitutes the discrete Laplacian symbol

$b$4

for the continuous $b$5, giving

$b$6

The point is not a new physical model but a better numerical match to sampled data. With a direct photon-counting detector, this yielded a $b$7 spatial-resolution improvement for single-material retrieval, $b$8 for two-material retrieval, and up to $b$9 in one dual-energy electron-density reconstruction (Pollock et al., 2021).

A third computational development is end-to-end learning-based quantitative phase retrieval under laboratory conditions. One reported implementation used a U-Net with two input channels, $\sigma$0 and $\sigma$1, one output channel, image size $\sigma$2, depth 3, $\sigma$3 kernels in encoder and decoder blocks, and $\sigma$4 kernels in skip connections. The target was the projected phase map at the contact plane corresponding to the peak energy of the spectrum, 24 keV. The method was trained on simulated data including polychromaticity, partial spatial coherence, blur, cone-beam geometry, and mixed Poisson-Gaussian noise; with correct flat-field scaling it was robust to at least 20% variation in propagation distance, and in multi-material phantoms it outperformed conventional analytic retrieval in NRMSE, SSIM, and PSNR (Deshpande et al., 2022).

4. Spectral, material-selective, and auxiliary contrast extensions

Spectral propagation-based imaging couples phase retrieval to material decomposition by using basis-material optical properties directly. For multi-energy propagated radiographs, the linearized TIE with a basis-material model gives

$\sigma$5

or, in matrix form, $\sigma$6. Solving

$\sigma$7

produces projected thickness maps $\sigma$8 of the selected materials directly from propagated images. In the two-material Al/PMMA simulations, using $\sigma$9 cm decreased $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$0 by a factor of approximately 3.15 relative to $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$1, while increasing $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$2 by a factor of 1.5. Experimentally, clean decomposition of a pure aluminium coin and PMMA spheres was obtained without residual phase contrast effects, and the same framework was used to generate a bone-free projected lung-air image from rabbit kitten chest radiographs acquired at 24 and 34 keV (Schaff et al., 2019).

Propagation-based data can also be mined for auxiliary contrast channels. A dark-field extraction method decomposes the propagated intensity into a TIE/Hom bright-field component and a higher-order diffraction component associated with weak Born perturbations. Under the homogeneous-object assumption $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$3, the object-plane intensity is recovered through the homogeneous Fresnel propagator

$\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$4

and the residual higher-order term is interpreted as a dark-field or SAXS-sensitive signal. In the homogeneous case this requires a single propagation-based image in 2D and one image per projection angle in CT. A preliminary breast CT example showed improved visualization of microcalcifications and resolution of clip structure in dark-field reconstructions (Gureyev et al., 2020).

These extensions change the role of phase retrieval. Instead of merely recovering a denoised thickness or phase map, they treat propagation as part of a joint inversion for materials, scattering signatures, or both. The recurring requirement is additional structure: a basis-material model, homogeneity, weak Born terms, or multiple energies.

5. Neutron propagation-based phase-contrast imaging

Neutron PBI follows the same propagation logic but uses neutron interaction parameters. In the homogeneous single-material model, the exit-surface intensity and phase are

$\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$5

with projected density $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$6. The TIE-based propagated intensity then becomes

$\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$7

which is inverted by a Paganin-type filter (Paganin et al., 2019).

A proof-of-concept neutron tomography experiment at ISIS/IMAT showed that this combination of propagation contrast and phase retrieval can overcompensate the flux loss due to collimation. The measured retrieval gain was $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$8 in SNR, corresponding to a net SNR gain of $\Delta\phi(x,y)=-k\int \delta_\omega(x,y,z)\,dz,$9 relative to conventional attenuation-based tomography after accounting for the collimation penalty, and to an effective brilliance gain of $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$0. The same paper derived the collimation condition

$$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$1

and the optimal divergence

$$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$2

which summarize the trade-off between flux and phase sensitivity (Paganin et al., 2019).

A later neutron study demonstrated that single-distance Paganin-style retrieval also remains practical with pink and white beams. For a polychromatic beam, the retrieval is written in terms of spectrum-averaged quantities, but in the reported experiments it was functionally almost identical to applying the monochromatic filter at the weighted mean wavelength. In a weakly absorbing Al/Zr near-pure-phase sample, no discernible signal was seen until about $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$3, whereas by $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$4 the sample was clearly visible in raw propagated images. In human cortical bone, phase retrieval significantly improved signal-to-noise ratio and interpretability, and in the bone plus $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$5 experiment it enabled separation of bone and $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$6 in a correlative imaging workflow. The practical conclusion was that propagation-based neutron phase contrast remains effective under pink and white beams, which is especially important for time-resolved and in situ experiments (Østergaard et al., 2022).

Neutron PB-PCI is therefore complementary to X-ray PB-PCI rather than derivative of it. The contrast mechanism is still free-space propagation of wavefront distortions, but the material dependence follows neutron scattering length and cross section rather than electron density, which makes isotope-sensitive problems such as $$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$7 transport particularly attractive (Østergaard et al., 2022).

6. Identifiability, model fidelity, and optimization

The inverse problem of PBI is not fundamentally ambiguous in the near field under the projection and paraxial approximations. A uniqueness theorem shows that a compactly supported complex object illuminated by a known plane wave or Gaussian beam is uniquely determined by a single intensity measurement at one detector distance and one wavelength. The key reason is that the unscattered part of the illumination acts as a known reference wave, so the measured Fresnel pattern is not ordinary modulus-only Fourier data (Maretzke, 2014).

Uniqueness, however, does not imply uniform stability. For the linearized contrast-transfer-function model with compact support, a stability estimate of the form

$$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$8

holds for arbitrary complex objects, but the constant decays nearly exponentially with Fresnel number,

$$\log_e\left[\frac{I_\omega(x,y,z=z_0)}{I_\omega(x,y,z=0)}\right] = -\int \mu_\omega(x,y,z)\,dz, \qquad \mu_\omega=2k\beta_\omega.$$9

For homogeneous objects and for two-distance measurements, the dependence becomes algebraic rather than nearly exponential, which is why monomorphous assumptions and multi-distance data remain so valuable in practice (Maretzke et al., 2016).

Finite detector size introduces an additional instability that is independent of the phase problem. Fresnel propagation has global support, so any compactly supported object produces a propagated field that is nonzero everywhere. With a finite field of view, reconstruction from truncated data is severely ill-posed even if the complex wavefield were known. Quantitative locality estimates nevertheless show a sharp resolution limit: the smallest resolvable length scale is $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$0 times the detector’s aspect length, where $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$1 is the Fresnel number associated with that scale. The limit is also position dependent, becoming more restrictive near detector boundaries; for real-valued or homogeneous objects the stable region is substantially larger except near corners (Maretzke, 2018).

Model mismatch enters at several levels. The projection approximation, ubiquitous in PB-CT, can fail when internal propagation inside the object becomes significant. In simulated synchrotron micro-CT of a 5 mm zebrafish phantom at 20 keV and 50 mm propagation, reconstructed images generated with the projection approximation and with a 10-slice multislice model showed little difference at $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$2 detector resolution, corresponding to $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$3, but a noticeable difference at $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$4 resolution, corresponding to $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$5, especially around fine details. The implication is that accurate phase retrieval at sub-micron resolution may require forward models that avoid the projection approximation (Jadick et al., 17 Aug 2025).

Optimization studies make these trade-offs quantitative. For 2D in-line imaging of a homogeneous edge feature, the magnification that maximizes contrast is

$$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$6

whereas the magnification that maximizes the biomedical quality metric is

$$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$7

In the 3D case of propagation-based phase-contrast tomography with Paganin phase retrieval, the optimal magnification $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$8 is close to $$\psi_\omega(x,y,z=\Delta)=\mathcal{D}_\Delta \psi_\omega(x,y,z=0),$$9. For the IMBL-like breast example, the reported optima were $$\mathcal{D}_{\Delta} = e^{ik\Delta}\, \mathcal{F}^{-1} \exp\left[ \frac{-i\Delta(k_x^2+k_y^2)}{2k} \right] \mathcal{F}.$$0 in 2D and $$\mathcal{D}_{\Delta} = e^{ik\Delta}\, \mathcal{F}^{-1} \exp\left[ \frac{-i\Delta(k_x^2+k_y^2)}{2k} \right] \mathcal{F}.$$1 in 3D; the same analysis gave practical transmission criteria of about $$\mathcal{D}_{\Delta} = e^{ik\Delta}\, \mathcal{F}^{-1} \exp\left[ \frac{-i\Delta(k_x^2+k_y^2)}{2k} \right] \mathcal{F}.$$2 bulk transmission for the 2D biomedical metric and about $$\mathcal{D}_{\Delta} = e^{ik\Delta}\, \mathcal{F}^{-1} \exp\left[ \frac{-i\Delta(k_x^2+k_y^2)}{2k} \right] \mathcal{F}.$$3 for the 3D tomographic case (Gureyev et al., 2 Oct 2025, Gureyev et al., 25 Jun 2025).

Taken together, these results locate the central tension of propagation-based phase-contrast imaging. Near-field propagation can encode materially stronger contrast than absorption alone, and under the right priors that encoding can be inverted with large SNR gains. At the same time, identifiability depends on known references and support, stability depends on object class and detector extent, and quantitative fidelity depends on whether approximations such as homogeneity, projection, and near-Fresnel transport remain valid at the desired spatial scale.