Dos filosofías para triangular una nube de puntos — Bowyer

El signo del circuncírculo voltea la diagonal, y el producto de forma y tamaño elige el siguiente triángulo del frente.

En 1934, Borís Delaunay extrajo un dual de un diagrama que llevaba el nombre de su maestro, Georgy Voronói. Unir los puntos en triángulos y dejar cualquier otro punto fuera de cada círculo circunscrito. Han pasado noventa años. Los generadores de malla siguen corriendo con esa regla.

Doce puntos con semilla 3 devuelven 17 triángulos. Antes de leer esa cuenta, conviene ver hacia dónde va la diagonal cuando cuatro puntos forman un cuadrilátero. Al arrastrar el vértice superior hacia arriba y hacia abajo, la diagonal se invierte en el instante en que cruza la arista compartida.

Drag vertex T through the shared edge. When T enters the circumcircle of the lower triangle, the diagonal snaps from LR to BT — a single Lawson flip restoring the empty-circumcircle rule.

El punto a mirar es dónde entra el vértice T en el círculo circunscrito del triángulo inferior. En ese único movimiento la diagonal pasa de LR a BT. Ese movimiento es un giro de Lawson, la edición mínima que restituye en lo local el círculo circunscrito vacío.

Dos filosofías para cubrir una nube de puntos#

Un mismo conjunto de puntos admite muchas triangulaciones. En CFD estorban las finas y afiladas. Un triángulo aplanado agranda el error de las derivadas. También empeora el número de condición del sistema lineal.

La triangulación de Delaunay maximiza el ángulo mínimo entre las opciones. La esquina más afilada queda tan roma como permiten las posiciones. Con solo las posiciones, la conectividad queda fija: volver a unir los mismos puntos después no reordena las aristas.

El frente de avance recorre el trabajo al revés. No dispersa puntos primero. La frontera es el arranque, y los triángulos se pegan hacia adentro de uno en uno. Un punto nuevo aparece solo cuando un triángulo lo necesita. Cada paso es avidez local: tomar el mejor elemento disponible en ese momento.

El círculo circunscrito vacío — la propiedad que define a Delaunay#

Una sola frase fija la malla. Ningún círculo circunscrito de un triángulo contiene a otro vértice en su interior. Es la condición del círculo circunscrito vacío. Que un punto $d$ caiga dentro del círculo circunscrito del triángulo $abc$ se decide con un determinante. Se ordena $abc$ en sentido antihorario y se calcula:

\begin{vmatrix} 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{vmatrix} > 0

Aquí $a,b,c$ son los tres vértices del triángulo y $d$ es el punto a examinar. Un determinante positivo significa que $d$ está dentro del círculo, negativo fuera, y cero justo sobre la circunferencia. Sin división ni raíz cuadrada. Solo importa el signo, así que la prueba es relativamente robusta frente al error de punto flotante. Cómo se rompe y se recupera ese signo en una inserción está en el círculo vacío y la inserción de Bowyer–Watson.

Bowyer–Watson: cavar una cavidad y volver a rellenarla#

Repetir giros de Lawson también construye una malla de Delaunay. El camino más recto es la inserción incremental que Adrian Bowyer y David Watson publicaron uno junto al otro en la misma revista en 1981.

Al colocar un punto nuevo $p$ en una malla que ya es Delaunay, todo triángulo cuyo círculo circunscrito contiene a $p$ queda invalidado. Al borrarlos queda un hueco poligonal, una cavidad. Conectar cada arista de frontera de esa cavidad con $p$ rellena el hueco.

Tres pasos.

Encontrar todos los triángulos (los triángulos malos) cuyo círculo circunscrito contiene al nuevo punto $p$ .
Extraer las aristas de frontera de la cavidad que forman. Una arista de frontera pertenece a un solo triángulo malo.
Borrar los triángulos malos y, por cada arista de frontera, formar un triángulo nuevo con $p$ .

Se arranca desde un enorme supertriángulo que encierra a todos los puntos de entrada. Tras la última inserción, se retira cualquier triángulo que aún toque un vértice del supertriángulo. Solo queda la malla de Delaunay final.

Un clic en el lienzo suelta puntos. Cada inserción borra los triángulos malos y vuelve a rellenar la cavidad.

Show circumcircles

Click on the canvas to insert a point. Each insertion rebuilds the Delaunay triangulation via Bowyer–Watson. Toggle the circumcircles: every one stays empty of other vertices — that is the Delaunay property.

Con los círculos circunscritos activados, ninguno debe contener otro vértice. Aunque los puntos caigan muy densos, la propiedad se mantiene. Eso es Delaunay.

Python — construir Delaunay incremental desde cero#

La prueba del círculo circunscrito y el núcleo de inserción, pasados a Python puro y sin librerías.

import numpy as np
 
def circumcircle(a, b, c):
    """Devuelve el centro y el radio al cuadrado del circuncírculo. None si colineal."""
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        return None
    a2, b2, c2 = ax * ax + ay * ay, bx * bx + by * by, cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return (ux, uy), (ax - ux) ** 2 + (ay - uy) ** 2
 
def in_circumcircle(p, tri, pts):
    """True si el punto p está dentro del circuncírculo del triángulo tri."""
    cc = circumcircle(pts[tri[0]], pts[tri[1]], pts[tri[2]])
    if cc is None:
        return False
    (ux, uy), r2 = cc
    return (p[0] - ux) ** 2 + (p[1] - uy) ** 2 < r2 - 1e-12
 
def bowyer_watson(points):
    """Construye una triangulacion de Delaunay por insercion incremental."""
    pts = list(points)
    xs = [p[0] for p in pts]; ys = [p[1] for p in pts]
    dmax = 20 * max(max(xs) - min(xs), max(ys) - min(ys))
    cx, cy = (min(xs) + max(xs)) / 2, (min(ys) + max(ys)) / 2
    base = len(pts)                                    # indice inicial de vertices del supertriangulo
    pts += [(cx - dmax, cy - dmax), (cx + dmax, cy - dmax), (cx, cy + dmax)]
    tris = [(base, base + 1, base + 2)]
 
    for pi in range(base):                             # insertar puntos de uno en uno
        p = pts[pi]
        bad = [t for t in tris if in_circumcircle(p, t, pts)]
        edge_count = {}                                # frontera de la cavidad = aristas vistas una vez
        for t in bad:
            for e in [(t[0], t[1]), (t[1], t[2]), (t[2], t[0])]:
                key = tuple(sorted(e))
                edge_count[key] = edge_count.get(key, 0) + 1
        boundary = [e for e, n in edge_count.items() if n == 1]
        tris = [t for t in tris if t not in bad]
        tris += [(a, b, pi) for a, b in boundary]      # conectar p a cada arista de frontera
 
    return [t for t in tris if all(i < base for i in t)]  # quitar el supertriangulo
 
if __name__ == "__main__":
    rng = np.random.default_rng(3)
    P = [tuple(xy) for xy in rng.random((12, 2))]
    T = bowyer_watson(P)
    print(f"{len(P)} puntos -> {len(T)} triangulos")
    print("indices de vertices del primer triangulo:", T[0])

La ejecución imprime 12 puntos -> 17 triangulos. El número de triángulos dentro de la envolvente convexa es aproximadamente $2n - 2 - h$ ( $n$ es el número de puntos, $h$ el número de puntos en la envolvente). El núcleo de inserción puro cabe en unas 60 líneas.

Con doce puntos, esa cuenta da $2\times 12-2=22$ . La salida es 17, así que cuando la fórmula es exacta la envolvente lleva $22-17=5$ puntos. La escala del supertriángulo es 20 veces el ancho de la nube. La muestra sale del cuadrado unidad, de modo que el ancho no pasa de 1, y los tres vértices quedan a esa escala del centro de la caja, con los 12 puntos dentro. El 20 es holgura para encerrar. No dice nada sobre lo afilado de un triángulo.

El supertriángulo del código tiene una base de dos veces la escala y una altura de dos veces la escala, así que el área es dos veces la escala al cuadrado. Con un ancho de nube igual a 1, la escala es 20 y el área 800. Frente al cuadrado unidad de área 1 que contiene los puntos, el supertriángulo es 800 veces mayor. Al retirar todo triángulo que aún toca un vértice del supertriángulo, esos 800 no quedan en el resultado. Lo que queda son los 17 triángulos hechos con los 12 puntos.

El determinante de arriba solo lee un signo. El código declara un punto interior solo cuando la distancia al cuadrado es menor que el radio al cuadrado menos 1e-12. Un punto casi sobre la circunferencia queda fuera. Tres puntos casi colineales chocan con el mismo 1e-12 y se abandona el centro. Del lado de la cavidad solo cuenta cuántas veces aparece una arista. Una arista vista dos veces es interior al hueco. Una arista vista una vez pasa a ser la base de un triángulo nuevo, y la malla gana un triángulo por cada una.

Cómo respetar la frontera — restricciones y la tubería#

Hasta aquí bastaban los puntos. Una geometría real de CFD no. La superficie de un perfil o la pared de un cilindro llevan líneas de frontera que deben existir como aristas de la malla. El Delaunay puro puede cortarlas.

La salida es una triangulación de Delaunay con restricciones. Para recuperar una línea de frontera perdida $AB$ , primero se localiza la banda de triángulos que $AB$ atraviesa, la tubería (pipe). Se parte de un triángulo pegado al punto $A$ , se camina por cada vecino que $AB$ perfora y se detiene en $B$ . La tubería queda completa.

Se retriangula el interior para que $AB$ vuelva como arista, por intercambio de diagonales (swapping of diagonals) o por divide y vencerás. La frontera sobrevive. Los triángulos vecinos quizá ya no sean perfectamente Delaunay. Se cambia fidelidad geométrica por calidad de elemento.

El frente de avance — apilar elementos con avidez#

La técnica de frente de avance (AFT) da la vuelta al problema de la frontera. La frontera es la línea de salida.

Se elige un segmento de referencia $AB$ de la lista de frontera. Se busca el punto $C$ que forma el mejor triángulo con $AB$ . Hay dos clases de candidato: un punto que ya está en el frente, o un punto interior nuevo $I$ . Se toma el de mayor calidad, se forma el triángulo y se actualiza el frente. Las aristas nuevas entran al frente. Las que ya estaban se dan por cerradas y se retiran.

La calidad es el producto de dos factores.

\lambda = \alpha \cdot \delta

Aquí $\alpha$ es el factor de forma del triángulo (1 para el equilátero, hacia 0 cuanto más se aplana) y $\delta$ mide cuánto coincide la longitud real de la arista con el tamaño de elemento pedido (1 cuando concuerdan). Cuanto mayor $\lambda$ , mejor el elemento.

El factor de forma se codifica así.

import math
 
def shape_quality(a, b, c):
    """Factor de forma alpha: 1 para equilatero, hacia 0 para triangulos aplanados."""
    l1 = math.dist(a, b); l2 = math.dist(b, c); l3 = math.dist(c, a)
    area = abs((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2
    denom = l1 * l1 + l2 * l2 + l3 * l3
    return 4 * math.sqrt(3) * area / denom if denom > 0 else 0.0

Cuando el frente se vacía, la frontera cerrada ha convergido hacia adentro y la malla queda terminada. Cada paso toma un óptimo local, así que la malla inicial de AFT suele salir más limpia que una de Delaunay. La diferencia es mayor en mallas muy direccionales, en particular en capas límite.

En un triángulo equilátero el área es el lado al cuadrado por $\sqrt{3}/4$ , y la suma de los lados al cuadrado es tres veces el lado al cuadrado. Al multiplicar el área por el $4\sqrt{3}$ del código y dividir por esa suma se obtiene $4\sqrt{3}\cdot(\sqrt{3}/4)/(3)=1$ . Si un lado se alarga y la altura baja, el numerador encoge y el denominador crece, de modo que $\alpha$ se va hacia 0. Aunque $\delta$ valga 1, el producto $\lambda$ encoge con él. Por eso un elemento aplanado no se queda en el frente.

La velocidad es el otro lado. Cada elemento nuevo debe contrastarse con las demás aristas del frente por intersección. Escrito de forma ingenua, eso es $O(n)$ por elemento y $O(n^2)$ en total. Si una malla de fondo o un quadtree no localizan la prueba, un frente largo gasta el tiempo en intersecciones y no en elementos.

Cuándo elegir cada uno#

Los dos métodos se reparten el trabajo más de lo que compiten.

Delaunay devuelve una respuesta única y rápida una vez fijados los puntos. Sirve para rellenar una envolvente convexa, interpolar datos dispersos y limpiar la malla a posteriori. No respeta una frontera por sí solo. Ese paso va aparte.

AFT es más fuerte en fidelidad de frontera, en la primera malla que produce y en la anisotropía. Cuesta más de implementar, y las estructuras de datos deciden el tiempo de ejecución.

Los generadores de producción los mezclan. AFT elige el orden de colocación. El criterio de Delaunay elige las conexiones. Ese es el reparto frente-Delaunay.

Una arista de perfil cortada termina la inserción pura#

Que doce puntos den 17 triángulos vale cuando los puntos ya están dentro del dominio. Si la superficie de un perfil queda cortada por una diagonal, parte de esos 17 deja de representar la forma. O se retriangula esa banda como tubería, o la arista cortada se deja en el frente como segmento de referencia y el siguiente punto $C$ se elige contra ella.

Apilar triángulos con un $\lambda = \alpha\delta$ grande no devuelve una frontera cortada. Dónde añadir puntos para que una arista restringida regrese está en recuperar la frontera cortada con puntos de Steiner. Primero se deja bien el núcleo de inserción de unas 60 líneas. Solo cuando la geometría pierde de verdad una arista se pasa a ese procedimiento.

Dos filosofías para triangular una nube de puntos — Bowyer–Watson y frente de avance