python - Encontrando vecinos más cercanos a una tesellación triangular.

performance graph-algorithm

Question

Tengo una teselación triangular como la que se muestra en la figura.

Dado el número N de triángulos en la teselación, tengo una matriz NX 3 X 3 que almacena las coordenadas (x, y, z) de los tres vértices de cada triángulo. Mi objetivo es encontrar para cada triángulo el triángulo vecino que comparte el mismo borde. La parte intrincada es toda la configuración de que no repito el recuento de vecinos. Es decir, si el triángulo j ya se contó como un vecino del triángulo i , entonces el triángulo i no debería contarse nuevamente como el vecino del triángulo j . De esta manera, me gustaría tener un mapa que almacene la lista de vecinos para cada triángulo de índice. Si comienzo con un triángulo en el índice i , entonces el índice i tendrá tres vecinos y todos los demás tendrán dos o menos. Como ilustración, supongamos que tengo una matriz que almacena los vértices del triángulo:

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

Supongamos que empiezo mi cuenta desde el índice de vértice 2 , es decir, el que tiene los vértices [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]] , entonces, me gustaría que mi salida ser algo como

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

Actualización: Siguiendo la respuesta de @ Ajax1234, creo que una buena forma de almacenar la salida es como @ Ajax1234 ha demostrado. Sin embargo, existe una ambigüedad en esa salida, en el sentido de que no es posible saber de quién es el vecino. Aunque la matriz de ejemplo no es buena, tengo vértices reales del icosaedro, entonces si comienzo con un triángulo dado, tengo la garantía de tener 3 vecinos para el primero y dos vecinos para el descanso (hasta que todos los recuentos de triángulos se agoten) . En este sentido, supongamos que tengo una matriz siguiente:

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

El algoritmo BFS que se muestra en la respuesta a continuación por @ Ajax1234 da la salida de

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

mientras que si simplemente cambio la posición del último elemento tal que

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

lo que da una salida de

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

Esto es un poco ambiguo, ya que las posiciones en la esquina no se han cambiado en absoluto, simplemente se intercambiaron. Por lo tanto, me gustaría tener una forma consistente en la búsqueda se lleva a cabo. Por ejemplo, la primera búsqueda de vecinos en el índice 2 proporciona [0, 1, 3] para ambos vertices1 y vertices2 , ahora me gustaría que la búsqueda estuviera en el índice 0, que no encuentra nada y, por lo tanto, vaya al siguiente elemento 1 que debería encontrar índice 7 para los vertices1 y el índice 5 para los vertices2 . Por lo tanto, la salida de corriente debe ser [0, 1, 3, 7] , [0, 1, 3, 5] para los vertices1 y vertices2 respectivamente. A continuación vamos al índice 3 , y así sucesivamente. Una vez que hayamos agotado toda la búsqueda, la salida final de la primera debería ser

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

y que para el segundo debe

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

¿Cuál sería la manera eficiente de lograr esto?

Answer 1

Puedes usar Trimesh

Las funciones están documentadas por comentarios en el código fuente. Una documentación normal sería realmente agradable. También lo encuentro no tan sencillo cuando lo usé por primera vez, pero si tienes lo básico, es un paquete poderoso y fácil de usar.

Supongo que el mayor problema es cómo llegar a una definición de malla limpia. Si solo tiene coordenadas de vértice (como en el formato stl), pueden producirse problemas, porque no está bien definido en qué puntos dos flotantes son iguales.

Ejemplo

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

Esto da una matriz 2d simple que las caras comparten un borde. Personalmente usaría este formato para cálculos adicionales. Si desea mantener su definición, crearía una matriz numpy nx3 (cada triángulo debería tener un máximo de 3 vecinos de borde). Los valores que faltan se pueden rellenar, por ejemplo, con NaNs o algo significativo.

Puedo agregar una manera eficiente, si realmente quieres hacer eso.

Answer 2

El proceso que está describiendo suena similar a una búsqueda de amplitud en primer lugar, que se puede utilizar para encontrar los triángulos que son vecinos. Sin embargo, la salida simplemente proporciona los índices, ya que no está claro cómo terminan las listas vacías en la salida final:

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

Salida:

[0, 1, 3, 4, 5]        

Answer 3

Me di cuenta de la respuesta, gracias a la orientación de @ Ajax1234. Había una pequeña complejidad, basada en cómo se comparan los elementos de la lista. Aquí hay un enfoque de trabajo:

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

Para una matriz de tamaño (3000, 3, 3) , el código tarda actualmente 18 segundos en ejecutarse. Si agrego @numba.jit(parallel = True, error_model=''numpy'') , entonces en realidad toma 30 segundos. Probablemente es porque la queue no es compatible con numba . Me alegraría si alguien pudiera sugerir cómo este código podría hacerse más rápido.

Actualizar

Hubo algunas redundancias en el código, que ahora se ha eliminado, y el código se ejecutó en 14 segundos, en lugar de 30 segundos, para un tamaño de d (4000 X 3 X 3) . Todavía no es estelar, pero un buen progreso, y el código parece más limpio ahora.

Answer 4

Si está dispuesto a usar la biblioteca networkx , puede aprovechar su rápida implementación de bfs. Lo sé, agregar otra dependencia es molesto, pero la ganancia de rendimiento parece enorme, ver más abajo.

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

Manifestación:

# 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