Whew, this was fun to solve! :-)
I used my idea of dilating the red station blobs to make sure they touch the roads they're related to, and @Micka's idea of a "wavefront" solver to figure out distances along the routes to form distance maps.
Following computation of those distance maps, looking up the distance from A to B is simply a matter of reading the value at A's position in B's map (or vice versa).
Code
The code is a bit on the long side, admittedly, but should be commented so you can figure out what's going on. :)
You can find the code which has additional Matplotlib stuff to generate diagnostics and images over here at GitHub.
from itertools import count, combinations
import cv2
import numpy as np
def inrange_thresh(image, color, thresh, binarize_thresh=None):
"""
Apply cv.inRange with a threshold near the given color, optionally threshold the final image.
"""
min_color = tuple(c - thresh for c in color)
max_color = tuple(c + thresh for c in color)
image = cv2.inRange(image, min_color, max_color)
if binarize_thresh is not None:
t, image = cv2.threshold(image, binarize_thresh, 255, cv2.THRESH_BINARY)
return image
def find_feature_bboxes(image):
"""
Find contours in the image and return their bounding boxes.
:return: Iterable of (x, y, w, h)
"""
cnts, *_ = cv2.findContours(image, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
for c in cnts:
yield cv2.boundingRect(c)
def distance_from_point(input_mask, start_point):
"""
Build a distance map following truthy paths in the input mask image, starting from start_point.
:return: Tuple of distance map matrix and infinity value for the matrix
"""
binary_mask = (input_mask > 127)
# Figure out a suitably large number to serve as "infinity" for the mask.
infinite_distance = max(binary_mask.shape) * 2
# Generate a distance map with unreachable points, then seed it with our start point.
dist_map = np.full_like(input_mask, infinite_distance, dtype="uint32")
dist_map[start_point[::-1]] = 0
# Precompute a structuring element we can use to dilate the "wavefront" to walk along the route with.
struct_el = cv2.getStructuringElement(cv2.MORPH_RECT, (3, 3))
for step in count(1):
# Compute a zero map for new neighboring pixels.
neighbor_map = np.full_like(dist_map, 0, dtype="uint8")
# Mask in all of the pixels that were filled in by the last step.
neighbor_map[dist_map == (step - 1)] = 255
# Dilate the map with the structuring element so new neighboring pixels would be included.
neighbor_map = cv2.dilate(neighbor_map, struct_el)
# Figure out which pixels in the dist map should be filled
new_dist_mask = (
(dist_map > step) & # must be more distant than we've already filled
(neighbor_map > 0) & # must be one of these new neighbors
binary_mask # must be walkable
)
if not np.any(new_dist_mask):
# If there are no new pixels, we're done.
break
dist_map[new_dist_mask] = step
return (dist_map, infinite_distance)
def main():
image = cv2.imread("hHwyu.png", cv2.IMREAD_COLOR)
marker_color = (239, 30, 40)[::-1] # RGB -> BGR
route_color = (0, 0, 0)[::-1]
# Grab a grayscale image of the markers only
markers = (inrange_thresh(image, marker_color, 5, 5) > 0).astype(np.uint8)
# Use the center of their bounding boxes as a station location
station_positions = [(int(x + w / 2), int(y + h / 2)) for x, y, w, h in find_feature_bboxes(markers)]
station_positions.sort(key=lambda pair: pair[1])
# Dilate the markers a bit so they'll always sit on the roads, then splat them on.
# We'll use this as the base map for the contour-walking algorithm so it's all connected.
markers_dilated = cv2.dilate(markers, cv2.getStructuringElement(cv2.MORPH_RECT, (5, 5)))
routes_binary = inrange_thresh(image, route_color, 25, 0)
routes_binary[markers_dilated > 0] = 255
station_infos = []
for station_id, start_point in enumerate(station_positions, 1):
print(f"Computing distance map for station {station_id} at {start_point}")
distmap, inf = distance_from_point(routes_binary, start_point)
station_infos.append((station_id, start_point, distmap, inf))
for (sa_id, sa_point, sa_map, sa_inf), (sb_id, sb_point, sb_map, sb_inf) in combinations(station_infos, 2):
distance = sa_map[sb_point[::-1]]
if distance >= sa_inf:
distance = np.inf
print(f"Distance between {sa_id} ({sa_point}) and {sb_id} ({sb_point}): {distance}")
if __name__ == '__main__':
main()
Output
The output (formatted as a distance matrix) is
B / A | #1 #2 #3 #4 #5 #6 #7 #8 #9 #10 #11
#1 | - 356 288 370 212 inf 574 304 inf 455 inf
#2 | 356 - 68 232 495 inf 436 587 inf 317 inf
#3 | 288 68 - 164 427 inf 368 519 inf 249 inf
#4 | 370 232 164 - 509 inf 379 601 inf 260 inf
#5 | 212 495 427 509 - inf 713 176 inf 594 inf
#6 | inf inf inf inf inf - inf inf inf inf inf
#7 | 574 436 368 379 713 inf - 805 inf 212 inf
#8 | 304 587 519 601 176 inf 805 - inf 686 inf
#9 | inf inf inf inf inf inf inf inf - inf 114
#10 | 455 317 249 260 594 inf 212 686 inf - inf
#11 | inf inf inf inf inf inf inf inf 114 inf -
given distance maps and stations as in the diagnostic image below.
Based on the fact that 6, 9 and 11 are disconnected from the rest of the network (and 6 is disconnected from everything!) and those tend to get inf distances in the output, I'd say this works. :-)
Illustrations
Distance maps for each station

Pretty animations
These sped-up gifs illustrate how the distance map and neighbor map "wavefront" work for the upper-left station
Distance map construction

Neighbor map construction
