Triangulación de Delaunay — la regla del circuncírculo vacío e inserción Bowyer–Watson

En 1934, Boris Delaunay definió una triangulación en honor a su mentor Georgy Voronoi. La regla cabía en una sola línea: ningún otro punto puede quedar dentro del circuncírculo de un triángulo. Esa única condición se volvió el corazón de casi todo generador de mallas no estructuradas que se usa hoy.

Este artículo muestra por qué esa regla del "circuncírculo vacío" produce buenas mallas. Después construimos una malla punto a punto con el algoritmo de Bowyer–Watson en 40 líneas de Python. Al final vemos cómo el punto flotante traiciona en silencio esta regla tan elegante.

Una malla definida por un solo círculo vacío#

Hay infinitas maneras de conectar un conjunto de puntos en triángulos. Delaunay elige exactamente una: la triangulación en la que el circuncírculo de cada triángulo (el círculo que pasa por sus tres vértices) está vacío.

Sea $C(T)$ el circuncírculo del triángulo $T = \{a, b, c\}$ . La condición se escribe así:

\forall\, p \notin T : \quad p \notin \mathrm{int}\, C(T)

Aquí $\mathrm{int}\,C(T)$ es el interior del circuncírculo y $p$ es cualquier otro punto de la malla. Ningún punto puede quedar dentro del circuncírculo de otro triángulo.

¿Por qué esta regla? Entre todas las triangulaciones de un conjunto de puntos, la de Delaunay maximiza el ángulo mínimo. Es decir, evita en lo posible los triángulos largos y delgados (slivers). Los triángulos finos arruinan la reconstrucción de gradientes y el número de condición de la matriz en métodos de volúmenes y elementos finitos. Por eso los generadores de mallas toman Delaunay como opción por defecto.

Como bono, el dual de una triangulación de Delaunay es justamente el diagrama de Voronoi. Al obtener una estructura, la otra llega gratis.

El determinante InCircle — probar "dentro" de un solo golpe#

Se puede responder "¿está el punto dentro del circuncírculo?" sin calcular nunca el centro ni el radio del círculo. Para un triángulo $(a, b, c)$ ordenado en sentido antihorario, el punto $d$ está dentro del circuncírculo si y solo si este determinante es positivo:

\det \begin{pmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \end{pmatrix} > 0

Cada fila traslada un vértice a coordenadas relativas a $d$ y luego añade su distancia al cuadrado. El signo de este valor depende de la orientación del triángulo: se invierte con un giro horario, así que el código real lo normaliza multiplicando por el signo de orientación.

En la figura de abajo, arrastra el vértice B hacia arriba con el control deslizante.

vertex B height (drag B up toward the diagonal)

The circle is the circumcircle of triangle A–B–C. Slide B upward: when D crosses into the circle the diagonal A–C becomes illegal and the Delaunay solution flips to B–D. Notice the larger of the two min-angle numbers always belongs to the chosen diagonal — Delaunay maximizes the smallest angle.

A medida que B sube, el circuncírculo del triángulo A–B–C crece. En el instante en que el punto D entra en él, la diagonal A–C se vuelve ilegal y la solución de Delaunay salta a B–D. Observa también que el mayor de los dos números de min-angle siempre corresponde a la diagonal elegida.

Bowyer–Watson — suelta un punto y rellena el hueco#

Ahora construimos la malla. Bowyer–Watson es un algoritmo incremental que inserta puntos uno a uno y mantiene la propiedad de Delaunay en cada paso. El flujo tiene cinco pasos.

Empezar con un supertriángulo lo bastante grande como para contener todos los puntos.
Para un punto nuevo $p$ , encontrar todos los triángulos cuyo circuncírculo contiene a $p$ — los llamamos triángulos malos.
Los triángulos malos siempre forman una sola cavidad en forma de estrella (star-shaped). Su frontera es el conjunto de aristas que pertenecen a exactamente un triángulo malo.
Borrar los triángulos malos y conectar $p$ con cada arista de la frontera, formando triángulos nuevos.
Tras insertar todos los puntos, quitar cualquier triángulo que comparta un vértice con el supertriángulo.

La clave es la propiedad del paso 3. Al insertar $p$ , los únicos triángulos que rompen la condición del circuncírculo vacío son "los triángulos cuyo circuncírculo contiene a $p$ ", y esos siempre forman una cavidad conexa. Basta con rellenar esa cavidad para que toda la malla vuelva a ser Delaunay.

Haz clic en el lienzo de abajo para agregar puntos por tu cuenta.

show circumcircles

Each thin cyan ring is the circumcircle of one triangle. In a valid Delaunay mesh no vertex (yellow dot) ever sits inside another triangle's ring — that is the empty-circumcircle property. Add points and watch the mesh repair itself.

Cada vez que se agrega un punto, la malla se reorganiza sola. Los anillos cian delgados son los circuncírculos. Si dejas activado show circumcircles, se ve directamente la propiedad del circuncírculo vacío: ningún punto amarillo cae jamás dentro del anillo de otro triángulo.

Volteo de Lawson — la elección local de una diagonal#

Donde Bowyer–Watson rellena toda una cavidad, el método de Lawson arregla una arista a la vez. Divide en tres el triángulo que contiene al punto nuevo y luego aplica la prueba in-circle a cada arista afectada. Si el punto opuesto cae dentro del circuncírculo, voltea la arista compartida (la diagonal). Repite hasta que no quede ninguna arista ilegal.

Un solo volteo siempre aumenta el ángulo mínimo de los dos triángulos implicados. Así que cuando el bucle termina, queda automáticamente una malla de Delaunay con el ángulo mínimo maximizado. El cambio A–C ↔ B–D que viste arriba es exactamente uno de esos volteos.

La misma herramienta resuelve la Delaunay con restricciones (constrained Delaunay), cuando hay que conservar una arista de frontera. Encuentra la tubería (pipe) de triángulos que cruza una arista de recuperación AB y vuelve a triangular esa franja intercambiando diagonales hasta que AB reaparezca como arista de la malla.

Python — 60 puntos aleatorios en una malla de Delaunay#

Toda la inserción incremental llevada a numpy. Corre sin ninguna biblioteca externa.

import numpy as np
 
def in_circle(a, b, c, d):
    """True si el punto d está estrictamente dentro del circuncírculo de a,b,c.
    Normalizado por el signo de orientación, así el winding no importa."""
    ax, ay = a[0] - d[0], a[1] - d[1]
    bx, by = b[0] - d[0], b[1] - d[1]
    cx, cy = c[0] - d[0], c[1] - d[1]
    det = np.linalg.det(np.array([
        [ax, ay, ax * ax + ay * ay],
        [bx, by, bx * bx + by * by],
        [cx, cy, cx * cx + cy * cy],
    ]))
    orient = (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])
    return det * np.sign(orient) > 1e-12
 
def super_triangle(points):
    pts = np.asarray(points)
    cx, cy = pts.mean(axis=0)
    r = np.max(np.linalg.norm(pts - [cx, cy], axis=1)) + 1.0
    return [(cx - 3 * r, cy - r), (cx + 3 * r, cy - r), (cx, cy + 3 * r)]
 
def bowyer_watson(points):
    st = super_triangle(points)
    tris = [(st[0], st[1], st[2])]
    for p in points:
        bad = [t for t in tris if in_circle(t[0], t[1], t[2], p)]
        edge_n = {}                              # cuenta cuántas veces aparece cada arista
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = frozenset((e[0], e[1]))
                edge_n[key] = edge_n.get(key, 0) + 1
        boundary = [tuple(k) for k, n in edge_n.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        for e in boundary:                       # rellena la cavidad con el punto p
            tris.append((e[0], e[1], p))
    s = set(st)
    return [t for t in tris if not (set(t) & s)]
 
rng = np.random.default_rng(7)
pts = [tuple(p) for p in rng.random((60, 2))]
mesh = bowyer_watson(pts)
print(f"points = {len(pts)}   triangles = {len(mesh)}")

La salida es, aproximadamente:

points = 60   triangles = 110

El número de triángulos no es casualidad. En una malla de Delaunay plana la cantidad de triángulos queda fija en $2N - 2 - h$ , donde $N$ es el número de puntos y $h$ el número de puntos sobre la envolvente convexa. Con $N = 60$ y 8 puntos en la envolvente, salen exactamente 110. Se puede verificar el conteo incluso antes de dibujar la malla.

Cuando los puntos casi caen sobre un mismo círculo#

La regla elegante tiene una trampa. Cuando cuatro puntos caen exactamente sobre un mismo círculo (cocirculares), el determinante in-circle vale cero. Cualquiera de las dos diagonales es Delaunay y la triangulación no es única. El peor caso es una malla cuadrada: las cuatro esquinas de cada cuadrado son cocirculares, así que la entrada más regular resulta ser la más ambigua.

El problema de verdad surge cerca de ese cero, por culpa del punto flotante. Cuando el determinante alcanza el límite de precisión de máquina, su signo se invierte al azar, la topología cae en contradicción y el código se queda en un bucle infinito o se cae. Por eso bibliotecas de producción como Triangle y CGAL usan predicados exactos adaptativos (adaptive exact predicates, Shewchuk) en lugar de un simple double. Calculan solo los dígitos necesarios para acertar el signo siempre.

En tres dimensiones aparece una trampa más. La condición de circunesfera vacía no prohíbe tetraedros casi planos (slivers). Por eso el mallado 3D agrega refinamiento de Delaunay (Ruppert, Chew), insertando puntos de Steiner para subir la calidad.

La última trampa práctica es la frontera. Una malla de Delaunay cruda no respeta la frontera geométrica que dibujaste. Para forzar que una arista de frontera sobreviva hace falta el paso de Delaunay con restricciones de arriba. "¿Por qué mi malla atraviesa la pared?" es la señal de que lo omitiste.

Tres líneas para anotar antes de tender una malla#

Delaunay = circuncírculo vacío = ángulo mínimo maximizado. Las tres formulaciones son lo mismo, todas apuntan a evitar los slivers.
Bowyer–Watson cabe en cinco líneas: "borra los triángulos malos, rellena el hueco con el punto". El truco es el teorema de que el hueco siempre tiene forma de estrella.
La entrada cocircular y el punto flotante rompen el signo de in-circle. Si hay una malla cuadrada de por medio, primero consigue un predicado exacto.