kotakotakota - 1 year ago 210
Python Question

# Implementing Discrete Gaussian Kernel in Python?

I'm looking to implement the discrete Gaussian kernel as defined by Lindeberg in his work about scale space theory.

It is defined as T(n,t) = exp(-t)*I_n(t) where I_n is the modified Bessel function of the first kind.

I am trying to implement this in Python using Numpy and Scipy but am running into some trouble.

``````def discrete_gaussian_kernel(t, n):
return math.exp(-t) * scipy.special.iv(n, t)
``````

I try plotting with:

``````import math
import numpy as np
import scipy
from matplotlib import pyplot as plt

def kernel(t, n):
return math.exp(-t) * scipy.special.iv(n, t)
ns = np.linspace(-5, 5, 1000)
y0 = discrete_gaussian_kernel(0.5, ns)
y1 = discrete_gaussian_kernel(1, ns)
y2 = discrete_gaussian_kernel(2, ns)
y3 = discrete_gaussian_kernel(4, ns)
plt.plot(ns, y0, ns, y1, ns, y2, ns, y3)
plt.xlim([-4, 4])
plt.ylim([0, 0.7])
``````

The output looks like:

From the Wikipedia article, it should look like:

I assume I'm making some really trivial mistake. :/ Any thoughts?
Thanks!

EDIT:
What I wrote is equivalent to
`scipy.special.ive(n, t)`
. I'm pretty sure it's supposed to be the modified Bessel function of the first kind, not the second, but can someone confirm?

If you want to get the wikipedia plot, replace

``````ns = np.linspace(-5, 5, 1000)
``````

with

``````ns = np.arange(-5, 5+1)
``````

Namely, the formula you use only makes sense at the integer points. The Bessel function as a function of negative order is an oscillating function, http://dlmf.nist.gov/10.27#E2 so the plot looks fine to me.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download