Normalizing derivatives of Gaussian kernels
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. Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.
On Tue, Jul 25, 2017 at 9:58 AM, Jaime Fernández del Río < jaime.frio@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
participants (2)
-
Jaime Fernández del Río -
Ralf Gommers