python - 如何使用 numpy 和 strides 来提高循环性能

Question

我是 numpy 的新手，刚刚听说过 strides，但我很难理解如何使用它们。我有一个非常慢的小循环，这似乎是一个可以用于跨步的好例子。

def semivariogram( ds, band, lag ):
    width = ds.RasterXSize
    height = ds.RasterYSize
    data = band.ReadAsArray( 0, 0, width, height ).astype(np.float)
    #print 'w: {}, h: {}'.format(width, height)

    sumw = 0.0
    sumh = 0.0
    for i in range(width-lag):
        for j in range(height-lag):
            sumw += data[i+lag,j]-data[i,j]
            sumh += data[i,j+lag]-data[i,j]

    Nh2 = 2.0*(width-lag)*(height-lag)

    return [sumw/Nh2, sumh/Nh2, (sumw/Nh2+sumh/Nh2)/2.0]

行:

data = band.ReadAsArray( 0, 0, width, height ).astype(np.float)

正在将图像的一个波段读取到 float 组中。 ds 是图像的句柄，例如:

ds = gdal.Open('test.tif')

Answer 1

好的，让我们一步一步来。

你有:

sumw = 0.0
sumh = 0.0
for i in range(width-lag):
    for j in range(height-lag):
        sumw += data[i+lag,j]-data[i,j]
        sumh += data[i,j+lag]-data[i,j]

为了清楚起见，让我们将其分成两个循环

sumw = 0.0
for i in range(width-lag):
    for j in range(height-lag):
        sumw += data[i+lag,j]-data[i,j]
sumh = 0.0
for j in range(height-lag):
    for i in range(width-lag):
        sumh += data[i,j+lag]-data[i,j]

我们可以写data[i,j+lag]如data[:-lag,lag:][i,j] ，和data[i,j]是 data[:-lag,:-lag][i,j] (对于我们的范围 i 和 j ，并假设 lag != 0)。所以我们的循环变成:

sumw = 0.0
for i in range(width-lag):
    for j in range(height-lag):
        sumw += data[lag:,:-lag][i,j]-data[:-lag,:-lag][i,j]
sumh = 0.0
for j in range(height-lag):
    for i in range(width-lag):
        sumh += data[:-lag,lag:][i,j]-data[:-lag,:-lag][i,j]

但现在我们注意到所有循环都只是迭代一个简单的索引 [i,j] ，这样我们就可以将它们完全展平:

sumw = (data[lag:,:-lag] - data[:-lag,:-lag]).sum()
sumh = (data[:-lag,lag:] - data[:-lag,:-lag]).sum()

通过发现差值之和就是总和之差，这可以变得更快:

sumw = data[lag:,:-lag].sum() - data[:-lag,:-lag].sum()
sumh = data[:-lag,lag:].sum() - data[:-lag,:-lag].sum()

此时，您可以直观地思考整个任务:

您想要sumw作为绿色轮廓上的总和减去红色轮廓上的总和，以及 sumh作为蓝色轮廓上的总和减去红色轮廓上的总和。

在这里，我们可以进行以下两种简化之一:

1. 我们注意到第一个和中的大部分元素都从第二个和中删除了。事实上，唯一不属于这种情况的元素是 [:lag]和[-lag:] 。所以这就变成了:

   sumw = data[-lag:,:-lag].sum() - data[:lag,:-lag].sum()
   sumh = data[:-lag,-lag:].sum() - data[:-lag,:lag].sum()
   

   如果 lag < len(data) / 2 的话，速度会更快。在这里，我们设置sumw = green - (red + magenta) ，和sumh = blue - (red + yellow)

2. 我们注意到第二个和是重复的，所以有

   sum_shared = data[:-lag,:-lag].sum()
   sumw = data[lag:,:-lag].sum() - sum_shared 
   sumh = data[:-lag,lag:].sum() - sum_shared