OK. In theory, convolution in the time domain is the same as multiplication in the frequency domain. And filtering a signal with a linear filter is mathematically the same as convolving by its impulse response.
So in the time domain, if you have an image $x_{n, m}$, and a filter with an impulse response $h_{n, m}$, then after filtering $x$ with $h$, you get a signal $y = x * h$. (Note that especially in image processing, a filter's impulse response is often called a kernel. Until you're smarter about this than I usually am, just treat "impulse response" and "kernel" as synonyms).
If you let $X = \mathcal F \left \lbrace x \right \rbrace$ denote the Fourier transform, then if $y = x * h$, $Y = H X$. This turns out to be a happy result from 1-dimensional Fourier analysis, and an even happier result is that if you're using Cartesian coordinates, it applies to 2-dimensional (on 3- or 4-, or 42-dimensional) Fourier analysis, too.
That's the theory. Basically, it's what you've found in your research*.
So if you have $h$ then you can compute $y$ as $y = x * h$ (or $y = h * x$ -- convolution commutes). But you can also compute $y = \mathcal F^{-1} \left \lbrace \mathcal F \left \lbrace \mathcal h \right \rbrace F \left \lbrace \mathcal x \right \rbrace \right \rbrace $. This seems to be dreadfully counter-productive, but for two things.
First, you can pre-compute $H = F \left \lbrace \mathcal h \right \rbrace$. If it's settled enough, you can even compute it during design time and just have a pre-defined array in code.
Second, there's this thing called the fast Fourier transform (FFT). From a mathematical perspective, a fast Fourier transform would be just an overly complicated way of doing a plain old discrete Fourier transform except that for large images it uses far fewer operations (hence -- Fast Fourier Transform).
Convolution of an $N \times M$ image by a $N_h \times M_h$ kernel takes -- if I'm getting my math right -- $N M N_h M_h$ multiply-add operations. Basically, it's the area of the kernel times the area of the image. Taking the FFT and inverse FFT of that same image area takes $k N M \log(N M)$ multiply-add operations**. So for a large enough filter kernel, using the FFT can be faster -- particularly because you can use a technique called "overlap-add" to speed things up.
So, theoretically there is no difference between doing your filtering using convolution in the spatial-domain or by multiplying in the frequency domain.
Practically, doing it in the frequency domain is always more complicated, but it can be faster. Also practically, when you use the FFT there can be some slight loss of numerical precision -- basically, the speed comes at the cost of a bit of extra numerical scrud due to rounding errors and whatnot. This numerical effect is almost certainly of no concern to you.
Finally -- practically --
- it's rare that you have useful filter kernels that are large enough that doing the computation in the frequency domain will speed things up.
- Moreover, if all you're interested in is Gaussian filtering, there are some fast almost-Gaussian filtering technique that apply a small (3x3, I think) kernel to an image repeatedly. Because of the Central Limit Theorem, if you choose a small kernel correctly and filter by it repeatedly, the resulting filtering will tend to approximate filtering by a Gaussian.
* Things may have been a bit clouded by the fact that there's a result that's specific to Gaussian filters, in that if $h_{n, m}$ is Gaussian, then so is $H(\omega_1, \omega_2)$ -- it's just that in general it has a different scaling. This is not always the case, but if you're looking at some page that doesn't fully explain it, the notion of multiplying by a Gaussian in the frequency domain being the same as convolving by a Gaussian in the spatial domain would look extra-weird.
** I don't remember the exact number, but I think for a naïve FFT implementation, $k=6$, and there's more clever implementations that get it smaller yet.
