Construir la BVH: partición, SAH y aplanado

Sin estructura de aceleración, cada rayo pregunta a todas las primitivas de la escena: coste lineal, y con dos millones de rayos por fotograma eso es inviable a partir de unos pocos miles de triángulos. Una jerarquía de volúmenes envolventes convierte esa pregunta en un descenso por un árbol y baja el coste a logarítmico. Construirla es un trabajo de CPU que se hace una vez, pero cada decisión que tomas ahí —dónde cortar, cuándo parar, en qué orden escribir los nodos en memoria— se paga o se cobra millones de veces por fotograma durante el recorrido.

Lineal frente a logarítmico, con números

Un modelo de un millón de triángulos a 1920 por 1080 son 2 073 600 rayos primarios. Por fuerza bruta, cada uno pregunta a un millón de triángulos: 2,07 billones de tests de intersección para una sola imagen sin un solo rebote. Con una jerarquía, cada rayo desciende por el árbol y visita del orden de treinta a sesenta cajas y prueba un puñado de triángulos en las hojas que atraviesa. Redondeando a cuarenta y cinco tests por rayo, son 93 millones de tests: un factor de veintidós mil. Ese es el algoritmo entero.

La forma correcta de leer ese número no es “la jerarquía acelera veintidós mil veces”. Es que el coste por rayo pasa de crecer proporcionalmente al número de primitivas a crecer con su logaritmo. Duplicar la escena, de un millón a dos millones de triángulos, duplica el coste por fuerza bruta y añade un nivel al árbol: un test de caja más por rayo. Por eso el trazado de rayos es el único algoritmo de visibilidad cuyo coste apenas nota la complejidad geométrica, y por eso las escenas de cine se trazan y no se rasterizan.

El coste no es gratis. Hay que pagarlo en tres monedas. Memoria: el árbol ocupa entre 16 y 64 MB para un millón de triángulos, según el tamaño de las hojas. Tiempo de construcción: unos segundos en JavaScript para ese mismo millón, así que se construye en un worker o directamente fuera de línea y se envía el buffer ya aplanado. Y latencia de memoria durante el recorrido, que es la moneda que más duele y de la que se ocupa el recorrido en GPU.

Particionar primitivas, no espacio

Una BVH es un árbol binario donde cada nodo guarda una caja alineada a los ejes que contiene a todas las primitivas de su subárbol, y cada hoja guarda un rango de primitivas. La caja alineada a los ejes gana a la esfera y a la caja orientada porque su test de intersección con un rayo son seis restas, seis multiplicaciones y cuatro mínimos: menos operaciones que un solo triángulo, que es exactamente la condición para que descartar salga a cuenta.

La palabra que separa la BVH de sus competidoras es partición. Un kd-tree y un octree parten el espacio: eligen un plano o un punto y dividen el volumen en dos o en ocho. La consecuencia es que una primitiva que cruza el plano de corte tiene que estar referenciada en las dos mitades, así que el número de referencias no está acotado por el número de primitivas y una escena mal dispuesta puede hacerlo explotar. A cambio, las celdas nunca se solapan y el espacio vacío se descarta con exactitud.

Una BVH parte el conjunto de primitivas: cada triángulo va a una de las dos ramas y a una sola. El número de nodos está acotado por 2n − 1, la memoria es predecible antes de empezar, y la construcción no tiene que decidir nada sobre primitivas partidas. Lo que pierde es la exclusividad: las cajas de dos hermanos pueden solaparse, y un rayo que entra en la zona común está obligado a bajar por las dos ramas. Todo el trabajo de la heurística de construcción consiste en minimizar ese solapamiento.

Hay una tercera razón por la que la BVH ganó y no aparece en los libros de texto antiguos: se puede refitear. Si la geometría se deforma pero conserva su topología —un personaje animado, tela, una malla que respira—, basta con recorrer el árbol de abajo arriba recalculando las cajas y dejando la estructura intacta. La calidad se degrada poco a poco, pero cuesta una fracción de una reconstrucción. Un kd-tree no admite eso: sus planos de corte dejan de ser válidos en cuanto la geometría se mueve.

La construcción de arriba abajo tiene tres decisiones. El eje: el más largo de la caja de los centroides, no de la caja del nodo. La distinción importa: un nodo con un triángulo enorme y cien pequeños tiene una caja alargada en la dirección del grande, pero los centroides pueden estar todos apretados en otro eje, y cortar por el eje equivocado deja una rama vacía. El punto de corte: la mediana da un árbol perfectamente equilibrado, de profundidad exactamente log₂ n, y es lo que hay que hacer si necesitas construir rápido. El criterio de parada: dejar de partir cuando queda una primitiva es lo peor que puedes hacer, porque duplica el número de nodos para ahorrar un test de triángulo que cuesta lo mismo que el de caja.

La mediana tiene un fallo famoso, el de la tetera en el estadio: una escena con un objeto denso y unos pocos triángulos gigantes lejanos. La mediana reparte por número, así que mete la mitad de la tetera en la misma rama que el estadio, y esa rama tiene una caja del tamaño del estadio con casi nada dentro. Cualquier rayo que la toque baja por ella para nada. La heurística del área de superficie existe precisamente para no hacer eso.

La heurística del área de superficie

La pregunta que resuelve es cuánto va a costar, en promedio, un rayo que llegue a este nodo, en función de dónde cortes. El coste esperado de un nodo interno es el del propio test más el coste de bajar por cada hijo, ponderado por la probabilidad de tener que bajar:

C = Ct + P(izq) · N_izq · Ci + P(der) · N_der · Ci

Ct es el coste de recorrer un nodo, Ci el de intersectar una primitiva, y N el número de primitivas de cada lado. Toda la heurística depende de esas dos probabilidades, y ahí es donde entra la geometría.

La probabilidad es la razón de áreas. El resultado viene de la geometría integral: la medida de las rectas que cortan un cuerpo convexo es proporcional a su área de superficie. La fórmula de Cauchy lo dice de forma más concreta —el área proyectada media de un convexo sobre todas las direcciones vale exactamente su área partida por cuatro—, y de ahí, para un convexo B contenido en otro convexo A, la probabilidad condicional de que una recta uniforme que corta A corte también B es Área(B) / Área(A). Sustituyendo:

C = Ct + ( A_izq · N_izq + A_der · N_der ) · Ci / A

Ese A del denominador es común a todos los cortes candidatos del mismo nodo, así que para comparar basta con minimizar el numerador. Y el coste de no partir, de dejarlo como hoja, es simplemente N · Ci: si ningún corte lo mejora, la hoja gana y el criterio de parada deja de ser un número arbitrario.

Sobre los valores: PBRT usa Ct = 1/8 con Ci = 1, otros motores usan proporciones más cercanas a 1/2. La diferencia importa menos de lo que parece porque el mínimo de la función de coste es muy plano; lo que sí importa es que Ct sea claramente menor que Ci, que es lo que empuja al árbol a hacerse profundo y a poner pocas primitivas por hoja.

Evaluar la SAH exacta significa ordenar las primitivas por cada eje y barrer todos los cortes posibles: O(n log n) por nodo, con constantes horribles. El binning la aproxima: proyecta los centroides en un número fijo de cubetas —doce o dieciséis es lo habitual—, acumula por cubeta el conteo y la caja, y evalúa solo los K − 1 cortes entre cubetas con dos barridos. El coste por nodo baja a O(n + K) y la calidad del árbol resultante queda a un uno o dos por ciento del óptimo.

const BINS = 12;
const COSTE_RECORRIDO = 0.125;   // Ct, con Ci = 1
const MAX_PRIMS_HOJA  = 8;

interface Caja { min: number[]; max: number[] }

const cajaVacia = (): Caja => ({
  min: [Infinity, Infinity, Infinity],
  max: [-Infinity, -Infinity, -Infinity],
});

function unirCaja(a: Caja, b: Caja): void {
  for (let k = 0; k < 3; k++) {
    if (b.min[k] < a.min[k]) a.min[k] = b.min[k];
    if (b.max[k] > a.max[k]) a.max[k] = b.max[k];
  }
}

function area(c: Caja): number {
  const dx = c.max[0] - c.min[0];
  const dy = c.max[1] - c.min[1];
  const dz = c.max[2] - c.min[2];
  if (dx < 0 || dy < 0 || dz < 0) return 0;          // caja vacia
  return 2 * (dx * dy + dy * dz + dz * dx);
}

interface NodoTmp {
  caja: Caja;
  primero: number;
  cuenta: number;              // 0 => nodo interno
  izq: NodoTmp | null;
  der: NodoTmp | null;
}

function construir(
  cajas: Caja[], centroides: Float32Array, orden: Uint32Array,
  ini: number, fin: number,
): NodoTmp {
  const caja = cajaVacia();
  const cc   = cajaVacia();                          // caja de los centroides
  for (let i = ini; i < fin; i++) {
    const p = orden[i];
    unirCaja(caja, cajas[p]);
    for (let k = 0; k < 3; k++) {
      const c = centroides[p * 3 + k];
      if (c < cc.min[k]) cc.min[k] = c;
      if (c > cc.max[k]) cc.max[k] = c;
    }
  }

  const n = fin - ini;
  const hoja = (): NodoTmp => ({ caja, primero: ini, cuenta: n, izq: null, der: null });
  if (n <= 1) return hoja();

  const ext = [cc.max[0] - cc.min[0], cc.max[1] - cc.min[1], cc.max[2] - cc.min[2]];
  const eje = ext[0] > ext[1] ? (ext[0] > ext[2] ? 0 : 2) : (ext[1] > ext[2] ? 1 : 2);
  if (ext[eje] <= 0) return hoja();                  // todos los centroides coinciden

  const escala = BINS / ext[eje];
  const cuentaBin = new Int32Array(BINS);
  const cajaBin: Caja[] = Array.from({ length: BINS }, cajaVacia);
  const cubeta = (p: number) =>
    Math.min(BINS - 1, ((centroides[p * 3 + eje] - cc.min[eje]) * escala) | 0);

  for (let i = ini; i < fin; i++) {
    const p = orden[i];
    const b = cubeta(p);
    cuentaBin[b]++;
    unirCaja(cajaBin[b], cajas[p]);
  }

  // Dos barridos: area y conteo acumulados por la izquierda y por la derecha.
  const aIzq = new Float64Array(BINS - 1), nIzq = new Int32Array(BINS - 1);
  const aDer = new Float64Array(BINS - 1), nDer = new Int32Array(BINS - 1);

  let acu = cajaVacia(), cnt = 0;
  for (let b = 0; b < BINS - 1; b++) {
    unirCaja(acu, cajaBin[b]); cnt += cuentaBin[b];
    aIzq[b] = area(acu); nIzq[b] = cnt;
  }
  acu = cajaVacia(); cnt = 0;
  for (let b = BINS - 1; b > 0; b--) {
    unirCaja(acu, cajaBin[b]); cnt += cuentaBin[b];
    aDer[b - 1] = area(acu); nDer[b - 1] = cnt;
  }

  const invA = 1 / area(caja);
  let mejor = n <= MAX_PRIMS_HOJA ? n : Infinity;    // coste de dejarlo como hoja
  let corte = -1;
  for (let b = 0; b < BINS - 1; b++) {
    if (nIzq[b] === 0 || nDer[b] === 0) continue;
    const coste = COSTE_RECORRIDO + invA * (aIzq[b] * nIzq[b] + aDer[b] * nDer[b]);
    if (coste < mejor) { mejor = coste; corte = b; }
  }
  if (corte < 0) return hoja();

  // Particion in situ del array de indices, estilo Hoare.
  let i = ini, j = fin - 1;
  while (i <= j) {
    const p = orden[i];
    if (cubeta(p) <= corte) { i++; }
    else { orden[i] = orden[j]; orden[j] = p; j--; }
  }
  if (i === ini || i === fin) return hoja();

  return {
    caja, primero: 0, cuenta: 0,
    izq: construir(cajas, centroides, orden, ini, i),
    der: construir(cajas, centroides, orden, i, fin),
  };
}

El array orden es una permutación de índices de primitiva que la construcción reordena in situ. Al terminar, las primitivas de cualquier hoja ocupan un rango contiguo de ese array, que es lo que permite que una hoja se describa con dos números: dónde empieza y cuántas hay. Ese array también viaja a la GPU.

Aplanar a 32 bytes

Un árbol de punteros no sirve en GPU. Hay que convertirlo en un array plano, y el diseño del nodo es una de esas decisiones que parecen cosméticas y no lo son.

struct Nodo {
  min    : vec3<f32>,   // offset  0, 12 bytes
  indice : u32,         // offset 12: hoja -> primera primitiva | interno -> hijo derecho
  max    : vec3<f32>,   // offset 16, 12 bytes
  cuenta : u32,         // offset 28: 0 => nodo interno
};

@group(0) @binding(0) var<storage, read> nodos      : array<Nodo>;
@group(0) @binding(1) var<storage, read> ordenPrims : array<u32>;

Son exactamente 32 bytes y no lleva ni un byte de relleno, porque vec3<f32> en WGSL tiene alineación 16 y tamaño 12: cada vector deja hueco justo para el u32 que le sigue. Eso no es casualidad, es el motivo por el que este layout es universal. Un nodo de 32 bytes son dos accesos de 16 bytes, y con líneas de caché de 64 o 128 bytes caben dos o cuatro nodos por línea. Estirarlo a 36 bytes para meter un campo más rompe la alineación y hace que un nodo pueda quedar a caballo entre dos líneas: puedes duplicar los fallos de caché del recorrido por añadir cuatro bytes.

El truco del aplanado está en el orden en profundidad. Si escribes el árbol recorriéndolo en preorden, el hijo izquierdo de un nodo cae siempre en la posición inmediatamente siguiente. Eso ahorra cuatro bytes, pero lo que de verdad compra es localidad: el nodo que vas a visitar a continuación en el caso común está pegado al que acabas de leer, muchas veces en la misma línea de caché. Solo el hijo derecho necesita un índice explícito, y va en el mismo campo que las hojas usan para la primera primitiva porque un nodo nunca es las dos cosas.

const contar = (n: NodoTmp): number =>
  n.cuenta > 0 ? 1 : 1 + contar(n.izq!) + contar(n.der!);

function aplanar(raiz: NodoTmp): ArrayBuffer {
  const buf = new ArrayBuffer(contar(raiz) * 32);
  const f32 = new Float32Array(buf);
  const u32 = new Uint32Array(buf);
  let siguiente = 0;

  const escribir = (nodo: NodoTmp): number => {
    const yo = siguiente++;
    const o  = yo * 8;                                    // 8 palabras de 32 bits
    f32[o + 0] = nodo.caja.min[0];
    f32[o + 1] = nodo.caja.min[1];
    f32[o + 2] = nodo.caja.min[2];
    f32[o + 4] = nodo.caja.max[0];
    f32[o + 5] = nodo.caja.max[1];
    f32[o + 6] = nodo.caja.max[2];

    if (nodo.cuenta > 0) {
      u32[o + 3] = nodo.primero;
      u32[o + 7] = nodo.cuenta;                           // mayor que 0 => hoja
    } else {
      escribir(nodo.izq!);                                // cae en yo + 1: implicito
      u32[o + 3] = escribir(nodo.der!);                   // el derecho, explicito
      u32[o + 7] = 0;                                     // 0 => nodo interno
    }
    return yo;
  };

  escribir(raiz);
  return buf;
}

const bufNodos = device.createBuffer({
  size: nodosPlanos.byteLength,
  usage: GPUBufferUsage.STORAGE | GPUBufferUsage.COPY_DST,
  label: 'bvh-nodos',
});
device.queue.writeBuffer(bufNodos, 0, nodosPlanos);

Cuenta la memoria antes de dar por buena una escena. Con hojas de hasta ocho primitivas y SAH, un millón de triángulos produce del orden de medio millón de nodos: 16 MB, cómodo. En el peor caso, una hoja por primitiva, son 2n − 1 nodos, casi dos millones, 64 MB: la mitad de los 128 MiB de maxStorageBufferBindingSize. Los triángulos son el otro consumidor grande. Guardados como tres vec4<f32> son 48 bytes cada uno, y 134 217 728 / 48 = 2 796 202: con un solo binding de storage no pasas de unos 2,8 millones de triángulos, y ese, no la memoria de la GPU, es el techo real de una escena en WebGPU con el layout ingenuo. Bajar a índices de 32 bits sobre un buffer de vértices compartido multiplica ese techo por tres o por cuatro.