Finding a path between Point in a RGB image

Question

I have a RGB image which is a map including roads and stations.

I am trying to find the closest station to each other that there is a road between those. The road is showing in black and the stations are Red. I added the name so it would be more understandable.

O -> M
A -> B
B -> C
E -> F
G -> H
G -> C
M -> N
.
. 
.

The distance between the two stations is calculated base on the length of the road between those.

Some ideas about to do solve this: I was thinking removing the station from the roads and we have disconnected roads, then use the contour to calculate the length of the road in pixels. But in some cases (like F or E) the station doesn't cover the road completely and do not discount the road to three different part.

Please let me know what you think, how you solve it.

By removing station roads would be like this:

Here is the main image without the labels.

And if it is possible, trace the path between two stations and create the path image. For example From Station A to C:

Answer 1

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