numpy 从图像中构造邻接矩阵的优化方法

pkmbmrz7  于 12个月前  发布在  其他
关注(0)|答案(1)|浏览(171)

我试图从高程栅格的像素构建邻接矩阵。栅格是一个GeoTIFF图像,在分水岭外具有指定的nodata值。
我现在的解决方案是有效的,但远远不是最佳的。
我所尝试的:

  • Scikit-learn:sklearn.feature_extraction.image.img_to_graph工作正常,但它只能得到4路邻居(车邻居),我需要它8路(皇后邻居)。
  • Scipy:我现在有了解决方案。在一个numpy 2d数组上搜索,在一个3x 3窗口上搜索,将中心与它的邻居进行比较,搜索非nodata像素,并馈送一个scipy稀疏矩阵。

问题是,我可以运行2000 x2000的图像,但我需要测试它与一个更大的图像(超过25000像素的每一方)。
有什么方法可以优化这个算法之前,我可以尝试它在一个更强大的计算机?
先谢了。
这是我现在的代码:

def image_to_adjacency_matrix(image_array, nodata = None):
    image_flat = image_array.flatten()

    height, width = image_array.shape
    adjacency_matrix = scipy.sparse.lil_matrix((height * width, height * width), dtype=int)

    for i in range(height):
        for j in range(width):
            current_pixel = i * width + j

            if image_flat[current_pixel] == nodata:
                continue

            for ni in range(i - 1, i + 2):
                for nj in range(j - 1, j + 2):
                    if 0 <= ni < height and 0 <= nj < width:
                        neighbor_pixel = ni * width + nj

                        if image_flat[neighbor_pixel] == nodata:
                            continue

                        if current_pixel != neighbor_pixel:
                            adjacency_matrix[current_pixel, neighbor_pixel] = 1

    return adjacency_matrix.tocsr()

字符串

q5iwbnjs

q5iwbnjs1#

尝试使用NumPy将其矢量化。
注意事项:我实现了你的问题的简化版本,它通过只创建1到width - 1的边而不是从0到width的边来避免环绕光栅的边缘。这简化了逻辑,足以让我用NumPy索引解决问题。

def image_to_adjacency_matrix_opt(image_array, nodata = None):
    image_flat = image_array.flatten()

    height, width = image_array.shape
    N = height * width
    image_has_data = image_flat != nodata
    index_dtype = np.int32 if N < 2 ** 31 else np.int64
    adjacents = np.array([
        -width - 1, -width, -width + 1,
        -1,                 1,
        width - 1,  width,  width + 1 
    ], dtype=index_dtype)
    row_idx, col_idx = np.meshgrid(
        np.arange(1, height - 1, dtype=index_dtype),
        np.arange(1, width - 1, dtype=index_dtype),
        indexing='ij'
    )
    row_idx = row_idx.reshape(-1)
    col_idx = col_idx.reshape(-1)
    pixel_idx = row_idx * width + col_idx
    pixels_with_data = image_has_data[pixel_idx]
    pixel_idx = pixel_idx[pixels_with_data]
    neighbors = pixel_idx.reshape(-1, 1) + adjacents.reshape(1, -1)
    neighbors_with_data = image_has_data[neighbors]
    row = np.repeat(pixel_idx, repeats=neighbors_with_data.sum(axis=1))
    col = neighbors[neighbors_with_data]
    data = np.ones(len(row), dtype='uint8')
    adjacency_matrix = scipy.sparse.coo_matrix((data, (row, col)), dtype=int, shape=(N, N))
    return adjacency_matrix.tocsr()

字符串
我发现对于一个500 x500的矩阵来说,这大约快了100倍,其中50%的条目随机设置为nodata。
测试数据集示例:

N = 500
image = np.random.randint(0, 100, size=(N, N))
image[np.random.rand(N, N) < 0.5] = -1
image = image.astype('int8')

相关问题