在骨架化的运河结构中寻找死角

Finding dead ends in a skeletonized canal structure

我对编码还很陌生,这是我正在处理的第一个真正的代码。我正在尝试编写代码,它在骨架化的 3d 结构中找到死角并设置标记和死角的位置。所以我只有值为 0 和 1 的像素。我现在的问题是,当运河分裂时,我只能跟踪 1 个分支,而第二个分支没有迭代,所以我错过了很多死胡同。 这是我用自建运河写的代码:


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()

我真的不知道我怎么能在分裂后抓住每条可能的运河并以这种方式找到每个死者。我的代码现在只采用它找到的第一条运河,然后遍历那条运河。 我愿意接受有关特定问题或一般代码的任何建议。 提前致谢!

欢迎来到 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.

我们可以将数据表示为网络(而不是 3D 图像)。在网络中,每个 non-zero 像素对应一个节点,如果像素在 3D 数组中相邻,则两个节点相连。在这个网络中,“死胡同”很容易识别,因为这些节点最多只连接到一个其他节点,因此度数为 1(so-called 图论术语中的叶节点)。

#!/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()