Auxiliary Variable Methods for Markov Chain Monte Carlo with Applications

Image quality of positron emission tomography (PET) reconstructions is degraded by subject motion occurring during the acquisition. Magnetic resonance (MR)-based motion correction approaches have been studied for PET/MR scanners and have been successful at capturing regular motion patterns, when used in conjunction with surrogate signals (e.g. navigators) to detect motion. However, handling irregular respiratory motion and bulk motion remains challenging. In this work, we propose an MR-based motion correction method relying on subspace-based real-time MR imaging to estimate motion fields used to correct PET reconstructions. We take advantage of the low-rank characteristics of dynamic MR images to reconstruct high-resolution MR images at high frame rates from highly undersampled k-space data. Reconstructed dynamic MR images are used to determine motion phases for PET reconstruction and estimate phase-to-phase nonrigid motion fields able to capture complex motion patterns such as irregular respiratory and bulk motion. MR-derived binning and motion fields are used for PET reconstruction to generate motion-corrected PET images. The proposed method was evaluated on in vivo data with irregular motion patterns. MR reconstructions accurately captured motion, outperforming state-of-the-art dynamic MR reconstruction techniques. Evaluation of PET reconstructions demonstrated the benefits of the proposed method in terms of motion artifacts reduction, improving the contrast-to-noise ratio by up to a factor 3 and achieveing a target-to-background ratio up to 90% superior compared to standard/uncorrected methods. The proposed method can improve the image quality of motion-corrected PET reconstructions in clinical applications.

Positron emission tomography (PET) reconstruction has become an ill-posed inverse problem due to low-count projection data, and a robust algorithm is urgently required to improve imaging quality. Recently, the deep image prior (DIP) has drawn much attention and has been successfully applied in several image restoration tasks, such as denoising and inpainting, since it does not need any labels (reference image). However, overfitting is a vital defect of this framework. Hence, many methods have been proposed to mitigate this problem, and DeepRED is a typical representation that combines DIP and regularization by denoising (RED). In this article, we leverage DeepRED from a Bayesian perspective to reconstruct PET images from a single corrupted sinogram without any supervised or auxiliary information. In contrast to the conventional denoisers customarily used in RED, a DnCNN-like denoiser, which can add an adaptive constraint to DIP and facilitate the computation of derivation, is employed. Moreover, to further enhance the regularization, Gaussian noise is injected into the gradient updates, deriving a Markov chain Monte Carlo (MCMC) sampler. Experimental studies on brain and whole-body datasets demonstrate that our proposed method can achieve better performance in terms of qualitative and quantitative results compared to several classic and state-of-the-art methods.

The Bayesian approach to inverse problems, in which the posterior probability distribution on an unknown field is sampled for the purposes of computing posterior expectations of quantities of interest, is starting to become computationally feasible for partial differential equation (PDE) inverse problems. Balancing the sources of error arising from finite-dimensional approximation of the unknown field, the PDE forward solution map and the sampling of the probability space under the posterior distribution are essential for the design of efficient computational Bayesian methods for PDE inverse problems. We study Bayesian inversion for a model elliptic PDE with an unknown diffusion coefficient. We provide complexity analyses of several Markov chain Monte Carlo (MCMC) methods for the efficient numerical evaluation of expectations under the Bayesian posterior distribution, given data δ. Particular attention is given to bounds on the overall work required to achieve a prescribed error level ε. Specifically, we first bound the computational complexity of ‘plain’ MCMC, based on combining MCMC sampling with linear complexity multi-level solvers for elliptic PDE. Our (new) work versus accuracy bounds show that the complexity of this approach can be quite prohibitive. Two strategies for reducing the computational complexity are then proposed and analyzed: first, a sparse, parametric and deterministic generalized polynomial chaos (gpc) ‘surrogate’ representation of the forward response map of the PDE over the entire parameter space, and, second, a novel multi-level Markov chain Monte Carlo strategy which utilizes sampling from a multi-level discretization of the posterior and the forward PDE. For both of these strategies, we derive asymptotic bounds on work versus accuracy, and hence asymptotic bounds on the computational complexity of the algorithms. In particular, we provide sufficient conditions on the regularity of the unknown coefficients of the PDE and on the approximation methods used, in order for the accelerations of MCMC resulting from these strategies to lead to complexity reductions over ‘plain’ MCMC algorithms for the Bayesian inversion of PDEs.

The authors studied certain characteristics of the expectation-maximization (EM) algorithm for image reconstruction in positron emission tomography (PET). They have attempted to speed up the algorithm by performing some time-consuming computation only once, at the initiation stage. They have studied the resulting memory size, disk space, and CPU requirements versus the number of pixels in the reconstructed image and the number of gamma -ray detectors in the tomograph. The convergence and image quality have been studied as a function of the number of iterations performed. Preliminary results are presented. >

In this article, we evaluate Parallel Level Sets (PLS) and Bowsher's method as segmentation-free anatomical priors for regularized brain positron emission tomography (PET) reconstruction. We derive the proximity operators for two PLS priors and use the EM-TV algorithm in combination with the first order primal-dual algorithm by Chambolle and Pock to solve the non-smooth optimization problem for PET reconstruction with PLS regularization. In addition, we compare the performance of two PLS versions against the symmetric and asymmetric Bowsher priors with quadratic and relative difference penalty function. For this aim, we first evaluate reconstructions of 30 noise realizations of simulated PET data derived from a real hybrid positron emission tomography/magnetic resonance imaging (PET/MR) acquisition in terms of regional bias and noise. Second, we evaluate reconstructions of a real brain PET/MR data set acquired on a GE Signa time-of-flight PET/MR in a similar way. The reconstructions of simulated and real 3D PET/MR data show that all priors were superior to post-smoothed maximum likelihood expectation maximization with ordered subsets (OSEM) in terms of bias-noise characteristics in different regions of interest where the PET uptake follows anatomical boundaries. Our implementation of the asymmetric Bowsher prior showed slightly superior performance compared with the two versions of PLS and the symmetric Bowsher prior. At very high regularization weights, all investigated anatomical priors suffer from the transfer of non-shared gradients.

Including accurate modeling of the point spread function (PSF) in positron emission tomography (PET) reconstruction algorithms results in improvements in image spatial resolution and contrast. In this work, we sampled the PSF in our first-generation radio-frequency brain dedicated PET insert for simultaneous PET/MR imaging using a 100 µCi NEMA standard 250 µm diameter Na-22 point source at 13 different positions within a subsection of the system field of view (FoV). The acquired list mode data was converted into the canonical sinogram format from which the spatial positioning of the source and standard deviations were calculated. The subset was then used to extrapolate the PSF for the full system FoV. This model was then fed as an input parameter into a graphical processing unit based ordered subset expectation maximization (OSEM) reconstruction algorithm and used to generate reconstructed images with and without spatially varying PSF modeling for the Na-22 point source and a Hoffman brain phantom. Results indicate that for point source reconstruction, the FWHM of the horizontal profile of the point source is smaller with spatially variant PSF especially closer to the edges. Effect of spatially varying PSF modeling is also presented with Hoffman phantom reconstruction and CNR value has increased with spatially varying PSF.

The promise of new PET image reconstruction

MRI-guided (MRIg) positron emission tomography (PET) reconstruction can potentially reduce noise and increase spatial resolution compared to standard clinical ordered-subsets expectation-maximization (OSEM) image quality. However, to adjust for the desired image quality, the balance between measured data and prior information usually requires manual tuning. This work presents an adaptive method to automatically control the influence of the magnetic resonance imaging (MRI) information on the PET emission data using maximum a posteriori (MAP) image reconstruction, robust against a wide range of counts. The method was evaluated on different static brain PET datasets using [18F]-FDG, [18F]-Florbetapir and [11C]-PiB, acquired in a simultaneous PET/MRI scanner and a PET/CT scanner, followed by an MRI scan. Noise in gray and white matter was measured for a wide range of statistics. Furthermore, noise and quantification accuracy were analyzed in different cortical and subcortical brain regions with different levels of tracer uptake, and at different levels of counts. Results demonstrated consistent improved image quality in terms of noise and spatial resolution with MRI-guided MAP PET (MRIg-MAP) reconstruction compared to OSEM. Additionally, it was shown that the number of collected counts could be reduced by ~1.6– <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">$2.3\times $ </tex-math></inline-formula> using MRIg-MAP reconstruction compared to OSEM, without increasing the noise significantly, either by reducing the scan time or injected activity. In conclusion, we presented a novel method to adaptively balance the influence of the anatomical information on the emission data for MRIg-MAP reconstruction, which showed image quality improvements compared to OSEM for different radiotracers, at different levels of counts, and applicable to simultaneous and sequential PET-MRI scans.

Where several different kinds of geophysical datasets have been acquired from a particular region, each of these can contain valuable information about the Earth, which may not be present in the other datasets. Jointly determining a common model, therefore, often gives a more thorough and more constrained description of the Earth structure than considering each dataset individually. For example, a seismic velocity inversion is only weakly constrained by first arrival seismic refraction data, but considering it alongside Magneto‐Telluric (MT) and gravity data can greatly assist in the constraint (Jegen‐Kulcsar et al., 2009). Strategies for joint inversion are therefore an active area of research. To date, most schemes for accomplishing this have been deterministic in nature. Using a deterministic technique often means that it is conceptually difficult to include prior beliefs about the system under determination, uncertainties both in measurement and the relationship between the different physical quantities (velocity, resistivity, density), and the discrepancy between the model and the real Earth. Statistical strategies such as MCMC (Markov Chain Monte Carlo) model searches exist for assessing this kind of problem, but the number of potentially computationally expensive forward model runs required to effectively sample the whole model space and thus achieve a meaningful result is normally prohibitively high (> 105), even for simple 1D models, so such schemes are not generally implemented. However, a technique known as emulation is used in various scientific fields eg. cosmology (Vernon and Goldstein, 2009), whereby computationally expensive forward modelling code (a simulator) is approximated by an uncertainty‐calibrated computationally cheap function. Here we apply emulation to the problem of stochastic joint model determination. We show that emulation can be used to quickly exclude large areas of implausible model space, allowing fast updating of beliefs about an Earth structure. It thus provides a means by which the input model space for a deterministic inversion or MCMC scheme can be greatly reduced. We also show how an emulator can, by itself, effectively constrain a region of the Earth. We demonstrate the concept using a 1D model.

We address the problem of parameter estimation for diffusion driven stochastic volatility models through Markov chain Monte Carlo (MCMC). To avoid degeneracy issues we introduce an innovative reparametrization defined through transformations that operate on the time scale of the diffusion. A novel MCMC scheme which overcomes the inherent difficulties of time change transformations is also presented. The algorithm is fast to implement and applies to models with stochastic volatility. The methodology is tested through simulation based experiments and illustrated on data consisting of US treasury bill rates.

This work is based on a cooperative co-evolution algorithm called ‘Fly Algorithm’, which is an evolutionary algorithm (EA) where individuals are called ‘flies’. It is a specific case of the ‘Parisian Approach’ where the solution of an optimisation problem is a set of individuals (e.g. the whole population) instead of a single individual (the best one) as in typical EAs. The optimisation problem considered here is tomography reconstruction in positron emission tomography (PET). It estimates the concentration of a radioactive substance (called a radiotracer) within the body. Tomography, in this context, is considered as a difficult ill-posed inverse problem. The Fly Algorithm aims at optimising the position of 3-D points that mimic the radiotracer. At the end of the optimisation process, the fly population is extracted as it corresponds to an estimate of the radioactive concentration. During the optimisation loop a lot of data is generated by the algorithm, such as image metrics, duration, and internal states. This data is recorded in a log file that can be post-processed and visualised. We propose using information visualisation and user interaction techniques to explore the algorithm’s internal data. Our aim is to better understand what happens during the evolutionary loop. Using an example, we demonstrate that it is possible to interactively discover when an early termination could be triggered. It is implemented in a new stopping criterion. It is tested on two other examples on which it leads to a 60% reduction of the number of iterations without any loss of accuracy.

Objective.This article explores a new graphics processing unit (GPU)-based techniques for efficient image reconstruction in organ-targeted positron emission tomography (PET) scanners with planar detectors.Approach.GPU-based reconstruction is applied to the Radialis low-dose organ-targeted PET technology, developed to overcome the issues of high exposure and limited spatial resolution inherent in traditional whole-body PET/CT (Computed Tomography) scans. The Radialis planar detector technology is based on four-side tileable sensor modules that can be seamlessly combined into a sensing area of the needed size, optimizing the axial field-of-view for specific organs, and maximizing geometric sensitivity. The article explores the transition from central processing unit-based maximum likelihood expectation maximization algorithms to a GPU-based counterpart, demonstrating a tenfold overall speedup in image reconstruction with a hundredfold improvement in iteration speed.Main results.Through standardized PET performance tests and clinical image analysis, this work demonstrates that GPU-based image reconstruction maintains diagnostic image quality while significantly reducing reconstruction times. The application of this technology, particularly in breast imaging using the Radialis low-dose positron emission mammography, significantly reduces exam times thus improving patient comfort and throughput in clinical settings.Significance.This study represents an important advancement in the clinical workflow of PET imaging, providing insights into optimizing reconstruction algorithms to effectively leverage the parallel processing capabilities of GPUs.

The preconditioned Crank--Nicolson (pCN) method is a Markov chain Monte Carlo (MCMC) scheme, specifically designed to perform Bayesian inferences in function spaces. Unlike many standard MCMC algorithms, the pCN method can preserve the sampling efficiency under the mesh refinement, a property referred to as being dimension independent. In this work we consider an adaptive strategy to further improve the efficiency of pCN. In particular we develop a hybrid adaptive MCMC method: the algorithm performs an adaptive Metropolis scheme in a chosen finite dimensional subspace and a standard pCN algorithm in the complement space of the chosen subspace. We show that the proposed algorithm satisfies certain important ergodicity conditions. Finally with numerical examples we demonstrate that the proposed method has competitive performance with existing adaptive algorithms.

An algorithm is described for three-dimensional whole-body (3DWB) image reconstruction in positron emission tomography. For whole-body applications, improvements to the popular fixed-axial-acceptance-angle technique are achieved by combining axially adjacent projection data available with a long-axis data set. Time-consuming reprojection of unmeasured oblique lines of response is reduced or eliminated by axial overlap of bed positions, while pixel variance and reconstructed axial resolution are made more uniform by the overlap. Improvements in noise and resolution uniformity are accompanied by gains in reconstruction efficiency, and may be optimized against increased acquisition time due to overlapping acquisition segments and reduced axial coverage. An 11-detector-ring overlap improves axial uniformity in coronal images of a long, uniform cylinder from 23% to 8% with uniform axial resolution. Associated with a 37% improvement in reconstruction time is a 34% reduction in axial coverage for four bed positions. A smaller degree of overlap is found to provide the best trade-off between image uniformity, total scan duration, and reconstruction time because of a proportionally greater reduction in reprojected lines of response. Using a sample optimization scheme, we find a three-ring overlap is best for a 60 cm axial field of view and a five-ring overlap for an 80 cm axial field of view.

In recent years, various approaches to overcome the problems associated with iterative statistical image reconstruction methods have been investigated, including issues such as slow convergence rate, costly computational time and nonuniform correction efficiency. A multiscale Bayesian image reconstruction algorithm based on a 1-D wavelet technique is proposed. In this approach, the authors first apply a wavelet transform to the sinogram measurements, which yields a multiresolution object representation. The image is reconstructed at the corresponding grid by the Bayesian image reconstruction algorithm for each scale level in the frequency space; the backprojector is also modified to match each corresponding resolution scale. After the low-frequency components of the image have been recovered sufficiently on the coarser scale level, the resulted image is used as the starting point for the finer level at a new iteration. An important feature in the authors' new approach is that the wavelet decomposition is carried out in the sinogram space, while multiple grids are also used in the reconstructed image space.