What it does. Build raymarched space phenomena in Three.js. Use for black-hole lensing, accretion disks, wormhole throat transits, curved-ray and null-geodesic integration, lensed celestial spheres, procedural star fields and galactic skies, relativistic-looking distortion, bounded volumetric structures, and GPU effects that need controlled numerical integration. Part of scottstts/Threejs-Awesome-Graphics-Agent-Skills (scottstts/Threejs-Awesome-Graphics-Agent-Skills).


Upstream	scottstts/Threejs-Awesome-Graphics-Agent-Skills
Skill file	skills/threejs-raymarched-space-effects/SKILL.md
License	MIT
Author	Scott Sun (scottstts)
Fetched	2026-09-10

Install

npx skills add scottstts/Threejs-Awesome-Graphics-Agent-Skills --skill threejs-raymarched-space-effects, or copy the skill folder into ~/.claude/skills/threejs-raymarched-space-effects/.
Raw file: curl -sL https://raw.githubusercontent.com/scottstts/Threejs-Awesome-Graphics-Agent-Skills/HEAD/skills/threejs-raymarched-space-effects/SKILL.md

SKILL.md (verbatim)

name: threejs-raymarched-space-effects
description: Build raymarched space phenomena in Three.js. Use for black-hole lensing, accretion disks, wormhole throat transits, curved-ray and null-geodesic integration, lensed celestial spheres, procedural star fields and galactic skies, relativistic-looking distortion, bounded volumetric structures, and GPU effects that need controlled numerical integration.

Raymarched Space Effects

Treat these effects as numerical renderers with explicit integration state. The visual character depends on coordinate choice, step policy, and how rays interact with emissive structures.

This skill contains exemplary examples and assets beyond descriptive guidance, they're worth studying, referencing, or even copying. Use them sufficiently when relevant and do NOT blindly skip them.

Workflow

Define the effect-space transform and camera ray.
Choose a physical, physically inspired, or purely artistic bending model.
Bound the integration domain.
Track ray position, direction, throughput, and accumulated radiance.
Detect crossings with disks, shells, throats, or event boundaries.
Sample the background only after integration terminates, filtered by the ray-bundle footprint rather than point sampled.
Add diagnostics for trajectory, step count, and termination reason.

Read references/curved-ray-integrators.md for the two wormhole throat models, observer transport through a throat, the artistic curved-ray accretion integrator, disk composition, and implementation defects.

Read references/lensed-celestial-spheres.md for footprint-filtered skies, flux-conserving star point spreads, dust extinction, bright sources under display clipping, progressive accumulation, and the bloom kernel a compact source needs.

Read the curved-ray accretion volume for the inverse-square steering loop, thin disk density, front-to-back accumulation, deterministic star environment, and integrator diagnostics.

Read the Schwarzschild geodesic black-hole effect for RK2 null-geodesic integration, interpolated equatorial disk crossings, Doppler/redshift disk emission, lensed procedural deep field, and HDR bloom composition.

Read the traversable wormhole transit for adaptive RK4 integration of the reduced geodesic system, an observer frame parallel-transported through the throat, exterior-region selection, progressive Halton accumulation, and the 13-tap bloom pyramid; its celestial spheres are a standalone GLSL chunk for the footprint-filtered galactic sky, resolved star layers, and the analytic ringed planet a lensed ray lands on.

Constraints

Do not call a UV swirl “gravitational lensing.”
Cap iterations and provide early termination.
Use continuous crossing tests for thin structures.
Keep numerical stability independent from frame rate.
Separate the integrator from shading of the accretion disk or wormhole interior.
Size the step from the local curvature scale, not from a global constant, once the metric has a narrow feature such as a lensing shoulder.
Filter the background by the ray-bundle footprint. A lensed map compresses solid angle enough that a point sample aliases where the image is most interesting.
Give capped rays a defined result. Their exit direction is arbitrary, so route them to a mean-radiance path rather than letting them speckle.
Carry a transported frame, not a Cartesian camera, when the domain has no global Cartesian chart.
Provide a cheaper approximation for non-hero views.

Routing boundary

Use $threejs-procedural-vfx for ordinary particles, trails, plasma, and event effects. This skill is for per-pixel numerical ray integration through curved or bounded space-effect domains.

Progressive accumulation and the bloom pyramid shipped with the wormhole example are owned by that renderer because the integration cost forces them. For a scene-wide HDR bloom pass over ordinary geometry use $threejs-bloom, and for exposure metering, tone-map ownership, and LUT grading use $threejs-exposure-color-grading.

Other files in this skill

references/curved-ray-integrators.md (verbatim)

Curved-ray numerical integrators

Use this reference for two wormhole throat integrators and an artistic inverse-square curved-ray accretion volume, including state reduction, observer transport, disk accumulation, background lensing, and numerical diagnostics. For metric-derived Schwarzschild null-geodesic integration with interpolated disk crossings and Doppler-shifted emission, read the schwarzschild-geodesic-black-hole example directly. For what a wormhole ray lands on once it escapes, read lensed-celestial-spheres.md.

Two throat models

Both models are spherically symmetric and ultrastatic, so both reduce every 3D ray to a planar problem with one conserved quantity. They differ in the shape function r(l) and therefore in step policy and cost:

Fixed-radius throat. r(l) = sqrt(l^2 + Rth^2): one smooth hyperbolic neck, no free lensing width. A uniform step is adequate everywhere. Cheap, and the distortion is broad and soft.
Shaped throat. A cylindrical neck of half-length a joined to lensing shoulders of width W, so the throat radius is exactly constant across the neck and the lensing is concentrated in a band you control independently of the throat size. This is what produces a razor-thin rim around a crystal-ball sphere instead of a soft fisheye. It needs an adaptive step, because the curvature scale collapses from ~r far away to M at the shoulder. The traversable-wormhole-transit example implements this model.

Constants below belong to one model or the other; do not mix them.

Fixed-radius throat: state reduction

Reduce each 3D ray to a two-dimensional integration state:

y.x = signed radial coordinate l
y.y = radial momentum pL
impact parameter b = length(cross(rayOrigin, rayDirection))
throat radius Rth = 1.2

Construct an orbital plane:

normal = normalize(cross(origin, direction))
u = normalize(origin)
v = cross(normal, u)

Near-radial rays use fallback axes to avoid a zero cross product.

Initial signed coordinate:

l = sqrt(max(length(origin)^2 - Rth^2, 0.001))
pL = dot(normalize(origin), direction)

Fixed-radius throat: RK4 integration

The derivative is:

r2 = l^2 + Rth^2
dl/ds = r2 * pL
dpL/ds = b^2 * l / r2

Run fourth-order Runge–Kutta with:

maximum iterations = 920
base step = 0.0042
per-ray step jitter = +/- 0.00045
escape distance = abs(l) > 40
azimuth accumulation = step * b

On escape:

finalDirection =
  normalize(u * cos(phi) + v * sin(phi))

The sign of final l selects which exterior universe is visible. Failure to escape renders a bright fallback color, making iteration-cap pixels observable.

This is materially stronger than a UV swirl because the final environment direction comes from numerical integration.

Shaped throat: metric and observer frame

The traversable-wormhole-transit example integrates

ds^2 = -dt^2 + dl^2 + r(l)^2 (dtheta^2 + sin^2 theta dphi^2)

with the shape function and its derivative:

x       = 2 (|l| - a) / (pi M)
r(l)    = rho                                          for |l| <= a
r(l)    = rho + M (x atan(x) - 0.5 ln(1 + x^2))        otherwise
dr/dl   = 0                                            for |l| <= a
dr/dl   = sign(l) (2/pi) atan(x)                       otherwise

l is proper radial distance and is negative on the far side: one global coordinate runs straight through the neck into the other exterior region, so there is no handoff and no second chart. The lensing parameter follows from the lensing width by W = 1.42953 M. The example fixes

rho = 1.0            (the unit of length)
W / rho = 0.05       ->  M = 0.05 / 1.42953 = 0.03497...
2a / rho = 0.01      ->  a = 0.005

Because the metric is ultrastatic and spherically symmetric, every null geodesic stays in the plane spanned by the observer's radial vector and the ray's tangential component. That reduction is exact, removes the polar coordinate singularity, and leaves three variables with one conserved quantity b = r^2 dpsi/dt:

y      = (l, psi, p_l)
dl/dt   = p_l
dpsi/dt = b / r^2
dp_l/dt = b^2 (dr/dl) / r^3

The observer cannot be a Cartesian camera, because no global Cartesian chart survives the neck. Carry a point on the sphere plus a tangent frame instead:

U = unit sphere point       A, B = tangent frame,  A x B = U
l = signed radial coordinate
yaw, pitch -> uCamRot, whose columns are the camera x/y/z axes
              written in (e_l, A, B) components

Per ray, seeded from that frame:

n     = uCamRot * normalize(vec3(px 2 tan(fov/2), py 2 tan(fov/2), -1))
m     = length(n.yz)                       tangential magnitude
tHat  = normalize(n.y A + n.z B)           (fall back to A when m < 1e-9)
b     = r(l_cam) * m                       conserved angular momentum
y     = (l_cam, 0, n.x)

Motion is a step in l plus a parallel transport of the frame along a great circle. For a displacement of length d whose components in (e_l, A, B) are (cl, ca, cb):

l += cl d
t  = normalize(ca A + cb B)
rotate U, A, B about (U x t) by  d |(ca, cb)| / max(r(l), 1e-4)
re-orthogonalise: U normalized, A -= (A.U) U, B = U x A

Speed must scale with the local sphere radius or the neck is impassably slow and the far field is impassably large:

v = speed dt clamp(0.35 r(l), 0.10, 14.0)      speed = 1.0, x4 while sprinting

Shaped throat: adaptive RK4 integration

The step is chosen from the local curvature scale: M at the shoulder, ~r far away, and unbounded inside the neck, where r is constant so the solution is linear and RK4 is exact for any h.

k = 0.15                                     step coefficient
|l| <= a : h = min((a - |l|) / max(|p_l|, 1e-4) + k M, 40 (a + rho))
|l| >  a : h = k min(r, M + 0.9 (|l| - a))
h = max(h, 1e-5)

Standard RK4 on y with the derivative above, capped at 1024 iterations. k = 0.15 holds the worst case — a ray grazing the lensing rim — within about 0.13 px of a 40000-step solution at 1080p and a 50 degree field of view, for roughly 25-98 steps per ray.

Termination, tested before each step:

escaped: r > 260 and (dr/dl) p_l > 0        asymptotically flat and receding
escaped: abs(l) > 1e7                       numerical runaway, treat as escaped
capped : the iteration limit was reached

The asymptotic direction is rebuilt in the frame carried around by the swept angle psi, from the radial component (dr/dl) p_l and the tangential component b / r:

uf =  U cos(psi) + tHat sin(psi)
tf = -U sin(psi) + tHat cos(psi)
D  = (dr/dl) p_l uf + (b / r) tf            normalize; fall back to uf if degenerate

Universe selection

The sign of l at termination selects the exterior region, and each region owns an independent directional field:

A cheap pairing is two procedural directional fields with five-octave FBM, animated coordinate drift, broad/fine structure, and different plane orientation. A small direction-hashed grain of amplitude 0.01 reduces visible gradient banding.
The traversable-wormhole-transit example instead evaluates one galaxy model twice, under different galactic frames, palettes and seeds, and filters it by the ray-bundle footprint. See lensed-celestial-spheres.md.

Whichever field is used, sample it only after integration. Lensing must change the lookup direction rather than distort an already rendered screen image.

Orient the far region deliberately. A ray fired straight down the neck exits along one fixed axis, so put that region's brightest structure on that axis; otherwise the throat reads as a black hole rather than a window.

Accretion-volume integration

The curved-ray-accretion-volume example is evaluated on a sphere surrounding the effect. Defaults:

iterations = 128
step = 0.0071
ray jitter = 0.01
bending power = 0.3
core radius = 0.13
disk half-width = 0.03

Per step:

r = length(rayPosition)
steerMagnitude = step * power / r^2
steerRange = remapClamped(r, 1 -> 0.5, 0 -> 1)
newDirection = normalize(direction - radial * steerMagnitude * steerRange)

The ray direction is bent toward the center only inside the configured range. Unlike the wormhole, this is an artistic inverse-square steering field, not a validated metric geodesic.

Disk density and color

Disk coordinates rotate around the local Z axis with radius and time:

rotation phase = radialDistance * 4.27 - time * 0.1
noise UV = rotatedPosition * 2

A repeated deep-noise texture is multiplied by a quadratic band across [-width, 0, +width]. Radial distance, noise value, and a nearby noise sample produce a ramp coordinate.

The three-point color ramp is:

white-hot at 0.06
gold at 0.33
dark amber at 1.0
emission scale 1.95
additional emission color (1.0, 0.72, 0.26)

The central sphere below originRadius is black. Disk opacity accumulates front-to-back into alphaAcc; color accumulates using remaining transmittance.

Background lensing

After the fixed loop, the final bent direction samples a deterministic equirectangular star texture generated from 5200 seeded stars on a sphere.

final = accumulated disk color
      + remaining transmittance * star environment

The deterministic star field is important: lensing motion can be compared without random stars moving between runs.

Observed defects and boundaries

The demonstrated accretion-volume loop advances rayPosition twice per iteration while its steering magnitude uses a single step. Treat the effective distance step as 2 * step, or remove the duplicate advance and retune the complete density-and-bending system explicitly.
The accretion volume has no early exit and no termination IDs; every pixel pays the full iteration count.
Its disk is detected by local band density at samples, not by a continuous plane-crossing test. A large step can skip a thin disk.
The artistic inverse-square steering must not be described as general relativity.
The fixed-radius throat uses a fixed high iteration cap and per-ray step jitter but has no CPU reference-ray tests.
Rays that wind the photon sphere at |l| ~ 0 never escape. Their exit direction is wherever the integrator happened to be pointing, which speckles pixel to pixel. They sample effectively the whole sphere, so hand them the maximum footprint and let filtering return the mean radiance instead of a point sample; a flat fallback color is a diagnostic, not a shipping image.
Raising the iteration cap does not remove those pixels. It only moves the boundary; the winding set is measure-nonzero at any cap.
The shaped throat's step policy is tuned to rho = 1. Changing rho, W, or a changes the curvature scale, so re-measure the step coefficient against a high-step solution instead of assuming 0.15 still holds.
The shaped throat is integrated in a plane per pixel with no CPU reference rays either; parity claims need an independent CPU integrator first.
Both accretion effects render on bounded proxy geometry; coordinate transforms must be verified if the proxy is nonuniformly scaled or moved far from the origin.
The star texture is finite-resolution and can alias under extreme magnification.

Diagnostics

Expose:

wormhole l and pL, or (l, psi, p_l) for the shaped throat
impact parameter / conserved b and the orbital-plane basis
RK4 step count and escaped/capped state
final exterior side and environment direction
ray-bundle footprint used to filter the celestial sphere
accretion-volume radius and steering magnitude
effective traveled distance
disk band, noise, ramp coordinate, and local alpha
accumulated alpha and remaining transmittance
core-hit mask
final bent background direction
NaN/invalid-state mask

Read them this way: a step-count view that saturates over a wide band means the step policy has collapsed, not that the scene is expensive; a capped-ray view should be a thin ring around the throat rim, and a filled disc means the escape radius is too large for the iteration cap.

Add CPU reference rays for either wormhole before claiming physical parity, and a continuous disk-crossing variant before increasing the accretion volume's step size.

references/lensed-celestial-spheres.md (verbatim)

Lensed celestial spheres

Use this reference for what a geodesic ray lands on after it escapes: a directional sky that stays stable when lensing compresses solid angle by orders of magnitude, plus the accumulation and bloom that a sky with 400x point sources in it demands. Constants are those of the traversable-wormhole-transit example, whose celestial-spheres.js is the sky and whose wormhole-effect.js owns the accumulation and bloom passes. For the integrator that produces the directions, read curved-ray-integrators.md.

The footprint contract

Lensing conserves radiance but not solid angle. Near a lensing rim a pixel's ray bundle can cover a huge patch of sky, and one sample of a point-sampled star field there is noise, not signal. Every field below is therefore a function of direction and of the bundle's angular width:

foot = clamp(0.5 (length(dFdx(D)) + length(dFdy(D))), 0.0, 0.06)

D is the escaped direction, so the derivatives are taken through the whole integration and automatically widen wherever the map is compressive. Rays that never escaped get foot = 0.06 — the maximum — so filtering hands back mean radiance instead of a speckled point sample.

Two consumers use it directly:

detail = 1.0 - smoothstep(0.0035, 0.028, foot)     fine nebula knots fade out
w      = max(foot, 2.0e-4)                          planet limb feather
e      = clamp(foot / 0.052, 0.0, 0.35)             ring-edge feather

Star lattice and flux-conserving point spread

Stars live on a cube-face lattice, one candidate per cell, so a layer costs a 3x3 cell neighbourhood and no texture. Given a direction D, pick the dominant axis, project to that face's uv, and scan the neighbourhood:

cellAng = 1 / cells
s  = clamp(max(S0, 0.62 foot), S0, 0.5 cellAng)     spread radius
k  = S0^2 / s^2                                      flux-conserving peak
I  = flux k (exp(-0.5 d^2/s^2) + 0.005 exp(-d^2/(18 s^2)))

with S0 = 0.00030 the intrinsic angular radius of a star. The second lobe is the PSF wing: 3x the width, carrying 4.5% of the flux. Because the peak falls as the width grows, total flux is invariant — compression at a lensing rim brightens the star instead of aliasing it.

Existence and position come from hashes of the cell, never from the pixel:

h  = hash33(vec3(cell, face 17 + seed))
sd = faceDirOf(face, (cell + h.xy) cellAng)
keep if h.z <= min(baseDens discDens(sd, ...), 1.0)

Evaluate the density at the star's own direction sd, not at the pixel's D, or stars flicker in and out as the camera moves. Reject anything below I < 2.0e-4: that is under an 8-bit step and under the bloom threshold.

Past roughly one cell per pixel, the lattice cannot be resolved and the correct answer is the layer's analytic mean radiance — the field's own top mip, and the physically right answer since lensing conserves radiance:

meanFlux = lum 2.78 1.036                    E[u^-2/3] with u clamped, plus giants
meanRad  = (0.98, 0.96, 1.00) meanDens meanFlux 1.045 2pi S0^2 cells^2
result   = mix(sum, meanRad, smoothstep(0.30, 1.25, foot cells))

Three layers share one disc profile and cover three decades of population:

cells =  30, density = 0.46 sDens, luminosity = 0.60  sLum, seed + 1
cells = 104, density = 0.34 sDens, luminosity = 0.18  sLum, seed + 11
cells = 300, density = 0.22 sDens, luminosity = 0.055 sLum, seed + 23

Luminosity law, temperature, and desaturation

Flux follows the number-count law N(<m) ~ 10^(0.6 m), which inverts to F ~ u^(-2/3) for uniform u:

u     = max(hash, 0.0016)                    clamp on the power law
giant = step(0.972, hash)                    ~3% evolved red giants
flux  = lum u^(-2/3) mix(1.0, 2.3, giant)

That heavy tail is the point: a handful of dominant stars over a dense faint floor, instead of the flat field a uniform random amplitude produces.

Temperature is correlated with flux, because the main sequence is — and the giants deliberately break the correlation, which is why the brightest stars in a real sky are not all blue:

ts = clamp(0.62 hash^6.5 + 0.38 (1 - u^0.30), 0.0, 1.0)
ts = mix(ts, 0.015 + 0.05 hash, giant)
colour = blackbody(mix(2700, 24000, ts))

blackbody(T) evaluates a piecewise fit to the Planck locus, then divides by its own Rec. 709 luminance so temperature sets hue only and the luminosity function alone sets brightness. Faint stars are mixed back toward white, because below the cone threshold the eye reads a point source as colourless:

tint = mix(vec3(1.0), colour, clamp(0.35 + 0.55 log2(1 + 6 flux), 0.35, 1.0))

No twinkle. Twinkle is atmospheric; there is no atmosphere in the scenes this sky is for.

Galactic disc, dust, and nebulae

One galaxy model is instanced per exterior region with a different frame, seed, and palette. A galactic frame is a north pole plus a core direction orthogonalised against it. Order the ingredients the way light reaches the camera:

sb    = dot(D, pole)                          sine of galactic latitude
lon   = 0.42 + 0.58 smoothstep(-0.55, 0.95, dot(D, core))
disc  = exp(-0.5 sb^2 / scaleH^2)
dpl   = exp(-0.5 sb^2 / (0.55 scaleH)^2)      dust is thinner than stars
bulge = max(dot(D, core), 0)^3.2 exp(-0.5 sb^2 / 0.11^2)

Dust is a column, not a paint layer. Build it from a ridged multifractal in a domain-warped coordinate (filaments, not blobs), then attenuate with a wavelength-dependent coefficient so the band is dimmed and reddened and the dark nebulae are carved out of the light behind them:

tau = dustAmp dpl lon (0.30 bulk + 1.25 smoothstep(0.34, 0.88, fil) (0.35 + 0.65 bulk))
ext = exp(-tau (1.00, 1.24, 1.52))            A ~ 1/lambda

Unresolved starlight is the milky band plus the bulge, tinted between a cool arm colour and a warm bulge colour and multiplied by ext. Emission and reflection nebulae are gated by a region mask so they read as a few discrete complexes strung along the arms rather than an even wash, and they sit inside the disc, so on average half the dust column lies in front of them — sqrt(ext), not ext. Resolved stars are attenuated by how much of the column is in front of them, which is a proxy for depth:

layer  30 (nearest): mix(vec3(1.0), ext, 0.40)
layer 104:           mix(vec3(1.0), ext, 0.75)
layer 300 (farthest): ext

The two regions in the example differ only by frame, seed and parameters — a younger, bluer, thinner-dusted galaxy at scaleH = 0.185 and an older, dustier, redder one at scaleH = 0.235 with 1.30x star density and 1.45x star luminosity.

Analytic planet and ring plane

A ringed gas giant is ray-traced analytically against the celestial sphere, so it lenses exactly like the star field instead of being composited afterwards. Sphere of angular radius 0.052, R = sin(ang):

bq   = dot(D, C)
tS   = bq - sqrt(max(bq^2 - (1 - R^2), 0))         sphere entry
aP   = 1 - smoothstep(ang - w, ang + w, acos(bq))   footprint-feathered coverage
tR   = dot(C, axis) / dot(D, axis)                  ring-plane crossing
rad  = length(tR D - C) / R                         ring radius in planet radii

Ring optical depth is three banded regions with a ringlet modulation and a gap:

C ring: smoothstep(1.12, 1.20+e) (1 - smoothstep(1.49, 1.53+e)) 0.30
B ring: smoothstep(1.52, 1.57+e) (1 - smoothstep(1.91, 1.95+e)) 0.95
A ring: smoothstep(2.02, 2.06+e) (1 - smoothstep(2.23, 2.27+e)) 0.66
ringlets: x (0.86 + 0.14 sin(173 rad))
gap:      x mix(0.18, 1.0, smoothstep(0.0, 0.006+e, abs(rad - 2.214)))

Composite by depth, not by fixed order: compare tR against tS, put the nearer one in front, and return premultiplied colour with coverage in alpha so the caller blends once. The ring shadow on the planet and the planet shadow on the rings both come from projecting onto the light direction; skipping them is the single most visible way this object stops reading as a solid body.

Bright sources under display clipping

A sun in frame is two orders of magnitude brighter than anything else and dominates the bloom, so its rendered image must be invariant under sub-pixel motion. Flux invariance alone is not enough, because the display clips: ACES sends everything above roughly 2.5 to white, so under a 420x peak only the bottom ~0.7% of a feather is visible and the tone-mapped limb of a plain disc lives in about 0.1 px. That hard sub-pixel edge pops as it crosses the pixel grid, and the bloom prefilter inherits the popping.

Give it a point spread wide enough for the display: flat core of radius R, Gaussian limb with sigma tied to the bundle footprint, peak normalised by the profile's own integral pi (R^2 + sqrt(2 pi) R sigma + 2 sigma^2):

sunR = 0.0050
sunS = 1.2 max(foot, 1.0e-5)
sunQ = max(length(D - sunDir) - sunR, 0.0)
c   += (1.00, 0.97, 0.92) 420 (sunR^2 / (sunR^2 + 2.5066283 sunR sunS + 2 sunS^2))
       exp(-0.5 sunQ^2 / sunS^2)
c   += (1.00, 0.95, 0.88) 0.55 exp(-lenD^2 / (2 0.0333^2))     the glare halo

Total flux is exactly 420 pi R^2 in every regime, so lensing can squeeze the disc below a pixel and it dims into a wider splat instead of point-sampling in and out of existence. Use the chord length(D - sunDir), not dot: near cos ~ 1 a dot comparison spends all its precision against the ~1e-6 direction noise surviving the integrator.

Progressive accumulation

A per-pixel geodesic sky is too expensive to supersample per frame and too high-contrast to leave at one sample, so converge it over time and reset on any observer change:

weight = accumCount > 0 ? 1/(accumCount + 1) : 1
out    = mix(previous, current, weight)          weight >= 0.999 replaces outright
accumCount = min(accumCount + 1, 512)

Replacing outright at full weight matters: it is what stops a stale or NaN history from being blended into the first frame after a resize.

The jitter sequence is not a plain Halton sequence:

accumCount == 0        -> jitter (0, 0)          pixel centre
accumCount <  512      -> Halton(2), Halton(3) at index accumCount + 1, minus 0.5
accumCount >= 512      -> white noise, minus 0.5

Frame 0 after a reset is the image while the observer moves, so it must sit at the pixel centre — the mean of the converged set. Starting at Halton index 1 puts every moving frame 1/6 px below the converged still, which snaps at every start and stop, and a drag that hesitates for one frame restarts at index 2 or 3 with blend weight 1/2 or 1/3 — a sub-pixel hop that a bright compact source shows plainly. In motion the sampling phase must be stationary, not well-distributed.

Bloom pyramid for compact sources

Every pyramid level decimates by two, so a sun with a display-space limb is one or two texels wide at the levels that carry the halo. A 2x2 box — a single bilinear tap — or a 4-tap corner box re-partitions a source that small between output texels as it crosses their grid. Energy is conserved, but the shape the coarse levels hand back up morphs with sub-texel phase and the halo breathes while the camera moves.

Use the 13-tap downsample instead: five overlapping 2x2 boxes, half the weight on the centre one, DC gain exactly 1.

inner four taps at +/-1 texel : 0.125 each
centre tap                    : 4 x 0.03125
edge taps at +/-2             : 2 x 0.03125 each
corner taps at +/-2           : 1 x 0.03125 each

The rest of the chain: prefilter with a soft-knee threshold at 0.85 and knee 0.55, 5 levels each halving, then a 3x3 tent upsample added into the level below with additive blending — which means the render target must not be cleared between passes. Composite adds bloom x 0.35 of level 0 to the accumulated scene, then exposure 1.15, ACES, vignette 0.34, sRGB transfer, and 0.03 of frame-indexed grain, in that order. Grain after the transfer is deliberate; grain before it is reshaped by the curve and stops masking banding.

Observed limits

The sky is a pure function of direction, so it has no parallax. Objects on it are infinitely far away and any apparent motion comes from the integrator.
The analytic mean radiance is a mean, not a blur: it is correct once cells are smaller than a pixel, and slightly flat in the transition band around foot cells in [0.30, 1.25].
The star lattice is a cube-face grid, so cell solid angle varies by about 3^(3/2) between face centre and corner. Density is uniform per cell, not per steradian.
Rendering at a fraction of the output resolution (0.65 in the example) widens the effective footprint, which is self-consistent — the derivatives are taken in the render target — but changes the scale at which layers cross into mean radiance.
Accumulation converges sampling noise, not integration error. Every sample runs the same step policy, so a step that is too coarse averages to the same wrong answer no matter how long the camera holds still.
At 512 samples the blend weight freezes at 1/513 and the jitter switches to white noise, so the result becomes a rolling average holding a constant noise level rather than continuing to converge.
Half-float targets are required. An 8-bit accumulation target loses the entire faint-star floor after the first blend.

Diagnostics

Expose:

ray-bundle footprint, and the mask where it saturates at 0.06
per-layer point-sample vs mean-radiance blend factor
star existence mask and per-cell density before the min(..., 1) clamp
dust column tau and the extinction it produces
region/HII/reflection masks separately from their colours
planet coverage, ring optical depth, and which surface won the depth test
accumulation sample count and the current blend weight
bloom prefilter output and each pyramid level

Read them this way: if the footprint view is saturated over a wide area the integrator is diverging, not the sky; if the blend-factor view shows a hard ring the layer's cells and the render scale disagree; if the bloom halo breathes while the camera moves, the downsample kernel is the suspect, not the threshold.

Compare a converged still against a single frame at the same pose to separate sampling noise from field error — they look identical in a screenshot and completely different in motion.

Back to scottstts/Threejs-Awesome-Graphics-Agent-Skills or Agent skills.

threejs-raymarched-space-effects skill (Threejs-Awesome-Graphics-Agent-Skills)

Install

SKILL.md (verbatim)

Raymarched Space Effects

Workflow

Constraints

Routing boundary

Other files in this skill

references/curved-ray-integrators.md (verbatim)

Curved-ray numerical integrators

Contents

Two throat models

Fixed-radius throat: state reduction

Fixed-radius throat: RK4 integration

Shaped throat: metric and observer frame

Shaped throat: adaptive RK4 integration

Universe selection

Accretion-volume integration

Disk density and color

Background lensing

Observed defects and boundaries

Diagnostics

references/lensed-celestial-spheres.md (verbatim)

Lensed celestial spheres

Contents

The footprint contract

Star lattice and flux-conserving point spread

Luminosity law, temperature, and desaturation

Galactic disc, dust, and nebulae

Analytic planet and ring plane

Bright sources under display clipping

Progressive accumulation

Bloom pyramid for compact sources

Observed limits

Diagnostics