检测 rings/circuits 个连接的体素
Detect rings/circuits of connected voxels
我有一个骨架化的体素结构，如下所示：
实际结构比这个例子大很多有没有办法找到结构中的闭环？
我尝试将其转换为图形并使用基于图形的方法，但它们都有一个问题，即图形没有节点位置的空间信息，因此图形可以具有多个同源环。
不可能找到所有环然后过滤掉感兴趣的环，因为图太大了。戒指的大小差异很大。
感谢您的帮助和贡献！
尽管我主要在 Python 和 Matlab 中工作，但欢迎任何语言方法和伪代码。
编辑：
不，该图不是平面的。
Graph cycle base 的问题与其他简单的基于图的方法相同。该图缺少任何空间信息，不同的空间配置可以具有相同的循环基础，因此循环基础不一定对应于图中的循环或空洞。
这是稀疏格式的邻接矩阵：
NodeID1 NodeID2 Weight
Pastebin with adjacency matrix
这里是图表节点对应的 X、Y、Z 坐标：
X Y Z
Pastebin with node coordinates
（实际结构明显比这个例子大）
首先，我通过将 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))
我有一个骨架化的体素结构，如下所示：
实际结构比这个例子大很多有没有办法找到结构中的闭环？ 我尝试将其转换为图形并使用基于图形的方法，但它们都有一个问题，即图形没有节点位置的空间信息，因此图形可以具有多个同源环。
不可能找到所有环然后过滤掉感兴趣的环，因为图太大了。戒指的大小差异很大。
感谢您的帮助和贡献！
尽管我主要在 Python 和 Matlab 中工作，但欢迎任何语言方法和伪代码。
编辑：
不，该图不是平面的。 Graph cycle base 的问题与其他简单的基于图的方法相同。该图缺少任何空间信息，不同的空间配置可以具有相同的循环基础，因此循环基础不一定对应于图中的循环或空洞。
这是稀疏格式的邻接矩阵：
NodeID1 NodeID2 Weight
Pastebin with adjacency matrix
这里是图表节点对应的 X、Y、Z 坐标：
X Y Z
Pastebin with node coordinates
（实际结构明显比这个例子大）
首先，我通过将 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))