当前位置: 首页 > 知识库问答 >
问题:

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

谢英光
2023-03-14

我已经从图像中加载了一个数字高程图图像(一个浮点高度图),我正在数组中的每个2x2平方子矩阵上迭代,并执行计算并对结果求和。

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

例如,如果我有以下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]

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

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

任何能给我指明正确方向的东西都很感谢。

该算法的目的是计算给定区域的表面积,给定高度阵列和像素之间的固定采样距离。该算法必须检查哪些像素组合不是“空”高度,并相应地调整计算(这就是位掩码所做的)。

共有1个答案

景国兴
2023-03-14

使用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×2窗口的左上角值)。

 类似资料:
  • 问题内容: 我有一个数字列表,它们表示另一个程序产生的矩阵或数组的展平输出,我知道原始数组的尺寸,并想将这些数字读回到列表列表或NumPy矩阵中。原始数组中的尺寸可能超过2个。 例如 将产生: [[0,2,7,6],[3,1,4,5]] 提前加油 问题答案: 用途: 您也可以直接分配给的属性,如果你想避免在内存中复制它:

  • 问题内容: 我正在尝试创建具有混合数据类型(字符串,整数,整数)的NumPy数组/矩阵(Nx3)。但是,当我通过添加一些数据来添加此矩阵时,出现错误: TypeError:无效的类型提升 。拜托,有人可以帮我解决这个问题吗? 当我用示例数据创建一个数组时,NumPy将矩阵中的所有列都转换为一种“ S”数据类型。而且我无法为数组指定数据类型,因为当我执行此操作时, res = np.array([“

  • 这得到了我想要的,但可能没有很好地扩展? 产量

  • 问题内容: 假设您具有以下数组: 您将如何将其转换为XML字符串,使其看起来像: 一种方法是通过类似如下的递归方法: 我正在寻找一种使用迭代的方法。 问题答案:

  • 问题内容: 我有以下内容: 如何在XYZ_2上执行与在XYZ_2上相同的操作?我会以某种方式首先重塑数组吗? 问题答案: 您似乎正在尝试的最后一个轴 与最后一个 。因此,您可以像这样使用- 相关帖子了解。 为了完整起见,在交换的最后两个轴后,我们当然也可以使用,例如- 这将不如一个高效。 运行时测试- 一般而言,涉及张量时,效率要高得多。由于的轴只有一个,因此我们可以通过重整,使用,获取结果并将其

  • 做一些类似的事情 使用多个内核,运行良好。 所以,如果我要做整数矩阵乘法,我得做下面的一个: 使用numpy慢得让人痛苦的并庆幸我可以保留8位整数。 使用Scipy的并使用4倍内存。 使用numpy的并且只使用2倍内存,但要注意的是,在float16数组上的速度要比在float32数组上慢得多,比int8慢得多。 为多线程整数矩阵乘法找到一个优化的库(其实Mathematica就是这么做的,但我更