Dahu pseudo-distance

This example shows some applications of the Dahu pseudo-distance [4] for grayscale images. These examples are derived from [5].

import matplotlib.pyplot as plt
import numpy as np
from skimage.graph import route_through_array
from numba import njit
from typing import Optional
from collections.abc import Sequence
from imageio.v3 import imread

import pylena as pln

noted that the ToS nodemap is returned in the Khalimsky domain, with a domain that is twice the size of the original one. Furthermore, an artificial border has been added to ensure that the root of the ToS is a shape that surround all the image.

img = imread("../../source/_static/images/eye.tif")
img = np.mean(img, axis=2).astype(np.uint8)
t = pln.morpho.tos(img, padding="median", subsampling="full")

Shortest path finding based on the Dahu pseudo-distance

First we define a function to compute the mask of the connected component on the path between two nodes

@njit
def _dahu_tos_mask_indices(parent: np.ndarray, n1: int, n2: int, depth: np.ndarray) -> np.ndarray:
    M = np.zeros_like(parent)

    while depth[n1] > depth[n2]:
        M[n1] = 1
        n1 = parent[n1]
    while depth[n2] > depth[n1]:
        M[n2] = 1
        n2 = parent[n2]
    while n1 != n2:
        M[n1] = 1
        M[n2] = 1
        n1 = parent[n1]
        n2 = parent[n2]
    M[n1] = 1

    return M


def dahu_tos_mask(t: pln.morpho.ComponentTree, n1: int, n2: int, depth: Optional[np.ndarray] = None) -> np.ndarray:
    if depth is None:
        depth = t.compute_depth()
    mask_nodes = _dahu_tos_mask_indices(t.parent, n1, n2, depth).astype(bool)
    return t.reconstruct(mask_nodes).astype(bool)

Then, we define a function that compute the shortest path in a 2D image between two points based on the Dahu pseudo-distance

def dahu_shortest_path(
    t: pln.morpho.ComponentTree,
    p1: Sequence[int],
    p2: Sequence[int],
    depth: Optional[np.ndarray] = None,
) -> np.ndarray:
    assert len(p1) == len(p2) == t.nodemap.ndim
    if depth is None:
        depth = t.compute_depth()
    p1 = tuple([p1[i] * 2 - 1 for i in range(len(p1))])
    p2 = tuple([p2[i] * 2 - 1 for i in range(len(p2))])

    n1 = t.nodemap[tuple(p1)]
    n2 = t.nodemap[tuple(p2)]
    ROI = dahu_tos_mask(t, n1, n2, depth)
    w = np.ones_like(t.nodemap, dtype=np.uint8) * 255
    w[ROI] = 1
    path, _ = route_through_array(w, p1, p2, geometric=False, fully_connected=False)

    return np.asarray(path)

We define the source and destination points in the image.

p1 = (312, 461)
p2 = (456, 287)

We compute the Dahu shortest path between p1 and p2

path = dahu_shortest_path(t, p1, p2)

We create a color image on which we display the resulting shortest path

mask1 = np.zeros((t.nodemap.shape[0], t.nodemap.shape[1], 3), dtype=np.uint8)
mask1[:, :, 0] = t.reconstruct()
mask1[:, :, 1] = mask1[:, :, 0]
mask1[:, :, 2] = mask1[:, :, 0]
mask1[path[:, 0], path[:, 1], :] = [255, 0, 0]

plt.figure(figsize=(10, 10))
plt.imshow(mask1)
plt.show()