Source-linked AI summary
Bayesian Joint Reconstruction and Uncertainty Quantification for Limited-Data Multi-Contrast X-Ray CT
Hassan Yazdanian, Henrik Mäkinen
TL;DR
Limited-data acquisition reduces scan costs but degrades the data, and posterior UQ for jointly reconstructed attenuation, phase, and dark-field images remains comparatively little explored. The study combines Bayesian reconstruction from separately retrieved sinograms with channel-wise TV and a JTV prior, and uses NUTS sampling to estimate pixel-wise uncertainty. JTV lowers experimental mean dark-field posterior standard deviation by 23–30% and raises median per-pixel ESS by 4–38% under limited data; it modestly improves synthetic accuracy and cleans experimental dark-field background.
Problem
Limited-data acquisition reduces scan costs but degrades the data, and posterior UQ for jointly reconstructed attenuation, phase, and dark-field images remains comparatively little explored.
Source: §I · pp. 1–2Method
The study combines Bayesian reconstruction from separately retrieved sinograms with channel-wise TV and a JTV prior, and uses NUTS sampling to estimate pixel-wise uncertainty.
Source: §I · p. 2Results
JTV lowers experimental mean dark-field posterior standard deviation by 23–30% and raises median per-pixel ESS by 4–38% under limited data; it modestly improves synthetic accuracy and cleans experimental dark-field background.
Source: §IV · p. 13, §V · p. 15Takeaways & Limitations
Both priors suppress undersampling artifacts relative to FBP, while JTV's clearest reported benefit is a cleaner experimental dark-field background.
Source: §V · p. 15Synthetic posterior spread understates boundary errors, and phase chains mix slowly under angular undersampling.
Source: §IV · p. 13, §V · p. 15
Abstract
from arXiv · showhide
Multi-contrast X-ray computed tomography (CT) combines attenuation contrast with phase and dark-field contrast, which are sensitive to refraction and ultra-small-angle scattering, respectively. Common implementations record several intensity measurements per projection angle by phase or mask stepping, which lengthens acquisition and, at a fixed exposure per measurement, raises the dose. Limited-data acquisition via sparse angular sampling or reduced stepping lowers these costs but makes the data less informative, calling for prior information and uncertainty quantification. We study joint Bayesian reconstruction of attenuation, phase, and dark-field images from separately retrieved sinograms. A joint total variation term in the prior couples edge locations across the images but not their values. Posterior sampling with the No-U-Turn Sampler complements maximum a posteriori estimates with pixel-wise uncertainty. We compare this joint prior with channel-wise total variation alone on synthetic and experimental Talbot-Lau interferometry data under full, sparse-angle, sparse-step, and combined acquisition (up to 24 times fewer exposures). Both maximum a posteriori reconstructions suppressed undersampling artifacts, with synthetic relative errors 1.6-11 times lower than those of filtered backprojection. The joint prior gave modest synthetic accuracy gains, a cleaner experimental dark-field background, 23-30 % lower mean dark-field posterior standard deviation, and 4-38% higher median effective sample size in experimental limited-data runs. In the synthetic study, posterior standard deviations matched or exceeded the errors in flat regions but underestimated them at boundaries, especially in phase. Experimental phase chains mixed slowly under angular undersampling. The formulation also applies to edge illumination with adapted forward and noise models.
I. INTRODUCTION
Multi-contrast CT combines attenuation, phase, and dark-field information, but limited-data acquisition makes reconstruction less informative and leaves questions about shared priors and uncertainty. This work combines Bayesian reconstruction and posterior sampling for separately retrieved sinograms, using a prior that favors shared edge locations while allowing contrast values to differ.
I. INTRODUCTION
- Sparse angular sampling and reduced stepping lower acquisition costs but produce underdetermined reconstructions and noisier retrieved sinograms, respectively.Regularization introduces prior information, but a regularized image alone does not show which features the measurements support.Source: §I · p. 1
- Shared edges motivate structural coupling, while prior approaches either regularize channels separately or impose prescribed relationships between channel values.Structural regularizers instead couple edge locations or orientations while leaving image values free.Source: §I · pp. 1–2
- Posterior UQ remains comparatively little explored for jointly reconstructed attenuation, phase, and dark-field images.The posterior characterizes uncertainty under assumed models and complements point estimates with pixel-wise uncertainty maps.Source: §I · p. 2
II. METHODS
- The proposed Bayesian reconstruction combines channel-wise TV with smoothed JTV and uses NUTS sampling to complement MAP estimates with pixel-wise uncertainty.The study compares the priors across full, sparse-angle, sparse-step, and combined acquisition using synthetic and experimental TLI data.Source: §I · p. 2
A. Forward Problem
- The forward problem predicts attenuation, differential phase, and dark-field sinograms from parameter images and known acquisition geometry.The DPC map applies a detector-coordinate derivative in the sensitivity direction; the dark-field model assumes visibility decays exponentially along each ray.Source: A. Forward Problem · p. 2
B. Bayesian Inverse Problem
The Bayesian inverse problem estimates cross-sectional parameters from noisy multi-contrast observations. It models parameters and noise as random vectors and uses their assumed independence to form the likelihood.
B. Bayesian Inverse Problem
- The inverse problem estimates cross-sectional parameters from noisy multi-contrast observations.Source: B. Bayesian Inverse Problem · p. 3
- The Bayesian formulation models the parameters and noise as random vectors, with observed data a realization of the resulting random data.Source: B. Bayesian Inverse Problem · p. 3
- Assuming noise is independent of the parameters, the additive observation model yields the likelihood from the prior and noise densities.Source: B. Bayesian Inverse Problem · p. 3
1) Likelihood:
The likelihood models the three sinogram channels with independent Gaussian noise. The prior combines channel-wise smoothed TV with a joint term that favors co-located edges without tying contrast values.
1) Likelihood:
- The attenuation, DPC, and dark-field noise blocks are modeled as mutually independent Gaussian vectors with known channel variances.The likelihood uses the Euclidean norm.Source: 1) Likelihood · p. 3
2) Prior:
- The channel-wise prior applies smoothed TV to scaled images, with channel weights setting penalty strength and an amplitude penalty making the prior proper.The images are scaled because the channels differ in physical units and magnitude.Source: 2) Prior · p. 3
- JTV groups gradients across channels at each pixel, favoring co-located edges without tying the size or sign of their jumps.The coupling weight αJ controls this term; αJ = 0 recovers the channel-wise prior.Source: 2) Prior · p. 3
3) Posterior:
The negative log-posterior combines the negative log-likelihood and the negative log-prior. The channel-wise terms factorize across channels, so only JTV couples them; strong convexity gives each posterior a unique mode.
3) Posterior:
- The negative log-posterior is the negative log-likelihood plus the negative log-prior, up to an additive constant.Source: 3) Posterior · p. 4
- The likelihood and channel-wise prior factorize over channels, so only the JTV term couples them.Source: 3) Posterior · p. 4
- The quadratic amplitude penalty makes the objective strongly convex, yielding log-concave posteriors with a unique mode.Source: 3) Posterior · p. 4
C. Inference
The inference procedure obtains MAP images by minimizing the negative log-posterior, then uses posterior sampling to quantify pixel-wise uncertainty. Synthetic and experimental reconstructions are evaluated across full and undersampled acquisitions.
1) MAP Estimate:
- MAP estimation minimizes the negative log-posterior, making the reconstruction a variational regularization problem whose prior weights act as regularization parameters.The optimization method and parameter settings are specified in the numerical experiments.Source: 1) MAP Estimate · p. 4, §III · p. 4
2) Sampling and UQ:
- NUTS post-warmup draws from independent chains are used to estimate the posterior mean and pixel-wise standard deviation.Uncertainty is conditional on the assumed forward model, noise model, and prior; convergence and sampling efficiency are assessed per pixel with R-hat and ESS.Source: 2) Sampling and UQ · p. 4
III. NUMERICAL EXPERIMENTS AND RESULTS
- The experiments compare FBP, channel-wise TV, and JTV on synthetic and experimental data under full, sparse-angle, sparse-step, and combined undersampling.Synthetic MAP and posterior UQ are evaluated against ground truth; experimental evaluation uses data consistency and full-data reconstructions as surrogate references.Source: §III · p. 4
- The reported reconstructions and posterior sampling use an NVIDIA RTX 6000 Ada GPU, and the code and experimental data are available online.The software stack includes Python, PyTorch, ASTRA Toolbox, and Pyro.Source: §III · p. 4
A. Synthetic Phantom
The synthetic study uses a multi-contrast phantom designed to test shared and channel-specific features, cross-talk, and uncertainty on small or low-contrast structures. Simulation and reconstruction settings define the noise, acquisition, and parameter tuning used for evaluation.
1) Phantom and Data Simulation:
- The phantom tests shared boundaries across channels, channel-specific contrast, cross-talk, and recovery of small or low-contrast features for uncertainty analysis.Inclusion F is dark-field dominant, P differs from the disc only in phase, and the phantom also contains a low-contrast ellipse, decreasing-size discs, and a one-pixel-wide dark-field crack.Source: 1) Phantom and Data Simulation · p. 4, Figure 1 · p. 5
- Full synthetic data use 180 views and 10 phase steps, while sparse settings use 30 views and 4 steps.The phantom is rendered at 1024^2, block-averaged to 512^2 for simulation, and reconstructed on a 256^2 grid; a common image scaling keeps max |DPC| at or below 0.8π.Source: 1) Phantom and Data Simulation · p. 4
- Independent zero-mean Gaussian noise is added to the retrieved attenuation, DPC, and dark-field sinograms.The supplied passages describe channel-specific relative noise levels and scaling by the number of phase steps.Source: 1) Phantom and Data Simulation · pp. 4–5
2) Reconstruction and Evaluation Settings: √
- Prior parameters are held fixed across acquisition conditions, with weights chosen by averaging relative-error performance over all four conditions.Tuning against ground truth on the evaluated noise realization gives best-case performance for each prior; tuning only on full data would over-regularize undersampled cases.Source: 2) Reconstruction and Evaluation… · p. 5
- Posterior sampling uses two NUTS chains of 1,000 post-warmup samples after 500 warmup iterations, with 147,192 active dimensions.Images are scored against ground truth using relative error, SSIM, and PSNR.Source: 2) Reconstruction and Evaluation… · p. 5
3) MAP Reconstructions:
Bayesian reconstructions substantially reduce synthetic undersampling error relative to FBP, while JTV provides modest accuracy gains and improves dark-field structure. Posterior uncertainty increases with undersampling, but errors at boundaries—especially in phase—are underestimated, and sparse-angle sampling can mix slowly.
3) MAP Reconstructions:
- 1.6–11 times lower synthetic relative errors than FBP are achieved by Bayesian reconstructions under combined undersampling, with channel REs of 0.066–0.071 in µ and δ and 0.111–0.146 in ϵ.FBP RE reaches 0.425 in µ and 1.56 in ϵ; its phase RE is 0.170.Source: 3) MAP Reconstructions · p. 5, Figure 2 · p. 6, Table I · p. 7
- JTV lowers RE in 11 of 12 channel–condition pairs, with the exception of δ under combined undersampling (+0.4 %).Gains are modest, reaching 3.1 % in full-data ϵ; dark-field SSIM increases by 0.008–0.013 in every condition, and channel-specific feature contrast is not transferred.Source: 3) MAP Reconstructions · pp. 5–6, Table I · p. 7
4) Posterior Sampling and UQ:
- Sparse-angle sampling leaves some coordinates slowly mixing, with minimum ESS of 5–13 at 30 views compared with 142–311 at 180 views.Maximum R-hat is 1.10–1.17 for 30 views and 1.02–1.04 for 180 views.Source: 4) Posterior Sampling and UQ · p. 7
- From full data to combined undersampling, TV mean posterior standard deviation rises by 30 %, 57 %, and 36 % in µ, δ, and ϵ; JTV reduces mean standard deviation by 12–14 % in µ and 23–27 % in ϵ.JTV changes the δ mean standard deviation by about 2 % relative to TV.Source: 4) Posterior Sampling and UQ · p. 7
- At edges, RMS error-to-standard-deviation ratios are 2.5–3.2 for µ, 5.3–10.5 for δ, and 1.2–1.8 for ϵ, indicating underestimated errors at boundaries.In flat regions, uncertainty is conservative for µ and ϵ and roughly matched for δ; the reported diagnostics are per-pixel marginals and do not assess joint posterior coverage.Source: 4) Posterior Sampling and UQ · p. 7, Table II · p. 7
B. Real Data: Talbot–Lau Interferometry
The experimental TLI study uses a sugar–water phantom and a multi-step, multi-angle acquisition, then retrieves attenuation, differential-phase, and dark-field sinograms for full and reduced-data conditions.
1) TLI Setup and Phantom:
- The TLI setup used gratings G0, G1, and G2, with G1 producing a π phase shift at 28 keV and G2 analyzing the interference pattern.G0 provides spatial coherence at the tube window; G1 and G2 are positioned downstream.Source: 1) TLI Setup and Phantom · pp. 7, 9
- Object and reference scans used 10 phase steps at each of 361 projection angles, with 18 s per step and about 18 h of exposure per scan.Total exposure was about 40 h including sample repositioning.Source: 1) TLI Setup and Phantom · p. 9
- The phantom was a 15 mm PMMA rod with three 4 mm cavities containing approximately 10, 20, and 30 m-% sugar–water solutions.It was submerged in water to reduce phase wrapping and beam-hardening artifacts; mean reference visibility was about 15%.Source: 1) TLI Setup and Phantom · p. 9
2) Data Retrieval and Preprocessing:
- Single-harmonic Fourier analysis of the stepping curves estimated mean intensity, visibility, and modulation phase from 10 steps, yielding attenuation, differential-phase, and dark-field sinograms.Attenuation and dark-field are line integrals, while differential phase is the detector-direction derivative of the phase line integral; Bayesian reconstruction uses differential phase directly.Source: 2) Data Retrieval and Preprocessing · p. 9
- Noise levels for each channel were estimated by Monte Carlo propagation through the retrieval pipeline, including nonlinear amplification at finite visibility.The five-step retrieval increased the estimated noise standard deviations by approximately √2 across attenuation, differential phase, and dark-field.Source: 2) Data Retrieval and Preprocessing · p. 9
- Sparse-angle, sparse-step, and combined acquisitions used 12, 2, and 24 times fewer object exposures than the full scan, respectively.Sparse-angle data used 30 projections; sparse-step data used five steps, and exposure time and dose fell by the same factors at fixed tube settings.Source: 2) Data Retrieval and Preprocessing · pp. 9–10
3) Reconstruction and Evaluation Settings:
Experimental reconstructions use retrieved attenuation, differential-phase, and dark-field sinograms within a calibrated fan-beam model. A round-trip test finds discrepancies near the estimated noise scale, with structured residuals near strong gradients.
3) Reconstruction and Evaluation Settings:
- The reconstruction used 842 detector bins and a 512² grid, with a circular support containing the rod and surrounding-water margin; the likelihood used the three retrieved sinograms and their estimated noise levels.The model used a nearest-neighbor projector, while FBP used a Ram–Lak filter.Source: 3) Reconstruction and Evaluation… · p. 10
- MAP estimates used zero-initialized L-BFGS, while posterior sampling used two NUTS chains of 1,000 samples with 200 warmup iterations.Without ground truth, assessment used noise-normalized sinogram residuals, visual agreement with full-data reconstructions, and comparisons of posterior means and standard deviations.Source: 3) Reconstruction and Evaluation… · p. 10
4) Forward-Model Adequacy:
- 0.022 was the relative intensity discrepancy in the round-trip test; attenuation, differential-phase, and dark-field sinogram discrepancies were 0.44, 0.69, and 1.19 versus noise references of 0.28, 0.42, and 0.76.These sinogram discrepancies were approximately 10%, 15%, and 11% above the nominal noise-only reference, with residual structure concentrated near strong gradients.Source: 4) Forward-Model Adequacy · p. 10
- 0.25 was the phase-image discrepancy, with SSIM 0.68 and an approximately 15% uniform offset inside the rod; attenuation and dark-field image comparisons were noise-dominated.Their discrepancies were 1.31 and 1.40, respectively.Source: 4) Forward-Model Adequacy · p. 10
- The round trip tests consistency among measured data, retrieval, geometry, and channel models; it does not establish the FBP images as ground truth.Source: 4) Forward-Model Adequacy · p. 10
5) Prior Parameter Selection:
Prior parameters were selected using held-out sinogram prediction and imperfect reconstruction references. In sparse-angle MAP results, TV and JTV suppressed FBP artifacts, while JTV’s clearest experimental difference was a cleaner dark-field background.
5) Prior Parameter Selection:
- Without ground truth, channel weights were selected by predicting every fifth held-out angle, with relative ℓ2 error of the predicted sinograms as the score.The selected weights were αµ = 5.62, αδ = 10, and αϵ = 1.78; phase held-out error was within 0.42% of its minimum across [10, 31.6].Source: 5) Prior Parameter Selection · p. 10
- β = 0.03 maximized phase SSIM against full-data FBP in sparse-angle data; λ0 = 10−2 was selected to add curvature in weakly constrained directions.Relative to λ0 = 10−6, it moved only the phase MAP 4% further from FBP in relative ℓ2 with unchanged SSIM.Source: 5) Prior Parameter Selection · p. 10
- Held-out error varied by less than 0.1% over αJ ∈ [0, 2.15]; αJ = 1 was selected because mean SSIM stayed nearly constant up to 1 and decreased above it.The comparison used sparse-angle JTV against full-data channel-wise TV as a surrogate reference, and parameter groups were selected separately rather than jointly.Source: 5) Prior Parameter Selection · p. 10
6) MAP Reconstruction Results:
- TV and JTV suppressed severe sparse-angle FBP artifacts and resolved all three cavities in attenuation and phase; JTV’s dark-field background standard deviation was 21–39% lower than TV’s.JTV also produced a more continuous dark-field rod boundary.Source: 6) MAP Reconstruction Results · p. 10
- Sinogram residuals decreased from 1.43–4.06 for FBP to 1.21–1.31 for TV and JTV, with TV residuals smaller throughout and JTV values up to 1.4% higher.Residuals measure data consistency, not reconstruction accuracy; values above unity can reflect prior-induced bias, model discrepancy, or limitations of the noise model.Source: 6) MAP Reconstruction Results · p. 11, Table III · p. 11
7) Posterior Sampling and UQ:
Experimental posterior sampling shows that JTV can reduce dark-field posterior spread and improve sampling efficiency in limited-data settings, while phase mixing remains slow under angular undersampling. The uncertainty maps and efficiency gains remain conditional on the assumed models, and boundary errors can be understated.
7) Posterior Sampling and UQ:
- JTV raises median ESS in every limited-data channel–condition pair, but these gains compare different target posteriors and leave slow phase mixing unresolved.The reported increase is 4–38 %.Source: 7) Posterior Sampling and UQ · p. 13
IV. DISCUSSION
- For experimental dark-field reconstructions, JTV reduces mean posterior standard deviation by 23–30%; across limited-data channel–condition pairs, it raises median ESS by 4–38%.Narrower posteriors alone do not establish better calibration.Source: §IV · p. 13, 7) Posterior Sampling and UQ · p. 13
- Posterior spread can substantially understate boundary error, especially for phase, despite matching or exceeding error in flat regions.The spatial summaries come from one noise realization and are not repeated-sampling coverage estimates.Source: §IV · p. 13
- Retrieved-channel inference assumes independent additive Gaussian noise, approximating nonlinear intensity-noise propagation and neglecting possible correlations.Experimental residuals exceed the estimated noise scale, with round-trip discrepancies about 10–15 % above the nominal independent-noise reference; their source is not uniquely identified.Source: §IV · pp. 13–14
- At fixed tube settings, combined undersampling corresponds to a 24-fold exposure-time and dose reduction, while phase mean posterior standard deviation rises by 89 %.These are retrospective estimates for a water-immersed phantom, without dosimetry or equal-image-quality validation.Source: §IV · p. 15
7) Posterior Sampling and UQ:
- Under combined undersampling, phase mixing remains slow: its 99th-percentile R̂ is approximately 1.18 for both priors, with no minimum-ESS improvement from JTV.Under sparse angles, JTV lowers phase 99th-percentile R̂ from 1.12 to 1.08 and raises minimum ESS from 2.7 to 5.1.Source: 7) Posterior Sampling and UQ · pp. 12–13, §IV · p. 15
IV. DISCUSSION
- A 361-view experimental chain takes about 9.5 h, and slow phase mixing limits quantitative uncertainty summaries.The authors recommend validating Monte Carlo error with longer or additional chains and assessing stability across noise realizations.Source: §IV · p. 15
- Three-dimensional extension will require more scalable inference, alongside intensity-domain noise models, repeated-sampling coverage studies, and further biomedical and industrial datasets.Source: §IV · p. 15
V. CONCLUSION
The study evaluates joint Bayesian reconstruction of attenuation, phase, and dark-field images under limited-data acquisition. JTV improves synthetic MAP accuracy modestly and produces a cleaner experimental dark-field background, but uncertainty calibration and slow phase mixing remain concerns.
- Across synthetic and experimental TLI studies, MAP reconstruction with either prior suppresses undersampling artifacts relative to FBP.Source: §V · p. 15
- JTV gives modest synthetic accuracy gains and a cleaner experimental dark-field background than channel-wise TV, while lowering experimental mean dark-field posterior standard deviation by 23–30 % and raising limited-data median per-pixel ESS by 4–38 %.Source: §V · p. 15
- These posterior-spread and sampling-efficiency gains do not establish calibration; synthetic spread understated boundary error, and phase mixing remains slow under angular undersampling.Source: §V · p. 15