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

查看:267
本文介绍了检测连接的体素的环/电路的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我有一个骨架化的体素结构,看起来像这样:



实际结构比该示例大得多在结构中找到闭合环的任何方法?
我尝试将其转换为图并使用基于图的方法,但是它们都存在一个问题,即图没有节点位置的空间信息,因此图可以有多个同源的环。



由于图形太大,因此不可能找到所有的环,然后过滤掉感兴趣的环。环的大小差别很大。



感谢您的帮助和贡献!



任何语言方法和伪方法-code是可以接受的,尽管我主要使用Python和Matlab。






编辑:



没有图形不是平面的。
图周期基础的问题与其他基于图的简单方法相同。该图缺少任何空间信息,并且不同的空间配置可以具有相同的循环基数,因此循环基数不一定与图中的循环或孔相对应。



此处是稀疏格式的邻接矩阵:

  NodeID1 NodeID2权重




由于某些原因,我无法正确确定周期的基础,但是我认为以下内容一定可以帮助您入门,也许其他人可以加入。



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



 将numpy导入为np 
从skimage.morphology导入为plt
的matplotlib.pyplot,从matplotlib.patches导入medial_axis,binary_closing
.patches导入Path,PathPatch
导入itertools
i mport networkx as nx

img = 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 ='gray')
已关闭= binary_closing(bw,selem = np.ones((50,50)))#连接断开的段
# plt.imshow(关闭,cmap ='灰色')
骨架= medial_axis(关闭)

无花果,ax = plt.subplots(1,1)
ax.imshow(骨骼,cmap ='灰色')
ax.set_xticks([])
ax.set_yticks([])

  def img_to_graph(binary_img, allowed_steps):

参数:
----------
binary_img-二维布尔数组,用于标记节点
的位置-允许的步骤列表;例如[[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-字典将节点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()
for node_pos中的(i,j):
for allow_steps中的(delta_i,delta_j):#如果(i + delta_i,j + pos_to_idx中的delta_j):#即目标节点也存在
g.add_edge(pos_to_idx [(i,j)],pos_to_idx [(i + delta_i,j + delta_j)])

idx_to_pos = dict(zip(zip(node_idx,node_pos))

return g,idx_to_pos,pos_to_idx

allowed_steps = set(itertools.product((-1,0,1),重复= 2))-设置([[0,0) ])
g,idx_to_pos,pos_to_idx = img_to_graph(骨架,允许步数)

图,ax = 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 = [节点的节点,如果度== 2,则度为g.degree()]
链= g.subgraph(is_chain)

#合约连接的组件(应为可变长度的链)到单节点
组件= list(nx.components.connected_component_subgraphs(链))
超节点= g.number_of_nodes()
超节点= []
超边= []
hypernode_to_nodes = dict()
false_alarms = []
用于组件中的组件:
如果component.number_of_nodes()> 1:

hypernodes.append(hypernode)
vs = [用于component.nodes()中节点的节点]
hypernode_to_nodes [hypernode] = vs

#从相邻节点创建新边链结束到超节点
component_edges = [e代表e in component.edges()]
用于v,如果不是,则在[e中用于g.edges(vs)中的e((component_edges中的e或component_edges中的(e [::-1]))]] b $ b如果v在组件中:
hyperedges.append([hypernode,w])
其他:
hyperedges.append([v,hypernode])

hypernode + = 1

else:#不会崩溃,因为组件中只有一个节点:$ b​​ $ b 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(超边缘)

return h,hypernode_to_nodes

h,hypernode_to_nodes =合同(g)

#将超级节点的位置设置为超级节点链的中心
的位置,hypernode_to_nodes.ite中的节点ms():
链= g.subgraph(nodes)
first,last = [节点的节点,chain.degree(),如果度== 1,则度数]
路径= nx。 shortest_path(链,第一个,最后一个)
center = path [len(path)/ 2]
idx_to_pos [hypernode] = idx_to_pos [centre]

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





< h3>查找周期基础

  cycle_basis = nx.cycle_basis(h)

图,ax = plt .subplots(1,1)
nx.draw(h,pos = idx_to_pos,node_size = 10,ax = ax)
在cycle_basis中循环:
个顶点= [idx_to_pos [idx]为idx in cycle]
path =路径(顶点)
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)表示循环中的循环数)#仅冻结集是可散列的
而new_basis!= old_basis :
old_basis = new_basis
for itertools.combinations(old_basis,2)中的cycle_a,cycle_b:如果len(frozenset.union(cycle_a,cycle_b))> = 2:#b $ b检查它们是否共享一条边
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(图,节点) new_basis中的节点]
返回ordered_cycles

def _symmetric_difference(graph,cycle_a,cycle_b):
#获取边
edgesa = list(graph.subgraph(cycle_a)。 edges())
edges_b = list(graph.subgraph(cycle_b).edges())

#也会得到反向边缘,因为图形是无向的
Edges_a + = [e [: :-1] for edges_a]中的e
edgesb + = [e [::-1] edges_b中e的e]

#查找在两个$ b中都没有的边缘$ b edge_c =设置(edges_a)^设置(edges_b)

cycle_c = Frozenset(nx.Graph(list(edges_c))。nodes())
return cycle_c

def _select_cycles(cycles,cost_function):

选择具有周期的节点的覆盖范围,以使与该覆盖范围内的所有周期相关的最大成本
最小化。

周期=清单(周期)
成本= [周期中的cost_function(周期)周期
订单= np.argsort(成本)

节点= Frozenset.union(* cycles)
被覆盖= set()
基础= []

#贪心;从ii的最低成本
开始顺序:
周期=周期[ii]
,如果周期< =被覆盖:
通过
否则:
base.append(cycle)
被覆盖| =循环
(如果已覆盖)==节点:$ b​​ $ b中断

返回集(基础)

def _get_cost(循环,hypernode_to_nodes):
成本= 0
对于周期中的节点:$ b​​ $ b,如果hypernode_to_nodes中的节点:$ b​​ $ b成本+ = len(hypernode_to_nodes [node])
否则:
成本+ = 1
退货成本

def _order_nodes_in_cycle(图,节点):
订单,= nx.cycle_basis(graph.subgraph(节点))
退货订单

个孔= f ind_holes(h,cost_function = partial(_get_cost,hypernode_to_nodes = hypernode_to_nodes))

图,ax = plt.subplots(1,1)
nx.draw(h,pos = idx_to_pos,node_size ii = 10,ax = ax)
,枚举(孔)中的孔:(len(hole)> 3):
个顶点= 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屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆