Cuando una bala atraviesa la pared — Detección continua de colisiones por bisección con funciones de inclusión

En una simulación que corre a 60fps, una bala veloz simplemente atraviesa la pared. En un fotograma está a la izquierda de la pared; en el siguiente, a la derecha. Nunca se solapa en ninguno de los dos instantes, así que el detector de colisiones informa que no pasó nada. Este artículo trata de inspeccionar ese "entremedio": la detección continua de colisiones (CCD, la técnica que halla el primer instante de contacto sobre toda la trayectoria lineal de dos cuerpos). Convertimos la colisión en una raíz polinómica mediante la condición de coplanaridad y luego implementamos una bisección con funciones de inclusión que jamás pierde una colisión, ni siquiera en coma flotante.

El tema toca todo código donde los cuerpos se mueven rápido entre fotogramas: interacción fluido-estructura (FSI), métodos de partículas y simulación de contacto basada en mallas.

La bala que se coló entre dos fotogramas#

La detección estática de colisiones solo ve una instantánea de un único momento. Pregunta si dos triángulos se solapan ahora mismo. Para cuerpos lentos basta. El problema es la velocidad.

Si un cuerpo se desplaza más que su propio grosor durante un paso de tiempo $\Delta t$, se cuela entre las instantáneas. Este fenómeno se llama túnel (tunneling). Y lo peor es que no hay recuperación. Una vez perdida la colisión, los dos cuerpos quedan interpenetrados, y en el siguiente paso ni siquiera queda claro qué cara está de qué lado.

Al subir point travel crece la distancia que el punto recorre en un paso. La barra gris es la pared.

Al subir point travel, la prueba discreta que solo mira los extremos cae pronto a MISS, porque ambos extremos quedan lejos de la pared. La prueba continua, en cambio, marca con un anillo rojo el instante exacto del primer contacto $t^\*$.

Por qué fuga cuando solo se revisan los extremos#

Revisar solo dos extremos significa muestrear una trayectoria continua en dos puntos. Lo que ocurrió entre las muestras es invisible.

La solución parte de un cambio de perspectiva. En vez de "¿se solapan ahora?", pregunta "¿para qué $t \in [0,1]$ se solapan?". Si cada vértice se mueve en línea recta durante un paso, su posición es una función lineal del tiempo.

Aquí $\mathbf{x}_i^{0}$ es la posición al inicio del paso y $\mathbf{x}_i^{1}$ la del final. Ahora la colisión se vuelve un problema de búsqueda de raíces en el tiempo.

Coplanaridad — convertir la colisión en una raíz polinómica#

La primera colisión entre dos triángulos en movimiento solo cae en dos casos. Un vértice toca una cara (vertex–face), o una arista toca otra arista (edge–edge). En ambos casos la observación geométrica clave es que cuatro puntos se vuelven coplanares en el instante de contacto.

La condición de coplanaridad para un vértice $\mathbf{p}$ y un triángulo $(\mathbf{a},\mathbf{b},\mathbf{c})$ es que el volumen (producto triple escalar) se anule.

Como cada $\mathbf{x}(t)$ es lineal en $t$, $f(t)$ es un polinomio cúbico en $t$. Es el enfoque univariado (univariate). Es rápido, pero tiene una trampa. Coplanar no es lo mismo que colisionando. Tras hallar una raíz $t^\*$ hay que comprobar aparte si el punto realmente cae dentro del triángulo. El vuelco de interior y exterior al caer un triángulo sobre una malla está contado en la paridad que voxeliza un STL. Perder una celda deja el mismo síntoma que una trayectoria que atravesó la pared.

Una vía más robusta es la forma multivariada (multivariate). Escribe la colisión directamente como el cero de un vector de holgura.

Hay colisión si existe un $(t,u,v)$ con $\mathbf{F}=\mathbf{0}$, $u,v \ge 0$ y $u+v \le 1$. Aquí $u,v$ son coordenadas baricéntricas sobre la cara. Esta forma no tiene casos de frontera y es puramente algebraica.

Las colisiones que la búsqueda numérica de raíces pierde#

Resolver la cúbica en coma flotante es lo más rápido, así que la mayoría de los simuladores hacen justo eso. Sin embargo, las pruebas de referencia a gran escala desde 2012 llegaron a una conclusión sorprendente: la mayoría de los solvers numéricos univariados producen falsos negativos (colisiones perdidas).

Las raíces perdidas salen de dos fallos distintos. Las raíces de una cúbica suelen ser irracionales, así que la coma flotante no las guarda exactas. Por otro lado, cuando dos cuerpos se deslizan en paralelo dentro del mismo plano, $f(t)$ es idénticamente cero. Un umbral numérico pensado para tirar ese caso de raíces infinitas se lleva también impactos reales.

Un falso positivo (informar una colisión que no existe) solo cuesta algo de trabajo extra y es recuperable. Pero un falso negativo lleva directo al túnel. Si el código de respuesta a colisiones evita la penetración con una búsqueda en línea, una sola pérdida hunde toda la simulación. Por eso la propiedad que CCD realmente necesita no es la velocidad, sino la conservatividad (conservativeness). Mejor informar un fantasma que perder una colisión.

Funciones de inclusión y bisección — acotar la raíz de forma conservadora#

La forma más antigua, y aún la más limpia, de garantizar la conservatividad es la bisección con funciones de inclusión de Snyder (1992).

La herramienta clave es la función de inclusión (inclusion function). Dada una caja de parámetros $B$, se calcula una caja envolvente alineada con los ejes $\Box\mathbf{F}(B)$ que encierra todos los valores que $\mathbf{F}$ puede tomar sobre $B$. Entonces se cumple un hecho simple: si $\mathbf{0} \notin \Box\mathbf{F}(B)$, no hay ninguna raíz dentro de $B$.

Por suerte, $\mathbf{F}$ es lineal en cada parámetro (multilineal). Una función multilineal siempre alcanza sus extremos en los vértices de una caja. Así que la caja de inclusión más ajustada es la envolvente de $\mathbf{F}$ evaluada en los cuatro (u ocho) vértices. No hace falta aritmética de intervalos por separado.

El algoritmo es divide y vencerás.

1. Inicializa la cola de candidatos con todo el dominio de parámetros $[0,1]^d$.
2. Extrae una caja $B$ de la cola y evalúa $\Box\mathbf{F}(B)$.
3. Si $\mathbf{0} \notin \Box\mathbf{F}(B)$, descártala (sin raíz garantizado).
4. Si el ancho de la caja es menor que $\delta$, acéptala como raíz.
5. Si no, parte por la mitad la dimensión más ancha y encola ambas mitades.
6. Para obtener el primer instante de contacto, ordena la cola por $t$.

En el cuadrado $(t,s)$, la celda cuyo recinto de inclusión no contiene el origen se descarta. Solo queda como raíz la celda más estrecha que $\delta$.

Python — primer contacto de un punto y un segmento en movimiento#

Implementamos el primer instante de contacto de un punto y un segmento en movimiento en 2D mediante bisección con funciones de inclusión. La función de holgura es $\mathbf{F}(t,s) = \mathbf{p}(t) - [\mathbf{a}(t) + s(\mathbf{b}(t)-\mathbf{a}(t))]$, y la raíz se busca en $(t,s)\in[0,1]^2$.

import heapq
import numpy as np
 
def gap_F(t, s, query):
    """F(t,s) = P(t) - [A(t) + s (B(t)-A(t))], trayectorias rectas (t in [0,1])."""
    P0, P1, A0, A1, B0, B1 = query
    P = (1 - t) * P0 + t * P1
    A = (1 - t) * A0 + t * A1
    B = (1 - t) * B0 + t * B1
    return P - (A + s * (B - A))
 
def inclusion_aabb(box, query):
    """Cota AABB más ajustada de F: envolvente de los cuatro vertices (exacta, multilineal)."""
    t0, t1, s0, s1 = box
    corners = [gap_F(t0, s0, query), gap_F(t0, s1, query),
               gap_F(t1, s0, query), gap_F(t1, s1, query)]
    lo = np.min(corners, axis=0)
    hi = np.max(corners, axis=0)
    return lo, hi
 
def origin_inside(lo, hi):
    return bool(np.all(lo <= 0) and np.all(hi >= 0))
 
def inclusion_bisection_toi(query, delta=1e-4, max_checks=200000):
    """Primer instante de contacto en [0,1], o None. Bisecccion de inclusion Snyder/Redon."""
    # cola de prioridad ordenada por t0 -> la primera caja convergida es el primer contacto
    pq = [(0.0, 0.0, 1.0, 0.0, 1.0)]   # (t0, t0, t1, s0, s1)
    checks = 0
    while pq and checks < max_checks:
        _, t0, t1, s0, s1 = heapq.heappop(pq)
        checks += 1
        lo, hi = inclusion_aabb((t0, t1, s0, s1), query)
        if not origin_inside(lo, hi):
            continue                              # sin colision en esta caja
        if max(t1 - t0, s1 - s0) < delta:
            return t0, checks                     # convergido: primer instante de contacto
        if (t1 - t0) >= (s1 - s0):                # partir la dimension mas ancha
            tm = 0.5 * (t0 + t1)
            heapq.heappush(pq, (t0, t0, tm, s0, s1))
            heapq.heappush(pq, (tm, tm, t1, s0, s1))
        else:
            sm = 0.5 * (s0 + s1)
            heapq.heappush(pq, (t0, t0, t1, s0, sm))
            heapq.heappush(pq, (t0, t0, t1, sm, s1))
    return None, checks
 
# Una bala que cruza una pared en un fotograma. Ambos extremos quedan lejos de la pared.
query = (np.array([0.5, 3.5]), np.array([9.5, 3.5]),   # punto P0 -> P1
         np.array([5.0, 1.0]), np.array([5.0, 1.0]),   # extremo de pared A (estatico)
         np.array([5.0, 6.0]), np.array([5.0, 6.0]))   # extremo de pared B (estatico)
 
toi, checks = inclusion_bisection_toi(query, delta=1e-4)
print(f"time-of-impact t* = {toi:.5f}   (cajas revisadas: {checks})")

Salida:

time-of-impact t* = 0.50000   (cajas revisadas: 53)

step avanza el corte un nivel, y $\delta$ es el ancho en el que una celda roja puede detenerse.

Las celdas rojas son la región de la raíz que bajó de $\delta$. Un $\delta$ mayor se detiene antes de las 53 revisiones y las celdas crecen. El $t^{*}=0.50000$ impreso es el borde izquierdo $t_0$ de esa caja, en cinco cifras.

Coincide exactamente con la solución analítica. El punto llega a $x=5$ en $t=0.5$, y allí $y=3.5$ cae dentro del rango del segmento $[1,6]$. Justo la colisión que la prueba de solo extremos respondió como MISS quedó acotada en 53 evaluaciones de caja. Conservadora y, aun así, barata.

El impacto es una división. $x$ va de 0.5 a 9.5, una distancia de 9, y la pared en $x=5$ queda a 4.5, así que $t=(5-0.5)/9=0.5$, el 0.50000 impreso. $y$ se queda en 3.5, y el segmento $[1,6]$ deja 2.5 de holgura en cada extremo. Bajar un eje desde ancho 1 hasta $\delta=10^{-4}$ pide 14 cortes, porque $2^{13}=8192$ y $2^{14}=16384$. La corrida revisó 53 cajas. Un barrido completo $14\times 14=196$ es mayor. La inclusión en los vértices tira las cajas vacías. El valor devuelto es el borde izquierdo $t_0$ de la caja aceptada. Imprimir 0.50000 significa que ese borde coincide con 0.5 en cinco cifras.

Intercambiar precisión por velocidad con un solo δ#

La virtud más práctica de la bisección con inclusión es que tiene un único mando de usuario, $\delta$.

Subir $\delta$ parte menos cajas y se detiene antes, a costa de una región de raíz mayor y más falsos positivos. Bajarlo afina el resultado, pero aumenta el número de revisiones. En simulación en tiempo real se fija además un tope de revisiones $m_I$ para acotar el peor tiempo de ejecución. Al alcanzar el tope, se devuelve de forma conservadora el último intervalo de colisión registrado: sigue sin haber pérdidas.

Añade encima la separación mínima (minimum separation) y podrás evitar también la penetración por tolerancias de fabricación o por redondeo. Un pequeño truco —cambiar la distancia de euclídea a Chebyshev (L∞)— simplifica el problema. En vez de $\mathbf{F}=\mathbf{0}$ se resuelve $\lVert\mathbf{F}\rVert_\infty \le d$, lo cual equivale a comprobar si un cubo de lado $2d$ se solapa con la caja de inclusión en lugar del origen. En código es una sola línea que agranda la caja en $d$.

Antes de llevarlo a código#

• Redondea siempre en el sentido conservador. Suma una cota de error hacia adelante a las evaluaciones de $\mathbf{F}$ en los vértices para agrandar un poco la caja de inclusión. Esa línea evita los falsos negativos.
• Cuidado con el atractivo de la cúbica univariada. Es rápida, pero se atasca en raíces infinitas bajo movimiento coplanar paralelo. Cuando necesites conservatividad, usa multivariado más bisección.
• Pon antes una fase amplia (broad phase). Filtra pares candidatos con un hash espacial o un árbol AABB antes de llamar a CCD. CCD es la costosa última puerta de la fase estrecha.
• Normaliza $\delta$ y $m_I$ a la escala de la escena. Átalos al tamaño del dominio y no a una longitud absoluta, para que sigan funcionando al cambiar velocidad y escala.

Subir δ hincha la celda y la colisión sigue dentro#

$\delta=10^{-4}$ pide 14 cortes en un eje. Con ese ancho la salida es $t^{*}=0.50000$ tras 53 revisiones. Un $\delta$ mayor agranda las celdas rojas y suma falsos positivos, pero no empuja el origen fuera del recinto de inclusión. Las pérdidas vienen de filtrar una cúbica coplanar con un umbral numérico, no del ancho.

Dos extremos llaman a esta bala MISS. En $t=0.5$ está en $x=5$, $y=3.5$, dentro de $[1,6]$, y las dos instantáneas quedan lejos de la pared. Cuando el fluido vive en la malla y el sólido cruza entre celdas, la frontera inmersa de Peskin no clava el contacto en un nodo. Si esa frontera cruza una celda entre fotogramas, la prueba de extremos vuelve a informar MISS.