Skip to the plate
FORMA PUBLIC DOMAIN GENERATIVE ATLAS / ED. 0.28
Plate 29, Voronoi Diagram: a still of the partition / nearest site plate as the atlas renders it, in the lattices accent.

PL. 29  ·  LATTICES / PARTITION / NEAREST SITE

Voronoi Diagram

Georgy Voronoy, 1908 · Descartes, 1644

OPEN THE LIVE PLATE ▸

DEFINITION

cell(sᵢ) = { p : ‖p − sᵢ‖ ≤ ‖p − sⱼ‖ ∀j }

NOTES

Divide the plane so every point belongs to whichever site is nearest. Descartes sketched the idea for the distribution of matter in the solar system; Voronoy formalised it. It is the dual of the Delaunay triangulation and shows up wherever territory is decided by proximity — cell walls, crystal grains, service areas, giraffe coats.

PROVENANCE

Origin
G. Voronoy, 1908; anticipated by Descartes (1644) and Dirichlet
Standing
Public domain
Constants
Sites drift on slow circular orbits here

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. 29 · VORONOI DIAGRAM — Georgy Voronoy, 1908 · Descartes, 1644
//   cell(sᵢ) = { p : ‖p − sᵢ‖ ≤ ‖p − sⱼ‖ ∀j }
// 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.66;    // this plate's own grid phase, 0..1
// FORMA's LATTICES accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.4185, 0.1389);
const vec3 u_pal_b = vec3(0.5, 0.4549, 0.151);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_sites = 26.0;        // sites · live 4 .. 80
const float p_edges = 0.55;        // edge emphasis · live 0 .. 1
const float p_sharp = 9.0;         // edge sharpness · live 2 .. 24

/* 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){
  float gridN = sqrt(p_sites * 2.5);
  vec2 st = uv * gridN;
  ivec2 ipos = ivec2(floor(st));
  vec2 fpos = fract(st);
  float m_dist = 8.0;
  float m_dist2 = 8.0;
  ivec2 m_tile = ipos;
  for (int y = -1; y <= 1; y++){
    for (int x = -1; x <= 1; x++){
      ivec2 tile = ipos + ivec2(x, y);
      float h1 = hash2(tile.x, tile.y);
      float h2 = hash2(tile.x + 13, tile.y + 37);
      float ph = h1 * 6.283 + u_phase * 6.283;
      /* Drift amplitude and edge gain both follow the JS path: 0.04 of the
         plate converts to cells by gridN, and the sharpness divides back out
         so both tabs draw the same seam width. The two used to disagree
         (0.08 cells and a gain of 4), which meant the tabs showed different
         diagrams. */
      vec2 point = vec2(float(x), float(y)) + vec2(h1 + cos(u_t * 0.13 + ph) * 0.04 * gridN, h2 + sin(u_t * 0.13 + ph) * 0.04 * gridN);
      float dist = length(point - fpos);
      if (dist < m_dist){
        m_dist2 = m_dist;
        m_dist = dist;
        m_tile = tile;
      } else if (dist < m_dist2){
        m_dist2 = dist;
      }
    }
  }
  float edge = clamp((m_dist2 - m_dist) * p_sharp / gridN, 0.0, 1.0);
  float val = hash2(m_tile.x, m_tile.y * 3);
  return ramp(val * 0.5 + 0.1 + (1.0 - edge) * p_edges * 0.6);
}

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. 29 · VORONOI DIAGRAM — Georgy Voronoy, 1908 · Descartes, 1644
//   cell(sᵢ) = { p : ‖p − sᵢ‖ ≤ ‖p − sⱼ‖ ∀j }
// 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_voronoi : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_sites; // sites · live 4 .. 80
  float p_edges; // edge emphasis · live 0 .. 1
  float p_sharp; // edge sharpness · live 2 .. 24

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_sites, "p_sites", 26.0f);
    defineParam(p_edges, "p_edges", 0.55f);
    defineParam(p_sharp, "p_sharp", 9.0f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.66f;    // this plate's own grid phase, 0..1
    // FORMA's LATTICES accent as cosine-gradient coefficients
    u_pal_a = float3(0.46f, 0.4185f, 0.1389f);
    u_pal_b = float3(0.5f, 0.4549f, 0.151f);
    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){
    float gridN = sqrt(p_sites * 2.5f);
    float2 st = uv * gridN;
    int2 ipos = int2(floor(st));
    float2 fpos = forma_fract(st);
    float m_dist = 8.0f;
    float m_dist2 = 8.0f;
    int2 m_tile = ipos;
    for (int y = -1; y <= 1; y++){
      for (int x = -1; x <= 1; x++){
        int2 tile = ipos + int2(x, y);
        float h1 = hash2(tile.x, tile.y);
        float h2 = hash2(tile.x + 13, tile.y + 37);
        float ph = h1 * 6.283f + u_phase * 6.283f;
        /* Drift amplitude and edge gain both follow the JS path: 0.04f of the
           plate converts to cells by gridN, and the sharpness divides back out
           so both tabs draw the same seam width. The two used to disagree
           (0.08f cells and a gain of 4), which meant the tabs showed different
           diagrams. */
        float2 point = float2(float(x), float(y)) + float2(h1 + cos(u_t * 0.13f + ph) * 0.04f * gridN, h2 + sin(u_t * 0.13f + ph) * 0.04f * gridN);
        float dist = length(point - fpos);
        if (dist < m_dist){
          m_dist2 = m_dist;
          m_dist = dist;
          m_tile = tile;
        } else if (dist < m_dist2){
          m_dist2 = dist;
        }
      }
    }
    float edge = forma_clamp((m_dist2 - m_dist) * p_sharp / gridN, 0.0f, 1.0f);
    float val = hash2(m_tile.x, m_tile.y * 3);
    return ramp(val * 0.5f + 0.1f + (1.0f - edge) * p_edges * 0.6f);
  }

  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);
  }
};