迭代多维numpy数组中的正方形子矩阵 [英] Iterating over square submatrices in multidimensional numpy array

问题描述

我已经从图像中加载了数字高程图图像(浮点高度图)，并且正在数组中的每个2x2正方形子矩阵上进行迭代，并进行计算和求和.

I have loaded a digital elevation map image (a floating-point height map) from an image and I am iterating over each 2x2 square submatrix within the array and performing a calculation and summing the results.

此操作非常慢，因为我正在使用的高程图可能非常大(16Kx16K)，并且纯Python循环方法比numpy或scipy慢得多(或者我读过).但是，我找不到有关如何迭代多维numpy数组块的任何具体信息.

This operation is incredibly slow because the elevation maps I am working with can be extremely large (16Kx16K), and a pure Python-loop approach is vastly slower than numpy or scipy would be (or so I read). However, I can't find any concrete information on how to iterate over blocks of multidimensional numpy arrays.

例如，如果我有以下3x3 numpy数组(请记住，这可能是一个NxM数组):

For example, if I had the following 3x3 numpy array (keeping in mind this could be an NxM array):

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

我想要一个快速迭代器，其结果如下:

I would want a fast iterator that would yield something like the following:

> [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]

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

The actual order of the values within the sub-matrix are not important so long as they are consistent (ie. counter-clockwise, clockwise, zig-zag etc.)

下面是相关的代码，并且不使用numpy.

The relevant bit of code is below, and does not use 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

任何能将我指向正确方向的东西都会受到赞赏.

Anything that can point my in the right direction is appreciated.

该算法的目的是在给定高度数组和像素之间的固定采样距离的情况下，计算给定区域的表面积.该算法必须检查哪些像素组合不是空"高度，并相应地调整计算(这就是位屏蔽的作用).

The purpose of the algorithm is to calculate the surface area of a given area, given an array of heights and a fixed sampling distance between pixels. The algorithm has to check what combination of pixels are not "null" heights, and adjust the calculation accordingly (which is what the bit-masking is doing).

推荐答案

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

Using scikit-image's view_as_windows is a possible way to go:

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

编辑

如果上述方法对您无效，请

If the approach above is not valid for you, scipy.ndimage.generic_filter might do the trick:

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×2窗口返回左上角的值).

Notice that you have to change the function surface_area so that it performs the computation correctly (in my toy example it simply returns the upper left value for each 2×2 window).

这篇关于迭代多维numpy数组中的正方形子矩阵的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！