Mathematical description of the ACF using fft2
When the signal is assumed to be ergodic, then its ACF can be computed using time averages which can also be computed using the following convolution of length $N$ sequence $x[n]$ by its conjugate symmetric version:

$$ \hat{r}_{xx}[m] = \frac{1}{N} ~~~x[m] ~~\star ~~ x^*[-m] ~~ $$

Taking the $2N-1$ point DFT of both sides yields the following:

$$ \hat{R}_{xx}[k] = \frac{1}{N} ~~~X[k] X^*[k] =  \frac{1}{N} |X[k]|^2  $$

In a MATLAB implementation, an fftshift is required to bring $r_{xx}[0]$ to the center position for convenience. Hence the ACF will be obtained by an $2N-1$ point inverse DFT/FFT as 

rxk = (1/N)*real(fftshift( ifft( abs(fft(x,2N-1)).^2 ,2N-1) ) );

The 2D case follows very similarly to this 1D outline...