CODEWITHSUNDEEP
pythonnumpyscipy
zwep
zwep

Reputation: 1340

convolution with gaussian vs gaussian filter

Simple task.. I want to smoothen out some vector with a Gaussian.. This is just a test case, later on I want to apply this to an image.

import numpy as np
import scipy.stats
import scipy.ndimage

m = 7  # size of the 'signal'
n = 7  # size of the filter
sgm = 2  # dev for standard distr
weight_conv = np.zeros(2*m*n).reshape(2*n, m)  # Weights for the convolution

input_signal = np.array(range(m))  # input signal..
x1 = np.linspace(-4*sgm, 4*sgm, n)  # x-values for the normal-dstr
input_filter = scipy.stats.norm.pdf(x1, loc=0, scale=sgm)

# create my own weight matrix
for i in range(weight_conv.shape[1]):
    weight_conv[i:(len(input_filter)+i), i] = input_filter

# My own way of calculating the convolution
np.sum(weight_conv * input_signal, axis=1)
# Convolution provided by numpy
np.convolve(input_signal, input_filter)
# Apply the scipy gaussian filter...
scipy.ndimage.filters.gaussian_filter(input_signal, sigma=sgm)
scipy.ndimage.filters.gaussian_filter1d(input_signal, sigma=sgm)

Now my idea is that these all should be similar. My method is does produce similar output as the numpy convolution, but the scipy method is different...

scipy.ndimage.filters.gaussian_filter(input_signal, sigma=sgm)
array([1, 1, 2, 3, 3, 4, 4])

Now it must be the case that scipy is doing something different. But WHAT? I dont know. I have checked the source, and there it seems they are just using a convolution with a gaussian kernel (which is what I am also doing). But the answers dont add up...

Anyone has a different idea?

Upvotes: 1

Views: 4594

Answers (2)

zwep
zwep

Reputation: 1340

With the help of @filippo and this SO-question I was able to reconstruct the scipy implementation. The way the method propagates the information to either side is crucial.

Below is a code snippet that shows how they are equal

import numpy as np
import scipy.stats
import scipy.ndimage
import matplotlib.pyplot as plt
np.set_printoptions(linewidth=160)


m_init = 7
# Create any signal here...
input_signal_init = []  
input_signal_init = np.arange(m_init)
input_signal_init = np.random.choice(range(m_init),m_init)
# Convert to float for better results in scipy.
input_signal_init = np.array(input_signal_init).astype(float)

# Simulating method='reflect'
input_signal = np.array([*input_signal_init[::-1], *input_signal_init, *input_signal_init[::-1]])

# Define new length of input signal
m = len(input_signal)  
# Properties of the Gaussian
sgm = 2  # dev for standard distr
radius = 4 * sgm
x = numpy.arange(-radius, radius+1)
n = len(x)
weight_conv = np.zeros(m*(n+m)).reshape(n+m, m)  

# Calculate the gaussian
p = np.polynomial.Polynomial([0, 0, -0.5 / (sgm * sgm)])
input_filter = numpy.exp(p(x), dtype=numpy.double)
input_filter /= input_filter.sum()

# Calculate the filter weights
for i in range(weight_conv.shape[1]):
    weight_conv[i:(len(input_filter)+i), i] = input_filter

# My own way of calculating the convolution
self_conv = np.sum(weight_conv * input_signal, axis=1)[(2*m_init+1):(3*m_init+1)]
# Convolution provided by numpy
numpy_conv = np.convolve(input_signal, input_filter)[(2*m_init+1):(3*m_init+1)]
# Convolution by scipy with method='reflect'
# !! Here we use t[![enter image description here][2]][2]he original 'input_signal_init'
scipy_conv = scipy.ndimage.filters.gaussian_filter(input_signal_init, sigma=sgm)

Plotting the results always give confidence that you did a good job... so

plt.plot(scipy_conv, 'r-')
plt.plot(self_conv, 'bo')
plt.plot(numpy_conv, 'k.-')
plt.show()

Gives the following image.

It can also be verified that setting the mode of the scipy filter to 'constant' will create identical convolutions as well.

Upvotes: 2

filippo
filippo

Reputation: 5294

Scipy multidimensional gaussian filter uses a bigger kernel. By default the kernel radius is truncated to 4 sigmas, which in your case should be somewhat similar to a 17x17 filter.

See _gaussian_kernel1d for the exact implementation. It also uses several 1d separable correlations but that shouldn't make much difference.

The other key difference is in the output vector size and precision. From ndimage docs:

The intermediate arrays are stored in the same data type as the output. Therefore, for output types with a lower precision, the results may be imprecise because intermediate results may be stored with insufficient precision. This can be prevented by specifying a more precise output type.

So in your case the ouput precision is limited by your input_signal.dtype. Try using a float input array or a different array for the ouput.

The output size and edge handling is a bit trickier, not sure if is there a way to get the same behavior from np.convolve and scipy.ndimage.gaussian_filter

Upvotes: 1

Related Questions