从图像构建邻接矩阵的最佳方法

Question

我正在尝试从高程栅格的像素构建邻接矩阵。栅格是一个 GeoTIFF 图像，在分水岭之外具有指定的节点数据值。

我现在的解决方案是有效的，但还远非最佳。

我尝试过的：

• Scikit-learn：

  sklearn.feature_extraction.image.img_to_graph

  工作正常，但它只能获得 4 路邻域（rook 邻域），而我需要它 8 路邻域（queen 邻域）。
• Scipy：我现在的解决方案。要迭代 numpy 2d 数组并迭代 3x3 窗口，将中心与其邻居进行比较，搜索非节点数据像素并输入 scipy 稀疏矩阵。

问题是，我可以运行 2000x2000 图像，但我需要用更大的图像（每边超过 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()

Answer 1

尝试使用 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
    image_has_data = image_flat != nodata
    adjacents = np.array([
        -width - 1, -width, -width + 1,
        -1,                 1,
        width - 1,  width,  width + 1 
    ])
    row_idx, col_idx = np.meshgrid(np.arange(1, height - 1), np.arange(1, width - 1), 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=(height * width, height * width))
    return adjacency_matrix.tocsr()

我发现对于 500x500 矩阵（其中 50% 的条目随机设置为 nodata）来说，速度大约快 100 倍。

测试数据集示例：

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