How to choose basis functions that contribute most efficiently per term to an approximation of an image f(x,y)?
I do not  understand why you need an infinite number of  terms for a Fourier approxiimation. Yet, if you are looking for tools on non-discretized 2D data,  I would suggest you to start with the concept of Matching pursuit of images. You take a dictionary of 2D atoms that is complete (could be complex exponentials, gaussian bumps). Then, you find the  atom that best matches you  data, subtract it with the necessary scale factor, and repeat.

You could also resort to anisotropic  triangulations, but  this would not suit your "need to be combined in an inner product fashion".

The  quality of approximation could be evaluated over some classes for regular functions ($C^p$ smooth inside $C^q$ curves), and are related to properties of  your  dictionary. You might find some references in section 4.1.1. Matching pursuits from A panorama on multiscale geometric representations, Signal Processing, 2011, L. Jacques et al.

I do suspect that the inner-product  constraint will limit the efficiency: optimal approximations in 1D do not extend generically  in 2D.