高效的复杂掩蔽移动窗口分析

Question

我正在 Python 中开发一个移动窗口算法，该算法将用于在大型 numpy 数组中滑动（例如，我的测试数组尺寸为 6349x9849）。我需要在 25 x 25 移动窗口中计算多个不同的统计数据，该移动窗口在每个像素位置被 9x9 窗口遮盖。

还有一个警告阻止了我使用 convolve（例如，这对于计算移动窗口中的平均值非常快）：如果移动窗口中的中心像素为 0，我将统计值设置为 -9999作为标志，或者如果 25x25 移动窗口包含超过一半的 0 值，我将统计值设置为 -1 作为标志。然后我可以稍后处理这些标志。

我已经编写了 Python 代码来执行此操作，但这是我在 Python 学习中第一次处理特别大量的数据，因此遇到了效率问题 - 我编写的代码需要很长时间已经写过（对于这个图像大小，每个统计大约需要 6 个小时......！）。

我想请教如何更有效地做到这一点。我想先优化代码，然后再投入更多的计算能力（我相信我可以使用多处理例如执行此操作的模块）。

我的代码，例如统计数据（标准偏差）如下。我对每个统计数据重复这一行（我有 6 个统计数据要计算）：

# Calculate the standard deviation of the masked moving window
stats_std = [-9999 if ds_array[row,column] == 0
else -1 if np.count_nonzero(ds_array[row-border_buff:row+border_buff+1,column-border_buff:column+border_buff+1]) < (outer_box**2)/2
else np.std([i for i in np.ma.compressed(np.ma.masked_array(ds_array[row-border_buff:row+border_buff+1,column-border_buff:column+border_buff+1],mask)) if i!=0])
for row in range(border_buff,m-border_buff)
for column in range(border_buff,n-border_buff)]

# Reshape the list into the image dimensions
stats_std = np.reshape(stats_std,(m-2*border_buff,n-2*border_buff))

如果需要，我可以提供一个示例子集数组和所需的输出以供试用，但不确定提供此功能的最佳方式，如果需要使上述内容更清晰，请告诉我。

附：我尝试了一种将 2d 数组转换为数组所有 25x25 子集的 3d 的方法，然后在每个子集的一个步骤中有效地计算 6 个统计数据，希望能够显着节省计算量。但这导致了 MemoryError。

Answer 1

pyvips 可让您在巨大的图像efficiently and using very little memory 上计算复杂的事物。它是 LGPL，可在 Linux、macOS 和 Windows 上运行，并且适用于 Python 的每个版本。大多数 linux 在包管理器中都有它。

当你在 pyvips 中进行a + b 之类的操作时，它实际上并没有进行任何处理，它只是将另一个节点添加到图像处理操作图中。当您最终将结果写入某处时，整个图形都会进行评估，并将图像以一组小块并行传输到您的系统中。

因为中间图像并不真正存在，你只需要少量的内存，而且因为它是并行的，所以速度很快。

例如，您可以像这样进行 sdev 计算：

import sys
import pyvips

# load the input image ... the access hint means we will only make a single
# top-to-bottom pass over image, and it can therefore be streamed 
image = pyvips.Image.new_from_file(sys.argv[1], access='sequential')

# our convolution ... total pixels in an M x M window
# it's a simple box filter, so we can use a seperable convolution 
# (two 1D filters at 90 degrees)
window_size = 25
size = window_size * window_size
sum_mask = pyvips.Image.new_from_array([1] * window_size)

# standard deviation ... sum and sum of squares
s = image.convsep(sum_mask)
s2 = (image * image).convsep(sum_mask)
sdev = ((s2 - (s * s / size)).abs() / (size - 1)) ** 0.5

# find all zero input pixels ... these become -9999 in the output
sdev = (image == 0).ifthenelse(-9999, sdev)

# find all pixels where more than half of the window is zero ... these become
# -1 in the output
# pyvips uses 255 for TRUE and 0 for FALSE
more_than_half_zero = (image == 0).convsep(sum_mask) > 255 * size / 2
sdev = more_than_half_zero.ifthenelse(-1, sdev)

sdev.write_to_file(sys.argv[2])

我可以这样运行：

$ vipsheader x.jpg 
x.jpg: 10000x10000 uchar, 1 band, b-w, jpegload
$ /usr/bin/time -f %M:%e python3 sdev.py x.jpg x.pfm
81432:4.11

从 10,000 x 10,000 像素的 jpg 图像制作 10,000 x 10,000 像素的 PFM（一种可以存储浮点值的简单格式），其中 PFM 中的每个像素都是相应 25 x 25 窗口的 sdev，加上您的规则为零。这个桌面需要4s，最大内存81mb。

我可以通过只使用一个线程来减少内存使用，但当然会慢很多：

$ VIPS_CONCURRENCY=1 /usr/bin/time -f %M:%e python3 sdev.py x.jpg x.pfm
54128:16.92

现在只有 54mb 的内存，但是 17s 的时间。

您也可以在 numpy 数组中读取和写入图像，请参阅chapter in the docs。