Python语言中的离散拉普拉斯(相当于del2)
Python语言中的离散拉普拉斯(相当于del2)
提问于 2011-01-14 22:41:05
我需要Python / Numpy等效于Matlab (Octave)离散拉普拉斯算子(函数) del2()。我尝试了几个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在Python上我尝试过
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表明它是拉普拉斯运算符。我是不是遗漏了什么？
回答 4
Stack Overflow用户
回答已采纳
发布于 2011-01-14 22:44:30
Stack Overflow用户
发布于 2011-01-18 22:25:11
您可以使用卷积通过将数组与适当的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')Stack Overflow用户
发布于 2011-01-15 22:06:00
基于这里的代码
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页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/4692196
复制相关文章
相似问题