Simultaneous hyperparameter estimation and Bayesian image reconstruction for PET

The authors present a new iterative algorithm for Bayesian PET image reconstruction that simultaneously estimates the PET image and the global hyperparameter β of a Gibbs prior. True maximum likelihood (ML) estimation of β is intractable for the PET reconstruction problem due to the complexity and high dimensionality of the probability densities involved. The new algorithm replaces the true likelihood function for the hyperparameter with an approximation in which the marginalization with respect to the image sample space is reduced to the product of a set of one dimensional integrals; one per image pixel. The approximation is closely related to the mean field theory of statistical mechanics. In essence, this reduction in complexity is achieved by approximating the influence of the neighbors of each pixel over their entire sample space with their estimated posterior modes. A preconditioned conjugate gradient algorithm is used to iteratively compute a MAP estimate of the image. At periodic intervals, the most recent image generated by this iterative procedure is used as an estimate of the posterior mode in the approximate marginalized log likelihood for the data given β, which in turn is used to update the ML estimate of β. The procedure is repeated until convergence of both the MAP image estimate and the ML estimate of β. Results of a validation study using Monte Carlo simulations are presented