Нахождение ближайших соседей треугольной мозаики

Question

У меня треугольная тесселяция, подобная той, что показана на рисунке.

Учитывая N количество треугольников в тесселяции, у меня есть массив N X 3 X 3, которыйхранит (x, y, z) координаты всех трех вершин каждого треугольника.Моя цель - найти для каждого треугольника соседний треугольник, имеющий одинаковое ребро.Это сложная часть всей установки, что я не повторяю количество соседей.То есть, если треугольник j уже был посчитан как сосед треугольника i, то треугольник i не должен снова считаться соседом треугольника j.Таким образом, я хотел бы иметь карту, хранящую список соседей для каждого индексного треугольника.Если я начну с треугольника в индексе i, то у индекса i будет три соседа, а у всех остальных будет два или меньше.В качестве иллюстрации предположим, что у меня есть массив, в котором хранятся вершины треугольника:

import numpy as np
vertices = np.array([[[2.0, 1.0, 3.0],[3.0, 1.0, 2.0],[1.2, 2.5, -2.0]],
                     [[3.0, 1.0, 2.0],[1.0, 2.0, 3.0],[1.2, -2.5, -2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[2.2, 2.0, 1.0]],
                     [[1.0, 2.0, 3.0],[2.2, 2.0, 1.0],[4.0, 1.0, 0.0]],
                     [[2.0, 1.0, 3.0],[2.2, 2.0, 1.0],[-4.0, 1.0, 0.0]]])

Предположим, я начинаю свой отсчет с индекса вершины 2, то есть с вершины [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]], тогда я бынапример, мой вывод будет выглядеть примерно так:

neighbour = [[], [], [0, 1, 3], [4, 5], [], []].

Обновление: После ответа от @ Ajax1234, я думаю, что хороший способ хранения вывода такой же, как продемонстрировал @ Ajax1234.Однако в этом выводе есть двусмысленность, в том смысле, что невозможно узнать, чей сосед какой.Хотя пример массива не очень хорош, у меня есть фактические вершины из икосаэдра, но если я начну с заданного треугольника, у меня гарантированно будет 3 соседа для первого и два соседа для отдыха (пока все треугольники не истощатся),В связи с этим предположим, что у меня есть следующий массив:

vertices1 = [[[2, 1, 3], [3, 1, 2], [1, 2, -2]], 
            [[3, 1, 2], [1, 2, 3], [1, -2, 2]], 
            [[1, 2, 3], [2, 1, 3], [3, 1, 2]], 
            [[1, 2, 3], [2, 1, 3], [2, 2, 1]],
            [[1, 2, 3], [2, 2, 1], [4, 1, 0]], 
            [[2, 1, 3], [2, 2, 1], [-4, 1, 0]],
            [[3, 1, 3], [2, 2, 1], [-4, 1, 0]],
            [[8, 1, 2], [1, 2, 3], [1, -2, 2]]]

Алгоритм BFS, показанный в ответе ниже @ Ajax1234, дает вывод

[0, 1, 3, 7, 4, 5, 6]

, а если я просто поменяю местамипозиция последнего элемента такая, что

vertices2 = [[[2, 1, 3], [3, 1, 2], [1, 2, -2]], 
            [[3, 1, 2], [1, 2, 3], [1, -2, 2]], 
            [[1, 2, 3], [2, 1, 3], [3, 1, 2]], 
            [[1, 2, 3], [2, 1, 3], [2, 2, 1]],
            [[1, 2, 3], [2, 2, 1], [4, 1, 0]], 
            [[8, 1, 2], [1, 2, 3], [1, -2, 2]],
            [[2, 1, 3], [2, 2, 1], [-4, 1, 0]],
            [[3, 1, 3], [2, 2, 1], [-4, 1, 0]]]

, которая дает вывод

[0, 1, 3, 4, 5, 6, 7].

Это несколько неоднозначно, так как позиции в гирде не изменились вообще, онибыли просто поменяны местами.Поэтому я хотел бы иметь последовательный способ поиска.Например, первый раз поиск соседей по индексу 2 дает [0, 1, 3] для vertices1 и vertices2, теперь я хотел бы, чтобы поиск был по индексу 0, который ничего не находит, и, следовательно, переход к следующему элементу 1 должен найти индекс7 для vertices1 и индекс 5 для vertices2.Таким образом, выходной ток должен быть [0, 1, 3, 7], [0, 1, 3, 5] для vertices1 и vertices2 соответственно.Далее мы идем к индексу 3 и так далее.После того, как мы исчерпали весь поиск, окончательный результат для первого должен быть

[0, 1, 3, 7, 4, 5, 6]

, а для второго -

[0, 1, 3, 5, 4, 6, 7].

Каким был бы эффективный способ добиться этого?

Answer 1

Вы можете использовать trimesh

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

Я предполагаю, что самая большая проблема заключается в том, как добраться доопределения чистой сетки.Если у вас есть только координаты вершин (как в stl-формате), могут возникнуть проблемы, потому что неясно, в каких точках равны два числа с плавающей точкой.

Пример

import trimesh
import numpy as np

vertices = np.array([[[2.0, 1.0, 3.0],[3.0, 1.0, 2.0],[1.2, 2.5, -2.0]],
                     [[3.0, 1.0, 2.0],[1.0, 2.0, 3.0],[1.2, -2.5, -2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[2.2, 2.0, 1.0]],
                     [[1.0, 2.0, 3.0],[2.2, 2.0, 1.0],[4.0, 1.0, 0.0]],
                     [[2.0, 1.0, 3.0],[2.2, 2.0, 1.0],[-4.0, 1.0, 0.0]]])

#generate_faces
# I assume here that your input format is N x NPoints x xyz
faces=np.arange(vertices.shape[0]*3).reshape(-1,3)
#reshape_vertices (nx3)
vertices=vertices.reshape(-1,3)

#Create mesh object
mesh=trimesh.Trimesh(vertices=vertices, faces=faces)

# Set the tolerance to define which vertices are equal (default 1e-8)
# It is easy to prove that int(5)==int(5) but is 5.000000001 equal to 5.0 or not?
# This depends on the algorithm/ programm from which you have imported the mesh....
# To find a proper value for the tolerance and to heal the mesh if necessary, will 
# likely be the most complicated part
trimesh.constants.tol.merge=tol

#merge the vertices
mesh.merge_vertices()

# At this stage you may need some sort of healing algorithm..

# eg. reconstruct the input
mesh.vertices[mesh.faces]

#get for example the faces, vertices
mesh.faces #These are indices to the vertices
mesh.vertices

# get the faces which touch each other on the edges
mesh.face_adjacency

Это дает простой 2d массив, грани которого имеют общий край.Я бы лично использовал этот формат для дальнейшего расчета.Если вы хотите придерживаться своего определения, я бы создал массив nx3 (у каждого треугольника должно быть не более 3 ребер-соседей).Пропущенные значения могут быть заполнены, например, NaN или чем-то значимым.

Я могу добавить эффективный способ, если вы действительно хотите это сделать.

Answer 2

Я разобрался с ответом, благодаря руководству @ Ajax1234.Была небольшая сложность, основанная на том, как вы сравниваете элементы списка.Вот один из рабочих подходов:

import numpy as np
from collections import deque
import time
d = [[[2, 1, 3], [3, 1, 2], [1, 2, -2]], 
     [[3, 1, 2], [1, 2, 3], [1, -2, 2]], 
     [[1, 2, 3], [2, 1, 3], [3, 1, 2]], 
     [[1, 2, 3], [2, 1, 3], [2, 2, 1]],
     [[1, 2, 3], [2, 2, 1], [4, 1, 0]],
     [[2, 1, 3], [2, 2, 1], [-4, 1, 0]],
     [[3, 1, 3], [2, 2, 1], [-4, 1, 0]]]
def bfs(d, start):
  queue = deque([d[start]])
  seen = [start]
  results = []
  while queue:
    _vertices = queue.popleft()
    current = [[i, a] for i, a in enumerate(d) if len([x for x in a if x in _vertices])==2 and i not in seen]
    if len(current)>0:
        current_array = np.array(current, dtype=object)
        current0 = list(current_array[:, 0])
        current1 = list(current_array[:, 1])
        results.extend(current0)
        queue.extend(current1)
        seen.extend(current0)
  return results

time1 = time.time()
xo = bfs(d, 2)
print(time.time()-time1)
print(bfs(d, 2))

Для массива размером (3000, 3, 3) код в настоящее время занимает 18 секунд для запуска.Если я добавлю @numba.jit(parallel = True, error_model='numpy'), то на самом деле это займет 30 секунд.Возможно, это связано с тем, что queue не поддерживается numba.Я был бы рад, если бы кто-нибудь мог предложить, как этот код может быть сделан быстрее.

Обновление

Были некоторые избыточности в коде, который теперь был удален, и код выполняется за 14 секунд вместо 30 секунд,для d размером (4000 X 3 X 3).Все еще не звездный, но хороший прогресс, и код теперь выглядит чище.

Answer 3

Если вы хотите использовать библиотеку networkx, вы можете воспользоваться ее быстрой реализацией bfs.Я знаю, добавление еще одной зависимости раздражает, но прирост производительности кажется огромным, см. Ниже.

import numpy as np
from scipy import spatial
import networkx

vertices = np.array([[[2.0, 1.0, 3.0],[3.0, 1.0, 2.0],[1.2, 2.5, -2.0]],
                     [[3.0, 1.0, 2.0],[1.0, 2.0, 3.0],[1.2, -2.5, -2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]],
                     [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[2.2, 2.0, 1.0]],
                     [[1.0, 2.0, 3.0],[2.2, 2.0, 1.0],[4.0, 1.0, 0.0]],
                     [[2.0, 1.0, 3.0],[2.2, 2.0, 1.0],[-4.0, 1.0, 0.0]]])


vertices1 = np.array([[[2, 1, 3], [3, 1, 2], [1, 2, -2]], 
                      [[3, 1, 2], [1, 2, 3], [1, -2, 2]], 
                      [[1, 2, 3], [2, 1, 3], [3, 1, 2]], 
                      [[1, 2, 3], [2, 1, 3], [2, 2, 1]],
                      [[1, 2, 3], [2, 2, 1], [4, 1, 0]], 
                      [[2, 1, 3], [2, 2, 1], [-4, 1, 0]],
                      [[3, 1, 3], [2, 2, 1], [-4, 1, 0]],
                      [[8, 1, 2], [1, 2, 3], [1, -2, 2]]])


def make(N=3000):
    """create a N random points and triangulate"""
    points= np.random.uniform(-10, 10, (N, 3))
    tri = spatial.Delaunay(points[:, :2])
    return points[tri.simplices]

def bfs_tree(triangles, root=0, return_short=True):
    """convert triangle list to graph with vertices = triangles,
    edges = pairs of triangles with shared edge and compute bfs tree
    rooted at triangle number root"""
    # use the old view as void trick to merge triplets, so they can
    # for example be easily compared
    tr_as_v = triangles.view(f'V{3*triangles.dtype.itemsize}').reshape(
        triangles.shape[:-1])
    # for each triangle write out its edges, this involves doubling the
    # data becaues each vertex occurs twice
    tr2 = np.concatenate([tr_as_v, tr_as_v], axis=1).reshape(-1, 3, 2)
    # sort vertices within edges ...
    tr2.sort(axis=2)
    # ... and glue them together
    tr2 = tr2.view(f'V{6*triangles.dtype.itemsize}').reshape(
        triangles.shape[:-1])
    # to find shared edges, sort them ...
    idx = tr2.ravel().argsort()
    tr2s = tr2.ravel()[idx]
    # ... and then compare consecutive ones
    pairs, = np.where(tr2s[:-1] == tr2s[1:])
    assert np.all(np.diff(pairs) >= 2)
    # these are the edges of the graph, the vertices are implicitly 
    # named 0, 1, 2, ...
    edges = np.concatenate([idx[pairs,None]//3, idx[pairs+1,None]//3], axis=1)
    # construct graph ...
    G = networkx.Graph(edges.tolist())
    # ... and let networkx do its magic
    res = networkx.bfs_tree(G, root)
    if return_short:
        # sort by distance from root and then by actual path
        sp = networkx.single_source_shortest_path(res, root)
        sp = [sp[i] for i in range(len(sp))]
        sp = [(len(p), p) for p in sp]
        res = sorted(range(len(res.nodes)), key=sp.__getitem__)
    return res

Демонстрация:

# OP's second example:
>>> bfs_tree(vertices1, 2)
[2, 0, 1, 3, 7, 4, 5, 6]
>>> 
# large random example
>>> random_data = make()
>>> random_data.shape
(5981, 3, 3)
>>> bfs = bfs_tree(random_data)
# returns instantly

Answer 4

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

from collections import deque
d = [[[2.0, 1.0, 3.0], [3.0, 1.0, 2.0], [1.2, 2.5, -2.0]], [[3.0, 1.0, 2.0], [1.0, 2.0, 3.0], [1.2, -2.5, -2.0]], [[1.0, 2.0, 3.0], [2.0, 1.0, 3.0], [3.0, 1.0, 2.0]], [[1.0, 2.0, 3.0], [2.0, 1.0, 3.0], [2.2, 2.0, 1.0]], [[1.0, 2.0, 3.0], [2.2, 2.0, 1.0], [4.0, 1.0, 0.0]], [[2.0, 1.0, 3.0], [2.2, 2.0, 1.0], [-4.0, 1.0, 0.0]]]
def bfs(d, start):
  queue = deque([d[start]])
  seen = [start]
  results = []
  while queue:
    _vertices = queue.popleft()
    exists_at = [i for i, a in enumerate(d) if a == _vertices][0]
    current = [i for i, a in enumerate(d) if any(c in a for c in _vertices) and i != exists_at and i not in seen]
    results.extend(current)
    queue.extend([a for i, a in enumerate(d) if any(c in a for c in _vertices) and i != exists_at and i not in seen])
    seen.extend(current)
  return results

print(bfs(d, 2))

Вывод:

[0, 1, 3, 4, 5]