Metaballs

DEFINITION

D(r) = b·e^(−a·r²)
F(p) = Σᵢ D(|p − cᵢ|)
draw the curves where F = threshold

NOTES

Blinn was rendering electron density maps for molecular models and needed surfaces that merged where atoms overlapped. Each centre contributes a Gaussian falloff, the field is simply their sum, and the surface is wherever that sum crosses a threshold — so two approaching blobs bulge toward one another and fuse without any special case for the join. He called them blobby models; everyone else calls them metaballs.

PROVENANCE

Origin

J. F. Blinn, "A Generalization of Algebraic Surface Drawing", ACM Transactions on Graphics 1(3), 1982

Standing

Public domain — a summation of exponentials

Contours

Extracted with marching squares (Lorensen & Cline, 1987, in its 2D form), which is why this field exports as real vector

Source

doi:10.1145/357306.357310

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. 63 · METABALLS — James Blinn, 1982
//   D(r) = b·e^(−a·r²)
//   F(p) = Σᵢ D(|p − cᵢ|)
//   draw the curves where F = threshold
// 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.5731;    // this plate's own grid phase, 0..1
// FORMA's FIELDS accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.3031, 0.11);
const vec3 u_pal_b = vec3(0.5, 0.3294, 0.1196);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_balls  = 6.0;         // centres · live 3 .. 9
const float p_tight  = 2.6;         // a — falloff · live 1.4 .. 5
const float p_levels = 3.0;         // contour levels · live 1 .. 5
const float p_iso    = 0.45;        // T — threshold · live 0.15 .. 0.9

/* 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;
}


/* Blinn's sum, centre for centre with the JS path: the same slow ellipses off
   the same clock and the same plate phase, the same 0.7 + 0.5*(i mod 3)/2
   weights, and the same aspect correction so a blob is round on a wide plate. */
float mbField(vec2 uv, float K, float ar){
  float f = 0.0;
  for (int i = 0; i < 9; i++){          // 9 is the centre slider's ceiling
    if (float(i) >= K) break;
    float ph = u_phase * 6.283 + float(i) * 6.283 / K;
    float cx = 0.5 + 0.30 * cos(u_t * 0.21 + ph * 1.7);
    float cy = 0.5 + 0.26 * sin(u_t * 0.17 + ph * 2.3);
    float cr = 0.7 + 0.5 * mod(float(i), 3.0) / 2.0;
    vec2 d = vec2(uv.x - cx, (uv.y - cy) * ar);
    f += cr * exp(-p_tight * 8.0 * dot(d, d));
  }
  return f;
}

vec3 plate(vec2 uv){
  /* Declared here rather than at file scope: this body is a composition
     operand, both operands may be this same specimen, and two globals of one
     name will not link. Only functions are namespaced. */
  vec3 ink = vec3(0.01569, 0.02353, 0.03922);   // #04060A
  float ar = u_res.y / u_res.x;
  float K = floor(p_balls + 0.5);
  /* The field is the same for every contour level, so it and its screen-space
     gradient are computed once, outside the loop — and outside any branch,
     which is where a derivative may honestly be taken. */
  float F = mbField(uv, K, ar);
  float grad = max(length(vec2(dFdx(F), dFdy(F))), 1e-7);
  /* The JS path strokes max(0.7, W/520) CSS pixels wide. Multiplied by the
     device ratio that is u_res.x/520 device pixels exactly, the ratio cancelling
     — the floor is all that has to be restated. */
  float lw = max(1.0, u_res.x / 520.0);
  float L = floor(p_levels + 0.5);
  vec3 col = ink;
  for (int l = 0; l < 5; l++){          // 5 is the level slider's ceiling
    if (float(l) >= L) break;
    float u = float(l) / max(1.0, L);
    /* The base threshold is the T the equation names; the levels above it are
       the same 0.42 fan the JS path spreads. */
    float iso = p_iso + 0.42 * u;
    float w = lw * (l == 0 ? 1.0 : 0.7);
    /* |F − iso| divided by how fast F changes per pixel is the distance to the
       contour measured in pixels, which is the only way a shader can hold a
       line to a constant width on a field whose slope varies by orders of
       magnitude between the centre of a blob and the space between two. */
    float cover = 1.0 - smoothstep(w * 0.5 - 0.5, w * 0.5 + 0.5, abs(F - iso) / grad);
    /* Laid down in the same order and at the same two alphas the JS path
       strokes them, so overlapping contours composite the same way. */
    col = mix(col, ramp(0.86 + 0.2 * u), (l == 0 ? 0.9 : 0.55) * cover);
  }
  return col;
}

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

AFTER EFFECTS · EXPRESSION

The same published mathematics as a Shape Layer path expression. Paste it onto a Path property; every constant is the published value plus a Slider Control named _tweak, so a bare paste already draws the figure and each slider moves one constant in its own units. Trim Paths is the comet.

// FORMA — PL. 63 · METABALLS — James Blinn, 1982
//   D(r) = b·e^(−a·r²)
//   F(p) = Σᵢ D(|p − cᵢ|)
//   draw the curves where F = threshold
// After Effects port — paste onto a Shape Layer's Path property
// (Contents › Shape › Path). Written from the published mathematics, not
// adapted from any code. Constants arrive at their published values; add a
// Slider Control (Effect › Expression Controls) named <k>_tweak and that
// constant moves in its own units, starting at 0 — the published figure.
// The plate's comet and its reveal are Trim Paths; the stroke colour is
// FORMA's FIELDS accent, #FFA83D. Animation runs on time.
// This plate draws 9 separate paths at its published constants:
// duplicate the group (Contents › Group) that many times and each copy draws
// its own part, read from its position in the layer. A Slider Control named
// "part" on the layer pins one instead.
// https://forma-gen.com/#plate=metaballs

// A missing slider reads 0, so a bare paste already draws the figure.
function forma_tweak(n){ try { return effect(n)("Slider"); } catch (e){ return 0; } }
var p_balls  = 6 + forma_tweak("balls_tweak");      // centres · live 3 .. 9
var p_tight  = 2.6 + forma_tweak("tight_tweak");    // a — falloff · live 1.4 .. 5
var p_levels = 3 + forma_tweak("levels_tweak");     // contour levels · live 1 .. 5
var p_iso    = 0.45 + forma_tweak("iso_tweak");     // T — threshold · live 0.15 .. 0.9

// The frame: the plate's W × H canvas is this comp, with the origin at the
// layer's anchor; canvas y already runs down, as After Effects' does.
var forma_W = thisComp.width, forma_H = thisComp.height, forma_t = time;
var forma_phase = 0.5731273817364126;   // this plate's own fixed phase, as the page has it
function forma_pt(x, y){ return [x - forma_W / 2, y - forma_H / 2]; }
function forma_partIndex(){
  try { return Math.round(effect("part")("Slider")); } catch (e){}
  try { return thisProperty.propertyGroup(3).propertyIndex - 1; } catch (e){ return 0; }
}
var forma_part = forma_partIndex();

// Blinn's metaballs as isolines: `balls` Gaussian sources drifting about the
// frame from this plate's phase, summed into a field, and `levels` contours
// of it from `iso` upward traced by marching squares on the page's own
// 96-wide lattice — the 16-case table, endpoints by linear interpolation
// along each cell edge, the two ambiguous saddles resolved the same way every
// time. A contour is a closed curve, so each is its own path: the cell
// segments are chained end to end (an edge shared by two cells interpolates
// to the same point in both), and the loops are ordered by level, then by
// length; one that leaves the frame is an open path. Their number changes
// as the blobs fuse and split — nine covers the published constants over
// their first seconds (nine at the first frame); a copy past the count draws
// nothing until a blob divides. The page's per-level tones are the
// stroke's.
// parts: 9
var nx = 96, ny = Math.max(24, Math.round(nx * forma_H / forma_W));
var K = Math.round(p_balls), cxs = [], cys = [], crs = [];
for (var i = 0; i < K; i++){
  var ph = forma_phase * 6.283 + i * 6.283 / K;
  cxs.push(0.5 + 0.30 * Math.cos(forma_t * 0.21 + ph * 1.7));
  cys.push(0.5 + 0.26 * Math.sin(forma_t * 0.17 + ph * 2.3));
  crs.push(0.7 + 0.5 * (i % 3) / 2);
}
var ar = forma_H / forma_W;
function forma_field(u, v){
  var f = 0;
  for (var b = 0; b < K; b++){
    var dx = u - cxs[b], dy = (v - cys[b]) * ar;
    f += crs[b] * Math.exp(-p_tight * 8 * (dx * dx + dy * dy));
  }
  return f;
}
var g = [];
for (var j = 0; j <= ny; j++) for (var i2 = 0; i2 <= nx; i2++) g.push(forma_field(i2 / nx, j / ny));
var L = Math.round(p_levels), loops = [];
function forma_chain(segs){
  // link segments end to end; endpoints keyed exactly. A contour that leaves
  // the frame is open, so those are chained first from their free end —
  // starting mid-way would cut one open contour into two — and what is left
  // closes on itself.
  var byStart = {}, byEnd = {};
  function keyOf(p){ return p[0].toFixed(6) + ',' + p[1].toFixed(6); }
  for (var s = 0; s < segs.length; s++){
    var k0 = keyOf(segs[s][0]), k1 = keyOf(segs[s][1]);
    (byStart[k0] = byStart[k0] || []).push(s);
    byEnd[k1] = 1;
  }
  var used = [], out = [];
  function follow(s0){
    var loop = [segs[s0][0]], cur = s0;
    while (cur >= 0 && !used[cur]){
      used[cur] = 1;
      var end = segs[cur][1];
      loop.push(end);
      var cands = byStart[keyOf(end)], next = -1;
      if (cands) for (var c = 0; c < cands.length; c++) if (!used[cands[c]]){ next = cands[c]; break; }
      cur = next;
    }
    var a = loop[0], b = loop[loop.length - 1];
    if (loop.length >= 3) out.push({ pts: loop, closed: a[0] === b[0] && a[1] === b[1] });
  }
  for (var s1 = 0; s1 < segs.length; s1++) if (!used[s1] && !byEnd[keyOf(segs[s1][0])]) follow(s1);
  for (var s2 = 0; s2 < segs.length; s2++) if (!used[s2]) follow(s2);
  return out;
}
for (var l = 0; l < L; l++){
  var iso = p_iso + 0.42 * (l / Math.max(1, L)), segs = [];
  for (var jj = 0; jj < ny; jj++) for (var ii = 0; ii < nx; ii++){
    var a = g[jj * (nx + 1) + ii], b2 = g[jj * (nx + 1) + ii + 1];
    var c2 = g[(jj + 1) * (nx + 1) + ii + 1], d2 = g[(jj + 1) * (nx + 1) + ii];
    var code = (a > iso ? 8 : 0) | (b2 > iso ? 4 : 0) | (c2 > iso ? 2 : 0) | (d2 > iso ? 1 : 0);
    if (code === 0 || code === 15) continue;
    var x0 = ii / nx * forma_W, x1 = (ii + 1) / nx * forma_W, y0 = jj / ny * forma_H, y1 = (jj + 1) / ny * forma_H;
    var mixT = x0 + (x1 - x0) * (iso - a) / (b2 - a), mixR = y0 + (y1 - y0) * (iso - b2) / (c2 - b2);
    var mixB = x0 + (x1 - x0) * (iso - d2) / (c2 - d2), mixL = y0 + (y1 - y0) * (iso - a) / (d2 - a);
    var top = [mixT, y0], right = [x1, mixR], bottom = [mixB, y1], left = [x0, mixL];
    // oriented so the inside (above iso) is always on the right hand of the
    // travel, which is what lets the chain close
    switch (code){
      case 1: segs.push([left, bottom]); break;   case 14: segs.push([bottom, left]); break;
      case 2: segs.push([bottom, right]); break;  case 13: segs.push([right, bottom]); break;
      case 3: segs.push([left, right]); break;    case 12: segs.push([right, left]); break;
      case 4: segs.push([right, top]); break;     case 11: segs.push([top, right]); break;
      case 6: segs.push([bottom, top]); break;    case 9: segs.push([top, bottom]); break;
      case 7: segs.push([left, top]); break;      case 8: segs.push([top, left]); break;
      case 5: segs.push([top, left], [bottom, right]); break;
      case 10: segs.push([left, bottom], [right, top]); break;
    }
  }
  var found = forma_chain(segs);
  found.sort(function (u, v){ return v.pts.length - u.pts.length; });
  for (var f = 0; f < found.length; f++) loops.push(found[f]);
}
var pts = [], pick = loops[Math.min(forma_part, loops.length - 1)] || { pts: [], closed: false };
var last = pick.closed ? pick.pts.length - 1 : pick.pts.length;
for (var m = 0; m < last; m++) pts.push(forma_pt(pick.pts[m][0], pick.pts[m][1]));
createPath(pts, [], [], pick.closed);