7

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? 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?

kotakotakota
  • 731
  • 8
  • 26
  • Are you mixing up the order or arguments with the n and t? Looks like it to me. Or maybe not. I need coffee. – cadolphs Nov 17 '16 at 19:49
  • @Lagerbaer I don't think so. I checked it over a couple of times, and have also tried it but to no avail... – kotakotakota Nov 17 '16 at 19:54

1 Answers1

5

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.

pv.
  • 33,875
  • 8
  • 55
  • 49
  • That would do it. I guess that would explain why it was so explicit about integer values and why the discrete Gaussian kernel can't be used for continuous cases. Thanks for the help! – kotakotakota Nov 17 '16 at 19:58