Understanding the Bilateral Filter (Image Filtering)

I asked this question already on StackOverflow and was told to ask it here instead, so I will just copy the content of my question. I do understand the general concept of bilateral filtering and I did read several sources on the topic (including some questions on StacOverflow). Most sources use this equation to describe the process (I think originally it is from Course Notes - A Gentle Introduction to Bilateral Filtering and its Applications source, but at this point I am not sure): With the "general" gaussian filter (that only considers the space) being represented by this The Gaussian Filter I did already implement successfully. It is only for gray-scale images and so will the bilateral filter (hopefully) be. It's in no way having a good performance, I am aware, but it does the job. In general I'd think I just need to add the multiplication with (Ip - Iq) to make it work as a bilateral filter. My issue now is: what exactly does (Ip - Iq) represents? Or asked differently, where do I retrieve those values?
I'm not sure if my code is needed but here's my implementation anyways. The kernel-algorithm I took from SO - How to calculate a Gaussian kernel matrix efficiently in numpy.

def gaussfiltering (img = [], kernelsize=3, sigma=1.9):
rows = len(img)
columns = len(img)
kernel = gaussian_kernel(kernelsize, sigma)
horizontallyfiltered = iterate_horizontally(kernel, img, rows, columns, len(kernel))
filtered = iterate_vertically(kernel, horizontallyfiltered, rows, columns, len(kernel))
return filtered

def gaussian_kernel(kernel_size=3, sigma=1.9):
x = numpy.linspace(-sigma, sigma, kernel_size+1)
y = st.norm.cdf(x)
kernel = numpy.diff(y)
return((kernel/kernel.sum()))

def iterate_horizontally (kernel=[], img=[], rows=0, columns=0, kernelsize=0):
horizontallyfiltered = img.copy()
for r in range(rows):
for c in range(columns):
halved = int(kernelsize/2)
if ((c - halved) >= 0 and (c + halved) < columns): #is in range for filter
indexOfImage = c - halved
summed = 0.0
for i in range(kernelsize):
summed += img[r][indexOfImage] * kernel[i]
indexOfImage+=1
horizontallyfiltered[r][c]=summed
return horizontallyfiltered

def iterate_vertically (kernel=[], img=[], rows=0, columns=0, kernelsize=0):
verticallyfiltered = img.copy()
for r in range(rows):
for c in range(columns):
halved = int(kernelsize/2)
if ((r - halved) >= 0 and (r + halved) < rows): #is in range for filter
indexOfImage = r - halved
summed = 0.0
for i in range(kernelsize):
summed += img[indexOfImage][c] * kernel[i]
indexOfImage+=1
verticallyfiltered[r][c]=summed
return verticallyfiltered

There are some good resources on our site:

In my opinion, the best way to understand the Bilateral Filter is looking at concise and simple implementation.

So in What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF I showed how to calculate the spatial and range weights.

function [ mW ] = CalcSpatialWeights( kernelRadius, spatialStd )

vW = exp(-(vX .* vX) / (2 * spatialStd * spatialStd));
mW = vW.' * vW;
mW = mW / sum(mW(:));

end

We can see that for any pixel coordination the weights are the same. This is since it is basically Gaussian Blur which is implemented in a convolution. If something is implemented by convolution it means it is shift / spatially invariant. Namely, the weights are independent of the pixels values, but only depend on the relative distance.

Now, the range weights:

function [ mW ] = CalcRangeWeights( mI, vRefPixlCoord, kernelRadius, rangeStd )

kernlLength = (2 * kernelRadius) + 1;

refPixelVal = mI(vRefPixlCoord(1), vRefPixlCoord(2));

mW = zeros(kernlLength, kernlLength);

for ii = 1:(kernlLength * kernlLength)
mW(ii) = exp(-((mP(ii) - refPixelVal) ^ 2) / (2 * rangeStd * rangeStd));
end

mW = mW / sum(mW(:));

end

Now we can see the values will depend on the pixel location (As it defines its neighborhood).
This is basically a Gaussian Kernel yet the difference is the difference of values of the pixels relative to the value of the reference pixel.

The Gaussian Bilateral Filters multiplies both weights so the complete effect is:

1. The weighing of a pixel decreases with it distance from the reference pixel.
2. The weighing of a pixel decreases with its value being different from the reference value.

Hence the Bilateral Filter forms a local filter with edge preserving filtration.

A simple filter (convolution) replace value of a pixel (for simplicity assume its the pixel positioned at the center of kernel) with a weighted sum of its neighbors and the weights are constant ( the spatial kernel only depend on the relative position of center pixel and the neighbor pixel). In a bilateral filter the weights (in addition to it's relative position) also depend on the intensity value of center pixel and the neighbor which is represented by range kernel.

In your code, you have to create a 2D discrete gaussain (either by writing the gaussian using exponential function or multiplying 2 gaussian P.D.F (why differentiating a C.D.F?) for the seperable kernel case). For the Range kernel also either give the difference of intensity values to a Gaussian function using pdf object or write the gaussian function yourself using exponentials.