Summary: Scattered data approximation refers to the computation of a multi-dimensional function from measurements obtained from scattered spatial locations. For this problem, the class of methods that adopt a roughness minimization are the best performing ones. These methods are called variational methods and they are capable of handling contrasting levels of sample density. These methods express the required solution as a continuous model containing a weighted sum of thin-plate spline or radial basis functions with centres aligned to the measurement locations, and the weights are specified by a linear system of equations. The main hurdle in this type of method is that the linear system is ill-conditioned. Further, getting the weights that are parameters of the continuous model representing the solution is only a part of the effort. Getting a regular grid image requires re-sampling of the continuous model, which is typically expensive. We develop a computationally efficient and numerically stable method based on roughness minimization. The method leads to an algorithm that uses standard regular grid array operations only, which makes it attractive for parallelization. We demonstrate experimentally that we get these computational advantages only with a little compromise in performance when compared with thin-plate spline methods.