Finding dead ends in a skeletonized canal structure

Question

I am very new to coding and this is the first real code I am working on. I am trying to write code, which finds dead ends in a skeletonized 3d structure and sets a marker and the position of the dead end. So I only have pixels with the value 0 and 1. My problem right now is when a canal splits I can only track 1 branch and the second branch is not iterated through so that I miss a lot of dead ends. This is the code I wrote with self constructed canals:


c = np.zeros([40,40,30], dtype=np.unit8)

c[5,5,0:4] = 1
c[5,4,4] = 1
c[5,3,4:10] = 1
c[4,4,10:25] = 1
c[5,5,25:30] = 1
c[5,6,4] = 1
c[5,7,4:6] = 1
c[6,7,6:15] = 1
c[6,6,15:24] = 1
c[5,5,24:25] = 1

c[15,15,0:15] = 1
c[16,15,15:16] = 1
c[17,16,16:18] = 1
c[16,16,18:20] = 1
c[15,17,20:25] = 1
c[18,17,18:20] = 1
c[19,17,20:30] = 1
c[14,14,15:16] = 1
c[13,13,16:20] = 1
c[13,12,20:21] = 1
c[14,12,21:23] = 1
c[15,11,23:25] = 1
c[12,13,20:21] = 1
c[11,13,21:25] = 1

c[30,30,20:30] = 1
c[29,29,18:20] = 1
c[28,28,15:18] = 1
c[27,27,14:15] = 1
c[26,26,13:14] =1
c[29,29,14:15] = 1
c[30,30,10:14] = 1
c[29,29,9:10] = 1
c[28,29,8:9] = 1
c[27,28,3:8] = 1
c[31,31,9:10] = 1
c[32,32,6:9] = 1

k = np.pad(c,pad_width=2)
b = np.array(k,dtype=np.uint8)

d = np.size(b,axis=0)
h = np.size(b,axis=1)
w = np.size(b,axis=2)

img = np.zeros_like(b,dtype=np.uint8)

occupied_pixels = np.where(b[:,:,3] == 1)
pixels_x = occupied_pixels[0]
pixels_y = occupied_pixels[1]


def is_deadend(x, y, z):
    pixel_found_in_cone = False
    for dz in range(1,3):
        for dx in [2, 1, 0, -1, -2]:   #for dx i2,n [-1, 0, 1]:
            for dy in [2, 1, 0, -1, -2]:   #for dy in [-1, 0, 1]
                if b[x+dx,y+dy,z+dz] == 1:
                    pixel_found_in_cone = True
                    position, bool_deadend = is_deadend(x+dx,y+dy,z+dz) 
                    if bool_deadend:
                        return position, bool_deadend   
    return [x,y,z], not pixel_found_in_cone 

markers_found = []

for i in range(len(pixels_x)): 
    occupied_pixel = [pixels_x[i],pixels_y[i],3] 
    position, bool_deadend = is_deadend(occupied_pixel[0],occupied_pixel[1],occupied_pixel[2])
    print(is_deadend(occupied_pixel[0],occupied_pixel[1],occupied_pixel[2]))
    if bool_deadend:
        markers_found.append(position)

for marker in markers_found:
    img[marker[0],marker[1],marker[2]] = 1

#for visualizing:
    
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
data = b
x, y, z = data.nonzero()
ax.scatter(x,y,z, c=z, alpha =1)
plt.show()


fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
data = img
x, y, z = data.nonzero()
ax.scatter(x,y,z, c=z, alpha =1)
plt.show()

I don´t really know how I can catch every possible canal after a splitup and find every dead in that way. My code right now only takes the first canal it finds and then iterates through that one. I am open for any suggestion regarding the specific problem or the code in general. Thanks in advance!

user_2912349 · Answer

Welcome to SO!

I don´t really know how I can catch every possible canal after a splitup and find every dead in that way. My code right now only takes the first canal it finds and then iterates through that one.

We can represent the data as a network (instead of a 3D image). In the network, each non-zero pixel corresponds to a node, and two nodes are connected if the pixels are next to each other in the 3D array. In this network, "dead ends" are easily identified by the fact that these nodes only connect to at most one other node, and hence have a degree of 1 (so-called leaf nodes in graph theory parlance).

#!/usr/bin/env python
import itertools
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx


def data_to_graph(data, allowed_steps):

    nodes = list(range(np.sum(data > 0)))
    positions = list(zip(*np.where(data)))
    node_positions = dict(zip(nodes, positions))
    inverse_mapping = {position : node for node, position in node_positions.items()}

    g = nx.Graph()
    for node, position in node_positions.items():
        for delta in allowed_steps:
            if tuple(np.array(position) + np.array(delta)) in inverse_mapping:
                g.add_edge(node, inverse_mapping[tuple(np.array(position) + np.array(delta))])

    return g, node_positions


if __name__ == '__main__':

    data = np.zeros([40,40,30], dtype=np.uint8)

    data[5,5,0:4] = 1
    data[5,4,4] = 1
    data[5,3,4:10] = 1
    data[4,4,10:25] = 1
    data[5,5,25:30] = 1
    data[5,6,4] = 1
    data[5,7,4:6] = 1
    data[6,7,6:15] = 1
    data[6,6,15:24] = 1
    data[5,5,24:25] = 1

    data[15,15,0:15] = 1
    data[16,15,15:16] = 1
    data[17,16,16:18] = 1
    data[16,16,18:20] = 1
    data[15,17,20:25] = 1
    data[18,17,18:20] = 1
    data[19,17,20:30] = 1
    data[14,14,15:16] = 1
    data[13,13,16:20] = 1
    data[13,12,20:21] = 1
    data[14,12,21:23] = 1
    data[15,11,23:25] = 1
    data[12,13,20:21] = 1
    data[11,13,21:25] = 1

    data[30,30,20:30] = 1
    data[29,29,18:20] = 1
    data[28,28,15:18] = 1
    data[27,27,14:15] = 1
    data[26,26,13:14] =1
    data[29,29,14:15] = 1
    data[30,30,10:14] = 1
    data[29,29,9:10] = 1
    data[28,29,8:9] = 1
    data[27,28,3:8] = 1
    data[31,31,9:10] = 1
    data[32,32,6:9] = 1

    # identify dead ends
    dx = (-1, 0, 1)
    dy = (-1, 0, 1)
    dz = (-1, 0, 1)
    allowed_steps = set(itertools.product(dx, dy, dz)) - set([(0, 0, 0)])
    g, node_positions = data_to_graph(data, allowed_steps)
    dead_ends = [node for node, degree in dict(g.degree).items() if degree == 1]
    dead_end_positions = [node_positions[node] for node in dead_ends]
    print(f"Pixels that are a dead end: {dead_end_positions}")

    # visualise
    fig = plt.figure(figsize=(8,8))
    ax = fig.add_subplot(111, projection='3d')
    x, y, z = np.array(list(node_positions.values())).transpose()
    node_color = ['tab:red' if node in dead_ends else 'tab:blue' for node in node_positions]
    ax.scatter(x, y, z, c=node_color, s=20)
    for (source, target) in g.edges:
        ax.plot(*zip(node_positions[source], node_positions[target]), color='gray', linewidth=1)

    plt.show()

Finding dead ends in a skeletonized canal structure

1 Answers