Generate a Gradient Operator for a Fourier Transform

Question

I came across this equation while trying to process an image $z$,

$$ \mathcal{F}\left(\mathbf{D}^T \mathbf{D}\right) \mathcal{F}\left(z\right), $$

where $\mathcal{F}$ is the 2-D Fourier transform, and $\mathbf{D}$ is the gradient operator, which I think is in matrix form. However, I am not sure on how to implement this on code; it is simple if $\mathbf{D}$ is acting on the image, but the Fourier transform only acts on the operators themselves.

I came across these two posts that generates the horizontal and vertical gradients and transposes the gradient, but I am not sure on how to apply the horizontal and vertical gradients to acquire $\mathbf{D}$ since they have different sizes on a non-square image (I have thought of cutting the image into square sections, but the $\mathcal{F}\left(z\right)$ accounts the whole image), while the example given in the transposition of the gradient is acting on the image itself (although I could act the transpose of the gradient to the gradient itself, but that leads me to the first problem).

My initial thought is to use $\mathbf{D}$ as,

$$ \mathbf{D} = \begin{bmatrix} D_h \\ D_v \end{bmatrix}, $$

where $D_h$ and $D_v$ are the horizontal and vertical gradients, respectively, but I am not sure what their sizes should be if I were to append zeros, given for example that my image is $n \times m$.

This equation was found on the paper Optics Temperature Dependent Non Uniformity Correction Via $ {l}_{0} $ Regularized Prior for Airborne Infrared Imaging Systems, at equation (7).

So, my question is how to generate the gradient operator $\mathbf{D}$. Thank you in advance.

Answer 1

There are 2 things to remember in the context you're asking:

1. Image as Column Stacked Vector
   It seems what's confuses you is that you think the operator is working on the image in its 2D form (Matrix). Yet in practice you should reshape it into a column vector. In my answer you linked to, have a look how it operates on the image in mIx = reshape(mDh * mI(:), numRows, numCols - 1);.

2. The Context of Your Paper
   In the paper above they use to apply the operators on the frequency domain to make the computation faster. So what you need to do is understand that $ {D}^{T} D $ is basically the Laplace Operator (See in Discrete Laplace Operator how it is implemented for 2D Discrete case). Then all needed is to use Replicate MATLAB's conv2() in Frequency Domain to build the operator in Fourier Transform and apply the equation.