¿Cómo determinar qué puntos están dentro de un polígono y cuáles no (gran cantidad de puntos)?

Question

Tengo un gran conjunto de puntos de datos (100.000+) almacenados en una matriz numpy bidimensional (1ra columna: coordenadas x, 2da columna: coordenadas y). También tengo varias matrices de 1 dimensión que almacenan información adicional para cada punto de datos. Ahora me gustaría crear gráficos de subconjuntos de estas matrices 1D que incluyen solo los puntos que están en un polígono dado.

me ocurrió con la siguiente solución que no es ni elegante ni rápido:

#XY is the 2D array. 
#A is one of the 1D arrays. 
#poly is a matplotlib.patches.Polygon 

mask = np.array([bool(poly.get_path().contains_point(i)) for i in XY]) 

matplotlib.pylab.hist(A[mask], 100) 
matplotlib.pylab.show() 

¿Me podría ayudar a mejorar este código? Intenté jugar con np.vectorize en lugar de la lista de comprensión, pero no pude lograr que funcionara.

Answer 1

Use matplotlib.nxutils.points_inside_poly, que implementa una prueba muy eficiente.

Ejemplos y explicación adicional de este algoritmo de 40 años en el matplotlib FAQ.

Actualización: Tenga en cuenta que points_inside_poly está en desuso desde la versión 1.2.0 de matplotlib. Use matplotlib.path.Path.contains_points en su lugar.

Answer 2

Me temo que no estoy familiarizado con las bibliotecas que estás utilizando, pero creo que tengo una idea razonable del algoritmo que podrías utilizar y voy a ver cómo implementarlo con vanilla python y entonces estoy seguro de que puedes mejorarlo e implementarlo usando estas bibliotecas. Además, no estoy afirmando que esta es la mejor manera de lograr esto, pero quería obtener mi respuesta con bastante rapidez así que aquí va.

---

Ahora, la idea proviene de la utilización de la producto cruzado de dos vectores en los algoritmos para encontrar el conjunto convexo de un conjunto de puntos, por ejemplo, Graham's Scan. Digamos que tenemos dos puntos p1 y p2, que definen los vectores de puntos p1 y p2, desde el origen (0,0) hasta (x1, y1) y (x2, y2) respectivamente. El producto cruzado de p1 x p2 da un tercer vector p3 que es perpendicular a ambos p1 y p2 y tiene magnitud dada por el área del paralelogramo delimitada por los vectores.

Un resultado muy útil es que el determinante de la matriz

/ x1, x2 \ 
\ y1, y2/

... que es x1 * y2 - x2 * y1 da la magnitud del vector p3 y el signo indica si p3 está "saliendo del avión" o "entrando" en él. El punto crucial aquí es que si esta magnitud es positiva, entonces p2 es "a la izquierda" de p1 y si es negativa, entonces p2 es "a la derecha" de p1.

suerte, este ejemplo de la técnica ascii ayudará:

. p2(4, 5) 
/
/
/
/_ _ _ _ _. p1(5, 0) 

x1 * y2 - x2 * Y1 = 5 * 4 - 0 * 5 = 20 y así p2 es "a la izquierda" de p1

¡Finalmente en por qué esto es útil para nosotros! Si tenemos una lista de vértices del polígono y un conjunto de los otros puntos en el gráfico, para cada borde del polígono podemos obtener el vector de ese borde. También podemos obtener los vectores que unen el vértice de inicio con todos los demás puntos del gráfico y al probar si se encuentran a la izquierda o a la derecha del borde podemos eliminar algunos puntos para cada borde. Todos los que no se eliminan al final del proceso son esos puntos dentro del polígono. De todos modos, vamos a un código para tener más sentido de esto.

---

Obtener una lista de los vértices de su polígono en el orden que visitar ellos si los estuviera dibujando en sentido contrario a las agujas del reloj, por ejemplo, algunos pentágono podría ser:

poly = [(1, 1), (4, 2), (5, 5), (3, 8), (0, 4)]

Obtenga un conjunto que contenga todos los otros puntos en el gráfico, gradualmente eliminaremos los puntos inválidos de este conjunto hasta que los que quedan al final del proceso sean exactamente esos puntos que están dentro del polígono.

points = set(['(3, 0), (10, -2), (3,3), ...])

El bit principal del código en sí mismo es en realidad bastante compacto para el tiempo que me ha llevado a escribir sobre cómo funciona. to_right toma dos tuplas que representan vectores y devuelve True si v2 se encuentra a la derecha de v1. Luego, los bucles simplemente pasan por todos los bordes del polígono y eliminan puntos del conjunto de trabajo si están a la derecha de cualquiera de los bordes.

def to_right(v1, v2): 
    return (v1[0]*v2[1] - v1[1]*v2[0]) < 0 

for i in range(len(poly)): 
    v1 = poly[i-1] 
    v2 = poly[i] 
    for p in points: 
     if(to_right(v2-v1, p-v1)): 
      points.remove(p) 

de edición: Para aclarar, el hecho de que se eliminan si están a la derecha frente a la izquierda está vinculada a la orden en el que se especifican los vértices del polígono. Si estuvieran en el sentido de las agujas del reloj, querría eliminar los puntos de la izquierda en su lugar. No tengo una solución especialmente buena para este problema en este momento.

---

De todos modos, es de esperar que tengo razón de estas cosas y es de alguna ayuda para alguien, incluso si no es la OP. La complejidad asintótica de este algoritmo es O (mn) donde n es el número de puntos en el gráfico ym es el número de vértices del polígono, ya que en el peor caso todos los puntos se encuentran dentro del polígono y debemos verificar cada punto para cada borde, sin quitar ninguno.