OpenCV filter2d gives incorrect result

Question

I am currently trying to filter an image with a Laplacian kernel that I have constructed myself. However, when filtering an input image with this kernel it gives an unexpected result compared to the implementation in SciPy.

The Laplacian kernel I have constructed should be verified by the following images

The code for filtering the image:

im = cv2.imread("test.png",0)
im = im.astype(np.float32)

def lkern(t=1.):
    ax = np.arange(np.round(-5*np.sqrt(t),0),np.round(5*np.sqrt(t),0)+1)
    xx, yy = np.meshgrid(ax, ax)

    kernel = -1/(np.sqrt(2*np.pi*t)*t)*np.exp(-(xx**2+yy**2)/(2*t))+
        (xx**2+yy**2)/(np.sqrt(2*np.pi*t)*t**2)*np.exp(-(xx**2+yy**2)/(2*t))


    return kernel.astype(np.float)

t = 25**2/2
l = lkern(t)

L = cv2.filter2D(im/255,-1,l)

plt.figure()
plt.imshow(L,cmap="gray")
plt.show()

which results in

Comparing to SciPy's ndimage.gaussian_laplace, the result should have been

which is very different, and I cannot figure out how to do this properly.

Answer 1

The code in the OP seems to take the equation for a 1D Laplace of Gaussian and use that to construct a 2D radially symmetric function. That is, along any of the diameters of the kernel the function looks like a 1D Laplace of Gaussian. This is not the correct way to create a 2D Laplace of Gaussian.

The Laplace of Gaussian is defined as the sum of the second order derivative of the Gaussian kernel along each of the axes. That is,

LoG = d²/dx² G + d²/dy² G

With G the Gaussian kernel.

With Numpy, you can construct this kernel as follows. I'm using the separability of the Gaussian to reduce the computational complexity.

s = 5;
x = np.arange(np.floor(-4*s),np.ceil(4*s)+1)
g = 1/(np.sqrt(2*np.pi)*s)*np.exp(-x**2/(2*s**2))
d2g = (x**2 - s**2)/(s**4) * g
log = g * d2g[:,None] + g[:,None] * d2g

Tricks here: g and d2g are 1D functions. g[:,None] turns the 1D function on its side, so that the multiplication causes broadcasting, leading to a 2D output.

I have written the kernel this way, rather than expressing the full 2D equation in one go, because this leads to large efficiencies in your code: The convolution of an image f with the kernel log can be written as:

conv(f, log) = conv(f, g * d2g[:,None] + g[:,None] * d2g)
             = conv(conv(f, g), d2g[:,None]) + conv(conv(f, g[:,None]), d2g)

That is, instead of one convolution with a large 2D kernel, we compute 4 convolutions with relatively small 1D kernels. Note that the actual order here does not matter:

• One applies a 1D kernel g and on the result a 1D kernel d2g along the other axis. These two operations can be reversed.
• Then one repeats this process changing the axes along which each of the operations is applied.
• Finally one adds the two results.

(It is OK to use cv2.filter2D where I wrote conv. conv just indicates any convolution function, but a correlation function like filter2D is fine because the kernels are all symmetric.)