Какой эффективный способ определить, лежит ли точка в выпуклой оболочке облака точек?

Question

У меня есть облако точек координат в NumPy. Для большого числа точек я хочу выяснить, лежат ли точки в выпуклой оболочке облака точек.

Я попытался Pyhull, но я не могу понять, как проверить, находится ли точка в ConvexHull:

hull = ConvexHull(np.array([(1, 2), (3, 4), (3, 6)]))
for s in hull.simplices:
    s.in_simplex(np.array([2, 3]))

повышает LinAlgError: массив должен быть квадратным.

Answer 1

Вот простое решение, которое требует только scipy:

def in_hull(p, hull):
    """
    Test if points in `p` are in `hull`

    `p` should be a `NxK` coordinates of `N` points in `K` dimensions
    `hull` is either a scipy.spatial.Delaunay object or the `MxK` array of the 
    coordinates of `M` points in `K`dimensions for which Delaunay triangulation
    will be computed
    """
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    return hull.find_simplex(p)>=0

Возвращает логический массив, где True значения указывают точки, которые лежат в данной выпуклой оболочке. Это можно использовать так:

tested = np.random.rand(20,3)
cloud  = np.random.rand(50,3)

print in_hull(tested,cloud)

Если у вас установлен matplotlib, вы также можете использовать следующую функцию, которая вызывает первую и отображает результаты. Только для 2D данных, заданных Nx2 массивы:

def plot_in_hull(p, hull):
    """
    plot relative to `in_hull` for 2d data
    """
    import matplotlib.pyplot as plt
    from matplotlib.collections import PolyCollection, LineCollection

    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)

    # plot triangulation
    poly = PolyCollection(hull.points[hull.vertices], facecolors='w', edgecolors='b')
    plt.clf()
    plt.title('in hull')
    plt.gca().add_collection(poly)
    plt.plot(hull.points[:,0], hull.points[:,1], 'o', hold=1)


    # plot the convex hull
    edges = set()
    edge_points = []

    def add_edge(i, j):
        """Add a line between the i-th and j-th points, if not in the list already"""
        if (i, j) in edges or (j, i) in edges:
            # already added
            return
        edges.add( (i, j) )
        edge_points.append(hull.points[ [i, j] ])

    for ia, ib in hull.convex_hull:
        add_edge(ia, ib)

    lines = LineCollection(edge_points, color='g')
    plt.gca().add_collection(lines)
    plt.show()    

    # plot tested points `p` - black are inside hull, red outside
    inside = in_hull(p,hull)
    plt.plot(p[ inside,0],p[ inside,1],'.k')
    plt.plot(p[-inside,0],p[-inside,1],'.r')

Answer 2

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

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

import numpy as np
from scipy.optimize import linprog

def in_hull(points, x):
    n_points = len(points)
    n_dim = len(x)
    c = np.zeros(n_points)
    A = np.r_[points.T,np.ones((1,n_points))]
    b = np.r_[x, np.ones(1)]
    lp = linprog(c, A_eq=A, b_eq=b)
    return lp.success

n_points = 10000
n_dim = 10
Z = np.random.rand(n_points,n_dim)
x = np.random.rand(n_dim)
print(in_hull(Z, x))

Для примера я решил задачу на 10000 точек в 10 измерениях. Время выполнения находится в диапазоне мс. Не хотел бы знать, сколько времени это займет с QHull.

Answer 3

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

1. создать точку, которая определенно находится за пределами вашего корпуса. Назовите это Y
2. создайте отрезок, соединяющий вашу точку (X) с новой точкой Y.
3. Обведите все сегменты вашего выпуклого корпуса. проверьте для каждого из них, пересекается ли сегмент с XY.
4. Если число пересечений, которое вы посчитали, является четным (включая 0), то X находится вне корпуса. В противном случае X находится внутри корпуса.
5. если это произойдет, XY пройдите через одну из ваших вершин на корпусе или непосредственно пересекаются с одной из кромок вашего корпуса, немного сдвиньте Y.
6. вышесказанное работает и для вогнутого корпуса. Вы можете видеть на рисунке ниже (Зеленая точка - это точка X, которую вы пытаетесь определить. Желтая точка обозначает точки пересечения.

Answer 4

Во-первых, получите выпуклую оболочку для вашего облака точек.

Затем зациклите все ребра выпуклой оболочки против часовой стрелки. Для каждого из ребер проверьте, находится ли ваша целевая точка слева от этого ребра. При этом рассматривайте края как векторы, направленные против часовой стрелки вокруг выпуклой оболочки. Если целевая точка находится "слева" от всех векторов, то она содержится в многоугольнике; в противном случае он лежит вне многоугольника.

Этот другой раздел переполнения стека включает решение для определения, на какой "стороне" линии находится точка: Определите, на какой стороне линии лежит точка

---

Сложная среда выполнения этого подхода (если у вас уже есть выпуклая оболочка) - это O(n), где n - количество ребер, которые имеет выпуклая оболочка.

Обратите внимание, что это будет работать только для выпуклых многоугольников. Но вы имеете дело с выпуклой оболочкой, поэтому она должна соответствовать вашим потребностям.

Похоже, у вас уже есть способ получить выпуклую оболочку для вашего облака точек. Но если вы обнаружите, что вам нужно реализовать свой собственный, в Википедии есть хороший список алгоритмов выпуклой оболочки: Алгоритмы выпуклой оболочки

Answer 5

Использовать equations атрибут ConvexHull:

def point_in_hull(point, hull, tolerance=1e-12):
    return all(
        (np.dot(eq[:-1], point) + eq[-1] <= tolerance)
        for eq in hull.equations)

Словом, точка находится в корпусе тогда и только тогда, когда для каждого уравнения (описывающего грани) произведение точек между точкой и вектором нормали (eq[:-1]) плюс смещение (eq[-1]) меньше или равно нулю. Вы можете сравнить с небольшой положительной константой tolerance = 1e-12 а не к нулю из-за проблем числовой точности (в противном случае вы можете обнаружить, что вершина выпуклой оболочки не находится в выпуклой оболочке).

Демонстрация:

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

for simplex in hull.simplices:
    plt.plot(points[simplex, 0], points[simplex, 1])

plt.scatter(*points.T, alpha=.5, color='k', s=200, marker='v')

for p in random_points:
    point_is_in_hull = point_in_hull(p, hull)
    marker = 'x' if point_is_in_hull else 'd'
    color = 'g' if point_is_in_hull else 'm'
    plt.scatter(p[0], p[1], marker=marker, color=color)

Answer 6

Просто для полноты, вот решение бедного человека:

import pylab
import numpy
from scipy.spatial import ConvexHull

def is_p_inside_points_hull(points, p):
    global hull, new_points # Remove this line! Just for plotting!
    hull = ConvexHull(points)
    new_points = numpy.append(points, p, axis=0)
    new_hull = ConvexHull(new_points)
    if list(hull.vertices) == list(new_hull.vertices):
        return True
    else:
        return False

# Test:
points = numpy.random.rand(10, 2)   # 30 random points in 2-D
# Note: the number of points must be greater than the dimention.
p = numpy.random.rand(1, 2) # 1 random point in 2-D
print is_p_inside_points_hull(points, p)

# Plot:
pylab.plot(points[:,0], points[:,1], 'o')
for simplex in hull.simplices:
    pylab.plot(points[simplex,0], points[simplex,1], 'k-')
pylab.plot(p[:,0], p[:,1], '^r')
pylab.show()

Идея проста: вершины выпуклой оболочки множества точек P не изменится, если вы добавите точку p который падает "внутри" корпуса; вершины выпуклой оболочки для [P1, P2, ..., Pn] а также [P1, P2, ..., Pn, p] подобные. Но если p падает "наружу", тогда вершины должны измениться. Это работает для n-измерений, но вы должны вычислить ConvexHull дважды.

Два примера графиков в 2-D:

Ложь:

Правда:

Answer 7

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

      import numpy as np
from scipy.spatial.qhull import _Qhull

def in_hull(points, queries):
    hull = _Qhull(b"i", points,
                  options=b"",
                  furthest_site=False,
                  incremental=False, 
                  interior_point=None)
    equations = hull.get_simplex_facet_array()[2].T
    return np.all(queries @ equations[:-1] < - equations[-1], axis=1)

# ============== Demonstration ================

points = np.random.rand(8, 2)
queries = np.random.rand(3, 2)
print(in_hull(points, queries))

Обратите внимание, что я использую нижний уровень _Qhull класс по эффективности.

Answer 8

Похоже, вы используете 2D-облако точек, поэтому я бы хотел направить вас к тесту включения для тестирования выпуклых многоугольников точка-в-многоугольнике.

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

Время выполнения этого подхода выглядит следующим образом:

• O (N log N) для построения выпуклой оболочки
• O(h) при предварительной обработке для вычисления (и сохранения) углов клина от внутренней точки
• O(log h) на запрос точка-полигон.

Где N - количество точек в облаке точек, а h - количество точек в облаках точек выпуклой оболочки.

Answer 9

Для тех, кто заинтересован, я сделал векторизованную версию ответа @charlie- brummit :

      def points_in_hull(p, hull, tol=1e-12):
    return np.all(hull.equations[:,:-1] @ p.T + np.repeat(hull.equations[:,-1][None,:], len(p), axis=0).T <= tol, 0)

куда pтеперь [N,2]множество. Это примерно в 6 раз быстрее, чем рекомендованное решение (@Sildoreth answer), и примерно в 10 раз быстрее, чем оригинал.

Answer 10

Чтобы откатить этот ответ , чтобы сразу проверить все точки в массиве numpy, это сработало для меня:

      import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([(1, 2), (3, 4), (3, 6), (2, 4.5), (2.5, 5)])
hull = ConvexHull(points)

np.random.seed(1)
random_points = np.random.uniform(0, 6, (100, 2))

# get array of boolean values indicating in hull if True
in_hull = np.all(np.add(np.dot(random_points, hull.equations[:,:-1].T),
                        hull.equations[:,-1]) <= tolerance, axis=1)

random_points_in_hull = random_points[in_hull]

Answer 11

Если вы хотите остаться со Сципионом, вы должны выпуклый корпус (вы так и сделали)

>>> from scipy.spatial import ConvexHull
>>> points = np.random.rand(30, 2)   # 30 random points in 2-D
>>> hull = ConvexHull(points)

затем построить список точек на корпусе. Вот код из документа для построения корпуса

>>> import matplotlib.pyplot as plt
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> for simplex in hull.simplices:
>>>     plt.plot(points[simplex,0], points[simplex,1], 'k-')

Итак, начиная с этого, я бы предложил для расчета списка точек на корпусе

pts_hull = [(points[simplex,0], points[simplex,1]) 
                            for simplex in hull.simplices] 

(хотя я не пробовал)

И вы также можете прийти со своим собственным кодом для вычисления оболочки, возвращая точки x,y.

Если вы хотите знать, находится ли точка из вашего исходного набора данных на корпусе, то все готово.

Я хочу знать, находится ли какая-либо точка внутри корпуса или снаружи, вы должны сделать немного больше работы. То, что вы должны будете сделать, может быть

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

• если точка находится ниже всех линий или выше всех линий, она находится за пределами корпуса

В качестве ускорения, как только точка оказывается выше одной линии и ниже друг друга, она оказывается внутри корпуса.

Answer 12

Основываясь на этом посте, вот мое быстрое и грязное решение для выпуклых областей с 4 сторонами (вы можете легко расширить его до большего)

def same_sign(arr): return np.all(arr > 0) if arr[0] > 0 else np.all(arr < 0)

def inside_quad(pts, pt):
    a =  pts - pt
    d = np.zeros((4,2))
    d[0,:] = pts[1,:]-pts[0,:]
    d[1,:] = pts[2,:]-pts[1,:]
    d[2,:] = pts[3,:]-pts[2,:]
    d[3,:] = pts[0,:]-pts[3,:]
    res = np.cross(a,d)
    return same_sign(res), res

points = np.array([(1, 2), (3, 4), (3, 6), (2.5, 5)])

np.random.seed(1)
random_points = np.random.uniform(0, 6, (1000, 2))

print wlk1.inside_quad(points, random_points[0])
res = np.array([inside_quad(points, p)[0] for p in random_points])
print res[:4]
plt.plot(random_points[:,0], random_points[:,1], 'b.')
plt.plot(random_points[res][:,0], random_points[res][:,1], 'r.')