Geometría computacional en Python: de la teoría a la aplicación

Cuando las personas piensan en geometría computacional, en mi experiencia, generalmente piensan en una de dos cosas:

1. Vaya, eso suena complicado.
2. Oh sí, casco convexo.

En esta publicación, me gustaría arrojar algo de luz sobre la geometría computacional, comenzando con una breve descripción general del tema antes de pasar a algunos consejos prácticos basados ​​en mis propias experiencias (sáltese si tiene un buen manejo del tema).

¿Qué es todo este alboroto?

Si bien los algoritmos de geometría computacional de casco convexo generalmente se incluyen en un curso introductorio de algoritmos, la geometría computacional es un tema mucho más rico que rara vez recibe suficiente atención del desarrollador / científico informático promedio (a menos que esté creando juegos o algo así).

Teóricamente intrigante…

Desde un punto de vista teórico, las cuestiones de la geometría computacional suelen ser sumamente interesantes; las respuestas, convincentes; y los caminos por los que se llega a ellos, variados. Estas cualidades por sí solas lo convierten en un campo que vale la pena estudiar, en mi opinión.

Por ejemplo, considere el problema de la galería de arte: somos dueños de una galería de arte y queremos instalar cámaras de seguridad para proteger nuestra obra de arte. Pero tenemos un presupuesto ajustado, por lo que queremos usar la menor cantidad de cámaras posible. ¿Cuántas cámaras necesitamos?

Cuando traducimos esto a notación geométrica computacional, el 'plano de planta' de la galería es solo un polígono simple. Y con un poco de esfuerzo, podemos demostrar que n/3 cámaras siempre son suficientes para un polígono en n vértices, sin importar lo complicado que sea. La prueba en sí usa gráficos duales, algo de teoría de gráficos, triangulaciones y más.

Aquí vemos una técnica de prueba inteligente y un resultado que es lo suficientemente curioso como para apreciarlo por sí solo. Pero si la relevancia teórica no es suficiente para ti…

E importante en la práctica

Como mencioné anteriormente, el desarrollo de juegos se basa en gran medida en la aplicación de geometría computacional (por ejemplo, la detección de colisiones a menudo se basa en calcular el casco convexo de un conjunto de objetos); al igual que los sistemas de información geográfica (SIG), que se utilizan para almacenar y realizar cálculos sobre datos geográficos; y la robótica también (p. ej., para problemas de visibilidad y planificación).

¿Por qué es tan difícil?

Tomemos un problema de geometría computacional bastante sencillo: dado un punto y un polígono, ¿el punto se encuentra dentro del polígono? (Esto se llama el problema del punto en el polígono o PIP).

PIP hace un gran trabajo al demostrar por qué la geometría computacional puede ser (engañosamente) difícil. Para el ojo humano, esta no es una pregunta difícil. Vemos el siguiente diagrama y es inmediatamente obvio para nosotros que el punto está en el polígono:

Incluso para polígonos relativamente complicados, la respuesta no se nos escapa por más de uno o dos segundos. Pero cuando alimentamos este problema a una computadora, es posible que vea lo siguiente:

 poly = Polygon([Point(0, 5), Point(1, 1), Point(3, 0), Point(7, 2), Point(7, 6), Point(2, 7)]) point = Point(5.5, 2.5) poly.contains(point)

Lo que es intuitivo para el cerebro humano no se traduce tan fácilmente al lenguaje informático.

De manera más abstracta (e ignorando la necesidad de representar estas cosas en código), los problemas que vemos en esta disciplina son muy difíciles de rigorizar ('hacer rigurosos') en un algoritmo de geometría computacional. ¿Cómo describiríamos el escenario de un punto en un polígono sin usar un lenguaje tan tautológico como 'Un punto está dentro de un polígono si está dentro del polígono'? Muchas de estas propiedades son tan fundamentales y tan básicas que es difícil definirlas concretamente.

Difícil, pero no imposible. Por ejemplo, podría endurecer el punto en el polígono con las siguientes definiciones:

• Un punto está dentro de un polígono si cualquier rayo infinito que comienza en el punto se cruza con un número impar de aristas del polígono (conocida como la regla par-impar).
• Un punto está dentro de un polígono si tiene un número de vueltas distinto de cero (definido como el número de veces que la curva que define el polígono se desplaza alrededor del punto).

A menos que haya tenido alguna experiencia con la geometría computacional, estas definiciones probablemente no serán parte de su vocabulario existente. Y tal vez eso es emblemático de cómo la geometría computacional puede empujarlo a pensar de manera diferente .

Presentamos CCW

Ahora que tenemos una idea de la importancia y la dificultad de los problemas de geometría computacional, es hora de mojarnos las manos.

En la columna vertebral del tema hay una operación primitiva engañosamente poderosa: en sentido contrario a las agujas del reloj, o 'CCW' para abreviar. (Te lo advierto ahora: CCW aparecerá una y otra vez).

CCW toma tres puntos A, B y C como argumentos y pregunta: ¿estos tres puntos componen un giro en sentido contrario a las agujas del reloj (frente a un giro en el sentido de las agujas del reloj)? En otras palabras, ¿es A -> B -> C un ángulo en sentido antihorario?

Por ejemplo, los puntos verdes son CCW, mientras que los puntos rojos no lo son:

Por qué importa CCW

CCW nos da una operación primitiva sobre la cual podemos construir. Nos da un lugar para comenzar a endurecer y resolver problemas de geometría computacional.

Para darle una idea de su poder, consideremos dos ejemplos.

Determinación de la convexidad

La primera: dado un polígono, ¿puedes determinar si es convexo? La convexidad es una propiedad invaluable: saber que sus polígonos son convexos a menudo le permite mejorar el rendimiento en órdenes de magnitud. Como ejemplo concreto: hay un algoritmo PIP bastante sencillo que se ejecuta en tiempo Log(n) para polígonos convexos, pero falla para muchos polígonos cóncavos.

Intuitivamente, esta brecha tiene sentido: las formas convexas son 'agradables', mientras que las formas cóncavas pueden tener bordes afilados que sobresalen hacia adentro y hacia afuera, simplemente no siguen las mismas reglas.

Un algoritmo de geometría computacional simple (pero no obvio) para determinar la convexidad es verificar que cada triplete de vértices consecutivos sea CCW. Esto requiere solo unas pocas líneas de código de geometría de Python (suponiendo que los points se proporcionen en el sentido contrario a las agujas del reloj; si los points están en el sentido de las agujas del reloj, querrá que todos los tripletes estén en el sentido de las agujas del reloj):

 class Polygon(object): ... def isConvex(self): for i in range(self.n): # Check every triplet of points A = self.points[i % self.n] B = self.points[(i + 1) % self.n] C = self.points[(i + 2) % self.n] if not ccw(A, B, C): return False return True

Pruebe esto en papel con algunos ejemplos. Incluso puede usar este resultado para definir la convexidad. (Para hacer las cosas más intuitivas, tenga en cuenta que una curva CCW de A -> B -> C corresponde a un ángulo de menos de 180, que es una forma ampliamente enseñada de definir la convexidad).

Intersección de línea

Como segundo ejemplo, considere la intersección de un segmento de línea, que también se puede resolver usando solo CCW:

 def intersect(a1, b1, a2, b2): """Returns True if line segments a1b1 and a2b2 intersect.""" return ccw(a1, b1, a2) != ccw(a1, b1, b2) and ccw(a2, b2, a1) != ccw(a2, b2, b1)

¿Por qué es este el caso? La intersección del segmento de línea también se puede expresar como: dado un segmento con extremos A y B, ¿los extremos C y D de otro segmento se encuentran en el mismo lado de AB? En otras palabras, si los giros de A -> B -> C y A -> B -> D están en la misma dirección, los segmentos no pueden intersecarse. Cuando usamos este tipo de lenguaje, queda claro que tal problema es el pan y la mantequilla de CCW.

Una definición rigurosa

Ahora que tenemos una idea de la importancia de CCW, veamos cómo se calcula. Dados los puntos A, B y C:

 def ccw(A, B, C): """Tests whether the turn formed by A, B, and C is ccw""" return (Bx - Ax) * (Cy - Ay) > (By - Ay) * (Cx - Ax)

Para entender de dónde viene esta definición, considere los vectores AB y BC. Si tomamos su producto cruzado, AB x BC, este será un vector a lo largo del eje z. Pero, ¿en qué dirección (es decir, +z o -z)? Resulta que, si el producto vectorial es positivo, el giro es en sentido contrario a las agujas del reloj; de lo contrario, es en el sentido de las agujas del reloj.

Esta definición parecerá poco intuitiva a menos que tenga una muy buena comprensión del álgebra lineal, la regla de la mano derecha, etc. Pero es por eso que tenemos abstracción: cuando piense en CCW, solo piense en su definición intuitiva en lugar de su cálculo. El valor será inmediatamente claro.

Mi inmersión en la geometría computacional y la programación usando Python

Durante el último mes, he estado trabajando en la implementación de varios algoritmos de geometría computacional en Python. Como los dibujaré en las próximas secciones, me tomaré un segundo para describir mis aplicaciones de geometría computacional, que se pueden encontrar en GitHub.

Nota: Mi experiencia es ciertamente limitada. Como he estado trabajando en esto durante meses en lugar de años, tome mi consejo con pinzas. Dicho esto, aprendí mucho en esos pocos meses, así que espero que estos consejos sean útiles.

Algoritmo de Kirkpatrick

El núcleo de mi trabajo fue una implementación del Algoritmo de Kirkpatrick para la ubicación de puntos. El enunciado del problema sería algo como: dada una subdivisión plana (un grupo de polígonos que no se superponen en el plano) y un punto P, ¿qué polígono contiene P? Piense en un punto en un polígono con esteroides: en lugar de un solo polígono, tiene un avión lleno de ellos.

Como caso de uso, considere una página web. Cuando un usuario hace clic con el mouse, la página web necesita averiguar en qué hizo clic el usuario lo más rápido posible. ¿Era el botón A? ¿Fue el enlace B? La página web está compuesta de polígonos que no se superponen, por lo que el Algoritmo de Kirkpatrick estaría bien posicionado para ayudar.

Si bien no discutiré el algoritmo en profundidad, puede obtener más información aquí.

Triángulo delimitador mínimo

Como subtarea, también implementé el algoritmo de O'Rourke para calcular un triángulo envolvente/limitador mínimo (es decir, encontrar el triángulo más pequeño que encierra un conjunto de puntos convexos) en tiempo lineal.

Nota: calcular el triángulo delimitador mínimo no ayuda ni perjudica el rendimiento asintótico del algoritmo de Kirkpatrick, ya que el cálculo en sí es de tiempo lineal, pero es útil para fines estéticos.

Consejos prácticos, aplicaciones e inquietudes

Las secciones anteriores se centraron en por qué la geometría computacional puede ser difícil de razonar rigurosamente.

En la práctica, tenemos que lidiar con toda una nueva serie de preocupaciones.

¿Recuerdas CCW? Como buena continuación, veamos otra de sus grandes cualidades: nos protege contra los peligros de los errores de coma flotante.

Errores de coma flotante: por qué CCW es el rey

En mi curso de geometría computacional, Bernard Chazelle, un estimado profesor que ha publicado más artículos de los que puedo contar, estableció como regla que no podíamos mencionar ángulos al intentar describir un algoritmo o una solución.

¿Por qué? Los ángulos están desordenados. Los ángulos están "sucios". Cuando tiene que calcular un ángulo, necesita dividir o usar alguna aproximación (cualquier cosa que involucre a Pi, por ejemplo) o alguna función trigonométrica.

Cuando tenga que calcular un ángulo en el código , casi siempre estará aproximando. Estará equivocado por un pequeño grado de precisión de punto flotante, lo cual es importante cuando está probando la igualdad. Puede resolver algún punto en el plano a través de dos métodos diferentes y, por supuesto, esperar que p1.x == p2.x and p1.y == p2.y . Pero, en realidad, esta verificación fallará a menudo . Además (y obviamente), estos puntos tendrán diferentes hashes.

Para empeorar las cosas, su grado de error aumentará a medida que sus pequeñas diferencias se propaguen a través de sus cálculos. (Para ver algunos ejemplos más científicos, este documento analiza lo que puede salir mal al calcular el casco convexo o la triangulación de Delaunay).

¿Entonces, qué podemos hacer al respecto?

casi igual

Parte del problema de la geometría computacional de Python es que requerimos exactitud en un mundo donde las cosas rara vez son exactas. Esto se convertirá en un problema más a menudo que cuando se manejan ángulos. Considera lo siguiente:

 # Define two random points p1 = RandomPoint() p2 = RandomPoint() # Take the line through them l1 = Line(p1, p2) # Shift both points up by sqrt(2) p1.y += sqrt(2) p2.y += sqrt(2) l2 = Line(p1, p2) # Slope 'should' be the same? if abs(l1.slope - l2.slope) > 0: print "Error!" # Error!

De hecho, este código imprimirá "¡Error!" aproximadamente el 70% del tiempo (empíricamente). Podemos abordar esta preocupación siendo un poco más indulgentes con nuestra definición de igualdad; es decir, sacrificando un grado de precisión.

Un enfoque que he usado (y visto, por ejemplo, en algunos módulos de OpenCV) es definir dos números como iguales si difieren solo en algún épsilon de valor pequeño. En Python, es posible que tenga:

 def almostEqual(x, y, EPSILON=1e-5): return abs(x - y) < EPSILON class Point(object): ... def __eq__(self, that): return (almostEqual(self.x, that.x) and almostEqual(self.y, that.y))

En la práctica, esto es muy útil. Rara vez, si acaso, calcularía dos puntos que difieren en menos de 1e-5 que en realidad se supone que son puntos diferentes. Recomiendo encarecidamente implementar este tipo de anulación. Se pueden usar métodos similares para las líneas, por ejemplo:

 class Line(object): ... def __eq__(self, that): return (almostEqual(self.slope, that.slope) and almostEqual(self.intercept, that.intercept))

Por supuesto, se han propuesto soluciones más avanzadas. Por ejemplo, la escuela de pensamiento del 'cómputo geométrico exacto' (descrita en este documento) tiene como objetivo que todas las rutas de decisión en un programa dependan únicamente del signo de algún cómputo, en lugar de su valor numérico exacto, eliminando muchas de las preocupaciones relacionadas a los cálculos de punto flotante. Nuestro enfoque de casi igualdad solo araña la superficie, pero a menudo será suficiente en la práctica.

CCW es el rey

En un nivel superior, es (posiblemente) problemático que incluso definamos nuestras soluciones en términos de cantidades computacionales tan exactas como ángulos o coordenadas de puntos. En lugar de abordar solo los síntomas (es decir, lavar los errores de coma flotante con almostEqual ), ¿por qué no abordar la causa? La solución: en lugar de pensar en términos de ángulos, piense en términos de CCW , lo que ayudará a abstraer las preocupaciones asociadas con el cálculo de coma flotante.

Aquí hay un ejemplo concreto: digamos que tienes un polígono P convexo, un vértice v y un punto u fuera del polígono. ¿Cómo puedes averiguar si la línea uv interseca a P por encima o por debajo de v , o no lo hace en absoluto, en tiempo constante?

La solución de fuerza bruta (además de ser de tiempo lineal, en lugar de constante) sería problemática ya que tendría que calcular algunos puntos de intersección de línea exactos.

Un enfoque de tiempo constante que he visto implica:

• Cálculo de algunos ángulos usando arctan2 .
• Convertir estos ángulos a grados multiplicando por 180/Pi.
• Examinar las relaciones entre estos diversos ángulos.

Afortunadamente, el autor usó la técnica almostEqual anterior para suavizar los errores de punto flotante.

En mi opinión, sería mejor evitar por completo el problema de los errores de coma flotante. Si se toma unos minutos para analizar el problema en papel, puede obtener una solución basada completamente en CCW. La intuición: si los vértices adyacentes a v están del mismo lado de uv , entonces la línea no se corta; de lo contrario, ver si u y v están en el mismo lado de la línea entre los vértices adyacentes y, según el resultado, comparar sus alturas.

Aquí está el código de Python para probar la intersección por encima de v (la intersección por debajo simplemente invierte la dirección de las comparaciones):

 def intersectsAbove(verts, v, u): """ Returns True if uv intersects the polygon defined by 'verts' above v. Assumes v is the index of a vertex in 'verts', and u is outside of the polygon. """ n = len(verts) # Test if two adjacent vertices are on same side of line (implies # tangency) if ccw(u, verts[v], verts[(v - 1) % n]) == ccw(u, verts[v], verts[(v + 1) % n]): return False # Test if u and v are on same side of line from adjacent # vertices if ccw(verts[(v - 1) % n], verts[(v + 1) % n], u) == ccw(verts[(v - 1) % n], verts[(v + 1) % n], verts[v]): return uy > verts[v].y else: return uy < verts[v].y

La solución no es inmediatamente obvia a simple vista, pero está en el lenguaje de un algoritmo de geometría computacional: 'mismo lado de la línea' es un elemento clásico de ese algoritmo confiable.

Hecho es mejor que perfecto

En la literatura geométrica computacional, a menudo hay una buena cantidad de magia involucrada en operaciones aparentemente simples. Esto le da una opción: puede hacer las cosas de la manera más difícil, siguiendo un documento que define una solución increíblemente avanzada para un problema no tan avanzado, o puede hacer las cosas de la manera más fácil con un poco de fuerza bruta.

Nuevamente, usaré un ejemplo: muestrear un punto interior aleatorio de un polígono arbitrario. En otras palabras, te doy un polígono simple y tú me das un punto aleatorio dentro de él (distribuido uniformemente a lo largo del polígono).

A menudo, se requieren puntos interiores para las pruebas. En ese caso, no tiene requisitos de tiempo de ejecución específicos en el algoritmo de geometría computacional que los produce (dentro de lo razonable). La solución rápida y sucia, que tarda ~2 minutos en implementarse, sería elegir un punto aleatorio dentro de un cuadro que contiene el polígono y ver si el punto en sí está dentro del polígono.

Por ejemplo, podemos fallar dos veces y encontrar una muestra válida solo en el tercer punto:

Aquí está el código:

 class Polygon(object): ... def interiorPoint(self): """Returns a random point interior point""" min_x = min([px for p in self.points]) max_x = max([px for p in self.points]) min_y = min([py for p in self.points]) max_y = max([py for p in self.points]) def x(): return min_x + random() * (max_x - min_x) def y(): return min_y + random() * (max_y - min_y) p = Point(x(), y()) while not self.contains(p): p = Point(x(), y()) return p def contains(self, p): for i in range(self.n): p1 = self.points[i] p2 = self.points[(i + 1) % self.n] p3 = self.points[(i + 2) % self.n] if not ccw(p1, p2, p3): return False return True

Esto se conoce como muestreo de rechazo : tome puntos al azar hasta que uno satisfaga sus criterios. Si bien puede requerir varias muestras para encontrar un punto que cumpla con sus criterios, en la práctica, la diferencia será insignificante para su conjunto de pruebas. Entonces, ¿por qué trabajar más duro? En resumen: no tengas miedo de tomar el camino sucio cuando la ocasión lo requiera.

Por cierto: si desea un algoritmo exacto para el muestreo aleatorio, aquí hay uno inteligente que he implementado a continuación. La esencia de esto:

1. Triangula tu polígono (es decir, divídelo en triángulos).
2. Elige un triángulo con probabilidad proporcional a su área.
3. Tome un punto aleatorio dentro del triángulo elegido (una operación de tiempo constante).

Tenga en cuenta que este algoritmo requiere que triangule su polígono, lo que impone inmediatamente un límite de tiempo de ejecución diferente en el algoritmo, así como la necesidad de tener una biblioteca para triangular polígonos arbitrarios (utilicé poly2tri con enlaces de Python).

 from p2t import CDT class Triangle(object): ... def area(self): return abs((Bx * Ay - Ax * By) + (Cx * By - Bx * Cy) + (Ax * Cy - Cx * Ay)) / 2 def interiorPoint(self): r1 = random() r2 = random() # From http://www.cs.princeton.edu/~funk/tog02.pdf return (1 - sqrt(r1)) * A + sqrt(r1) * (1 - r2) * B + r2 * sqrt(r1) * C class Polygon(object): ... def triangulate(self): # Triangulate poly with hole cdt = CDT(poly.points) triangles = cdt.triangulate() def convert(t): A = Point(tax, tay) B = Point(tbx, tby) C = Point(tcx, tcy) return Triangle(A, B, C) return map(convert, triangles) def interiorPoint(self): # Triangulate polygon triangles = self.triangulate() areas = [t.area() for t in triangles] total = sum(areas) # Calculate normalized areas probabilities = [area / total for area in areas] weighted_triangles = zip(triangles, probabilities) # Sample triangles according to area r = random() count = 0 for (triangle, prob) in weighted_triangles: count += prob # Take random point from chosen triangle if count > r: return triangle.interiorPoint()

Con suerte, el esfuerzo adicional es evidente en el código. Recuerda: como dicen en Facebook, “más vale hecho que perfecto”. Lo mismo ocurre con los problemas de geometría computacional.

Pruebas visuales y automatizadas

Dado que muchos de los problemas en los que trabaja en geometría computacional se definen en términos de cualidades o cantidades fácilmente visualizables, las pruebas visuales son particularmente importantes, aunque insuficientes por sí solas. El conjunto de pruebas ideal tendrá una combinación de pruebas automatizadas visuales y aleatorias.

Nuevamente, procedemos con el ejemplo. Considere probar nuestra implementación del Algoritmo de Kirkpatrick. En un paso, el algoritmo necesita delimitar el polígono dado por un triángulo y triangular la región entre el polígono y el triángulo exterior. Aquí hay un ejemplo visual, donde la línea verde continua define el polígono inicial y las líneas discontinuas definen la región triangulada:

Confirmar que esta triangulación se ha ejecutado correctamente es muy difícil de verificar a través del código, pero es inmediatamente evidente para el ojo humano. Nota: Sugiero encarecidamente usar Matplotlib para ayudar en sus pruebas visuales; hay una buena guía aquí.

Posteriormente, querremos comprobar que el algoritmo localiza correctamente los puntos. Un enfoque aleatorio y automatizado sería generar un montón de puntos interiores para cada polígono y asegurarnos de que devolvamos el polígono deseado. En codigo:

 class TestLocator(unittest.TestCase): ... def runLocator(self, polygons): # Pre-process regions l = Locator(polygons) # Ensure correctness for polygon in polygons: # Test 100 random interior points per region for k in range(100): target = polygon.interiorPoint() target_polygon = l.locate(target) self.assertEqual(polygon, target_polygon) self.assertTrue(target_polygon.contains(target))

Entonces podríamos usar el método runLocator en diferentes conjuntos de polígonos, dándonos un conjunto de pruebas bien diversificado.

Soluciones de código abierto

La geometría computacional tiene un buen conjunto de bibliotecas y soluciones de código abierto disponibles independientemente del lenguaje de programación que elija (aunque las bibliotecas de C++ parecen surgir en una cantidad desproporcionada).

Los beneficios de usar soluciones de código abierto existentes (como con la computación científica en Python) son bien conocidos y se han discutido extensamente, por lo que no me extenderé al respecto aquí. Pero pensé en mencionar algunos recursos centrados en Python que encontré útiles:

• poly2tri: una gran biblioteca para triangulaciones rápidas de polígonos. También admite (y esto suele ser crucial) polígonos con agujeros . Escrito en C++, poly2tri también tiene enlaces de Python y fue bastante fácil de poner en marcha. Vea mi método triangulate anterior para probar las llamadas a funciones.
• scipy.spatial: incluye funciones para calcular cascos convexos, triangulaciones de Delaunay y más. Rápido (como siempre), confiable, etc. Nota: me resultó útil usar mi propio tipo de datos Point con un método toNumpy : def np(self): return [self.x, self.y] . Entonces, podría llamar fácilmente a los métodos scipy.spatial, por ejemplo: scipy.spatial.ConvexHull(np.array(map(lambda p: p.np()), points)) .
• OpenCV: la biblioteca de visión artificial de código abierto tiene algunos buenos módulos de geometría computacional independientes. En particular, usé su función de triángulo envolvente mínimo por un tiempo antes de implementarla yo mismo.

Conclusión

Espero que esta publicación le haya dado una idea de la belleza de la geometría computacional como desarrollador de Python, un tema rico en problemas fascinantes y aplicaciones igualmente fascinantes.

En la práctica, las implementaciones geométricas computacionales presentan desafíos únicos que lo empujarán a ejercitar nuevas y emocionantes habilidades para resolver problemas.

Si está interesado en aprender más o tiene alguna pregunta para mí, me puede contactar en [email protected].