Sympy中迭代相交的线段......有更好的方法吗？

Posted 2023-02-23

【中文标题】Sympy中迭代相交的线段......有更好的方法吗？

【英文标题】：iteratively intersecting line segments in Sympy... is there a better way?

【发布时间】：2015-06-03 16:37:02

【问题描述】：

好的。我有构成多边形边界的点。我想（a）使用 Sympy 的几何模块从任何一对点之间的所有可能的线段中确定哪些线段不跨越周边。该结果将是允许在 (b) Networkx 中的最短距离分析中使用的“边”。我的最终目标是通过许多形状来迭代这个过程，但在这个例子中我只硬编码了 1 个形状的坐标。

import numpy
import networkx as nx
from sympy import geometry
from itertools import combinations
from matplotlib import pyplot as plot
arr_bou = numpy.array([[-542.62545014,  961.34455209],
   [-544.45425379,  961.34455209],
   [-544.45425379,  962.25895392],
   [-547.19745928,  962.25895392],
   [-547.19745928,  963.17335575],
   [-549.02626294,  963.17335575],
   [-549.02626294,  964.08775758],
   [-550.85506659,  964.08775758],
   [-550.85506659,  961.34455209],
   [-552.68387025,  961.34455209],
   [-552.68387025,  962.25895392],
   [-553.59827208,  962.25895392],
   [-553.59827208,  965.91656123],
   [-552.68387025,  965.91656123],
   [-552.68387025,  967.7453649 ],
   [-551.76946842,  967.7453649 ],
   [-551.76946842,  968.65976672],
   [-550.85506659,  968.65976672],
   [-550.85506659,  967.7453649 ],
   [-548.11186111,  967.7453649 ],
   [-548.11186111,  965.91656123],
   [-547.19745928,  965.91656123],
   [-547.19745928,  964.08775758],
   [-546.28305745,  964.08775758],
   [-546.28305745,  965.00215941],
   [-543.53985197,  965.00215941],
   [-543.53985197,  963.17335575],
   [-542.62545014,  963.17335575],
   [-542.62545014,  964.08775758],
   [-540.79664648,  964.08775758],
   [-540.79664648,  963.17335575],
   [-539.88224465,  963.17335575],
   [-539.88224465,  962.25895392],
   [-542.62545014,  962.25895392],
   [-542.62545014,  961.34455209]])

boundXY = []
for i in arr_bou:
    boundXY.append((i[0],i[1]))

points     = [geometry.Point(i) for i in boundXY]
poly       = geometry.Polygon(*points) # use the * first to unpack the points (necessary to avoid errors)


G          = nx.Graph()
positions  =                                   # build a dictionary
for i in xrange(len(boundXY)):                   # that contains coordinates
    positions[i] = boundXY[i]                    # of each node on the graph's perimeter
G.add_path(positions.keys())# add nodes to graph w/ boundary edges
G.add_path([min(G.nodes()),max(G.nodes())])    combos_o   = list(combinations(positions.keys(),2))
combos     = [i for i in combos_o if i not in G.edges()]

keepcombos = []
for combo in combos:
    pt1 = positions[combo[0]]
    pt2 = positions[combo[1]]
    line = geometry.Polygon(pt1,pt2)
    # there are 4 polygon sides that do not count as intersections
    # because 2 sides will intersect a point on each end
    test = True
    for side in poly.sides:
        if side.p1 != geometry.Point(pt1) and side.p1 != geometry.Point(pt2):
            if side.p2 != geometry.Point(pt1) and side.p2 != geometry.Point(pt2):
                if geometry.intersection(line,side):
                    test = False
                    break
                else:
                    try:
                        if poly.encloses(line.midpoint):
                            pass
                        else:
                            test = False
                            break
                    except NotImplementedError:
                        pass
    if test == True:
        keepcombos.append(combo)
G.add_edges_from(keepcombos)

我已经让它适用于小多边形（14 个顶点），但这需要永远，即使只有 35 个顶点，其他多边形仍然会比这更大。

是否有更有效的方法来查找所有多边形内节点对？

谢谢！！

【参考方案1】：

如果凸包是定义包含所有其他多边形的凸多边形的点集

>>> coords = [[-542.62545014,  961.34455209],
...    [-544.45425379,  961.34455209],
...    [-544.45425379,  962.25895392],
...    [-547.19745928,  962.25895392],

...    [-542.62545014,  962.25895392],
...    [-542.62545014,  961.34455209]]
>>> from sympy import *
>>> pts = [Point(*i) for i in coords]
>>> h = convex_hull(*pts)

那么您是否只对不在外围的所有点感兴趣？一旦知道了这些点，您是否只想生成这些点的所有成对组合？

inside = [p for p in pts if p not in h.vertices]
from sympy.utilities.iterables import combinations
pairs = combinations(inside, 2)

看着https://www.desmos.com/calculator/up4k2qkxy9，我发现有几个点似乎落在了边缘，所以也许它们不应该被包含在“内部”点中。但希望你能明白我回答的要点。

【讨论】：

感谢输入，但凸包内的点不是我感兴趣的。“坐标”中的点是多边形的周长，所以里面没有点。我感兴趣的是找到连接完全在多边形内的 2 个点的所有线段。我想我找到了一个更有效的解决方案，我现在将发布我的新代码作为答案。再次感谢您！

【参考方案2】：

我找到了一个解决方案，可以将处理速度提高大约 13 倍（对于具有 35 个点的多边形（如上面列出的数据），问题中代码中的旧方法需要大约 4 小时才能找到多边形内的所有线段. 这个新方法花了 18 分钟。）

上面我遍历了这些点，并且在每次迭代时分别查看每个边界（“边”）以查看线是否相交。我将其更改为将线与整个多边形相交。如果它在内部或边缘上，应该只有2个相交点，所以如果相交的长度>2，我把组合扔掉

for combo in combos:
    pt1  = geometry.Point(positions[combo[0]],evaluate=False)
    pt2  = geometry.Point(positions[combo[1]],evaluate=False)
    line = geometry.Polygon(pt1,pt2)
    try:
        if poly.encloses(line.midpoint):
            pass
        else:
            continue
    except NotImplementedError:
        continue
    intersect = geometry.intersection(line,poly)
    if len(intersect)>2:
        continue
    keepcombos.append(combo)

此列表“keepcombos”现在包含我希望包含在 Dijkstra 路径分析中的所有行（或 networkx“边”）