Typical methods for image segmentation, or labeling, formulate and solve an optimization problem to produce a single optimal solution. For applications in clinical decision support relying on automated medical image segmentation, it is also desirable for methods to inform about (i) the uncertainty in label assign- ments or object boundaries or (ii) alternate close-to-optimal solutions. However, typical methods fail to do so. To estimate uncertainty, while some Bayesian methods rely on simplified prior models and approximate variational inference schemes, others rely on sampling segmentations from the associated posterior model using (i) traditional Markov chain Monte Carlo (MCMC) methods based on Gibbs sampling or (ii) approximate perturbation models. However, in such typical approaches, in practice, the resulting inference or generated sample set are approximations that deviate significantly from those indicated by the true posterior. To estimate uncertainty, we propose the modern paradigm of exact MCMC sampling to sample multi-label segmentations from generic Bayesian Markov random field (MRF) models, in finite time for exact inference. Fur- thermore, for exact sampling in generic Bayesian MRFs, we extend the the- ory underlying Fill’s algorithm to generic MRF models by proposing a novel bounding-chain algorithm. On several classic problems in medical image analy- sis, and several modeling and inference schemes, results on simulated data and clinical brain magnetic resonance images show that our uncertainty estimates gain accuracy over several state-of-the-art inference methods.