Mach Bands

DEFINITION

L(u) = j∕(N−1),  j = ⌊x(u)·N⌋,  x(u) triangular in u, period 1
DoG(d) = Gc(d) − Gs(d),  Gσ(d) = exp(−d²∕2σ²)∕(σ√2π)     Rodieck 1965
response(u) = (N−1)·(L * DoG)(u)
perceived(u) = L(u) + gain·response(u)

NOTES

A staircase of flat grey bars, each one shade brighter than the last, does not look like a staircase of flat grey bars: every bar reads slightly darker near the edge where a darker bar sits beside it, and slightly lighter near the edge where a lighter one does, so a thin bright seam runs along the light side of every step and a thin dark seam along the shadowed side. Photograph the same staircase and the seams are gone, because nothing in the light changes there — only two flat greys sit side by side, and a photometer reports exactly the two numbers printed on the bars. Mach reported the seams to the Vienna Academy in 1865, before anyone had a working model of why a retina would invent contrast that is not in the light at all. Rodieck supplied the model in 1965: a retinal ganglion cell answers for a small patch of the image with one number, built from two overlapping pools, a narrow one that excites and a wider, weaker one around it that inhibits — the difference of two Gaussian sensitivity profiles, normalised so a uniformly lit patch of any brightness reports nothing, because only a change in the light moves the two pools out of balance. Convolve the staircase against that kernel in closed form and the seams stop being a description: a positive bump appears on the light side of every riser, a negative dip on the dark side, and both fall to exactly zero across any plateau wide enough that no edge reaches it, which a smooth ramp of the same overall range never has, since a bandpass kernel answers only to curvature and a straight ramp has none. The top register is the staircase as a photometer would read it: flat bars, hard edges, nothing invented. The bottom draws two curves over the same bars — a faint flat line at each bars true grey, and riding past it at every edge, the bright curve the receptor field actually reports, by an amount the ratio and gain sliders set. More bars in the same strip means a smaller true jump at each one, so the response is measured in units of one risers own height before either slider touches it, otherwise a finer staircase would quietly bury its own exhibit by shrinking the very thing being shown. The whole strip drifts sideways without end, one infinite corrugation sliding behind a fixed window, so the bands on screen now are never quite the ones a moment ago.

PROVENANCE

Origin

E. Mach, "Ueber die Wirkung der raeumlichen Vertheilung des Lichtreizes auf die Netzhaut", Sitzungsberichte der mathematisch-naturwissenschaftlichen Classe der kaiserlichen Akademie der Wissenschaften (Wien) 52, 1865, 303-322. Pre-DOI — no resolvable identifier exists for a paper this old, and none is offered. Verified against the German Wikipedia entry for Machsche Streifen and independently against michaelbach.de/ot/lum-MachBands, which both give the same title, volume and page range

The receptive-field model

R. W. Rodieck, "Quantitative analysis of cat retinal ganglion cell response to visual stimuli", Vision Research 5(11), 1965, 583-601 (the src link above). This is the paper that proposed a ganglion cell receptive field is well described by a difference of two Gaussian sensitivity profiles — the DoG this plate runs, unmodified in shape, just built at whatever size the ratio slider names rather than fitted to one cats retina. Verified on Crossref against title, author, year and venue before the doi went in

What is a model, not a measurement

Mach measured a real perceptual effect in 1865. Rodieck proposed a mathematical account of the retinal mechanism a century later, fitted to specific cells, not to Machs own plates. Running Rodiecks kernel over Machs staircase, as this plate does, is a reconstruction — it shows what the published account of lateral inhibition predicts a bar of that kind would provoke, not a photograph of any one nerve firing. Every DoG plate on the market makes this same substitution; the honest thing is to say so rather than let the bottom register read as a measurement

Standing

Public domain. The staircase is quantised directly from its own definition, the kernel is Rodiecks published difference of two Gaussians evaluated in closed form, and the two are combined by an ordinary convolution written out as a discrete sum. Nothing is taken from any vision-science or imaging codebase

Constants

steps (N) sets how many flat bars fill the strip; the strip is a mirrored (triangular) staircase, climbing across the first half of its own period and descending back across the second, so every riser is the same height, including the one where the drifting domain wraps — a plain sawtooth staircase would have one riser N-1 times taller than the rest exactly at the seam, which would dwarf every other band and is not what is being shown. ratio is the surround sigma divided by the centre sigma in Rodiecks kernel; both sigmas are set internally as a fixed fraction (0.22) of the local half-step width, so the kernel scales itself to the staircase rather than needing a fourth slider, and the relative geometry (how many step-widths the surround spans) stays the same at every step count by construction. gain multiplies the response after it has already been divided by one risers own height (N-1 lots of the raw convolution) — without that normalisation the raw response shrinks roughly as 1/(N-1) simply because each riser gets smaller as more of them fill the same strip, and a finer staircase would look like a weaker exhibit for no reason connected to lateral inhibition at all. Internal (unexposed) constants: TAPS = 64 taps for the discrete convolution and EXT_K = 3.2 surround-sigmas of half-width, both fixed rather than dialled, chosen from the quadrature check below

Checked, not just plotted

Quadrature: the discrete convolution at TAPS = 64 was checked against a 2000-tap reference at the worst corner in the declared box (steps = 9, ratio = 3.5, where the centre Gaussian is narrowest relative to the tap spacing) and differs by 2.1% in the responses own peak — a shape number, not a visibility one. Scale invariance: normalising by (N-1) makes the response nearly step-count-independent, as the underlying mathematics of an isolated edge says it should be — measured at ratio = 2.0, the peak response in riser-height units runs 0.1647 at steps = 3 to 0.1729 at steps = 9, a 5% spread across the whole range rather than the roughly 3x spread the un-normalised response shows. Liveness: 960 random all-constants tuples (8 seeds of 120, steps drawn from 3..9, ratio from 1.4..3.5, gain from 1.0..2.2) all cleared 0.09 of the 2.0-unit plotted axis, worst case 0.098 — no tuple came back dark. One-at-a-time, holding the other two dials at their published defaults (steps = 6, ratio = 2.0, gain = 1.5): steps alone runs the overshoot from 0.2470 to 0.2586 across its full range (by design, nearly flat — the normalisation is doing its job); ratio alone runs it from 0.1301 to 0.4930; gain alone from 0.1683 to 0.3702. Every dial has a live, graduated range on its own and in combination

Source

doi:10.1016/0042-6989(65)90033-7

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. 130 · MACH BANDS — Ernst Mach, 1865 · Robert W. Rodieck, 1965
//   L(u) = j∕(N−1),  j = ⌊x(u)·N⌋,  x(u) triangular in u, period 1
//   DoG(d) = Gc(d) − Gs(d),  Gσ(d) = exp(−d²∕2σ²)∕(σ√2π)     Rodieck 1965
//   response(u) = (N−1)·(L * DoG)(u)
//   perceived(u) = L(u) + gain·response(u)
// 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.8096;    // this plate's own grid phase, 0..1
// FORMA's COLOUR accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.2002, 0.3896);
const vec3 u_pal_b = vec3(0.5, 0.2176, 0.4235);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);

const float p_steps = 6.0;         // N — steps in the staircase · live 3 .. 9
const float p_ratio = 2.0;         // surround / centre sigma ratio · live 1.4 .. 3.5
const float p_gain  = 1.5;         // response gain (per riser height) · live 1 .. 2.2

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


/* Same construction as the JS path: a mirrored (triangular) luminance
   staircase convolved against a difference-of-Gaussians receptive field,
   Rodieck 1965. The JS path caches the response in a table keyed on step
   count and ratio, because it does not depend on time at all; a fragment
   shader has nowhere to put that table, so every fragment in the curve
   register runs the same 64-tap convolution the table was built from — the
   numbers are identical, the two paths are one picture sampled twice. */

/* The staircase, mirrored the same way the JS path's triLevel is. */
float mach_triLevel(float u, float n){
  float uu = fract(u);
  float x = uu < 0.5 ? uu * 2.0 : (1.0 - uu) * 2.0;
  float j = clamp(floor(x * n), 0.0, n - 1.0);
  return j / (n - 1.0);
}

/* Rodieck 1965, the same two normalised Gaussians as the JS path. */
float mach_dog(float d, float sc, float ss){
  float gc = exp(-0.5 * (d * d) / (sc * sc)) / (sc * 2.5066282746310002);
  float gs = exp(-0.5 * (d * d) / (ss * ss)) / (ss * 2.5066282746310002);
  return gc - gs;
}

/* Constant loop bound — 64, matching TAPS in the JS table build exactly, so
   both paths integrate the same quadrature. It does not depend on any
   parameter, so there is no live ceiling to reach for and nothing to break
   on. Normalised by (n-1), one riser's own height, for the reason the JS
   path's Constants note gives: without it a finer staircase would quietly
   shrink the very exhibit being shown. */
float mach_response(float u, float n, float sc, float ss){
  float ext = 3.2 * ss;
  float dd = 2.0 * ext / 64.0;
  float sum = 0.0;
  for (int k = 0; k < 64; k++){
    float d = -ext + (float(k) + 0.5) * dd;
    float uu = u - d;
    sum += mach_triLevel(uu, n) * mach_dog(d, sc, ss) * dd;
  }
  return sum * (n - 1.0);
}

float mach_mark(float e0, float e1, float x){
  float s = clamp((x - e0) / (e1 - e0), 0.0, 1.0);
  return 1.0 - s * s * (3.0 - 2.0 * s);
}

vec3 plate(vec2 uv){
  const float PHYS_T = 0.07, PHYS_B = 0.42, CURVE_T = 0.52, CURVE_B = 0.93;
  const float VMIN = -0.5, VMAX = 1.5;
  const float CFRAC = 0.22;
  float asp = u_res.y / u_res.x;

  float n = floor(p_steps + 0.5);
  float drift = cycle(u_t, 20.0);
  float xi = fract(uv.x + drift);
  float L = mach_triLevel(xi, n);

  /* The staircase as a photometer would read it — cheap, so it costs
     nothing to evaluate before branching. */
  if (uv.y > PHYS_T && uv.y < PHYS_B) return vec3(L);

  vec3 ink = vec3(4.0, 6.0, 10.0) / 255.0;
  if (uv.y > CURVE_T && uv.y < CURVE_B){
    /* Both sigmas are a fixed fraction of the local half-step width, so the
       kernel scales itself to the staircase — matching the JS path exactly. */
    float stepW = 0.5 / n, sc = CFRAC * stepW, ss = p_ratio * sc;
    float P = clamp(L + p_gain * mach_response(xi, n, sc, ss), VMIN, VMAX);
    float vL = CURVE_B - (L - VMIN) / (VMAX - VMIN) * (CURVE_B - CURVE_T);
    float vP = CURVE_B - (P - VMIN) / (VMAX - VMIN) * (CURVE_B - CURVE_T);
    const float LW = 0.007;
    vec3 col = ink;
    float dL = abs(uv.y - vL) * asp;
    col = mix(col, vec3(0.30, 0.32, 0.38), mach_mark(LW * 0.5, LW, dL) * 0.5);
    /* The order accent, in the ramp's bright lobe — a sparse mark on a dark
       ground, exactly as the JS path keeps it. */
    float dP = abs(uv.y - vP) * asp;
    col = mix(col, ramp(0.92), mach_mark(LW * 0.5, LW, dP));
    return col;
  }
  return ink;
}

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. 130 · MACH BANDS — Ernst Mach, 1865 · Robert W. Rodieck, 1965
//   L(u) = j∕(N−1),  j = ⌊x(u)·N⌋,  x(u) triangular in u, period 1
//   DoG(d) = Gc(d) − Gs(d),  Gσ(d) = exp(−d²∕2σ²)∕(σ√2π)     Rodieck 1965
//   response(u) = (N−1)·(L * DoG)(u)
//   perceived(u) = L(u) + gain·response(u)
// 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_mach : ImageComputationKernel<ePixelWise>
{
  Image<eWrite> dst;

param:
  float u_t;             // seconds; 0 is the still frame
  float p_steps; // N — steps in the staircase · live 3 .. 9
  float p_ratio; // surround / centre sigma ratio · live 1.4 .. 3.5
  float p_gain;  // response gain (per riser height) · live 1 .. 2.2

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_steps, "p_steps", 6.0f);
    defineParam(p_ratio, "p_ratio", 2.0f);
    defineParam(p_gain, "p_gain", 1.5f);
  }

  void init(){
    u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
    u_phase = 0.8096f;    // this plate's own grid phase, 0..1
    // FORMA's COLOUR accent as cosine-gradient coefficients
    u_pal_a = float3(0.46f, 0.2002f, 0.3896f);
    u_pal_b = float3(0.5f, 0.2176f, 0.4235f);
    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;
  }


  /* Same construction as the JS path: a mirrored (triangular) luminance
     staircase convolved against a difference-of-Gaussians receptive field,
     Rodieck 1965.f The JS path caches the response in a table keyed on step
     count and ratio, because it does not depend on time at all; a fragment
     shader has nowhere to put that table, so every fragment in the curve
     register runs the same 64-tap convolution the table was built from — the
     numbers are identical, the two paths are one picture sampled twice. */

  /* The staircase, mirrored the same way the JS path's triLevel is. */
  float mach_triLevel(float u, float n){
    float uu = forma_fract(u);
    float x = uu < 0.5f ? uu * 2.0f : (1.0f - uu) * 2.0f;
    float j = forma_clamp(floor(x * n), 0.0f, n - 1.0f);
    return j / (n - 1.0f);
  }

  /* Rodieck 1965, the same two normalised Gaussians as the JS path. */
  float mach_dog(float d, float sc, float ss){
    float gc = exp(-0.5f * (d * d) / (sc * sc)) / (sc * 2.5066282746310002f);
    float gs = exp(-0.5f * (d * d) / (ss * ss)) / (ss * 2.5066282746310002f);
    return gc - gs;
  }

  /* Constant loop bound — 64, matching TAPS in the JS table build exactly, so
     both paths integrate the same quadrature. It does not depend on any
     parameter, so there is no live ceiling to reach for and nothing to break
     on. Normalised by (n-1), one riser's own height, for the reason the JS
     path's Constants note gives: without it a finer staircase would quietly
     shrink the very exhibit being shown. */
  float mach_response(float u, float n, float sc, float ss){
    float ext = 3.2f * ss;
    float dd = 2.0f * ext / 64.0f;
    float sum = 0.0f;
    for (int k = 0; k < 64; k++){
      float d = -ext + (float(k) + 0.5f) * dd;
      float uu = u - d;
      sum += mach_triLevel(uu, n) * mach_dog(d, sc, ss) * dd;
    }
    return sum * (n - 1.0f);
  }

  float mach_mark(float e0, float e1, float x){
    float s = forma_clamp((x - e0) / (e1 - e0), 0.0f, 1.0f);
    return 1.0f - s * s * (3.0f - 2.0f * s);
  }

  float3 plate(float2 uv){
    float3 forma_r = float3(0.0f, 0.0f, 0.0f);
    for (int forma_once = 0; forma_once < 1; forma_once++){
      const float PHYS_T = 0.07f;
      const float PHYS_B = 0.42f;
      const float CURVE_T = 0.52f;
      const float CURVE_B = 0.93f;
      const float VMIN = -0.5f;
      const float VMAX = 1.5f;
      const float CFRAC = 0.22f;
      float asp = u_res.y / u_res.x;

      float n = floor(p_steps + 0.5f);
      float drift = cycle(u_t, 20.0f);
      float xi = forma_fract(uv.x + drift);
      float L = mach_triLevel(xi, n);

      /* The staircase as a photometer would read it — cheap, so it costs
         nothing to evaluate before branching. */
      if (uv.y > PHYS_T && uv.y < PHYS_B) { forma_r = float3(L); break; }

      float3 ink = float3(4.0f, 6.0f, 10.0f) / 255.0f;
      if (uv.y > CURVE_T && uv.y < CURVE_B){
        /* Both sigmas are a fixed fraction of the forma_local half-step width, so the
           kernel scales itself to the staircase — matching the JS path exactly. */
        float stepW = 0.5f / n;
        float sc = CFRAC * stepW;
        float ss = p_ratio * sc;
        float P = forma_clamp(L + p_gain * mach_response(xi, n, sc, ss), VMIN, VMAX);
        float vL = CURVE_B - (L - VMIN) / (VMAX - VMIN) * (CURVE_B - CURVE_T);
        float vP = CURVE_B - (P - VMIN) / (VMAX - VMIN) * (CURVE_B - CURVE_T);
        const float LW = 0.007f;
        float3 col = ink;
        float dL = fabs(uv.y - vL) * asp;
        col = lerp(col, float3(0.30f, 0.32f, 0.38f), mach_mark(LW * 0.5f, LW, dL) * 0.5f);
        /* The order accent, in the ramp's bright lobe — a sparse mark on a dark
           ground, exactly as the JS path keeps it. */
        float dP = fabs(uv.y - vP) * asp;
        col = lerp(col, ramp(0.92f), mach_mark(LW * 0.5f, LW, dP));
        { forma_r = col; break; }
      }
      { forma_r = ink; 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);
  }
};