Procedural Lightning

An R&D pipeline in progress - Houdini propagator, Alembic handoff, Blender render

Status: R&D, not a finished technique. This is an active exploration, not a battle-tested production pipeline - the values below are starting points I'm still validating, not defaults proven across shots. Treat the code as a base to adapt, not a drop-in solution.

Pipeline Overview

The split is deliberate: Houdini handles the procedural simulation of the bolt's path - geometry and attributes only - and Blender handles the final render - emissive shading, bloom, compositing. The two talk through Alembic (.abc, Ogawa format), which carries curve geometry and custom point attributes across reliably.

Three Ways to Generate a Bolt

L-System fractal - recursive rewriting rules, simple to drive through Houdini's native L-System SOP. A good teaching starting point, but the branching reads too regular for production realism.

Pyro/VDB simulation - a probabilistic path lit from within by an emissive pyro shader. Gives good halos and volumetric glow but is hard to control precisely and expensive to compute - less relevant here since rendering is delegated to Blender anyway.

DLA / Dielectric Breakdown (the approach used below) - simulating an electric field across a grid, with "leaders" advancing toward the path of least resistance. The most physically accurate of the three, implementable in VEX over a potential-field grid - this is the logic behind SideFX's own Labs Lightning SOP.

The Real Structure of a Lightning Bolt

Worth understanding the physics before writing the algorithm. The stepped leader - the main channel - advances in discrete jumps rather than a continuous line, each step a probabilistic decision steered by the electric field. Dead branches are bifurcations that stop for lack of potential to keep going, and they stay visible as the characteristic residue of a real strike. The return stroke is the bright flash when the leader reaches its target and the channel lights up back along its length - very brief, one to three frames. Streamers are the glowing filaments that persist after the return stroke, fading exponentially over 8-15 frames. Splits happen when the main channel divides into two comparably strong branches, when two paths of similar potential coexist.

Spatializing the Path in 3D

Four techniques, usually combined. A VDB potential field models the electric field between source and target as a float VDB; the bolt descends the gradient with fractal noise for local deviation - the most physically correct method. A guide curve with attraction defines the general trajectory, with a Point Attract or velocity field keeping branches inside a volume around it - well suited to shot-specific art direction. SDF obstacle avoidance samples scene geometry's SDF at each step via volumesample(), letting the bolt avoid solids or seek them out (a lightning rod, a conductive surface). Target-surface scatter places weighted points on a target surface (normal + distance + noise) for surface arcs - lightning crawling across a rock, striking a specific object. In production, a guide curve for macro artistic control combined with a VDB field for micro physical behavior is the strongest combination.

The Node Network

The full chain, source through export:

[GEO_SOURCE]   [GEO_TARGET]
      |               |
  [scatter]       [scatter]
      |               |
[wrangle_init_source]  [wrangle_init_target]
           |
      [merge]  <- source + target into one geo
           |
   [for_each: wrangle_propagator]  <- iterative DLA loop
           |
   [wrangle_tortuosity]  <- high-frequency noise on position
           |
   [add (by attribute: parent_id)]  <- reconnects points into curves
           |
   [wrangle_width]  <- radius from energy + branch_depth
           |
   [resample]  <- uniform segment length
           |
   [polywire]  <- meshes the tubes from width
           |
   [normal]  <- clean point normals for Blender
           |
   [rop_alembic]

Two source points feed the system: a single scattered point on the source geometry (the initial leader - full energy, depth 0) and a single scattered point on the target (a passive sentinel the propagator searches for, flagged with branch_depth = -1). Merging them into one geo lets a single wrangle loop handle both ends.

The Propagator Wrangle

This runs inside a For-Each block (Block Begin set to Fetch Input, Block End set to Feedback Each Iteration - without that setting the loop isn't cumulative, and points created in iteration N simply don't exist for iteration N+1). Each active tip, per iteration: finds the target position, blends 60% of its previous direction with 40% attraction toward the target (inertia against sudden direction changes), adds fractal noise deviation, advances one step, decays its energy, and rolls a probability check to spawn a branch before deactivating itself as a tip.

if (i@is_tip == 0) return;

float step_len    = chf("step_length");    // default 0.15
float chaos       = chf("chaos");          // default 0.4
float freq        = chf("noise_freq");     // default 1.2
float p_branch1   = chf("prob_branch1");   // depth 0->1, default 0.25
float p_branch2   = chf("prob_branch2");   // depth 1->2, default 0.18
float energy_main = chf("energy_decay_main");   // default 0.88
float energy_b1   = chf("energy_decay_b1");     // default 0.60
float energy_b2   = chf("energy_decay_b2");     // default 0.55
float min_energy  = chf("min_energy");     // stop threshold, default 0.08
int   max_depth   = chi("max_depth");      // default 2 (3 levels: 0,1,2)
int   seed_offset = chi("seed");           // default 42

vector target_pos = {0, -5, 0};
int npts = npoints(0);
for (int k = 0; k < npts; k++) {
    if (point(0, "branch_depth", k) == -1) {
        target_pos = point(0, "P", k);
        break;
    }
}

vector to_target = normalize(target_pos - v@P);
vector prev_dir = normalize(v@dir);
vector base_dir = normalize(lerp(prev_dir, to_target, 0.4));

vector noise_offset = set(
    noise(v@P * freq + set(seed_offset * 0.1, 0, 0)),
    noise(v@P * freq + set(0, seed_offset * 0.1, 0)),
    noise(v@P * freq + set(0, 0, seed_offset * 0.1))
);
noise_offset = fit(noise_offset, set(0,0,0), set(1,1,1), set(-1,-1,-1), set(1,1,1));

vector dir = normalize(base_dir + noise_offset * chaos);
v@dir = dir;

if (f@energy < min_energy) { i@is_tip = 0; return; }

vector new_pos = v@P + dir * step_len;
int new_pt = addpoint(0, new_pos);

int   new_depth  = i@branch_depth;
float new_energy = f@energy;
if      (new_depth == 0) new_energy *= energy_main;
else if (new_depth == 1) new_energy *= energy_b1;
else                      new_energy *= energy_b2;

setpointattrib(0, "branch_depth", new_pt, new_depth);
setpointattrib(0, "is_tip",       new_pt, 1);
setpointattrib(0, "parent_id",    new_pt, i@ptnum);
setpointattrib(0, "energy",       new_pt, new_energy);
setpointattrib(0, "dir",          new_pt, dir);
i@is_tip = 0;

if (new_depth < max_depth) {
    float p_branch = (new_depth == 0) ? p_branch1 : p_branch2;
    float rng = rand(i@ptnum * 1973 + seed_offset + @Time * 100);
    if (rng < p_branch * f@energy) {
        vector up = set(0, 1, 0);
        if (abs(dot(dir, up)) > 0.95) up = set(1, 0, 0);
        vector side = normalize(cross(dir, up));
        float twist = fit01(rand(i@ptnum * 3571 + seed_offset), -1.2, 1.2);
        vector branch_dir = normalize(side * cos(twist) + cross(dir, side) * sin(twist));
        branch_dir = normalize(branch_dir + to_target * 0.3);
        int branch_pt = addpoint(0, v@P + branch_dir * step_len * 0.85);
        float branch_energy = (new_depth == 1) ? f@energy * energy_b2 : f@energy * energy_b1;
        setpointattrib(0, "branch_depth", branch_pt, new_depth + 1);
        setpointattrib(0, "is_tip",       branch_pt, 1);
        setpointattrib(0, "parent_id",    branch_pt, i@ptnum);
        setpointattrib(0, "energy",       branch_pt, branch_energy);
        setpointattrib(0, "dir",          branch_pt, branch_dir);
    }
}

The Block End's iteration count (40-80 in early tests) sets how far the leader travels and how many branch generations appear. For true splits rather than one-sided branching, the propagator can run simultaneously from source and target and connect where the two fronts meet - available in the Labs Lightning SOP's advanced mode, not yet built into this version.

Width, Tortuosity and Tubes

The macro skeleton comes out of the propagator; a separate high-frequency noise pass on point position afterward adds the characteristic fast micro-jitter of a real bolt, with branches noisier than the main channel (tort_amp_br around 0.035 against tort_amp_main around 0.02). Width is derived from energy and depth - a thicker main channel, thinner sub-branches, tapering toward zero as a branch's energy runs out, which is what sells a dead branch as actually dead rather than just stopping abruptly:

float base_w = (i@branch_depth == 0) ? chf("width_main")
             : (i@branch_depth == 1) ? chf("width_b1") : chf("width_b2");
f@width = base_w * fit(f@energy, 0.0, 1.0, 0.1, 1.0);
if (i@branch_depth == 0 && f@energy > 0.6) f@width *= 1.25; // thicker at forks

Points reconnect into curves through an Add SOP set to build polygons by the parent_id attribute - no manual wiring needed, though a noisy or disconnected result usually means a stray parent_id value slipped in outside the source/target sentinels. A Resample pass (uniform segment length, subdivision-curve treatment) cleans the curve up before PolyWire meshes it into tubes, reading radius straight from the width point attribute - PolyWire handles per-point radius and fork joins more cleanly than a Sweep here. Point normals get stripped before export and left for Blender to recompute, since imported Alembic normals tend to produce artifacts on the Blender side.

The Alembic Handoff to Blender

Export settings that matter: Ogawa format (lighter and faster than HDF5), Single Partition for one object (or partition by branch_depth for separate Blender objects per level), Build Hierarchy on, and an explicit point-attribute list - energy life is_main branch_depth width - rather than "Save All," which drags along internal Houdini attributes (pscale, id, is_tip, parent_id) that do nothing in Blender and just bloat the file.

What crosses cleanlyWhat doesn't
width - read natively as curve radiusString attributes - ignored entirely
Float/int point attributes - via Geometry Nodes' Named AttributePrimitive (not point) attributes - sometimes dropped
Curve hierarchy, if built cleanly in HoudiniAttribute names with special characters

On the Blender side, the Alembic import lands as a curve object; a Geometry Nodes setup reads the attributes back through Named Attribute nodes (matching the exact name set in Houdini) and reconstructs the tubes, ready for a multi-layer emissive shader.

What's Still Unvalidated

This is the current state of an active exploration, not a closed loop - worth being upfront about what hasn't been proven yet rather than presenting it as finished:

  • True bi-directional splits (propagating from both source and target) aren't implemented in this version
  • Return-stroke timing (an animated life attribute driving the flash) hasn't been built or tested
  • The full Alembic Ogawa → Blender 4.x attribute handoff hasn't been run end to end yet
  • The Blender-side Geometry Nodes setup for tube reconstruction from the imported curves doesn't exist yet
  • The multi-layer emissive shader (white core, blue corona, outer halo) is still just a plan

Before trusting a first test: confirm the Block End is actually set to Feedback Each Iteration (the whole loop silently does nothing cumulative without it), color points by branch_depth before the Add SOP to sanity-check the tree structure, check the Spreadsheet for stray parent_id values before connecting curves, and export a single frame to Alembic before committing to a full animated sequence.