Path tracing progresivo: de un rayo a la integral

Un rayo por píxel te da visibilidad: qué se ve desde dónde. La iluminación es otra cosa, porque el brillo de un punto no depende de una dirección sino de todas: es una integral sobre el hemisferio, y esa integral no tiene solución analítica en cuanto la escena tiene más de dos objetos. El path tracing la estima disparando caminos aleatorios y promediando, con la garantía de que el promedio converge al valor exacto y con la penalización de que converge despacio, con la raíz cuadrada del número de muestras. Toda la ingeniería de esta lección sale de esas dos frases.

La ecuación y su estimador

La radiancia que sale de un punto en una dirección es la que ese punto emite más la integral, sobre el hemisferio alrededor de su normal, de la radiancia que le llega por la BRDF y por el coseno del ángulo de incidencia:

L_o(x, ω_o) = L_e(x, ω_o) + ∫  f_r(x, ω_i, ω_o) · L_i(x, ω_i) · (n · ω_i) dω_i
                             Ω

Es recursiva —L_i es la L_o de otro punto— e imposible de resolver a mano. Monte Carlo la convierte en una suma: si eliges N direcciones al azar con una densidad de probabilidad p, la media de f_r · L_i · cosθ / p converge a la integral. Ese cociente por la densidad es lo que hace el estimador insesgado: cada muestra se compensa por lo probable que era elegirla.

La elección de p decide cuánto ruido tendrás. Muestrear el hemisferio uniformemente funciona, pero desperdicia muestras en direcciones rasantes que el coseno va a aplastar a casi cero. Muestrear proporcionalmente al coseno —p(ω) = cosθ / π— reparte las muestras donde pesan, y con una BRDF lambertiana f_r = ρ/π produce una cancelación preciosa:

f_r · cosθ / p  =  (ρ/π) · cosθ · (π / cosθ)  =  ρ

Todo el estimador de un rebote difuso se reduce a multiplicar por el albedo. Ni divisiones, ni cosenos, ni el número pi. Esa cancelación es la razón de que un path tracer básico quepa en cuarenta líneas, y también un aviso: en cuanto cambies la BRDF por algo especular, la cancelación desaparece y vuelven todos los términos.

El bucle de rebotes

La formulación es recursiva y WGSL no tiene recursión, así que hay que aplanarla. El truco es llevar un acumulador de atenuación: en lugar de multiplicar al volver de la llamada recursiva, multiplicas al bajar. La radiancia final es la suma, sobre todos los rebotes, de la emisión encontrada por la atenuación acumulada hasta ese punto.

const MAX_REBOTES : u32 = 8u;

fn caminar(rayoInicial: Rayo) -> vec3<f32> {
  var radiancia  = vec3<f32>(0.0);
  var atenuacion = vec3<f32>(1.0);
  var rayo = rayoInicial;

  for (var rebote = 0u; rebote < MAX_REBOTES; rebote = rebote + 1u) {
    let h = trazar(rayo, 0.0, 1e30);

    if (h.id == SIN_IMPACTO) {                       // el camino se escapa al cielo
      radiancia = radiancia + atenuacion * cielo(rayo.dir);
      break;
    }

    let mat = materiales[matDePrim[h.id]];
    radiancia  = radiancia  + atenuacion * mat.emision;
    atenuacion = atenuacion * mat.albedo;            // la cancelacion del coseno

    // Ruleta rusa a partir del tercer rebote.
    if (rebote >= 2u) {
      let p = clamp(max(atenuacion.r, max(atenuacion.g, atenuacion.b)), 0.05, 1.0);
      if (aleatorio() > p) { break; }
      atenuacion = atenuacion / p;
    }

    let punto = rayo.origen + h.t * rayo.dir;
    let n = normalize(select(h.normal, -h.normal, dot(h.normal, rayo.dir) > 0.0));
    rayo = Rayo(desplazarOrigen(punto, n), muestrearCoseno(n));
  }
  return radiancia;
}

El for con límite constante no es una concesión estética: le dice al compilador cuál es el tiempo de vida máximo de cada variable y le permite decidir el reparto de registros. El tMin de trazar vale cero porque el origen ya viene desplazado en unidades de precisión por el mismo mecanismo que evita la autointersección.

La ruleta rusa merece su párrafo porque casi todo el mundo la aplica y muy poca gente sabe justificar por qué no ensucia el resultado. Sea X la contribución que le quedaba a este camino. En lugar de calcularla, lánzala a cara o cruz con probabilidad p de continuar, y si continúa multiplica el resultado por 1/p. El valor esperado del nuevo estimador es p · (X/p) + (1−p) · 0 = X: exactamente el mismo. Terminar caminos así no pierde energía, solo la reparte en menos muestras más brillantes, lo que sube un poco la varianza a cambio de bajar mucho el coste. Usar la componente máxima de la atenuación como p es lo estándar: los caminos que ya se han oscurecido mueren pronto, los que siguen aportando sobreviven. El clamp inferior a 0.05 evita que una superficie casi negra produzca un factor 1/p de veinte mil y una mota blanca permanente.

El límite duro de MAX_REBOTES sí introduce sesgo: la energía de los caminos más largos simplemente se pierde. Cuánta, se calcula. Con un albedo de 0,8, la atenuación tras ocho rebotes vale 0,8⁸ = 0,168, así que estás recortando una cola que como mucho aporta un diecisiete por ciento, y en la práctica bastante menos porque esa cola ya la estaba diezmando la ruleta. En una escena blanca de albedo cercano a uno, en cambio, el mismo límite produce un interior visiblemente más oscuro de lo que debería. Las dos técnicas conviven: ruleta rusa para terminar sin sesgo, límite duro como red de seguridad contra bucles infinitos entre dos espejos.

El muestreo del hemisferio con coseno se hace por el método de Malley —un punto uniforme del disco unidad, elevado al hemisferio— y necesita una base ortonormal alrededor de la normal. La construcción sin ramas de Duff y compañía la resuelve en seis líneas y sin ningún caso especial salvo el signo:

fn muestrearCoseno(n: vec3<f32>) -> vec3<f32> {
  let r1 = aleatorio();
  let r2 = aleatorio();

  let r   = sqrt(r1);                        // radio en el disco
  let phi = 6.2831853071795864 * r2;
  let x = r * cos(phi);
  let y = r * sin(phi);
  let z = sqrt(max(0.0, 1.0 - r1));          // altura: proyeccion al hemisferio

  // Base ortonormal sin ramas alrededor de n (n debe estar normalizada).
  let s = select(-1.0, 1.0, n.z >= 0.0);
  let a = -1.0 / (s + n.z);
  let b = n.x * n.y * a;
  let t1 = vec3<f32>(1.0 + s * n.x * n.x * a, s * b, -s * n.x);
  let t2 = vec3<f32>(b, s + n.y * n.y * a, -n.y);

  return normalize(t1 * x + t2 * y + n * z);
}

El s que salta de +1 a −1 según el signo de n.z es lo que evita la singularidad de a cuando la normal apunta hacia −z. Es el mismo motivo por el que la construcción ingenua con un producto vectorial contra un eje fijo se rompe: falla justo cuando la normal es paralela a ese eje.

Aleatoriedad en el shader

Cada invocación necesita su propio flujo de números, distinto del de sus vecinas y distinto del del fotograma anterior. Sin estado global compartido y sin poder usar un generador con siembra centralizada, la solución es un hash entero con estado local.

var<private> semilla : u32;

// PCG de 32 bits: un paso de LCG y una permutacion de salida (xorshift-multiply-xorshift).
fn pcgHash(v: u32) -> u32 {
  let estado  = v * 747796405u + 2891336453u;
  let palabra = ((estado >> ((estado >> 28u) + 4u)) ^ estado) * 277803737u;
  return (palabra >> 22u) ^ palabra;
}

fn sembrar(indicePixel: u32, indiceFrame: u32) {
  semilla = pcgHash(pcgHash(indicePixel) + indiceFrame);
}

fn siguiente() -> u32 {
  semilla = semilla * 747796405u + 2891336453u;
  let palabra = ((semilla >> ((semilla >> 28u) + 4u)) ^ semilla) * 277803737u;
  return (palabra >> 22u) ^ palabra;
}

fn aleatorio() -> f32 {
  return f32(siguiente()) * 2.3283064365386963e-10;   // 1 / 2^32
}

La siembra en dos niveles no es paranoia. Considera qué pasa si sirves los números directamente desde un LCG sembrado con el índice del píxel. La primera salida de un LCG es una función afín de la semilla: x₁ = a·x₀ + c módulo dos elevado a treinta y dos. Píxeles adyacentes tienen semillas que se diferencian en uno, así que sus primeras salidas se diferencian en exactamente a, siempre. Convertido a coma flotante, eso es una rampa perfectamente regular que recorre la imagen: bandas diagonales nítidas, no ruido. Y como el patrón está atado a la posición del píxel, no se promedia al acumular; se queda ahí. Hay un segundo fallo que agrava el primero: en un LCG de módulo potencia de dos, el bit k tiene periodo 2^(k+1), de modo que los bits bajos son casi constantes. Si sacas números tomando los bits bajos, no tienes un generador, tienes un contador.

La permutación de salida de PCG —desplazamiento dependiente de los propios datos, multiplicación por un primo grande y otro desplazamiento con xor— existe exactamente para destruir esa estructura afín, y por eso el mismo LCG se vuelve utilizable. El hash exterior de la siembra hace lo propio con la combinación de píxel y fotograma.

Una recomendación negativa que ahorra horas: el clásico fract(sin(dot(uv, vec2(12.9898, 78.233))) * 43758.5453) que arrastran mil shaders de la web depende de la precisión de sin con argumentos enormes, que es lo peor especificado de toda la aritmética de GPU. Funciona en tu portátil y produce bandas visibles o repeticiones en otro adaptador. Un hash entero cuesta lo mismo y es determinista en todas partes.

Acumular, y el coste real

El acumulador es un buffer de un vec4<f32> por píxel. Los tres primeros canales guardan la media de las muestras, no la suma, y el cuarto el número de muestras. Que sea la media importa: sumar diez mil muestras de magnitud cercana a uno lleva el acumulador a 10⁴, donde el escalón de un f32 vale 2⁻¹⁰ ≈ 0,00098, así que cualquier muestra por debajo de la mitad de eso deja de sumar y las zonas oscuras se congelan. Con la media incremental el acumulador se queda siempre en la magnitud de la imagen y ese problema no aparece.

@group(0) @binding(0) var<uniform> cam : Uniformes;
@group(0) @binding(1) var<storage, read_write> acumulador : array<vec4<f32>>;
@group(0) @binding(2) var salida : texture_storage_2d<rgba32float, write>;

@compute @workgroup_size(8, 8)
fn principal(@builtin(global_invocation_id) gid: vec3<u32>) {
  let dims = textureDimensions(salida);
  if (gid.x >= dims.x || gid.y >= dims.y) { return; }

  let indice = gid.y * dims.x + gid.x;
  sembrar(indice, cam.frame);

  let jitter = vec2<f32>(aleatorio(), aleatorio());
  let color  = caminar(generarRayo(gid.xy, dims, jitter));

  // Media incremental. WebGPU inicializa los buffers a cero, asi que en la
  // primera muestra n vale 1 y mix devuelve el color tal cual.
  let n      = f32(cam.muestras) + 1.0;
  let previo = acumulador[indice].rgb;
  let media  = mix(previo, color, 1.0 / n);

  acumulador[indice] = vec4<f32>(media, n);
  textureStore(salida, gid.xy, vec4<f32>(media, 1.0));
}

El reinicio se controla desde JavaScript con dos contadores que no son el mismo, y confundirlos es un fallo sutil y difícil de ver:

let muestras = 0;              // se reinicia al mover la camara
let numeroDeFrame = 0;         // no se reinicia nunca: es la semilla temporal
const vistaPrevia = new Float32Array(16);

function fotograma() {
  const vista = camara.matrizVista();
  let igual = true;
  for (let i = 0; i < 16; i++) {
    if (vista[i] !== vistaPrevia[i]) { igual = false; break; }
  }
  if (!igual) { muestras = 0; vistaPrevia.set(vista); }

  u[11] = numeroDeFrame++;
  u[15] = muestras;
  device.queue.writeBuffer(bufCamara, 0, datos);
  // ...dispatch de ceil(ancho/8) por ceil(alto/8) workgroups...
  muestras++;
}

Si reinicias también la semilla temporal, tras cada movimiento de cámara vuelves a disparar exactamente la misma secuencia de números aleatorios y las primeras muestras se repiten: el ruido deja de promediar y se congela un patrón. La comparación de la matriz es exacta, sin tolerancia, y hay que reiniciar por todo lo que invalide lo acumulado: la vista, la proyección, el tamaño del lienzo, un material, una luz.

Sobre el formato: rgba32float como storage texture con acceso write está en el núcleo de WebGPU y no necesita ninguna feature. Lo que sí la necesita es filtrarla con un sampler en la pasada de presentación —para escalar o desenfocar—: eso exige float32-filterable, que hay que solicitar en requestDevice y que no todos los adaptadores ofrecen. Si presentas a resolución nativa con textureLoad, no hay sampler y no hay problema. Vigila también el recuento de bindings: nodos, orden de primitivas, triángulos, acumulador, materiales e índices de material ya son seis storage buffers, y maxStorageBuffersPerShaderStage son ocho por defecto.

Ahora los números, que es donde esto se pone honesto. A 1920 por 1080 el acumulador ocupa 2 073 600 · 16 = 31,6 MiB, holgado frente a los 128 MiB de maxStorageBufferBindingSize. Cada muestra dispara un rayo primario más los rebotes: con cuatro rebotes efectivos son unos diez millones de rayos por muestra. Un recorrido de BVH bien escrito en compute, sin unidades de intersección en hardware, mueve del orden de cien millones de rayos por segundo en una GPU de escritorio de gama media si son coherentes, y bastante menos cuando la mayoría son rebotes difusos. Eso sitúa el resultado en torno a diez muestras por segundo a 1080p, y en una integrada o en un móvil hay que dividir por cinco o por diez.

Y sobre eso actúa la ley de la raíz cuadrada. El error estándar de un estimador de Monte Carlo es σ/√N, así que para reducir el ruido a la mitad hacen falta cuatro veces más muestras, y para dividirlo por diez, cien veces más. La progresión en la práctica: 16 muestras es grano grueso, 256 es cuatro veces menos ruidoso, 4096 otras cuatro veces menos. Si una escena de interior necesita mil muestras para quedar limpia y consigues diez por segundo, son cien segundos de imagen quieta.

Traducido a tiempo real: el presupuesto de un fotograma a sesenta por segundo son 16,6 milisegundos, y a diez muestras por segundo llevas 0,17 muestras por fotograma. Ni una muestra por píxel entera. Esto no da sesenta fotogramas por segundo a 1080p y no hay forma de escribirlo mejor para que los dé. Lo que sí funciona es bajar la resolución —a 960 por 540 tienes la cuarta parte de píxeles y por tanto cuatro veces las muestras— y aceptar la naturaleza del sistema: una imagen que se limpia sola mientras el usuario no toca la cámara. Con ese contrato, esto es perfectamente útil para configuradores de producto, previsualización de iluminación, horneado de mapas de luz y render de imagen fija en el navegador. Sin unidades de intersección en hardware, que WebGPU no expone, y sin un filtro de ruido, que hay que escribir aparte y cuesta varios milisegundos por fotograma, esa es la frontera honesta.