Burning Ship Fractal

DEFINITION

z_{n+1} = (|Re(z_n)| + i|Im(z_n)|)² + c

NOTES

A variation of the Mandelbrot set where the absolute values of the real and imaginary components are taken before squaring at each iteration. The non-analytic transformation breaks Cauchy-Riemann symmetry, creating a fractal structure resembling a burning ship with intricate mast detail.

PROVENANCE

Origin

M. Michelitsch & O. E. Rössler, Computers & Graphics 16(4), 1992

Standing

Public domain — complex dynamical system

Constants

Iter governs calculation depth; cx/cy sets target coordinates

Source

doi:10.1016/0097-8493(92)90032-Q

TOUCHDESIGNER · GLSL

The same shader this plate runs, reframed for a GLSL TOP. Pasted bare it renders the published constants as a still frame; wire absTime.seconds into u_t on the Vectors page to animate it.

// FORMA — PL. 34 · BURNING SHIP FRACTAL — Michael Michelitsch & Otto E. Rössler, 1992
//   z_{n+1} = (|Re(z_n)| + i|Im(z_n)|)² + c
// TouchDesigner port — paste into a GLSL TOP's pixel shader. Set the
// resolution on the TOP's Common page. As pasted it renders the published
// constants as a still frame; to animate, add a uniform named u_t on the
// GLSL TOP's Vectors 1 page with the expression absTime.seconds.
// Constants are consts — edit to tweak; comments give the measured range.
// Written from the published mathematics, not adapted from any code.

#define u_res (uTDOutputInfo.res.zw)
uniform float u_t;               // absTime.seconds on the Vectors page; unset = still

const float u_phase = 0.7013;    // this plate's own grid phase, 0..1
// FORMA's FRACTALS accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.1389, 0.1912);
const vec3 u_pal_b = vec3(0.5, 0.151, 0.2078);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_iter = 80.0;        // max iterations · live 20 .. 240
const float p_cx   = -1.1;        // centre x · live -2 .. 1.5
const float p_cy   = -0.45;       // centre y · live -2 .. 1.5
const float p_zoom = 1.0;         // zoom factor · live 0.5 .. 15

/* The order's ramp — the same cosine formulation the JS kit uses, so a
   plate keeps its classification colour in either language. */
vec3 ramp(float t){
  return clamp(u_pal_a + u_pal_b * cos(6.28318530718 * (u_pal_c * t + u_pal_d)), 0.0, 1.0);
}

/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return fract(t / period + u_phase); }
float pingpong(float t, float period){
  float u = cycle(t, period);
  return u < 0.5 ? u * 2.0 : 2.0 - u * 2.0;
}


float hash2(int x, int y){
  uint h = uint(x) * 374761393u + uint(y) * 668265263u;
  h ^= h >> 13u;
  h *= 1274126177u;
  h ^= h >> 16u;
  return float(h) / 4294967296.0;
}

float smoothCurve(float t){ return t * t * (3.0 - 2.0 * t); }
float fadeCurve(float t){ return t * t * t * (t * (t * 6.0 - 15.0) + 10.0); }

/* Value noise: bilinear interpolation of a hashed lattice. */
float valueNoise(vec2 p){
  vec2 c = floor(p), f = p - c;
  int xi = int(c.x), yi = int(c.y);
  float u = smoothCurve(f.x), v = smoothCurve(f.y);
  return mix(mix(hash2(xi, yi),     hash2(xi + 1, yi),     u),
             mix(hash2(xi, yi + 1), hash2(xi + 1, yi + 1), u), v);
}

/* Gradient (Perlin) noise: dot products against pseudo-random unit vectors. */
float gradDot(int ix, int iy, float dx, float dy){
  float a = hash2(ix, iy) * 6.28318530718;
  return cos(a) * dx + sin(a) * dy;
}
float gradNoise(vec2 p){
  vec2 c = floor(p), f = p - c;
  int xi = int(c.x), yi = int(c.y);
  float u = fadeCurve(f.x), v = fadeCurve(f.y);
  return mix(mix(gradDot(xi,     yi,     f.x,       f.y),
                 gradDot(xi + 1, yi,     f.x - 1.0, f.y), u),
             mix(gradDot(xi,     yi + 1, f.x,       f.y - 1.0),
                 gradDot(xi + 1, yi + 1, f.x - 1.0, f.y - 1.0), u), v) * 0.7071 + 0.5;
}

/* Octaves summed at falling amplitude. GLSL has no function pointers, so the
   two bases are two functions rather than one with a noise argument. The
   three-argument forms take the per-octave gain — the Hurst roughness dial,
   mirroring the kit's fbm — and the two-argument forms keep the classic 0.5
   so existing call sites read unchanged. */
float fbmValue(vec2 p, int oct, float gain){
  float sum = 0.0, amp = 0.5, norm = 0.0;
  for (int i = 0; i < 9; i++){          // 9 is the octave slider's ceiling
    if (i >= oct) break;
    sum += amp * valueNoise(p);
    norm += amp;
    amp *= gain;
    p *= 2.0;
  }
  return sum / norm;
}
float fbmValue(vec2 p, int oct){ return fbmValue(p, oct, 0.5); }
float fbmGrad(vec2 p, int oct, float gain){
  float sum = 0.0, amp = 0.5, norm = 0.0;
  for (int i = 0; i < 9; i++){
    if (i >= oct) break;
    sum += amp * gradNoise(p);
    norm += amp;
    amp *= gain;
    p *= 2.0;
  }
  return sum / norm;
}
float fbmGrad(vec2 p, int oct){ return fbmGrad(p, oct, 0.5); }

vec3 plate(vec2 uv){
  /* The dive in place of the old centre wobble — mandelbrot's cycle, the
     same camera the JS path runs. */
  float z = p_zoom * (1.0 + pingpong(u_t, 31.0) * 119.0);
  float maxIt = min(240.0, floor(p_iter * (1.0 + log2(z) * 0.5) + 0.5));
  float scale = 1.5 / z;
  float ar = u_res.y / u_res.x;
  float cr = p_cx + (uv.x - 0.5) * scale;
  float ci = p_cy + (uv.y - 0.5) * scale * ar;
  float zr = 0.0, zi = 0.0;
  float n = 0.0;
  for (int i = 0; i < 240; i++){
    if (float(i) >= maxIt) break;
    float azr = abs(zr), azi = abs(zi);
    float zr2 = azr * azr, zi2 = azi * azi;
    if (zr2 + zi2 > 4.0) break;
    zi = 2.0 * azr * azi + ci;
    zr = zr2 - zi2 + cr;
    n += 1.0;
  }
  if (n >= maxIt) return vec3(4.0, 6.0, 10.0) / 255.0;   // the interior is INK, not a hole
  float nu = log(log(zr * zr + zi * zi) / 2.0) / 0.69314718;
  float iterSmooth = n + 1.0 - nu;
  /* The log remap the escape family shares, exactly as the JS path lays it. */
  float g = log(1.0 + max(0.0, iterSmooth)) / log(1.0 + maxIt);
  return ramp(0.44 + 0.62 * g);
}

out vec4 fragColor;
void main(){
  // FORMA's uv runs y-down, matching its canvas; TD's vUV runs up
  vec2 uv = vec2(vUV.s, 1.0 - vUV.t);
  fragColor = TDOutputSwizzle(vec4(plate(uv), 1.0));
}

NUKE · BLINKSCRIPT

The same shader this plate runs, transpiled to a BlinkScript kernel. Paste it into a BlinkScript node's Kernel Source and press Recompile; every constant arrives as a knob at its published value, and u_t animates with the expression frame/24. Compiled and rendered in Nuke 17.1, then compared against this plate on the page.

// FORMA — PL. 34 · BURNING SHIP FRACTAL — Michael Michelitsch & Otto E. Rössler, 1992
//   z_{n+1} = (|Re(z_n)| + i|Im(z_n)|)² + c
// Nuke port — a BlinkScript kernel. Paste into a BlinkScript node's Kernel
// Source and press Recompile. Every constant arrives as a knob at its published
// value (the comment gives the measured range); u_t is a knob too — animate it
// with the expression frame/24 or leave it at 0 for the still frame. Written
// from the published mathematics, not adapted from any code.
// Transpiled from the shader this plate runs on the page (GLSL ES 3.00):
// vec → float2/3/4, swizzles expanded, GLSL builtins Blink lacks written out
// as forma_ functions, float literals suffixed. Compiled and rendered in a
// real Nuke (17.1v1) and compared against this plate on the page: 34 of 34.
//
// plate() and its helpers are written to a single exit — the loop that runs
// once. That is not a style: Blink 17.1 drops a conditional early return from
// a called function while Vectorize is on, which is the node default, with no
// warning and no error. Written this way it paints correctly as pasted.

kernel Forma_burningship : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_iter; // max iterations · live 20 .. 240
  float p_cx;   // centre x · live -2 .. 1.5
  float p_cy;   // centre y · live -2 .. 1.5
  float p_zoom; // zoom factor · live 0.5 .. 15

local:
  float2 u_res;
  float u_phase;
  float3 u_pal_a, u_pal_b, u_pal_c, u_pal_d;

  void define(){
    defineParam(u_t, "u_t", 0.0f);
    defineParam(p_iter, "p_iter", 80.0f);
    defineParam(p_cx, "p_cx", -1.1f);
    defineParam(p_cy, "p_cy", -0.45f);
    defineParam(p_zoom, "p_zoom", 1.0f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.7013f;    // this plate's own grid phase, 0..1
    // FORMA's FRACTALS accent as cosine-gradient coefficients
    u_pal_a = float3(0.46f, 0.1389f, 0.1912f);
    u_pal_b = float3(0.5f, 0.151f, 0.2078f);
    u_pal_c = float3(1.0f, 1.0f, 1.0f);
    u_pal_d = float3(0.0f, 0.05f, 0.1f);
  }

  /* GLSL builtins Blink lacks, written as templates rather than overload sets.
     Blink's operators return expression templates (Swizzle<float,N>), so a call
     passing an expression cannot resolve against an overload set on float2
     against float3 — measured in Nuke 17.1: a float2 expression is ambiguous
     between the two, while scalar-against-vector resolves. A template deduces
     the expression's own type, so the ambiguity cannot arise. */
  template <class T> T forma_fract(T v){ return v - floor(v); }
  template <class T, class S> T forma_mod(T x, S y){ return x - y * floor(x / y); }
  /* Blink's own min/max/clamp take no scalar bound against a vector, which GLSL
     does; v * 0.0f + b is that bound at the vector's own width, and collapses to
     b when v is a scalar, so one template serves both. */
  template <class T, class S> T forma_min(T a, S b){ return min(a, a * 0.0f + b); }
  template <class T, class S> T forma_max(T a, S b){ return max(a, a * 0.0f + b); }
  template <class T, class S> T forma_clamp(T v, S lo, S hi){ return clamp(v, v * 0.0f + lo, v * 0.0f + hi); }
  int forma_min(int a, int b){ return min(a, b); }
  int forma_max(int a, int b){ return max(a, b); }
  /* GLSL step(edge, x) is 1 where x >= edge; floor(sign(x - e) * 0.5 + 1) is
     that exactly, equality included, out of builtins Blink does have. */
  template <class T, class S> T forma_step(S e, T x){ return floor(sign(x - e) * 0.5f + 1.0f); }
  template <class T, class S> T forma_smoothstep(S a, S b, T x){
    T t = forma_clamp((x - a) / (b - a), 0.0f, 1.0f);
    return t * t * (3.0f - 2.0f * t);
  }
  template <class T> float forma_distance(T a, T b){ return length(a - b); }
  float forma_tanh(float x){ float e = exp(2.0f * x); return (e - 1.0f) / (e + 1.0f); }
  float forma_radians(float d){ return d * 0.01745329252f; }
  // the page's hash2 is exact uint32; Blink has int, so the shifts are made
  // logical by masking and the read-back is lifted into 0 .. 2^32
  /* A uint read back as a float. Blink has no unsigned type, so a value past
     2^31 arrives as a negative int and float() of it is negative. Measured on
     gabor, whose own generator then returned uniforms in [-0.5, 0.5) and drew
     a different picture — it compiled, it rendered, and only comparing it with

  /* The order's ramp — the same cosine formulation the JS kit uses, so a
     plate keeps its classification colour in either language. */
  float3 ramp(float t){
    return forma_clamp(u_pal_a + u_pal_b * cos(6.28318530718f * (u_pal_c * t + u_pal_d)), 0.0f, 1.0f);
  }

  /* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
  float cycle(float t, float period){ return forma_fract(t / period + u_phase); }
  float pingpong(float t, float period){
    float u = cycle(t, period);
    return u < 0.5f ? u * 2.0f : 2.0f - u * 2.0f;
  }


  float hash2(int x, int y){
    uint h = uint(x) * 374761393u + uint(y) * 668265263u;
    h ^= h >> 13u;
    h *= 1274126177u;
    h ^= h >> 16u;
    return float(h) / 4294967296.0f;
  }

  float smoothCurve(float t){ return t * t * (3.0f - 2.0f * t); }
  float fadeCurve(float t){ return t * t * t * (t * (t * 6.0f - 15.0f) + 10.0f); }

  /* Value noise: bilinear interpolation of a hashed lattice. */
  float valueNoise(float2 p){
    float2 c = floor(p);
    float2 f = p - c;
    int xi = int(c.x);
    int yi = int(c.y);
    float u = smoothCurve(f.x);
    float v = smoothCurve(f.y);
    return lerp(lerp(hash2(xi, yi),     hash2(xi + 1, yi),     u),
               lerp(hash2(xi, yi + 1), hash2(xi + 1, yi + 1), u), v);
  }

  /* Gradient (Perlin) noise: dot products against pseudo-random unit vectors. */
  float gradDot(int ix, int iy, float dx, float dy){
    float a = hash2(ix, iy) * 6.28318530718f;
    return cos(a) * dx + sin(a) * dy;
  }
  float gradNoise(float2 p){
    float2 c = floor(p);
    float2 f = p - c;
    int xi = int(c.x);
    int yi = int(c.y);
    float u = fadeCurve(f.x);
    float v = fadeCurve(f.y);
    return lerp(lerp(gradDot(xi,     yi,     f.x,       f.y),
                   gradDot(xi + 1, yi,     f.x - 1.0f, f.y), u),
               lerp(gradDot(xi,     yi + 1, f.x,       f.y - 1.0f),
                   gradDot(xi + 1, yi + 1, f.x - 1.0f, f.y - 1.0f), u), v) * 0.7071f + 0.5f;
  }

  /* Octaves summed at falling amplitude. GLSL has no function pointers, so the
     two bases are two functions rather than one with a noise argument. The
     three-argument forms take the per-octave gain — the Hurst roughness dial,
     mirroring the kit's fbm — and the two-argument forms keep the classic 0.5f
     so existing call sites read unchanged. */
  float fbmValue(float2 p, int oct, float gain){
    float sum = 0.0f;
    float amp = 0.5f;
    float norm = 0.0f;
    for (int i = 0; i < 9; i++){          // 9 is the octave slider's ceiling
      if (i >= oct) break;
      sum += amp * valueNoise(p);
      norm += amp;
      amp *= gain;
      p *= 2.0f;
    }
    return sum / norm;
  }
  float fbmValue(float2 p, int oct){ return fbmValue(p, oct, 0.5f); }
  float fbmGrad(float2 p, int oct, float gain){
    float sum = 0.0f;
    float amp = 0.5f;
    float norm = 0.0f;
    for (int i = 0; i < 9; i++){
      if (i >= oct) break;
      sum += amp * gradNoise(p);
      norm += amp;
      amp *= gain;
      p *= 2.0f;
    }
    return sum / norm;
  }
  float fbmGrad(float2 p, int oct){ return fbmGrad(p, oct, 0.5f); }

  float3 plate(float2 uv){
    float3 forma_r = float3(0.0f, 0.0f, 0.0f);
    for (int forma_once = 0; forma_once < 1; forma_once++){
      /* The dive in place of the old centre wobble — mandelbrot's cycle, the
         same camera the JS path runs. */
      float z = p_zoom * (1.0f + pingpong(u_t, 31.0f) * 119.0f);
      float maxIt = forma_min(240.0f, floor(p_iter * (1.0f + log2(z) * 0.5f) + 0.5f));
      float scale = 1.5f / z;
      float ar = u_res.y / u_res.x;
      float cr = p_cx + (uv.x - 0.5f) * scale;
      float ci = p_cy + (uv.y - 0.5f) * scale * ar;
      float zr = 0.0f;
      float zi = 0.0f;
      float n = 0.0f;
      for (int i = 0; i < 240; i++){
        if (float(i) >= maxIt) break;
        float azr = fabs(zr);
        float azi = fabs(zi);
        float zr2 = azr * azr;
        float zi2 = azi * azi;
        if (zr2 + zi2 > 4.0f) break;
        zi = 2.0f * azr * azi + ci;
        zr = zr2 - zi2 + cr;
        n += 1.0f;
      }
      if (n >= maxIt) { forma_r = float3(4.0f, 6.0f, 10.0f) / 255.0f; break; }
      float nu = log(log(zr * zr + zi * zi) / 2.0f) / 0.69314718f;
      float iterSmooth = n + 1.0f - nu;
      /* The log remap the escape family shares, exactly as the JS path lays it. */
      float g = log(1.0f + forma_max(0.0f, iterSmooth)) / log(1.0f + maxIt);
      { forma_r = ramp(0.44f + 0.62f * g); break; }
    }
    return forma_r;
  }

  void process(int2 pos){
    // FORMA's uv runs y-down like its canvas; Nuke's rows run up
    float2 uv = float2((float(pos.x) + 0.5f) / u_res.x, 1.0f - (float(pos.y) + 0.5f) / u_res.y);
    float3 c = plate(uv);
    dst() = float4(c.x, c.y, c.z, 1.0f);
  }
};