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?

Leastwise

In the interior, the operators are all the same (Matlab apparently divides by 4 where Python does not). On the boundary, you can make the two Python versions the same by also providing mode="wrap" to laplace(). But by just looking at the Matlab result, I have no idea what Matlab does on the boundaries. – Kile

Actually it does cubic extrapolation on the edges: mathworks.it/it/help/matlab/ref/del2.html So if you try the final example with laplace() there's no way to get the right result on the boundaries too. – Sentience

You could also try this scipy.sparse.csgraph.laplacian – Eyeleteer

github.com/ipython-books/cookbook-2nd/blob/master/… has one implementation. – Halloween

Answer 1

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

Kile

I tested this function and compared to del2 output, it's different. – Leastwise

Of cause there can be differences in the discretisation, for example on the boundaries. Could you be more specific how you tested it and how it is different? – Kile

I tested a = [3 4 6;7 8 9;1 3 3]; disp(del2(a)). On the Python side I called the function you mentioned, the results are completely different, on each cell. – Leastwise

Ah, I tested a bigger matrix, and yes, the differences are on the boundaries. How can I have the boundaries come out the same as del2? – Leastwise

@user423805: First find out what del2 does on the boundary (I don't know), then look up in the linked documentation how to get the same behaviour in Python. (But if you don't even know what del2 does, why is it important to you to get the same in Python?) – Kile

Oddly enough, after porting del2, I saw problems were elsewhere, fixing those, it turned out laplace() call also works in my application. – Leastwise

@becko Fixed. (Googling for the function name still yielded the desired docs.) – Kile

Answer 2

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')

Terris

is e2 the name of the input array (2D)? is dL spacing between points of the grid? – Gayden

Answer 3

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

Leastwise

I guess the first if-clause should be mc >= 3. Also I'm not sure if it returns the correct values when dx != 1 or dy != 1 – Verena

Answer 4

# 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

Bonnybonnyclabber

Answer 5

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.

Sheaff