如何快速判断 A - (A ∩ (ΣBi)) 集合是否存在并从中选取一个随机值?

问题描述 投票:0回答:1

我试图使我的代码随机并快速找到A光盘中的一个点(xy),如果该点确实存在,则该点不属于任何其他现有光盘。我发现的唯一方法是一种漫长而费力的方法,因为它随机选择一个点并检查它是否符合我的标准。然而,如果超过 99% 的 A 圆盘与其他 Bi 圆盘相交,这显然几乎是不可能的。

我似乎找不到更好的解决方案。手工操作很简单,但在算法的离散领域则更难。大家有什么想法吗?

这里是问题的说明,其中蓝色部分代表可以随机选择点的区域: illustration of problem

python math random geometry-surface
1个回答
1
投票

我们可以使用

shapely
来近似圆并找到设定的差异。 为了从结果差异中进行采样,我们可以使用此线程中的答案:https://codereview.stackexchange.com/questions/69833/ 并将其应用到乘差中的每个多边形上,并按每个多边形的面积加权以分布采样点。

这里我们设置一个基本的圆形对象并调整采样方法以确保每个点都在多边形内。

import random
from collections import Counter
import shapely
from shapely.affinity import affine_transform
from shapely.geometry import Point, Polygon, MultiPolygon
from shapely.ops import triangulate


class Circle:
    def __init__(self, x: float, y: float, radius: float):
        self.x = x
        self.y = y
        self.radius = radius
        self._geom = Point(x, y).buffer(radius, quad_segs=90)
    
    def __repr__(self) -> str:
        return f'<Circle(x={self.x}, y={self.y}, radius={self.radius})'
    
    def difference(self, other: 'Circle') -> Polygon|MultiPolygon:
        return shapely.difference(self._geom, other._geom).normalize()
    
    def difference_all(self, *others: 'Circle') -> Polygon|MultiPolygon:
        return shapely.difference(
            self._geom, 
            shapely.union_all([o._geom for o in others])
        ).normalize()
    

# adapted from https://codereview.stackexchange.com/questions/69833/
def random_points_in_polygon(polygon: Polygon, k: int) -> list[Point]:
    "Return list of k points chosen uniformly at random inside the polygon."
    areas = []
    transforms = []
    for t in triangulate(polygon):
        areas.append(t.area)
        (x0, y0), (x1, y1), (x2, y2), _ = t.exterior.coords
        transforms.append([x1 - x0, x2 - x0, y2 - y0, y1 - y0, x0, y0])
    points = []
    for transform in random.choices(transforms, weights=areas, k=k):
        x, y = [random.random() for _ in range(2)]
        if x + y > 1:
            p = Point(1 - x, 1 - y)
        else:
            p = Point(x, y)
        p_transformed = affine_transform(p, transform)
        if polygon.contains(p_transformed):
            points.append(p_transformed)
    if len(points) < k:
        points.extend(random_points_in_polygon(polygon, k-len(points)))
    return points

这里是我们处理结果差异的地方,它可以是多边形或多多边形。 如果是多多边形,

k
样本点将根据面积比率分布到每个子多边形。

def sample_difference(ps: Polygon|MultiPolygon, k: int) -> list[Point]:
    """
    Returns k samples from the difference of multiple Circles.
    """
    if isinstance(ps, Polygon):
        ps = MultiPolygon([ps])
    poly_indices = list(range(len(ps.geoms)))
    sample_probs = [p.area / ps.area for p in ps.geoms]
    poly_index_k_map = Counter(random.choices(poly_indices, weights=sample_probs, k=k))

    points = []
    for ix, p in enumerate(ps.geoms):
        k_poly = poly_index_k_map.get(ix, 0)
        if not k_poly:
            continue
        poly_points = random_points_in_polygon(p, k_poly)
        points.extend(poly_points)
    return points

这是一个例子:

a = Circle(0, 0, 5.1)
b_list = (
    [Circle(1, i, 1.5) for i in range(-5, 6)]
    + [Circle(-3, i, 1.5) for i in range(-5, 6)]
)
# the difference cuts a into 3 segments
a_blist_diff = a.difference_all(*b_list)
sampled_points = sample_difference(a_blist_diff, 100)

这就是

a
b_list
重叠的样子:

enter image description here

这是

a_blist_diff

enter image description here

© www.soinside.com 2019 - 2024. All rights reserved.