title: Implementing Discrete Gaussian Kernel in Python?
body: <p>I'm looking to implement the <a href="https://en.wikipedia.org/wiki/Scale_space_implementation#The_discrete_Gaussian_kernel" rel="noreferrer">discrete Gaussian kernel</a> as defined by Lindeberg in his work about scale space theory.</p>
  
  <p>It is defined as T(n,t) = exp(-t)*I_n(t) where I_n is the <a href="https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions" rel="noreferrer">modified Bessel function of the first kind</a>.</p>
  
  <p>I am trying to implement this in Python using Numpy and Scipy but am running into some trouble.</p>
  
  <pre><code>def discrete_gaussian_kernel(t, n):
      return math.exp(-t) * scipy.special.iv(n, t)
  </code></pre>
  
  <p>I try plotting with:</p>
  
  <pre><code>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])
  </code></pre>
  
  <p>The output looks like:</p>
  
  <p><a href="https://i.stack.imgur.com/syKLI.png" rel="noreferrer"><img src="https://i.stack.imgur.com/syKLI.png" alt="python impl of discrete gaussian. it&#39;s not right! o.o"></a></p>
  
  <p>From the Wikipedia article, it should look like:</p>
  
  <p><a href="https://i.stack.imgur.com/iZFgs.png" rel="noreferrer"><img src="https://i.stack.imgur.com/iZFgs.png" alt="enter image description here"></a></p>
  
  <p>I assume I'm making some really trivial mistake. :/ Any thoughts?
  Thanks!</p>
  
  <p>EDIT:
  What I wrote is equivalent to <code>scipy.special.ive(n, t)</code>. I'm pretty sure it's supposed to be the modified Bessel function of the first kind, not the second, but can someone confirm?</p>