Finding nearest neighbours of a triangular tesellation

Question

I have a triangular tessellation like the one shown in the figure.

Given N number of triangles in the tessellation, I have a N X 3 X 3 array which stores (x, y, z) coordinates of all three vertices of each triangle. My goal is to find for each triangle the neighbouring triangle sharing the same edge. The is an intricate part is the whole setup that I do not repeat the neighbour count. That is if triangle j was already counted as a neighbour of triangle i, then triangle i should not be again counted as neighbour of triangle j. This way, I would like to have a map storing list of neighbours for each index triangle. If I start with a triangle in index i, then index i will have three neighbours, and all others will have two or less. As an illustration suppose I have an array which stores vertices of the triangle:

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

Suppose I start my count from vertex index 2, i.e. the one with the vertices [[1.0, 2.0, 3.0],[2.0, 1.0, 3.0],[3.0, 1.0, 2.0]], then, I would like my output to be something like:

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

Update: Following the answer from @Ajax1234, I think a good way of storing the output is just like how @Ajax1234 has demonstrated. However, there is ambiguity in that output, in a sense that it is not possible to know whose neighbour is which. Although the example array are not good, I have an actual vertices from icosahedron, then if I start with a given triangle, I am guaranteed to have 3 neighbours for the first one, and two neighbours for rest (until all the triangle counts deplete). In this regard, suppose I have a following array:

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

The BFS algorithm shown in the answer below by @Ajax1234 gives the output of

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

while if I just swap the position of the last element such that

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

which gives an output of

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

This is kind of ambiguous, as the positions in the gird have not been changed at all, they were just swapped. Therefore, I would like to have a consistent way the search is carried. For example, first time search of neighbours at index 2 gives [0, 1, 3] for both vertices1 and vertices2, now I would like the search to be at index 0, which finds nothing and thus go to next element 1 should find index 7 for vertices1 and index 5 for vertices2. Thus the current output should be [0, 1, 3, 7], [0, 1, 3, 5] for vertices1 and vertices2 respectively. Next we go to index 3, and so on. After we have exhausted all the search, the final output for the first should be

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

and that for the second should

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

What would be the efficient way to achieve this?

Intrepid

What is the rule for the inclusion of [] in the neighbour list? – Beamer

It means that the particular index triangle has no neighbours. – Intrepid

You can do that with trimesh github.com/mikedh/trimesh. In general I would convert your definition of the mesh to vertices and faces, which is a lot more stable. – Laccolith

@Laccolith I was also looking at the same github package. Got a bit confused on exactly which part of the package achieves that, and how. But, thanks. – Intrepid

Answer 1

If you are willing to use the networkx library, you can take advantage of its fast bfs implementation. I know, adding another dependency is annoying, but the performance gain seems huge, see below.

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

Demo:

# 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

Appropriation

Thanks Paul. I think your code is meant to be be Python 3, but I am running 2.7, and it does not support .view(f'...') type code. Still trying to figure how to make this compatible in Python 2.7. – Intrepid

Thanks. Sorry to bother you again, it throws one more error TypeError: data type "V{3*triangles.dtype.itemsize}" not understood. Could you please suggest me how to fix this? – Intrepid

Actually, the quotation mark was just in the error message. The full message is File "test.py", line 25, in bfs_tree tr_as_v = triangles.view('V{3*triangles.dtype.itemsize}').reshape(triangles.shape[:-1]) TypeError: data type "V{3*triangles.dtype.itemsize}" not understood. – Intrepid

Paul, this works perfectly. But I have a question, is it possible to modify this code such that I get a desired output of what I posted in my updated question for vertices1. In particular, for that case, your code gives an output of [(1, 7), (2, 0), (2, 1), (2, 3), (3, 4), (3, 5), (5, 6)], which is accurate, but I want it in a linear sequence like [0, 1, 3, 7, 4, 5, 6]. – Intrepid

@Intrepid Currently the code is just going with the format networkx is using. One would have to manually convert. Problem is, I don't understand how your format works. For example, how do you know that 7 connects to 1 and not to 0? – Appropriation

the current example is not so good in that sense that I do not know if 7 connects to 1 or 0. But in the real code I am dealing with, e very triangle is guaranteed to have two neighbours (except for the edge ones), so it is non-ambiguous in my real example case. – Intrepid

@Intrepid I think I have implemented your format. Only difference is that the very first element is the root node. – Appropriation

Wow, that's really cool. Actually I am getting an error, though. Again, should be an issue of python version. res = sorted(range(len(res.nodes)), key=sp.__getitem__) TypeError: object of type 'instancemethod' has no len(). Could you please suggest me how to fix this? – Intrepid

@Intrepid looks like in your setup res.nodes is a function. What happens if you replace res.nodes with res.nodes()? – Appropriation

Answer 2

I figured out the answer, thanks to the guidance of @Ajax1234. There was a small intricacy, based on how you compare the list elements. Here is one working approach:

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

For an array of size (3000, 3, 3), the code currently takes 18 seconds to run. If I add @numba.jit(parallel = True, error_model='numpy'), then it actually takes 30 seconds. It is probably because queue is not supported by numba. I would be pleased if any one could suggest how this code could be made faster.

Update

There were some redundancies in the code, which has now been removed, and the code run in 14 seconds, instead of 30 seconds, for a of d of size (4000 X 3 X 3). Still not stellar, but a good progress, and the code looks cleaner now.

Intrepid

Answer 3

The process you are describing sounds similar to a breadth-first search, which can be utilized to find the triangles which are neighbors. The output merely gives the indices, however, since it is unclear as to how the empty lists end up in the final output:

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

Output:

[0, 1, 3, 4, 5]        

Beamer

The empty index suggest that the given index has no neighbour. In the particular example, triangle 0, triangle 1, triangle 5 and triangle 6 have no neighbours. – Intrepid

In your code, how do you know which are the neighbours of a given index? It is pretty clear for now that 0, 1, 3 are neighbours of index 2, but how do I know if 4, 5 are neighbours of index 0 or index 1 or index 3? – Intrepid

Answer 4

You can use trimesh

The functions are documented by comments in the source code. A normal documentation would be realy nice. I also find it not that straight forward when I used it the first time, but if you have the basics it is a powerful and easy to use package.

I assume that the biggest problem is how to get to a clean mesh definitions. If you have only vertex coordinates (like in the stl-format) problems may occur, because it is not well defined at which points two floats are equal.

Example

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

This gives a simple 2d array which faces share a edge. I would personally use this format for further calculation. If you want to stick to your definition I would create a nx3 numpy array (each triangle should have max. 3 edge-neighbours). The missing values can be filled for example with NaNs or something meaningfull.

I can add an efficient way, if you really want to do that.

Laccolith

Thanks. Personally, I would like the output to be of the format of answer of Paul Panzer. – Intrepid