Curve-Shortening Flow

DEFINITION

∂γ/∂t = κ N = ∂²γ/∂s²
dL/dt = −∮ κ² ds ≤ 0
dA/dt = −∮ κ ds = −2π  (any simple closed curve)
4πA/L² → 1

NOTES

Give every point of a closed loop a velocity equal to its own curvature, aimed at the centre of the circle that best fits the loop there. Written against arclength the rule is nothing but the heat equation applied to position, so the curve diffuses: bulges flatten, dents fill, and a wiggle decays fastest exactly where it is sharpest. Gage and Hamilton proved in 1986 that a convex curve treated this way stays convex, grows circular, and shrinks away in finite time to what they called a round point. Grayson closed the gap a year later, and it is the harder half: any embedded curve at all, meaning any closed loop that does not cross itself, first becomes convex — after which the earlier theorem finishes the job. So every seeded squiggle on this plate ends the same way, and the ending is not a modelling choice, it is a theorem. The seed here is a polar graph with a strictly positive radius. That is deliberate rather than convenient: a radial graph cannot cross itself, so the starting curve is embedded by construction and sits inside the hypothesis Grayson actually proved rather than merely near it. The bright arcs are the vertices where the polygon turns the wrong way, and they are the only mark on the plate that means the curve is not yet convex. Watch them go out — that is Grayson, live, and by the theorem they never come back. Two numbers are worth following. Area falls at exactly 2π per unit time whatever shape the curve is in, because the rate of change of area is the integral of curvature around the loop and that integral is the total turning of a simple closed curve. This plate uses that as its own clock, advancing a fixed slice of area every frame, and banks a faint outline at equal decrements of area — which by the same identity is equal intervals of flow time, making the nested family a stroboscopic record at a fixed exposure interval rather than a stack of arbitrary snapshots. The other number is the isoperimetric ratio 4πA/L², which is 1 for a circle and less for everything else. From the default seed it starts at 0.736, passes 0.972 at the instant the last concave arc disappears, and reaches 0.9995 with the curve at a third of the length it began with. Nothing is rescaled to keep it in frame: the shrink drawn is the shrink the equation produces, because vanishing to a point is the statement.

PROVENANCE

Origin

M. Gage and R. S. Hamilton, "The heat equation shrinking convex plane curves", Journal of Differential Geometry 23(1), 1986, doi:10.4310/jdg/1214439902 — verified on Crossref against title, both authors, year, journal, volume and issue. This is the convex case: a convex embedded curve stays convex, becomes asymptotically circular, and contracts to a point in finite time

Any embedded curve

M. A. Grayson, "The heat equation shrinks embedded plane curves to round points", Journal of Differential Geometry 26(2), 1987, doi:10.4310/jdg/1214441371 — verified on Crossref the same way. Grayson supplies the missing half: an embedded curve becomes convex in finite time without ever developing a self-intersection, so the 1986 result then applies to it. Both papers sit in the same journal a year apart, and the pair is the reason this plate can seed any squiggle at all rather than only convex ones. src points at Grayson because that is the theorem the plate exhibits; the 1986 paper is the one it ends with

What the registry does and does not carry

Neither deposited record carries a page range — the DOIs resolve to the articles, and the Crossref metadata stops at volume and issue. So no page range is quoted above, on this atlas own rule that a citation which is nearly right is worse than one that stops where the evidence does. The commonly given ranges are easy to find and would have looked more authoritative; they would not have been checked

Standing

Public domain — a partial differential equation and two theorems about it, implemented here from the published statements. Nothing was adapted from any existing implementation

Checked, not just plotted

The identity dA/dt = −2π is independent of everything the scheme does, so it is the right thing to measure against. On a circle the discrete flow is exact to the last digit: the vertex Laplacian of an n-gon inscribed in a circle of radius r evaluates to 1/r exactly, so the polygon obeys dA/dt = −n sin(2π/n), which is −6.28157 at 160 nodes against −2π = −6.28319, a gap of 0.026% that is the polygon and not the method. Measured over a full run from a seeded squiggle at 160 nodes with the redistribution switched off, dA/dt comes out −6.27895, which is 0.067% from −2π and 0.042% from the exact n-gon figure. Radius against the closed form is the second check: r(t) = sqrt(r0² − 2t) is reproduced to 0.098% at 160 nodes over a shrink to three tenths of the starting radius, and the residual halves as the node count doubles, which is the forward-Euler truncation behaving as it should

The time step is measured, not assumed

The vertex Laplacian on a uniform polygon has its most negative eigenvalue at −4/h², the alternating zigzag mode, so forward Euler multiplies that mode by 1 − 4λ each step where λ = dt/h², and the scheme is stable exactly while λ ≤ 1/2. That bound was measured rather than taken on trust: seeding a zigzag of one part in a thousand onto a 160-gon and running 400 steps, the amplitude comes back multiplied by 4e−12 at λ = 0.2, by 3e−12 at 0.45, by exactly 1.000 at λ = 0.500 — the marginal case, where the factor is −1 — and by 4.9 at 0.501 and 1.1e3 at 0.55. This plate runs at λ = 0.25, half the bound, with dt recomputed from the shortest edge every step because the edges shrink with the curve

Why the nodes are redistributed

Under this flow arclength decays at rate κ², so nodes crowd into the sharp arcs and the flat arcs starve. Both halves of that are fatal and the plate would not announce either: the time step is bounded by the SHORTEST edge, so crowding buys nothing but arithmetic, while the stretched edges quietly stop resolving the curve. Measured on a modes 8, wobble 0.75 seed over 6,000 steps, left alone the spacing ratio reaches 20.6 and the run covers 15.9 times less flow time than a redistributed one. The fix is to resample to uniform arclength — but a resampled polygon is inscribed in the one it replaces, so it can only lose area, and doing it every step biases dA/dt by 0.83%. Firing it only when the spacing ratio passes 1.25 holds the spacing inside 26% for the price of 0.22%, at seven resamples in twelve hundred steps. Those two numbers are the whole trade and both are on the record

Embedded, and it stays embedded

Grayson proves an embedded curve never touches itself under the flow, and a discrete scheme is free to violate that even when the continuum does not — a neck closing between two nodes would simply pass through. Checked by testing every pair of non-adjacent edges for crossing, sampled through whole runs at five corners of the constant box including the two-lobed dumbbell that necks hardest: zero crossings anywhere, with the closest approach between nodes more than six apart never falling below 4.4 mean edge lengths

What is drawn, and what is not

The true shrink, at true scale. Normalising the area would hold the curve at a constant size and make the plate legible for longer, and it would also delete the theorem — Gage and Hamilton prove it vanishes, so the vanishing is the picture. The nested outlines are what makes that readable: each is banked at an equal decrement of area, which is an equal interval of flow time, so the family thickens toward the end exactly as the self-similar final phase demands. The seed outline is kept faintly for the whole run, because within a few seconds there is nothing else left to say where the flow began

Constants

nodes is resolution and cost: the same seed curve sampled more finely, at steps growing as nodes², and locked from REGENERATE for that reason — measured mean 4.7 steps a frame at 90 nodes against 25.7 at 260, for cycles of 258 and 269 frames. rate is a percentage of the seed area consumed per frame rather than a step count, which is only possible because dA/dt is constant: it holds the cycle at the same length in frames across the whole nodes range, where a fixed step budget would have varied by eight to one. Measured at the extremes it runs 498 frames at 0.15 and 91 at 0.9, and both are worth having. modes and wobble shape the seed; wobble is the sum of the harmonic amplitudes, so the floor of 0.22 is where a seed is still reliably non-convex and the ceiling of 0.75 leaves a quarter of the radius, keeping the polar graph positive and the curve embedded. Swept the way REGENERATE jitters — all five constants at once, 40 tuples across 5 seeds — all 200 runs closed a cycle, the longest taking 403 frames, the latest convexity landing at frame 109, and the most expensive single frame costing 66 steps and 0.79 ms at 900 px

Source

doi:10.4310/jdg/1214441371

HOUDINI · VEX

The same published mathematics as a Detail Wrangle body. Paste it into a Wrangle with Run Over set to Detail; every constant is the published value plus a tweak channel, so Create Spare Parameters gives a slider that starts where the paper does.

// FORMA — PL. 140 · CURVE-SHORTENING FLOW — Michael Gage, Richard Hamilton, 1986 · Matthew Grayson, 1987
//   ∂γ/∂t = κ N = ∂²γ/∂s²
//   dL/dt = −∮ κ² ds ≤ 0
//   dA/dt = −∮ κ ds = −2π  (any simple closed curve)
//   4πA/L² → 1
// Paste into a Detail Wrangle (Run Over: Detail), no inputs needed.
// Written from the published mathematics, not adapted from any code.
// Constants arrive at their published values. Press the node's Create
// Spare Parameters button and every tweak becomes a slider — starting
// at 0, the published figure, and moving in the constant's own units.
// https://forma-gen.com/#plate=curvatureflow

float p_nodes  = 160 + chf('nodes_tweak');      // nodes on the curve · live 90 .. 260
float p_modes  = 5 + chf('modes_tweak');        // harmonics in the seed · live 2 .. 8
float p_wobble = 0.5 + chf('wobble_tweak');     // seed amplitude (Σ|aₖ|) · live 0.22 .. 0.75
float p_ghosts = 16 + chf('ghosts_tweak');      // outlines banked · live 6 .. 28

// The plate's own colour: FORMA's CURVES accent as a cosine ramp,
// brightest near t = 0 and t = 1, near-black around t = 0.5.
vector forma_ramp(float t){
  return set(
    0.11 + 0.1196 * cos(6.28318530718 * (t + 0)),
    0.46 + 0.5 * cos(6.28318530718 * (t + 0.05)),
    0.2453 + 0.2667 * cos(6.28318530718 * (t + 0.1)));
}

// Curve-shortening flow cooked to its own stopping test in a single pass.
// The page advances the flow one slice of area per frame against a frame
// clock and shows the curve moving; a cook has no clock, so this runs the
// same scheme to the same two exits — the isoperimetric ratio reaching the
// round-point threshold, or the area falling past the floor — and emits
// what the page has accumulated on screen by then: the seed outline, the
// nested outlines banked at equal decrements of area, and the curve it
// stopped on. Equal area is equal flow time, because dA/dt = -2*PI for any
// simple closed curve, so the nest is a stroboscopic record at a fixed
// interval rather than a stack of arbitrary snapshots.
//
// The scheme is the one the plate runs and the numbers are the measured
// ones. Velocity is the second derivative of position with respect to
// arclength on the three-point stencil for unequal spacing, which IS the
// curvature vector kappa*N — it points at the centre of the osculating
// circle, so no normal and no sign has to be derived separately. The step
// is 0.25*hmin^2 against a stability bound measured at 0.5*hmin^2. The
// nodes are redistributed to uniform arclength whenever the spacing ratio
// passes 1.25, because arclength decays at rate kappa^2 and left alone the
// nodes crowd into the sharp arcs while the flat arcs starve.
//
// Deterministic: the seed harmonics are random(counted seed) on the seed
// the page uses for its first epoch, so every cook gives the same curve.
// The page also turns the seed by the per-plate PHASE, which exists to
// keep a grid of plates out of step with each other and has nothing to do
// in a scene, so it is left out here. Domain is the 560-unit reference
// canvas the corpus shares, and canvas y is negated onto the Houdini up
// axis.
// waived: rate — it paces the plate against a frame clock and a cook has
// no clock; nodes, modes, wobble and ghosts all reach the body

// One closed outline. Every input arrives as a parameter — a VEX function
// cannot read a snippet-scope variable.
function void forma_outline(float xs[]; float ys[]; int n; float w; float hh; vector col; float a){
    int prim = addprim(0, "polyline");
    for (int i = 0; i <= n; i++){
        int k = (i == n) ? 0 : i;              // last vertex repeats the first
        int pt = addpoint(0, set(xs[k] - w / 2.0, (hh / 2.0) - ys[k], 0.0));
        setpointattrib(0, "Cd", pt, col);
        setpointattrib(0, "Alpha", pt, a);
        addvertex(0, prim, pt);
    }
}

// The seed, with the vertices that turn the wrong way marked. Those arcs
// are the only thing on the plate that means the curve is not yet convex,
// and Grayson is the theorem that says they all go and none come back.
function void forma_seedline(float xs[]; float ys[]; int mark[]; int n; float w; float hh;
                    vector col; vector hi; float a){
    int prim = addprim(0, "polyline");
    for (int i = 0; i <= n; i++){
        int k = (i == n) ? 0 : i;
        int pt = addpoint(0, set(xs[k] - w / 2.0, (hh / 2.0) - ys[k], 0.0));
        setpointattrib(0, "Cd", pt, mark[k] ? hi : col);
        setpointattrib(0, "Alpha", pt, mark[k] ? 0.7 : a);
        addvertex(0, prim, pt);
    }
}

float W = 560.0, H = 560.0;
float TAU = 2.0 * M_PI;
int   forma_n     = int(rint(p_nodes));
int   forma_modes = int(rint(p_modes));
int   forma_g     = int(rint(p_ghosts));
int   forma_max   = 9000;      // the cook ceiling; see the note at the exit
float forma_cfl   = 0.25;      // dt = cfl*hmin^2, half the measured bound
float forma_space = 1.25;      // redistribute past this spacing ratio
float forma_round = 0.9995;    // the 4*PI*A/L^2 that counts as a round point
float forma_floor = 0.06;      // or stop here, as a fraction of the seed area

// ---- the seed ----------------------------------------------------------
// A polar graph r(theta) on top of a circle. The amplitudes are normalised
// to sum to wobble, so r stays at or above R*(1 - wobble) and the ceiling of
// 0.75 still leaves a quarter of the radius. A radial graph with a positive
// radius cannot cross itself, so the seed is EMBEDDED BY CONSTRUCTION —
// which is the hypothesis Grayson proves the theorem under, rather than
// something merely hoped for.
float amp[], ph[];
resize(amp, forma_modes);
resize(ph, forma_modes);
int rc = 1987;                             // the page seed, counted upward
float asum = 0.0;
for (int m = 0; m < forma_modes; m++){
    amp[m] = 0.35 + 0.65 * random(rc);  rc++;
    asum += amp[m];
}
for (int m = 0; m < forma_modes; m++) amp[m] *= p_wobble / asum;
for (int m = 0; m < forma_modes; m++){ ph[m] = random(rc) * TAU;  rc++; }

float x[], y[], h[], vx[], vy[], ox[], oy[];
resize(x, forma_n);   resize(y, forma_n);   resize(h, forma_n);
resize(vx, forma_n);  resize(vy, forma_n);
resize(ox, forma_n);  resize(oy, forma_n);

// Sampled at unit radius, then fitted to the frame from its own measured
// extent: the harmonics only reach the full wobble where their phases happen
// to agree, so a radius sized for the worst case measured 59% to 71% of the
// frame across the wobble range instead of the 88% it was sized for. Scale
// is free — the time step goes as the square of the size and so does the
// area budget — so fitting costs the mathematics nothing.
float minx = 1e9, maxx = -1e9, miny = 1e9, maxy = -1e9;
for (int i = 0; i < forma_n; i++){
    float th = TAU * float(i) / float(forma_n);
    float r = 1.0;
    for (int m = 0; m < forma_modes; m++) r += amp[m] * cos(float(m + 2) * th + ph[m]);
    float ux = r * cos(th), uy = r * sin(th);
    ox[i] = ux;  oy[i] = uy;
    minx = min(minx, ux);  maxx = max(maxx, ux);
    miny = min(miny, uy);  maxy = max(maxy, uy);
}
float sc = min(0.86 * W / (maxx - minx), 0.86 * H / (maxy - miny));
float fx = W / 2.0 - (minx + maxx) * 0.5 * sc;
float fy = H / 2.0 - (miny + maxy) * 0.5 * sc;
for (int i = 0; i < forma_n; i++){ x[i] = fx + ox[i] * sc;  y[i] = fy + oy[i] * sc; }

// ---- the flow ----------------------------------------------------------
// The seed measurements are taken on the first pass through the loop below
// rather than in a block of their own, so there is exactly one copy of the
// redistribution in this file. The polar graph is sampled at equal angle and
// therefore does not arrive at equal arclength, which is why step zero
// redistributes whatever its spacing ratio happens to be.
float A0 = 0.0, nextBank = 0.0;
int   banked = 0;
int   steps = 0;
int   mark[];
resize(mark, forma_n);

for (steps = 0; steps <= forma_max; steps++){
    float L = 0.0, hmin = 1e9, hmax = 0.0;
    for (int i = 0; i < forma_n; i++){
        int j = (i + 1 == forma_n) ? 0 : i + 1;
        h[i] = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
        L += h[i];
        hmin = min(hmin, h[i]);
        hmax = max(hmax, h[i]);
    }

    // Redistribute to uniform arclength. This is a purely tangential motion —
    // it slides the vertices along the curve, not the curve — but the polygon
    // it produces is inscribed in the one it replaces, so it loses a little
    // area every time it runs. Measured on the page, firing it every step
    // biases dA/dt by 0.83% and firing it on this spacing test costs 0.22%,
    // which is why it is a test and not a schedule.
    if ((steps == 0 || hmax > forma_space * hmin) && L > 0.0){
        for (int i = 0; i < forma_n; i++){ ox[i] = x[i];  oy[i] = y[i]; }
        float seg = L / float(forma_n);
        int e = 0;
        float acc = 0.0;
        for (int i = 1; i < forma_n; i++){
            float target = float(i) * seg;
            while (e < forma_n - 1 && acc + h[e] < target){ acc += h[e];  e++; }
            float u = (h[e] > 0.0) ? (target - acc) / h[e] : 0.0;
            int j = (e + 1 == forma_n) ? 0 : e + 1;
            x[i] = ox[e] + (ox[j] - ox[e]) * u;
            y[i] = oy[e] + (oy[j] - oy[e]) * u;
        }
        L = 0.0;  hmin = 1e9;  hmax = 0.0;
        for (int i = 0; i < forma_n; i++){
            int j = (i + 1 == forma_n) ? 0 : i + 1;
            h[i] = sqrt((x[j] - x[i]) * (x[j] - x[i]) + (y[j] - y[i]) * (y[j] - y[i]));
            L += h[i];
            hmin = min(hmin, h[i]);
            hmax = max(hmax, h[i]);
        }
    }

    // Twice the signed area, taken about the first vertex rather than about
    // the origin. The two are algebraically the same; on this side they are
    // not, because VEX floats carry seven digits and a region of a hundred
    // thousand square units sitting three hundred units out spends two of
    // them on cancellation — the lesson the Lloyd port paid for. The sign
    // carries the winding, which the concave test needs, so it is kept.
    float b2 = 0.0;
    for (int i = 0; i < forma_n; i++){
        int j = (i + 1 == forma_n) ? 0 : i + 1;
        b2 += (x[i] - x[0]) * (y[j] - y[0]) - (x[j] - x[0]) * (y[i] - y[0]);
    }
    float A = abs(b2) * 0.5;
    float iso = 4.0 * M_PI * A / (L * L);

    // Step zero is the seed, now evened out: its area is what every later
    // measurement is a fraction of, and its concave vertices are the arcs
    // Grayson says will all go and none come back. The page marks them in
    // HILITE on the live curve; a still can only show them where they start.
    if (steps == 0){
        A0 = A;
        nextBank = A0 - A0 / float(forma_g + 1);
        float wind = (b2 >= 0.0) ? 1.0 : -1.0;
        for (int i = 0; i < forma_n; i++){
            int a = (i == 0) ? forma_n - 1 : i - 1;
            int b = (i + 1 == forma_n) ? 0 : i + 1;
            float ux = x[i] - x[a], uy = y[i] - y[a];
            float wx = x[b] - x[i], wy = y[b] - y[i];
            mark[i] = ((ux * wy - uy * wx) * wind < 0.0) ? 1 : 0;
        }
        forma_seedline(x, y, mark, forma_n, W, H, forma_ramp(0.84), set(1.0, 1.0, 1.0), 0.16);
    }

    // Banked at equal decrements of area, which is equal flow time. The
    // brightness rises with how far the flow has run rather than with the
    // index of the outline, so an outline reads the same whatever the
    // stopping test turns out to be.
    if (A <= nextBank && banked < forma_g){
        float u = 1.0 - A / A0;
        forma_outline(x, y, forma_n, W, H, forma_ramp(0.86 + 0.14 * u), 0.12 + 0.24 * u);
        nextBank -= A0 / float(forma_g + 1);
        banked++;
    }

    // The two exits the plate uses, and a ceiling that is the cook talking
    // rather than the mathematics: the slowest corner of the declared box —
    // 260 nodes on a two-lobed seed at full wobble — takes about 10,800
    // steps on the page, and a wrangle that runs that far is a wrangle
    // nobody waits for. Raise forma_max above if that corner is the one
    // you want.
    if (iso > forma_round) break;
    if (A < forma_floor * A0) break;
    if (steps == forma_max) break;

    float dt = forma_cfl * hmin * hmin;
    for (int i = 0; i < forma_n; i++){
        int a = (i == 0) ? forma_n - 1 : i - 1;
        int b = (i + 1 == forma_n) ? 0 : i + 1;
        float hm = h[a], hp = h[i];
        float k = 2.0 / (hm + hp);
        vx[i] = k * ((x[b] - x[i]) / hp - (x[i] - x[a]) / hm);
        vy[i] = k * ((y[b] - y[i]) / hp - (y[i] - y[a]) / hm);
    }
    for (int i = 0; i < forma_n; i++){ x[i] += dt * vx[i];  y[i] += dt * vy[i]; }
}

// The curve it stopped on, brightest.
forma_outline(x, y, forma_n, W, H, forma_ramp(1.02), 0.92);

// And the limit point: the centroid of what is left, which is where the
// theorem says the curve was going all along.
float gx = 0.0, gy = 0.0;
for (int i = 0; i < forma_n; i++){ gx += x[i];  gy += y[i]; }
gx /= float(forma_n);
gy /= float(forma_n);
int lim = addpoint(0, set(gx - W / 2.0, (H / 2.0) - gy, 0.0));
setpointattrib(0, "Cd", lim, set(1.0, 1.0, 1.0));
setpointattrib(0, "Alpha", lim, 0.95);