Python 不强制周期性边界的二维卷积

Question

我正在为一个疾病问题建模，其中二维景观中的每个人都具有由（径向基）核函数描述的传递率。我的目标是将内核与人口密度进行卷积，以便输出捕获整个景观的传播风险。

我使用 NumPy 的 2D FFT 和反 FFT 函数执行了卷积。但是，这会在结果中强制使用 periodic/wrapped 边界条件，这不适合我的模型。有没有一种方法可以在原始固定边界的上下文中进行卷积？

import numpy as np
import random
from math import *
import matplotlib.pyplot as plt

''' Landscape parameters ''' 
L = 10.
nx = 100
dx = L/nx
hs = .5 * dx  # half-step
ulist = np.linspace(hs, L-hs, nx)

''' Radial Basis Function Kernel '''
alpha = 1.
i, j = ulist.reshape(nx,1), ulist.reshape(1,nx)
r = np.minimum(i-ulist[0], L-i+ulist[0])**2 + np.minimum(j-ulist[0], L-j+ulist[0])**2
rbf = sqrt(1 / (2 * alpha ** 2))
ker = np.exp(-(rbf * r) ** 2)
ker = ker/np.sum(ker)

''' Population Density '''
ido = np.random.randint(nx, size=(1000,2)).astype(np.int)
og = np.zeros((nx,nx))
np.add.at(og, (ido[:,0], ido[:,1]), 1)

''' Convolution via FFT and inverse-FFT '''
v1 = np.fft.fft2(ker)
v2 = np.fft.fft2(og)
v0 = np.fft.ifft2(v1*v2)
dd = np.abs(v0)

plt.plot(ido[:,1], ido[:,0], 'ko', alpha=.5)
plt.imshow(dd, origin='origin')
plt.show()

Answer 1

为了做你想做的，你需要用零填充og，并相应地扩展ker（因为它已经循环移位，你只需要扩展它更多以适应og).

的大小

填充 og (300x300)：原始 ker (100x100)：展开 ker (300x300)：

代码：

pad = 100
og = np.pad(og, pad, mode='constant')
new_ker = np.zeros_like(og)
new_ker[:nx//2,:nx//2] = ker[:nx//2,:nx//2]
new_ker[:nx//2,-nx//2:] = ker[:nx//2,-nx//2:]
new_ker[-nx//2:,:nx//2] = ker[-nx//2:,:nx//2]
new_ker[-nx//2:,-nx//2:] = ker[-nx//2:,-nx//2:]
ker = new_ker

# doesn't change
v1 = np.fft.fft2(ker)
v2 = np.fft.fft2(og)
v0 = np.fft.ifft2(v1*v2)

# unpadding the result
v0 = v0[pad:-pad,pad:-pad]

(展开比较乱，直接生成内核就可以了，我只是想把展开的部分分开。)

原结果：填充结果：

如果您不介意使用 scipy，它可以为您完成这一切（包括边界条件的其他变体）。有关更多选项，请参阅 fftconvolve or convolve2d。请注意，您不必必须 pre-roll 您的内核才能使用这些函数（也就是说，最大值应该在中心，而不是在角落）。示例代码（产生与上面手动填充相同的结果）：

# unroll the kernel (again, can be done directly when generating `ker`)
ker = np.roll(ker, 50, 0)
ker = np.roll(ker, 50, 1)

from scipy.signal import fftconvolve
dd = fftconvolve(og, ker, mode='same')

Answer 2

每当你使用基于 fft 的方法时，你总是在问题中引入周期性，本质上，根据定义。

如果你不想要这个，你应该首先在计算级别之前在数学建模级别定义你想要的是什么。

记住一些根据我的经验在很多情况下都非常有用的技巧：

-- 将您的问题嵌入更大的 2d space 以减少周期性影响。换句话说，尽可能地推动周期性边界以模拟 'free-space' 情况。

-- 使用离散余弦变换代替。余弦变换具有使边界条件对称化的作用。它的工作原理就好像在每个边界之外都有您的数据的镜像反射。在某些问题中，这用于模拟反射边界条件。

Python 不强制周期性边界的二维卷积

Python 2D convolution without forcing periodic boundaries

python

numpy

fft

convolution