Synthetic magnetic resonance (MR) imaging is an approach suggested in the literature to predict MR images at different design parameter settings from at least three observed MR scans. However, performance is poor when no regularization is used in the estimation and otherwise computationally impractical to implement for three-dimensional imaging methods. We propose a method which accounts for spatial context in MR images by the imposition of a Gaussian Markov Random Field (MRF) structure on a transformation of the spin-lattice relaxation time, the spin-spin relaxation time and the proton density at each voxel. The MRF structure is specified through a Matrix Normal distribution. We also model the observed magnitude images using the more accurate but computationally challenging Rice distribution. A One-Step-Late Expectation-Maximization approach is adopted to make our approach computationally practical. We evaluate predictive performance in generating synthetic MR images in a clinical setting: our results indicate that our suggested approach is not only computationally feasible to implement but also shows excellent performance.