Convolution of an Image with a Kernel That Is a Product of Two Functions Suppose that $G(i,j)$ is a Gaussian decay function on the distance between points $i$ and $j$ of an image. In addition, $D(i,j)$ is the difference between the VALUES of the image at those points. 
Now, at every point $i$ of an image I need to compute the summation of the neighboring $G(i,j) D(i,j) $ , this is, the product of both functions. This is a convolution-like operation: for every point I add the contributions of a function at every neighboring point.
However, I do not know if that can be expressed as a convolution because G and D are functions on different variables: the distance (a raw gaussian filter) and the image intensity difference. 
I need help at expressing the summation as a combination of convolution operations in order to easily code it in a program without the need for loops. In addition, it would allow me to possibly do it in the frequency domain.
Any help will be appreciated!
It is equivalent to a bilateral filter. You can move the kernel to each local part of image to do the dot product, then sum them up to get the result in a loop. There is also a fast implementation with O(1) in matlab programming for your reference. And the relevant publications are:


K.N. Chaudhury, D. Sage, and M. Unser, "Fast O(1) bilateral  filtering using 
trigonometric range kernels," IEEE Transactions on Image Processing, vol. 20, 
no. 11, 2011.
K.N. Chaudhury, "Acceleration of the shiftable O(1) algorithm for bilateral filtering 
and non-local means," arXiv:1203.5128v1. 


I assume the Gaussian kernel is position shift-invariant. So normalization means the sum of elements inside the kernel equals to 1. For a 3*3 Gaussian kernel, it is

1/16  2/16  1/16
2/16  4/16  2/16
1/16  2/16  1/16


The bilateral filter can also be normalized with the K factor of integral of the product of  Gaussian kernel and intensity kernel. You can try both to observe the effect of normalization on your image data.