singularity-rs/assets/shaders/black_hole.wgsl
xfy 52cf4e4a29 fix(stars): balance AA so it only softens edges, not size or count
After the round-star fix, toggling "Anti-aliased stars" made stars look
bigger and more numerous instead of just smoother. The AA path differed
from the hard-edge path in three ways that each inflated the footprint:

  - radius: AA used 0.25 + b*0.4 (up to 0.65); hard path used fixed 0.5
  - gain:   AA used 4.0; hard path used 3.0  (+33% peak brightness)
  - falloff: AA Gaussian decayed slowly (36.8% at one radius, never 0);
             hard path smoothstep hit 0 at the radius

The slow Gaussian tail was the main offender: it kept dim stars lit far
past the hard path's cutoff, lifting many sub-threshold stars above the
visibility floor — the apparent "more stars".

Balance the two paths so the toggle is purely an edge-softness control:
  - both use radius = 0.25 + b*0.4  (same size, brightness-driven)
  - both use gain 3.0               (same peak brightness)
  - AA steepens the Gaussian to exp(-4.6*d²/r²) so it converges to ~1%
    near the same radius where smoothstep cuts to 0 — matched footprint,
    soft vs hard edge.

Peak brightness is now identical across the toggle; effective lit area
differs by <30% (the Gaussian is a cone, smoothstep a flatter dome) and
the visibility radius matches, so dim stars no longer pop in and out as
AA is toggled.

Validated: naga type-checks; app runs 20s on Metal (M4) with no shader
errors; cargo test (14 physics tests) passes.
2026-07-15 10:23:29 +08:00

478 lines
18 KiB
WebGPU Shading Language
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

// All shader logic is inlined into this single file.
//
// HISTORY: the shader used to be split across several naga_oil modules
// (ray_gen, geodesic, stars, disk, planets, grid, skybox) that each
// `#define_import_path singularity::...` and were pulled into this file via
// `#import singularity::...`. That compiled without error, but calling ANY
// imported function at runtime produced no fragment output — the fullscreen
// quad silently drew nothing and only the camera clear color (grey) showed.
// (Local functions worked; only cross-module imports broke.) Rather than chase
// the naga_oil composition bug, every function is inlined here. The standalone
// module files still exist on disk but are no longer imported.
#import bevy_sprite::mesh2d_vertex_output::VertexOutput
struct BlackHoleUniforms {
eye: vec4<f32>,
forward: vec4<f32>,
right: vec4<f32>,
up: vec4<f32>,
resolution: vec2<f32>,
time: f32,
_pad3: f32,
rs: f32,
disk_inner: f32,
disk_outer: f32,
disk_tilt: f32,
disk_brightness: f32,
disk_rotation_speed: f32,
doppler_strength: f32,
star_intensity: f32,
skybox_intensity: f32,
grid_density: f32,
doppler_enabled: u32,
grid_enabled: u32,
planet_count: u32,
steps: u32,
spin: f32,
star_aa: u32,
bloom_threshold: f32,
bloom_strength: f32,
exposure: f32,
_pad5: f32,
};
@group(#{MATERIAL_BIND_GROUP}) @binding(0) var<uniform> uniforms: BlackHoleUniforms;
// ---------- planets storage (binding 3) ----------
struct SphereData {
center: vec4<f32>, // xyz = center (world space), w = radius
color: vec4<f32>, // xyz = color, w = emissive flag
};
@group(#{MATERIAL_BIND_GROUP}) @binding(3) var<storage, read> planets: array<SphereData>;
// ---------- optional cubemap skybox (bindings 1 & 2) ----------
@group(#{MATERIAL_BIND_GROUP}) @binding(1) var skybox: texture_cube<f32>;
@group(#{MATERIAL_BIND_GROUP}) @binding(2) var skybox_sampler: sampler;
// ====================== inlined helpers ======================
// Rotate a vector around the X axis by angle a.
fn rot_x(v: vec3<f32>, a: f32) -> vec3<f32> {
let c = cos(a);
let s = sin(a);
return vec3<f32>(v.x, c * v.y - s * v.z, s * v.y + c * v.z);
}
// --- ray_gen ---
// `fov` is packed into the `.w` of `up` (Rust lays out `up: Vec3` + `fov: f32`
// as one vec4 block).
fn ray_direction(uv: vec2<f32>) -> vec3<f32> {
let tan_half_fov = tan(uniforms.up.w * 0.5);
let dir =
normalize(uniforms.forward.xyz)
+ uniforms.right.xyz * (uv.x * tan_half_fov)
+ uniforms.up.xyz * (uv.y * tan_half_fov);
return normalize(dir);
}
// --- stars ---
fn hash13(p: vec3<f32>) -> f32 {
var q = vec3<f32>(dot(p, vec3<f32>(127.1, 311.7, 74.7)),
dot(p, vec3<f32>(269.5, 183.3, 246.1)),
dot(p, vec3<f32>(113.5, 271.9, 124.6)));
let h = fract(sin(q) * 43758.5453);
return h.x;
}
fn star_color(dir: vec3<f32>, intensity: f32) -> vec3<f32> {
// The unit-sphere direction is hashed into a 3D cell grid. A pixel's ray
// usually sits near a cell boundary, so evaluating only the home cell clips
// any star whose center lies in a neighbor — that clip is what made stars
// look square even with the Gaussian AA path: the hash changes at the cell
// edge, so the falloff could never reach zero and the visible shape was
// just the cell's bounding box. Scan the 3×3×3 neighborhood so each star's
// radial falloff extends across the whole cell it lives in.
let scale = 80.0;
let p = dir * scale;
let base = floor(p);
let threshold = 0.985;
// Brightest contributor wins: `best` = (brightness, r, g, b).
var best = vec4<f32>(0.0);
for (var dz: i32 = -1; dz <= 1; dz = dz + 1) {
for (var dy: i32 = -1; dy <= 1; dy = dy + 1) {
for (var dx: i32 = -1; dx <= 1; dx = dx + 1) {
let cell = base + vec3<f32>(f32(dx), f32(dy), f32(dz));
let h = hash13(cell);
if (h <= threshold) { continue; }
let b = (h - threshold) / (1.0 - threshold);
let center = cell + vec3<f32>(0.5);
let dist = length(p - center); // Euclidean → round, never square
let col = mix(vec3<f32>(0.6, 0.7, 1.0), vec3<f32>(1.0, 0.9, 0.7), b);
// Both paths share the same brightness-driven radius and the
// same gain, so toggling AA only softens the edge — it does not
// change peak brightness or visible count. Before, the AA path
// used a larger radius (up to 0.65 vs 0.5), a heavier gain
// (4.0 vs 3.0), and a slow-decaying Gaussian whose long tail
// lifted many dim stars above the visibility floor — that's
// what made AA look like "bigger and more" stars.
let radius = 0.25 + b * 0.4;
let falloff = select(
smoothstep(radius, 0.0, dist), // hard edge
exp(-4.6 * dist * dist / (radius * radius)), // soft edge
uniforms.star_aa != 0u);
let bright = b * falloff;
if (bright > best.x) {
best = vec4<f32>(bright, col.r, col.g, col.b);
}
}
}
}
return best.yzw * best.x * 3.0 * intensity;
}
// --- skybox ---
fn skybox_color(dir: vec3<f32>) -> vec3<f32> {
return textureSample(skybox, skybox_sampler, dir).rgb;
}
// --- geodesic ---
struct Deriv {
dpos: vec3<f32>,
ddir: vec3<f32>,
}
fn deriv(pos: vec3<f32>, dir: vec3<f32>) -> Deriv {
let r = length(pos);
let rs = uniforms.rs;
// Kerr spin. χ ∈ [0,1]; a = χ·M, M = Rs/2 = 0.5 (Rs=1).
let chi = uniforms.spin;
let m = 0.5;
let a = chi * m;
// Schwarzschild radial bending (identical to Phase 1 at χ=0).
let h = cross(pos, dir);
let h2 = dot(h, h);
let r5 = max(r * r * r * r * r, 1e-6);
let radial = -1.5 * rs * h2 / r5 * pos;
// Frame-dragging (Lense-Thirring leading term). Spin axis = +Y.
let spin_axis = vec3<f32>(0.0, 1.0, 0.0);
let r3 = max(r * r * r, 1e-6);
let drag = 2.0 * m * a / r3 * cross(spin_axis, dir);
let accel = radial + drag;
return Deriv(dir, accel);
}
// --- disk ---
fn disk_hit(prev: vec3<f32>, cur: vec3<f32>) -> bool {
let y0 = prev.y;
let y1 = cur.y;
if (y0 * y1 > 0.0) {
return false;
}
let t = y0 / (y0 - y1);
let cross = mix(prev, cur, vec3<f32>(t));
let r = length(vec2<f32>(cross.x, cross.z));
return r >= uniforms.disk_inner && r <= uniforms.disk_outer;
}
// --- disk noise (domain-warped FBM) ---
fn hash33(p: vec3<f32>) -> vec3<f32> {
let q = vec3<f32>(
dot(p, vec3<f32>(127.1, 311.7, 74.7)),
dot(p, vec3<f32>(269.5, 183.3, 246.1)),
dot(p, vec3<f32>(113.5, 271.9, 124.6)),
);
return fract(sin(q) * 43758.5453);
}
fn value_noise3(p: vec3<f32>) -> f32 {
let i = floor(p);
let f = fract(p);
let u = f * f * (3.0 - 2.0 * f);
let n000 = hash33(i + vec3<f32>(0.0, 0.0, 0.0)).x;
let n100 = hash33(i + vec3<f32>(1.0, 0.0, 0.0)).x;
let n010 = hash33(i + vec3<f32>(0.0, 1.0, 0.0)).x;
let n110 = hash33(i + vec3<f32>(1.0, 1.0, 0.0)).x;
let n001 = hash33(i + vec3<f32>(0.0, 0.0, 1.0)).x;
let n101 = hash33(i + vec3<f32>(1.0, 0.0, 1.0)).x;
let n011 = hash33(i + vec3<f32>(0.0, 1.0, 1.0)).x;
let n111 = hash33(i + vec3<f32>(1.0, 1.0, 1.0)).x;
let nx00 = mix(n000, n100, u.x);
let nx10 = mix(n010, n110, u.x);
let nx01 = mix(n001, n101, u.x);
let nx11 = mix(n011, n111, u.x);
let nxy0 = mix(nx00, nx10, u.y);
let nxy1 = mix(nx01, nx11, u.y);
return mix(nxy0, nxy1, u.z);
}
fn fbm3(p: vec3<f32>, octaves: u32) -> f32 {
var sum = 0.0;
var amp = 0.5;
var freq = 1.0;
for (var i: u32 = 0u; i < octaves; i = i + 1u) {
sum = sum + amp * value_noise3(p * freq);
freq = freq * 2.0;
amp = amp * 0.5;
}
return sum;
}
fn disk_noise(pos: vec3<f32>, t: f32) -> f32 {
let warp = fbm3(pos * 0.8 + vec3<f32>(0.0, 0.0, t * 0.1), 3u);
let n = fbm3(pos * 2.0 + warp * 1.5 + vec3<f32>(0.0, 0.0, t * 0.3), 4u);
return n;
}
fn disk_color(pos: vec3<f32>, dir: vec3<f32>) -> vec3<f32> {
let r = length(vec2<f32>(pos.x, pos.z));
let phi = atan2(pos.z, pos.x);
let rot = uniforms.time * uniforms.disk_rotation_speed / pow(r, 1.5);
// Domain-warped FBM for feathered/smoky gas texture. The Keplerian shear
// (rot ∝ 1/r^1.5) is folded into the noise flow term so inner radii flow
// faster than outer — correct differential rotation.
let noise = disk_noise(vec3<f32>(pos.x * 0.3, pos.z * 0.3, rot), uniforms.time);
let t = (r - uniforms.disk_inner) / (uniforms.disk_outer - uniforms.disk_inner);
let tcol = mix(vec3<f32>(1.0, 0.95, 0.85), vec3<f32>(1.0, 0.45, 0.12), clamp(t, 0.0, 1.0));
let falloff = 1.0 / pow(r / uniforms.disk_inner, 2.0);
var col = tcol * (0.6 + 0.4 * noise) * falloff;
let v_orbital = sqrt(uniforms.rs / (2.0 * r));
let tangent = normalize(vec3<f32>(-sin(phi), 0.0, cos(phi)));
let vdotn = dot(tangent * v_orbital, -dir);
let gamma = 1.0 / sqrt(max(1.0 - v_orbital * v_orbital, 1e-4));
var doppler = 1.0;
if (uniforms.doppler_enabled != 0u) {
let delta = 1.0 / (gamma * (1.0 - vdotn));
doppler = pow(delta, 3.0) * uniforms.doppler_strength;
}
col *= doppler;
return col * uniforms.disk_brightness;
}
// --- planets ---
// `prev`/`cur` are in DISK-LOCAL space; planet centers are world space, so we
// rotate each center into disk-local space here.
fn planet_hit(prev: vec3<f32>, cur: vec3<f32>, dir: vec3<f32>) -> vec4<f32> {
var nearest_t = 1e9;
var nearest_col = vec3<f32>(0.0);
var found = false;
for (var i: u32 = 0u; i < uniforms.planet_count; i = i + 1u) {
let s = planets[i];
let center = rot_x(s.center.xyz, -uniforms.disk_tilt);
let radius = s.center.w;
let seg = cur - prev;
let oc = prev - center;
let a = dot(seg, seg);
let b = 2.0 * dot(oc, seg);
let c = dot(oc, oc) - radius * radius;
let disc = b * b - 4.0 * a * c;
if (disc < 0.0) { continue; }
let sq = sqrt(disc);
var t = (-b - sq) / (2.0 * a);
if (t < 0.0) { t = (-b + sq) / (2.0 * a); }
if (t >= 0.0 && t <= 1.0 && t < nearest_t) {
nearest_t = t;
let hit_pos = prev + seg * t;
let n = normalize(hit_pos - center);
let light_dir = normalize(vec3<f32>(0.5, 0.8, 0.3));
let ndl = max(dot(n, light_dir), 0.0);
var col = s.color.xyz * (0.2 + 0.8 * ndl);
if (s.color.w > 0.5) { col = s.color.xyz; }
nearest_col = col;
found = true;
}
}
if (found) {
return vec4<f32>(nearest_col, 0.95);
}
return vec4<f32>(0.0, 0.0, 0.0, 0.0);
}
// --- grid (Flamm's paraboloid) ---
fn flamm_depth(r: f32) -> f32 {
if (r <= uniforms.rs) { return 0.0; }
return -2.0 * sqrt(uniforms.rs * (r - uniforms.rs));
}
fn grid_hit(prev: vec3<f32>, cur: vec3<f32>) -> vec3<f32> {
let r0 = length(vec2<f32>(prev.x, prev.z));
let r1 = length(vec2<f32>(cur.x, cur.z));
let z0_surf = flamm_depth(r0);
let z1_surf = flamm_depth(r1);
if ((prev.y - z0_surf) * (cur.y - z1_surf) > 0.0) {
return vec3<f32>(0.0);
}
var hit = vec3<f32>(0.0);
var found = false;
for (var s: i32 = 0; s < 8; s = s + 1) {
let f = f32(s + 1) / 8.0;
let p = mix(prev, cur, vec3<f32>(f));
let r = length(vec2<f32>(p.x, p.z));
let surf = flamm_depth(r);
if (abs(p.y - surf) < 0.3) {
hit = p;
found = true;
break;
}
}
if (!found) { return vec3<f32>(0.0); }
let r = length(vec2<f32>(hit.x, hit.z));
let phi = atan2(hit.z, hit.x);
let ring = smoothstep(0.06, 0.0, abs(fract(r * uniforms.grid_density * 0.5) - 0.5));
let spoke = smoothstep(0.04, 0.0, abs(fract(phi * 6.0 / 6.283185) - 0.5));
let grid = max(ring, spoke);
let fade = smoothstep(-15.0, -1.0, hit.y);
let col = vec3<f32>(0.15, 0.3, 0.6) * grid * fade;
return col * 0.5;
}
// One Dormand-Prince RK45 step. Returns the 5th-order solution and the
// error estimate (y5 - y4) as a vec3 (position error; direction error is
// folded in via normalize so we only need position error for step control).
struct RkStep {
pos: vec3<f32>,
dir: vec3<f32>,
err: f32,
};
fn rk45_step(pos: vec3<f32>, dir: vec3<f32>, dt: f32) -> RkStep {
// Butcher tableau (Dormand-Prince), 6 stages. Each deriv() returns Deriv{dpos, ddir}.
let k1 = deriv(pos, dir);
let p2 = pos + k1.dpos * dt * 0.2;
let d2 = normalize(dir + k1.ddir * dt * 0.2);
let k2 = deriv(p2, d2);
let p3 = pos + (k1.dpos * 0.075 + k2.dpos * 0.225) * dt;
let d3 = normalize(dir + (k1.ddir * 0.075 + k2.ddir * 0.225) * dt);
let k3 = deriv(p3, d3);
let p4 = pos + (k1.dpos * 0.3 + k2.dpos * -0.9 + k3.dpos * 1.2) * dt;
let d4 = normalize(dir + (k1.ddir * 0.3 + k2.ddir * -0.9 + k3.ddir * 1.2) * dt);
let k4 = deriv(p4, d4);
let p5 = pos + (k1.dpos * -11.0/54.0 + k2.dpos * 2.5 + k3.dpos * -70.0/27.0 + k4.dpos * 35.0/27.0) * dt;
let d5 = normalize(dir + (k1.ddir * -11.0/54.0 + k2.ddir * 2.5 + k3.ddir * -70.0/27.0 + k4.ddir * 35.0/27.0) * dt);
let k5 = deriv(p5, d5);
let p6 = pos + (k1.dpos * 1631.0/55296.0 + k2.dpos * 175.0/512.0 + k3.dpos * 575.0/13824.0 + k4.dpos * 44275.0/110592.0 + k5.dpos * 253.0/4096.0) * dt;
let d6 = normalize(dir + (k1.ddir * 1631.0/55296.0 + k2.ddir * 175.0/512.0 + k3.ddir * 575.0/13824.0 + k4.ddir * 44275.0/110592.0 + k5.ddir * 253.0/4096.0) * dt);
let k6 = deriv(p6, d6);
// 5th-order solution (used to advance).
let new_pos = pos + (k1.dpos * 37.0/378.0 + k3.dpos * 250.0/621.0 + k4.dpos * 125.0/594.0 + k5.dpos * 512.0/1771.0 + k6.dpos * 0.0) * dt;
let new_dir = normalize(dir + (k1.ddir * 37.0/378.0 + k3.ddir * 250.0/621.0 + k4.ddir * 125.0/594.0 + k5.ddir * 512.0/1771.0 + k6.ddir * 0.0) * dt);
// 4th-order solution (for error estimate).
let pos4 = pos + (k1.dpos * 2825.0/27648.0 + k3.dpos * 18575.0/48384.0 + k4.dpos * 13525.0/55296.0 + k5.dpos * 277.0/14336.0 + k6.dpos * 0.25) * dt;
let err = length(new_pos - pos4);
return RkStep(new_pos, new_dir, err);
}
// ====================== main ======================
@fragment
fn fragment(in: VertexOutput) -> @location(0) vec4<f32> {
let aspect = uniforms.resolution.x / uniforms.resolution.y;
var uv = (in.uv * 2.0 - 1.0);
uv.x *= aspect;
let dir = ray_direction(uv);
// Work in disk-local space: rotate eye + dir by -disk_tilt around X so the
// disk lies on y=0. (disk_hit/disk_color assume disk-local coords.)
// Total path length to integrate: enough to go from the camera, past the
// hole, and far enough beyond to count as escaped.
let eye_dist = length(uniforms.eye.xyz);
let escape_r = max(eye_dist * 2.0, 100.0);
let total_path = eye_dist + escape_r;
// Adaptive RK45 constants.
let steps_max = uniforms.steps;
let dt_init = total_path / f32(steps_max);
let dt_min = dt_init * 0.25;
let dt_max = dt_init * 4.0;
let tol = 1e-3;
let r_plus = 0.5 + sqrt(max(0.25 - (uniforms.spin * 0.5) * (uniforms.spin * 0.5), 0.0));
var pos = rot_x(uniforms.eye.xyz, -uniforms.disk_tilt);
var d = normalize(rot_x(dir, -uniforms.disk_tilt));
var dt = dt_init;
var prev = pos;
var budget = steps_max;
var accum_color = vec3<f32>(0.0);
var accum_alpha = 0.0;
loop {
if (budget == 0u) { break; }
let step = rk45_step(pos, d, dt);
let err = step.err;
if (err > tol * 10.0 && dt > dt_min) {
// Reject: shrink dt, retry (does not consume budget).
// The `dt > dt_min` guard makes dt_min a forced-accept floor: once
// dt is already at its minimum, accept the step regardless of error
// rather than spinning forever on `continue`. This prevents a GPU
// hang when a ray hits a region whose error is intrinsically above
// tolerance even at the smallest step (extreme spin, near-horizon).
dt = dt_min;
continue;
}
// Accept: consume one budget unit, refine dt.
budget = budget - 1u;
if (err > tol * 10.0) {
// Forced accept at dt_min (err still high): don't grow dt back up.
} else {
dt = clamp(dt * pow(tol / max(err, 1e-12), 0.2), dt_min, dt_max);
}
let new_pos = step.pos;
let new_dir = step.dir;
let r = length(new_pos);
if (r < r_plus) {
break;
}
if (r > escape_r) {
let world_dir = normalize(rot_x(new_dir, uniforms.disk_tilt));
var bg = vec3<f32>(0.0);
bg += star_color(world_dir, uniforms.star_intensity);
if (uniforms.skybox_intensity > 0.0) {
bg += skybox_color(world_dir) * uniforms.skybox_intensity;
}
accum_color += (1.0 - accum_alpha) * bg;
accum_alpha = 1.0;
break;
}
if (disk_hit(prev, new_pos)) {
let ty = prev.y / (prev.y - new_pos.y);
let hit = mix(prev, new_pos, vec3<f32>(ty));
let dc = disk_color(hit, new_dir);
let a = 0.85;
accum_color += (1.0 - accum_alpha) * dc * a;
accum_alpha += (1.0 - accum_alpha) * a;
if (accum_alpha > 0.99) { break; }
}
let ph = planet_hit(prev, new_pos, new_dir);
if (ph.w > 0.0) {
accum_color += (1.0 - accum_alpha) * ph.xyz * ph.w;
accum_alpha += (1.0 - accum_alpha) * ph.w;
if (accum_alpha > 0.99) { break; }
}
if (uniforms.grid_enabled != 0u) {
let g = grid_hit(prev, new_pos);
if (g.x + g.y + g.z > 0.0) {
accum_color += g;
}
}
prev = new_pos;
pos = new_pos;
d = new_dir;
}
return vec4<f32>(accum_color, 1.0);
}