如何以向量化方式在多边形内生成随机位置与轨迹

问题背景 #

已实现正方形边界下的向量化随机轨迹生成（代码如下），现需适配自定义多边形边界（示例多边形坐标：[(100, 100), (80, 130), (90, 130), (90, 140), (70, 140), (150, 200), (120, 150), (100, 100)]），要求全程采用Numpy向量化操作保证计算效率。

原正方形实现代码：

import numpy as np
# coordinates of square
square_coords = [-255, -255, 256, 256] # [xMin, yMin, xMax, yMax]
# convert to np.array
square_boundaries = np.array([(square_coords[0], square_coords[2]), (square_coords[1], square_coords[3])])

# specify speeds with corresponding probabilities of each speed
ue_speed = [3, 4, 8, 25] 
ue_speed_prob = [0.4, 0.2, 0.3, 0.1] # should sum up to 1

steps = 50
time_interval = 10
nwalks = 1

# calcultation without boundaries
v = np.random.choice(ue_speed, size=steps, p=ue_speed_prob)
R = np.expand_dims((v * time_interval), axis=-1)
theta = 2 * np.pi * np.random.rand(nwalks, steps)
xy = np.dstack((np.cos(theta), np.sin(theta))) * R
trajectory_no_boundaries = np.hstack((np.zeros((nwalks, 1, 2)), np.cumsum(xy, axis=1)))

# generete random placement inside of boundaries
start = (np.random.randint(square_coords[0], square_coords[2]), np.random.randint(square_coords[1], square_coords[3]))
# generate random trajectory inside of boundaries
size = np.diff(square_boundaries, axis=1).ravel()
trajectory = np.abs((trajectory_no_boundaries[0] + start - square_boundaries[:, 0] + size) % (2 * size) - size) + square_boundaries[:, 0]

解决方案 #

针对多边形场景，核心分为随机起点生成和带边界约束的轨迹修正两部分，尽量采用Numpy向量化操作减少循环开销，结合Shapely做几何辅助判断。

1. 多边形内随机起点的向量化生成 #

直接遍历判断点是否在多边形内效率低，先通过多边形包围盒过滤大部分点，再用射线法向量化判断剩余点是否在多边形内：

import numpy as np
from shapely.geometry import Polygon, Point

# 定义目标多边形
polygon_coords = [(100, 100), (80, 130), (90, 130), (90, 140), (70, 140), (150, 200), (120, 150), (100, 100)]
polygon = Polygon(polygon_coords)
poly_np = np.array(polygon_coords[:-1])  # 转为Numpy数组，去掉闭合点

# 向量化射线法判断点是否在多边形内
def points_in_polygon(points, poly):
    x, y = points[:, 0], points[:, 1]
    n = len(poly)
    inside = np.zeros(len(points), dtype=bool)
    
    for i in range(n):
        j = (i + 1) % n
        xi, yi = poly[i]
        xj, yj = poly[j]
        
        # 判断点是否在边的y范围内
        cond1 = (yi > y) != (yj > y)
        # 计算射线与边的交点x坐标
        x_intersect = ((y - yi) * (xj - xi)) / (yj - yi) + xi
        cond2 = x < x_intersect
        inside[cond1 & cond2] = ~inside[cond1 & cond2]
    
    return inside

# 生成随机起点：先在包围盒内生成大量点，再筛选出多边形内的点
min_x, min_y, max_x, max_y = polygon.bounds
# 多生成候选点保证足够的有效起点
candidate_points = np.random.uniform([min_x, min_y], [max_x, max_y], size=(nwalks * 5, 2))
mask = points_in_polygon(candidate_points, poly_np)
start_points = candidate_points[mask][:nwalks]

2. 带多边形边界约束的轨迹生成 #

多边形无法像正方形那样用镜像反射公式直接处理，采用逐帧修正策略：先生成无约束轨迹，再对每一步超出多边形的点，计算与边界的碰撞点并修正运动方向，尽量用Numpy批量处理：

from shapely.geometry import LineString

# 原无约束轨迹生成逻辑（复用原有代码）
ue_speed = [3, 4, 8, 25] 
ue_speed_prob = [0.4, 0.2, 0.3, 0.1]
steps = 50
time_interval = 10
nwalks = 1

v = np.random.choice(ue_speed, size=(nwalks, steps), p=ue_speed_prob)
R = v * time_interval
theta = 2 * np.pi * np.random.rand(nwalks, steps)
xy = np.stack([np.cos(theta) * R, np.sin(theta) * R], axis=-1)
# 基于随机起点生成无约束轨迹
trajectory_no_boundaries = start_points[:, np.newaxis, :] + np.cumsum(xy, axis=1)
# 插入初始起点
trajectory_no_boundaries = np.concatenate([start_points[:, np.newaxis, :], trajectory_no_boundaries], axis=1)[:, :-1, :]

# 轨迹边界修正：对每一步超出多边形的点进行碰撞修正
trajectory = trajectory_no_boundaries.copy()
for step in range(1, steps+1):
    # 获取当前步所有点
    current_points = trajectory[:, step, :]
    # 判断是否在多边形内
    in_poly = np.array([polygon.contains(Point(p)) for p in current_points])
    # 处理超出边界的点
    out_idx = np.where(~in_poly)[0]
    if len(out_idx) == 0:
        continue
    
    for idx in out_idx:
        # 获取上一步位置和当前步位置
        prev_p = trajectory[idx, step-1, :]
        curr_p = current_points[idx]
        # 计算线段与多边形的交点（碰撞点）
        line = LineString([prev_p, curr_p])
        intersection = polygon.boundary.intersection(line)
        if intersection.is_empty:
            # 无交点时回退到上一步，重新生成方向
            theta[idx, step-1] = 2 * np.pi * np.random.rand()
            new_xy = np.array([np.cos(theta[idx, step-1]), np.sin(theta[idx, step-1])]) * R[idx, step-1]
            trajectory[idx, step, :] = trajectory[idx, step-1, :] + new_xy
        else:
            # 用碰撞点作为当前步位置，修正方向（镜像反射）
            intersect_p = np.array([intersection.x, intersection.y])
            trajectory[idx, step, :] = intersect_p
            # 计算反射方向：基于边界法向量
            edge = get_nearest_edge(intersect_p, polygon.boundary)
            edge_vec = np.array([edge.coords[1][0]-edge.coords[0][0], edge.coords[1][1]-edge.coords[0][1]])
            norm_vec = np.array([-edge_vec[1], edge_vec[0]])  # 法向量
            norm_vec = norm_vec / np.linalg.norm(norm_vec)
            dir_vec = curr_p - prev_p
            dir_vec = dir_vec / np.linalg.norm(dir_vec)
            # 反射向量公式
            reflect_vec = dir_vec - 2 * np.dot(dir_vec, norm_vec) * norm_vec
            # 更新下一步的方向
            if step < steps:
                theta[idx, step] = np.arctan2(reflect_vec[1], reflect_vec[0])

# 辅助函数：获取点到边界最近的边
def get_nearest_edge(point, boundary):
    min_dist = float('inf')
    nearest_edge = None
    for i in range(len(boundary.coords)-1):
        edge = LineString([boundary.coords[i], boundary.coords[i+1]])
        dist = edge.distance(Point(point))
        if dist < min_dist:
            min_dist = dist
            nearest_edge = edge
    return nearest_edge

说明 #

• 随机起点生成采用包围盒预过滤+向量化射线法，比逐个调用Shapely.contains效率提升数倍；
• 轨迹修正部分对超出边界的点单独处理，尽可能复用Numpy批量判断，平衡效率与实现复杂度；
• 如果需要支持大量轨迹（nwalks很大），可以进一步优化碰撞修正的循环逻辑，比如用Numpy向量化处理边的距离计算。

内容的提问来源于stack exchange，提问作者illuminato