[SciPy-Dev] Normalizing derivatives of Gaussian kernels

[Ralf Gommers]

On Tue, Jul 25, 2017 at 9:58 AM, Jaime Fernández del Río <
jaime.frio at gmail.com> wrote:

> There is an interesting issue that was posted yesterday: #7644
> <https://github.com/scipy/scipy/issues/7644>
>
> in a nutshell, when ndimage does an N-dimensional convolution with a
> Gaussian kernel, we separate it into N 1-dimensional convolutions, one
> along each axis. If the border is set to 'constant' mode, we need to
> replace that constant with the result of convolving the 1D kernel with that
> constant value, something we are currently not doing.
>
> In an ideal world that would mean leaving the constant alone if the order
> of the Gaussian kernel is zero, or setting it to zero if the order is
> larger than zero, because that's what the integrals from -inf to +inf of
> those kernels are supposed to be.
>
> In real life we get things right for zero order (because we normalize the
> kernel values to add up to 1.0) and the odd-ordered higher orders, (because
> they are antisymmetric, so they naturally add up to zero).
>
> But for even-ordered Gaussian kernels, which are symmetric, truncation
> means that we get a slightly negative sum, and that e.g. convolution of a
> constant field with a second order kernel, which should return all zeros,
> returns a small negative error.
>
> There are two ways we can solve this:
>
>    1.  Make our separable approach to Gaussian filtering equivalent to
>    what one would get without treating the filter as separable, by multiplying
>    the constant value by the sum of the 1-D kernel after each 1-D convolution.
>    This will at least make ndimage consistent with itself.
>    2. Normalize the Gaussian kernels of higher order, so that their total
>    sums are equal to zero. I'm not sure if there's is a standard way of doing
>    this in the literature, my searches haven't found anything. But e.g. all
>    kernels you find as being an approximation to the Laplacian of Gaussian do
>    add up to exactly zero.
>
> Thoughts and comments are very welcome.
>

It sounds like (2) is the more correct approach. As far as I can tell
either approach will change the numerical results for even-ordered kernels
in a similar way, and making the change can be considered a bug fix rather
than a backwards compatibility change. Is that right?

Ralf
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20170727/1402a1e1/attachment.html>