The problem of image reconstruction from incomplete data can be formulated as a linear inverse problem and is usually approached using optimization theory tools. Total variation (TV) regularization has been widely applied in this framework, due to its effectiveness in capturing spatial information and availability of elegant, fast algorithms. In this paper we show that significant improvements can be gained by extending this approach with a Markov Random Field (MRF) model for image gradient magnitudes. We propose a novel method that builds upon the Chambolle's fast projected algorithm designed for solving TV minimization problem. In the Chambolle's algorithm, we incorporate a MRF model which selects only a subset of image gradients to be effectively included in the algorithm iterations. The proposed algorithm is especially effective when a large portion of image data is missing. We also apply the proposed method to demosacking where algorithm shows less sensitivity to the initial choice of the tuning parameter and also for its wide range of values outperformes the method without the MRF model.