Обнаружение колец / цепей подключенных вокселей

У меня есть скелетированная структура вокселей, которая выглядит следующим образом:  пример воксельных колец

Фактическая структура значительно больше, чем в этом примере. Есть ли способ найти замкнутые кольца в структуре? Я попытался преобразовать его в граф и использовать подходы, основанные на графах, но все они имеют проблему, заключающуюся в том, что граф не имеет пространственной информации о положении узла и, следовательно, граф может иметь несколько колец, которые являются гомологичными.

Невозможно найти все кольца, а затем отфильтровать те, которые представляют интерес, поскольку граф слишком велик. Размер колец значительно различается.

Спасибо за вашу помощь и вклад!

Приветствуются любые языковые подходы и псевдокод, хотя я работаю в основном на Python и Matlab.


РЕДАКТИРОВАТЬ:

Нет, график не плоский. Проблема с базой цикла Graph такая же, как и с другими простыми подходами, основанными на графах. В графе отсутствует какая-либо пространственная информация, и разные пространственные конфигурации могут иметь одну и ту же базу цикла, следовательно, база цикла не обязательно соответствует циклам или дырам в графе.

Вот матрица смежности в разреженном формате:

NodeID1 NodeID2 Weight

Pastebin с матрицей смежности

А вот соответствующие координаты X, Y, Z узлов графика:

X Y Z

Pastebin с координатами узла

(Фактическая структура значительно больше, чем в этом примере)


person DarkCell    schedule 26.04.2017    source источник
comment
Вы ищете не кольца, а дырочки.   -  person Paul Brodersen    schedule 26.04.2017
comment
Можете ли вы опубликовать пример структуры вокселей в виде разреженной матрицы?   -  person Paul Brodersen    schedule 26.04.2017
comment
Граф плоский? Другими словами, можно ли проецировать координаты x, y, z на двумерное многообразие без изменения топологии графа? Какую структуру мы смотрим? Это какая-то складчатая мембрана?   -  person Paul Brodersen    schedule 26.04.2017


Ответы (1)


Сначала я значительно уменьшаю размер проблемы, сжимая соседние узлы степени 2 в гиперузлы: каждая простая цепочка в графе заменяется одним узлом.

Затем я нахожу базис цикла, для которого максимальная стоимость циклов в базисном наборе составляет минимальный.

Для центральной части сети решение может быть легко построено, так как оно является плоским:

введите здесь описание изображения

По какой-то причине мне не удается правильно определить основу цикла, но я думаю, что нижеследующее определенно должно помочь вам начать работу, и, возможно, кто-то еще сможет вмешаться.

Восстановить данные из опубликованного изображения (так как 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)

введите здесь описание изображения

NB: это не красные линии, это множество красных точек, соответствующих узлам на графике.

График контракта

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

ДЕЛАТЬ:

Найдите правильную основу цикла (я могу не понять, что такое базис цикла или 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))
person Paul Brodersen    schedule 27.04.2017
comment
Большое спасибо за очень подробный ответ, но, к сожалению, с ним есть некоторые проблемы: Показанный вами подход я уже реализовал. Я могу создать очень чистый трехмерный график из моих исходных данных вокселей. (Возможно, это было не так ясно, и я приложил данные к своему вопросу). Проблема с подходами, основанными на графах, заключается в том, что им не хватает какой-либо пространственной информации об узлах. Следовательно, такие вещи, как основа цикла, на самом деле не определены должным образом. Рассмотрим следующие графики. Они представляют собой графы одинаково! - person DarkCell; 27.04.2017
comment
Конечно, но вы можете просто отслеживать, какой узел соответствует какой позиции, как это делаю я в своем коде со словарем idx_to_pos. - person Paul Brodersen; 27.04.2017
comment
да, вы можете отслеживать их, но код для вычисления основы цикла не использует эту информацию в каком-либо смысле, и никакой алгоритм чистого графа никогда не будет (вы можете получить часть этой информации, имея веса ребер, но не полную информацию). Когда вы посмотрите на мой рисунок, основа цикла будет одинаковой в обоих случаях, хотя структура, которая меня интересует, будет другой. Вот почему базис цикла, который вы получаете из networkx, кажется, не соответствует кольцам / дырам, хотя они являются правильным основанием цикла. - person DarkCell; 27.04.2017
comment
Если вы отслеживаете позиции узлов, вы можете просто (-ish) наложить штраф на базисные наборы цикла, которые вам не нравятся. См. Обновленный ответ. - person Paul Brodersen; 27.04.2017
comment
Ну, прежде всего: большое спасибо за этот отличный подробный ответ. Я пробовал делать dfs и bfs с некоторым штрафом за пройденное расстояние, но мне никогда не приходило в голову собрать базу для велосипеда с наименьшими затратами. Я не пробовал это на моем полном графике, но я поиграю в ближайшие дни. Это выглядит многообещающе, но да, это определенно кроличья нора :) минимизация окружности не идеальна, поскольку отверстия могут быть вогнутыми, но у меня в голове уже есть некоторые функции стоимости, которые могут сработать. Я мог бы вернуться к вам, но пока я собираюсь принять этот отличный ответ. Есть новые идеи :) - person DarkCell; 03.05.2017
comment
Да, я подумал, что вы можете использовать другую функцию стоимости, которая больше связана с площадью колец / отверстий. Кроме того, если какой-либо вариант этого окажется в исследовательской статье, я хочу знать об этом! - person Paul Brodersen; 03.05.2017
comment
Я играю как с реализацией в Matlab, так и с вашей реализацией на Python. Мне было интересно, верен ли ваш критерий для прерывания цикла выбора. В случае, если у меня есть только узлы со степенью ›2 (без узлов на линиях), разрыв, если все узлы посещены, оставляет ребра, которые были посещены ранее. Если я изменю это на все посещенные ребра, которые были посещены до того, как я посещу все ребра и, следовательно, все узлы, но у меня все равно будет меньше циклов, чем раньше. Но я думал, что количество базовых компонентов цикла никогда не должно меняться? Как я знаю, что нашел подходящую базу для цикла? - person DarkCell; 10.05.2017
comment
Строка new_basis = new_basis.union([cycle_c]) добавляет кучу циклов сверх базового набора. Функция select_cycles сжимает, что снова устанавливается на базисный набор. Разрыв в select_cycles гарантирует, что у нас есть правильный базовый набор. - person Paul Brodersen; 11.05.2017
comment
Да, я полностью понимаю. Но вы не можете просто проверить, все ли узлы были посещены. После вчерашнего разведения я уверен, что вам нужно проверить, является ли цикл линейно независимым, прежде чем добавлять. Допустим, вы начинаете с самого дешевого цикла, добавляете второй, а затем добавляете третий, который по-прежнему дешев, но на самом деле представляет собой линейную комбинацию первого. Пример: представьте себе квадрат с крестом внутри: ☒ У вас есть 5 узлов и 8 ребер, поэтому 8-5 + 1 = 4 цикла в качестве основы. Если вы добавите левый треугольник и правый треугольник, у вас будут посещены все узлы, но не база полного цикла. - person DarkCell; 11.05.2017
comment
Ты прав. Не могу придумать решение сегодня, но на выходных подумаю. - person Paul Brodersen; 11.05.2017
comment
Думаю, я решил проблему: вместо того, чтобы отслеживать узлы для циклов, я создаю вектор инцидентности для каждого цикла. У меня m ребер. Цикл теперь является вектором типа bool размером 1xm с единицами, если ребро находится в цикле, и нулями в противном случае. Тогда симметричная разность - это побитовое xor, и я могу использовать логическое гауссовское исключение после добавления нового цикла и проверить, получил ли в результате нулевой вектор. Если это так, цикл линейно зависит от уже найденной базы и будет отклонен, а затем я могу выйти из цикла, когда будут найдены N циклов (или перебрать остальные, чтобы проверить) - person DarkCell; 12.05.2017
comment
Это сработает, и это намного проще и, вероятно, быстрее, чем моя идея. Я собирался изменить код таким образом, чтобы, получив начальную основу цикла, я отслеживал все возможности в своего рода мультивселенной решений, а затем использовал Витерби / Баума-Велча, чтобы найти оптимальный набор циклов. Вроде как вы ищете оптимальную последовательность в скрытой марковской модели. - person Paul Brodersen; 12.05.2017
comment
На самом деле это может быть полезно для людей. Нам нужно создать репозиторий на github. Вы пишете код Matlab, я пишу код Python, а затем мы выбираем общий API. - person Paul Brodersen; 12.05.2017
comment
Да, мы должны это сделать. В настоящее время я много оптимизирую, так как мне это понадобится для работы с большим графом, и, поскольку это довольно грубая сила, есть некоторая оптимизация в отношении памяти и типов данных, которые необходимы. Возможно, версии на python это не понадобится, раз уж у вас есть наборы, но мы можем это увидеть :) - person DarkCell; 16.05.2017
comment
В моем профиле указан мой веб-сайт, на котором указан мой адрес электронной почты. Отправь мне электронное письмо, когда захочешь поговорить. - person Paul Brodersen; 16.05.2017
comment
Очень интересно написать! Оказалось ли это в базе кода, как предложил @Paul? Я хотел бы взглянуть и / или внести свой вклад, как только я протестирую это с моим вариантом использования - person Tom Hemmes; 04.07.2018
comment
@TomHemmes Нет, извините, мы так и не разобрались с этим полностью. Мне пришлось закончить докторскую степень, и в процессе это выпало из моего поля зрения. DarkCell также обнаружила некоторые проблемы с моим подходом, поэтому ваш опыт может отличаться. Я опубликую его электронную почту в комментариях ниже. Я думаю, что мой подход по-прежнему является хорошей эвристикой, но вы должны внимательно смотреть на свои результаты, поскольку нет никаких гарантий. - person Paul Brodersen; 05.07.2018
comment
Основная проблема, которую я до сих пор обнаружил, - это вычислительные затраты на определение того, является ли кольцо линейно независимым до их добавления, а другая проблема заключается в том, что ваш подход не обязательно сходится или дает оптимальное решение. Т.е. может случиться так, что если вы возьмете 3 круга со стоимостью a, b, c, новое кольцо a и b будет иметь более высокую стоимость, чем a, b и c, но новое кольцо (a & b) и c может иметь более низкую стоимость, чем a, б и в. Но этого никогда не происходит, потому что кольцо a и b исключено из-за его более высокой стоимости. - person Paul Brodersen; 05.07.2018
comment
Основная проблема здесь в том, что для каждой функции стоимости, кроме функции с минимальным весом, проблема NP-трудна. При тщательной инициализации подход все же дает хорошее приближение, поэтому я думаю, что стоит продолжить. - person Paul Brodersen; 05.07.2018