For PIL.ImageFilter.GaussianBlur how what kernel is used and does the radius parameter relate to standard deviation?

Question

After reading an image with PIL I usually perform a Gaussian filter using scipy.ndimage as follow

import PIL
from scipy import ndimage

PIL_image = PIL.Image.open(filename)
data = PIL_image.getdata()
array = np.array(list(data)).reshape(data.size[::-1]+(-1,))
img = array.astype(float)
fimg = ndimage.gaussian_filter(img, sigma=sigma, mode='mirror', order=0)

There is Gaussian blur function within PIL as follows (from this answer), but I don't know how it works exactly or what kernel it uses:

from PIL import ImageFilter
fimgPIL = PIL_image.filter(ImageFilter.GaussianBlur(radius=r)

This documentation does not provide details.

Questions about PIL.ImageFilter.GaussianBlur:

1. What exactly is the radius parameter; is it equivalent to the standard deviation σ?
2. For a given radius, how far out does it calculate the kernel? 2σ? 3σ? 6σ?

---

This comment on an answer to Gaussian Blur - standard deviation, radius and kernel size says the following, but I have not found information for PIL yet.

OpenCV uses kernel radius of (sigma * 3) while scipy.ndimage.gaussian_filter uses kernel radius of int(4 * sigma + 0.5)

Answer 1

From the source code, it looks like PIL.ImageFilter.GaussianBlur uses PIL.ImageFilter.BoxBlur. But I wasn't able to figure out how the radius and sigma are related.

I wrote a script to check the difference between scipy.ndimage.gaussian_filter and PIL.ImageFilter.GaussianBlur.

import numpy as np
from scipy import misc
from scipy.ndimage import gaussian_filter
import PIL
from PIL import ImageFilter
import matplotlib.pyplot as plt


# Load test color image
img = misc.face()

# Scipy gaussian filter
sigma = 5
img_scipy = gaussian_filter(img, sigma=(sigma,sigma,0), mode='nearest')

# PIL gaussian filter
radius = 5
PIL_image = PIL.Image.fromarray(img)
img_PIL = PIL_image.filter(ImageFilter.GaussianBlur(radius=radius))
data = img_PIL.getdata()
img_PIL = np.array(data).reshape(data.size[::-1]+(-1,))
img_PIL = img_PIL.astype(np.uint8)

# Image difference
img_diff = np.abs(np.float_(img_scipy) - np.float_(img_PIL))
img_diff = np.uint8(img_diff)

# Stats
mean_diff = np.mean(img_diff)
median_diff = np.median(img_diff)
max_diff = np.max(img_diff)

# Plot results
plt.subplot(221)
plt.imshow(img_scipy)
plt.title('SciPy (sigma = {})'.format(sigma))
plt.axis('off')

plt.subplot(222)
plt.imshow(img_PIL)
plt.title('PIL (radius = {})'.format(radius))
plt.axis('off')

plt.subplot(223)
plt.imshow(img_diff)
plt.title('Image difference \n (Mean = {:.2f}, Median = {:.2f}, Max = {:.2f})'
          .format(mean_diff, median_diff, max_diff))
plt.colorbar()
plt.axis('off')

# Plot histogram
d = img_diff.flatten()
bins = list(range(int(max_diff)))

plt.subplot(224)
plt.title('Histogram of Image difference')

h = plt.hist(d, bins=bins)
for i in range(len(h[0])):
    plt.text(h[1][i], h[0][i], str(int(h[0][i])))

Output for sigma=5, radius=5:

Output for sigma=30, radius=30:

The outputs of scipy.ndimage.gaussian_filter and PIL.ImageFilter.GaussianBlur are very similar and the difference is negligible. More than 95% of difference values are <= 2.

PIL version: 7.2.0, SciPy version: 1.5.0

Answer 2

This is a supplementary answer to @Nimal's accepted answer.

Basically the radius parameter is like sigma. I won't dig too deep, but I think the Gaussian kernel is slightly different internally in order to preserve normalization after rounding back to integers, since the PIL method returns 0 to 255 integer levels.

The script below generates an image that is 1 on the left and 0 on the right, then does a sigma = 10 pixel blur with both methods, then plots the center horizontal lines through each, plus their differences. I do difference twice since log can only display positive differences.

The first panel is the difference between PIL and the SciPy float results, the second is the truncated integer SciPy results, and the third is rounded SciPy.

import numpy as np
import matplotlib.pyplot as plt
import PIL
from scipy.ndimage import gaussian_filter
from PIL import ImageFilter

import PIL

sigma = 10.0
filename = 'piximg.png'

# Save a PNG with a central pixel = 1
piximg = np.zeros((101, 101), dtype=float)
piximg[:, :50] = 1.0
plt.imsave(filename, piximg, cmap='gray')

# Read with PIL
PIL_image = PIL.Image.open(filename)

# Blur with PIL
img_PIL = PIL_image.filter(ImageFilter.GaussianBlur(radius=sigma)) 
data = img_PIL.getdata()
img_PIL = np.array(list(data)).reshape(data.size[::-1]+(-1,))
g1 = img_PIL[..., 1]

# Blur with SciPy
data = PIL_image.getdata()
array = np.array(list(data)).reshape(data.size[::-1]+(-1,))
img = array.astype(float)
fimg = gaussian_filter(img[...,:3], sigma=sigma, mode='mirror', order=0)
g2 = fimg[..., 1]
g2u = np.uint8(g2)
g2ur = np.uint8(g2+0.5)

if True:
    plt.figure()
    plt.subplot(3, 1, 1)
    plt.plot(g1[50])
    plt.plot(g2[50])
    plt.plot(g2[50] - g1[50])
    plt.plot(g1[50] - g2[50])
    plt.yscale('log')
    plt.ylim(0.1, None)
    plt.subplot(3, 1, 2)
    plt.plot(g1[50])
    plt.plot(g2u[50])
    plt.plot(g2u[50] - g1[50])
    plt.plot(g1[50] - g2u[50])
    plt.yscale('log')
    plt.ylim(0.1, None)
    plt.subplot(3, 1, 3)
    plt.plot(g1[50])
    plt.plot(g2ur[50])
    plt.plot(g2ur[50] - g1[50])
    plt.plot(g1[50] - g2ur[50])
    plt.yscale('log')
    plt.ylim(0.1, None)
    plt.show()