kotakotakota - 1 year ago 112

Python Question

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)`

Answer Source

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.