Calcular y mostrar un casco convexo en Python

Un casco convexo es el polígono convexo más pequeño que contiene un conjunto de puntos. En Python, ConvexHull de SciPy lo calcula y Matplotlib puede dibujar su contorno; una pequeña implementación propia de Quickhull resulta útil para estudiar el algoritmo.

¿Qué es un casco convexo?

Para puntos en un plano, el casco convexo es el polígono convexo más pequeño que contiene todos los puntos. Imagina una goma elástica tensa alrededor del conjunto: los puntos que toca son los vértices del casco, mientras que los que quedan dentro no forman parte del contorno. Los lados del polígono conectan esos vértices exteriores en orden.

Esta definición se refiere a un conjunto de puntos, no al límite de un objeto que ya es convexo. Los cascos convexos son útiles en geometría, para delimitar colisiones y en el análisis espacial, aunque descartan deliberadamente cómo se distribuyen los puntos en el interior.

Calcular y representar un casco convexo con SciPy

Instala NumPy, SciPy y Matplotlib en el entorno de Python donde ejecutarás el script; por ejemplo, con python -m pip install numpy scipy matplotlib. ConvexHull recibe una matriz cuyas filas son puntos y cuyas dos columnas contienen las coordenadas x e y. Su atributo vertices proporciona los índices de los puntos del contorno en orden; al conectar esa secuencia y repetir el primer punto, se cierra el polígono.

Ejecuta este ejemplo:

import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import ConvexHull

points = np.array([
    [0.0, 0.0], [1.0, 0.0], [2.0, 0.0],
    [2.0, 1.0], [2.0, 2.0], [1.0, 2.0],
    [0.0, 2.0], [0.0, 1.0], [1.0, 1.0],
])
hull = ConvexHull(points)
boundary = points[hull.vertices]
closed_boundary = np.vstack((boundary, boundary[0]))

plt.scatter(points[:, 0], points[:, 1], label="Points")
plt.plot(closed_boundary[:, 0], closed_boundary[:, 1], "r-", label="Hull")
plt.gca().set_aspect("equal", adjustable="box")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

El punto interior forma parte de los datos de entrada, pero no del casco. En este caso, es preferible usar hull.vertices en vez de representar hull.simplices arista por arista: los vértices ordenados forman un único contorno cerrado. El ejemplo abre una ventana de Matplotlib; en un cuaderno, la figura aparece en el área de salida. Las etiquetas del gráfico dentro del ejemplo inmutable permanecen en inglés.

Implementar Quickhull en Python

Quickhull es un método recursivo para encontrar un casco convexo en 2D. Empieza con los puntos que tienen las coordenadas x menor y mayor. Esos extremos dividen el contorno candidato en dos lados. Para cada lado dirigido, el algoritmo selecciona el punto con el mayor producto vectorial positivo (su valor es proporcional a la distancia perpendicular máxima) y después busca en los dos subconjuntos más pequeños situados a cada lado de ese punto. Si ningún candidato queda a la izquierda de un segmento, ese segmento forma parte del contorno. La prueba estricta de lado excluye cada punto seleccionado de las llamadas recursivas posteriores, por lo que la recursión termina para un conjunto finito de puntos distintos.

La implementación siguiente acepta un iterable de pares de coordenadas 2D y devuelve los vértices del polígono en el orden del recorrido, sin repetir el primero. Elimina las coordenadas repetidas antes de iniciar la búsqueda recursiva. Si hay menos de tres puntos distintos o todos están en una misma línea, devuelve los extremos disponibles. El ejemplo solo usa la biblioteca estándar de Python, así que puedes ejecutarlo como un único script y pasar después la lista devuelta a una biblioteca de gráficos.

Ejecuta este ejemplo:

def cross(a, b, point):
    """Return the signed 2D cross product of AB and A-to-point."""
    return (b[0] - a[0]) * (point[1] - a[1]) - (b[1] - a[1]) * (point[0] - a[0])


def quickhull(points):
    """Return the vertices of a 2D convex hull in boundary order."""
    unique = sorted(set(tuple(point) for point in points))
    if len(unique) <= 2:
        return unique

    left = min(unique, key=lambda point: (point[0], point[1]))
    right = max(unique, key=lambda point: (point[0], point[1]))

    def chain(start, end, candidates):
        candidates = [point for point in candidates if cross(start, end, point) > 0]
        if not candidates:
            return []
        farthest = max(candidates, key=lambda point: cross(start, end, point))
        first = [point for point in candidates if cross(start, farthest, point) > 0]
        second = [point for point in candidates if cross(farthest, end, point) > 0]
        return chain(start, farthest, first) + [farthest] + chain(farthest, end, second)

    upper = chain(left, right, unique)
    lower = chain(right, left, unique)
    return [left] + upper + [right] + lower


points = [
    (0, 0), (1, 0), (2, 0), (2, 1), (2, 2),
    (1, 2), (0, 2), (0, 1), (1, 1),
]
hull_points = quickhull(points)
closed_hull = hull_points + hull_points[:1]
print(hull_points)

El signo del producto vectorial indica en qué lado de una recta dirigida se encuentra un punto; su magnitud permite ordenar los candidatos de ese lado. Las llamadas recursivas conservan únicamente los puntos que quedan estrictamente fuera de los dos segmentos nuevos. La lista devuelta sigue el contorno del polígono, así que añade el primer punto al final para dibujar una línea cerrada. Esta implementación educativa está pensada para puntos 2D y no sustituye a una biblioteca geométrica probada en producción. Además, su recursión directa puede requerir bastante trabajo con ciertas distribuciones de puntos.

Elegir un método y tratar los casos límite

Usa SciPy si necesitas una rutina de biblioteca mantenida, los índices del casco o compatibilidad con más de dos dimensiones. Usa Quickhull propio si quieres entender la partición recursiva o adaptar una rutina pequeña para 2D. No elijas suponiendo que será más rápido: el rendimiento depende de los datos y del entorno, y este artículo no presenta resultados de benchmarks.

La forma y la geometría de los datos importan. Para formar un polígono 2D, SciPy necesita al menos tres puntos no colineales; con muy pocos puntos o con puntos colineales no se puede definir ese polígono, y Qhull informa de un error geométrico. Si tu aplicación debe gestionar esa situación, captura scipy.spatial.QhullError o valida los datos antes de llamar a ConvexHull. Los puntos duplicados no añaden vértices al contorno. La función propia anterior elimina duplicados y devuelve los extremos disponibles si los puntos son colineales, pero presupone pares de coordenadas numéricas finitas; valida los datos que no sean de confianza antes de pasárselos. Para cascos en 3D o en dimensiones superiores, utiliza la implementación de biblioteca de SciPy en vez de este algoritmo propio para 2D.

Preguntas frecuentes

¿Cómo encuentra Quickhull el contorno?

Empieza por los puntos extremos, selecciona el candidato más alejado en cada lado y examina recursivamente los conjuntos de puntos que quedan fuera de los segmentos nuevos. Si un lado no tiene candidatos restantes, aporta una arista al casco.

¿Se puede calcular un casco convexo en 3D?

Sí. ConvexHull de SciPy admite matrices de puntos con más de dos columnas de coordenadas. La implementación propia de este artículo usa productos vectoriales 2D y no admite datos 3D.

¿Qué ocurre si los puntos son colineales?

No forman un polígono 2D. SciPy informa de un error geométrico de Qhull; la función propia devuelve los dos extremos (o menos si hay menos puntos distintos).

¿Qué método debería usar?

Elige SciPy para código de aplicación y datos de dimensiones superiores. Elige la versión propia para estudiar el algoritmo 2D; antes de adaptarla a un sistema más grande, añade validación de entradas y pruebas.