About this technique →
quaternion-julia/sketch.js
// Quaternion Julia Set
//
// The Julia set extended to quaternion arithmetic (4D hypercomplex numbers).
// Iterate q = q² + c where q and c are quaternions; the bounded set is a 3D
// solid sliced from the full 4D shape. Rendered via distance-estimated raymarching.
//
// Quaternion multiplication: if q = (r, i, j, k), then
//   q² = (r²-i²-j²-k², 2ri, 2rj, 2rk)
// The derivative estimate ||dz/dr|| enables a smooth distance estimator.

export const PARAMS = {
  cR:         { value: -0.4,  min: -1,   max: 1,   step: 0.01, label: "c.r",                folder: "Constant" },
  cI:         { value: 0.6,   min: -1,   max: 1,   step: 0.01, label: "c.i",                folder: "Constant" },
  cJ:         { value: 0.0,   min: -1,   max: 1,   step: 0.01, label: "c.j",                folder: "Constant" },
  cK:         { value: 0.0,   min: -1,   max: 1,   step: 0.01, label: "c.k",                folder: "Constant" },
  iterations: { value: 8,     min: 4,    max: 16,  step: 1,    label: "Max Iterations",     folder: "Quality" },
  speed:      { value: 0.2,   min: 0,    max: 1.5, step: 0.01, label: "Camera Orbit Speed", folder: "Behavior" },
  hue:        { value: 280,   min: 0,    max: 360, step: 1,    label: "Hue",                folder: "Appearance" },
  glow:       { value: 0.5,   min: 0,    max: 1.5, step: 0.01, label: "Glow",               folder: "Appearance" },
};

export const SHARE = { bookmarked: ["cR", "cI", "cJ", "cK", "iterations", "hue", "glow"] };

export function shaderUniforms(params) {
  return {
    u_cR:         { value: params.cR },
    u_cI:         { value: params.cI },
    u_cJ:         { value: params.cJ },
    u_cK:         { value: params.cK },
    u_iterations: { value: params.iterations | 0 },
    u_speed:      { value: params.speed },
    u_hue:        { value: params.hue },
    u_glow:       { value: params.glow },
  };
}

export function shaderAnimate(uniforms, params) {
  uniforms.u_cR.value         = params.cR;
  uniforms.u_cI.value         = params.cI;
  uniforms.u_cJ.value         = params.cJ;
  uniforms.u_cK.value         = params.cK;
  uniforms.u_iterations.value = params.iterations | 0;
  uniforms.u_speed.value      = params.speed;
  uniforms.u_hue.value        = params.hue;
  uniforms.u_glow.value       = params.glow;
}

export function fragmentShader() {
  return `
    uniform float u_time;
    uniform vec2  u_resolution;
    uniform float u_cR;
    uniform float u_cI;
    uniform float u_cJ;
    uniform float u_cK;
    uniform int   u_iterations;
    uniform float u_speed;
    uniform float u_hue;
    uniform float u_glow;

    // Quaternion multiply: (a) * (b)
    vec4 qmul(vec4 a, vec4 b) {
      return vec4(
        a.x*b.x - a.y*b.y - a.z*b.z - a.w*b.w,
        a.x*b.y + a.y*b.x + a.z*b.w - a.w*b.z,
        a.x*b.z - a.y*b.w + a.z*b.x + a.w*b.y,
        a.x*b.w + a.y*b.z - a.z*b.y + a.w*b.x
      );
    }

    // Quaternion Julia distance estimator
    // p is the 3D point; w=0 slice of 4D quaternion space
    float qjulia(vec3 p) {
      vec4 c = vec4(u_cR, u_cI, u_cJ, u_cK);
      vec4 z  = vec4(p, 0.0);
      vec4 zp = vec4(1.0, 0.0, 0.0, 0.0);  // dz/dr for derivative
      for (int i = 0; i < 16; i++) {
        if (i >= u_iterations) break;
        zp = 2.0 * qmul(z, zp);
        z  = qmul(z, z) + c;
        if (dot(z, z) > 4.0) break;
      }
      float r = length(z);
      float dr = length(zp);
      if (dr < 1e-6) return 0.0;
      return 0.5 * log(r) * r / dr;
    }

    // Tetrahedron-epsilon normal (4 samples, no branch)
    vec3 calcNormal(vec3 p) {
      const float h = 0.001;
      const vec2 k = vec2(1.0, -1.0);
      return normalize(
        k.xyy * qjulia(p + k.xyy * h) +
        k.yyx * qjulia(p + k.yyx * h) +
        k.yxy * qjulia(p + k.yxy * h) +
        k.xxx * qjulia(p + k.xxx * h)
      );
    }

    // Cosine palette — matches Mandelbulb style
    vec3 palette(float t, float baseHue) {
      float h = baseHue / 360.0 + t * 0.4;
      vec3 a = vec3(0.5);
      vec3 b = vec3(0.5);
      vec3 c = vec3(1.0);
      vec3 d = vec3(0.00, 0.10, 0.20) + h;
      return a + b * cos(6.2831853 * (c * t + d));
    }

    void main() {
      vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution) / u_resolution.y;

      // Orbiting camera — slow vertical bob to see the shape from angles
      float a = u_time * u_speed;
      vec3 ro = vec3(2.8 * cos(a), 0.8 * sin(a * 0.37), 2.8 * sin(a));
      vec3 ta = vec3(0.0);
      vec3 ww = normalize(ta - ro);
      vec3 uu = normalize(cross(vec3(0.0, 1.0, 0.0), ww));
      vec3 vv = cross(ww, uu);
      vec3 rd = normalize(uv.x * uu + uv.y * vv + 1.5 * ww);

      // Raymarch
      float t = 0.0;
      float steps = 0.0;
      bool hit = false;
      for (int i = 0; i < 200; i++) {
        vec3 p = ro + rd * t;
        float d = qjulia(p);
        if (d < 0.0005) { hit = true; break; }
        if (t > 8.0) break;
        t += max(d, 0.001);
        steps += 1.0;
      }

      vec3 col = vec3(0.02, 0.02, 0.04);
      if (hit) {
        vec3 p  = ro + rd * t;
        vec3 n  = calcNormal(p);
        vec3 ld = normalize(vec3(1.5, 2.0, 1.0));
        float diff = max(0.0, dot(n, ld));
        float spec = pow(max(0.0, dot(reflect(-ld, n), -rd)), 48.0);
        vec3 base = palette(0.5 + 0.5 * cos(t * 0.8), u_hue);
        col = base * (0.2 + 0.8 * diff) + vec3(spec) * 0.4;
      }
      col += u_glow * vec3(steps / 200.0) * palette(steps / 200.0, u_hue);
      col = pow(max(col, vec3(0.0)), vec3(1.0 / 2.2));
      gl_FragColor = vec4(col, 1.0);
    }
  `;
}