Python 中的区域交集

我有一个代码，它采用条件 C 作为输入，并将问题的解决方案计算为 (x,y) 空间上的“允许区域”A。该区域由多个“管”组成，这些“管”由两条永远不会交叉的线定义。

我要寻找的最终结果必须满足 k 个条件 {C1, .., Ck}，因此是 k 个区域 {A1, .., Ak} 之间的交集 S。

这是一个有 2 个条件的示例（A1：绿色，3 管。A2：紫色，1 管）；溶液S为红色。

当我处理 4 个区域（每个区域大约有 10 个管子）时，如何找到 S？ （最后的剧情太糟糕了！）

我需要能够绘制它，并找到 S 中点的平均坐标和方差（每个坐标的方差）。 [如果有一种有效的方法来知道点 P 是否属于 S，我将使用蒙特卡罗方法]。

理想情况下，我还希望能够实现从 S 中删除的“禁止管道”[这可能比将 S 与禁止区域的外部相交要复杂一些，因为来自同一区域的两个管道可以交叉（即使定义管子的线永远不会交叉）]。

Note:

• 该代码还存储线的弧长。

• 这些线存储为点数组（每条线大约 1000 个点）。定义管子的两条线不一定具有相同数量的点，但 Python 可以在 1 秒内将所有这些线作为其弧长的函数进行插值。

• 这些线是参数函数（即我们不能写 y = f(x)，因为线允许垂直）。

• 使用绘图编辑绘图以获得右侧的结果...效率不高！

Edit:

• 我不知道如何使用 plt.fill_ Between 进行多重交叉（我可以在此处执行两个条件，但当眼睛判断的行太多时，我需要代码自动执行此操作）。

• 现在我只是生成线条。我没有写任何东西来寻找最终的解决方案，因为我绝对不知道哪种结构最适合这个。 [但是，以前版本的代码能够找到两个不同管子的线之间的交点，我计划将它们作为多边形传递给 shapely，但这意味着其他几个问题..]

• 我不认为我能做到sets：以所需精度扫描整个 (x,y) 区域代表大约 6e8 个点...[由于可变步长（适应曲率），这些线只有 1e3 个点，但整个问题相当大]

问题用Shapely解决了！

我将每个管子定义为Polygon，面积 A 是MultiPolygon作为其管子的联合而构建的对象。

The intersection然后方法计算我正在寻找的解决方案（所有区域之间的重叠）。

整个事情几乎是瞬间发生的。我不知道 Shapely 对于大型物体如此有效 [每管大约 2000 个点，每个区域 10 个管，4 个区域]。

感谢您的帮助！ :)

Edit:

一个工作示例。

import matplotlib.pyplot as plt
import shapely
from shapely.geometry import Polygon
from descartes import PolygonPatch
import numpy as np

def create_tube(a,height):
    x_tube_up = np.linspace(-4,4,300)
    y_tube_up = a*x_tube_up**2 + height
    x_tube_down = np.flipud(x_tube_up)          #flip for correct definition of polygon
    y_tube_down = np.flipud(y_tube_up - 2)

    points_x = list(x_tube_up) + list(x_tube_down)
    points_y = list(y_tube_up) + list(y_tube_down)

    return Polygon([(points_x[i], points_y[i]) for i in range(600)])

def plot_coords(ax, ob):
    x, y = ob.xy
    ax.plot(x, y, '+', color='grey')


area_1 = Polygon()          #First area, a MultiPolygon object
for h in [-5, 0, 5]:
    area_1 = area_1.union(create_tube(2, h))

area_2 = Polygon()
for h in [8, 13, 18]:
    area_2 = area_2.union(create_tube(-1, h))

solution = area_1.intersection(area_2)      #What I was looking for

##########  PLOT  ##########

fig = plt.figure()
ax = fig.add_subplot(111)

for tube in area_1:
    plot_coords(ax, tube.exterior)
    patch = PolygonPatch(tube, facecolor='g', edgecolor='g', alpha=0.25)
    ax.add_patch(patch)

for tube in area_2:
    plot_coords(ax, tube.exterior)
    patch = PolygonPatch(tube, facecolor='m', edgecolor='m', alpha=0.25)
    ax.add_patch(patch)

for sol in solution:
    plot_coords(ax, sol.exterior)
    patch = PolygonPatch(sol, facecolor='r', edgecolor='r')
    ax.add_patch(patch)

plt.show()

和情节：