我有一个问题，我需要将一个非常大的二维数组（磁盘上的一个文件）与一个适合内存的较小数组进行卷积。 scipy.signal.fftconvolve 在数组适合内存时很棒，但在不适合时无济于事。除了遍历每个数组中的所有点来手动计算卷积之外，还有其他合理的方法吗？我的数学不太好，但我想知道fftconvolve 是否可以分成几部分并通过一点重叠重新组合？还有什么？
【问题讨论】：
标签： scipy numerical-methods convolution
我可以建议你两种不同的方法（虽然我不会冒险提供一些示例代码，希望你不会介意弄清楚）：
1) 使用numpy.memmap； “内存映射文件用于访问磁盘上大文件的小片段，而不是将整个文件读入内存。(...) memmap 对象可以在任何接受 ndarray 的地方使用。”
2) 将大数组拆分为瓦片，与mode='full' 进行卷积，并将结果叠加。对于每个图块，您将在图块周围获得一个“边框”，其宽度与卷积核的宽度相同。
可以将这两种方法结合起来（例如，从一个 memmapped 文件中读取切片，并将结果叠加到另一个 memmapped 文件中，即结果）。
【讨论】：
fftconvolve 直接应用于完整的 memmap 数组？我已经在 C# 中实现了卷积，并且算法本身不会随着基本数组的增加而增加内存使用量（是的，它取决于内核大小）。所以我相信应该可以只使用memmap，而不需要平铺/分割大数组......
fftconvolve，小心地将结果添加回现有数组。它很光滑。
为了帮助 heltonbiker 的响应，快速完成它很重要。就像重新组装你的“瓷砖”一样。如果您无法将大数组加载到内存中，则需要将其作为 memmapped 文件加载...但要将其作为 memmapped 文件，您需要先创建一个。
这是一个粗略的伪代码。
#you need to know how big your data matrix is
#todo: assign nrows and ncols based on your data.
fp = np.memmap('your_new_memmap_filename.dat', dtype='float32', mode='w+', shape=(nrows, ncols))#create file with appropriate dimensions
data = open('yourdatafile.txt', 'r')
i = 0
for line in data:
arr = map(float, line.strip().split(',')) #comma delimited? header row?
fp[i, :] = arr
i += 1
del fp #see documentation...del will flush the output and close the file.
现在处理..可以继续或新脚本
convolve_matrix = somenumpyarray
fp_result = np.memmap('your_results_memmap_filename.dat', dtype='float32', mode='w+', shape=(nrows, ncols))
#open big file read only
fp = np.memmap('your_new_memmap_filename.dat', dtype='float32', mode='r', shape=(nrows, ncols))
chunksize = 10000 #?
for i in range(int(nrows/chunksize) - 1): #don't forget the remainder at the end
chunk = fp[i * chunksize: (i + 1) * chunksize, :]
res = fftconvolve(chunk, convolve_matrix)
fp_result[i * chunksize: (i + 1) * chunksize, :] = res
#deal with remainder
del fp_result
del fp
请注意，此伪代码不重叠，您需要填补一些空白。此外，一旦您开始使用拼贴，请确保您使用 Joblib 并并行处理拼贴。 https://pythonhosted.org/joblib/parallel.html 抱歉，我不能提供更多代码，我有一个为 gis 制作的 2-d tiler/reassembler，但它不在这台计算机上。它甚至可能没有太大帮助，因为您的 tiler 不会返回实际的切片，而是切片列表，可能是几个获取切片的位置（在大数组上）的列表，将其放在结果中的位置（大结果数组）以及切片的位置切片的结果（从大数组中抓取的切片的卷积结果）......遍历切片列表，处理将很容易。但是制作切片函数会很棘手。
for source_slice, result_slice, mini_slice in zip(source_slice, result_slice, mini_slice):
matrix2convolve = big_fp[source_slice[0]:source_slice[1], :]
convolve_result = fftconvolve(matrix2convolve, convolve_matrix)
big_result_fp[result_slice[0]:result_slice[1], :] = convolve_result[mini_slice[0]:mini_slice[1], :]
【讨论】：
=)。或者您可以只分配卷积核（顺便说一下，它应该更快、更容易并提供相同的结果）。但是你需要包含参数fftconvolve(mode='same')。
这听起来很幼稚，但如果你移动结果像素，它是 100% 正确的。
卷积通常是错误的。实际操作完全不需要内存。您可以简单地通读文件，执行卷积，然后将其写回同一位置。
算法在执行时的缺陷是它要求我们有一个中心像素。如果您宁愿将结果像素放在左上角，请执行卷积。您可以简单地将卷积作为扫描线操作来阅读。由于此操作没有改变向下或向右的像素，因此您的结果是正确的。
一旦打破了对先前像素的这种完全任意且毫无意义的依赖，您就可以做一些很酷的事情，例如组合卷积内核，并在当前正在从中读取它的相同内存占用空间（或文件）中执行操作。 http://godsnotwheregodsnot.blogspot.com/2015/02/combining-convolution-kernels.html
在此处查看源代码http://pastebin.com/bk0A2Z5D
将像素放在中心的原因是它们不会移动。但是，这不值得。如果您对结果周围的垃圾行相同，您真的可以将它们移回同一内存中的正确位置。事实上，如果你在右下角执行下一个卷积，并向后迭代数组，你最终会回到它们或多或少开始的地方或内核为：
0,0
0,1
简而言之，您的问题的唯一原因是前一段时间有人决定应该有一个中心像素。当你放弃那个愚蠢的想法时，你总是知道结果像素在哪里，你可以做一些有趣的操作，比如预编译卷积矩阵，从而使用组合内核一次执行所有卷积，对你来说，纯粹通过从磁盘读取来执行操作。试想一下，你可能是第一个拥有 NASA 仙女座星系图像的模糊副本的人。
如果有人知道谁首先做出了这个决定有一个中心像素的愚蠢决定的历史，我有点想知道。因为没有它，卷积只是逐像素扫描线的操作。
【讨论】：