多维numpy阵列中平方子矩阵的迭代
我已经从图像中加载了数字高程图(浮点高度图)，我正在对数组中的每个2x2平方子矩阵进行迭代，并执行计算并对结果进行求和。
这个操作非常慢，因为我正在使用的高程映射可能非常大(16Kx16K)，而纯Python循环方法要比numpy或can慢得多(至少我读过)。但是，我找不到任何关于如何迭代多维numpy数组块的具体信息。
例如，如果我有以下3x3numpy数组(请记住，这可能是一个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
Stack Overflow用户
发布于 2019-03-05 19:15:12
使用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窗口的左上角值)。
https://stackoverflow.com/questions/55009855
复制相似问题
腾讯云开发者
Copyright © 2013 - 2026 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有
深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 
粤公网安备44030502008569号
腾讯云计算（北京）有限责任公司 京ICP证150476号 | 京ICP备11018762号