Unexpected Gaussian filtering behavior with Scipy
Considering that I have an image loaded f(x,y)
, for example
I want to compute the Gaussian derivative of an ∂/∂x ∂/∂y G*f
image f
, where G
is the Gaussian filter and *
denotes convolution. This can be easily done with Scipy:
from scipy.ndimage.filters import gaussian_filter
imshow(gaussian_filter(g, sigma, order=1))
With sigma=50
this yields the following result:
Now, for application reasons, I need to perform the calculation using mode='constant'
:
imshow(gaussian_filter(g, sigma, order=1, mode='constant', cval=0))
However, the result looks reasonable:
However, note that the intensity of the background image is the image 1
, not 0
. Hence, it is wise to use cval=1
:
imshow(gaussian_filter(g, sigma, order=1, mode='constant', cval=1))
Now this is unexpected! This result doesn't make sense, does it?
For the record, I also checked the partial differentials ∂/∂x G*f
and ∂/∂y G*f
. While
imshow(gaussian_filter(g, sigma, order=[0, 1], mode='constant', cval=1)
looks reasonable
other
imshow(gaussian_filter(g, sigma, order=[1, 0], mode='constant', cval=1)
:
Why is this?
source to share
There gaussian_filter
is an error in B that appears when both order
and cval
are nonzero. Specifically, it's here :
for axis, sigma, order, mode in axes:
gaussian_filter1d(input, sigma, axis, order, output, mode, cval, truncate)
input = output
The filter re-convolves 1d and every time it passes in cval
to 1d of the filter. The problem is that if any derivatives were taken, then cval
it should be set to 0, since the derivative of any constant is zero. This is why the result does not match order=[1, 0]
, but not with order=[0, 1]
. Without testing (no SciPy dev environment), I think it would be correct:
for axis, sigma, order, mode in axes:
gaussian_filter1d(input, sigma, axis, order, output, mode, cval, truncate)
if order > 0:
cval = 0.0
input = output
Bypass
Without non-zero, cval
it can be emulated by subtracting it from the image before filtering (and adding back after filtering only if the order is zero). Example:
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import gaussian_filter
g = np.ones((500, 500))
g[200:300, 200:300] = 2
sigma = 50
cval = 1
gf = gaussian_filter(g-cval, sigma, order=1, mode='constant')
plt.matshow(gf)
plt.show()
returns
what is the expected result. (My original image is slightly different from yours and I'm using a different render tool.)
source to share