Fitting regression models can be challenging when regression coefficients are high-dimensional. Especially when large spatial or temporal effects need to be taken into account the limits of computational capacities of normal working stations are reached quickly. The analysis of images with several million pixels, where each pixel value can be seen as an observation on a new spatial location, represent such a situation. A Markov chain Monte Carlo (MCMC) framework for the applied statistician is presented that allows to fit models with millions of parameters with only low to moderate computational requirements. The method combines a modified sampling scheme with novel accomplishments in iterative methods for sparse linear systems. This way a solution is given that eliminates potential computational burdens such as calculating the log-determinant of massive precision matrices and sampling from high-dimensional Gaussian distributions. In an extensive simulation study with models of moderate size it is shown that this approach gives results that are in perfect agreement with state-of-the-art methods for fitting structured additive regression models. Furthermore, the method is applied to two real world examples from the field of medical imaging.