Space Colonisation

DEFINITION

each attractor pulls its nearest node
nodes step toward the mean pull
attractors die inside the kill radius

NOTES

Scatter points where a plant could grow, then let the branches compete for them: every attractor tugs on its nearest node, each node steps toward the average tug, and any attractor approached closely enough is consumed. That is the whole algorithm, and it produces leaf venation, tree crowns and root systems depending only on how the attractors are scattered. Growth claims territory — hence colonisation — and when the field is spent, the plate scatters a fresh one.

PROVENANCE

Origin

A. Runions, B. Lane & P. Prusinkiewicz, "Modeling Trees with a Space Colonization Algorithm", Eurographics Workshop on Natural Phenomena, 2007

Standing

Public domain — implemented from the published description

Constants

The attractor count is locked: it is the cost of every growth step

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. 62 · SPACE COLONISATION — Adam Runions, Brendan Lane & Przemysław Prusinkiewicz, 2007
//   each attractor pulls its nearest node
//   nodes step toward the mean pull
//   attractors die inside the kill radius
// 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=venation

float p_pts  = 300 + chf('pts_tweak');      // attractors · live 150 .. 500
float p_step = 6 + chf('step_tweak');       // growth step · live 4 .. 10
float p_kill = 1.8 + chf('kill_tweak');     // kill radius (× step) · live 1 .. 4

// The plate's own colour: FORMA's LATTICES 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.46 + 0.5 * cos(6.28318530718 * (t + 0)),
    0.4185 + 0.4549 * cos(6.28318530718 * (t + 0.05)),
    0.1389 + 0.151 * cos(6.28318530718 * (t + 0.1)));
}

// Runions, Lane & Prusinkiewicz's space colonisation, grown to its finish
// in one cook: every attractor tugs on its nearest node, each node steps
// toward the average tug, and any attractor an advancing tip comes close
// enough to is consumed. The plate iterates twice a frame; this iterates
// until the field is spent or the tips stall against each other, which is
// the finished tree the plate holds on screen before rescattering.
//
// The domain is the plate's own canvas at its 560-unit reference width
// (the width its line weights were tuned against), root at the bottom
// centre, y negated onto Houdini's up so the tree grows upward. Branches
// are two-point polyline prims, parent to child, coloured by the plate's
// three growth bands with its stroke alpha carried as Alpha; the attractors
// the finished tree never claimed are dim points carrying the plate's own
// underlay alpha. Deterministic: the scatter is random(counted seed) on the
// plate's first-epoch seed, so every cook is the same tree.
float W = 560.0, H = 560.0;
int forma_maxnodes = 3600;

int rc = 2007;                   // seeded(2007 + epoch·4241) at epoch 0
int npts = int(rint(p_pts));
float ax[], ay[];
for (int i = 0; i < npts; i++){
    push(ax, W * (0.08 + 0.84 * random(rc)));  rc++;
    push(ay, H * (0.06 + 0.68 * random(rc)));  rc++;
}
float nx[], ny[];
int   npar[];
push(nx, W / 2.0);  push(ny, H * 0.94);  push(npar, -1);

float stp = p_step, kill2 = (stp * p_kill) * (stp * p_kill);
int stall = 0;
while (len(ax) > 0 && stall <= 8 && len(nx) < forma_maxnodes){
    int ns = len(nx);
    float px[], py[];
    int   pc[];
    resize(px, ns);  resize(py, ns);  resize(pc, ns);
    // every attractor pulls its nearest node, unit-normalised — the mean of
    // directions, not of offsets, which is what keeps thin twigs straight
    for (int a = 0; a < len(ax); a++){
        int bi = 0;
        float bd = 1e30;
        for (int i = 0; i < ns; i++){
            float dx = ax[a] - nx[i], dy = ay[a] - ny[i];
            float d2 = dx * dx + dy * dy;
            if (d2 < bd){ bd = d2;  bi = i; }
        }
        float d = sqrt(bd);
        if (d == 0.0) d = 1.0;
        px[bi] += (ax[a] - nx[bi]) / d;
        py[bi] += (ay[a] - ny[bi]) / d;
        pc[bi]++;
    }
    int grew = 0;
    for (int i = 0; i < ns; i++){
        if (!pc[i]) continue;
        float m = sqrt(px[i] * px[i] + py[i] * py[i]);
        if (m == 0.0) m = 1.0;
        push(nx, nx[i] + px[i] / m * stp);
        push(ny, ny[i] + py[i] / m * stp);
        push(npar, i);
        grew = 1;
    }
    if (grew){
        // attractors die inside the kill radius of the nodes just grown
        float kx[], ky[];
        for (int a = 0; a < len(ax); a++){
            int dead = 0;
            for (int i = ns; i < len(nx); i++){
                float dx = ax[a] - nx[i], dy = ay[a] - ny[i];
                if (dx * dx + dy * dy < kill2){ dead = 1;  break; }
            }
            if (!dead){ push(kx, ax[a]);  push(ky, ay[a]); }
        }
        ax = kx;  ay = ky;
        stall = 0;
    } else stall++;
}

// the field the finished tree never claimed, as the plate leaves it showing
for (int a = 0; a < len(ax); a++){
    int pt = addpoint(0, set(ax[a] - W / 2.0, (H / 2.0) - ay[a], 0.0));
    setpointattrib(0, "Cd", pt, forma_ramp(0.9));
    setpointattrib(0, "Alpha", pt, 0.3);
}
// the tree: one prim per branch segment, banded through the growth's thirds
int N = len(nx);
for (int i = 1; i < N; i++){
    int band = min(2, i * 3 / N);
    vector c = forma_ramp(0.84 + 0.1 * float(band));
    int pr = npar[i];
    int p0 = addpoint(0, set(nx[pr] - W / 2.0, (H / 2.0) - ny[pr], 0.0));
    int p1 = addpoint(0, set(nx[i]  - W / 2.0, (H / 2.0) - ny[i],  0.0));
    setpointattrib(0, "Cd", p0, c);
    setpointattrib(0, "Cd", p1, c);
    setpointattrib(0, "Alpha", p0, 0.85);
    setpointattrib(0, "Alpha", p1, 0.85);
    addprim(0, "polyline", p0, p1);
}