In the context of image processing, given a $k$-th order, homogeneous and linear differential operator with constant coefficients, we study a class of variational problems whose regularizing terms depend on the operator. Precisely, the regularizers are integrals of spatially inhomogeneous integrands with convex dependence on the differential operator applied to the image function. The setting is made rigorous by means of the theory of Radon measures and of suitable function spaces modeled on $BV$. We prove the lower semicontinuity of the functionals at stake and existence of minimizers for the corresponding variational problems. Then, we embed the latter into a bilevel scheme in order to automatically compute the space-dependent regularization parameters, thus allowing for good flexibility and preservation of details in the reconstructed image. We establish existence of optima for the scheme and we finally substantiate its feasibility by numerical examples in image denoising. The cases that we treat are Huber versions of the first and second order total variation with both the Huber and the regularization parameter being spatially dependent. Notably the spatially dependent version of second order total variation produces high quality reconstructions when compared to regularizations of similar type, and the introduction of the spatially dependent Huber parameter leads to a further enhancement of the image details.