Generar rayos: del píxel a la dirección

La rasterización lleva treinta años respondiendo a una pregunta muy concreta: dado este triángulo, qué píxeles cubre. Es una pregunta cómoda porque se contesta proyectando y recorriendo un área, y porque los píxeles que toca están juntos en memoria. El trazado de rayos hace la pregunta contraria —dado este píxel, qué geometría veo— y esa inversión lo cambia todo: la estructura de datos, el patrón de acceso a memoria y el modelo de coste. Aquí montamos la primera mitad del problema, que es la barata: convertir dos coordenadas enteras de píxel en un origen y una dirección.

Dispersar frente a recolectar

Un rasterizador es un bucle sobre triángulos. Para cada uno proyecta sus tres vértices, calcula el rectángulo que los envuelve en pantalla y recorre los píxeles de dentro comprobando el signo de tres funciones de arista. Es un algoritmo de dispersión: un dato de entrada escribe en muchas salidas, y esas salidas están contiguas. El buffer de profundidad resuelve la visibilidad a posteriori, comparando y descartando. El coste es proporcional al número de triángulos más el número de fragmentos generados, y el acceso a memoria es casi perfecto porque un triángulo pequeño toca un puñado de píxeles vecinos que caben en el mismo tile de la caché de color.

Un trazador es un bucle sobre píxeles. Para cada uno construye un rayo y pregunta cuál es la primera superficie que corta. Es un algoritmo de recolección: cada salida lee de todas las entradas. Sin estructura de aceleración el coste es el producto de píxeles por primitivas, y el acceso a memoria es un desastre porque dos rayos vecinos pueden acabar leyendo triángulos que están en extremos opuestos del buffer.

Pon números. Un modelo de un millón de triángulos a 1920 por 1080 son 2 073 600 píxeles. La rasterización hace un millón de configuraciones de triángulo y del orden de diez millones de fragmentos. El trazado por fuerza bruta hace 2 073 600 por 1 000 000 = 2,07 billones de tests de intersección. Un test de triángulo cuesta unas veinte operaciones en coma flotante, así que son cuatro por diez elevado a trece operaciones para un solo fotograma: en una GPU de cinco teraflops, ocho segundos por imagen suponiendo eficiencia perfecta. Por eso la estructura de aceleración no es una optimización opcional sino la mitad del algoritmo.

Entonces, ¿por qué molestarse? Porque la rasterización solo sabe contestar la pregunta desde un punto de vista a la vez, y solo desde el que ha proyectado. El trazado la contesta desde cualquier punto y en cualquier dirección. Esa generalidad es exactamente lo que hace falta para sombras duras y blandas, reflejos de superficies curvas, refracciones y luz indirecta: todos son la misma consulta —qué ve este punto mirando hacia allá— hecha desde sitios que no son la cámara. Los trucos que la rasterización usa para fingir esos efectos (mapas de sombra, sondas de reflejo, oclusión ambiental en espacio de pantalla) son aproximaciones a una consulta que el trazado responde de forma exacta.

Y WebGPU no tiene nada de esto en hardware. No existe ninguna API de aceleración de rayos en la especificación: ni estructuras de aceleración gestionadas, ni pipelines de rayos, ni unidades de intersección expuestas. Todo lo que hagas aquí es cómputo puro sobre storage buffers, exactamente igual que una simulación de partículas.

Del píxel a la dirección

Un rayo son seis números: origen y dirección. El origen de un rayo primario de cámara pinhole es siempre el mismo punto, la posición de la cámara. La dirección es lo único que hay que calcular, y hay dos formas de hacerlo.

La primera es la geométrica, con la base ortonormal de la cámara. Necesitas derecha, arriba y adelante como vectores unitarios ortogonales, la tangente de la mitad del campo de visión vertical y la relación de aspecto. La construcción es directa: pasas la coordenada de píxel a un rango normalizado de menos uno a uno, la escalas por el semiancho del plano de imagen a distancia uno, y combinas.

struct Camara {
  origen      : vec3<f32>,
  tanMitadFov : f32,        // tan(fovY / 2)
  derecha     : vec3<f32>,
  aspecto     : f32,        // ancho / alto
  arriba      : vec3<f32>,
  frame       : u32,
  adelante    : vec3<f32>,
  muestras    : u32,
};

struct Rayo {
  origen : vec3<f32>,
  dir    : vec3<f32>,
};

@group(0) @binding(0) var<uniform> cam : Camara;

fn generarRayo(px: vec2<u32>, dims: vec2<u32>, jitter: vec2<f32>) -> Rayo {
  // uv en [0,1): el jitter reparte la muestra dentro del pixel.
  let uv  = (vec2<f32>(px) + jitter) / vec2<f32>(dims);
  // ndc en [-1,1). La fila 0 de la textura es la de arriba: y se invierte.
  let ndc = vec2<f32>(uv.x * 2.0 - 1.0, 1.0 - uv.y * 2.0);

  let x = ndc.x * cam.aspecto * cam.tanMitadFov;
  let y = ndc.y * cam.tanMitadFov;

  return Rayo(cam.origen, normalize(cam.derecha * x + cam.arriba * y + cam.adelante));
}

Esa struct Camara ocupa exactamente 64 bytes y no lleva relleno explícito porque los vec3<f32> de WGSL tienen alineación 16 y tamaño 12: cada uno deja hueco justo para el escalar que le sigue. El lado de JavaScript escribe las dos vistas sobre el mismo ArrayBuffer:

const datos = new ArrayBuffer(64);
const f = new Float32Array(datos);
const u = new Uint32Array(datos);

f.set(origen, 0);    f[3]  = Math.tan(fovY * 0.5);
f.set(derecha, 4);   f[7]  = ancho / alto;
f.set(arriba, 8);    u[11] = frame;
f.set(adelante, 12); u[15] = muestras;

device.queue.writeBuffer(bufCamara, 0, datos);

La segunda forma es con la inversa de la matriz de proyección-vista, y tiene la ventaja de que reutiliza exactamente la misma matriz que usa tu pasada de rasterización, así que los dos caminos coinciden píxel a píxel. Aquí importa un detalle de WebGPU: su volumen de recorte va de cero a uno en profundidad, no de menos uno a uno como el de WebGL. El plano cercano es z = 0 y el lejano z = 1.

@group(0) @binding(1) var<uniform> invViewProj : mat4x4<f32>;

fn generarRayoMatriz(uv: vec2<f32>) -> Rayo {
  let ndc = vec2<f32>(uv.x * 2.0 - 1.0, 1.0 - uv.y * 2.0);

  let hCerca = invViewProj * vec4<f32>(ndc, 0.0, 1.0);   // z = 0 -> plano cercano
  let hLejos = invViewProj * vec4<f32>(ndc, 1.0, 1.0);   // z = 1 -> plano lejano

  let pCerca = hCerca.xyz / hCerca.w;
  let pLejos = hLejos.xyz / hLejos.w;

  return Rayo(pCerca, normalize(pLejos - pCerca));
}

Con una proyección de plano lejano infinito el segundo punto degenera —w tiende a cero— y hay que desproyectar dos profundidades finitas, por ejemplo 0.0 y 0.5, en lugar de los dos planos. Con una proyección normal funciona tal cual.

Normalizar la dirección no es opcional aunque parezca que sí. Todo el código que viene después interpreta t como una distancia en unidades del mundo: el intervalo de búsqueda, el desplazamiento del origen para evitar la autointersección, la atenuación por distancia. Si la dirección no es unitaria, t es una distancia en unidades de “longitudes de dirección” y cada rayo usa una escala distinta.

El dispatch bidimensional

Una imagen es una rejilla, así que el dispatch es bidimensional. @workgroup_size(8, 8) da 64 invocaciones por grupo, cómodamente por debajo del límite de 256 de maxComputeInvocationsPerWorkgroup, y con una propiedad que importa mucho más que el número: las 64 invocaciones cubren un bloque cuadrado de pantalla. Los 64 rayos de ese bloque salen del mismo punto con direcciones casi paralelas, así que atraviesan casi los mismos nodos del árbol y leen casi la misma memoria. Un @workgroup_size(64) unidimensional cubriría una tira de 64 píxeles de ancho por uno de alto, con mucha menos coherencia vertical y peor comportamiento de caché.

Ojo con el matiz de cómo se reparten esas invocaciones entre los grupos SIMD del hardware. El índice lineal dentro del workgroup es y * 8 + x, así que en una máquina con grupos de 32 carriles el primero cubre las filas 0 a 3 y el segundo las filas 4 a 7: dos bloques de 8 por 4, que siguen siendo compactos en 2D. Ese es todo el argumento a favor del 8 por 8, y es suficiente.

@group(0) @binding(2) var salida : texture_storage_2d<rgba32float, write>;

@compute @workgroup_size(8, 8)
fn generar(@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;
  let jitter = jitterDe(indice, cam.frame);
  let rayo   = generarRayo(gid.xy, dims, jitter);

  // Placeholder verificable: un degradado vertical de cielo.
  let a = 0.5 * (rayo.dir.y + 1.0);
  let color = mix(vec3<f32>(1.0, 1.0, 1.0), vec3<f32>(0.5, 0.7, 1.0), a);

  textureStore(salida, gid.xy, vec4<f32>(color, 1.0));
}

El número de workgroups sale de redondear hacia arriba, y ahí aparece el problema que la comprobación de límites resuelve:

const pass = enc.beginComputePass({ label: 'rayos-primarios' });
pass.setPipeline(pipeline);
pass.setBindGroup(0, bindGroup);
pass.dispatchWorkgroups(Math.ceil(ancho / 8), Math.ceil(alto / 8));
pass.end();

A 1920 por 1080 la división es exacta —240 por 135 grupos— y no sobra nada. A 1366 por 768 no: 1366 / 8 = 170,75, así que lanzas 171 grupos y cubres 1368 columnas. Sobran dos columnas por 768 filas, 1536 invocaciones fantasma. Si escribes a una storage texture, WebGPU descarta esas escrituras y no pasa nada visible. Si escribes a un storage buffer con el índice lineal gid.y * ancho + gid.x, el desastre es silencioso: la invocación (1366, 5) calcula 5 * 1366 + 1366 = 8196, que es el índice del primer píxel de la fila 6. Un índice perfectamente válido, el píxel equivocado, y una carrera de escritura contra la invocación que de verdad le corresponde. El resultado es una columna de píxeles corruptos en el borde izquierdo que aparece y desaparece según la GPU. Escribe la comprobación siempre.

La alternativa a la storage texture es un storage buffer de acumulación, que es lo que necesitas en cuanto sumas muestras entre fotogramas:

@group(0) @binding(3) var<storage, read_write> acumulador : array<vec4<f32>>;

A 1920 por 1080 con vec4<f32> son 2 073 600 · 16 = 33 177 600 bytes, unos 31,6 MiB. Cabe de sobra en los 128 MiB de maxStorageBufferBindingSize y en los 256 MiB de maxBufferSize. A 3840 por 2160 son 126,6 MiB y ya rozas el límite de binding con un solo buffer: a partir de ahí toca dividirlo en varios o bajar a rgba16float.

El jitter que sale gratis

Sin jitter, acumular mil fotogramas de la misma cámara da exactamente la misma imagen mil veces. Todos los rayos salen del centro geométrico del píxel, todos cortan la misma geometría, y el borde de un triángulo sigue siendo el mismo escalón duro. Con jitter, cada fotograma dispara el rayo por un punto aleatorio distinto dentro del píxel, y la media de las muestras converge a la integral de la radiancia sobre el área del píxel. Eso es antialiasing por supermuestreo, y sale gratis porque la acumulación ya estaba ahí.

Hace falta un número aleatorio distinto por píxel y por fotograma. Para el jitter basta un hash entero de buena calidad:

// lowbias32: avalancha completa, tres multiplicaciones y tres desplazamientos.
fn hash(v: u32) -> u32 {
  var x = v;
  x = x ^ (x >> 16u);
  x = x * 0x7feb352du;
  x = x ^ (x >> 15u);
  x = x * 0x846ca68bu;
  x = x ^ (x >> 16u);
  return x;
}

fn aFloat(u: u32) -> f32 { return f32(u) * 2.3283064365386963e-10; }  // 1 / 2^32

fn jitterDe(indice: u32, frame: u32) -> vec2<f32> {
  let h = hash(indice ^ hash(frame));
  return vec2<f32>(aFloat(h), aFloat(hash(h)));
}

Fíjate en el sembrado: hash(indice ^ hash(frame)). El hash del fotograma descorrelaciona el tiempo, el xor con el índice descorrelaciona el espacio, y el hash exterior mezcla los dos. Sembrar con indice + frame a secas y esperar que el hash arregle el resto funciona con este hash concreto porque tiene avalancha completa, pero deja de funcionar en cuanto lo cambias por algo más barato, y el síntoma —franjas diagonales que se mueven— tarda en atribuirse a la causa. El generador con estado que necesita el path tracing es una versión más elaborada de esta misma idea.

Un jitter uniforme en el cuadrado del píxel equivale a filtrar con una caja, que es el filtro de reconstrucción más pobre que existe: deja pasar frecuencias por encima de Nyquist y produce un ligero aliasing residual en patrones muy finos. Un filtro de tienda o gaussiano se implementa ponderando cada muestra según su distancia al centro del píxel en lugar de promediarlas a peso igual. Empieza por la caja; el salto de calidad grande está en pasar de una muestra a muchas, no en el filtro.