Esta es la continuación espiritual de shaders en corto. Aquí describiré los pasos para renderizar Mandelbulb en GLSL.

Qué es Mandelbulb? un fractal en 3-dimensiones. De nada.

En el post previo se muestra Mandelbrot (un fractal en 2-dimensiones) como se muestra a continuación:

Mandelbulb tiene la intención de mostrar como sería Mandelbrot en 3-dimensiones. Y si me equivoco, pues me equivoqué.

Trayectoria:

• Contexto: El objetivo es renderizar Mandelbulb en GLSL, por diversión 👍
• Limitaciones: Conocimientos de ray marching prematuros y de álgebra lineal no muy remotos. (lado positivo: algunas cosas son caja negra)
• Retroalimentación: No hay 👍

La intención es construir una versión propia, separada por partes (a mi manera (es una dictadura)).

Ray marching en GLSL

El ray marching es la técnica que se usará para renderizar Mandelbulb. Existen otras, es probable que también sirvan, pero parto de la asunción de que es posible ya que existen muchas implementaciones.

Esta técnica es iterativa, busca la intersección entre un supuesto rayo de luz y un objeto. En otras palabras, se lanza un rayo de luz desde un punto de la pantalla hasta encontrar un objeto.

Para aplicar la técnica en GLSL la parte importante es el fragment shader.

El fragment shader es el siguiente:

precision highp float;
varying vec2 uv;
uniform float uTime;

mat3 rotationX(float angle) {
  float c = cos(angle);
  float s = sin(angle);
  return mat3(
    1.0, 0.0, 0.0,
    0.0, c, -s,
    0.0, s, c
  );
}

mat3 rotationY(float angle) {
  float c = cos(angle);
  float s = sin(angle);
  return mat3(
    c, 0.0, s,
    0.0, 1.0, 0.0,
    -s, 0.0, c
  );
}

float sdCube(vec3 p) {
  vec3 boxSize = vec3(1.0, 1.0, 1.0);
  vec3 d = abs(p) - boxSize;
  return length(max(d, 0.0)) + min(max(d.x, max(d.y, d.z)), 0.0);
}

float distanceField(vec3 p) {
  float angle = uTime;
  p = rotationY(angle) * p;
  p = rotationX(angle) * p;
  return sdCube(p);
}

vec3 rayMarch(vec3 ro, vec3 rd) {

  // Menor valor oculta la escena de de atrás hacia adelante
  const float maxDistance = 100.0;

  // Menor valor oculta las esquinas de los objetos
  const float minDistance = 0.001;

  // Mayor valor aumenta la fidelidad de la imágen a mayor coste computacional
  const int maxIter = 100;

  float totalDistance = 0.0;
  vec3 p;
  for (int i = 0; i < maxIter; i++) {
    p = ro + totalDistance * rd;
    float d = distanceField(p);
    if (d < minDistance) {
      return vec3(p); // Superficie
    }
    totalDistance += d;
    if (totalDistance > maxDistance) {
      break;
    }
  }

  return vec3(0.0, 0.0, 0.0); // Fondo
}

void main() {
  vec3 ro = vec3(0.0, 0.0, 4.0); // Posición del observador
  vec3 rd = normalize(vec3(uv, -1.0)); // Dirección del rayo

  vec3 color = rayMarch(ro, rd);

  gl_FragColor = vec4(color, 1.0);
}

Olvida todo!

Son muchas cosas pero lo importante es rayMarch:

vec3 rayMarch(vec3 ro, vec3 rd) {

  // Menor valor oculta la escena de de atrás hacia adelante
  const float maxDistance = 100.0;

  // Menor valor oculta las esquinas de los objetos
  const float minDistance = 0.001;

  // Mayor valor aumenta la fidelidad de la imágen a mayor coste computacional
  const int maxIter = 100;

  float totalDistance = 0.0;
  vec3 p;
  for (int i = 0; i < 100; i++) {
    p = ro + totalDistance * rd;
    float d = distanceField(p);
    if (d < minDistance) {
      return vec3(p); // Superficie
    }
    totalDistance += d;
    if (totalDistance > maxDistance) {
      break;
    }
  }

  return vec3(0.0, 0.0, 0.0); // Fondo
}

A diferencia de otros métodos de renderizado, este es estimado. La parte más relevante es la que nos permite obtener la coordenada p. Se obtiene de forma iterativa, cuando la distancia cumple con la condición de ser menor a minDistance el punto p es considerado parte de la escena (de alguna superficie). Con ello se puede asignar un valor al color.

El resto del código son utilidades (álgebra lineal: posiciones, rayos de luz, rotaciones, etc...). Para los objetos existen muchas funciones con las que se puede calcular la distancia a uno de ellos, puedes revisar:

• https://iquilezles.org/articles/distfunctions/ ✨

En distanceField es donde se calculan las distancias a cada objeto, en este caso hace rotaciones y calcula la distancia a un cubo.

Mandelbulb

La referencia principal se toma del siguiente recurso:

• https://editor.p5js.org/Taxen99/sketches/47CDg5-nV by Taxen99

También sirve ver la referencia en wikipedia:

• https://wikipedia.org/wiki/Mandelbulb

Pasos:

1. Evaluar el punto en un área de 3-dimensiones.
2. Calcular si pertenece al conjunto dado un criterio.

Suena simple pero el criterio es extenso. Primero es necesario pasar la coordenada p a coordenadas esféricas:

vec3 toSpherical(vec3 p) {
  float r = length(p);
  float theta = acos(p.z / r);
  float phi = atan(p.y, p.x);
  return vec3(r, theta, phi);
}

No puedo decir mucho, la fórmula de White y Nylander hacen posible el cálculo del conjunto de Mandelbulb. ✨

Para no entrar en mucho detalle, empecemos del main:

void main() {
  vec3 ro = vec3(0.0, 0.0, 1.6);
  vec3 rd = normalize(vec3(uv, -1.0));
  vec4 color = rayMarch(ro, rd);
  gl_FragColor = color;
}

Es lo mismo que en ejemplos anteriores, se modificaron los valores para mostrar una visualización más apreciable.

rayMarch:

vec4 rayMarch(vec3 ro, vec3 rd) {
  const float maxDistance = 100.0;
  const float minDistance = 0.001;
  const int maxIter = 100;
  float totalDistance = 0.0;
  vec3 p;
  for (int i = 0; i < maxIter; i++) {
    p = ro + totalDistance * rd;
    vec2 dfVec = distanceField(p);
    float d = dfVec.x;
    if (d < minDistance) {
      float effort = dfVec.y;
      float g = effort;
      float b = 1.0 - effort;
      return vec4(0.0, g, b, 1.0);
    }
    totalDistance += d;
    if (totalDistance > maxDistance) {
      break;
    }
  }
  return vec4(0.0);
}

Se añadió un parámetro extra (effort) que representa cuantas iteraciones se realizaron para saber si una coordenada pertenece al conjunto (verde: muchas iteraciones, azul: pocas iteraciones).

distanceField:

vec2 distanceField(vec3 p) {
  float angle = uTime / 10.0;
  p = rotationY(angle) * p;
  p = rotationX(PI / 2.0) * p;
  return sdMandelbulb(p, 8.0);
}

Para apreciar el mejor el fractal, se hacen ligeras rotaciones.

sdMandelbulb: 💀

vec2 sdMandelbulb(vec3 p, float power) {
  const int maxIter = 12;

  vec3 zeta = p;
  float dr = 1.0;
  float r = 0.0;
  int iterations = 0;
  for (int _ = 0; _ < maxIter; _++) {
    iterations++;

    vec3 spherical = toSpherical(zeta);
    r = spherical.x;
    if (r > 2.0) {
      break;
    }

    dr = pow(r, power - 1.0) * power * dr + 1.0;

    float theta = spherical.y * power;
    float phi = spherical.z * power;
    float sinTheta = sin(theta);

    float powx = pow(r, power) * sinTheta * cos(phi);
    float powy = pow(r, power) * sinTheta * sin(phi);
    float powz = pow(r, power) * cos(theta);

    zeta.x = powx + p.x;
    zeta.y = powy + p.y;
    zeta.z = powz + p.z;
  }
  float d = 0.5 * log(r) * r / dr;
  float effort = float(iterations) / float(maxIter);
  return vec2(d, effort);
}

Esta función estima la distancia a una coordenada que pertenezca al conjunto de Mandelbulb. Es una mezcla de las referencias de los recursos mencionados y lo que pude entender de este fractal.

Hay algunas cosas que me generar ruido como dr o el cálculo de d. No las comprendo muy bien, pero entiendo que se busca normalizar la distancia estimada.

Adicionalmente, se añadió effort como una utilidad para distinguir entre los puntos del conjunto.

En este ejemplo se usa una potencia de 8, sin embargo, lo genial viene de probar distintos valores.

Experimentos

n = min(time, 8)

Este es muy similar al del post anterior.

n = 8, rotateXY(time)

Portada

Eso, el código con el que se creó la portada está en ese iframe.

Recursos

• Shaders en corto
• https://jamie-wong.com/2016/07/15/ray-marching-signed-distance-functions/
• https://wikipedia.org/wiki/Mandelbulb
• https://editor.p5js.org/Taxen99/sketches/47CDg5-nV

Nota: A nivel personal me queda pendiente experimentar con estados previos para calcular el siguiente, así como en el juego de la vida del post anterior. Igual lo dejo again (copiar es gratis):