多维numpy数组中方子矩阵的迭代

Question

我从一个图像加载了一个数字高程图图像（一个浮点高度图），我在数组中的每个2x2平方子矩阵上迭代，并执行计算和结果汇总。在

这个操作非常慢，因为我使用的立面图可能非常大（16Kx16K），而纯Python循环方法比numpy或scipy慢得多（或者我读过）。但是，我找不到任何关于如何迭代多维numy数组块的具体信息。在

例如，如果我有以下3x3 numpy数组（记住这可能是NxM数组）：

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

我想要一个快速迭代器，它可以产生如下结果：

^{pr2}$

子矩阵中值的实际顺序并不重要，只要它们是一致的（即逆时针、顺时针、之字形等）

下面是相关的代码位，不使用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

任何能为我指明正确方向的东西都是值得感激的。在

编辑

该算法的目的是计算给定区域的表面积，给定一个高度数组和像素之间的固定采样距离。算法必须检查哪些像素组合不是“空”高度，并相应地调整计算（这就是位掩蔽的作用）。在

Answer 1

使用scikit image的^{}是一种可能的方法：

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.]])

编辑

如果上述方法对您无效，^{}可能会起作用：

^{pr2}$

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