检测连接体素的环/电路 [英] Detect rings/circuits of connected voxels

问题描述

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

实际结构明显比这个例子大.有什么方法可以找到结构中的闭环吗?我尝试将其转换为图并使用基于图的方法，但它们都有一个问题，即图没有节点位置的空间信息，因此图可以有多个同源环.

不可能找到所有环然后过滤掉感兴趣的环，因为图表太大了.环的大小差异很大.

感谢您的帮助和贡献！

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

<小时>

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

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

NodeID1 NodeID2 权重

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

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

将 numpy 导入为 np将 matplotlib.pyplot 导入为 plt从 skimage.morphology 导入 medial_axis，binary_closure从 matplotlib.patches 导入路径，PathPatch导入迭代工具将 networkx 导入为 nximg = plt.imread("tissue_skeleton_crop.jpg")# plt.hist(np.mean(img, axis=-1).ravel(), bins=255) # 找到一个好的截止点bw = np.mean(img, axis=-1) <200# plt.imshow(bw, cmap='灰色')closed = binary_closing(bw, selem=np.ones((50,50))) # 连接断开的段# plt.imshow(关闭，cmap='灰色')骨架 = 内侧轴(关闭)无花果，斧头 = plt.subplots(1,1)ax.imshow(骨架，cmap='灰色')ax.set_xticks([])ax.set_yticks([])

def img_to_graph(binary_img, allowed_steps):"""论据:----------binary_img -- 标记节点位置的二维布尔数组allowed_steps -- 允许的步骤列表；例如[(0, 1), (1, 1)] 表示来自位置为 (i, j) 的节点 位置为 (i, j+1) 的节点和 (i+1, j+1) 是可访问的，回报:--------g -- networkx.Graph() 实例pos_to_idx -- 将 (i, j) 位置映射到节点 idx(用于测试路径是否存在)idx_to_pos -- dict 将节点 idx 映射到 (i, j) 位置(用于绘图)"""# 将数组索引映射到节点索引，反之亦然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))# 创建图表g = nx.Graph()对于 node_pos 中的 (i, j):for (delta_i, delta_j) in allowed_steps: # 尝试在所有允许的方向上步进if (i+delta_i, j+delta_j) in pos_to_idx: # 即目标节点也存在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))返回 g，idx_to_pos，pos_to_idxallowed_steps = set(itertools.product((-1, 0, 1), repeat=2)) - set([(0,0)])g，idx_to_pos，pos_to_idx = img_to_graph(骨架，allowed_steps)无花果，斧头 = plt.subplots(1,1)nx.draw(g, pos=idx_to_pos, node_size=1, ax=ax)

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

合约图

def 合约(g):"""将度数为 2 的相邻顶点的链收缩到一个超节点中.论据:----------g -- networkx.Graph 或 networkx.DiGraph 实例回报:--------h -- networkx.Graph 或 networkx.DiGraph 实例收缩图hypernode_to_nodes -- dict: int hypernode ->[v1, v2, ..., vn]将超节点映射到节点的字典"""# 创建所有度数为 2 的节点的子图is_chain = [node for node, degree in g.degree() if degree == 2]链 = g.subgraph(is_chain)# 将连接的组件(应该是可变长度的链)收缩到单个节点中组件=列表(nx.components.connected_component_subgraphs(链))超节点 = g.number_of_nodes()超节点 = []超边 = []hypernode_to_nodes = dict()假警报 = []对于组件中的组件:如果 component.number_of_nodes() >1:hypernodes.append(超节点)vs = [component.nodes() 中节点的节点]hypernode_to_nodes [hypernode] = vs# 创建从链端的邻居到超节点的新边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))]:如果 v 在组件中:hyperedges.append([超节点，w])别的:hyperedges.append([v, 超节点])超节点 += 1else: # 没有什么可折叠的，因为组件中只有一个节点:false_alarms.extend([component.nodes() 中节点的节点])# 用所有其他节点初始化新图not_chain = [如果不是 is_chain 中的节点，则 g.nodes() 中的节点的节点]h = g.subgraph(not_chain + false_alarms)h.add_nodes_from(超节点)h.add_edges_from(超边)返回 h，hypernode_to_nodesh, hypernode_to_nodes = contract(g)# 将超节点的位置设置为链中心的位置对于超节点，hypernode_to_nodes.items() 中的节点:链 = g.subgraph(节点)first, last = [node for node, degree in chain.degree() if degree==1]路径 = nx.shortest_path(链，第一，最后)中心 = 路径[len(路径)/2]idx_to_pos[超节点] = idx_to_pos[中心]无花果，斧头 = plt.subplots(1,1)nx.draw(h, pos=idx_to_pos, node_size=20, ax=ax)

寻找循环基础

cycle_basis = nx.cycle_basis(h)无花果，斧头 = plt.subplots(1,1)nx.draw(h, pos=idx_to_pos, node_size=10, ax=ax)对于cycle_basis中的循环:vertices = [idx_to_pos[idx] for idx in cycle]路径 = 路径(顶点)ax.add_artist(PathPatch(path, facecolor=np.random.rand(3)))

待办事项:

找到正确的循环基础(我可能会混淆

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

def find_holes(graph, cost_function):"""找到循环基，它使基组中循环的最大单个成本最小化."""# 获取循环基础周期 = nx.cycle_basis(graph)# 找到最小化最大成本的新基组old_basis = set()new_basis = set(frozenset(cycle) for cycle in cycles) # 只有frozensets是可散列的而 new_basis != old_basis:旧基础 = 新基础对于 itertools.combinations(old_basis, 2) 中的 cycle_a、cycle_b:if len(frozenset.union(cycle_a, cycle_b)) >= 2: # 也许应该检查它们是否共享一条边cycle_c = _symmetric_difference(图，cycle_a，cycle_b)new_basis = new_basis.union([cycle_c])new_basis = _select_cycles(new_basis, cost_function)ordered_cycles = [new_basis 中节点的 order_nodes_in_cycle(graph, nodes)]返回有序循环def _symmetric_difference(图，cycle_a，cycle_b):# 获取边edge_a = list(graph.subgraph(cycle_a).edges())edge_b = list(graph.subgraph(cycle_b).edges())# 也得到反向边作为无向图edge_a += [e[::-1] for e in edges_a]edge_b += [e[::-1] for e in edges_b]# 找出在其中一个中但不在两者中的边边缘_c = 集合(边缘_a) ^ 集合(边缘_b)cycle_c = freezeset(nx.Graph(list(edges_c)).nodes())返回周期_cdef _select_cycles(周期，成本函数):"""选择具有最小化最大成本的循环的节点覆盖与封面中的所有循环相关联."""周期=列表(周期)成本 = [周期成本函数(周期)]order = np.argsort(costs)节点 = freezeset.union(*cycles)覆盖=设置()基础= []# 贪婪的;从最低成本开始对于 ii 按顺序:周期 = 周期[ii]如果循环 <= 覆盖:经过别的:基础.追加(循环)覆盖 |= 循环如果覆盖 == 节点:休息返回集(基础)def _get_cost(周期，hypernode_to_nodes):成本 = 0对于循环中的节点:如果 hypernode_to_nodes 中的节点:成本 += len(hypernode_to_nodes[node])别的:成本 += 1退货成本def _order_nodes_in_cycle(图，节点):顺序，= nx.cycle_basis(graph.subgraph(nodes))退货单孔= find_holes(h，cost_function=partial(_get_cost，hypernode_to_nodes=hypernode_to_nodes))无花果，斧头 = plt.subplots(1,1)nx.draw(h, pos=idx_to_pos, node_size=10, ax=ax)对于 ii，枚举中的孔(孔):如果 (len(hole) > 3):vertices = np.array([idx_to_pos[idx] for idx in hole])路径 = 路径(顶点)ax.add_artist(PathPatch(path, facecolor=np.random.rand(3)))xmin, ymin = np.min(顶点, 轴=0)xmax, ymax = np.max(顶点, 轴=0)x = xmin + (xmax-xmin)/2.y = ymin + (ymax-ymin)/2.# ax.text(x, y, str(ii))

I have a skeletonized voxel structure that looks like this:

The actual structure is significantly larger than this exampleIs there any way to find the closed rings in the structure? I tried converting it to a graph and using graph based approaches but they all have the problem that a graph has no spatial information of node position and hence a graph can have multiple rings that are homologous.

It is not possible to find all the rings and then filter out the ones of interest since the graph is just too large. The size of the rings varies significantly.

Thanks for your help and contribution!

Any language approaches and pseudo-code are welcomed though I work mostly in Python and Matlab.

---

EDIT:

No the graph is not planar. The problem with the Graph cycle base is the same as with other simple graph based approaches. The graph lacks any spatial information and different spatial configurations can have the same cycle base, hence the cycle base does not necessarily correspond to the cycles or holes in the graph.

Here is the adjacency matrix in sparse format:

NodeID1 NodeID2 Weight

Pastebin with adjacency matrix

And here are the corresponding X,Y,Z coordinates for the Nodes of the graph:

X Y Z

Pastebin with node coordinates

(The actual structure is significantly larger than this example)

解决方案

First I reduce the size of the problem considerably by contracting neighbouring nodes of degree 2 into hypernodes: each simple chain in the graph is substituted with a single node.

Then I find the cycle basis, for which the maximum cost of the cycles in the basis set is minimal.

For the central part of the network, the solution can easily be plotted as it is planar:

For some reason, I fail to correctly identify the cycle basis but I think the following should definitely get you started and maybe somebody else can chime in.

Recover data from posted image (as OP wouldn't provide some real data)

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)

NB: These are not red lines, these are lots of red dots corresponding to nodes in the graph.

Contract Graph

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)

Find cycle basis

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

TODO:

Find the correct cycle basis (I might be confused what the cycle basis is or networkx might have a bug).

EDIT

Holy crap, this was a tour-de-force. I should have never delved into this rabbit hole.

So the idea is now that we want to find the cycle basis for which the maximum cost for the cycles in the basis is minimal. We set the cost of a cycle to its length in edges, but one could imagine other cost functions. To do so, we find an initial cycle basis, and then we combine cycles in the basis until we find the set of cycles with the desired property.

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

这篇关于检测连接体素的环/电路的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！