在 2D 三角形网格中查找对齐和几乎对齐的线段
Find aligned and nearly-aligned segments in a 2D triangular mesh
我有一个结构为 (x, y, index)
的二维点云
示例：
pts = [[0,0,1],[0,1,2],[0,2,3],[1,0,4],[1,1,5],[1,2,6],[2,0,7],[2,1,8],[2,2,9]]
和一系列具有结构 (index1, index2) 的线段，将这些点 2 乘 2 连接起来：
segs = [[1,2],[1,4],[1,5],[2,3],[2,5],[3,5],[3,6],[4,5],[4,7],[4,8],[5,6],[5,8],[5,9],[6,9],[7,8],[8,9]]
线段总是只连接两个点，它们永远不会经过第三个点。
当我应该只使用 seg 中的段一次时，如何检测对齐的段并使用 n 索引制作更长的段？
aligned = [[1,2,3],[1,4,7],[1,5,9],[2,5,8],[3,6,9],[3,5],[4,5,6],[4,8],[7,8,9]]
现在说点没有完全对齐：
pts = [[0.002,0.1,1],[0.0003,1.003,2],[-0.01,2.11,3],[1.01,0.101,4],[1.0005,1.2,5],[1.25,2.01,6],[2.007,-0.12,7],[1.996,1.1,8],[2.03,2.1,9]]
如果线段之间的角度几乎为零，我想连接几乎对齐的线段，由 tolerance 定义，例如，3°
我试过很多方法都没有成功。这必须在 10k+ 点的大量网络上完成。
此代码重现了我的示例：
import numpy as np
import matplotlib.pyplot as plt
#pts = np.array([[0,0,1],[0,1,2],[0,2,3],[1,0,4],[1,1,5],[1,2,6],[2,0,7],[2,1,8],[2,2,9]])
pts = np.array([[0.002,0.1,1],[0.0003,1.003,2],[-0.01,2.11,3],[1.01,0.101,4],[1.0005,1.2,5],[1.25,2.01,6],[2.007,-0.12,7],[1.996,1.1,8],[2.03,2.1,9]])
indices = pts[:,2].astype(int)
seg = np.array([[1,2],[1,4],[1,5],[2,3],[2,5],[3,5],[3,6],[4,5],[4,7],[4,8],[5,6],[5,8],[5,9],[6,9],[7,8],[8,9]])
for s in seg:
i1, i2 = s
segment = np.array((pts[indices == i1][0], pts[indices == i2][0]))
plt.plot(segment[:,0],segment[:,1],color="black")
plt.axis("equal")
plt.show()
aligned = np.array([[1,2,3],[1,4,7],[1,5,9],[2,5,8],[3,6,9],[3,5],[4,5,6],[4,8],[7,8,9]],dtype=object)
for a in aligned:
p = []
for index in a:
p += [pts[indices == index][0]]
p = np.array(p)
print(p)
plt.plot(p[:, 0], p[:, 1])
plt.axis("equal")
plt.show()
注意：这个问题的表述方式与我之前在此处提出的问题不同且更普遍 我希望它不会被视为重复问题。
根据这些一般假设，不一定所有索引都存在。
我会尝试这样的算法：(python/pseudocode)
point_dict = dictionary( index : point )
seg_set = set( seg )
seg_endpoints = dictionary( index of point : list of segments containing this point )
result = []
while not empty(seg_set):
segment = seg_set.pop()
result.append(segment)
left, right = segment.left, segment.right
# going to the left
while True:
possible_lefts = seg_endpoints[left]
additional_left = get_line_to_left(possible_lefts)
if additional_left is None:
break
remove_additional_left_from_seg_set()
remove_additional_left_from_seg_endpoints()
left = additional_left.left
results[-1].prepend(left) # not real syntax offcource
# going to the right
while True:
possible_rights = seg_endpoints[left]
additional_right = get_line_to_right(possible_rights)
if additional_right is None:
break
remove_additional_right_from_seg_set()
remove_additional_right_from_seg_endpoints()
left = additional_right.right
results[-1].append(right) # this is real syntax :)
我遗漏了很多实现细节，但是既然你提到你尝试了不同的方法但没有取得太大的成功，我希望你至少可以从这个算法中得到一些想法。
一个重要的实现细节是如何将段存储在 seg_endpoints 中以及如何从该段中删除特定段。
一种可能性是将段转换为数据类，然后您可以 link 转换为这些段。另一个是用一些额外的索引来跟踪所有这些段的位置。
无论如何，对于数万个点的网络，我当然认为制作这些字典会为你的算法节省很多时间。
祝你好运！
这是一种可能的方法。
为了有一些测试用例，让我们从生成问题中描述的类型的随机图的代码开始。函数generate_graph() returns 一个元组(points, edges)
其中 points 是具有图顶点坐标的元组列表，edges 是具有通过边连接的顶点索引的元组列表。图顶点在具有 nrow 行和 ncol 列的网格中均匀分布。参数 p 给出了连接两个相邻顶点的边被创建的概率。 p 的值越高，图形的边数就越多。 distortion 参数可用于扭曲顶点网格。
from itertools import product
import numpy as np
def generate_graph(nrow, ncol, p=0.5, distortion=0):
X, Y = np.meshgrid(np.arange(ncol), np.arange(nrow))
if distortion:
X = X + np.random.random(X.shape)*distortion
Y = Y + np.random.random(Y.shape)*distortion
pts = list(zip(X.ravel(), Y.ravel()))
edges = []
for i, j in product(range(nrow-1), range(ncol-1)):
rands = np.random.random(4)
if rands[0] < p:
edges.append(np.ravel_multi_index([[i, i+1], [j, j]], (nrow, ncol)))
if rands[1] < p:
edges.append(np.ravel_multi_index([[i, i], [j, j+1]], (nrow, ncol)))
if rands[2] < 1 - (1 - p)**2:
if rands[3] > 0.5:
edges.append(np.ravel_multi_index([[i, i+1], [j, j+1]], (nrow, ncol)))
else:
edges.append(np.ravel_multi_index([[i, i+1], [j+1, j]], (nrow, ncol)))
edges = [tuple(t) for t in edges]
return pts, edges
这是一个示例输出：
np.random.seed(333)
pts, edges = generate_graph(nrow=5, ncol=10, p=0.7)
print(f"sample vertex coordinates: {pts[:5]}")
print(f"sample edges: {edges[:5]}")
这给出：
sample vertex coordinates: [(0, 0), (1, 0), (2, 0), (3, 0), (4, 0)]
sample edges: [(0, 10), (1, 10), (1, 11), (1, 2), (2, 11)]
请注意，与问题的陈述不同，顶点坐标不包括顶点索引，因为顶点已经根据它们在 pts 列表中的位置进行了索引。
下一个函数绘制结果图：
import matplotlib.pyplot as plt
def plot_graph(pts, edges):
fs = 10
plt.figure(figsize=(fs, fs))
plt.axes().set_aspect('equal')
for s in edges:
plt.plot(*zip(pts[s[0]], pts[s[1]]), 'g-', alpha=0.5)
plt.scatter(*zip(*pts), s=50, c='r', marker='o', zorder=10)
for i, p in enumerate(pts):
t = plt.text(*(np.array(p) + 0.1), str(i))
plt.axis('off')
plt.show()
对于上面创建的图表，我们得到：
plot_graph(pts, edges)
下一个函数搜索与顶点相邻的共线边。对于给定的图顶点，函数 colinear() 找到所有以该顶点作为其端点之一的边，并且 returns 是共线的边对列表。搜索共线对是通过首先计算每条边与正 x 轴之间的角度，然后比较任意两条边的角度来完成的。共线边对将具有相等或相隔 180 度的角度。该函数采用三个参数：center 是要考虑的图顶点的索引，pts 是包含图顶点的列表。seg 最初是图边的列表 - 但更多关于它下面。
def colinear(center, pts, segs, tol=10**(-9)):
edges = [i for i in range(len(segs)) if center in [segs[i][0], segs[i][-1]]]
if len(edges) < 2:
return []
vertices = [segs[i][1] if segs[i].index(center) == 0 else segs[i][-2] for i in edges]
vects = np.array([pts[i] for i in vertices]) - np.array(pts[center])
vects = vects[:, 0] + vects[:, 1] * 1j
angles = np.angle(vects)
dirs = np.abs(angles - angles.reshape(-1, 1))
dirs = np.minimum(dirs, np.abs(dirs - np.pi))
pairs = np.array(edges)[np.argwhere(np.tril(dirs < tol, k=-1))]
pairs = [[segs[i][:], segs[j][:]] for [i, j] in pairs]
return pairs
例如，我们可以找到与顶点编号 26 相邻的共线边：
colinear(26, pts, edges)
这给出：
[[(26, 35), (17, 26)], [(26, 36), (16, 26)]]
这意味着连接顶点 26 到顶点 35 和 17 的边是共线的，同样对于连接顶点 26 到顶点 36 和 16 的边。看上面的图可以确认这些都是共线边对在这种情况下。
我们需要的最后一个函数是concat()。对于共享一个端点的给定路径对（即由边连接的图顶点序列），它会生成一条新路径，将两条路径在公共端点处连接在一起。
def concat(pair):
center = (set(pair[0]) & set(pair[1])).pop()
new_pair = []
for i in range(2):
if pair[i].index(center) != 0:
new_pair.append(pair[i][::-1])
else:
new_pair.append(pair[i][:])
bds = [new_pair[i][-1] for i in range(2)]
new_seg = new_pair[0][::-1] + new_pair[1][1:]
if new_seg[0] > new_seg[1]:
new_seg = new_seg[::-1]
return new_seg
我们现在可以找到图中所有由共线边组成的路径。路径将存储在列表 segs 中。一开始，这个列表由图的所有边组成。然后，对于每个图形顶点，我们找到与该顶点相邻的共线路径对，并将每一对连接起来形成更长的路径：
segs = [sorted(s) for s in edges]
for i in range(len(pts)):
colins = colinear(i, pts, segs)
for pair in colins:
new_seg = concat(pair)
for i in range(2):
segs.remove(pair[i])
segs.append(new_seg)
最后segs包含所有由共线边组成的最长路径。这是一个打印输出，显示了至少包含 5 个顶点的所有此类路径：
for s in segs:
if len(s) >= 5:
print(s)
它给出：
[3, 4, 5, 6, 7]
[14, 15, 16, 17, 18]
[0, 10, 20, 30, 40]
[1, 11, 21, 31, 41]
[8, 17, 26, 35, 44]
[6, 16, 26, 36, 46]
[8, 18, 28, 38, 48]
可以将这些结果与上图进行比较。
添加失真
接下来，让我们创建一个顶点网格在一定程度上扭曲的图形：
np.random.seed(200)
pts, edges = generate_graph(nrow=5, ncol=10, p=0.8, distortion=0.2)
plot_graph(pts, edges)
函数 colinear() 有一个 tol 参数，可用于为图边的共线性计算增加一些容差。但是，当前编写的代码仅在 tol 与边之间的角度相比较小时才有效。更改此设置需要进行一些修改。无论如何，让我们将其应用于扭曲图，将公差设置为 5 度：
tol = 5*(2*np.pi/360)
segs = [sorted(s) for s in edges]
for i in range(len(pts)):
colins = colinear(i, pts, segs, tol=tol)
for pair in colins:
new_seg = concat(pair)
for i in range(2):
segs.remove(pair[i])
segs.append(new_seg)
我们可以打印长度至少为 4 的近似共线路径：
for s in segs:
if len(s) >= 4:
print(s)
它给出：
[0, 10, 20, 30]
[4, 13, 22, 31]
[30, 31, 32, 33]
[35, 36, 37, 38]
最后说明。以上代码可以优化。特别是，使用字典记录每个顶点是哪些路径的端点会更有效，这样函数 colinear() 就不需要在每次调用时都搜索这些路径。
我有一个结构为 (x, y, index)
示例：
pts = [[0,0,1],[0,1,2],[0,2,3],[1,0,4],[1,1,5],[1,2,6],[2,0,7],[2,1,8],[2,2,9]]
和一系列具有结构 (index1, index2) 的线段，将这些点 2 乘 2 连接起来：
segs = [[1,2],[1,4],[1,5],[2,3],[2,5],[3,5],[3,6],[4,5],[4,7],[4,8],[5,6],[5,8],[5,9],[6,9],[7,8],[8,9]]
线段总是只连接两个点，它们永远不会经过第三个点。
当我应该只使用 seg 中的段一次时，如何检测对齐的段并使用 n 索引制作更长的段？
aligned = [[1,2,3],[1,4,7],[1,5,9],[2,5,8],[3,6,9],[3,5],[4,5,6],[4,8],[7,8,9]]
现在说点没有完全对齐：
pts = [[0.002,0.1,1],[0.0003,1.003,2],[-0.01,2.11,3],[1.01,0.101,4],[1.0005,1.2,5],[1.25,2.01,6],[2.007,-0.12,7],[1.996,1.1,8],[2.03,2.1,9]]
如果线段之间的角度几乎为零，我想连接几乎对齐的线段，由 tolerance 定义，例如，3°
我试过很多方法都没有成功。这必须在 10k+ 点的大量网络上完成。
此代码重现了我的示例：
import numpy as np
import matplotlib.pyplot as plt
#pts = np.array([[0,0,1],[0,1,2],[0,2,3],[1,0,4],[1,1,5],[1,2,6],[2,0,7],[2,1,8],[2,2,9]])
pts = np.array([[0.002,0.1,1],[0.0003,1.003,2],[-0.01,2.11,3],[1.01,0.101,4],[1.0005,1.2,5],[1.25,2.01,6],[2.007,-0.12,7],[1.996,1.1,8],[2.03,2.1,9]])
indices = pts[:,2].astype(int)
seg = np.array([[1,2],[1,4],[1,5],[2,3],[2,5],[3,5],[3,6],[4,5],[4,7],[4,8],[5,6],[5,8],[5,9],[6,9],[7,8],[8,9]])
for s in seg:
i1, i2 = s
segment = np.array((pts[indices == i1][0], pts[indices == i2][0]))
plt.plot(segment[:,0],segment[:,1],color="black")
plt.axis("equal")
plt.show()
aligned = np.array([[1,2,3],[1,4,7],[1,5,9],[2,5,8],[3,6,9],[3,5],[4,5,6],[4,8],[7,8,9]],dtype=object)
for a in aligned:
p = []
for index in a:
p += [pts[indices == index][0]]
p = np.array(p)
print(p)
plt.plot(p[:, 0], p[:, 1])
plt.axis("equal")
plt.show()
注意：这个问题的表述方式与我之前在此处提出的问题不同且更普遍
根据这些一般假设，不一定所有索引都存在。
我会尝试这样的算法：(python/pseudocode)
point_dict = dictionary( index : point )
seg_set = set( seg )
seg_endpoints = dictionary( index of point : list of segments containing this point )
result = []
while not empty(seg_set):
segment = seg_set.pop()
result.append(segment)
left, right = segment.left, segment.right
# going to the left
while True:
possible_lefts = seg_endpoints[left]
additional_left = get_line_to_left(possible_lefts)
if additional_left is None:
break
remove_additional_left_from_seg_set()
remove_additional_left_from_seg_endpoints()
left = additional_left.left
results[-1].prepend(left) # not real syntax offcource
# going to the right
while True:
possible_rights = seg_endpoints[left]
additional_right = get_line_to_right(possible_rights)
if additional_right is None:
break
remove_additional_right_from_seg_set()
remove_additional_right_from_seg_endpoints()
left = additional_right.right
results[-1].append(right) # this is real syntax :)
我遗漏了很多实现细节，但是既然你提到你尝试了不同的方法但没有取得太大的成功，我希望你至少可以从这个算法中得到一些想法。 一个重要的实现细节是如何将段存储在 seg_endpoints 中以及如何从该段中删除特定段。
一种可能性是将段转换为数据类，然后您可以 link 转换为这些段。另一个是用一些额外的索引来跟踪所有这些段的位置。
无论如何，对于数万个点的网络，我当然认为制作这些字典会为你的算法节省很多时间。
祝你好运！
这是一种可能的方法。
为了有一些测试用例，让我们从生成问题中描述的类型的随机图的代码开始。函数generate_graph() returns 一个元组(points, edges)
其中 points 是具有图顶点坐标的元组列表，edges 是具有通过边连接的顶点索引的元组列表。图顶点在具有 nrow 行和 ncol 列的网格中均匀分布。参数 p 给出了连接两个相邻顶点的边被创建的概率。 p 的值越高，图形的边数就越多。 distortion 参数可用于扭曲顶点网格。
from itertools import product
import numpy as np
def generate_graph(nrow, ncol, p=0.5, distortion=0):
X, Y = np.meshgrid(np.arange(ncol), np.arange(nrow))
if distortion:
X = X + np.random.random(X.shape)*distortion
Y = Y + np.random.random(Y.shape)*distortion
pts = list(zip(X.ravel(), Y.ravel()))
edges = []
for i, j in product(range(nrow-1), range(ncol-1)):
rands = np.random.random(4)
if rands[0] < p:
edges.append(np.ravel_multi_index([[i, i+1], [j, j]], (nrow, ncol)))
if rands[1] < p:
edges.append(np.ravel_multi_index([[i, i], [j, j+1]], (nrow, ncol)))
if rands[2] < 1 - (1 - p)**2:
if rands[3] > 0.5:
edges.append(np.ravel_multi_index([[i, i+1], [j, j+1]], (nrow, ncol)))
else:
edges.append(np.ravel_multi_index([[i, i+1], [j+1, j]], (nrow, ncol)))
edges = [tuple(t) for t in edges]
return pts, edges
这是一个示例输出：
np.random.seed(333)
pts, edges = generate_graph(nrow=5, ncol=10, p=0.7)
print(f"sample vertex coordinates: {pts[:5]}")
print(f"sample edges: {edges[:5]}")
这给出：
sample vertex coordinates: [(0, 0), (1, 0), (2, 0), (3, 0), (4, 0)]
sample edges: [(0, 10), (1, 10), (1, 11), (1, 2), (2, 11)]
请注意，与问题的陈述不同，顶点坐标不包括顶点索引，因为顶点已经根据它们在 pts 列表中的位置进行了索引。
下一个函数绘制结果图：
import matplotlib.pyplot as plt
def plot_graph(pts, edges):
fs = 10
plt.figure(figsize=(fs, fs))
plt.axes().set_aspect('equal')
for s in edges:
plt.plot(*zip(pts[s[0]], pts[s[1]]), 'g-', alpha=0.5)
plt.scatter(*zip(*pts), s=50, c='r', marker='o', zorder=10)
for i, p in enumerate(pts):
t = plt.text(*(np.array(p) + 0.1), str(i))
plt.axis('off')
plt.show()
对于上面创建的图表，我们得到：
plot_graph(pts, edges)
下一个函数搜索与顶点相邻的共线边。对于给定的图顶点，函数 colinear() 找到所有以该顶点作为其端点之一的边，并且 returns 是共线的边对列表。搜索共线对是通过首先计算每条边与正 x 轴之间的角度，然后比较任意两条边的角度来完成的。共线边对将具有相等或相隔 180 度的角度。该函数采用三个参数：center 是要考虑的图顶点的索引，pts 是包含图顶点的列表。seg 最初是图边的列表 - 但更多关于它下面。
def colinear(center, pts, segs, tol=10**(-9)):
edges = [i for i in range(len(segs)) if center in [segs[i][0], segs[i][-1]]]
if len(edges) < 2:
return []
vertices = [segs[i][1] if segs[i].index(center) == 0 else segs[i][-2] for i in edges]
vects = np.array([pts[i] for i in vertices]) - np.array(pts[center])
vects = vects[:, 0] + vects[:, 1] * 1j
angles = np.angle(vects)
dirs = np.abs(angles - angles.reshape(-1, 1))
dirs = np.minimum(dirs, np.abs(dirs - np.pi))
pairs = np.array(edges)[np.argwhere(np.tril(dirs < tol, k=-1))]
pairs = [[segs[i][:], segs[j][:]] for [i, j] in pairs]
return pairs
例如，我们可以找到与顶点编号 26 相邻的共线边：
colinear(26, pts, edges)
这给出：
[[(26, 35), (17, 26)], [(26, 36), (16, 26)]]
这意味着连接顶点 26 到顶点 35 和 17 的边是共线的，同样对于连接顶点 26 到顶点 36 和 16 的边。看上面的图可以确认这些都是共线边对在这种情况下。
我们需要的最后一个函数是concat()。对于共享一个端点的给定路径对（即由边连接的图顶点序列），它会生成一条新路径，将两条路径在公共端点处连接在一起。
def concat(pair):
center = (set(pair[0]) & set(pair[1])).pop()
new_pair = []
for i in range(2):
if pair[i].index(center) != 0:
new_pair.append(pair[i][::-1])
else:
new_pair.append(pair[i][:])
bds = [new_pair[i][-1] for i in range(2)]
new_seg = new_pair[0][::-1] + new_pair[1][1:]
if new_seg[0] > new_seg[1]:
new_seg = new_seg[::-1]
return new_seg
我们现在可以找到图中所有由共线边组成的路径。路径将存储在列表 segs 中。一开始，这个列表由图的所有边组成。然后，对于每个图形顶点，我们找到与该顶点相邻的共线路径对，并将每一对连接起来形成更长的路径：
segs = [sorted(s) for s in edges]
for i in range(len(pts)):
colins = colinear(i, pts, segs)
for pair in colins:
new_seg = concat(pair)
for i in range(2):
segs.remove(pair[i])
segs.append(new_seg)
最后segs包含所有由共线边组成的最长路径。这是一个打印输出，显示了至少包含 5 个顶点的所有此类路径：
for s in segs:
if len(s) >= 5:
print(s)
它给出：
[3, 4, 5, 6, 7]
[14, 15, 16, 17, 18]
[0, 10, 20, 30, 40]
[1, 11, 21, 31, 41]
[8, 17, 26, 35, 44]
[6, 16, 26, 36, 46]
[8, 18, 28, 38, 48]
可以将这些结果与上图进行比较。
添加失真
接下来，让我们创建一个顶点网格在一定程度上扭曲的图形：
np.random.seed(200)
pts, edges = generate_graph(nrow=5, ncol=10, p=0.8, distortion=0.2)
plot_graph(pts, edges)
函数 colinear() 有一个 tol 参数，可用于为图边的共线性计算增加一些容差。但是，当前编写的代码仅在 tol 与边之间的角度相比较小时才有效。更改此设置需要进行一些修改。无论如何，让我们将其应用于扭曲图，将公差设置为 5 度：
tol = 5*(2*np.pi/360)
segs = [sorted(s) for s in edges]
for i in range(len(pts)):
colins = colinear(i, pts, segs, tol=tol)
for pair in colins:
new_seg = concat(pair)
for i in range(2):
segs.remove(pair[i])
segs.append(new_seg)
我们可以打印长度至少为 4 的近似共线路径：
for s in segs:
if len(s) >= 4:
print(s)
它给出：
[0, 10, 20, 30]
[4, 13, 22, 31]
[30, 31, 32, 33]
[35, 36, 37, 38]
最后说明。以上代码可以优化。特别是，使用字典记录每个顶点是哪些路径的端点会更有效，这样函数 colinear() 就不需要在每次调用时都搜索这些路径。