python - ¿Cuál es una forma eficiente de descubrir si un punto se encuentra en el casco convexo de una nube de puntos?

numpy convex-hull

Question

Tengo una nube de puntos de coordenadas en numpy. Para un número elevado de puntos, quiero saber si los puntos se encuentran en el casco convexo de la nube de puntos.

Intenté Pyhull pero no puedo averiguar cómo comprobar si hay un punto en el ConvexHull :

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

plantea LinAlgError: la matriz debe ser cuadrada.

Answer 1

Aquí hay una solución fácil que solo requiere 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

Devuelve una matriz booleana donde los valores True indican puntos que se encuentran en el casco convexo dado. Se puede usar así:

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

print in_hull(tested,cloud)

Si tiene instalado matplotlib, también puede usar la siguiente función que llama a la primera y traza los resultados. Solo para datos 2D, dados por matrices 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

Basado en this post, aquí está mi solución rápida y sucia para regiones convexas con 4 lados (puede extenderlo fácilmente a más)

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.'')

Answer 3

Hola, no estoy seguro de cómo usar su biblioteca de programas para lograr esto. Pero hay un algoritmo simple para lograr esto descrito en palabras:

1. Crea un punto que está definitivamente fuera de tu casco. Llámalo Y
2. produzca un segmento de línea que conecte su punto en cuestión (X) con el nuevo punto Y.
3. Haga un bucle alrededor de todos los segmentos del borde de su casco convexo. verifique si cada uno de ellos se cruza con XY.
4. Si el número de intersecciones que contó es par (incluido 0), X está fuera del casco. De lo contrario X está dentro del casco.
5. Si ocurre así, pase XY a través de uno de sus vértices en el casco, o si se superpone directamente con uno de los bordes de su casco, mueva Y un poco.
6. Lo anterior funciona también para el casco cóncavo. Puede ver en la siguiente ilustración (el punto verde es el punto X que está tratando de determinar. El amarillo marca los puntos de intersección.

Answer 4

No utilizaría un algoritmo de casco convexo, ya que no necesita calcular el casco convexo, solo desea verificar si su punto se puede expresar como una combinación convexa del conjunto de puntos de los cuales un subconjunto define un casco convexo. Además, encontrar el casco convexo es computacionalmente costoso, especialmente en dimensiones más altas.

De hecho, el mero problema de descubrir si un punto se puede expresar como una combinación convexa de otro conjunto de puntos se puede formular como un problema de programación lineal.

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))

Para el ejemplo, resolví el problema por 10000 puntos en 10 dimensiones. El tiempo de ejecución está en el rango de ms. No querría saber cuánto tiempo llevaría esto con QHull.

Answer 5

Parece que estás usando una nube de puntos 2D, por lo que me gustaría dirigirte a la prueba de inclusión para las pruebas de punto en polígono de polígonos convexos.

El algoritmo de casco convexo de Scipy permite encontrar cascos convexos en 2 o más dimensiones, lo que es más complicado de lo que debe ser para una nube de puntos 2D. Por lo tanto, recomiendo usar un algoritmo diferente, como este . Esto se debe a que, como realmente necesita una prueba de punto en polígono de un casco convexo, es la lista de puntos de casco convexo en el sentido de las agujas del reloj y un punto que está dentro del polígono.

El tiempo de rendimiento de este enfoque es el siguiente:

• O (N log N) para construir el casco convexo
• O (h) en el preprocesamiento para calcular (y almacenar) los ángulos de cuña desde el punto interior
• O (registro h) por consulta de punto en polígono.

Donde N es el número de puntos en la nube de puntos y h es el número de puntos en el casco convexo de las nubes de puntos.

Answer 6

Primero, obtenga el casco convexo para su nube de puntos.

Luego pase por todos los bordes del casco convexo en el sentido contrario a las agujas del reloj. Para cada uno de los bordes, compruebe si su punto de destino se encuentra a la "izquierda" de ese borde. Al hacer esto, trate los bordes como vectores que apuntan en sentido contrario a las agujas del reloj alrededor del casco convexo. Si el punto objetivo está a la "izquierda" de todos los vectores, entonces está contenido por el polígono; De lo contrario, queda fuera del polígono.

Este otro tema de desbordamiento de pila incluye una solución para encontrar en qué "lado" de una línea está el punto: determinar qué lado de una línea miente un punto

La complejidad del tiempo de ejecución de este enfoque (una vez que ya tiene el casco convexo) es O (n), donde n es el número de bordes que tiene el casco convexo.

Tenga en cuenta que esto funcionará solo para polígonos convexos. Pero se trata de un casco convexo, por lo que debe adaptarse a sus necesidades.

Parece que ya tienes una forma de obtener el casco convexo para tu nube de puntos. Pero si encuentra que tiene que implementar el suyo propio, Wikipedia tiene una buena lista de algoritmos de casco convexo aquí: Algoritmos de casco convexo

Answer 7

Sólo para completar, aquí está la solución de un hombre pobre:

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()

La idea es simple: los vértices del casco convexo de un conjunto de puntos P no cambiarán si agrega un punto p que cae "dentro" del casco; Los vértices del casco convexo para [P1, P2, ..., Pn] y [P1, P2, ..., Pn, p] son los mismos. Pero si p cae "afuera", entonces los vértices deben cambiar. Esto funciona para n-dimensiones, pero debe calcular el ConvexHull dos veces.

Dos diagramas de ejemplo en 2-D:

Falso:

Cierto:

Answer 8

Si quieres mantenerte con scipy, tienes que hacer un casco convexo (lo hiciste)

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

luego construir la lista de puntos en el casco. Aquí está el código del doc para trazar el casco.

>>> 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-'')

Así que a partir de eso, propondría para calcular la lista de puntos en el casco

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

(aunque no lo intenté)

Y también puede venir con su propio código para calcular el casco y devolver los puntos x, y.

Si desea saber si un punto de su conjunto de datos original está en el casco , entonces está listo.

Lo que quieres es saber si un punto está dentro o fuera del casco , debes hacer un poco más de trabajo. Lo que deberías hacer podría ser

• para todos los bordes que unen dos sencillos de su casco: decida si su punto está arriba o abajo

• si el punto está debajo de todas las líneas, o sobre todas las líneas, está fuera del casco

A medida que aumenta la velocidad, tan pronto como un punto ha estado por encima de una línea y por debajo de otro, está dentro del casco.

Answer 9

Utilice el atributo de equations de 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)

En palabras, un punto está en el casco si y solo si para cada ecuación (que describe las facetas) el producto puntual entre el punto y el vector normal ( eq[:-1] ) más el desplazamiento ( eq[-1] ) es menor o igual que cero. Es posible que desee comparar con una tolerance = 1e-12 constante pequeña y positiva tolerance = 1e-12 lugar de cero debido a problemas de precisión numérica (de lo contrario, es posible que un vértice del casco convexo no esté en el casco convexo).

Demostración:

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)