Computing a 2D Fourier transform on a disk

I am looking for a reliable and fast way of evaluating an integral like $$F(r, \phi)= \int_0^1 \int_0^{2\pi} f(\rho, \theta) e^{2\pi i \rho r \cos(\theta - \phi)}\rho\, d\theta\,d\rho,$$ where $f$ is a complex-valued function defined on the unit disk $\mathbb{D}=\{|z|\leq 1 \}$, for many values of $(r, \phi)$ and for several functions $f$. So, quadrature-based methods shouldn't work well here. Apparently neither the approximation by a discrete Fourier transform (evaluated by the FFT) works well. There is a so-called Extended Nijboer-Zernike (ENZ) Analysis (http://www.nijboerzernike.nl/) that tries to do it semi-analytically, approximating $f$ by Zernike polynomials, but still it is not totally satisfactory.

So, my question is: what is the state of the art of the numerical computation of such type of integrals?

-
Does your integral actually depend on $\phi$? –  Carl Nov 8 '13 at 9:04
@Carl I think so, $\phi$ is sitting in the cosine (argument of the exponential). –  Andrei MF Nov 8 '13 at 15:33
Sorry - you're right, though it might still simplify if you express $f$ as a Fourier series in $\theta$. Why are the Zernike polynomials unsatisfactory? –  Carl Nov 8 '13 at 16:45
@Carl The formulas that the ENZ theory gets using Zernike polynomials involve evaluation of infinite series of Bessel functions times certain binomial coefficients. Moreover, the terms of these sums are sign-changing and occasionally of similar size, thus cancellation errors might become a problem. –  Andrei MF Nov 8 '13 at 17:28