Summary: The problem of uncertainty quantification (UQ) for inverse problems has become of significant recent interest. However, UQ requires more than the classical methods for computing solutions of inverse problems. In this paper, we take a Bayesian approach for the solution of ill-posed deconvolution problems with a symmetric convolution kernel and Neumann boundary conditions. The prior is modeled as a Gaussian Markov random field (GMRF) with the same boundary conditions and symmetry assumptions. These assumptions yield better results in certain instances and also allow for the use of the discrete cosine transform for fast computations. Moreover, we use a hierarchical model for the noise precision (inverse-variance) and prior precision parameters. This leads to an a posterior density function from which we can compute samples using a basic Markov chain Monte Carlo (MCMC) method. The resulting samples can then be used for both estimation (using, e.g., the sample mean) and uncertainty quantification (using, e.g., histograms, the sample variance, or a movie created from the image samples). We provide a numerical experiment showing that the method is effective, computationally efficient, and that for certain problems, the boundary conditions can yield significantly better results than if a periodic boundary is assumed. The novelty in the work lies in the combination of the MCMC method, Neumann boundary conditions, GMRF priors, and in the use of a movie to visualize uncertainty in the unknown image.