python - 从栅格图像中提取街道网络

Question

我有一张街道网格的 512x512 图像:

我想为此图像中的每条街道提取折线(大蓝点 = 交叉点，小蓝点 = 沿折线的点):

我尝试了一些技巧!一个想法是从 skeletonize 开始将街道压缩到 1px 宽的线:

from skimage import morphology
morphology.skeletonize(streets_data))

不幸的是，这有一些差距会破坏街道网络的连通性；我不完全确定为什么，但我猜这是因为有些街道在某些地方窄 1px，而在其他地方宽 1px。 (更新:间隙不是真实的；它们完全是我展示骨架的人工制品。有关悲伤的故事，请参阅 this comment。骨架连接良好。)
我可以使用 binary_dilation 修补这些, 代价是让街道的宽度再次发生变化:

out = morphology.skeletonize(streets_data)
out = morphology.binary_dilation(out, morphology.selem.disk(1))

使用重新连接的网格，我可以运行 Hough 变换来查找线段:

import cv2
rho = 1  # distance resolution in pixels of the Hough grid
theta = np.pi / 180  # angular resolution in radians of the Hough grid
threshold = 8  # minimum number of votes (intersections in Hough grid cell)
min_line_length = 10  # minimum number of pixels making up a line
max_line_gap = 2  # maximum gap in pixels between connectable line segments

# Run Hough on edge detected image
# Output "lines" is an array containing endpoints of detected line segments
lines = cv2.HoughLinesP(
    out, rho, theta, threshold, np.array([]),
    min_line_length, max_line_gap
)

line_image = streets_data.copy()
for i, line in enumerate(lines):
    for x1,y1,x2,y2 in line:
        cv2.line(line_image,(x1,y1),(x2,y2), 2, 1)

这会产生一整个困惑的重叠线段，以及一些间隙(看看右侧的 T 交叉点):

在这一点上，我可以尝试去重复重叠的线段，但我不太清楚这是通往解决方案的途径，特别是考虑到这个差距。
是否有更直接的方法可用于获取我正在寻找的折线网络？特别是，有哪些方法可以用于:

寻找交叉点(四向交叉和 T 交叉)。

将街道缩小到 1px 宽，允许可能有一些可变宽度。

查找交叉点之间的折线。

Answer 1

如果你想改善你的“骨架化”，你可以尝试以下算法来获得“1像素宽的街道”:

import imageio
import numpy as np
from matplotlib import pyplot as plt
from scipy.ndimage import distance_transform_edt
from skimage.segmentation import watershed

# read image
image_rgb = imageio.imread('1mYBD.png')

# convert to binary
image_bin = np.max(image_rgb, axis=2) > 0

# compute the distance transform (only > 0)
distance = distance_transform_edt(image_bin)

# segment the image into "cells" (i.e. the reciprocal of the network)
cells = watershed(distance)

# compute the image gradients
grad_v = np.pad(cells[1:, :] - cells[:-1, :], ((0, 1), (0, 0)))
grad_h = np.pad(cells[:, 1:] - cells[:, :-1], ((0, 0), (0, 1)))

# given that the cells have a constant value,
# only the edges will have non-zero gradient
edges = (abs(grad_v) > 0) + (abs(grad_h) > 0)

# extract points into (x, y) coordinate pairs
pos_v, pos_h = np.nonzero(edges)

# display points on top of image
plt.imshow(image_bin, cmap='gray_r')
plt.scatter(pos_h, pos_v, 1, np.arange(pos_h.size), cmap='Spectral')

该算法适用于“街区”而不是“街道”，请查看 cells图片: