Ir al contenido

Introducción al Raymarching

He estado siguiendo la masterclass de Iñigo Quilez «Live Coding “Happy Jumping”», donde implementa en ShaderToy un shader con raymarching que genera una escena con un monstruo saltarín paseando por un terreno peculiar. Para él era un ejercicio de animación, ya que aplica principios como el squash and stretch.

Para mí, en cambio, es una oportunidad para aprender cómo se construyen estas escenas tridimensionales usando shaders, sin modelos 3D, solo con matemáticas y código. Con este artículo intento expandir algunos conceptos para que queden más claros y profundizar un poco más en el tema.

Ver el shader “Happy Jumping” en ShaderToy

Raymarching

Según Iñigo Quilez, el raymarching puede entenderse como una técnica para explorar una escena tridimensional mediante rayos. Cada píxel de la pantalla genera un rayo que se adentra en el espacio virtual y busca las superficies presentes en la escena.

Esquema de un rayo explorando una escena tridimensional
Con raymarching, la escena se describe mediante distancias en lugar de calcular de forma exacta la intersección entre un rayo y la geometría (ray tracing). El algoritmo avanza de forma iterativa a lo largo del rayo utilizando la distancia estimada hasta la superficie más cercana.

En cada iteración se evalúa una Signed Distance Function (SDF), una función que devuelve la distancia firmada hasta la superficie más cercana. Este valor determina cuánto puede avanzar el rayo de forma segura sin atravesar ningún objeto. Cuando la superficie está lejos, los pasos pueden ser grandes; a medida que el rayo se aproxima a un objeto, los incrementos se reducen progresivamente, aumentando la precisión del cálculo.

Avance iterativo del rayo usando la distancia estimada a la superficie

De GPU Gems 2: Chapter 8.

Gracias a este mecanismo, es posible representar escenas complejas definidas únicamente mediante funciones matemáticas, sin necesidad de mallas, vértices o estructuras geométricas tradicionales.

// This normalize the coordinates from (-aspectRatio,-1) to (aspectRatio, 1)
vec2 p = (2.0*(fragCoord+offset)-iResolution.xy) / iResolution.y;

vec3 ro = vec3(0.0,0.0,2.0);
vec3 rd = normalize(vec3(p,-1.5));

// if ray doesn't hit anything value -1
vec4 res = vec4(-1.0, -1.0, 0.0,1.0);

float t = 0.0;
for(int i = 0; i < 128; ++i)
{
	vec3 pos = ro + t* rd;
	// map function will return a distance that is secure to advance
	vec4 h = map(pos,atime);
	t+=h.x;
}

Pseudocódigo del raymarching

Ahora bien, ¿cómo se calcula esa distancia máxima que se puede avanzar sin atravesar los distintos objetos? No hay una respuesta única: dependiendo del tipo de objeto existen distintas fórmulas. Una esfera seria un caso simple:

Cálculo de la distancia h hasta una esfera
Demostración del cálculo de h

La distancia máxima (h) a la que podemos avanzar sin atravesar la esfera es la distancia entre el origen del rayo y el centro de la esfera, menos el radio de esta. Con esto ya podemos empezar a pintar algo: en cuanto sabemos que hemos dado con la esfera, podemos mostrarla en pantalla e incluso empezar a iluminarla.

Primer render de una esfera con raymarching
Primer render de una esfera con raymarching

Calcular normales

Para poder calcular la iluminación de la escena, primero debemos calcular las normales de los objetos con los que colisiona el rayo. Cuando sabemos a qué distancia t se produce la colisión, podemos hallar el punto exacto de la superficie, y a partir de él calcular la normal.

vec3 pos = ro + t*rd;
vec3 nor = calcNormal( pos, time );

Como la función map() es continua, la dirección en la que la distancia cambia más rápido (su gradiente) coincide con la normal de la superficie. Para aproximar ese gradiente, tomamos el punto de colisión y lo desplazamos ligeramente en cada eje (x, y, z), calculando la diferencia de distancia entre desplazar el punto hacia un lado y hacia el otro (diferencias finitas centradas). El resultado, normalizado, nos da la dirección de la normal.

vec3 calcNormal( in vec3 pos, float time ) 
{
	float e = 0.0001; // pequeño desplazamiento
	
	float dx = map( pos + vec3(e,0,0), time ).x - map( pos - vec3(e,0,0), time ).x;
	float dy = map( pos + vec3(0,e,0), time ).x - map( pos - vec3(0,e,0), time ).x;
	float dz = map( pos + vec3(0,0,e), time ).x - map( pos - vec3(0,0,e), time ).x;
	return normalize( vec3(dx, dy, dz) ); 
}

Diferencia finita centrada: en vez de comparar el punto de colisión con un único punto desplazado, comparamos dos puntos situados simétricamente a ambos lados (pos+e y pos-e). Esto da una aproximación más precisa de la derivada que si solo miráramos hacia un lado, ya que promedia la variación en ambas direcciones.

Cálculo de luz

Sun Diffuse

Para calcular cuánta luz directa del sol recibe un punto de la superficie, usamos la normal del objeto y la dirección hacia la luz (el sol, en este caso). Calculamos el producto escalar (dot product) entre ambos vectores; como los dos son vectores unitarios (normalizados), el resultado es directamente el coseno del ángulo que forman.

Producto escalar entre la normal y la dirección de la luz

Este coseno es lo que se conoce como ley del coseno de Lambert: cuanto más alineada esté la normal con la dirección de la luz (ángulo cercano a 0°), más luz recibe la superficie (coseno cercano a 1); cuanto más de lado reciba la luz (ángulo cercano a 90°), menos luz le llega (coseno cercano a 0). Esto simula cómo la misma cantidad de luz se “reparte” sobre más o menos superficie según el ángulo de incidencia.

vec3 sun_dir = normalize(vec3(0.6, 0.35, 0.5));
float sun_dif = clamp(dot(sun_dir, normal), 0.0, 1.0);

Hacemos clamp(..., 0.0, 1.0) porque, si el ángulo supera los 90° (la superficie mira en dirección opuesta al sol), el coseno se vuelve negativo — y una cantidad de luz negativa no tiene sentido físico. Acotando el valor entre 0 y 1 nos aseguramos de que esas superficies simplemente no reciban luz directa del sol (quedan en sombra propia), en vez de “restar” luz de forma incorrecta.

Sky Diffuse & Bounce Diffuse

Para darle una iluminación del cielo y del suelo (la luz que rebota), utilizamos la misma técnica que con la luz del sol.

Hacemos el dot product de la normal con el vector vertical: positivo (0,1,0) para el cielo, negativo (0,-1,0) para el rebote del suelo. También hacemos un remapeo para que los valores vayan de 0 a 1, ya que el dot product devuelve un valor entre -1 y 1. En sky_dif usamos 0.5 + 0.5*x, que reparte el rango de forma simétrica (mitad y mitad); en bou_dif usamos 0.1 + 0.9*x, una mezcla distinta que da más peso al dot product (rebote más intenso cuando la normal mira directamente al suelo) y menos peso a un valor base constante. Por último, el resultado se clampea para que no se salga del rango [0.0, 1.0].

float sky_dif = clamp(0.5 + 0.5*dot(normal, vec3(0.0, 1.0, 0.0)), 0.0, 1.0);
float bou_dif = clamp(0.1 + 0.9*dot(normal, vec3(0.0, -1.0, 0.0)), 0.0, 1.0);

Esto se puede simplificar, ya que multiplicar por un vector unitario que solo tiene una componente activa simplemente “extrae” esa componente. Por tanto, podemos usar normal.y directamente en vez de hacer el dot product:

float sky_dif = clamp(0.5 + 0.5*normal.y, 0.0, 1.0);
float bou_dif = clamp(0.1 - 0.9*normal.y, 0.0, 1.0);

Hard Shadows

La idea de calcular sombras en raymarching es bastante sencilla. Reutilizamos la misma función castRay que usamos para el rayo de la cámara, pero ahora lanzamos el rayo desde la superficie del objeto (ligeramente desplazado a lo largo de la normal, para no chocar consigo mismo por errores de precisión) en dirección a la luz. Si ese rayo golpea algo antes de llegar a la fuente de luz, el punto original está en sombra.

6_Raymarching.png

float sun_sha = step(castRay(pos + 0.001*nor, sun_lig, time), 0.0);

step(valor, 0.0) — la función step(edge, x) devuelve 0.0 si x < edge, y 1.0 si x >= edge. Aquí se usa “al revés” de como se suele leer: step(castRay(...), 0.0) compara 0.0 contra el resultado de castRay. Como en esta versión simplificada castRay devuelve -1.0 cuando no hay impacto:

  • Si no hay impacto → castRay(...) = -1.0step(-1.0, 0.0) = 1.0 (0.0 ≥ -1.0) → sin sombra, luz plena.
  • Si hay impacto (choca con algo) → castRay(...) es una distancia positiva (por ejemplo 3.5) → step(3.5, 0.0) = 0.0 (0.0 < 3.5) → en sombra total.

Calculo de la oclusión

Esta función simula un efecto muy sutil pero importante para el realismo: las zonas donde hay geometría cercana rodeando un punto (rincones, pliegues, huecos entre formas) reciben menos luz ambiental, porque parte de esa luz queda bloqueada por los propios objetos vecinos. Es lo que hace que las axilas de un personaje, o el hueco entre dos dedos, se vean ligeramente más oscuros aunque no haya una sombra proyectada por el sol.

float h = 0.01 + 0.11*float(i)/4.0; 
vec3 opos = pos + h*nor;
float d = map( opos, time ).x;
occ += (h-d)*sca;

Para cada una de las 5 iteraciones, tomamos el punto de la superficie (pos) y nos movemos una distancia h a lo largo de la normal (nor) — es decir, nos alejamos perpendicularmente hacia afuera del objeto. Si no hubiera nada más cerca (superficie totalmente convexa y despejada), la distancia real hasta cualquier otro objeto en ese punto (d, obtenida de map()) sería exactamente h — porque no hay ningún obstáculo interfiriendo.

Pero si hay geometría vecina cerca (un rincón, un pliegue), la distancia real d será menor que h, porque hay algo estorbando antes de esa distancia esperada. La diferencia h - d mide justo eso: cuánto “más cerca de lo esperado” hay superficie ahí — cuanta más geometría vecina haya, mayor será esta diferencia.

return clamp( 1.0 - 2.0*occ, 0.0, 1.0 );

occ acumula las diferencias (h-d) ponderadas: cuanto mayor sea occ, más “encerrado” está el punto por geometría vecina. Al hacer 1.0 - 2.0*occ, invertimos la escala: si occ es alto (mucha oclusión), el resultado se acerca a 0 (poca luz ambiental); si occ es bajo o cero (superficie despejada), el resultado se acerca a 1 (luz ambiental completa). El clamp asegura que el resultado quede siempre en el rango [0, 1], ya que la multiplicación por 2.0 es un factor arbitrario para amplificar el efecto (ajuste artístico, no físico) y podría sacar el valor fuera de rango.

Camera (Look-at)

Para construir una cámara que siempre mira a un punto fijo. Se define un target, y a partir de ahí la posición de la cámara (que también será el origen del rayo, ro).

Con el target y el origen podemos definir el vector “adelante” (ww), normalizado: es simplemente la dirección desde ro hacia target.

Para generar el vector de la “derecha” (uu), hacemos el producto vectorial del vector adelante y el vector vertical del mundo. El producto vectorial (cross) entre dos vectores da un tercer vector perpendicular a ambos. Al hacer cross(ww, vec3(0,1,0)) (adelante × arriba-del-mundo), obtenemos un vector perpendicular tanto a “hacia dónde mira la cámara” como al eje vertical del mundo — que es exactamente la dirección “derecha” de la cámara.

Una vez tenemos los vectores adelante y derecha, repetimos el proceso para obtener el vector “arriba” real de la cámara (vv). No podemos simplemente reutilizar vec3(0,1,0) como arriba, porque si la cámara mira hacia arriba o hacia abajo, ese “arriba del mundo” ya no sería perpendicular a los otros dos ejes. Por eso lo recalculamos con cross(uu, ww), asegurando que los tres vectores formen una base perfectamente ortogonal entre sí.

        vec3 target = vec3(0.0, 0.75, 0.4);
        // ray origin
        vec3 ro = target + vec3(1.5 * sin(an),-.1,1.5 * cos(an));
        
        vec3 ww = normalize(target - ro);
        vec3 uu = normalize(cross(ww, vec3(0,1,0)));
        vec3 vv = normalize(cross(uu,ww));

Estos vectores forman la base para construir la dirección de cada rayo por píxel:

// This normalize the coordinates from [-aspectRatio,-1] to [aspectRatio, 1]
vec2 p = (2.0*(fragCoord+offset) - iResolution.xy) / iResolution.y;
vec3 rd = normalize(p.x*uu + p.y*vv + 1.8*ww);

Aquí p.x y p.y son las coordenadas normalizadas del píxel en pantalla: p.y va de -1 a 1, mientras que p.x va de -aspectRatio a aspectRatio (donde aspectRatio es el ancho entre el alto de la pantalla). Estas coordenadas se combinan con los ejes uu (derecha), vv (arriba) y ww (adelante) de la cámara para construir la dirección del rayo en el espacio del mundo — es decir, “cuánto me desvío a la derecha/arriba desde el centro de la vista, y cuánto avanzo hacia adelante”. El 1.8 que multiplica a ww controla el campo de visión (FOV): cuanto mayor sea ese número, más “zoom” (campo de visión más estrecho); cuanto menor, más angular (FOV más amplio).

Escena simple

Este es el resultado final: el shader de Iñigo Quilez con una escena básica que reúne los elementos que hemos ido explicando — el bucle de raymarching, el cálculo de normales, y la iluminación combinando la luz del sol, del cielo, del rebote y las sombras.

Shader de Iñigo Quilez ShaderToy