У меня есть облако точек координат в 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: массив должен быть квадратным.
Вот простое решение, которое требует только 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')
Привет, я не уверен, как использовать вашу библиотеку программ для достижения этой цели. Но есть простой алгоритм для достижения описанного в словах:
Сначала получим выпуклую оболочку для вашего облака точек.
Затем контур по всем ребрам выпуклой оболочки в порядке против часовой стрелки. Для каждого из ребер проверьте, лежит ли ваша целевая точка на "левом" краю. При этом обработайте края как векторы, указывающие против часовой стрелки вокруг выпуклой оболочки. Если целевой точкой является "левый" всех векторов, то он содержится в многоугольнике; иначе он лежит вне многоугольника.
В этом другом разделе "Переполнение стека" есть решение найти, какая "сторона" строки содержит точку: Определите, на какой стороне линии лежит точка
Обратите внимание, что это будет работать только для выпуклых многоугольников. Но вы имеете дело с выпуклым корпусом, поэтому он должен соответствовать вашим потребностям.
Похоже, у вас уже есть способ получить выпуклый корпус для вашего облака точек. Но если вы обнаружите, что вам нужно реализовать свои собственные, у Википедии есть хороший список выпуклых алгоритмов корпуса здесь: Алгоритмы выпуклых корпусов
Просто для полноты, вот решение для бедных:
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:
False
True
Похоже, вы используете 2D-облако точек, поэтому я хотел бы направить вас на тест включения для point-in- полигональное тестирование выпуклых многоугольников.
Скриптовый алгоритм выпуклого корпуса позволяет находить выпуклые оболочки в 2 или более измерениях, что является более сложным, чем это необходимо для двумерного точечного облака. Поэтому я рекомендую использовать другой алгоритм, например этот. Это потому, что, поскольку вам действительно нужно для точечного полигона тестирования выпуклой оболочки, это список выпуклых точек корпуса по часовой стрелке и точка, находящаяся внутри многоугольника.
Время выполнения этого подхода выполняется следующим образом:
Где N - число точек в облаке точек, а h - количество точек в выпуклой оболочке точечных облаков.
Используйте атрибут 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)
Я бы не использовал алгоритм выпуклого корпуса, потому что вам не нужно вычислять выпуклую оболочку, вы просто хотите проверить, может ли ваша точка выражаться как выпуклая комбинация множества точек, из которых подмножество определяет выпуклую оболочку корпус. Более того, обнаружение выпуклой оболочки является дорогостоящим вычислительным, особенно в более высоких измерениях.
В самом деле, простая задача выяснить, может ли точка быть выражена в виде выпуклой комбинации другого набора точек, можно сформулировать как задачу линейного программирования.
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_[Z.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 измерениях. Время выполнения находится в диапазоне ms. Не хотел бы знать, сколько времени это займет с QHull.
Если вы хотите держаться с scipy, вы должны иметь выпуклый корпус (вы сделали это)
>>> 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 точки.
Если вы хотите узнать, находится ли точка из вашего исходного набора данных на корпусе, то вы закончили.
Я хочу, чтобы узнать, находится ли какая-либо точка внутри корпуса или снаружи, вы должны сделать немного больше работы. То, что вам нужно сделать, может быть
для всех ребер, соединяющих два симплекса вашего корпуса: решите, находится ли ваша точка выше или ниже
если точка находится ниже всех строк или выше всех строк, она находится вне корпуса
Как ускорение, как только точка была выше одной линии и ниже друг друга, она находится внутри корпуса.
Основываясь на этом сообщении, вот мое быстрое и грязное решение для выпуклых областей с четырьмя сторонами (вы можете легко расширить его до более)
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.')
