Qual è un modo efficace per scoprire se un punto si trova nello scafo convesso di una nuvola di punti?

Question

Ho una nuvola di punti di coordinate in numpy. Per un numero elevato di punti, voglio scoprire se i punti si trovano nello scafo convesso della nuvola di punti.

ho cercato pyhull ma riesco a capire come verificare se un punto è nella ConvexHull:

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

solleva solo LinAlgError: Array deve essere quadrata.

Answer 1

Se si desidera mantenere con SciPy, bisogna convesso (avete fatto così)

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

poi costruire l'elenco dei punti sullo scafo. Ecco il codice da doc per tracciare scafo

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

Quindi, a partire da quella, io proporrei per la lista dei punti di calcolo sullo scafo

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

(anche se non ho provato)

E puoi anche venire con il tuo codice per calcolare lo scafo, restituendo i punti x, y.

Se vuoi sapere se un punto dal set di dati originale è sullo scafo, poi si è fatto.

Io quello che vuoi è sapere se un qualsiasi punto è all'interno dello scafo o all'esterno, devi fare un po 'di lavoro in più. Quello che dovrete fare potrebbe essere

• per tutti i bordi di giunzione due simplessi dello scafo: decidere se il vostro punto è al di sopra o sotto

• se il punto è inferiore a tutte le linee, o al di sopra tutte le linee, è esternamente allo scafo

come accelerare, non appena un punto è stato sopra una linea e una sotto l'altra, è all'interno dello scafo.

Answer 2

Sembra che si sta utilizzando una nuvola di punti in 2D, quindi mi piacerebbe di indirizzarvi verso la inclusion test per i test point-in-poligono di poligoni convessi.

L'algoritmo di scafo convesso di Scipy consente di trovare scafi convessi in 2 o più dimensioni che è più complicato di quanto non debba essere per una nuvola di punti 2D. Pertanto, ti consiglio di utilizzare un algoritmo diverso, ad esempio this one. Questo perché, come realmente necessario per il test point-in-poligono di uno scafo convesso, c'è l'elenco dei punti scafo convessi in senso orario e un punto che si trova all'interno del poligono.

Le prestazioni tempo di questo approccio è come segue:

• O (N log N) per costruire convesso
• O (h) in preelaborazione per calcolare (e memorizzare) gli angoli di cuneo da il punto interno
• O (log h) per query point-in-poligono.

Dove N è il numero di punti nella nuvola di punti eh è il numero di punti nelle nuvole di punti convesso dello scafo.

Answer 3

Ecco una soluzione semplice che richiede solo 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 

restituisce una matrice booleana dove True valori indicano punti che giacciono nel dato convesso. Può essere utilizzato in questo modo:

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

print in_hull(tested,cloud) 

Se avete matplotlib installato, è possibile anche utilizzare la seguente funzione che chiama il primo e traccia i risultati. Per i dati 2D solo, fornite dal Nx2 array:

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 4

primo luogo, ottenere il convesso per la nuvola di punti.

Quindi avvolgere su tutti i bordi dello scafo convesso in senso antiorario. Per ciascuno dei bordi, controlla se il tuo obiettivo si trova sulla "sinistra" di quel bordo. Quando lo fai, tratti i bordi come vettori che puntano in senso antiorario attorno allo scafo convesso. Se il punto target è sulla "sinistra" di tutti i vettori, allora è contenuto dal poligono; altrimenti, si trova al di fuori del poligono.

Questo altro argomento Stack Overflow include una soluzione per trovare quale "collaterali" di una linea di un punto è: Determine Which Side of a Line a Point Lies

---

La complessità di esecuzione di questo approccio (una volta che si dispone già di lo scafo convesso) è O (n) dove n è il numero di bordi che ha lo scafo convesso.

Si noti che questo funzionerà solo per i poligoni convessi. Ma hai a che fare con uno scafo convesso, quindi dovrebbe essere adatto alle tue esigenze.

Sembra che tu abbia già un modo per ottenere lo scafo convesso per la tua nuvola di punti. Ma se trovi che devi implementare il tuo, Wikipedia ha una bella lista di algoritmi di convessità dello scafo qui: Convex Hull Algorithms

Answer 5

Ciao Non sono sicuro di come usare la libreria del tuo programma per raggiungere questo obiettivo. Ma c'è un semplice algoritmo per ottenere ciò descritto in parole:

1. creare un punto che è decisamente fuori dallo scafo. Chiamalo Y
2. produce un segmento di linea che collega il tuo punto in questione (X) al nuovo punto Y.
3. loop su tutti i segmenti del tuo scafo convesso. controlla per ognuna di esse se il segmento si interseca con XY.
4. Se il numero di intersezione contato è pari (compreso 0), X è al di fuori dello scafo. Altrimenti X è all'interno dello scafo.
5. in tal caso XY passa attraverso uno dei tuoi vertici sullo scafo, o si sovrappone direttamente con uno dei bordi dello scafo, sposta Y un po '.
6. i lavori di cui sopra per lo scafo concavo pure. Puoi vedere sotto l'illustrazione (il punto verde è il punto X che stai cercando di determinare.Il giallo segna i punti di intersezione.

Answer 6

Solo per completezza, ecco la soluzione uomo di un povero:

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

L'idea è semplice: i vertici del guscio convesso di un insieme di punti P non cambia se si aggiunge un punto p che cade "all'interno" dello scafo; i vertici dello scafo convesso per [P1, P2, ..., Pn] e [P1, P2, ..., Pn, p] sono gli stessi. Ma se p cade "fuori", allora i vertici devono cambiare. Questo funziona per n-dimensioni, ma è necessario calcolare il ConvexHull due volte.

Due esempi di trame in 2-D:

Falso:

Vero:

Answer 7

utilizzare l'attributo equations di 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) 

In parole, un punto è nello scafo se e solo se per ogni equazione (che descrive le sfaccettature) il prodotto scalare tra il punto e il vettore normale (eq[:-1]) più l'offset (eq[-1]) è minore o uguale a zero. Potresti voler confrontare una piccola costante positiva tolerance = 1e-12 piuttosto che zero a causa di problemi di precisione numerica (altrimenti potresti scoprire che un vertice dello scafo convesso non è nello scafo convesso).

Dimostrazione:

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 8

Sulla base di this posta, ecco la mia soluzione veloce-e-sporco per le regioni convesse con 4 lati (si può facilmente estendere a più)

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 9

Non userei un algoritmo di scafo convesso, perché non hai bisogno di calcolare lo scafo convesso, vuoi solo verificare se il tuo punto può essere espresso come una combinazione convessa dell'insieme di punti di cui un sottoinsieme definisce uno scafo convesso. Inoltre, trovare lo scafo convesso è costoso dal punto di vista computazionale, specialmente nelle dimensioni più elevate.

In effetti, il semplice problema di scoprire se un punto può essere espresso come una combinazione convessa di un altro insieme di punti può essere formulato come un problema di programmazione lineare.

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

Per l'esempio, ho risolto il problema per 10000 punti in 10 dimensioni. Il tempo di esecuzione è nel range ms. Non vorrei sapere quanto tempo ci vorrà con QHull.