Finding Nearest Neighbours Of A Triangular Tesellation
Solution 1:
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.
Solution 2:
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]]])
defmake(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]
defbfs_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 inrange(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
Solution 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 notin 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 notin seen])
seen.extend(current)
return results
print(bfs(d, 2))
Output:
[0, 1, 3, 4, 5]
Solution 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.
Post a Comment for "Finding Nearest Neighbours Of A Triangular Tesellation"