如何在图形中找到无弦环?
输入的数据是一个包含两维点的数组,这些点组成了多条线:
coordinates = [(x, y, x1, y1, x2, y2, etc.), (x3, y3, x4, y4, etc.)].
这里是所有的坐标:https://pastebin.com/raw/ascwXxKY。
任务是找出图形中所有没有弦的循环。你需要注意,弦不仅可以在图的点之间画,还可以,比如说,从一条线的中间到另一条线的中间。
这个循环可以不接触其他的线。
你还要记得,弦可以由多条线组成,就像下面的图所示。
我的图中总共有15个这样的循环:
我尝试使用networkx中的无弦循环方法,但到目前为止还没有找到适合我情况的方法。
下面是一段代码,可以帮助绘制我上面指定的坐标的所有线:
# Create coordinate lists for each pair of points
x_coords = []
y_coords = []
all_connections = []
for coord_set in coordinates:
x_pair = []
y_pair = []
connections_pair = []
for i in range(len(coord_set) // 2 - 1):
x1, y1 = coord_set[2*i], coord_set[2*i+1]
x2, y2 = coord_set[2*i+2], coord_set[2*i+3]
x_pair.extend([x1, x2])
y_pair.extend([y1, y2])
connections_pair.append(((x1, y1), (x2, y2)))
all_connections.append(((x1, y1), (x2, y2)))
x_coords.append(x_pair)
y_coords.append(y_pair)
# Plot the graph
for i in range(len(x_coords)):
plt.plot(x_coords[i], y_coords[i], marker='o')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Line Chart')
plt.grid(True)
plt.show()
3 个回答
假设输入中已经给出了每一条线交点的坐标:
- 找到所有的多边形
- Construct undirected graph.
( Vertices are specified in the input. Edges are lines between vertices. )
- Construct empty list of polygons
- Loop over every vertex specified in the input
- Check if it is a member of a polygon
( is there a path from the vertex that loops back again to the starting vertex.
Any decent graph theory library will give you a method to do this. )
- If polygon is unique, add to list of polygons
- 找到空的多边形
- Loop over vertices
- Loop over list of polygons
- IF vertex NOT part of polygon
- IF vertex INSIDE polygon
( algorithm https://ravenspoint.wordpress.com/2010/06/27/in-or-out/ )
- REMOVE polygon from list of polygons
( careful not to invalidate loop iterator )
- 输出空多边形的列表
注意,如果线可以交叉而输入中没有给出坐标,那么你需要检查属于同一个多边形的每一对线是否相交(任何计算几何的书籍都会提供这个算法),并且要去掉那些重叠的多边形。
这是我对你之前提问的同一主题的回答的改编版本。你需要将线条“对齐”到合适的点上,以处理你数据中已经存在的交点误差。要注意的是,如果你没有现成的交点,你需要手动找到这些点并把它们加到线条上。这个解决方案并不会寻找交点,而是假设这些点已经存在,除了浮点数误差之外。
import numpy as np
import networkx as nx
from shapely import geometry, ops, union_all, intersection
eps = 1.e-8
def add_intersections(coords):
all_points = [geometry.Point(pair) for line in coords
for pair in zip(line[::2],line[1::2])]
point_lists = [[geometry.Point(pair) for pair in zip(line[::2],line[1::2])]
for line in coords]
lines = [geometry.LineString(line) for line in point_lists]
for line_no in range(len(lines)):
for point in all_points:
lines[line_no] = ops.snap(lines[line_no], point, eps)
return [[val for coord in line.coords for val in coord] for line in lines]
def get_cycles(coords):
point_lists = [[pair for pair in zip(line[::2],line[1::2])] for line in coords]
edges = [segment for point_list in point_lists
for segment in zip(point_list[:-1], point_list[1:])]
G = nx.from_edgelist(edges)
return nx.chordless_cycles(G)
def is_elementary(cycle, cycles):
union = union_all(cycles)
intersection_ = intersection(union, cycle)
if intersection_.area > eps:
return False
else:
return True
def elementary_filter(cycles):
cycles = [geometry.Polygon(cycle) for cycle in cycles]
cycles.sort(key=lambda polygon: polygon.area)
elementary_cycles = [cycles[0],]
elementary_indices = [0,]
for i in range(1, len(cycles)):
if is_elementary(cycles[i], elementary_cycles):
elementary_cycles.append(cycles[i])
elementary_indices.append(i)
return elementary_cycles, elementary_indices
def get_elementary_cycles(coords):
coords = add_intersections(coords)
cycles = get_cycles(coords)
elementary_polygons, indices = elementary_filter(cycles)
return elementary_polygons
# Make sure coords is already defined
basis_cycles = get_elementary_cycles(coords)
print('The number of basis cycles:', len(basis_cycles))
如果你不想处理Polygon
对象,可以在get_elementary_cycles
中使用indices
来索引原始列表。我还添加了一些简单的函数来绘制线条和多边形:
import matplotlib.pyplot as plt
from itertools import cycle
def plot_lines(coords):
line_array = [np.array([pair for pair in zip(line[::2],line[1::2])]).T
for line in coords]
marker_list = ['o', 's', 'D', '^', 'v', '*', 'H', 'P', 'd']
markers = cycle(marker_list)
fig, ax = plt.subplots()
for line in line_array:
ax.plot(line[0], line[1], marker=next(markers), mfc='none', alpha=0.5)
ax.grid()
ax.set_title('Lines')
return fig
def plot_polygons(polygons):
fig, ax = plt.subplots()
for polygon in polygons:
x, y = polygon.exterior.xy
ax.plot(x, y, alpha=0.5)
ax.set_title('Polygons')
return fig
from matplotlib.collections import PatchCollection
from shapely.plotting import patch_from_polygon
colors = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
def plot_filled_polygons(polygons):
fig, ax = plt.subplots()
patches = [patch_from_polygon(polygon) for polygon in polygons]
color = [next(colors) for patch in patches]
patches = PatchCollection(patches, alpha=0.5, fc=color)
patches.set(alpha=0.5)
ax.add_collection(patches)
ax.set_xlim([85, 210])
ax.set_ylim([55, 110])
ax.set_title('Filled polygons')
return fig
fig1 = plot_lines(coords)
fig2 = plot_polygons(basis_cycles)
fig3 = plot_filled_polygons(basis_cycles)
plt.show()
因为数据集中有一些边相交(除了在它们的端点),所以你首先需要把它整理成一个没有这些交叉的图,也就是说,把这些边分开。
在你给出的某些例子中,原本应该相接的边实际上并没有相接,这是因为输入数据的精度有限。为了应对这个问题,你可以把那些“死胡同”的边稍微调整一下,让它们变得长一点。
接下来,你可以通过总是选择最左转(或最右转)的邻近边来进行操作。我建议先复制这些边,这样它们就可以双向存在。一旦访问了某个方向的边,就可以把它从图中移除。
会有一条路径沿着图的外部走。这条路径可以排除,因为它是以相反的方向完成循环的(相对角度的总和与其他真实的区域路径不同)。
这里有一段不使用图形库的代码:
from collections import defaultdict
coords = [(2536.948596571963, 1495.6001088264184, 2537.706596571963, 1491.4401088264185, 2538.720596571963, 1485.6001088264184, 2539.1375965719635, 1483.5001088264185, 2540.6675965719633, 1475.6001088264184, 2541.3665965719633, 1471.9401088264185, 2542.566596571963, 1465.6001088264184, 2543.209596571963, 1462.5601088264184, 2544.727596571963, 1455.6001088264184, 2546.9675965719634, 1446.4701088264185, 2547.1915965719636, 1445.6001088264184, 2549.0705965719635, 1436.4801088264185, 2549.2515965719635, 1435.6001088264184, 2549.3065965719634, 1435.6001088264184, 2549.3115965719635, 1435.6001088264184, 2552.0625965719632, 1435.6001088264184, 2552.3485965719633, 1435.6001088264184, 2554.729596571963, 1435.6001088264184, 2555.294596571963, 1435.6001088264184, 2557.3605965719635, 1435.6001088264184, 2559.0515965719633, 1435.6001088264184, 2559.9605965719634, 1435.6001088264184, 2561.3875965719635, 1435.6001088264184, 2562.5125965719635, 1435.6001088264184, 2563.7235965719633, 1435.6001088264184, 2564.8165965719636, 1435.6001088264184, 2566.1095965719633, 1435.6001088264184, 2567.5525965719635, 1435.6001088264184, 2568.874596571963, 1435.6001088264184, 2570.901596571963, 1435.6001088264184, 2572.4065965719633, 1435.6001088264184, 2573.961596571963, 1435.6001088264184, 2575.3975965719633, 1435.6001088264184, 2576.2905965719633, 1435.6001088264184, 2578.081596571963, 1435.6001088264184, 2579.1065965719636, 1435.6001088264184, 2579.2315965719636, 1435.6001088264184, 2580.4365965719635, 1435.6001088264184, 2580.4465965719633, 1435.6001088264184, 2580.545596571963, 1435.6001088264184, 2580.8465965719633, 1435.6001088264184, 2581.147596571963, 1435.6001088264184, 2581.440596571963, 1435.6001088264184, 2581.723596571963, 1435.6001088264184, 2581.913596571963, 1435.6001088264184, 2582.232596571963, 1435.6001088264184, 2582.502596571963, 1435.6001088264184, 2582.857596571963, 1435.6001088264184, 2583.247596571963, 1435.6001088264184, 2583.612596571963, 1435.6001088264184, 2584.059596571963, 1435.6001088264184, 2584.4775965719627, 1435.6001088264184, 2584.848596571963, 1435.6001088264184, 2585.337596571963, 1435.6001088264184, 2585.7265965719625, 1435.6001088264184, 2586.103596571963, 1435.6001088264184, 2586.4585965719625, 1435.6001088264184, 2587.013596571963, 1435.6001088264184, 2587.371596571963, 1435.6001088264184, 2587.684596571963, 1435.6001088264184, 2587.946596571963, 1435.6001088264184, 2588.185596571963, 1435.6001088264184, 2588.770596571963, 1435.6001088264184, 2589.032596571963, 1435.6001088264184, 2592.299596571963, 1435.6001088264184, 2592.6515965719627, 1435.6001088264184, 2592.938596571963, 1435.6001088264184, 2593.640596571963, 1435.6001088264184, 2593.805596571963, 1435.6001088264184, 2593.962596571963, 1435.6001088264184, 2594.5305965719626, 1435.6001088264184, 2594.683596571963, 1435.6001088264184, 2595.199596571963, 1435.6001088264184, 2595.330596571963, 1435.6001088264184, 2595.876596571963, 1435.6001088264184, 2595.988596571963, 1435.6001088264184, 2596.5435965719635, 1435.6001088264184, 2596.6505965719634, 1435.6001088264184, 2597.1995965719634, 1435.6001088264184, 2597.304596571963, 1435.6001088264184, 2597.8335965719634, 1435.6001088264184, 2598.327596571963, 1435.6001088264184, 2598.433596571963, 1435.6001088264184, 2598.8915965719634, 1435.6001088264184, 2599.298596571963, 1435.6001088264184, 2599.4085965719632, 1435.6001088264184, 2599.769596571963, 1435.6001088264184, 2600.071596571963, 1435.6001088264184, 2600.186596571963, 1435.6001088264184, 2600.431596571963, 1435.6001088264184, 2612.4455965719626, 1440.3401088264186, 2612.5935965719627, 1440.3401088264186, 2613.1735965719627, 1440.4201088264185, 2613.7825965719626, 1440.5001088264185, 2613.9165965719626, 1440.5001088264185, 2614.5255965719625, 1440.5801088264186, 2614.6675965719623, 1440.5801088264186, 2615.2515965719626, 1440.6601088264188, 2635.4905965719627, 1442.8301088264186, 2635.9525965719627, 1442.8901088264186, 2636.3505965719623, 1442.9401088264185, 2636.6685965719626, 1442.9801088264185, 2636.897596571963, 1443.0101088264184, 2637.0205965719624, 1443.0301088264184, 2637.1845965719626, 1443.0501088264184, 2637.3545965719627, 1443.0701088264184, 2637.4935965719624, 1443.0901088264184, 2637.5665965719627, 1443.1001088264184, 2647.2875965719627, 1445.6001088264184, 2647.322596571963, 1445.6001088264184, 2647.366596571963, 1445.6001088264184, 2647.4235965719627, 1445.6001088264184, 2647.484596571963, 1445.6001088264184, 2647.563596571963, 1445.6001088264184, 2650.223596571963, 1445.6001088264184, 2650.350596571963, 1445.6001088264184, 2650.531596571963, 1445.6001088264184, 2654.813596571963, 1445.6001088264184, 2654.870596571963, 1445.6001088264184, 2654.9325965719627, 1445.6001088264184, 2655.116596571963, 1445.6001088264184, 2655.2675965719627, 1445.6001088264184, 2655.830596571963, 1445.6001088264184, 2657.420596571963, 1445.6001088264184, 2666.8145965719627, 1449.8101088264184, 2670.3905965719628, 1451.9401088264185, 2673.974596571963, 1453.8601088264184, 2676.533596571963, 1455.6001088264184, 2676.675596571963, 1455.7901088264184, 2676.857596571963, 1456.0301088264184), (2548.531751224003, 1435.0954731549964, 2559.114905876043, 1431.6001088264184), (2548.531751224003, 1431.0954731549964, 2559.114905876043, 1431.6001088264184), (2548.531751224003, 1423.0954731549964, 2559.114905876043, 1419.6001088264184), (2548.531751224003, 1419.0954731549964, 2559.114905876043, 1415.6001088264184), (2548.531751224003, 1411.0954731549964, 2559.114905876043, 1407.6001088264184), (2559.114905876043, 1431.6001088264184, 2567.984493888275, 1431.6001088264184), (2559.114905876043, 1431.6001088264184, 2567.984493888275, 1423.6001088264184), (2559.114905876043, 1419.6001088264184, 2567.984493888275, 1419.6001088264184), (2559.114905876043, 1415.6001088264184, 2567.984493888275, 1415.6001088264184), (2559.114905876043, 1407.6001088264184, 2567.984493888275, 1403.6001088264184), (2567.984493888275, 1431.6001088264184, 2577.4785637051095, 1431.6001088264184), (2567.984493888275, 1423.6001088264184, 2577.4785637051095, 1423.6001088264184), (2567.984493888275, 1419.6001088264184, 2577.4785637051095, 1415.6001088264184), (2567.984493888275, 1415.6001088264184, 2577.4785637051095, 1411.6001088264184), (2567.984493888275, 1403.6001088264184, 2577.4785637051095, 1403.6001088264184), (2577.4785637051095, 1431.6001088264184, 2588.2013160893976, 1431.6001088264184), (2577.4785637051095, 1423.6001088264184, 2588.2013160893976, 1423.6001088264184), (2577.4785637051095, 1423.6001088264184, 2588.2013160893976, 1415.6001088264184), (2577.4785637051095, 1415.6001088264184, 2588.2013160893976, 1407.6001088264184), (2577.4785637051095, 1411.6001088264184, 2588.2013160893976, 1407.6001088264184), (2577.4785637051095, 1403.6001088264184, 2588.2013160893976, 1399.6001088264184), (2588.2013160893976, 1431.6001088264184, 2597.798842058678, 1431.6001088264184), (2588.2013160893976, 1423.6001088264184, 2597.798842058678, 1427.6001088264184), (2588.2013160893976, 1415.6001088264184, 2597.798842058678, 1423.6001088264184), (2588.2013160893976, 1407.6001088264184, 2597.798842058678, 1407.6001088264184), (2588.2013160893976, 1407.6001088264184, 2597.798842058678, 1403.6001088264184), (2588.2013160893976, 1399.6001088264184, 2597.798842058678, 1395.6001088264184), (2597.798842058678, 1431.6001088264184, 2607.842362905752, 1434.5239503797861), (2597.798842058678, 1427.6001088264184, 2607.842362905752, 1426.5239503797861), (2597.798842058678, 1423.6001088264184, 2607.842362905752, 1426.5239503797861), (2597.798842058678, 1407.6001088264184, 2607.842362905752, 1410.5239503797861), (2597.798842058678, 1403.6001088264184, 2607.842362905752, 1406.5239503797861), (2597.798842058678, 1395.6001088264184, 2607.842362905752, 1386.5239503797861), (2607.842362905752, 1434.5239503797861, 2617.0717529170083, 1436.8552636892455), (2607.842362905752, 1426.5239503797861, 2617.0717529170083, 1432.8552636892455), (2607.842362905752, 1426.5239503797861, 2617.0717529170083, 1424.8552636892455), (2607.842362905752, 1410.5239503797861, 2617.0717529170083, 1404.8552636892455), (2607.842362905752, 1406.5239503797861, 2617.0717529170083, 1400.8552636892455), (2607.842362905752, 1386.5239503797861, 2617.0717529170083, 1388.8552636892455), (2617.0717529170083, 1436.8552636892455, 2627.8066294831897, 1438.006243586899), (2617.0717529170083, 1432.8552636892455, 2627.8066294831897, 1434.006243586899), (2617.0717529170083, 1424.8552636892455, 2627.8066294831897, 1430.006243586899), (2617.0717529170083, 1404.8552636892455, 2627.8066294831897, 1410.006243586899), (2617.0717529170083, 1400.8552636892455, 2627.8066294831897, 1410.006243586899), (2617.0717529170083, 1388.8552636892455, 2627.8066294831897, 1402.006243586899), (2627.8066294831897, 1438.006243586899, 2637.720881676667, 1439.1397871271859), (2627.8066294831897, 1434.006243586899, 2637.720881676667, 1427.1397871271859), (2627.8066294831897, 1430.006243586899, 2637.720881676667, 1423.1397871271859), (2627.8066294831897, 1410.006243586899, 2637.720881676667, 1411.1397871271859), (2627.8066294831897, 1402.006243586899, 2637.720881676667, 1403.1397871271859), (2637.720881676667, 1439.1397871271859, 2648.126410539704, 1441.6001088264184), (2637.720881676667, 1427.1397871271859, 2648.126410539704, 1433.6001088264184), (2637.720881676667, 1423.1397871271859, 2648.126410539704, 1429.6001088264184), (2637.720881676667, 1411.1397871271859, 2648.126410539704, 1393.6001088264184), (2637.720881676667, 1403.1397871271859, 2648.126410539704, 1385.6001088264184), (2648.126410539704, 1441.6001088264184, 2657.704484528568, 1441.7273355985399), (2648.126410539704, 1433.6001088264184, 2657.704484528568, 1429.7273355985399), (2648.126410539704, 1429.6001088264184, 2657.704484528568, 1421.7273355985399), (2648.126410539704, 1393.6001088264184, 2657.704484528568, 1393.7273355985399), (2648.126410539704, 1385.6001088264184, 2657.704484528568, 1377.7273355985399), (2657.704484528568, 1441.7273355985399, 2667.161761375854, 1446.0168932314207), (2657.704484528568, 1429.7273355985399, 2667.161761375854, 1430.0168932314207), (2657.704484528568, 1421.7273355985399, 2667.161761375854, 1426.0168932314207), (2657.704484528568, 1393.7273355985399, 2667.161761375854, 1398.0168932314207), (2657.704484528568, 1377.7273355985399, 2667.161761375854, 1390.0168932314207), (2548.531751224003, 1439.0954731549964, 2548.531751224003, 1435.0954731549964), (2548.531751224003, 1435.0954731549964, 2548.531751224003, 1431.0954731549964), (2548.531751224003, 1431.0954731549964, 2548.531751224003, 1423.0954731549964), (2548.531751224003, 1423.0954731549964, 2548.531751224003, 1419.0954731549964), (2548.531751224003, 1419.0954731549964, 2548.531751224003, 1411.0954731549964), (2667.161761375854, 1450.0168932314207, 2667.161761375854, 1446.0168932314207), (2667.161761375854, 1446.0168932314207, 2667.161761375854, 1430.0168932314207), (2667.161761375854, 1430.0168932314207, 2667.161761375854, 1426.0168932314207), (2667.161761375854, 1426.0168932314207, 2667.161761375854, 1398.0168932314207), (2667.161761375854, 1398.0168932314207, 2667.161761375854, 1390.0168932314207)]
import math
def get_edges(coords):
return [
((x1, y1), (x2, y2))
for polyline in coords
for x1, y1, x2, y2 in zip(polyline[::2], polyline[1::2], polyline[2::2], polyline[3::2])
]
def cut_intersecting_edges(edge1, edge2):
a, b = edge1
c, d = edge2
x1, y1 = a
x2, y2 = b
x3, y3 = c
x4, y4 = d
if (max(x1, x2) < min(x3, x4) or max(x3, x4) < min(x1, x2)
or max(y1, y2) < min(y3, y4) or max(y3, y4) < min(y1, y2)):
# Edges are outside eachother's bounding box
return
tnum = (x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4)
unum = (y1 - y2) * (x1 - x3) - (x1 - x2) * (y1 - y3)
denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4)
if (0 <= tnum <= denom or 0 >= tnum >= denom) and (0 <= unum <= denom or 0 >= unum >= denom):
if denom == 0:
print("overlap")
# overlap
polyline = sorted([*edge1, *edge2])
return list(zip(polyline, polyline[1:]))
# intersection
t = tnum / denom
p = (x1 + t * (x2 - x1), y1 + t * (y2 - y1))
print("split at", p)
return [(a, p),(p, b),(c, p),(p, d)]
def normalize_edges(edgeset):
edges = list(edgeset)
for i, edge1 in enumerate(edges):
if not edge1:
continue
a, b = edge1
if a == b: # A zero-length edge
edges[i] = None # Mark as deleted
continue
for j, edge2 in enumerate(edges):
if j >= i:
break
if not edge2 or a in edge2 or b in edge2:
# Consider adjacent edges as OK
continue
newedges = cut_intersecting_edges(edge1, edge2)
if newedges:
# Mark the intersecting edges as deleted
edgeset.remove(edge1)
edgeset.remove(edge2)
edges[i] = edges[j] = None
# ...and replace them with non-intersecting ones:
edges.extend(newedges)
edgeset.update(newedges)
break
def direct_edges(edges):
for i in range(len(edges)):
edges.append(edges[i][::-1])
def create_adjacency_list(edges):
adj = defaultdict(list)
for edge in edges:
adj[edge[0]].append(edge)
return adj
def normalize_angle(angle):
return (angle + math.tau) % math.tau
def get_areas(edges):
adj = create_adjacency_list(edges)
edgeset = set(edges)
while edgeset:
edge = next(iter(edgeset))
path = []
turn = 0
while edge in edgeset:
edgeset.remove(edge)
a, b = edge
angle = math.atan2(b[1] - a[1], b[0] - a[0])
path.append(b)
# Find index of next edge
min_angle = 2*math.pi
choice = 0
for i, (_, c) in enumerate(adj[b]):
if c == a: # Not interested in the back-edge
continue
angle2 = math.atan2(b[1] - c[1], b[0] - c[0])
diff = normalize_angle(angle2 - angle)
if diff < min_angle:
min_angle = diff
choice = i
# choose that edge
edge = adj[b][choice]
turn += min_angle - math.pi
if turn < 0: # Not an outside "area"
yield path
def get_dead_ends(adj):
for p, neighbors in adj.items():
if len(neighbors) == 1: # dead end
yield neighbors[0]
def lengthen_dead_ends(edgeset, grow_factor):
adj = create_adjacency_list(edgeset)
for short_edge in get_dead_ends(adj):
p, q = short_edge
# Make this edge a bit longer (data patching)
p = (p[0] + (p[0] - q[0]) * grow_factor,
p[1] + (p[1] - q[1]) * grow_factor)
edgeset.remove(short_edge)
edgeset.remove(short_edge[::-1])
edgeset.add((p, q))
edgeset.add((q, p))
edges = get_edges(coords)
direct_edges(edges)
edgeset = set(edges)
lengthen_dead_ends(edgeset, 1e-10)
normalize_edges(edgeset)
areas = list(get_areas(edgeset))
print("found", len(areas), "areas")
每个areas
中的条目表示一条路径,形式是一个(x, y)的元组列表。