У меня есть облако точек координат в пакете 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: массив должен быть квадратным.
Вот простое решение, которое требует лишь составляющей:
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
Она возвращает логическое значение, массив, где "истинные" значения указывают точки, которые лежат в заданной выпуклой оболочки. Он может быть использован такой:
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 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.
Во-первых, получить выпуклую оболочку для облака точек.
Затем цикл по всем краям выпуклой оболочки в порядке обхода против часовой стрелки. Для каждого из ребер, проверьте, есть ли ваша целевая точка лежит в "левый" из этого края. При этом, обработать кромки, как векторы, указывающие против часовой стрелки вокруг выпуклой оболочки. Если проектное положение точки находится в "левый" из всех векторов, то оно содержится в многоугольник; в противном случае, это лежит вне многоугольника.
Другой переполнение стека теме включает в себя решение найти что-то "стороне" линии точка на: Определить, какой стороне линии точка лежит
<ч /> Сложность выполнения такого подхода (когда у вас уже есть выпуклая оболочка) является О(Н) где n-количество ребер, которые выпуклого корпуса.
Обратите внимание, что это будет работать только для выпуклых многоугольников. Но вы'ре дело с выпуклым днищем, поэтому она должна удовлетворить ваши потребности.
Похоже, у вас уже есть способ, чтобы получить выпуклую оболочку для облака точек. Но если вы обнаружите, что вам придется реализовать свой собственный, в Википедии есть хороший список алгоритмов выпуклой вот КАСКО: Выпуклое Алгоритмы Корпус
Просто для полноты, вот бедный'ы человек решение:
в
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()
Идея проста: вершины выпуклой оболочки набора точек П
выиграл't изменить, если вы добавляете точку "П", который падает на "внутри" в корпусе; вершины выпуклой оболочки для [Р1, Р2, ..., рп]
и [Р1, Р2, ..., РП, Р]
являются одинаковыми. Но если P
падает на "снаружи" и вершины должны меняться.
Это работает для n-размеров, но вы должны вычислить ConvexHull
дважды.
Два примера сюжетов в 2-Д:
Ложные:
Правда:
Использовать атрибут уравнения
о 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)
В словах, смысл в каско, если и только если для каждого уравнения (с описанием граней) скалярное произведение между точкой и нормальным вектором (экв[:-1]
) и смещения (экв[-1]
) меньше или равна нулю. Можно сравнить с небольшим, толерантности положительная константа ` = 1е-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)
Похоже, вы используете 2Д облака точек, поэтому я'd, как, чтобы направить вас на [включение] тест1 для точка-в-полигон для испытания выпуклых многоугольников.
Составляющей'алгоритм выпуклой оболочки S позволяет для нахождения выпуклых оболочек в 2 или более измерений, который является более сложным, чем он должен быть для 2D-облако точек. Поэтому я рекомендую использовать другой алгоритм, такой как этот. Это потому, что вам действительно нужно для точка-в-полигон для испытания выпуклой оболочки-это список выпуклыми точками корпуса по часовой стрелке, а точка внутри многоугольника.
Время этот подход как следовать:
Где n-число точек в облаке и H-число точек в облака точек выпуклой оболочки.
Если вы хотите сохранить с составляющей, вы должны выпуклую оболочку (вы сделали это)
>>> 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 точки.
Если вы хотите знать, если точки из исходного набора данных на КАСКО, то вы сделали.
Я то, что вы хотите знать, если какая-либо точка внутри корпуса или снаружи, вы должны сделать немного больше работы. Что вам придется сделать, может быть
на все кромки соединения двух симплексов из вашего корпуса: решить, будет ли ваша точка находится выше или ниже
Как ускорить, как только точка была выше и ниже друг друга, он находится внутри корпуса.
Основываясь на этом пост, вот мой быстрый и грязный решение для выпуклых областей с 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.')