Why isn't bilateral filtering shift-invariant? I have a source (that I can't share) that states that bilateral filtering isn't shift-invariant (that is, it depends on the input signal), hence we cannot use convolution and the Fourier convolution theorem (that is, convolution is the spatial domain corresponds to multiplication in the spectral one). 

Why isn't bilateral filtering shift-invariant? Why exactly can't we apply convolution and the convolution theorem to the bilateral filtering of an image? Why wouldn't a Gaussian filter also be signal dependent?

I understood that bilateral filtering combines domain and range filtering, that is, when updating a pixel of an image, it weights the contribution of other pixels in the image using both their Euclidean distance (to the pixel being updated) and their intensity (in case of grayscale) value. However, I don't get why we couldn't also use the convolution theorem to filter an image using BF.
The bilateral filter can be written as:

$$y[i, j] = \frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k, j+l]\,f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k,j+l] - x[i,j]\Big)\,g\Big(k^2+l^2\Big)}\tag{1}$$

Shift invariance can be defined as any shift $m,n$ of the input $x$ resulting in the corresponding shift in the output $y$:

$$x[i,j] \to x[i-m,j-n]\quad\Rightarrow\quad y[i,j]\to y[i-m,j-n].\tag{2}$$

By doing the first substitution on all occcurrences of $x$ on the right side of Eq. 1 we get as the output:

$$\begin{gather}\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i+k-m, j+l-n]\,f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i+k-m,j+l-n] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}\\
=\frac{\displaystyle\sum_{k,l\in \mathbb Z^2}x[i-m+k, j-n+l]\,f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)}{\displaystyle\sum_{k,l\in \mathbb Z^2}f\Big(x[i-m+k,j-n+l] - x[i-m,j-n]\Big)\,g\Big(k^2+l^2\Big)},\end{gather}\tag{3}$$

where we can recognize $i-m$ and $j-n$ as $i$ and $j$ of Eq. 1, so we have:

$$= y[i-m, j-n],\tag{4}$$

showing that the bilateral filter is shift-invariant.

However, these are shifts with integer $m$ and $n$. The non-linearity of the bilateral filter will make it not invariant to fractional shifts carried out by (sinc) interpolation. Perhaps this is what your source was referring to.