Дискретный лапласиан (эквивалент del2) в Python

Мне нужен эквивалент Python/Numpy дискретного лапласианского оператора (функции) del2() в Matlab (Octave). Я попробовал пару решений Python, ни одно из которых не соответствует выводу del2. На октаве у меня

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

это дает результат

   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

На питоне пробовал

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

что дает результат

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

я тоже пробовал

scipy.ndimage.filters.laplace(image)

Это дает результат

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

Таким образом, ни один из выходов не соответствует друг другу. Октавный код del2.m предполагает, что это оператор Лапласа. Я что-то упускаю?


python numpy scipy scientific-computing
person BBSysDyn    schedule 14.01.2011    source источник
comment
Внутри все операторы одинаковы (Matlab, по-видимому, делит на 4, а Python - нет). На границе вы можете сделать две версии Python одинаковыми, также предоставив mode="wrap" для laplace(). Но, просто глядя на результат Matlab, я понятия не имею, что Matlab делает на границах.   -  person Sven Marnach    schedule 15.01.2011
comment
На самом деле он выполняет кубическую экстраполяцию по краям: mathworks.it/it/help /matlab/ref/del2.html Так что, если вы попробуете последний пример с laplace(), вы не сможете получить правильный результат и на границах.   -  person astrojuanlu    schedule 08.04.2013
comment
Вы также можете попробовать этот scipy .sparse.csgraph.laplacian   -  person sh139u1032    schedule 02.07.2018
comment
github.com/ipython-books/cookbook-2nd/ blob/master/ имеет одну реализацию.   -  person Dilawar    schedule 13.07.2018

Ответы (3)


arrow_upward
10
arrow_downward

Возможно, вы ищете scipy.ndimage.filters.laplace()< /а>.

person Sven Marnach    schedule 14.01.2011
comment
Я протестировал эту функцию и по сравнению с выводом del2, она другая. - person BBSysDyn; 14.01.2011
comment
Конечно, могут быть различия в дискретизации, например, на границах. Не могли бы вы подробнее рассказать, как вы его тестировали и чем он отличается? - person Sven Marnach; 14.01.2011
comment
Я проверил a = [3 4 6;7 8 9;1 3 3]; дисп(дел2(а)). На стороне Python я назвал функцию, которую вы упомянули, результаты совершенно разные для каждой ячейки. - person BBSysDyn; 14.01.2011
comment
А, я протестировал матрицу побольше, и да, отличия на границах. Как сделать так, чтобы границы были такими же, как у del2? - person BBSysDyn; 15.01.2011
comment
@ user423805: Сначала узнайте, что del2 делает на границе (я не знаю), а затем посмотрите в связанной документации, как добиться такого же поведения в Python. (Но если вы даже не знаете, что делает del2, почему вам важно получить то же самое в Python?) - person Sven Marnach; 15.01.2011
comment
Как ни странно, после портирования del2 я увидел проблемы в другом месте, исправив их, оказалось, что вызов laplace() работает и в моем приложении. - person BBSysDyn; 16.01.2011
comment
@becko исправлено. (Поиск имени функции по-прежнему дал желаемые документы.) - person Sven Marnach; 31.01.2016

arrow_upward
7
arrow_downward

Вы можете использовать convolve для вычисления лапласиана путем свертки массива с соответствующим трафаретом:

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')
person 2daaa    schedule 18.01.2011

arrow_upward
3
arrow_downward

На основе кода здесь

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

Я попытался написать эквивалент Python. Кажется, это работает, любые отзывы будут оценены.

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
person BBSysDyn    schedule 15.01.2011
comment
Я предполагаю, что первое предложение if должно быть mc >= 3. Также я не уверен, возвращает ли он правильные значения, когда dx != 1 or dy != 1 - person Daniel; 24.02.2021