11

I'm wondering what would be the easiest way to generate a 1D gaussian kernel in python given the filter length. I think that the idea is to evaluate the normal distribution for the values of the vector [-filter-length,...,filter_length], is it correct?

So far, I've done this, but I don't know why it is not correct:

result = np.zeros( filter_length )

mid = filter_length/2
result=[(1/(sigma*np.sqrt(2*np.pi)))*(1/(numpy.exp((i**2)/(2*sigma**2)))) for i in range(-mid,mid+1)]  

return result

where sigma is the standard deviation, which is a parameter. filter-length is also a parameter.

It's incorrect because I get, for example, for length=3 and sigma=math.sqrt(1.0/2/math.log(2))

[0.23485931967491286, 0.46971863934982572, 0.23485931967491286]

And it should be:

[0.25, 0.5, 0.25]

So, is there any problem of rounding? I don't know what is going on...

Edit I think that I should truncate somehow

Problem Solved The problem was that I wasn't normalizing. I had to divide the vector by the sum of all its components.

5
  • what value is sigma for the results you give? Commented Feb 17, 2013 at 1:04
  • 1
    for sigma=1 i think the values should be [0.24,0.40,0.24]. i am unsure why you expect [0.25,0.5,0.25]. Commented Feb 17, 2013 at 1:12
  • no, it's for sigma=math.sqrt(1.0/2/math.log(2))
    – synack
    Commented Feb 17, 2013 at 8:05
  • I have added what I think is the correct sigma to my answer, at least with the Gaussian formula you use. Why did you choose sigma=math.sqrt(1.0/2/math.log(2))?
    – HugoRune
    Commented Feb 17, 2013 at 10:21
  • [0.25, 0.5, 0.25] is a triangular kernel, not a Gaussian. You cannot make a Gaussian kernel this short, it will never share the good properties of the Gaussian kernel. You need more samples. See here for details: crisluengo.net/archives/695 Commented Mar 19 at 13:33

2 Answers 2

5

I am not very firm with numpy syntax, but if you convolve a kernel with a dirac impulse, you get the same kernel as output.

So you could simply use the inbuild scipy.ndimage.filters.gaussian_filter1d function, and use this array as input: [ 0, 0, 0, ... 0, 1, 0, ...0, 0, 0]

The output should be a gaussian kernel, with a value of 1 at its peak. (replace 1 with the maximum you want in your desired kernel)

So in essence, you will get the Gaussian kernel that gaussian_filter1d function uses internally as the output. This should be the simplest and least error-prone way to generate a Gaussian kernel, and you can use the same approach to generate a 2d kernel, with the respective scipy 2d function. Of course if the goal is to do it from scratch, then this approach is only good as a reference

In regards to your equation:
to get [..., 0.5, ...] as the output with your formula, you need to solve
(1/(sigma*np.sqrt(2*np.pi)) = 0.5
so the correct sigma should be
sigma = math.sqrt(2*1/np.pi)

3
  • reading those docs, you'd want mode='constant' and also it's not clear to me whether the amplitude will be 1, or whether the filter has an integrated value of 1. so you might need to divide by max(...). Commented Feb 17, 2013 at 1:00
  • what do you mean that you get the same kernel?
    – synack
    Commented Feb 17, 2013 at 8:17
  • 1
    Yes, you get the same kernel as output that the gaussian_filter1d function uses internally. I am pretty sure that this is the simplest way to generate a 1D Gaussian kernel. Of course, if you want to generate the kernel from scratch as an exercise, you will need a different approach. But you could at least use this method as a reference to compare to your output
    – HugoRune
    Commented Feb 17, 2013 at 9:56
0

For those that like to have a fully working copy/paste code example:

import numpy as np

filter_length = 3
sigma=math.sqrt(1.0/2/math.log(2))
result = np.zeros( filter_length )

mid = int(filter_length/2)
result=[(1/(sigma*np.sqrt(2*np.pi)))*(1/(np.exp((i**2)/(2*sigma**2)))) for i in range(-mid,mid+1)]  
sumresult = np.sum(result)
print(result/sumresult)

[0.25 0.5 0.25]

1
  • If you’re going to normalize, you can leave the 1/(sigma*np.sqrt(2*np.pi)) computation out. It’s no longer doing anything. Commented Mar 19 at 13:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.