【问题标题】:How to evaluate the sum of values within array blocks如何评估数组块内的值的总和
【发布时间】:2016-04-03 08:16:59
【问题描述】:

我有数据数组,形状为 100x100。我想把它分成 5x5 的块,每个块有 20x20 的网格。我想要的每个块的值是其中所有值的总和。

有没有更优雅的方式来完成它?

x = np.arange(100)
y = np.arange(100)
X, Y = np.meshgrid(x, y)
Z = np.cos(X)*np.sin(Y)
Z_new = np.zeros((5, 5))
for i in range(5):
  for j in range(5):
    Z_new[i, j] = np.sum(Z[i*20:20+i*20, j*20:20+j*20])

这是基于index的,如果基于x呢?

x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
X, Y = np.meshgrid(x, y)
Z = np.cos(X)*np.sin(Y)
x_new = np.linspace(0, 1, 15)
y_new = np.linspace(0, 1, 15)

Z_new?

【问题讨论】:

标签: python numpy


【解决方案1】:

只需 reshape 将这两个轴中的每一个拆分为两个形状为 (5,20) 的轴,以形成一个 4D 数组,然后沿长度为 20 的轴求和减少,就像这样 -

Z_new = Z.reshape(5,20,5,20).sum(axis=(1,3))

np.einsum 的功能相同,但可能更快的选项 -

Z_new = np.einsum('ijkl->ik',Z.reshape(5,20,5,20))

通用块大小

扩展到一般情况 -

H,W = 5,5 # block-size
m,n = Z.shape
Z_new = Z.reshape(H,m//H,W,n//W).sum(axis=(1,3))

einsum 变成 -

Z_new = np.einsum('ijkl->ik',Z.reshape(H,m//H,W,n//W))

要计算块之间的平均值/平均值,请使用 mean 而不是 sum 方法。

通用块大小和缩减操作

扩展为使用具有ufuncsreduction 操作支持多个axes 参数和axis 进行缩减,这将是-

def blockwise_reduction(a, height, width, reduction_func=np.sum):
    m,n = a.shape
    a4D = a.reshape(height,m//height,width,n//width)
    return reduction_func(a4D,axis=(1,3))

因此,要解决我们的具体情况,应该是:

blockwise_reduction(Z, height=5, width=5)

对于逐块平均计算,它将是 -

blockwise_reduction(Z, height=5, width=5, reduction_func=np.mean)

【讨论】:

    【解决方案2】:

    您可以执行以下操作。

    t = np.eye(5).repeat(20, axis=1)
    Z_new = t.dot(Z).dot(t.T)
    

    这是正确的,因为Z_new[i, j] = t[i, k] * Z[k, l] * t[j, l]

    而且这似乎比 Divakar 的解决方案更快。

    【讨论】:

      【解决方案3】:

      这样的问题非常适合像scipy.ndimage.measurements.sum 这样的函数,因为它允许“分组”和“标记”术语。你会得到你想要的东西:

      labels = [[20*(y//5) + x//5 for x in range(100)] for y in range(100)]
      s = scipy.ndimage.measurements.sum(Z, labels, range(400))
      

      (未测试,但就是这个想法)。

      【讨论】:

        猜你喜欢
        • 2016-12-25
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多