In recent years, many investigators have proposed Gibbs prior models to regularize images reconstructed from emission computed tomography data. Unfortunately, hyperparameters used to specify Gibbs priors can greatly influence the degree of regularity imposed by such priors and, as a result, numerous procedures have been proposed to estimate hyperparameter values, from observed image data. Many of these, procedures attempt to maximize the joint posterior distribution on the image scene. To implement these methods, approximations to the joint posterior densities are required, because the dependence of the Gibbs partition function on the hyperparameter values is unknown. Here, the authors use recent results in Markov chain Monte Carlo (MCMC) sampling to estimate the relative values of Gibbs partition functions and using these values, sample from joint posterior distributions on image scenes. This allows for a fully Bayesian procedure which does not fix the hyperparameters at some estimated or specified value, but enables uncertainty about these values to be propagated through to the estimated intensities. The authors utilize realizations from the posterior distribution for determining credible regions for the intensity of the emission source. The authors consider two different Markov random field (MRF) models-the power model and a line-site model. As applications they estimate the posterior distribution of source intensities from computer simulated data as well as data collected from a physical single photon emission computed tomography (SPECT) phantom.