Как определить, какие точки находятся внутри многоугольника, а какие нет (большое количество точек)?

У меня есть большой набор точек данных (100 000+), хранящихся в 2-мерном массиве numpy (1-й столбец: координаты x, координаты 2-го столбца: y). У меня также есть несколько одномерных массивов, сохраняющих дополнительную информацию для каждой точки данных. Теперь я хотел бы создавать графики из подмножеств этих 1D-массивов, которые включают только точки, находящиеся в заданном полигоне.

Я придумал следующее решение, которое не является ни элегантным, ни быстрым:

#XY is the 2D array.
#A is one of the 1D arrays.
#poly is a matplotlib.patches.Polygon

mask = np.array([bool(poly.get_path().contains_point(i)) for i in XY])

matplotlib.pylab.hist(A[mask], 100)
matplotlib.pylab.show()

Не могли бы вы помочь мне улучшить этот код? Я пытался играть с np.vectorize вместо понимания списка, но не смог заставить его работать.

Ответы

Ответ 1

Используйте matplotlib.nxutils.points_inside_poly, который реализует очень эффективный тест.

Примеры и дальнейшее объяснение этого 40-летнего алгоритма в matplotlib часто задаваемых вопросов.

Обновление: Обратите внимание, что points_inside_poly устарел с версии 1.2.0 из matplotlib. Вместо этого используйте matplotlib.path.Path.contains_points.

Ответ 2

Я боюсь, что не знаком с библиотеками, которые вы используете, но я думаю, что у меня есть разумная идея для алгоритма, который вы могли бы использовать, и я просто займусь тем, как я буду реализовывать это с помощью vanilla python, а затем Я уверен, что вы можете улучшить его и реализовать с помощью этих библиотек. Кроме того, я не утверждаю, что это лучший способ добиться этого, но я хотел получить ответ в разумной степени, так что вот здесь.

Теперь идея исходит из использования кросс-произведения двух векторов в алгоритмах для нахождения выпуклого множества множества точек, например. Graham Scan. Скажем, мы имеем две точки p1 и p2, которые определяют точечные векторы p1 и p2, начиная с начала (0,0) до (x1, y1) и (x2, y2) соответственно. Перекрестное произведение p1 x p2 дает третий вектор p3, который перпендикулярен как p1, так и p2 и имеет величину, заданную площадью параллелограмма, ограниченного векторами.

Очень полезный результат состоит в том, что определитель матрицы

/ x1, x2 \
\ y1, y2 /

..., который является x1 * y2 - x2 * y1, дает величину вектора p3, а знак указывает, является ли p3 "выходящим" из плоскости или "входить" в него. Ключевым моментом здесь является то, что если эта величина положительная, то p2 находится "слева" от p1, а если она отрицательная, то p2 вправо " p1.

Надеюсь, этот пример искусства ascii поможет:

    . p2(4, 5)
   /
  /
 /
/_ _ _ _ _. p1(5, 0)

x1 * y2 - x2 * y1 = 5 * 4 - 0 * 5 = 20, и поэтому p2 находится "слева" от p1

Наконец, почему это полезно для нас! Если у нас есть список вершин многоугольника и множество других точек графа, то для каждого ребра многоугольника мы можем получить вектор этого ребра. Мы также можем получить векторы, соединяющие стартовую вершину со всеми остальными точками графа и путем проверки того, лежат ли они слева или справа от края, мы можем исключить некоторые точки для каждого ребра. Все те, которые не удалены в конце процесса, это те точки внутри многоугольника. Во всяком случае, на какой-то код, чтобы сделать это более понятным!

Получите список вершин вашего многоугольника в том порядке, в котором вы их посещали, если бы вы рисовали их против часовой стрелки, например, какой-то пятиугольник мог бы быть:

poly = [(1, 1), (4, 2), (5, 5), (3, 8), (0, 4)]

Получите набор, содержащий все остальные точки на графике, мы постепенно удалим недопустимые точки из этого набора, пока те, которые остались в конце процесса, не будут точно такими точками, которые находятся внутри многоугольника.

points = set(['(3, 0), (10, -2), (3,3), ...])

Основной бит самого кода на самом деле довольно компактен для того, как долго мне приходилось писать о том, как он работает. to_right принимает два кортежа, представляющих векторы, и возвращает True, если v2 лежит справа от v1. Затем петли проходят через все края многоугольника и удаляют точки из рабочего набора, если они находятся справа от любого из ребер.

def to_right(v1, v2):
    return (v1[0]*v2[1] - v1[1]*v2[0]) < 0

for i in range(len(poly)):
    v1 = poly[i-1]
    v2 = poly[i]
    for p in points:
        if(to_right(v2-v1, p-v1)):
            points.remove(p)

edit: Чтобы прояснить, тот факт, что они удалены, если они справа, а не слева, связаны с порядком, в котором указаны вершины многоугольника. Если бы они были в порядке по часовой стрелке, вам нужно было бы исключить левые точки вместо этого. На данный момент у меня нет особого решения этой проблемы.

В любом случае, надеюсь, что я прав по этому поводу, и это может помочь кому-то, даже если не OP. Асимптотическая сложность этого алгоритма равна O (mn), где n - количество точек в графе, а m - количество вершин многоугольника, так как в худшем случае все точки лежат внутри многоугольника, и мы должны проверять каждую точку для каждого ребра, при этом никто не удаляется.