Motivated by the fractal-like behavior of natural images, we propose a new smoothing technique that uses a regularization functional which is a fractional iterate of the Laplacian. This type of functional has previously been introduced by Duchon in the context of radial basis functions (RBFs) for the approximation of non-uniform data. Here, we introduce a new solution to Duchon's smoothing problem in multiple dimensions using non-separable fractional polyharmonic B-splines. The smoothing is performed in the Fourier domain by filtering, thereby making the algorithm fast enough for most multi-dimensional real-time applications.