在2D numpy阵列中的子诊断上迭代

Question

在图像的浮点高度图上，我在数组中的每个2x2 Square sub-matrix上迭代，执行计算，总结结果。这很慢，因为高程图很大（16k x 16k），并且环比Numpy或Scipy慢。如何迭代多维数阵列的块？

这个3x3 numpy数组可以是nxm：

[[0.0, 1.0, 2.0],
 [3.0, 4.0, 5.0],
 [6.0, 7.0, 8.0]]

我想要一个快速的迭代器，可以产生：

> [0.0, 1.0, 3.0, 4.0]
> [1.0, 2.0, 4.0, 5.0]
> [3.0, 4.0, 6.0, 7.0]
> [4.0, 5.0, 7.0, 8.0]

sub-matrix内的值顺序只要它们保持一致（逆时针，顺时针，Zig-Zag等）就不重要。我的代码不使用numpy：

shape_dem_data = shape_dem.getdata() # shape_dem is a PIL image
for x in range(width - 1):
    for y in range(height - 1):
        i = y * width + x
        z1 = shape_dem_data[i]
        z2 = shape_dem_data[i + 1]
        z3 = shape_dem_data[i + width + 1]
        z4 = shape_dem_data[i + width]
        # Create a bit-mask indicating the available elevation data
        mask = (z1 != NULL_HEIGHT) << 3 |
               (z2 != NULL_HEIGHT) << 2 |
               (z3 != NULL_HEIGHT) << 1 |
               (z4 != NULL_HEIGHT) << 0
        if mask == 0b1111:
            # All data available.
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (gsd, gsd, z3)))
            surface_area += area_of_triangle(((0, 0, z1), (gsd, gsd, z3), (0, gsd, z4)))
            pass
        elif mask == 0b1101:
            # Top left triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (0, gsd, z4)))
        elif mask == 0b0111:
            # Bottom right triangle
            surface_area += area_of_triangle(((gsd, 0, z2), (gsd, gsd, z3), (0, gsd, z4)))
        elif mask == 0b1011:
            # Bottom left triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, gsd, z3), (0, gsd, z4)))
        elif mask == 0b1110:
            # Top right triangle
            surface_area += area_of_triangle(((0, 0, z1), (gsd, 0, z2), (gsd, gsd, z3)))
return surface_area

目的是计算高度阵列的表面积（给定像素之间的固定采样距离）。位掩模检查哪种像素的组合不是" null"高度（并相应地调整计算）。

Answer 1

使用scikit-image的 view_as_windows是一种可能的方法：

In [55]: import numpy as np
In [56]: from skimage.util import view_as_windows
In [57]: wrows, wcols = 2, 2
In [58]: img = np.arange(9).reshape(3, 3).astype(np.float64)
In [59]: img
Out[59]: 
array([[0., 1., 2.],
       [3., 4., 5.],
       [6., 7., 8.]])
In [60]: view_as_windows(img, window_shape=(wrows, wcols), step=1).reshape(-1, wrows*wcols)
Out[60]: 
array([[0., 1., 3., 4.],
       [1., 2., 4., 5.],
       [3., 4., 6., 7.],
       [4., 5., 7., 8.]])

编辑

如果上述方法对您无效，则scipy.ndimage.generic_filter可能会解决：

In [77]: from scipy.ndimage import generic_filter
In [78]: def surface_area(block):
    ...:     z1, z2, z3, z4 = block
    ...:     # YOUR CODE HERE
    ...:     return z1
    ...: 
    ...: 
In [79]: generic_filter(img, function=surface_area, 
    ...:                size=(wrows, wcols), mode='constant', cval=np.nan)
    ...: 
Out[79]: 
array([[nan, nan, nan],
       [nan,  0.,  1.],
       [nan,  3.,  4.]])

请注意，您必须更改函数surface_area，以便它正确执行计算（在我的玩具示例中，它只需返回2＆times; 2窗口的左上值）。