Python中的区域交集

Question

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

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

以下是一个示例，其中有2个条件（A1：绿色，3个管道。A2：紫色，1个管道）；解决方案S为红色。

我如何在处理大约10个管道的4个区域时找到S？（最终图形很糟糕！）

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

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

---

注意：

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

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

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

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

---

编辑：

• 我不知道如何使用plt.fill_between进行多重交点（我可以在这里处理2个条件，但当有太多线条需要眼睛判断时，我需要代码自动执行它）。

• 目前我只是生成了这些线。我没有编写任何内容来找到最终解决方案，因为我绝对不知道哪种结构最适合此问题。[然而，代码的早期版本能够找到两个不同管道的线之间的交点，并且我计划将它们作为多边形传递给shapely，但这涉及到其他几个问题..]

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

- Lixons

---

根据您的最后一条留言，您可以使用plt.fill_between（或plt.fill）来完成此操作。 - Phlya

你能否在问题中添加一些代码，以展示你在实现过程中的具体进度，这样那些愿意帮助你的人就能更清楚地了解你可能卡在哪里，以及他们最擅长的领域是什么？谢谢。 - Dilettant

只是一个初步的想法，如果你能将每个管道计算为一组点，那么你可以对每对管道进行一次“集合”交集运算，以获得有助于你预期结果的信息。 - machine yearning

请查看http://matplotlib.org/examples/pylab_examples/fill_between_demo.html上的示例。 - Serenity

请查看 Shapely 库，它可以帮助解决这些问题：https://pypi.python.org/pypi/Shapely - tfv

如果您可以将管道的方程参数化（例如，管道1-> y = x + x ^ 2-5），则可以使用标准优化工具解决问题。如果您有由点表示的线条，则可以通过多项式近似线条，然后使用优化来找到面积。虽然不是非常容易的问题，但是可以完成。 - Imanol Luengo

Answer 1

使用Shapely解决了问题！

我将每个管道定义为一个Polygon，并将面积A作为MultiPolygon对象构建，其由其管道的联合组成。

然后intersection方法计算出了我要找的解决方案（所有区域之间的重叠部分）。

整个过程几乎是瞬间完成的。我不知道Shapely在处理大型对象时如此出色[每个管道约2000个点，每个区域10个管道，共4个区域]。

谢谢你的帮助！ :)

编辑：

一个可行的示例。

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()

而剧情是：

- Lixons

---

你可以发布一个解决方案的可用示例吗？顺便说一下，你可以将你的答案选为最佳解决方案。 - nicoguaro