Are scipy second-order Gaussian derivatives correct?

Question

For an edge detection algorithm, I need to compute second-order derivatives of an image, and I do this with use of Gaussian derivatives. I assumed that the scipy.ndimage.gaussian_filter implementation of this principle would work well, but I get a confusing result. My question probably revolves around whether this is erroneous behaviour (i.e., a bug) on the side of scipy, or whether there is a fundamental reason for it.

The test array is

import numpy as np
import matplotlib.pyplot as plt
A = np.zeros((30,30))
A[:,15] = 1
plt.imshow(A)
plt.show()

This has a non-zero second-order derivative in the 'row' direction (which the scipy code appears to compute appropriately). However, its second-order derivative in the 'column' direction should be zero, as the values are constant in the column direction. However, this is not what I find when computing it as follows:

from scipy.ndimage import gaussian_filter
B = gaussian_filter(A, sigma=1, order=[2, 0], mode='reflect')
plt.imshow(B)
plt.colorbar()
plt.show()

As you can see, values on the order of 1e-5 are somehow present. This doesn't quite make sense, as (1) the solution should be just 0 everywhere, (2) the second derivative shouldn't be varying in the row direction.

Interestingly, we /do/ get the correct result if we apply two successive first-derivative operations. To make sure we use the appropriate $\sigma$ value, we can use $\sigma_\text{total} = \sqrt{\sigma_\text{operation 1}^2 + \sigma_\text{operation 2}^2}$, so I'll choose $\sigma_\text{operation 1}=\sigma_\text{operation 2}=\sqrt{1/2}$. That gives the following code:

im = gaussian_filter(A, sigma=np.sqrt(1/2), order=[1, 0], mode='reflect')
im = gaussian_filter(im, sigma=np.sqrt(1/2), order=[1, 0], mode='reflect')
plt.imshow(im)
plt.colorbar()
plt.show()

This result simply has zeroes everywhere, which is what I would expect.

My question: am I right in assuming the scipy operation contains a bug somehow? What would be the likely error? It is not the same as https://stackoverflow.com/questions/45255265/unexpected-behavior-of-gaussian-filtering-with-scipy, which is specifically related to using a cval incompatible with the array. I use mode='reflect' to avoid this issue.

Answer 1

Ndimage generates a Gaussian kernel by sampling a Gaussian and normalizing it to 1. The derivative of this kernel is generated by modifying that normalized kernel according to the chain rule to compute the derivative, and this modification is applied repeatedly to obtain higher order derivatives. Here is the relevant source code.

This indeed leads to a kernel that produces imprecise second order derivatives, as described by the OP in the question and the very nice answer.

Indeed, one can not just sample a derivative of Gaussian to obtain a convolution kernel, because the Gaussian function is not band-limited, and so sampling causes aliasing. The Gaussian function is nearly bandlimited, sampling with a sample spacing of $\sigma$ leads to less than 1% of the energy being aliased. But as the order of the derivative increases, so does the bandlimit, meaning that the higher the derivative order, the more sampling error we get. [Note that if we had no sampling error, convolution with the sampled kernel would yield the same result as a convolution in the continuous domain.]

So we need some tricks to make the Gaussian derivatives more precise. In DIPlib (disclosure: I'm an author) the second order derivative is computed as follows (see the relevant bit of source code):

1. Sample the second order derivative of the Gaussian function.
2. Subtract the mean, to ensure that the response to a contant function is 0.
3. Normalize such that the kernel, multiplied by a unit parabola ($x^2$), sums to 2, to ensure that the response to a parabolic function has the right magnitude.

As we can see, the error in this case is within numerical precision of the double-precision floating-point values:

import diplib as dip
import numpy as np
import matplotlib.pyplot as plt

A = np.zeros((30,30))
A[:,15] = 1

B = dip.Gauss(A, [1], [0,2])  # note dimension order is always [x,y] in DIPlib

plt.imshow(B)
plt.colorbar()
plt.show()

The DIPlib code to generate a 1D second order derivative of the Gaussian is equivalent to the following Python code:

import numpy as np

sigma = 2.0
radius = np.ceil(4.0 * sigma)

x = np.arange(-radius, radius+1)
g = np.exp(-(x**2) / (2.0 * sigma**2)) / (np.sqrt(2.0 * np.pi) * sigma)
g2 = g * (x**2 / sigma**2 - 1) / sigma**2
g2 -= np.mean(g2)
g2 /= np.sum(g2 * x**2) / 2.0

Answer 2

I think I have figured out the answer. The results have to do with the difference between the integral of a continuous Gaussian function, versus the sum of a discretized Gaussian.

In the ideal case, continuous and discrete sums of Gaussians already differ

As one illustrative example, consider the integral of a Gaussian versus the discrete summation of a Gaussian with unit steps is $$ \int_{-\infty}^\infty e^{-x^2}\!\text{d}x=\sqrt{\pi}\approx 1.77245,\qquad \sum_{x\in\mathbb{Z}} e^{-x^2} = \theta_3(0,e^{-1})\approx 1.77264.$$ Apparently, the two expressions differ in the details (after 4 decimals). Such things matter when we compare a continuous convolution ($(f*g)(x)=\int_\tau f(x-\tau)g(\tau)\text{d}\tau $) against a discrete convolution ($(f*g)[x]=\sum_\tau f[x-\tau]g[\tau]$) using Gaussian kernels.

Symmetric versus asymmetric functions

A Gaussian is a symmetric function, i.e., $G(x)=G(-x)$. The first derivative of a Gaussian is antisymmetric, i.e., $\text{d} G(x)/\text{d} x=-\text{d} G(-x)/\text{d} x$. The second derivative of a Gaussian is symmetric again, and so on.

Implication for convolutions

We consider truncated (to $2T+1$ terms, although we include $T\to\infty$ in this description) convolutions with Gaussians $G_0$ or their $n$-th derivatives ($G_n$) with the following computation: $$ (f*G_n)[x]= \sum_{\tau=-T}^T f[x-\tau] G_n[\tau]. $$ Knowing that Gaussian functions are either anti-symmetric or symmetric, that means we can 'simplify' the summation to only consider non-negative integer values for $\tau$, $$ (f*G_n)[x] = f[x]G_n[0] + \sum_{\tau=1}^T \left(f[x-\tau]+(-1)^n f[x+\tau]\right)G_n[\tau].$$

We now note that for an odd (i.e., first-order, third-order, ...) derivative of a Gaussian, the antisymmetry also implies that $G_n[0]=0$. This is not true for an even derivative of a Gaussian, where $G_n[0]\neq 0$.

Specific example for constant functions, $f=1$

We can thus symbolically write down what happens when convolving a Gaussian derivative with a constant function $f=1$. For a first derivative, with $n=1$, we compute $$ (f*G_1)[x] = \underbrace{G_1[0]}_{=0} + \sum_{\tau=1}^T \underbrace{(1-(-1)^11)}_{=0}G_1[\tau] = 0, $$ that is, it evaluates to 0, which we expect from differentiating a constant function.

Conversely, for a second derivative, $n=2$, it results in $$ (f*G_2)[x] = G_2[0] + \sum_{\tau=1}^T\underbrace{(1+(-1)^2 1)}_{=2} G_2[\tau] = G_2[0] + 2\sum_{\tau=1}^T G_2[\tau].$$ This expression only equals the expected result of 0, if the truncated and discrete kernel $G_2$ sums to 0, which is alternatively expressed as the expression above. This strict requirement is not met for discrete sums of a Gaussian derivative, even without truncation (i.e., when letting $T\to\infty$ the sum is still not 0)!

Explaining the observed behavior

So, I can now answer the question posed above (which was 'why do two successive 1st-order derivatives give the correct result, but a single 2nd-order derivative not?'): after applying the first derivative operation, only zeroes remain (per the result for $G_1$ described above), which leads to field of zeroes for the second convolution operation. Conversely, an immediate application of $G_2$ simply returns the sum of the kernel $G_2$ multiplied with the local value. That second definition should be 0 in the ideal and continuous case -- but because we're working with a discrete application of convolutions it ceases to hold. Even if we use a very large truncation (e.g., specify truncate=1000 when calling the Gaussian filter in Python), it will not converge to 0. Instead, it will converge to $G_0[0]\sum_{\tau=-T}^T G_2[\tau]=-8.42695\times10^{-8}$ (the term $G_0[0]$ comes from the Gaussian which is applied in the horizontal direction), which we can confirm with our code:

import numpy as np
import matplotlib.pyplot as plt
A = np.zeros((30,30))
A[:,15] = 1
from scipy.ndimage import gaussian_filter
B = gaussian_filter(A, sigma=1, order=[2, 0], mode='reflect', truncate=10000)
plt.imshow(B)
plt.colorbar()
plt.show()
print(B[0,15]) # >> -8.426948379372423e-08

So, to conclude, the scipy implementation is correct, and it does what it is supposed to do, from a discrete point of view. However, from a continuous point of view, the results are slightly unexpected and seem wrong. I'll ponder what the most appropriate fix is. Perhaps, we must simply ensure that the summation convention holds, by slightly modifying the kernel $G_2[0]$, after which the discrete kernel will mimic the continuous kernel, at least for constant functions!