kotakotakota kotakotakota - 10 months ago 83
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:

python impl of discrete gaussian. it's not right! o.o

From the Wikipedia article, it should look like:

enter image description here

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

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?

pv. pv.
Answer Source

If you want to get the wikipedia plot, replace

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


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.