寻找三角形镶嵌的最近邻居
Finding nearest neighbours of a triangular tesellation
我有一个如图所示的三角形镶嵌。
给定镶嵌中的 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]]]
@Ajax1234 在下面的答案中显示的 BFS 算法给出了
的输出
[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 处搜索邻居时，vertices1 和 vertices2 都会得到 [0, 1, 3]，现在我希望搜索在索引 0 处，它什么也找不到，因此转到下一个元素 1 应该为 vertices1 找到索引 7，为 vertices2 找到索引 5。因此，对于 vertices1 和 vertices2，当前输出应该分别为 [0, 1, 3, 7]、[0, 1, 3, 5]。接下来我们转到索引 3，依此类推。在我们用尽所有搜索之后，第一个的最终输出应该是
[0, 1, 3, 7, 4, 5, 6]
第二个应该
[0, 1, 3, 5, 4, 6, 7].
实现此目标的有效方法是什么？
您描述的过程听起来类似于广度优先搜索，可用于查找相邻的三角形。然而，输出仅给出索引，因为不清楚空列表如何在最终输出中结束：
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]
感谢@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 秒。可能是因为numba不支持queue。如果有人可以建议如何使这段代码更快，我将很高兴。
更新
代码中有一些冗余，现在已被删除，代码 运行 在 14 秒内，而不是 30 秒，持续 [=19] =] 大小 (4000 X 3 X 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
你可以使用 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
这给出了一个简单的二维数组，其面共享一条边。我个人会使用这种格式进行进一步计算。如果你想坚持你的定义，我会创建一个 nx3 numpy 数组（每个三角形应该有最多 3 个边缘邻居）。缺失值可以用例如 NaN 或其他有意义的东西来填充。
如果你真的想这样做，我可以添加一个有效的方法。
我有一个如图所示的三角形镶嵌。
给定镶嵌中的 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]]]
@Ajax1234 在下面的答案中显示的 BFS 算法给出了
的输出[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 处搜索邻居时，vertices1 和 vertices2 都会得到 [0, 1, 3]，现在我希望搜索在索引 0 处，它什么也找不到，因此转到下一个元素 1 应该为 vertices1 找到索引 7，为 vertices2 找到索引 5。因此，对于 vertices1 和 vertices2，当前输出应该分别为 [0, 1, 3, 7]、[0, 1, 3, 5]。接下来我们转到索引 3，依此类推。在我们用尽所有搜索之后，第一个的最终输出应该是
[0, 1, 3, 7, 4, 5, 6]
第二个应该
[0, 1, 3, 5, 4, 6, 7].
实现此目标的有效方法是什么？
您描述的过程听起来类似于广度优先搜索，可用于查找相邻的三角形。然而，输出仅给出索引，因为不清楚空列表如何在最终输出中结束：
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]
感谢@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 秒。可能是因为numba不支持queue。如果有人可以建议如何使这段代码更快，我将很高兴。
更新
代码中有一些冗余，现在已被删除，代码 运行 在 14 秒内，而不是 30 秒，持续 [=19] =] 大小 (4000 X 3 X 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
你可以使用 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
这给出了一个简单的二维数组，其面共享一条边。我个人会使用这种格式进行进一步计算。如果你想坚持你的定义，我会创建一个 nx3 numpy 数组（每个三角形应该有最多 3 个边缘邻居）。缺失值可以用例如 NaN 或其他有意义的东西来填充。
如果你真的想这样做，我可以添加一个有效的方法。