Python中的区域交集 [英] Area intersection in Python

问题描述

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

I have a code that takes a condition C as an input, and computes the solution to my problem as an 'allowed area' A on the (x,y) space. This area is made of several 'tubes', which are defined by 2 lines that can never cross.

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

The final result I'm looking for must satisfy k conditions {C1, .., Ck}, and is therefore an intersection S between k areas {A1, .. , Ak}.

这是一个有 2 个条件的例子(A1:绿色，3 管.A2:紫色，1 管)；溶液 S 为红色.

Here is an example with 2 conditions (A1: green, 3 tubes. A2: purple, 1 tube); the solution S is in red.

当我处理 4 个区域，每个区域大约有 10 个管时，如何找到 S?(最后的剧情太烂了！)

How can I find S when I'm dealing with 4 areas of around 10 tubes each? (The final plot is awful!)

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

I would need to be able to plot it, and to find the mean coordinate and the variance of the points in S (variance of each coordinate). [If there is an efficient way of knowing whether a point P belongs to S or not, I’ll just use a Monte Carlo method].

理想情况下，我还希望能够实现我将从 S 中移除的禁管" [这可能比将 S 与我的禁区外部相交更复杂一些，因为来自同一个的两个管区域可以交叉(即使定义管的线永远不会交叉)].

Ideally, I’d also like to be able to implement "forbidden tubes" that I would remove from S [it might be a bit more complicated than intersecting S with the outside of my forbidden area, since two tubes from the same area can cross (even if the lines defining a tube never cross)].

注意:

• 代码还存储了线条的弧长.

• The code also stores the arc length of the lines.

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

The lines are stored as arrays of points (around 1000 points per line). The two lines defining a tube do not necessarily have the same number of points, but Python can interpolate ALL of them as a function of their arc length in 1 second.

线是参数函数(即我们不能写 y = f(x)，因为线可以是垂直的).

The lines are parametric functions (i.e. we cannot write y = f(x), since the lines are allowed to be vertical).

用油漆编辑了绘图以获得正确的结果......效率不高！

The plot was edited with paint to get the result on the right... Not very efficient!

• 我不知道如何使用plt.fill_between进行多个交叉点(我可以在这里为2个条件做，但我需要代码在行数过多时自动进行肉眼判断).

• I don't know how I can use plt.fill_between for a multiple intersection (I can do it here for 2 conditions, but I need the code to do it automatically when there are too many lines for eye judgement).

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

For now I just generate the lines. I didn’t write anything for finding the final solution since I absolutely don’t know which structure is the most adapted for this. [However, a previous version of the code was able to find the intersection points between the lines of 2 different tubes, and I was planning to pass them as polygons to shapely, but this implied several other problems..]

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

I don't think I can do it with sets: scanning the whole (x,y) area at required precision represents around 6e8 points... [The lines have only 1e3 points thanks to a variable step size (adapts to the curvature), but the whole problem is quite large]

推荐答案

Shapely 解决了问题！

Problem solved with Shapely!

我将每个管定义为一个 Polygon，区域 A 是一个 MultiPolygon 对象，构建为其管的并集.

I defined each tube as a Polygon, and an area A is a MultiPolygon object built as the union of its tubes.

intersection 方法然后计算我正在寻找的解决方案(所有区域之间的重叠).

The intersection method then computes the solution I was looking for (the overlap between all areas).

整个过程几乎是瞬间完成的.我不知道匀称对大物体这么好[每个管大约 2000 个点，每个区域 10 个管，4 个区域].

The whole thing is almost instantaneous. I didn't know shapely was so good with large objects [around 2000 points per tube, 10 tubes per area, 4 areas].

感谢您的帮助！:)

一个有效的例子.

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

情节:

这篇关于Python中的区域交集的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！