I came across multiple papers (e.g. here and here) where an inversion of the 2D Laplace operator $\Delta = \partial_x^2+\partial_y^2$ was obtained via Fourier transform. This was applied to obtain images for situations where only the gradient of the image was known. I understand that we can write the Laplacian of any function in terms of the Fourier transform like
$$\Delta f(\vec{x}) = -\mathcal{F}^{-1}\left[ ||\vec{q}||^2 \mathcal{F}\left[ f(\vec{x})\right](\vec{q}) \right](\vec{x})$$
if we define the 2D Fourier transform as $\mathcal{F}=\int \int f(\vec{x}) \exp(-i\vec{q}\cdot\vec{x})\text{d}x \text{d}y $. I also understand that we can invert the equation above and that the inversion might not lead to the integral we might expect, because any function $g(\vec{x})$ with $\Delta g(\vec{x})=0$ could be added to $f$ to obtain the same Laplacian.
My question is this: when I know the gradients $\partial_x f$ and $\partial_yf$ why would I differentiate them again and then use the inversion of the Laplace operator to integrate them as they did in this paper? Why could I not define an operator $G = \partial_x + \partial_y$ and express it in terms of the Fourier transform like so :
$$G f(\vec{x}) = \mathcal{F}^{-1}\left[ i(q_x+q_y)\mathcal{F}\left[f(\vec{x})\right](\vec{q})\right](\vec{x})$$
I could then also express $G^{-1}$ in terms of the Fourier Transform and apply it to the sum of the gradients, couldn't I? Sure, the "constant" of integration would be any function $h$ with $\partial_x h = -\partial_y h$ but is that any worse than what I get with the Laplace inversion?
I am glad for all help, because surely there is something that I am not seeing. Has it maybe to do with the fact that we are using the DFT on discrete images? Does it have something to do with the boundary conditions we are implicitly imposing?