Hamiltonian dynamics for stochastic reconstruction in emission tomography

Abstract

The AMIAS/RISE framework formulates emission tomography as a probabilistic inverse problem in which reconstructed images are sampled from a distribution defined by the measurement model and counting statistics. In this work we present a stochastic reformulation of this approach based on gradient-driven optimization combined with Hamiltonian Monte Carlo (HMC) sampling directly in high-dimensional voxel space. This formulation enables practical ensemble generation for tomographic image reconstructions and provides direct access to image fluctuations within the sampled ensemble. Beyond point reconstruction, we introduce a spatially resolved operator-weighted diagnostic-the sampled data-visible variance-which quantifies how image fluctuations propagate through the imaging operator and thereby probes the local conditioning of the inverse problem under realistic acquisition physics. Using controlled software phantoms, experimental anthropomorphic phantom measurements, and a clinical DATSCAN SPECT acquisition, we demonstrate that while point-estimate accuracy under ideal conditions is comparable to deterministic reconstruction methods, the stochastic formulation provides additional physically interpretable insight. In particular, the ensemble analysis distinguishes uncertainty arising from intrinsic ill-posedness of the inverse problem from that associated with forward-model inadequacy. The clinical example is included to illustrate methodological applicability under realistic acquisition statistics rather than to assess diagnostic performance. The results establish the stochastic reconstruction framework as a practical ensemble-based approach for uncertainty quantification and forward-model validation in emission tomography.

Figures (10)

• Figure 1: Key steps of the emission tomography process: radiotracer administration and photon propagation, detection and projection data formation, and the final inverse reconstruction problem.
• Figure 2: Overview of the stochastic reconstruction workflow and the concept of data-visible variance. Starting from the measured detector sinogram (top left), a probability density over images is defined as $\Pi(\mathbf{x}) \propto \exp(-\chi^2(\mathbf{x})/2)$ (top middle). A fast gradient-based reconstruction (SGD) is first used to locate a high-probability region of this distribution, providing initialization for Hamiltonian Monte Carlo (HMC) sampling. HMC then generates an ensemble of reconstructed images; the top-right panel illustrates the resulting $\chi^2$ distribution, and the two example images (bottom right) correspond to representative samples from this ensemble. From the ensemble, voxel-wise statistics can be computed directly in image space. In particular, the data-visible variance is obtained by propagating ensemble fluctuations through the data-weighted projection--backprojection operator $\mathbf{H}=\mathbf{P}^{\top}\mathbf{W}\mathbf{P}$ (bottom row), thereby isolating the component of variability that remains coupled to the measured data. Comparison of these maps under different forward-model assumptions (e.g. without and with attenuation correction) illustrates how ensemble diagnostics can be used to assess forward-model adequacy beyond global residual measures.
• Figure 3: (a) True activity distribution (normalized). (b) MLEM reconstruction and corresponding local $\chi^2$ map. (c) SGD reconstruction and corresponding local $\chi^2$ map. Both methods converge to visually and statistically comparable solutions under idealized conditions. The $\chi^2$ map metric is analyzed in detail in metrics.
• Figure 4: Voxel-intensity histograms for the true image (blue), MLEM (orange), and SGD (green). This metric is analyzed in detail in metrics.
• Figure 5: Application of the HMC framework to a realistic 3D voxelized software phantom. Left column: mean reconstruction. Middle column: spatially resolved $\chi^2$ map. Right column: sampled data-visible variance $\sigma(H\delta x)$. Rows correspond to forward-model configurations: NC; Model 1 (uniform attenuation, ideal collimator); Model 2 (attenuation correction with realistic collimator); Model 3 (incomplete attenuation correction). Structured residuals are progressively suppressed from NC to Model 2, after which $\chi^2$ becomes largely noise-like. In contrast, $\sigma(H\delta x)$ decreases and homogenizes with model refinement, while Model 3 reintroduces spatially structured variance, demonstrating sensitivity to attenuation-model mismatch. All $\sigma(H\delta x)$ maps are normalized by the mean reconstructed activity within the object mask (dimensionless).