computer-vision - 检测连接体素的环/电路

Question

我有一个骨架化的体素结构，如下所示:

实际结构比这个例子大很多有没有办法找到结构中的闭环？我尝试将其转换为图形并使用基于图形的方法，但它们都存在图形没有节点位置的空间信息的问题，因此图形可以具有多个同源环。

不可能找到所有环然后过滤掉感兴趣的环，因为图形太大了。环的大小变化很大。

感谢您的帮助和贡献!

尽管我主要使用 Python 和 Matlab 工作，但欢迎使用任何语言方法和伪代码。

---

编辑:

不，图形不是平面的。Graph cycle base 的问题与其他简单的基于图的方法相同。该图缺少任何空间信息，不同的空间配置可以具有相同的循环基础，因此循环基础不一定对应于图中的循环或空洞。

这是稀疏格式的邻接矩阵:

NodeID1 NodeID2 Weight

Pastebin with adjacency matrix

下面是图表节点对应的 X、Y、Z 坐标:

X Y Z

Pastebin with node coordinates

(实际结构比这个例子大很多)

Answer 1

首先，我通过将度数为 2 的相邻节点收缩为超节点来显着减小问题的规模:图中的每个简单链都被单个节点替代。

然后我找到 cycle basis ，其中基组循环的最大成本是最小的。

对于网络的中心部分，解决方案很容易绘制，因为它是平面的:

出于某种原因，我未能正确识别周期基础，但我认为以下内容绝对可以帮助您入门，也许其他人可以插话。

从发布的图像中恢复数据(因为 OP 不会提供一些真实数据)

import numpy as np
import matplotlib.pyplot as plt
from skimage.morphology import medial_axis, binary_closing
from matplotlib.patches import Path, PathPatch
import itertools
import networkx as nx

img = plt.imread("tissue_skeleton_crop.jpg")
# plt.hist(np.mean(img, axis=-1).ravel(), bins=255) # find a good cutoff
bw = np.mean(img, axis=-1) < 200
# plt.imshow(bw, cmap='gray')
closed = binary_closing(bw, selem=np.ones((50,50))) # connect disconnected segments
# plt.imshow(closed, cmap='gray')
skeleton = medial_axis(closed)

fig, ax = plt.subplots(1,1)
ax.imshow(skeleton, cmap='gray')
ax.set_xticks([])
ax.set_yticks([])

def img_to_graph(binary_img, allowed_steps):
    """
    Arguments:
    ----------
    binary_img    -- 2D boolean array marking the position of nodes
    allowed_steps -- list of allowed steps; e.g. [(0, 1), (1, 1)] signifies that
                     from node with position (i, j) nodes at position (i, j+1)
                     and (i+1, j+1) are accessible,

    Returns:
    --------
    g             -- networkx.Graph() instance
    pos_to_idx    -- dict mapping (i, j) position to node idx (for testing if path exists)
    idx_to_pos    -- dict mapping node idx to (i, j) position (for plotting)
    """

    # map array indices to node indices and vice versa
    node_idx = range(np.sum(binary_img))
    node_pos = zip(*np.where(np.rot90(binary_img, 3)))
    pos_to_idx = dict(zip(node_pos, node_idx))

    # create graph
    g = nx.Graph()
    for (i, j) in node_pos:
        for (delta_i, delta_j) in allowed_steps: # try to step in all allowed directions
            if (i+delta_i, j+delta_j) in pos_to_idx: # i.e. target node also exists
                g.add_edge(pos_to_idx[(i,j)], pos_to_idx[(i+delta_i, j+delta_j)])

    idx_to_pos = dict(zip(node_idx, node_pos))

    return g, idx_to_pos, pos_to_idx

allowed_steps = set(itertools.product((-1, 0, 1), repeat=2)) - set([(0,0)])
g, idx_to_pos, pos_to_idx = img_to_graph(skeleton, allowed_steps)

fig, ax = plt.subplots(1,1)
nx.draw(g, pos=idx_to_pos, node_size=1, ax=ax)

注意:这些不是红线，这些是与图中节点对应的许多红点。

合约图

def contract(g):
    """
    Contract chains of neighbouring vertices with degree 2 into one hypernode.

    Arguments:
    ----------
    g -- networkx.Graph or networkx.DiGraph instance

    Returns:
    --------
    h -- networkx.Graph or networkx.DiGraph instance
        the contracted graph

    hypernode_to_nodes -- dict: int hypernode -> [v1, v2, ..., vn]
        dictionary mapping hypernodes to nodes

    """

    # create subgraph of all nodes with degree 2
    is_chain = [node for node, degree in g.degree() if degree == 2]
    chains = g.subgraph(is_chain)

    # contract connected components (which should be chains of variable length) into single node
    components = list(nx.components.connected_component_subgraphs(chains))
    hypernode = g.number_of_nodes()
    hypernodes = []
    hyperedges = []
    hypernode_to_nodes = dict()
    false_alarms = []
    for component in components:
        if component.number_of_nodes() > 1:

            hypernodes.append(hypernode)
            vs = [node for node in component.nodes()]
            hypernode_to_nodes[hypernode] = vs

            # create new edges from the neighbours of the chain ends to the hypernode
            component_edges = [e for e in component.edges()]
            for v, w in [e for e in g.edges(vs) if not ((e in component_edges) or (e[::-1] in component_edges))]:
                if v in component:
                    hyperedges.append([hypernode, w])
                else:
                    hyperedges.append([v, hypernode])

            hypernode += 1

        else: # nothing to collapse as there is only a single node in component:
            false_alarms.extend([node for node in component.nodes()])

    # initialise new graph with all other nodes
    not_chain = [node for node in g.nodes() if not node in is_chain]
    h = g.subgraph(not_chain + false_alarms)
    h.add_nodes_from(hypernodes)
    h.add_edges_from(hyperedges)

    return h, hypernode_to_nodes

h, hypernode_to_nodes = contract(g)

# set position of hypernode to position of centre of chain
for hypernode, nodes in hypernode_to_nodes.items():
    chain = g.subgraph(nodes)
    first, last = [node for node, degree in chain.degree() if degree==1]
    path = nx.shortest_path(chain, first, last)
    centre = path[len(path)/2]
    idx_to_pos[hypernode] = idx_to_pos[centre]

fig, ax = plt.subplots(1,1)
nx.draw(h, pos=idx_to_pos, node_size=20, ax=ax)

寻找周期基础

cycle_basis = nx.cycle_basis(h)

fig, ax = plt.subplots(1,1)
nx.draw(h, pos=idx_to_pos, node_size=10, ax=ax)
for cycle in cycle_basis:
    vertices = [idx_to_pos[idx] for idx in cycle]
    path = Path(vertices)
    ax.add_artist(PathPatch(path, facecolor=np.random.rand(3)))

待办事项:

找到正确的循环基础(我可能会混淆 cycle basis 是什么，或者 networkx 可能有错误)。

编辑

我靠，这是一场绝技。我不应该深入这个兔子洞。

所以现在的想法是我们要找到循环基础，其中循环的最大成本是最小的。我们将循环的成本设置为其边长，但可以想象其他成本函数。为此，我们找到一个初始循环基础，然后我们在基础中组合循环，直到找到具有所需属性的循环集。

def find_holes(graph, cost_function):
    """
    Find the cycle basis, that minimises the maximum individual cost of the cycles in the basis set.
    """

    # get cycle basis
    cycles = nx.cycle_basis(graph)

    # find new basis set that minimises maximum cost
    old_basis = set()
    new_basis = set(frozenset(cycle) for cycle in cycles) # only frozensets are hashable
    while new_basis != old_basis:
        old_basis = new_basis
        for cycle_a, cycle_b in itertools.combinations(old_basis, 2):
            if len(frozenset.union(cycle_a, cycle_b)) >= 2: # maybe should check if they share an edge instead
                cycle_c = _symmetric_difference(graph, cycle_a, cycle_b)
                new_basis = new_basis.union([cycle_c])
        new_basis = _select_cycles(new_basis, cost_function)

    ordered_cycles = [order_nodes_in_cycle(graph, nodes) for nodes in new_basis]
    return ordered_cycles

def _symmetric_difference(graph, cycle_a, cycle_b):
    # get edges
    edges_a = list(graph.subgraph(cycle_a).edges())
    edges_b = list(graph.subgraph(cycle_b).edges())

    # also get reverse edges as graph undirected
    edges_a += [e[::-1] for e in edges_a]
    edges_b += [e[::-1] for e in edges_b]

    # find edges that are in either but not in both
    edges_c = set(edges_a) ^ set(edges_b)

    cycle_c = frozenset(nx.Graph(list(edges_c)).nodes())
    return cycle_c

def _select_cycles(cycles, cost_function):
    """
    Select cover of nodes with cycles that minimises the maximum cost
    associated with all cycles in the cover.
    """
    cycles = list(cycles)
    costs = [cost_function(cycle) for cycle in cycles]
    order = np.argsort(costs)

    nodes = frozenset.union(*cycles)
    covered = set()
    basis = []

    # greedy; start with lowest cost
    for ii in order:
        cycle = cycles[ii]
        if cycle <= covered:
            pass
        else:
            basis.append(cycle)
            covered |= cycle
            if covered == nodes:
                break

    return set(basis)

def _get_cost(cycle, hypernode_to_nodes):
    cost = 0
    for node in cycle:
        if node in hypernode_to_nodes:
            cost += len(hypernode_to_nodes[node])
        else:
            cost += 1
    return cost

def _order_nodes_in_cycle(graph, nodes):
    order, = nx.cycle_basis(graph.subgraph(nodes))
    return order

holes = find_holes(h, cost_function=partial(_get_cost, hypernode_to_nodes=hypernode_to_nodes))

fig, ax = plt.subplots(1,1)
nx.draw(h, pos=idx_to_pos, node_size=10, ax=ax)
for ii, hole in enumerate(holes):
    if (len(hole) > 3):
        vertices = np.array([idx_to_pos[idx] for idx in hole])
        path = Path(vertices)
        ax.add_artist(PathPatch(path, facecolor=np.random.rand(3)))
        xmin, ymin = np.min(vertices, axis=0)
        xmax, ymax = np.max(vertices, axis=0)
        x = xmin + (xmax-xmin) / 2.
        y = ymin + (ymax-ymin) / 2.
        # ax.text(x, y, str(ii))