How does the skimage.filters.laplace function work in Python?

Question

I am trying to figure out the kernel being used in skimage.filters's laplace function. I know that a Laplacian filter is based on matrix convolution, but I just can't seem to make sense of the values produced by skimage.filters's laplace function.

This is an example:

>>> import numpy as np
>>> import skimage
>>> from skimage.filters import laplace
>>> x = [[2,3,2],[5,3,6],[3,7,3]]
>>> x = np.array(x)
>>> x
array([[2, 3, 2],
       [5, 3, 6],
       [3, 7, 3]])
>>> laplace(x, ksize=3)
array([[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
       [ 0.00000000e+00, -9.75781955e-19,  0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00]])

If skimage.filters's laplace function uses a kernel/operator of

[[ 0,  1, 0],
 [ 1, -4, 1],
 [ 0,  1, 0]]

then according to matrix convolution, it should have produced the output of

[[ 0,  5, -1],
 [12, -9, 16],
 [ 0, 19, -1]]

instead of

[[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00],
 [ 0.00000000e+00, -9.75781955e-19,  0.00000000e+00],
 [ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00]]

I am very confused on what kernel/operator skimage.filters's laplace function is using to have almost every output value so close to zero, such as -9.75781955e-19. I honestly don't think any reasonable kernel/operator could produce this output, so maybe I am just not understanding how Python's skimage.filters's laplace function is working...

Any help/comments/suggestions/insights to this question would be greatly appreciated. Thank you.

Answer 1

welcome to the scikit-image thread of Stack Overflow! The reason for this strange behaviour is that the dtype of x is int64, and the scikit-image laplace function calls img_as_float in order to do the computation in float numbers, but when casting the dtype it also divides the array by the maximum value of the origin dtype (here 2^63 - 1), hence the very small values. If you want to avoid this problem you can convert the image to float before passing it to laplace:

>>> x = x.astype(np.float)                                                
>>> filters.laplace(x)                                                    
array([[-4.,  2., -5.],
       [ 7., -9., 10.],
       [-6., 12., -7.]])

(the function uses the default boundary-condition mode of scipy.ndimage.convolve which is 'reflect').

Note that this behaviour (dividing by the max value of the dtype) is likely to change with scikit-image 1.0, precisely because we've noticed that it can be confusing for users as in your case.