Discrete Laplacian (del2 equivalent) in Python

Question

I need the Python / Numpy equivalent of Matlab (Octave) discrete Laplacian operator (function) del2(). I tried couple Python solutions, none of which seem to match the output of del2. On Octave I have

image = [3 4 6 7; 8 9 10 11; 12 13 14 15;16 17 18 19]
del2(image)

this gives the result

   0.25000  -0.25000  -0.25000  -0.75000
  -0.25000  -0.25000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000
   0.25000   0.25000   0.00000   0.00000

On Python I tried

import numpy as np
from scipy import ndimage
import scipy.ndimage.filters
    
image =  np.array([[ 3,  4,  6,  7],
                   [ 8,  9, 10, 11],
                   [12, 13, 14, 15],
                   [16, 17, 18, 19]])
stencil = np.array([[0,  1, 0],
                    [1, -4, 1],
                    [0,  1, 0]])
print ndimage.convolve(image, stencil, mode='wrap')

which gives the result

[[ 23  19  15  11]
 [  3  -1   0  -4]
 [  4   0   0  -4]
 [-13 -17 -16 -20]]

I also tried

scipy.ndimage.filters.laplace(image)

That gives the result

[[ 6  6  3  3]
 [ 0 -1  0 -1]
 [ 1  0  0 -1]
 [-3 -4 -4 -5]]

So none of the outputs seem to match eachother. Octave code del2.m suggests that it is a Laplacian operator. Am I missing something?

Answer 1

I was also looking for a function to compute the Laplacian in Python. Therefore, I made a comparison with a Laplacian computed as suggested by Sven using scipy.ndimage.laplace, and a "custom" version made by iterating the use of numpy.gradient a couple of times.

I came out with the following piece of code

import numpy as np
import scipy.ndimage as sn

import matplotlib.pylab as pl
from mpl_toolkits.axes_grid1 import make_axes_locatable


def lapl(mat, dx, dy):
    """
    Compute the laplacian using `numpy.gradient` twice.
    """
    grad_y, grad_x = np.gradient(mat, dy, dx)
    grad_xx = np.gradient(grad_x, dx, axis=1)
    grad_yy = np.gradient(grad_y, dy, axis=0)
    return(grad_xx + grad_yy)
        
# The geometry of the example
dx = 0.1
dy = 0.1
lx = 20
ly = 15

x = np.arange(0, lx, dx)
y = np.arange(0, ly, dy)

x0 = 3*lx/4
y0 = 3*ly/4 

xx, yy = np.meshgrid(x, y)
zz = np.exp(-((xx-x0)**2+(yy-y0)**2))

#
# Plot
#
fig, ax = pl.subplots(1,3, sharex=True, sharey=True)
# Plot data
ax[0].imshow(zz)
# Plot Laplacian "scipy" way
im = ax[1].imshow(sn.laplace(zz)/(dx*dy))
divider = make_axes_locatable(ax[1])
cax = divider.append_axes("right", size="5%", pad=0.05)
pl.colorbar(im, cax=cax)
# Plot Laplacian "numpy.gradient twice" way
im = ax[2].imshow(lapl(zz, dx, dy))
divider = make_axes_locatable(ax[2])
cax = divider.append_axes("right", size="5%", pad=0.05)
pl.colorbar(im, cax=cax)

pl.show()

This seems somehow simpler than some of the solutions posted here and exploits for the computation some numpy functions, therefore I expect it to be a faster.

Answer 2

Maybe you are looking for scipy.ndimage.filters.laplace() (which in 2023 is scipy.ndimage.laplace).

Answer 3

# Laplacian operator (2nd order cetral-finite differences)
# dx, dy: sampling, w: 2D numpy array

def laplacian(dx, dy, w):
    """ Calculate the laplacian of the array w=[] """
    laplacian_xy = np.zeros(w.shape)
    for y in range(w.shape[1]-1):
        laplacian_xy[:, y] = (1/dy)**2 * ( w[:, y+1] - 2*w[:, y] + w[:, y-1] )
    for x in range(w.shape[0]-1):
        laplacian_xy[x, :] = laplacian_xy[x, :] + (1/dx)**2 * ( w[x+1,:] - 2*w[x,:] + w[x-1,:] )
    return laplacian_xy

Answer 4

You can use convolve to calculate the laplacian by convolving the array with the appropriate stencil:

from scipy.ndimage import convolve
stencil= (1.0/(12.0*dL*dL))*np.array(
        [[0,0,-1,0,0], 
         [0,0,16,0,0], 
         [-1,16,-60,16,-1], 
         [0,0,16,0,0], 
         [0,0,-1,0,0]])
convolve(e2, stencil, mode='wrap')

Answer 5

Based on the code here

http://cns.bu.edu/~tanc/pub/matlab_octave_compliance/datafun/del2.m

I attempted to write a Python equivalent. It seems to work, any feedback will be appreciated.

import numpy as np

def del2(M):
    dx = 1
    dy = 1
    rows, cols = M.shape
    dx = dx * np.ones ((1, cols - 1))
    dy = dy * np.ones ((rows-1, 1))

    mr, mc = M.shape
    D = np.zeros ((mr, mc))

    if (mr >= 3):
        ## x direction
        ## left and right boundary
        D[:, 0] = (M[:, 0] - 2 * M[:, 1] + M[:, 2]) / (dx[:,0] * dx[:,1])
        D[:, mc-1] = (M[:, mc - 3] - 2 * M[:, mc - 2] + M[:, mc-1]) \
            / (dx[:,mc - 3] * dx[:,mc - 2])

        ## interior points
        tmp1 = D[:, 1:mc - 1] 
        tmp2 = (M[:, 2:mc] - 2 * M[:, 1:mc - 1] + M[:, 0:mc - 2])
        tmp3 = np.kron (dx[:,0:mc -2] * dx[:,1:mc - 1], np.ones ((mr, 1)))
        D[:, 1:mc - 1] = tmp1 + tmp2 / tmp3

    if (mr >= 3):
        ## y direction
        ## top and bottom boundary
        D[0, :] = D[0,:]  + \
            (M[0, :] - 2 * M[1, :] + M[2, :] ) / (dy[0,:] * dy[1,:])

        D[mr-1, :] = D[mr-1, :] \
            + (M[mr-3,:] - 2 * M[mr-2, :] + M[mr-1, :]) \
            / (dy[mr-3,:] * dx[:,mr-2])

        ## interior points
        tmp1 = D[1:mr-1, :] 
        tmp2 = (M[2:mr, :] - 2 * M[1:mr - 1, :] + M[0:mr-2, :])
        tmp3 = np.kron (dy[0:mr-2,:] * dy[1:mr-1,:], np.ones ((1, mc)))
        D[1:mr-1, :] = tmp1 + tmp2 / tmp3

    return D / 4